так как вдоль рассматриваемой границы Li = 0. Теперь для тепло вого потока получаем следующее выражение:
(17.26)
Фиг. 17.1. Поток тепла на границе элемента.
Подстановка (17.26) в (17.22) дает
j [W]r
-fj- d°z=h J [Л^]7' [W] {Ф} <№— j IN]T
(17.27)
26
26
26
где N =[0 L2 L3]. В ы п о л н я я
интегрирование в (17.27) с
помощью
плоских L-координат, приходим к результатам, идентичным тем,
что получены в гл. 8.
Использование метода
Галёркина непосредственно
приводит
к слагаемым,
которые в
вариационной формулировке
должны
быть добавлены к функционалу, чтобы учесть граничные условия. Метод Галёркина применяется также при решении двумерных и трехмерных задач теории упругости [5]. В результате получается система уравнений, подобная той, которая соответствует вариаци онной формулировке этих задач.
17.4. Задача Коши
Метод Галёркина может быть использован при решении зада чи Коши, а также при решении переходной задачи, обсуждавшей ся в гл. 11. В этом разделе будет рассмотрена задача Коши для одного дифференциального уравнения, а затем проведено обобще ние на случай системы дифференциальных уравнений первого по рядка.
Рассмотрим дифференциальное уравнение
^ г + 4 у = 0
(17.28)
с начальными условиями у(0 )= 0 и у '(0 )= 4 . Это
уравнение име
ет очевидное решение
y = 4s\n2t.
Приближенное решение этого уравнения можно получить числен но. Чтобы проиллюстрировать применение метода Галёркина, ис пользуем именно этот способ.
Подставляя (17.28) в (17.1), получаем
t
j W (чрг + *у) dt= ° -
07.29)
о
Мы уже указывали на необходимость разбиения интеграла на сумму и преобразования интеграла, содержащего производную высокого порядка, в интеграл от первой производной. Подобное преобразование уже рассматривалось по отношению к dzy/dxz. Учитывая формулу (17.9), получаем
W ! ? - < # = №
- J
1 - « ,
где Те— шаг по времени (длина) отдельного элемента. Уравнение (17.29) может быть записано теперь в виде
7 -4 " /=0
£
%
- 4 [ № ‘Г
y ) d t = 0 , (17.30)
e=i
rt
где Ti
и Tj — значения времени, соответствующие узлам i и / эле
мента, a R — число элементов.
Применение формулы (17.30) будет проиллюстрировано с по
мощью
линейного интерполяционного
полинома
для у:
Соотношение (17.31) определено относительно местной системы координат с началом в i-uузле, что соответствует пределам ин тегрирования Ti = 0 и Tj = Te. Подстановка выражения (17.31) в (17.30) дает
п е
d [N (e)] T
d [N (e)]
- т т- % , - . - 2 К
[ У ) ~
dt
dt
e= l 0
— 4 [7V<e>]r [ЛД£>] {F }jd /=00.
(17.32)
Для первого
элемента получаем следующие уравнения:
■{о)
т [ _ 1
— 1
YA + ^
: у в - а
(17.33)
m
v j
1
6
Любому другому элементу соответствуют уравнения
— 1
{
^
Н
2
1 1 0 = 1 “ !.
07.34,
М - 5
1
:
1ш
\
так как во всех других узлах угол поворота не определен. Объ единяя уравнения для отдельных элементов и предполагая длину элементов одинаковой, получаем
—4'
1
— 1
Уг
0
— 1
2
— 1
У,
0
1
— 1
2
— 1
У3 +
г.
“
— 1
2 — 1
У4
2
1
Ух]
0 '
1
4
1
У2
0
4Т е
1
4
1
у3.
= . 0
(17.35)
1
4 1
у*
0
Заметим, что все уравнения, кроме первого, идентичны. Система уравнений (17.35) может быть записана в виде
-
4 - ^ - ( П
- П
)
+ 4 1
(2K1 + F2) = 0 ,
(17.36а)
( - У п.г +
2Уп- У п+г) +
- [
-
(Уп-1 +
4 К Л + Г п+1) = 0 ,
п ^ 2 , ( 1 7 . 3 6 6 )
где п — произвольный узел.
Зная Тб1 из этих двух уравнений мож
но определить все значения
{У}.
Пример
171. Требуется решить дифференциальное уравнение (17.28), считая шаг по времени равным 1/16 с. Начальные условия сле дующие: 1/(0) = 0 и у'(0 )= 4 . Сравнить результаты расчетов с ана литическим решением.
Рекуррентные соотношения, соответствующие рассматриваемо му уравнению, даны в (17.36). Заметим, что условие г/(0) = 0 = позволяет решить уравнение (17.36а) относительно У2:
_ 4 _ 1 6 ( - Г 2) + А ( - 1 - ) Г 2= 0 ,
ИЛИ
Уравнение (17.366) может быть решено относительно Уп+ь Полу чаем соотношение
у
, — 764
. у
у
,
1 п+1
385
1 п
1 п—1»
откуда находим
764
„
764
96
-о,
385
2 - У х — 385
*
385
У 3 = 0 , 4 9 4 8 .
Повторяя эту процедуру для каждого временного шага, вычисляем последовательно все узловые значения. Ниже приведены значе ния У, полученные численным методом, а также значения, соот ветствующие точному решению.
t
Метод конечных
Точное
t
Метод конечных
Точное
элементов
решение
элементов
решение
0
0 , 0
0 , 0
V ie
1,8051
1,8045
V ie
0,2494
0,2493
10/ю
1,8987
1,8980
2/ l 0
0,4948
0,4948
n /io
1,9627
1 , 9 6 1 8
3/ ю
0,7326
0,7325
12/io
1,9962
1,9950
V ie
0,9589
0,9589
I3/l0
1,9985
1,9971
V u
1,1703
1,1702
14/lO
1,9697
1,9680
o/io
1,3635
1,3633
15/l0
1,9101
1,9082
7 i o
1,5354
1,5351
‘ V io
1,8209
1,8186
7 i „
1,6833
1,6829
В случае использования метода Галёркина при решении задачи Коши получаются уравнения с двумя замечательными особенно стями. Шаг по времени может изменяться, если в этом есть не-
обходимость; могут варьироваться и функции формы, входящие в [ЛДе)]. В случае большой величины шага по времени можно ис пользовать элементы высокого порядка.
Изменение шага по времени вызовет модификацию системы уравнений (17.35). Эта модификация будет выражаться в появле нии более одной пары рекуррентных соотношений типа (17.36). Некоторые из этих соотношений будут включать как новые, так и старые приращения времени.
Если вместо линейного интерполирования (17.31) применить функции формы для квадратичного элемента, вместо двух будут получены три уравнения. Первые два уравнения используются для определения У2 и Уз. Третье соотношение рекуррентное, оно выра жает последовательно одно из узловых значений через три пре дыдущих:
Такая ситуация всегда возникает при решении задачи Коши с по мощью метода Галёркина; всегда имеется достаточное число урав нений, чтобы можно было вычислить значения У, требуемые для проведения вычислений по рекуррентным формулам.
17.5. Система дифференциальных уравнений первого порядка
Система дифференциальных уравнений первого порядка вида
[ С ] 4 { ф } + [^П ф ) + И = о
(17.37)
обсуждалась в гл. 11, где было дано конечно-разностное решение этих уравнений. Теперь решим эту систему методом Галёркина. В результате для вычисления значений {Ф} получим матричное уравнение, которое несколько отличается от уравнения (11.23).
Используя для интегрирования неизвестной величины ф линей ную модель, можно записать
Фр) =Ы^Фи +
ф р = № ? Ф и +№ *Ф 2},
(17.38)
Ф<«)=ЛГМФ„+ т е)фг.г
где индексы i и / используются для обозначения двух узловых зна чений, разделенных во времени на величину шага длины Те, а г обозначает общее число узлов. Соотношение (17.38) в матричном виде имеет вид
Здесь {Ф }— матрица размера n X l, а ЛД*> и N<г)— функции фор
мы. Подстановка выражения (17.39) в (17.1) дает два матричных уравнения:
те
+[7С1{Ф) 4{ F ) y t = 0
Jл Ц [ С ] - ^ }
(17.40а)
* е
J Nj ([С] - 4 М
+ 1К\ {Ф} 4- [F)) d t = 0.
(17.406)
Дифференцируя по времени соотношение (17.39), получаем
{Ф}. _
(ф).
(фЬ.
(17.41)
Подставляя (17.41) и (17.38) в (17.40а) и выполняя интегриро вание, имеем
Последние два уравнения могут быть объединены в одно:
(-4 ic ]+ -^ !*j)(4 -[c i+ -£ -ro ) (I®},)
1|(f|,|
( — g - L C l H —
[ / С ] ) ^ 4 _ [ С ] + - ^
[ / С ] )
2 ( И Л
Это уравнение для элемента,
соответствующее одному шагу
по времени. Для получения системы уравнений, определяющей значения {Ф}ь {ФЬ и {Ф}з, ...» оно должно быть объединено с аналогичными уравнениями для соседних временных шагов. Пусть всем элементам соответствуют одни и те же приращения времени; объединяя уравнения для первых двух шагов по времени, по-
лучим
' ( - - r t c i + ^ W )
( - - 5-1С]+тЧК|) (4 " №)
( - - J-lC]+-£-[Kl)
о
( 4 - i c i + ^ m ) X
( 4 - [ С ] + - ^ ( Л ] )
[(Ф)х1
4 - 1^ь
(17.45)
X . {ф}2
{F\t
(Ф}з
.
Шз
Так как {Ф} 1 известно, то первое уравнение в (17.45) можно ис
пользовать для вычисления {Ф}г- Все остальные уравнения, начи ная со второго, идентичны между собой и могут быть записаны в общей форме:
( - - Y [С] + [ffJ) {Ф}«-1+ - ^ [ К ] [ Ф ) п +
+ (4 1 с1 + ^ [/(1)1фи = -
^
п > 2-
(17.46)
Это соотношение позволяет вычислить все требуемые
значения
{Ф Ь для любого п больше двух.
иметь
в виду
При использовании формулы (17.46) следует
две особенности, которые отличают ее от результатов, полученных
разностным методом. Для
вычисления {Ф} п+1 необходимо знать
два
вектор-столбца
{Ф} п- 1
и {Ф}п, а кроме того, требуется пом
нить
три матрицы
размера
[С]. Последнее требование приводит
к значительному загружению машинной памяти и представляет определенный недостаток при решении систем большого порядка.
Для реализации метода Галёркина во временной области мож но применить другую процедуру: рассматривая шаг по времени как отдельный элемент, вычислить {Ф};, используя {Ф}*. При этом вектор-столбец {Ф}; может быть найден либо из уравнения
(17.42), либо из
(17.43). Уравнение (17.43), видимо, более широ
ко применяется,
хотя (17.42) будет давать для {Ф};- значения,
почти идентичные тем, что
получаются с помощью (17.43).
Как показано
в работе
[1], использование уравнения (17.43)
дает для {Ф}; значения, которые меньше осциллируют, чем в слу чае решения исходной системы (17.37) разностным методом, рас смотренным в гл. 11.
Вектор {F} в (17.37) может изменяться со временем. В этом случае его нужно вычислять для каждого временного шага и уравнения (17.42) и (17.43) должны быть видоизменены.
17.6. Заключение
Сочетание метода Галёркина с кусочной аппроксимацией ме тода конечных элементов является чрезвычайно эффективным спо собом решения многих дифференциальных уравнений. Несколько примеров, связанных с техническими расчетами, были обсуждены в этой главе. Метод Галёркина, безусловно, получит широкое рас пространение благодаря тому, что он позволяет обходиться без вариационной формулировки задачи.
Задачи
172. При осевом нагружении некоторого элемента конструкции его перемещения описываются следующим дифференциальным уравнением:
du
Р
О,
dx
АЕ
где и — перемещение, см;
Р — осевая
нагрузка, Я; А — площадь
сечения, см2; Е — модуль
упругости,
Н/см2. Получите матрицы
элемента для этого дифференциального уравнения, считая величи ну Р/АЕ постоянной.
173. С помощью матриц элемента, полученных в задаче 172, вычислите узловые перемещения для детали конструкции, изобра женной ниже. Используйте два элемента равной длины.
I
S000H
---------
120см
—
>1
Е=2*107 Н/см А=3см2
Кзадаче 173.
174.Определите матрицы элемента для дифференциального уравнения в задаче 172, если величина Р/АЕ меняется линейно от
одного узла до другого. Используйте результаты решения зада чи 108, рассмотренной в гл. 12.
175. Выведите уравнения для элемента, соответствующие диф ференциальному уравнению (17.28). Используйте элемент с тремя узлами (квадратичный). С помощью выведенных уравнений реши те исходную задачу Коши для интервала O ^ f ^ l .
22*
176. Используя линейный интерполяционный полином, получите уравнения для элемента, соответствующие дифференциальному уравнению первого порядка:
Н-яу + £ = 0 ,
где а и b — заданные константы, а у(0) известно.
177. Используя результаты, полученные в задаче 176, решите одно из следующих дифференциальных уравнений. Сравните ре
зультаты с теми, что дает аналитическое решение.
а)
У' + З у = 0 ,
#(0) = 4 ,
1;
б) У' — 4 0 + 2 = 0 , 0 ( 0 ) = 1 , 0 <
1 ;
в) у ' + 2 у — 1 = 0 ,
0(0)= 3,
0 < г < 1 / 2 ;
г ) у '— 6 у — 6 = 0 , 0 ( 0 ) = 0 ,
1 .
178. Решите
задачу 176,
применяя
следующие элементы:
а) квадратичный интерполяционный полином; б) кубичный интер поляционный полином. Решите одно из уравнений в задаче 177.
179. Выведите уравнения для элемента, соответствующие диф ференциальному уравнению второго порядка:
- Ж + а - Ж + ь У + с = ° >
где а, b и с — заданные константы. Для у используйте линейный интерполяционный полином.
180. Используя результаты, полученные в задаче 179, решите
следующие уравнения:
а)
у"—Зу' + 4г/=0, у(0)=у’(0) = 1, 0 < / <
1/2;
б)
у " - 2 у + 1 = 0 ,
у{0 )= 2 ,
г/'(0)=0,
0 < / <
1/2;
в)
у" + у ' ~ 6 = 0 ,
у(0)=0,
у'(0 )= 3 ,
0 < / <
1/2.
ЛИТЕРАТУРА
1.Donea J., On the Accuracy of Finite Element Solutions to the Transient HeatConduction Equation, Intern. J. for Numerical Methods in Engineering, 8 , 103— 110 (1974).
2.Martin H. C., Carey G. F.f Introduction to Finite Element Analysis, McGraw-Hill, N. Y., 1973.
3.Norrie D. H., deVries G., The Finite Element Method, Academic Press, N. Y., 1973.
4.Sokonikoff I. S., Mathematical Theory of Elasticity, 2-nd ed., McGraw-Hill, N. Y., 1956.
5.Szabo B. A., Lee G. C., Derivation of Stiffness Matrices for Problems in Plane Elasticity by Galerkin’s Method, Intern. J. for Numerical Methods in Engineering,
1, 301—310 (1969).
6 . Zienkiewicz О. C., The Finite Element Method in Engineering Science, McGrawHill, London, 1971; есть русский перевод: Зенкевич О., Метод конечных элемен тов в технике, изд-во «Мир», М., 1975.
УЧЕБНЫЕ ПРОГРАММЫ
Во всех предыдущих главах подчеркивалась необходимость машинной реализации метода конечных элементов. Очевидно, что метод конечных элементов не пригоден для проведения расчетов вручную. В этой главе будут рассмотрены некоторые программы, которыми следует пользоваться при изучении материала, пред ставленного в гл. 2, 6 и 8— 12. Настоящая глава не должна рас сматриваться как последняя глава этой книги. Изложенный здесь материал нужно использовать в качестве приложения при обсуж дении конкретных применений метода конечных элементов.
Программы, приведенные в этой главе, далеки от тех сложных программ, которые могут решать самые различные задачи. Они не предназначены для того, чтобы конкурировать с имеющимися стандартными программами. Эти программы преследуют только учебные цели. Такие несложные программы весьма желательны с учебной точки зрения, так как они сокращают до минимума класс ное время, требуемое для объяснения ее работы. Перечислим не которые характерные особенности, которые приводят к упрощени ям в программе:
1.Используются элементы только одного типа — линейные треугольники.
2.При разбиении области на элементы характеристики мате риала каждого элемента предполагаются одинаковыми.
3. Главные оси инерции параллельны координатным осям
х, у.
4. Каждая программа решает только однотипные задачи, на пример задачу о кручении упругого стержня, двумерный случай переноса тепла, двумерные течения жидкости или двумерные за дачи теории упругости.
5.Редко встречающиеся варианты счета исключаются.
6.Программы мало отличаются между собой в отношении тре буемых исходных данных и организации их ввода.
Даже в том случае, когда разнообразие вариантов счета в программе сведено к минимуму, методика программирования включает такие характеристики, которые делают эффективным ис пользование этих вариантов. Все коэффициенты системы уравне ний [К]!{Ф} = {Г} хранятся в машинной памяти в виде отдельного вектор-столбца. Такой способ исключает необходимость заранее проставлять размеры отдельных компонент [К], {Ф} и {F} и зна
