CC BY
183
УДК 624.04:519.6
А.В. Игнатьев
ФГБОУВПО «ВолгГАСУ»
ОСНОВНЫЕ ФОРМУЛИРОВКИ МЕТОДА КОНЕЧНЫХ ЭЛЕМЕНТОВ В ЗАДАЧАХ СТРОИТЕЛЬНОЙ МЕХАНИКИ
Часть 2*
Предложена классификация формулировок метода конечных элементов (МКЭ), позволяющая ориентироваться в огромном количестве опубликованных и продолжающих публиковаться работ по проблеме повышения эффективности этого самого распространенного численного метода. Во второй части статьи продолжено рассмотрение прямых формулировок МКЭ в форме метода перемещений, метода сил и в форме классического смешанного метода. На примере одномерной системы-балки рассмотрен вопрос о сходимости решения по МКЭ в форме классического смешанного метода при сгущении конечно-элементной сетки.
Ключевые слова: метод конечных элементов, классификация формулировок МКЭ, строительная механика.
Рассмотрим прямоугольный конечный элемент (КЭ) er в виде изгибаемой прямоугольной пластинки с 12-ю степенями свободы (КЭ Клафа) (рис. 1).
В классической теории пластин для упрощения задачи делаются предположения о линейном законе изменения деформаций и напряжений по нормали к плоскости пластины (гипотезы Кирхгофа — Лява).
Р-
У
уН
у
0,М.
-е-
Рис. 1. Изгибаемая прямоугольная пластинка с 12-ю степенями свободы (КЭ Клафа)
При принятых допущениях в каждом узле элемента имеем по три неизвестных перемещения: прогиб срединной плоскости и два угла поворота. Положительные направления осей координат, перемещений и сил показаны на рис. 6. Обозначим V перемещение вдоль оси г; 0х и 0у — углы поворота относительно осей х и у соответственно; р^ — поперечная нагрузка; М и Мх — соответствующие внешние моменты.
Вектор узловых перемещений для рассматриваемого элемента имеет вид
Ы =[q,, , qk, q ] r
где [q ] =
w
{') a{t) a{t)'
x ' y
, {t = i, j, k, l).
(1)
Таким образом, как и матрица жесткости, он имеет 12-й порядок.
Часть 1 опубликована в журнале «Вестник МГСУ», 2014, № 11.
Функции, аппроксимирующие линейные и угловые перемещения точек серединной поверхности, в соответствии с количеством степеней свободы узловых перемещений элемента также должны содержать 12 констант. Особенность здесь заключается в том, что в силу принятого допущения о линейном изменении деформаций и напряжений по нормали к плоскости пластины углы поворота 9x и 9y не являются независимыми величинами и определяются известными соотношениями
0w 0w
9y =—, 9Х = —. (2)
ox oy
Это означает, что деформированное состояние пластины полностью описывается одной величиной — прогибом w срединной поверхности пластины. Следовательно, функция, аппроксимирующая поверхность прогибов элемента (функция формы), должна содержать 12 констант.
Для изотропной пластинки, например, в качестве функции формы может быть принят неполный бикубический полином
w (x, y) = a1 + a2 x + a3 y + a4 x2 +a5 y2 + a6 xy +
2 2 3 3 3 3 (3)
+a7xy + a8xy +a9x +a10y +a11 xy + a12xy , удовлетворяющий однородному дифференциальному уравнению ее чистого изгиба
04w + 2 04w + 04w 0 (4)
dx4 0x 2dy2 0y4
так как рассматривается только узловая нагрузка на элемент.
Тогда в соответствии с (2)
9y =a2 + 2a4 x + a6 y + 2a7 xy + a8 y2 + 3a9 x2 + 3a11 x2 y + a12 y3;
9x =a3 +2a5 y + a6 x + a7 x2 +2a8 xy + 3a10 y2 +an x3 +3a12 xy2.
Полином (3) удобен во многих отношениях. В частности, он удовлетворяет, как отмечено выше, уравнению равновесия и в то же время условиям неразрывности деформаций по прогибам на границе двух элементов. Полагая в (3), например, x = const, получаем полином третьей степени относительно у, для определения четырех произвольных постоянных которого используются четыре общие для смежных элементов константы: два угла поворота относительно оси x и два прогиба в узловых точках на концах границы. Поскольку коэффициенты полинома единственным образом определяются этими четырьмя константами, то функция перемещений w (и, как следствие этого, функция 9x) вдоль всей границы между элементами является единственной и общей для смежных элементов.
Хуже обстоит дело с углами поворота по нормали — 9 y = .
При x = const для 9y получаем полином третьей степени относительно y, для определения четырех произвольных постоянных которого имеется только две константы: два угла поворота 9y на концах границы между элементами. Следствием этого является неоднозначность решения для коэффициентов полинома. Поэтому в общем случае угол наклона 9y для границы x = const может оказаться разрывным.
(5)
Таким образом, функция формы (3) является несогласованной. Погрешность в определении углов 9 на линиях х = const и 9x на линиях у = const является неустранимой при любом разбиении рассматриваемой области. КЭ с несогласованными функциями формы (КЭ Клафа) называются несовместными. Выполненные теоретические исследования показали, что использование КЭ Клафа обеспечивает сходимость решения задачи, и поэтому его широко используют в расчетах благодаря физической наглядности его степеней свободы. С помощью согласованных функций, например, рассмотренных в [1—3], можно получить более корректное, но, как правило, менее точное решение. Поэтому для практического применения рекомендуются функции типа (3), т.е. удовлетворяющие условию чистого изгиба (4).
Выражения (3) и (5) могут быть представлены для элемента е как и для стержневого КЭ, в виде
(6)
{U (x, y )} = [ Ф (x, y )]{а}
где
« (x,y) " Ф1 (x, y)
{U (x, У )} «2 (x. У) ■; [Ф(x,y)] = Ф2 (x, У)
«3 (x,У) _ Ф3 (x, У)_
\ (x,У) = w(x,y) = [ф (x,y)]{а};
д
'h (x y) =—[w (x y )] = [ф2 (x у )]{а} = O;
д
ъ (x y) = dy[w (x' у )] = [фз (x, у )]{а} = ox; } = [
(7)
¡а
а
а12 ] >
[ф (x, У )] =
Ф1 (x, У ) Ф2 (x, У)
Ф3 (x,y)
Так как вектор узловых перемещений {q} для элемента ег может быть получен из (6) путем подстановки координат узлов, то связь между векторами {q} и {а} определяется выражением
{q} = V ]{«}, (8)
~ Ф (' y iГ Ф (('yj)
Ф (, Ук)
_Ф (, yt )_
Решая уравнение (8) относительно {а.} и подставляя полученный результат в (6), находим аналогично (9)—(13) в [4].
U (x,y) = w (x,y) = [£/]{q}, (9)
где [U ] = [ Ф! (x, y )][V ]-1. (10)
в котором [V ] =
Проектирование и конструирование строительных систем. Проблемы механики в строительстве УЕБТЫНС
_мвви
По аппроксимирующей матрице перемещений найдем вектор {%}, характеризующий кривизны срединной поверхности элемента:
д? [ - (х у )] Ь (Х,У )]
2 ~дЬу [- (х,уЯ
С учетом (7) в [4] это выражение может быть представлено в виде (13) [4], т.е.
М={№ М (12)
М =
(11)
где [B] = {A}[U], {A} = -
д2 __д!
-2 ' д..2 '
— матрица дифференцирования.
дх Су дхду
Для осуществления перехода от вектора деформаций (кривизн) к вектору усилий (моментов) {Щ введем к рассмотрению матрицу упругих свойств материала [Щ, которая для ортотропного материала имеет вид
[D ] =
ЦА D
y y 0
0 0
0 2D
(13)
где Бх = ^/12 (1 -ц Ц ), Ву = Е2И3/12 (1 ), Бк = ЦИ3/б, Бх ц = Бу ц = = Бц; Е1, Е2 — соответствующие модули упругости; цх, ц — коэффициенты Пуассона.
Для изотропного материала
1 ц 0
ЕИ
(14)
[D ] =
ц 1 0 0 0 1 -ц
12 (1 -Ц2 )'
Используя теперь известную связь между деформациями (кривизнами) и усилиями (моментами) в пластинке, получим
{N } =
Mx
My
= [D ]{Х},
(15)
где Мх = | схгс1г — погонный изгибающий момент, который представляет
- и 2
собой интегральную характеристику напряжений в направлении оси х;
и 2
Му = | с— погонный изгибающий момент, т.е. интегральная характери-
2
стика напряжений в направлении оси у; Мху = | т— погонный крутящий
и 2
момент, т.е. интегральная характеристика касательных напряжений.
Процесс построения матрицы жесткости остается таким же, как и стержневого КЭ, т.е. стандартным. Поэтому в окончательном результате (23) [4] следует лишь интегрирование по длине заменить интегрированием по площади рассматриваемого элемента е .
и и
[К] = {{[Б]Т [Б][Б]хёу.
(16)
Приведение нагрузки, действующей на элемент, к узловой выполняется так же, как и в одномерном случае
(17)
При равномерно распределенной нагрузке вектор эквивалентных сил, прикладываемых в узлах элемента, имеет следующий вид:
{/У |_ /' /}' /' /1
= 4 раЬ
(18)
1 Ь а 1 Ь а 1 Ь а 1 Ь а _4' 12'12'4'12'12'4' 12' 12'4'12' 12 Рассмотрим теперь прямоугольный КЭ с 16-ю степенями свободы (КЭ Богнера — Фокса — Шмита) [5].
Для такого элемента в каждом узле принимается по четыре степени свободы. Дополнительно к трем степеням свободы узла КЭ Клафа учитывается деформация закручивания в узле. Вектор узловых перемещений для рассматриваемого КЭ (рис. 2) с номером г имеет вид
ы=
где {ъ I =
Я.
ул
ху ,1
( = 1,],к,1 или 11, 21 ' 22, 12).
Рис. 2. Вектор узловых перемещений для КЭ с номером г
ь
Функция формы для этого КЭ записывается в явном виде
^(х, у)=у у Ги> Э(1 > (х, у)+0(1)Э(2) (х,у) + е(уЬ(3} (х, у)+
\ ' у ) /,/ А зп зп \ > у / зп зп\ ) зп зп \ 'У /
з-1 п-1 (19)
+00у )э(4](х, у)
где Э(1))(х,у) и ™тп — составленная на основе произведений полиномов Эрмита (26) [4] координатная функция (функция формы) и обобщенная координата, соответствующие вертикальному перемещению в узле тп (рис. 3); Э32 (х, у) и е« — те же, соответствующие угловому перемещению относись Э(3) / \ а(у) С™
тельно оси х, т.е. —; Эзп (х,у) и етп — те же, относительно оси у, т.е. —;
ду дх '
/ ч п(х, у) д 2 ™
Э3п (х,у) и ет„ — те же, соответствующие закручиванию в узле тп, т.е.-.
дхду
Рис. 3. Координатная функция (функция формы) и обобщенная координата, соответствующие вертикальному перемещению в узле тп
p« (x,у) = Э0>и (x)Э0;Я (у); Э(2> (x, y) = 31>m (x)30>„ (у);
{Э(3] (x,у)(x,у) = Э^ (x)Э1,„ (у); Э!4] (x,у) = Э^ (x)Э^ (у);
{Э0,т-1 (x) = "7 (2x2 - 3ax2 + a ); Э0,и-2 (x) = (2x2 -3ax2);
a 4 y a x '
{Э1,т-1 (x) = Л (( - 2axl + a2x) ; Э1,т-2 (x) = Л (( - ax2 ) ■
a 4 y a x '
Полиномы Э0п1(у), Э0п2(у), Э1п1(у), Э1п2(у), получаются из соответствующих полиномов Эп(х) путем замены в них a на b и х на у.
Особенностью этих функций форм, определившей их выбор, является то, что они обеспечивают полную совместность деформаций со смежными КЭ и имеют нулевые значения во всех узлах, кроме узла и номера связи (степени свободы), к которой они относятся.
В этом нетрудно убедиться если, например, в функцию f = Э01(х) Э01(у) подставить координаты узлов, то получим
f = 1 f = 0 f = 0 f = 0
•/1,11 '-'1,21 '-'1,22 '-'1,12
Так как для рассматриваемого типа КЭ все степени свободы имеют физический смысл перемещений, то элемент K.. его матрицы жесткости — это усилие (реакция), возникающее в связи по направлению i-й степени свободы
'м; ХХ
М={^ н Му , {е}={х} = < % у
% ху
вестник 12/2014
от единичного перемещения связи по направлению у'-й степени свободы при условии, что перемещения связей по направлениям всех остальных степеней свободы равны нулю (д. = 0, 1 фу).
На основании принципа возможных перемещений
Ку = {с у е^П, (20)
п
где с . — напряжения по области П КЭ от перемещений д . = 1; е . — деформации по всей области КЭ от перемещений д . = 1.
При расчете пластин в качестве физических соотношений принимаются зависимости между интегральными характеристиками напряжений (погонными изгибающими и крутящим моментом) и кривизнами срединной поверхности:
{с} = те}, (21)
где [Б] — матрица упругих свойств (13),
(22)
вектор усилий (напряжений) по области КЭ, составляющими которого являются погонные изгибающие и крутящие моменты и вектор деформаций по области срединной поверхности КЭ, составляющими которого являются кривизны (5), (9).
С учетом (21) и (22) выражение (20) может быть записано для всех коэф -фициентов матрицы жесткости следующим образом:
[К Н[хГ [Б ][х]<Ю, (23)
п
где [х] — матрица кривизны (12); [Б] — матрица упругих характеристик КЭ (13).
Выполнение в [6] с использованием этого типа КЭ расчеты пластинок и оболочек в [7] показали, что при одинаковых сетках использование КЭ Богнера — Фокса — Шмита обеспечивает лучшую сходимость к точному решению по сравнению с КЭ Клафа.
Из приведенных примеров видно, что в общем случае (даже в случае пластинки) трудно подобрать функцию, аппроксимирующую поле перемещений и обеспечить тем самым необходимые условия для сходимости получаемого решения к теоретически точному. Это наиболее четко проявилось в [7], где в качестве узловых неизвестных КЭ тонкостенной оболочки приняты перемещения, их первые и смешанные производные.
Для преодоления проблем, связанных со сходимостью, разработаны и разрабатываются другие (не прямые) подходы, основанные на вариационных принципах с менее жесткими требованиями к аппроксимирующим функциям [3, 8, 9].
При использовании МКЭ в перемещениях матрица жесткости КЭ получается путем минимизации функционала полной потенциальной энергии системы. В вариационной постановке достаточным условием сходимости МКЭ являются критерии полноты и слабой согласованности [10].
В расчетах сложных конструкций часто используются, особенно при локальном сгущении, так называемые «несогласованные сетки» с независимой аппроксимацией перемещений или усилий в области КЭ [3, 11].
Для выполнения условий совместности на межэлементных границах используются различные зачастую сложные в реализации алгоритмы: переходные элементы с промежуточными узлами на сторонах и более сложными функциями формы; построение элементов с переменным количеством узлов или введение специальных промежуточных элементов.
Условия непрерывности перемещений или усилий (напряжений) на межэлементных границах принудительно удовлетворяются с помощью множителей Лагранжа или штрафных методов [12, 13].
МКЭ в форме метода сил. Формулировка метода сил аналогична формулировке метода перемещений [14]. В качестве неизвестных в основной системе метода сил для КЭ принимаются статические силы или напряжения. Поле сил (напряжений) внутри области КЭ, аппроксимируется с помощью интерполяционных функций. Предполагается, что должны быть выполнены условия равновесия каждого из КЭ, а также в узлах и на границах между элементами.
Эти условия стали непреодолимым препятствием на пути формулировки МКЭ в форме классического метода сил. Проблема состоит в том, что матрица податливости КЭ может быть получена только для двух случаев. В первом из них КЭ должен быть статически определимой системой, во втором — должен находиться под воздействием самоуравновешенной системы сил.
Ни по одному из этих вариантов получить уравнения, связывающие матрицу податливости КЭ с узловой нагрузкой невозможно.
Поэтому существующие варианты МКЭ с узловыми неизвестными в виде напряжений (усилий) основаны на модифицированном функционале дополнительной энергии.
Имеются несколько таких вариантов.
1. Решение в напряжениях строится на основе функционала дополнительной энергии всей конструкции в целом.
Для аппроксимации напряжений по области КЭ используются постоянные функции, что позволяет уменьшить количество неизвестных параметров и соответственно уменьшить вычислительные затраты [14].
Дополнительная энергия всей конструкции получается суммированием энергий всех КЭ
где [о] — глобальный вектор неизвестных для всей конструкции, состоящий из векторов неизвестных напряжений КЭ; [5] — блочно-диагональная матрица податливости конструкции, состоящая из блоков-матриц податливо-
Дополнительно (на основе принципа возможных перемещений) составляется система алгебраических уравнений равновесия для узлов конечно-элементной сетки конструкции:
(24)
стей КЭ.
К Г Ь} + ри = 0,0' = х, у, г), (25)
где {С. ' Г — матрицы податливости; {с {| — вектор неизвестных, состоящий
из напряжений в КЭ, примыкающих к узлу 7; {р . | — свободные (грузовые) коэффициенты, состоящие из узловых и объемных сил КЭ, примыкающих к
узлу 7.
Для получения решения задачи находится минимум функционала (24) при наличии ограничений в виде системы линейных уравнений равновесия (25):
П* = ±{а} [В]{а} + УУа({Сг,,} (а,} + Р1,1 | ^ шт.
(26)
Для минимизации квадратичного функционала (64) используется метод штрафных функций.
Необходимо отметить, что при измельчении сетки, получаемые по найденным напряжениям перемещения, сходятся к точным значениям сверху. Это означает, что они могут быть использованы для оценки точности получаемых решений вместе с решениями по методу перемещений.
2. Если для решения задачи использовать метод множителей Лагранжа, то получим функционал
п*=2 (} [в ]{аМ1«> ы
2 ! 1
С, 1} {} + ри
^тш,
(27)
где п.. — множители Лагранжа, имеющие физический смысл действительного
V
перемещения узла 7 по направлению.'.
Для минимизации функционала (27) приравниваются к нулю его частные производные по неизвестным напряжениям (усилиям) и перемещениям. В результате получается система уравнений
[В](а} + [СГ М = 0; [С ]{а} + [Р ] = 0,
(28)
где {п} — вектор неизвестных узловых перемещений.
Эту модификацию следует рассматривать уже не как метод сил, а как смешанную формулировку МКЭ, так как в ней принимаются за неизвестные как напряжения (усилия), так и перемещения.
При представлении конструкции объемными КЭ ее применение оказывается эффективным [14].
3. Более приемлемым в смысле возможностей применения является вариант, когда за основные неизвестные принимаются функция напряжений и ее производные [14].
Для простейшего КЭ прямоугольной формы с узлами только в вершинах (рис. 4) за основные неизвестные в узлах элемента принимаются соответствующие значения функции напряжений Ф и ее производных Ф, , Ф, , Ф, .
1 г-^ х у ху
Рис. 4. Простейший КЭ прямоугольной формы с узлами только в вершинах
По аналогии с МКЭ в перемещениях получаем КЭ с 16 степенями свободы. Вектор основных неизвестных в узлах элемента е можно представить в виде
где
Ф} = [< ! Ф 2 Ф 3 Ф 4 ], (29)
(<: } = [ф(;) <(? Ф(ХУ ], I = 1,2,3,4. (30)
Так как напряжения ох, оу, тху являются производными от функции напряжений Эри, т.е.
МГ =[Ф,УУ <Ф» ФхУ ], (31)
то для обеспечения постоянного поля напряжений в элементе и совместности на межэлементных границах функция напряжений Ф(х, у) в элементе аппроксимируется, как и в методе перемещений полиномами Эрмита:
< (х, У ) = Э (х, У), (32)
где (х,у) = Э (х(у). (33)
Выражениями (32) и (33) определено 16 интерполяционных функций, устанавливающих связь между функцией напряжений Ф е = Ф (х, у) и ее значениями в узлах элемента е.
Подставив в (31) выражение (32), получим
(а)Т = [Э,УУ Э,хх Эху ] = [Н ]е Ф• (34)
Дополнительная энергия деформации элемента е будет определяться выражением
и* = 2XI |К^Х^ (35)
2 '=1 ]=1 ое
Т V ~Т
где {в} =|^вх ву уJ — вектор относительных деформаций. С учетом (34) выражение (35) может быть записано в виде
и*=2 к к (ф^Ф J), (36)
2 :=1 J=1
где Бу е = | [итв НХ^О — податливость элемента на воздействие единичного
Ое
неизвестного Ф :1.
V
Система разрешающих уравнений, с помощью которых описывается напряженно-деформированное состояние всей конструкции, как системы КЭ, получается на основе уравнений для отдельных КЭ и дополнительных уравнений в виде силовых условий на контурах элементов. В этих уравнениях неизвестные — значения функции напряжений и ее производных в узлах
ЕФе - ре = 0, (37)
где Е — матрица коэффициентов при неизвестных; ре — вектор свободных членов, зависящих только от нагрузки.
Для всей конструкции в целом получаем
ЕФ - Р = 0. (38)
На основе выражений (30) и (31) формируется модифицированный функционал дополнительной энергии
иМод =ФТВФ -Хр + ФТЕТX, (39)
где В — матрица податливости всей конструкции; X — вектор неопределенных множителей Лагранжа (дополнительных неизвестных) с помощью которых вводятся условия — ограничения.
Примеры использования этого варианта МКЭ можно найти в [14].
МКЭ в форме смешанного метода. Построение матрицы откликов и вектора воздействий для КЭ — стержня. Рассмотрим изгибаемый КЭ — стержень плоской рамы (рис. 5, а).
Рис. 5. Изгибаемый КЭ
В основной системе смешанного метода для этого КЭ за неизвестные принимаются одновременно как усилия, так и перемещения [15, 16]. Самый эффективный вариант основной системы показан на рис. 5, б. В ней примем следующие обозначения неизвестных: д2, д3, д5, д6 — перемещения по направлению реакций Я1, Я2, Л3, Л4, Л5, Л6, во введенных связях; q1, — неизвестные усилия Q, М, N в удаленных связях.
Ввиду простоты основной системы элементы матрицы откликов для рассматриваемого КЭ получаются непосредственно из физического смысла задачи, т.е. из условий неразрывности, статических и кинематических условий.
Матрицу откликов КЭ представим в блочном виде г Г
(40)
[а ] =
5 5
Для определения коэффициентов, входящих в блок 5, матрицы откликов (40) используются известные формулы для определения перемещений:
=
|ММ„
Е1
ах
0
ЕЕ
■ж
ОЕ
ах.
0 0 0 Выполняя перемножение соответствующих эпюр усилий от единичных силовых воздействий ql = 1, q8 = 1, q9 = 1, находим (без учета сдвига)
б
а
§77 §78 §79 " l/EF 0 0
[§] = §87 §88 §89 = 0 13/l2El 0
_597 §98 §99 _ 0 0 l/El
(41)
Так как при смещении введенных в основной системе связей (см. рис. 5) реакции в них не возникают, то все коэффициенты блока г будут нулевыми
(г* = гп = 0).
Таким образом [г] — квадратная нулевая матрица шестого порядка. Вычисление коэффициентов матрицы откликов, входящих в блок 5, выполняется на основе геометрических соотношений между смещениями введенных связей и вызванными ими перемещениями по направлениям неизвестных усилий в месте разреза стержня (рис. 5, б):
§71 §72 §76
p]= §81 §82 §86 =
§92 §96 __
cos a sin a 0 - sin a - cos a 0
= sin a - cos a ll 2 cos a - sin a l/ 2
0 0 -1 0 0 1
(42)
Коэффициенты блока г находятся на основе теоремы о взаимности реакций и перемещений, т.е. Г* = —5 *
[ Г ] = —[5] Г. (43)
Расчет стержневой системы по МКЭ в форме смешанного метода выполняется не только на воздействие узловой нагрузки, но и на нагрузку в поле КЭ.
От этой нагрузки в основной системе КЭ возникают отклики {гР } и {5Р }. Элементы подвектора |гР } = Р, г2 Р, г3р, г4Р, г5Р, г6Р J находятся как реакции в связях основной системы КЭ от нагрузки (как узловой, так и в поле КЭ), а элементы подвектора |5Р | = |^57Р, 58Р, 59Р ] находятся как перемещения
от той же нагрузки, но по направлению силовых неизвестных в основной системе.
Прямоугольный КЭ изгибаемой пластинки. Рассмотрим прямоугольный КЭ ортотропной плиты. Основная система смешанного метода для этого элемента показана на рис. 6.
За неизвестные примем вертикальные линейные смещения Z узловых точек КЭ и по два изгибающих момента в направлении осей ох и оу соответственно в этих же угловых точках ((Мх, Му) [15, 17]. Необходимо при этом помнить, что на элемент в общем случае узловой нагрузки будут действовать по кромкам распределенные моменты, являющиеся линейными функциями координат х,у. Значения этих моментов в узловых точках — М'х, М'у, (^ = 7, ], к, I) принимаются за неизвестные. Направления этих моментов соответствуют растяжению нижних волокон КЭ. Итого для рассматриваемого КЭ будем иметь 12 неизвестных:
= [ q2... <?п ]
Рис. 6. Основная система смешанного метода для прямоугольного КЭ ортотроп-ной плиты
В этом векторе компоненты qt ( = 1, 2, 3, 4) относятся к неизвестным линейным смещениям ^ q2 q3 q4 ] = [zi 2. 2к 21 ], а компоненты ^ = 5, 6,..., 12) — к силовым неизвестным
|> йв... ] = [мх му мх му мкх мку мх му].
При узловой нагрузке функцию прогибов КЭ 2 = w( х, у) примем в том же виде, что и при расчете по МКЭ в перемещениях, т.е. в виде неполного бикубического полинома с 12 произвольными параметрами:
2 = w(x, у) = [ х, у )]{а}, (44)
[ФДх,у)] = [1 х у х2 у2 ху х2у ху2 х3 у3 х3у ху3 ], {а} = [а1 а2 ...а12].
Вектор кривизны срединной поверхности КЭ определяется выражением {х} = [В]{а}, (45)
где
где [В ] = (8}[ФД х, у)], {8} =
82 82
— матрица дифференцирования.
8х2 8у2 8х8у
Вектор усилий и вектор деформаций (моментов и кривизн в рассматриваемом случае) связаны зависимостью
\ Мх ~
Му[ = -[С ]{х} = -[С ][ В ]{а}, (46)
М_.
{N } =
где [С] — матрица упругих свойств материала. Для ортотропного материала
[С ] =
А А 0
в в
ц у
0
0 0 2В,
(41)
D =
Eh
D, = -
Eyh3
12 (1 -^ц )' у 12 (1 -ц^у)
— изгибные жесткости КЭ, в направле-
нии осей ox и oy соответственно; D = цxDx = ц D , Dk =
Gh3
-; ц x, ц у — ко-
[У ] =
ц г х х к
6
эффициенты Пуассона; Л — толщина пластинки.
В матрице [С] последняя строка умножена на два, чтобы учесть крутящие моменты М и М = М .
ух ху ух
Так как вектор {д} может быть получен из (44) и (45) путем подстановки координат узлов, то связь между векторами {д} и {а} определяется выражением
{д} = [^]{а}, (48)
в котором в соответствии с принятой нумерацией неизвестных в основной системе
А (х, у, У Ф1 ((, У_)
Ф1 (, Ук)
_Ф1 (я у
Ф2 (, У, ) _Ф_3_(_Х£ у Ф2 ( ,У;~)
Ф_3_((_ ,У_) Ф2_(хл,Ук ) Ф_3 _(х ) Ф2 (,У/) Ф3 (х/, У/у
Ф2 (x, У У= Ох ^ [Ф1 (x, У)] + В ^ [Ф1 (x, У)],
Ф3 (x, У У= Вц ^ [Ф1 (X, У)] + ВУ ^ [Ф1 (^ У)]. Из (48) следует {а} = [У1 ]-1 {д}.
Подставляя это выражение в (44), (45), (46), получаем: г = м>(х х, У) = [и (х, У )]{д};
М = [В ]{д}; {М } = -[С ][ В ]{д}, где [[ ] = {5}[ы(х,У)], [и(х,у)] = [ФДх,у)][[] или [В,] = [В][Уг]
(49)
(50)
(51)
1-г
вестник 12/2014
Процесс построения матрицы откликов КЭ остается в принципе таким же, как и при построении матрицы жесткости.
В линейно-деформируемых системах любой узловой отклик в смешанной форме МКЭ можно представить в виде линейной однородной функции от узловых воздействий (силовых, деформационных, температурных и т.д.)
{^ } = [<! ]{д}, (52)
где — вектор узловых откликов КЭ; [с1] — матрица откликов КЭ на единичные воздействия; {д} — вектор узловых воздействий (вектор смешанных неизвестных).
Для установления этой взаимно однозначной связи применим к рассматриваемому элементу одновременно принцип возможных перемещений и принцип возможных изменений напряженного состояния, принимая в качестве возможных изменений вектор смешанных неизвестных {д} вариации его действительных компонентов, а в качестве действительного состояния — вектор узловых откликов {^} = [Fi ¥у ¥к¥г ] ,
где = ЛХ Лу ]Т , X = 7, у, к, I.
Тогда возможная работа внешних воздействий может быть записана формально следующим образом:
5А = [5{д}]Т {^}. (53)
Аналогично находится возможная работа внутренних сил по всей области элемента:
а Ь
5Ж = — |{[5{х}]Т {Ы}йхйу. (54)
0 0
Подставляя в (48), (49) векторы {к} и {Ы} по (46), получаем
5Ж = — [5{д}]Т ({ { [Я, ]Т [С ][Я ] dxdy 1{д}. (55)
V оо у
Так как 5А + 5Ж = 0, то
а Ь
И =
{{[А ]Т [С ][Д ]dxdy
{д}. (56)
Сравнивая это выражение с (47), получаем
а Ь
[ ] = {{[Я ]Т [С ][Я ]dxdy. (57)
0 0
Подстановка матриц [Я,] и [С] в это выражение и интегрирование по площади дает значения элементов матрицы откликов КЭ, относящихся к блокам г и 5.
Непосредственное использование выражения (57) для нахождения элементов матрицы откликов КЭ, относящихся к блокам г и 5, является, как это показано в [15], некорректным.
Поэтому для нахождения коэффициентов блока г примем в качестве возможных состояний функции формы КЭ при единичных смещениях линейных связей = 1 (X = 1,2,3,4) в основной системе метода перемещений для данного КЭ.
Эти функции формы определяются выражением, используемым при получении матрицы жесткости КЭ:
[ (х, у )] = !>, (х, у )][Г ]-1 {д}. (58)
Параметры всех функций форм определяются матрицей [V] , из которой берутся лишь столбцы, соответствующие линейным смещениям (1, 4, 7, 10) [15, 18, 19].
Кривизны срединной поверхности КЭ, соответствующие этим функциям форм, определяются выражением
[Х( х, у )] = [ В2 ] [д], (59)
где [ В2 ] = [В ][¥ ]-1.
В соответствии с этим выражение для определения коэффициентов блока г следует записать в виде
[d ] = {{[ B][C ][ B2 ] dxdy = [?].
(60)
В качестве возможных состояний могут быть приняты также функции формы на основе полиномов Эрмита, используемых иногда при расчете по МКЭ в перемещениях.
Аналогично можно получить матрицу откликов для прямоугольного элемента с четырьмя неизвестными в каждом узле.
В этом случае в каждом узле элемента будут четыре неизвестных — прогиб срединной плоскости ^ и три усилия (изгибающие моменты в направлении осей X—М и У — М и обобщенный крутящий момент — Н = М = М ).
х у ^ ху ух'
Сходимость МКЭ в форме смешанного метода. Сходимость МКЭ в форме смешанного метода рассмотрим на том же примере, который использован для исследования сходимости МКЭ в перемещениях. Основная система смешанного метода для рассматриваемой балки с тем же количеством конечных элементов приведена на рис. 7. За неизвестные в основной системе приняты вертикальные перемещения узлов г. и опорные моменты М.
Рис. 7. Основная система смешанного типа
Система разрешающих уравнений МКЭ совпадает в данном случае с уравнениями смешанного метода, и поэтому можно сразу записать уравнения совместности деформаций и равновесия для .-го узла и двух КЭ (см. рис. 7) системы
(61)
r = яЛ + яЛ = 0, еЛ-еП = 0.
вестник 12/2014
В случае равномерно распределенной нагрузки, т.е. q = const имеем
Mtl + M J z t - zql
3EJ 6EJ l 24EJ
Mil M>J Z+1 - z ql 3
' 3Е/ 6EJ I 24EJ кл = Мы-Щ_+ & кп = М1+! -М, + & ' I 2' ' I 2 '
Подставив эти выражения в (61) получим после преобразований систему конечно-разностных уравнений МКЭ в смешанной форме:
M-1 + 4 M + M+i + zi-1- 2 zi + z+1 _ 0.
6EJ l2 '
Rл _ M-i - 2M + M+i + q R _-1-+ q.
(62)
При предельном переходе (/ ^ dx), получаем:
r z, j-2z, + zi+1 » Mt 1 + 4M, + Mt+1
lim —-i-— = z, , lim——-'--— = M■,
i^dx 12 - ' i^dx 6EJ
.. Mr i -2M- + M-+j »
lim——---— = M .
i ^dx i -
Таким образом, система уравнений (62) преобразуется к виду
M + EJz = 0, M" + q = 0. (63)
Система дифференциальных уравнений (63) разрешима. Это означает, что при предельном переходе система алгебраических уравнений (62) не вырождается и при сгущении сетки КЭ счет будет устойчивым. Следовательно, при решении задачи по МКЭ в смешанной форме математическая модель в виде системы уравнений (62) сходится к точному решению.
Это позволяет сделать вывод о том, что вырождение системы алгебраических уравнений МКЭ при предельном переходе не является специфическим свойством этого метода вообще, а проявляется лишь в некоторых конкретных случаях. В то же время полученные результаты свидетельствуют о том, что МКЭ в смешанной форме обеспечивает получение более устойчивых результатов при сгущении сетки конечных элементов, чем МКЭ в форме метода перемещений. Вполне очевидно, что те же свойства будут проявляться и в двумерных системах.
Библиографический список
1. Городецкий А.С., Заворицкий В.И., Лантух-Лященко А.И., Рассказов А.О. Метод конечных элементов в проектировании транспортных сооружений. М. : Транспорт, 1981. 143 с.
2. Постнов В.А., Хархурим И.Я. Метод конечных элементов в расчетах судовых конструкций. Л. : Судостроение, 1974. 344 с.
3. Секулович М. Метод конечных элементов / пер. с серб. Ю.Н. Зуева ; под ред. В.Ш. Барбакадзе. М. : Стройиздат, 1993. 664 с.
4. Игнатьев А.В. Основные формулировки метода конечных элементов в задачах строительной механики. Часть 1 // Вестник МГСУ 2014. № 11. С. 37—57.
5. Bogner F., Fox R., Schmit L. A cylindrical shell discrete element. AIAA J. 1967. Vol. 5. No. 4. Pp. 745—750.
6. Олман Д.Дж. Треугольные конечные элементы для расчета изгибаемых пластин при постоянных и линейно распределенных изгибающих моментах // Расчет упругих конструкций с использованием ЭВМ / пер. с англ. ; под ред. А.П. Филина. Л. : Судостроение, 1974. С. 80—101.
7. Клочков Ю.В. Развитие теории линейного и нелинейного деформирования оболочек на основе МКЭ с учетом смещения как жесткого целого и изменения толщины : дисс. ... д-ра техн. наук. Волгоград : Волгоградская ГСХА, 2001. 326 с.
8. Батэ К., Вилсон Е. Численные методы анализа и метод конечных элементов / пер. с англ. А.С. Алексеева, О.О. Андреева, В.Н. Сидорова ; под ред. А.Ф. Смирнова. М. : Стройиздат, 1982. 448 с.
9. Бате К. Методы конечных элементов / пер. с англ. В.П. Шидловского ; под ред. Л.И. Турчака. М. : Физматлит, 2010. 1024 с.
10. Цыбенко А.С. Применение треугольных трехузловых несогласованных элементов для решения осесимметричных задач теории упругости // Проблемы пластичности. 1986. № 3. С. 79—83.
11. Semenov V.A., SemenovP.Yu. Hybrid finite elements for analysis of shell structures // Proc. International Congress ICSS—98, 22—26 June 1998, Moscow, Russia. Moscow, 1998. Vol. 1. Pp. 244—251.
12. Bathe K.J. Finite Element Procedures. Prent. Hall, Englewood Cliffs. 1996. 1036 p.
13. Fraeijs de Veubeke B., Sander G. An equilibrium model for plate bending // International J. Solids and Structures. 1968. Vol. 4. No. 4. Pр. 447—468.
14. Тюкалов Ю.А. Решение задач строительной механики методом конечных элементов в напряжениях на основе функционала дополнительной энергии и принципа возможных перемещений : дисс. ... д-ра техн. наук. Киров : ВятГУ, 2006. 314 с.
15. Игнатьев В.А., Игнатьев А.В., Жиделев А.В. Смешанная форма метода конечных элементов в задачах строительной механики. Волгоград : ВолгГАСУ, 2006. 172 с.
16. Игнатьев А.В., Габова В.В. Алгоритм статического расчета плоских стержневых систем по методу конечных элементов в смешанной форме // Вестник Волгоградского государственного архитектурно-строительного университета. Серия: Естественные науки. 2007. Вып. 6 (23). С. 72—77.
17. Рекунов С.С., ВоронковаГ.В. Особенности расчета пластинок по методу конечных элементов в смешанной форме // Вестник Волгоградского государственного архитектурно-строительного университета. Серия: Строительство и архитектура. 2007. Вып. 7 (26). С. 74—77.
18. Масленников А.М. Расчет строительных конструкций численными методами. Л. : ЛГУ, 1987. 224 с.
19. Покровский А.А. Смешанная форма МКЭ в расчетах стержневых систем и сплошной среды : дисс. ... д-ра техн. наук. Пенза : ПГАСА, 2000. 308 с.
Поступила в редакцию в октябре 2014 г.
Об авторе: Игнатьев Александр Владимирович — кандидат технических наук, доцент, доцент кафедры строительной механики, Волгоградский государственный архитектурно-строительный университет (ФГБОУ ВПО «ВолгГАСУ»), 400074, г. Волгоград, ул. Академическая, д. 1, [email protected].
Для цитирования: Игнатьев А.В. Основные формулировки метода конечных элементов в задачах строительной механики. Часть 2 // Вестник МГСУ. 2014. № 12. С. 40—59.
becthmk
MfCY_12/2014
A.V. Ignat'ev
MAIN FORMULATIONS OF THE FINITE ELEMENT METHOD FOR THE PROBLEMS OF STRUCTURAL MECHANICS Part 2
The author offers a classification of Finite Element formulations, which allows orienting in a great number of the published and continuing to be published works on the problem of raising the efficiency of this widespread numerical method. The second part of the article offers examination of straight formulations of FEM in the form of displacement approach, area method and classical mixed-mode method. The question of solution convergence according to FEM in the form of classical mixed-mode method is considered on the example of single-input single-output system of a beam in case of finite element grid refinement.
The author draws a conclusion, that extinction of algebraic equations system of FEM in case of passage to the limit is not a peculiar feature of this method in general, but manifests itself only in some particular cases. At the same time the obtained results prove that FEM in mixed-mode form provides obtaining more stable results in case of finite element grid refinement in comparison with FEM in the form of displacement approach. It is quite obvious that the same qualities will appear also in two-dimensional systems.
Key words: finite elements method, Finite Element formulations classification, structural mechanics.
References
1. Gorodetskiy A.S., Zavoritskiy V.I., Lantukh-Lyashchenko A.I., Rasskazov A.O. Metod konechnykh elementov v proektirovanii transportnykh sooruzheniy [Finite Element Method in Transport Constructions Design]. Moscow, Transport Publ., 1981, 143 p. (In Russian)
2. Postnov V.A., Kharkhurim I.Ya. Metod konechnykh elementov v raschetakh sudovykh konstruktsiy [Finite Element Method in Ship Structures Calculation]. Leningrad, Sudostroenie Publ., 1974, 344 p. (In Russian)
3. Sekulovich M. Metod konechnykh elementov [Finite Element Method]. Translated from Serbian Yu.N. Zueva, editor V.Sh. Barbakadze. Moscow, Stroyizdat Publ., 1993, 664 p. (In Russian)
4. Ignat'ev A.V. Osnovnye formulirovki metoda konechnykh elementov v zadachakh stroitel'noy mekhaniki. Chast' 1 [Essential FEM Statements Applied to Structural Mechanics Problems. Part 1]. Vestnik MGSU [Proceedings of Moscow State University of Civil Engineering]. 2014, no. 11, pp. 37—57. (In Russian)
5. Bogner F., Fox R., Schmit L. A Cylindrical Shell Discrete Element. AIAA. 1967, vol. 5, no. 4, pp. 745—750. DOI: http://dx.doi.org/10.2514/3.4056.
6. Allman D.J. Treugol'nye konechnye elementy dlya rascheta izgibaemykh plastin pri postoyannykh i lineyno raspredelennykh izgibayushchikh momentakh [Trigonal Finite Elements for Bending Plates Calculation in Case of Permanent and Linearly Distributed Bending Moments]. Raschet uprugikh konstruktsiy s ispol'zovaniem EVM [Calculation of Elastic Structures Using Computer]. Translated from English, editor A.P. Filin. Leningrad, Sudostroenie Publ., 1974, pp. 80—101. (In Russian)
7. Klochkov Yu.V. Razvitie teorii lineynogo i nelineynogo deformirovaniya obolochek na osnove MKE s uchetom smeshcheniya kak zhestkogo tselogo i izmeneniya tolshchiny [Development of the Theory of Linear and Non-linear Deformation of Shells Basing on FEM with Account for the Displacement as Stiff Entire and Change of the Width]. Dissertation of the Doctor of Technical Sciences. Volgograd, Volgogradskaya GSKhA Publ., 2001, 326 p. (In Russian)
8. Bathe K.-J., Wilson E.L. Numerical Methods in Finite Element Analysis, New Jersey, Prentice-Hall, 1976, 528 p.
9. Bathe K. Metody konechnykh elementov [Finite Elements Methods]. Transl. from English by V.P. Shidlovskiy. Moscow, FIZMATLIT Publ., 2010, 1024 p. (In Russian)
10. Tsybenko A.S. Primenenie treugol'nykh trekhuzlovykh nesoglasovannykh elementov dlya resheniya osesimmetrichnykh zadach teorii uprugosti [Application of Trigonal Nonconforming Elements for Solving Axisymmetric Tasks of Elasticity Theory]. Problemy plastichnosti [Elasticity Problems]. 1986, no. 3, pp. 79—83.
11. Semenov V.A., Semenov P.Yu. Hybrid Finite Elements for Analysis of Shell Structures. Proc. International Congress ICSS—98, 22—26 June 1998, Moscow, Russia. Moscow, 1998, vol. 1, pp. 244—251.
12. Bathe K.J. Finite Element Procedures. Prent. Hall, Englewood Cliffs, 1996, 1036 p.
13. Fraeijs de Veubeke B., Sander G. An Equilibrium Model for Plate Bending. International J. Solids and Structures. 1968, vol. 4, no. 4, pp. 447—468. DOI: http://dx.doi. org/10.1016/0020-7683(68)90049-8.
14. Tyukalov Yu.A. Reshenie zadach stroitel'noy mekhaniki metodom konechnykh el-ementov v napryazheniyakh na osnove funktsionala dopolnitel'noy energii i printsipa voz-mozhnykh peremeshcheniy [Solving the Tasks of Structural Mechanics by Finite Element Method in Strains Basing on Additional Energy Functional and Principle of Possible Sisplace-ments]. Dissertation of the Doctor of Technical Sciences. Kirov, VyatGY Publ., 2006, 314 p. (In Russian)
15. Ignat'ev V.A., Ignat'ev A.V., Zhidelev A.V. Smeshannaya forma metoda konechnykh elementov v zadachakh stroitel'noy mekhaniki [Mixed Form of Finite Element Method in Problems of Structural Mechanics]. Volgograd, VolgGASU Publ., 2006, 172 p. (In Russian). (In Russian)
16. Ignat'ev A.V., Gabova V.V. Algoritm staticheskogo rascheta ploskikh sterzhnevykh sistem po metodu konechnykh elementov v smeshannoy forme [Algorithm of Static Analysis of Flat Truss Structures Using Finite Element Method in Mixed Mode Form]. Vestnik Volgo-gradskogo gosudarstvennogo arkhitekturno-stroitel'nogo universiteta. Seriya: Estestvennye nauki [Proceedings of Volgograd State University of Architecture and Civil Engineering. Series: Natural Sciences]. 2007, no. 6 (23), pp. 72—77. (In Russian)
17. Rekunov S.S., Voronkova G.V. Osobennosti rascheta plastinok po metodu konechnykh elementov v smeshannoy forme [Features of Calculating Plates Using Finite Elements Method in Mixed-Mode Form]. Vestnik Volgogradskogo gosudarstvennogo arkhitekturno-stroitel'nogo universiteta. Seriya: Stroitel'stvo i arkhitektura [Proceedings of Volgograd State University of Architecture and Civil Engineering. Series: Construction and Architecture]. 2007, no. 7 (26), pp. 74—77. (In Russian)
18. Maslennikov A.M. Raschet stroitel'nykh konstruktsiy chislennymi metodami [Calculation of Building Structures Using Numerical Methods]. Leningrad, LGU Publ., 1987, 224 p. (In Russian)
19. Pokrovskiy A.A. Smeshannaya forma MKE v raschetakh sterzhnevykh sistem i sploshnoy sredy [Mixed-Mode Form of FEM in Calculation of Truss Systems and Continuous Medium]. Dissertation of the Doctor of Technical Sciences. Penza, PGASA Publ., 2000, 308 p. (In Russian)
About the author: Ignat'ev Aleksandr Vladimirovich — Candidate of Technical Sciences, Associate Professor, Department of Structural Mechanics, Volgograd State University of Architecture and Civil Engineering (VSUACE), 1 Akademicheskaya str., Volgograd, 400074, Russian Federation; [email protected].
For citation: Ignat'ev A.V. Osnovnye formulirovki metoda konechnykh elementov v za-dachakh stroitel'noy mekhaniki. Chast' 2 [Main Formulations of the Finite Element Method for the Problems of Structural Mechanics. Part 2]. Vestnik MGSU [Proceedings of Moscow State University of Civil Engineering]. 2014, no. 12, pp. 40—59. (In Russian)
