лекция. глава4_5_6_2018. Макро и микромоделирование глава 4 Метод конечных элементов
 Скачать 1.17 Mb.
|
|
1 2 3 8 7 6 5 4 V 5 U 5 X Y 0 Рисунок 4.5 - Изопараметрический элемент Для определения перемещений произвольной точки внутри элемента задаются некоторой функцией формы, которая определяет их однозначно в зависимости от известных узловых перемещений. Эту зависимость пред- ставим как e N f } ]{ [ } { , (4.44) где T uvw f } { } { - вектор-столбец перемещений некоторой точки внутри элемента; {δ} e - вектор-столбец узловых перемещений элемента, } { } { 1 m e - матрица формы; [N] - число степеней свободы элемента. Каждая компонента матрицы [N] есть функция координат точек внутри элемента и равняется нулю за пределами данного элемента. Опишем процесс построения интерполирующих полиномов, которыми аппроксимируется искомая функция {f} в объеме конечного элемента. 131 1 8 X Y 0 5 5 2 3 4 1 Рисунок 4.6 - Изопараметрический элемент второго порядка Треугольный симплекс-элемент /5, 7/. Допустим, что перемещения то- чек внутри элемента выражаются полиномом первого порядка от их коор- динат x и y. y x v y x u 6 5 4 3 2 1 (4.45) Тогда выражение для перемещений произвольной точки внутри эле- мента имеет вид e k j i N I N I N I f } ]{ [ } { , (4.46) где I – единичная матрица размерности 2х2, а S y c x b a N i i i i 2 (4.47) и так далее. Здесь S – площадь элемента, j k k i i y x y x a , k i i y y b , j k i x x c 132 Выбранная функция перемещений автоматически гарантирует непре- рывность перемещений между смежными элементами. Изопараметрический квадратичный элемент [5,6]. Введем в изопара- метрическом четырехугольном элементе (рисунок 4.7) локальную коорди- натную систему ξ, η, которая удовлетворяет условиям 1 1 , 1 1 Функциям формы для интерполяции перемещений по их узловым зна- чениям на основе [5,6] имеют вид: ). 1 )( 1 ( 2 1 ), 1 )( 1 ( 2 1 ), 1 )( 1 ( 4 1 ) 1 )( 1 ( 4 1 ) 1 )( 1 ( 2 1 2 0 0 2 2 0 0 2 0 0 i i i N N N (4.48) Здесь i – номер узла элемента (i = 1,2,…,8), а i 0 , i 0 , ξ i , η i – значения локальных координат в i –м узле. (X 1 ,Y 1 ) X Y 0 4 e 1 (X 2, Y 2 ) 2 3 (X 3, Y 3 ) (X 4, Y 4 ) Рисунок 4.7 - Четырехугольный элемент Тетраэдральный элемент [6] . Перемещения любой точки определяются тремя компонентами u, v, w (рисунок 4.7). Тогда вектор перемещений име- ет вид T z y x w z y x v z y x u f )} , , ( ) , , ( ) , , ( { } { 133 По аналогии с (4.31) можно записать, например, z y x u 4 3 2 1 (4.49) Перемещения произвольной точки в объеме элемента можно предста- вить в виде e m k j i N I N I N I N I f } ]{ [ } { , (4.50) где I – единичная матрица размерности 3х3, а v z d y c x b a N i i i i i 6 (4.51) и так далее. Здесь V – объем конечного элемента; а коэффициенты a i , b i , c i и d i определяются через значения координат узловых точек элемента. Трехмерный изопараметрический элемент с линейным полем переме- щений [5,6]. Удобнее всего выразить функции формы через координаты узлов на границе элемента и использовать локальные координаты. Они по- казаны на рисункe 1.6 и выбраны так, чтобы стороны конечного элемента совпадали с координатными линиями ±1. Тогда получим следующие функции формы: , } ]{ 8 2 1 [ } { e N I N I N I f (4.52) где I – единичная матрица размерности 3х3, а ) 1 )( 1 )( 1 ( 8 1 0 0 0 i N , i 0 , i 0 , i 0 , (4.53) 134 где ξ i , η i , ζ i – значения локальных координат в i – м узле (i = 1, … ,8). Отметим, что объемный элемент с 24-мя cтепенями свободы позволяет существенно увеличить число неизвестных параметров в выражении для компонентов перемещений (4.51) и тем самым обеспечить лучшую ап- проксимацию действительного НДС в пределах объема конечного элемен- та. Описанные интерполирующие полиномы для каждого конечного эле- мента обеспечивают непрерывность функции U(x,y,z) и ее производных до (m-1) - го порядка включительно. Здесь 2m - порядок определяющего диф- ференциального уравнения. Опишем процедуры построения матриц жесткости рассматриваемых. Треугольный симплекс-элемент. Для того чтобы записать закон Гука для плоской задачи в матричной форме, представим матрицу упругости в виде 2 ) 1 ( 0 0 0 1 0 1 1 ] [ 2 E D (4.54) Такая запись подразумевает, что деформации объединены в вектор в следующем порядке: } { } { xy y x Матрица [B] получается при помощи дифференцирования соотноше- ний (4.45) и ее можно представить в виде ] [ ] [ k j i B B B B (4.55) или k k j j i i k j i k j i B C B C B C C C C B B B S B 0 0 0 0 0 0 2 1 ] [ (4.56) 135 Отметим, что в этом случае матрица [B] не зависит от координат точек внутри элемента и, следовательно, деформациив нем постоянны. Подставив теперь матрицы [B] и [D] в основное уравнение матрицы жесткости и узловых сил и проинтегрировав по объему конечного элемен- та, получим формулу для подсчета матрицы жесткости треугольного эле- мента st B D B K T e ] ][ [ ] [ ] [ , (4.57) где S - площадь элемента; t - толщина элемента, и силы, обусловленной начальной деформацией, st D B R T e } ]{ [ ] [ } { 0 (4.58) Тетраэдральный симплекс-элемент. Выражение для матрицы упругости в трехмерном случае выглядит как 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 ] [D , (4.59) где ) 2 1 )( 1 ( ) 1 ( E , ) 1 ( , 2 ) 1 ( ) 2 1 ( Такая запись подразумевает, что деформации объединены в вектор в следующем порядке T zx yz xy z y x T } { } { Выражение для матрицы [B] получается дифференцированием соотно- шений (4.52) и состоит из четырех блоков ], [ ] [ m k j i B B B B B (4.60) 136 где x N z N y N z N x N y N z N y N x N B i i i i i i i i i i 0 0 0 0 0 0 0 0 0 ] [ (4.61) Выражение для матрицы жесткости, определяемой соотношением, можно проинтегрировать точно, так как компоненты деформаций и напря- жений постоянны внутри элемента. Подматрица матрицы жесткости имеет размерность 3х3 и определяется соотношением e s T r e rs V B D B K ] ][ [ ] [ ] [ (4.62) где V e – объем тетраэдра; (r =1,…,4; s =1,…,4). Узловые силы, обусловленные начальной деформацией, записываются в виде e T e V D B R } ]{ [ ] [ } { 0 , (4.63) Трехмерный изопараметрический элемент. Особенностью определения матрицы жесткости данного элемента является разделение геометрических параметров элемента и физических параметров его материала. Формула для вычисления блока матрицы жесткости равна [4, 5] V s T r e rs dV B D B K ] ][ [ ] [ ] [ , (4.64) где матрица [B r ] определяется соотношением (4.61), a [D] - выражением (4.59); r =1,2,...,8. Преобразуем (4.56) к следующему виду: 137 3 1 3 1 ] [ ] [ n m V m s n r nm e rs dV x N x N D K (4.65) Здесь x n (n = 1,2,3) - глобальные координаты; [D nm ] - матрицы размер- ности 3х3, полученные из матрицы [D] за счет коэффициентов, стоящих на пересечении n-х строк и m-х столбцов. Причем индекс n пробегает после- довательно три значения, соответствующие расположению n i x N в каж- дой строке матрицы [B r ] T , а m пробегает значения, соответствующие рас- положению m i x N в столбцах матрицы [B s ]. Интеграл вида V m s n r nm rs dV x N x N K ] [ , (4.66) вычисляем, переходя к локальным координатам ξ, η, ζ 1 1 1 1 1 1 ] [ d d d I x N x N K m s n r nm rs , (4.67) где z z z y y y x x x I ] [ - матрица Якоби, (4.68) где а[I]=det([I]) – якобиан матрицы. Изопараметрический квадратичный элемент. При определении выра- жения для матрицы жесткости используется подход, описаний выше. В этом случае формулу для вычисления блока матрицы жесткости размерно- сти 2х2 можно представить в виде 2 1 2 1 ] [ ] [ n m V m s n r nm e rs dV x N x N D K (4.69) 138 Здесь x n - декартовы координаты; [D nm ] - матрицы размерности 2х2, полученные из матрицы [D]. При вычислении интегралов типа (4.67) в выражении (4.69) использу- ются локальные координаты ξ и η, а матрица Якоби имеет вид y y x x I ] [ (4.70) Одно из главных преимуществ формул (4.66) и (4.56) проявляется при решении физически нелинейных задач, когда изменяется только матрица упругости. Интегралы типа (4.67) вычисляются численно с использовани- ем квадратур Гаусса [5]: - для плоского элемента 1 1 1 1 )] , ( [ d d G I , (4.71) - для объемного элемента 1 1 1 1 1 1 )] , , ( [ d d d G I , (4.72) где G – некоторая функция, выражения для которой зависят от типа эле- мента, числа точек интегрирования и даны в [5]. Отметим, что, используя формулы (4.66) и (4.69), не представляет тру- да рассчитать матрицы жесткости и для элементов более высокого поряд- ка. Имея матрицы жесткости отдельных элементов, можно получить гло- бальную матрицу жесткости рассматриваемой области [K e ]. Для получения используется процедура объединения по элементам [5,7], то есть элемент- ная матрица [К] е прибавляется к системной матрице [К е ] сразу же после вычисления. Вектор-столбец узловых сил получается аналогичным обра- зом. 139 Обычное решение i i i i i i i i i d d d } { } { } { } { } { } { } { } { } { 1 1 1 , полученное МКЭ, является приближением к истинному или точному решению теории упругости. Сформулируем критерии сходимости, которые непосредственно сле- дуют из сказанного выше, и состоят в следующем: 1. Для того чтобы выполнялось требование сходимости, необходимо, чтобы представления перемещения и любой его производной, появляю- щейся в функционале, стремились к их точным значениям на каждом эле- менте по мере стремления к нулю размера элемента. 2. Условие сходимости состоит в том, что по мере стремления к нулю размера элемента члены с производными и функцией в функционале должны стремиться к функции той же гладкости, что и точное решение (предполагается непрерывность точного решения). Критерии 1 (полноты) и 2 (согласованности) являются достаточными условиями сходимости вариационного метода конечных элементов. При условии выполнения критериев полноты и согласованности путем доста- точного уменьшения размера элемента теоретически можно достичь лю- бой требуемой точности решения. 4.4 Другие формулировки МКЭ Несмотря на то, что приближенная минимизация функционала - самый распространенный способ подхода к МКЭ, это никоим образом не означа- ет, что такой подход является единственно возможным. Наиболее популярными из других вариантов МКЭ являются методы, позволяющие получить основные соотношения МКЭ непосредственно из дифференциальных уравнений задачи; возможные преимущества таких методов состоят в том, что: а) исчезает необходимость искать функциональный эквивалент извест- ным ДУЧП; б) эти методы могут быть распространены на задачи, для которых функционал - либо вообще не существует, либо пока еще не получен. 140 Рассмотрим задачи приближенного решения системы ДУЧП, которой должна удовлетворять неизвестная функция в области V . Запишем основ- ное уравнение в виде 0 ) ( L , (4.73) а граничное условие на границе С как 0 ) ( C (4.74) Если пробная функция, удовлетворяющая граничным условиям, запи- сана в общей форме (4.3) N ˆ , (4.75) где [N]- матрица базисных функций; { } - основные неизвестные задачи, то в общем случае 0 ˆ R L , (4.76) где ) ˆ ( ) ( L L R - невязка решения. Наилучшим решением будет то, которое дает во всех точках области V наименьшую невязку R . Очевидно, что решение можно получить, использовав то обстоятель- ство, что если невязка R тождественно равна нулю всюду в области, то 0 v WRdV , (4.77) где W – любая функция координат. Если число неизвестных параметров { } равно n, то, выбрав n линейно независимых функций W i , запишем со- ответствующую систему уравнений 141 0 ) ( dV N L W WRdV v i v , (4.78) из которой может быть найдена сеточная функция { }. Описанный процесс называется методом взвешенных невязок, a W i - весовой функцией. Выбор различных W i приводит к различным классиче- ским вариантам МКЭ. Наиболее популярным из них является метод Галер- кина. В этом случае W I =N I , т.е. в качестве весовой функции выбирается функция формы, с помощью которой аппроксимируется решение { }. Этот метод обычно приводит к наилучшим результатам. Отметим один недостаток метода взвешанных невязок. В этом методе дифференциальный оператор L, содержит производные более высоких по- рядков, чем вариационный функционал F. Поэтому необходимо обеспе- чить выполнение условий непрерывности функций формы более высокого порядка. 4.5 Классификация конечных элементов Наиболее очевидная классификация элементов делит их: - на одномерные; - двумерные; - трехмерные. Эти группы могут разделяться в зависимости от того, включают ли уз- ловые переменные только значение функций (лагранжевые элементы) или также и значение производных (эрмитовы элементы). 4.6 Базисные функции конечных элементов и выбор конечного элемента Произвольная пробная функция, определенная на элементе, записыва- ется в линейной форме e e e U N U ˆ , (4.79) где [N]-матрица базисных функций; {U} e -узловой вектор элемента. 142 Для лагранжева элемента в узле имеется только одна степень свободы- значение функции, и, следовательно, уравнение (4.79) может быть записа- но в виде s k s k U U U U N N N N U 2 1 2 1 ˆ , (4.80) где 1,…, K ,…, S-номера узлов, а S-общее количество узлов элемента. Для эрмитова элемента каждая из компонент U k должна записываться как столбец, включающий и производные указанной функции, которые в этом случае также являются узловыми переменными. Если каждый из S узлов элемента e имеет q степеней свободы, то каждая компонента U k в уравнении (4.80) представляет собой столбец kq k k k U U U U 2 1 , (4.81) где k=1,2,…,S и аналогичные базисные функции [N k ] должны записываться как строки kq k k k N N N N 2 1 , (4.82) Таким образом, уравнение (4.80) при соответствующем выборе переменных является общим для случаев лагранжевых эрмитовых эле- ментов. Ранее мы с вами установили, что базисные функции являются функци- ями независимых переменных (X,Y) и узловых координат. Для получения базисных функций любого отдельного элемента применимы два подхода: -в одном из них используются обобщенные координаты; 143 -в другом интерполяционные формулы. 4.6.1 Получение базисных функций в обобщенных координатах Этот подход особенно удобен для простых элементов, использующих полные полиномы низкого порядка. В случае более сложных элементов указанный подход становится неэффективным, и вместо него обычно при- меняют метод интерполяции. Рассмотренный подход ранее был проиллюстрирован на примере тре- угольного-симплекс-элемента. Применим его для прямоугольного элемен- та со сторонами параллельными осям глобальной системы координат OXY (рисунок 4.7). Простейшая пробная функция для элемента содержит только четыре неизвестных параметраL i |