Метод конечных элементов

Введение Метод конечных элементов (МКЭ) – численный метод решения дифференциальных уравнений, широко использующийся в различных областях техники (ракето- и самолетостроение, кораблестроение, строительство и др.). Основоположником теории МКЭ считается Курант (1943 г.). М. Тернер, Х. Мартин и др. внедрили МКЭ в строительную механику и механику сплошных сред (конец пятидесятых – начало шестидесятых годов двадцатого века). Существенно расширили область применения МКЭ Б. Сабо, О. Зенкевич и др. (конец шестидесятых – начало семидесятых годов), показав, что его можно использовать для решения любых дифференциальных уравнений. Большой вклад в развитие МКЭ внесли отечественные ученые Л. Розин, В. Корнеев, В. Постнов и др. Развитие МКЭ неразрывно связано с совершенствованием вычислительной техники, ускоряющей сложные численные расчеты. Соответственно совершенствовались вычислительные программы, реализующие этот метод. Наиболее распространенными программами расчета конструкций на основе МКЭ являются в настоящее время COSMOS, ЛИРА, STARK (строительные конструкции). Основная идея МКЭ. Преимущества метода Основная идея МКЭ состоит в том, что любую непрерывную в некоторой области величину (например, внутреннее усилие в фундаментной балке, перемещение в плите перекрытия и т.п.) можно аппроксимировать дискретной моделью, которая создается из множества кусочно-непрерывных функций, определенных в конечном числе подобластей (элементов). Обычно такими функциями являются полиномы – линейные, квадратичные, кубичные, … . Кусочно-непрерывные функции строятся с помощью значений непрерывной величины в точках соединения элементов (в узлах). Таким образом, чтобы определить неизвестную непрерывную величину, нужно определить ее значения в узлах. Основные этапы создания дискретной модели неизвестной величины следующие: 1. В исследуемой области задается конечное число точек (узлов). 2. Значения непрерывной величины в каждом узле считаются неизвестными, которые должны быть определены. 3. Исследуемая область разбивается на конечное число подобластей (элементов), имеющих общие точки (узлы). 4. Непрерывная величина в каждом элементе аппроксимируется полиномом, который определяется с помощью узловых значений этой величины: для каждого элемента определяется свой полином, но его коэффициенты подбираются так, чтобы сохранялась непрерывность величины на каждой границе элемента. Основную идею МКЭ иллюстрирует следующий пример. Рассматриваются прогибы v в стержне (рис. 1). q х v(x) l Рис. 1 Непрерывная величина – функция прогиба  (x) . Ее область определения (исследуе- мая область) – стержень длиной l. Задается пять точек (узлов). Фиксируются прогибы в каждом узле: v1, v2 ,..., v5 (см. рис. 2). v1 v2 1 2 v3 3 v4 v5 4 5 Рис. 2 Аппроксимирующая функция – линейный по х полином, так как на каждый элемент приходится по два узла. Окончательная аппроксимация  (x) - четыре кусочно-линейные функции, каждая из которых определена на отдельном элементе (рис. 3). v v1 1 v2 2 v3 v4 v5 x 3 2 Рис. 3 4 5 Неизвестные узловые значения  (x) должны быть отрегулированы таким образом, чтобы приближение к истинной функции  (x) было наилучшим. Это осуществляется минимизацией некоторой величины, связанной с физической сущностью задачи. Процесс минимизации сводится к решению систем линейных алгебраических уравнений относительно узловых значений  (x) . Если неизвестная непрерывная величина  определена в двух – или трехмерной области, аппроксимирующими функциями являются функции от х и у или от х, у и z соответственно. Двумерная область разбивается обычно на элементы в форме треугольника или четырехугольника, трехмерная область – на элементы в форме тетраэдра или параллелепипеда. Аппроксимирующие функции изображаются в таком случае плоскими (см. рис. 4, а) или криволинейными (см. рис. 4, б) поверхностями (двумерная область).   у у х х б) а) Рис. 4 Из сказанного выше следует, что основными преимуществами МКЭ являются следующие: 1. Возможность исследовать тела (конструкции), составленные из нескольких материалов (так как свойства материалов соседних элементов могут быть разными). 2. Возможность исследовать области (конструкции) любой формы (так как криволинейная область аппроксимируется прямолинейными элементами или точно описывается криволинейными элементами). 3. Возможность учета различных граничных условий: с разрывной нагрузкой, смешан3 ных. 4. Возможность составления общих методик и программ для решения различных по физике задач одного определенного вида (например, программа осесимметричной задачи о распространении тепла может быть использована для решения любой задачи данного типа: о распределении напряжений в осесимметричной конструкции и т.п.). Дискретизация области Разбиение области на подобласти – первый этап в решении задачи МКЭ. Разбиение области требует инженерных навыков и опыта. Неудачное разбиение приведет к ошибочным результатам решения задачи. При разбиении области необходимо уже иметь некоторые общие представления о результатах решения задачи, чтоб уменьшить размер элементов в тех частях области, где ожидаемый результат может резко меняться, и увеличить размеры элементов в тех частях, где ожидаемый результат близок к постоянному. Вообще, при разбиении области всегда идет поиск «золотой середины»: с одной стороны, элементы должны быть достаточно малыми, чтоб получить результаты необходимой точности; с другой стороны, чем крупнее элементы, тем меньше вычислительной работы. а. Одномерные элементы Одномерный элемент – это стержневой элемент. Он используется при расчете стержневых конструкций (фермы, балки, рамы, …). Одномерные элементы могут быть с двумя узлами (см. рис. 5, а), тремя (квадратичные) (рис. 5, б), четырьмя (кубические) (рис. 5, в). а) б. Двумерные элементы б) Рис. 5 Основные виды двумерных элементов – треугольные и четырехугольные (рис. 6). 4 Рис. 6 Толщина элементов может быть постоянна или является функцией координат. Такие элементы используются при расчете различных пластин (плиты перекрытий, стеновые панели и т.п.). в. Трехмерные элементы Наиболее часто используются трехмерные элементы в виде тетраэдра или параллелепипеда (рис. 7). а) б) в) Рис. 7 При рассмотрении конструкций специфической формы (например, осесимметричных) используются специальные элементы (рис. 8). z  r Рис. 8 Разбиение одномерной конструкции на элементы не представляет трудностей. При разбиении двумерной области чаще используются треугольные элементы. 5 Сначала область делится на треугольные и четырехугольные подобласти (зоны): границы между зонами определяются изменением геометрии области, изменением нагрузки, свойств материалов. Затем зоны разбиваются на элементы. Наиболее легко разбить треугольную зону, выбрав определенное количество узлов вдоль каждой стороны зоны и соединив соответствующие узлы линиями (см. рис. 9). а) б) Рис. 9 Если на каждой стороне такой зоны выбрано по п узлов, число полученных треуголь2 ных элементов - п  1 . Четырехугольные зоны обычно разбиваются соединением узлов на противоположных сторонах (рис. 10). Если число узлов на двух противоположных сторонах такой зоны одинаково и равно п или т для двух пар противоположных сторон, число прямоугольных элементов этой зоны будет 2n  1m  1 . Желательное разбиение а) Нежелательное разбиение б) Рис. 10 Четырехугольные элементы в дальнейшем можно разбивать на треугольные: проведением более короткой диагонали в четырехугольнике. Разбиение более короткой диагональю дает элементы, наиболее близкие к равностороннему треугольнику, что ведет к более точ6 ным результатам. Равномерное разбиение, когда все элементы имеют одинаковые размеры, используется редко. Обычно из-за ожидаемой концентрации напряжений размеры элементов варьируются, и эта возможность является важным достоинством метода конечных элементов. Наиболее простой способ изменения размеров элементов – применение черырехугольных зон с различным числом узлов на противоположных сторонах (рис. 11). Рис. 11 Порядок нумерации узлов влияет на эффективность вычислений метода конечных элементов. Реализация метода конечных элементов приводит к решению системы алгебраических уравнений, большое число коэффициентов которой нулевые. В матрице коэффициентов этой системы все ненулевые коэффициенты заключены между линиями, параллельными главной диагонали матрицы (рис. 12). Ширина полосы a a  a  0 0  0  0 a a a a a a a a a a a a a a a 0 0  0  0 a  a a  Рис. 12 Расстояние между главной диа- гональю и одной из этих линий называ- ется шириной полосы матрицы (см. рис. 12) и является показателем эффективности вычис7 лений: чем ширина полосы уже, тем меньше размер требуемой машинной памяти и время вычислений. Ширина полосы определяется порядком нумерации узлов: B  R  1Q , где В – ширина полосы, R – наибольшая для исследуемой области разница между номерами узлов в элементе, Q - число неизвестных в каждом узле. Таким образом, при нумерации узлов необходимо стремиться к тому, чтобы разница между номерами узлов в элементах области была как можно меньше. На примере простой прямоугольной области это проиллюстрировано рис. 13. 1 9 17 25 2 3 4 5 6 7 10 11 12 13 14 15 16 18 19 20 21 22 23 24 8 26 27 28 29 30 31 32 9 5 1 13 17 21 25 29 2 6 10 14 18 22 26 30 3 7 11 15 19 23 27 31 8 12 16 20 24 28 32 4 б) а) Рис. 13 Вариант а) менее эффективен: наибольшая разница между номерами узлов в элементе 8. Вариант б) более эффективен: соответствующая разница – 4. Расчет стержневых систем Рассматривается пространственная стержневая система. Предполагается, что материал стержней идеально упругий. Система является линейно-деформируемой. Стержневая система разбивается на конечное число элементов, соединенных с соседними в узлах. Соединение в узле может быть жестким (рис. 14, а - для плоской стержневой системы) или шарнирным (рис. 14, б - также для плоской стержневой системы). а) б) Рис. 14 8 Жесткий узел обеспечивает равенство всех перемещений концевых сечений элементов, примыкающих к узлу (и линейных перемещений, и углов поворота). Шарнирный узел обеспечивает равенство только линейных перемещений концов элементов. За неизвестные принимаются перемещения узлов стержневой системы. Рассмотрим отдельный конечный элемент с шарнирными узлами i, j (рис. 15). Такие элементы называются ферменными. Свяжем с ним систему координат ху (рис. 15). Такая система координат, привязанная к конкретному конечному элементу, называется местной. y l j i ui ´ x vj vi i j ´ uj Рис. 15 Под действием внешних нагрузок стержневая система деформируется, в том числе деформируется рассматриваемый элемент. Его узлы переместятся в новые положения i , j  (рис. 15). Матрицы-векторы перемещений узлов i и j соответственно u j  ui  qi    , q j    v j  vi  ui    qi  vi  q  u  q j   j  v j  Теперь рассмотрим усилия, действующие на элемент, - это внутренние усилия, заменяющие действие остальной части стержневой системы на наш элемент (рис. 16). y Rix i Рис. 16 Riy j Rjx x Rjy 9 Они образуют матрицы-векторы внутренних усилий для узлов i и j  Rix   R jx  Ri    , R j     Riy   R jy   Rix     Ri   Riy  R   R j   R jx     R jy  Если нагрузка на такой прямолинейный элемент приложена только в узлах, Riy  R jy  0 . Действительно, из условия равновесия элемента  M i  0  R jyl  0  R jy  0  y  0  Riy  R jy  0  R iy x  0  R iх 0  R jx  0  R jx   Riх  Rx Это означает, что продольные усилия одинаковы по величине и противонаправлены. На рис. 17 представлена соответствующая картина с растягивающими усилиями. y Rx j i Rx x Рис. 17 Удлинение элемента l  u j  ui Rx  EA  1  EA  0 R l  1   0 10  l EA  u j  ui l l 0 1 0 0 1 0 0  0 0 q 0  0  Rix   1 R    iy   EA  0  R jx  l  1     0  R jy  0 ui    0 vi   0 u j     0 v j  0 1 1 Обозначим  1  EA  0 KM  l  1   0 0 1 0 0 1 0 0 0 0  0  0 R  KM q Теперь рассмотрим отдельный конечный элемент с жесткими узлами i, j (рис. 18). y j j i vi ´ j i x uj ui z vj i ´ l Рис. 18 Такие элементы называются балочными. В результате деформации стержневой системы под действием внешних нагрузок узлы рассматриваемого элемента переместятся в новые положения i, j  , кроме этого, концевые сечения повернутся (см. рис. 18). Векторы перемещений узлов 11 u j   ui    qi   vi  , q j   v j  ,    j   i    а вектор узловых перемещений элемента  ui  v   i   i  q   . u j  vj     j  Усилия, действующие на элемент, - на рис. 19: y j Rjx Mjz x Miz Rix i z Rjy Riy Рис. 19 Векторы внутренних усилий для узлов i, j:  R jx   Rix    Ri   Riy  , Ri   R jy  ,    M jz   M iz    а для всего элемента  Rix  R   iy   M iz  R . R jx    R jy    M  jz  12 Матрицу жесткости КM балочного элемента получим следующим образом. Будем отдельно рассматривать перемещения узлов по оси х, по оси y, повороты узловых сечений относительно оси z для каждого узла (рис. 20, а – е): y y uj ui i i x  z j i z а) y x б) y  vj i j j vi x j i z z в) y г) y i x j i j x x j i z j i z д) е) Рис. 20 Узловые усилия можно определить, воспользовавшись табличными результатами для расчета рам методом перемещений (рис. 21, а, б): A i A B M RA B RB а) A 12i ; l2 EJ l M A  4i ; M B  2i ; A = 1 A R A  RB  i l M =1 M 6i ; l 6i MB  ; l MA  l i Рис. 21 B 13 R A = 1 R M A б) B B Учитывая направления перемещений узлов, а также то, что смещения узлов нашего балочного элемента не единичные, а vi , v j , i ,  j , получим усилия, соответствующие отдельным смещениям узлов, показанные на рис. 22, а – г. y 6 EJ vi l2 i 12 EJ vi l3 vi 12 EJ vi l3 i 6 EJ vi l2 x j z а) y 12 EJ vj l3 6 EJ vj l2 j vi 6 EJ vj l2 x i j 12 EJ vj l3 z б) y i 4 EJ i l z 2 EJ i l i j 6 EJ i l2 6 EJ i l2 в) Рис. 22 14 x y j 4 EJ j l 2 EJ j l i j 6 EJ j l2 z x 6 EJ j l2 г) Окончание рис. 22 Тогда узловые усилия: Riy  0  ui  12 EJ 6 EJ 12 EJ 6 EJ vi  i  0  u j  vj  j, l3 l2 l3 l2 R jy  0  ui  12 EJ 6 EJ 12 EJ 6 EJ vi  i  0  u j  vj  j, l3 l2 l3 l2 M iz  0  ui  6 EJ 4 EJ 6 EJ 2 EJ vi  i  0  u j  vj  j, l l l2 l2 M jz  0  ui  6 EJ 2 EJ 6 EJ 4 EJ vi  i  0  u j  vj  j. l l l2 l2 Выражения для узловых усилий Rix и R jx получаются с использованием формулы (1): Rix  EA EA ui  0  vi  0  i  uj  0vj  0 j, l l R jx   EA EA ui  0  vi  0  i  uj  0vj  0 j. l l В матричной форме 15  EA  l  12 EJ  Rix   0 R   l3 iy    6 EJ  M iz   0 l2     EA  R jx    R jy   l 12 EJ     M  jz   l3 6 EJ   0 l2  6 EJ l2 4 EJ l EA l 12 EJ l3 6 EJ  l2 EA l 6 EJ l2 2 EJ l     6 EJ   ui   l2   v  2 EJ   i    l    i   u 0   j  v  6 EJ   j    l 2   j  4 EJ  l  12 EJ l3 6 EJ  l2 или  A 12 J   Rix   R  l2 iy    6J  M iz  E  0 l     A R  jx  l   0  12 J  R jy   l2    6 EJ  M jz   0 l2  6J l A 4J 6J  l A 2J 12 J  l2 6J  l 12 J l2 6J  l 0  6 J   ui  l  v    i  2J        i . 0  u j  6J   vj  l    4 J   j   l  В общем виде это уравнение выглядит так же, как для ферменного элемента (3): R  KM q , где  A 12 J   l2  6J  E l KM   l  A  0  12 J  l2  6 EJ  0 l2  16 6J l A 4J 6J  l A 2J 12 J  l2 6J  l 12 J l2 6J  l 0  6J  l   2J  . 0  6J   l  4J   l  (4) - матрица жесткости балочного элемента. Ранее рассматривались стержневые элементы, на которые действует только узловая нагрузка (рис. 16, 18). В действительности внешняя нагрузка на стержневые конструкции может быть и распределенной (рис. 23, а). В таких случаях перед решением задачи МКЭ проводится так называемое приведение внешних нагрузок к узловым усилиям. Так как рассматриваются линейно-деформируемые упругие стержневые системы, можно использовать принцип суперпозиции (принцип независимости действия сил). В соответствие с этим принципом результат действия группы сил равен сумме результатов, полученных от действия каждой силы в отдельности. Представим решение нашей задачи в виде суммы решений для той же стержневой системы от нескольких внешних воздействий, равных в сумме исходному воздействию. Закрепим все узлы системы от всех возможных смещений (рис. 23, б): в данном случае это означает введение жесткой заделки в узел I и шарнирно-неподвижной опоры в узел j. В результате qi  0 , q j  0 . Получившаяся при этом система – это не связанные друг с другом стержневые элементы (рис. 23, в). Определить усилия в такой системе можно относительно просто: каждый отдельный стержень – дважды или трижды статически неопределимая система, которая может быть рассчитана методом сил. В случае простейших внешних нагрузок можно воспользоваться готовыми решениями (рис. 23, г – эпюра изгибающих моментов M изг в отдельных стержнях). q q q q i j h F h/2 F l 17 3 ql 8 ql 2 0,125Fh 8 q q 5 ql 8 F 2 F F 2 0,125Fh в) г) Рис. 23 г) Теперь рассмотрим заданную систему (рис. 23, а) под действием узловых усилий, равных тем, которые приходились в отдельных стержнях на их опоры (см. рис. 24).  ql 2     , 125 Fh  8    5 ql 8 F 2 3 ql 8 F 2 0,125Fh Рис. 24 По принципу суперпозиции решение исходной задачи (рис. 23, а) можно представить в виде суммы решений задачи рис. 23,б и задачи рис. 24. 18 Таким образом, расчет стержневой системы на действие произвольной внешней нагрузки состоит из расчета отдельных ее элементов на действие внешней нагрузки, относящейся к ним, и расчета всей стержневой системы в целом на действие только узловых усилий. Первый расчет, как было сказано выше, относительно прост. Таким образом, задача сводится в основном к расчету стержневой системы на узловые нагрузки. Описанная операция и называется приведением внешних нагрузок к узловым усилиям. Рассмотрим ферменный элемент в некоторой системе координат ху, единой для всей стержневой системы, состоящей из ферменных элементов (рис. 25). Такая система координат называется глобальной. xr yr - местная система координат данного элемента. y yr EA Rcos i R R xr Rsin j  Rcos l Rsin x Рис. 25 В результате деформации узлы элемента переместятся (рис. 26). j yr j i vi vj y xr i uj ui x Рис. 26 Вектор узловых перемещений элемента в глобальной системе координат 19  ui  v  i q . u j     v j  Вектор узловых усилий  Rix   R cos  R    iy    R sin    R  .  R jx   R cos        R jy   R sin   Удлинение элемента     l  u j  ui cos   v j  vi sin  . Тогда R      EA EA l  u j  ui cos   v j  vi sin  . l l Вектор узловых усилий       u j  ui  Rix   R  iy  EA   u j  ui  R   R jx  l  u j  ui     R jy   u j  ui cos2   v j  vi sin  cos cos sin   v j  vi sin 2   cos2   v j  vi sin  cos cos sin   v j  vi sin 2   ui cos2   vi sin  cos   u j cos2   v j sin  cos    EA  ui sin  cos   vi sin 2   u j sin  cos   v j sin 2    l  ui cos2   vi sin  cos   u j cos2   v j sin  cos    2 2   ui sin  cos   vi sin   u j sin  cos   v j sin  20  cos2  sin  cos   cos2   sin  cos   ui      sin 2   sin  cos   sin 2   vi  EA  sin  cos    . l   cos2   sin  cos  cos2  sin  cos  u j       sin 2  sin  cos  sin 2   v j   sin  cos  Матрицу  cos2  sin  cos   cos2   sin  cos    sin 2   sin  cos   sin 2  EA  sin  cos  K l   cos2   sin  cos  cos2  sin  cos     sin 2  sin  cos  sin 2   sin  cos  называют матрицей жесткости ферменного элемента в глобальной системе координат. Коротко R  Kq . Нетрудно проверить, что эта матрица жесткости связана с матрицей жесткости ферменного элемента в местной системе координат K M (2) соотношением K  CKM C , (5) где С - так называемая матрица направляющих косинусов осей xr yr местной системы координат относительно осей ху глобальной системы координат: 0   cos   sin   sin  cos  0   . C  0 cos   sin     sin  cos    0 Рассмотрим балочный элемент в глобальной системе координат (рис. 27). Вектор узловых перемещений элемента в глобальной системе координат 21 ui  vi     i  q  u  . j v   j   j    j j uj y vj yr i xr ui i vi i EJ  j l zr x z Рис. 27 Вектор узловых усилий  Rix   Riy    M  iz  R  R  . jx   R jy   M   jz  Подобно тому, как выше была получена матрица жесткости ферменного элемента в глобальной системе координат, может быть получена матрица жесткости балочного элемента в глобальной системе координат: R  Kq , где 22   12 J 12 J  6J 12 J 12 J  6J     2 2  A  sin    A cos 2   sin 2    sin    sin  cos  A  sin  cos    A cos   2 sin   2 2 2 l l l l  l l         12 J 6 J 12 J 12 J 6 J          A  12 J  sin  cos  2 2 2 2 cos   cos    A sin   cos   cos     A sin   A  sin  cos     2 2 2 2 l l l  l l  l        6J 6J 6J 6J    sin  sin  4 J sin   cos  2 J E l l l l  K  12 J 12 J 6 J 12 J 12 J 6 J        l   2 2 sin  sin 2   A sin   A  sin  cos   A cos 2    sin  cos    A cos   2 sin   2 2 2 l l l l  l l          12 J 6 J 12 J 12 J 6 J         A  12 J  sin  cos  2 2 2 2   A sin   cos    cos  A cos    cos    sin  cos   A sin     2 2 2 2 l l l  l l  l        6 J 6 J 6 J 6 J 4J    sin  cos  2J sin   cos    l l l l l Эта матрица также связана с матрицей жесткости балочного элемента в местной системе координат (4) соотношением (5), при этом  cos   sin   sin  cos    0 C  0  0   0 1 0 cos   sin  0  sin  cos  0 0  0 . 0 0  1 Выше были введены матрицы, характеризующие отдельный стержневой элемент, и получены соотношения между ними. Теперь получим характеристики совокупности стержневых элементов, образующих стержневую систему. Первоначально будем считать элементы не связанными друг с другом. Общее число элементов системы – т.  q1       Вектор узловых перемещений системы q   q r  .       m q   R1       Вектор узловых усилий системы R   R r  .       m R  K1 н у    Матрица жесткости системы K   Kr    у л н и    .     m и K  л Далее рассмотрим элементы, связанные в единую систему (например, раму см. рис. 23,а). Чтобы заставить m несвязанных элементов работать так же, как работает соответствующая система связанных элементов, нужно приравнять перемещения соответствующих узлов этих двух систем. В матричной форме эти соотношения выглядят так: q  Hq , где q – вектор узловых перемещений системы несвязанных элементов (6), Н – так называемая матрица соединения (или логическая матрица), q - вектор узловых перемещений системы связанных элементов. Матрица соединения состоит из единичных матриц Е. Для системы ферменных элементов 1 0 E . 1   Для системы балочных элементов 1 0 0  E  0 1 0  .   0 0 1 Матрица жесткости системы связанных элементов K  H  KH . 25