Курс лекций - Математическое моделирование технических объектов (1075784), страница 8
Текст из файла (страница 8)
Однако условие полной непрерывности между элементаминалагается только в узлах, а не на всем протяжении границ раздела. Диапазон применимостиМКЭ, его высокая эффективность и сравнительная простота, с которой могут быть учтеныреальные граничные условия, делают МКЭ серьезным соперником для любого конкурирующего метода.К недостаткам метода следует отнести: (а) в основе МКЭ лежит дискретизация всегоИТО, что неизбежно ведет к очень большому количеству конечных элементов (особенно втрехмерных задачах); (б) МКЭ часто приводит к нереальным разрывам значений физическихвеличин между смежными элементами.При решении задач методом КЭ используются одномерные, 2 и 3-мерные КЭ.
Одномерный КЭ показан на рисунке 9.1. Площадь поперечного сечения одномерного КЭ можетизменяться по длине, но в большинстве практических задач ее считают постоянной. Наиболее часто такой элемент используется в одномерных задачах распространения тепла и в задачах строительной механики.Рис. 9.1Для построения дискретной модели двумерной области используют треугольники и 4хугольники, которые могут иметь как прямолинейные, так и криволинейные стороны.Собственно процесс дискретизации ИТО может быть разделен на два этапа: (а) разбиение ИТО на КЭ и (б) нумерация элементов и узлов.
Последний этап может существенноповлиять на эффективность вычислений.Требование простоты КЭ связано с тем, что при моделировании области должно бытьиспользовано большое число элементов, поэтому деление области на треугольники являетсянаилучшим способом разбиения.МКЭ основан на аппроксимации непрерывной функции (температура, давление, перемещение и др.) дискретной моделью, которая строится на множестве кусочно-непрерывныхфункций, определенных на конечном числе подобластей, называемых элементами.
В качестве функции, действующей внутри границ (и на границах) элемента обычно применяется полином, порядок которого и определяет тип элемента.На практике используются три типа элементов: симплекс-элемент, комплекс-элемент имультиплекс-элемент.31Симплекс - элементам соответствуют полиномы, содержащие константу и линейныечлены. Число коэффициентов в таком полиноме на 1 больше размерности координатногопространства.
Например, полином: = 1 +2 x +3 yпредставляет собой симплекснуюфункцию для двумерного треугольного элемента. Этот полином линеен по X и Y и содержиттри коэффициента, потому что треугольник имеет три узла.Комплекс – элементам соответствуют полиномиальные функции, содержащие константу, линейные члены, а также члены второго, третьего и более высоких порядков, еслиэто необходимо.
Форма комплекс – элементов может быть такой же как у симплекс - элементов, но с дополнительными граничными (и даже внутренними) узлами. Число узлов в комплекс – элементе должно быть больше размерности координатного пространства + 1. Интерполяционный полином для 2-мерного треугольного комплекс – элемента имеет вид: = 1 + 2x + 3y + 4x2 + 5 xy + 6y2(9.2)Это соотношение включает шесть коэффициентов, поэтому рассматриваемый элемент должен иметь шесть узлов.Для мультиплекс – элементов также используются полиномы, содержащие члены высокого порядка, Однако, для обеспечения непрерывности при переходе от одного мультиплекс – элемента к другому границы элементов должны быть параллельны осям координат.Примером мультиплекс – элемента является прямоугольный элемент.5.2.
Функции формы.Рис. 9.2Одномерный симплекс – элемент представляет собой прямолинейный отрезок длины Lс двумя узлами – по одному на каждом конце отрезка (рисунок 9.2). Узлы обозначаются индексами i и j, значения функции в узлах – через Фi иФj соответственно.Начало системы координат располагаетсявне КЭ.
Полиномиальная функция для скалярной величины (например, температуры – Т илидавления – Р) такова: = 1 +2 x(9.3)Коэффициенты 1 и 2 определяются с помощью условий в узловых точках: = Фi при x = Xi и = Фj при x = Xj.Эти узловые условия приводят к системе двух уравнений:Фi = 1 +2 XiФj = 1 +2решение которой дает: 1= (Фi Xj - Фj Xi)/L; 2 = (Фj - Фi )/LПодставляя найденные значения 1 и 2 в формулу (9.3), получим: = (ФiXj-ФjXi)/L +{(Фj-Фi)/L}xДанное уравнение может быть переписано в виде: = [(Xj-x)/L]Фi+[(x-Xi)/L]Фj(9.4)Линейные функции от х в формуле (9.5) называются функциями формы (ФФ) или интерполяционными функциями. Далее эти функции обозначаются через N. Каждая ФФ долж-32на быть снабжена нижним индексом для обозначения узла, к которому она относится.
Произвольную ФФ будем обозначать через N. В формулу (8.5) входят следующие ФФ:Ni =Xj-xL;иx-XiLNj =Используя эти ФФ, запишем выражение (9.5) в матричной форме:Фi = NiФi + NjФj = [N]{Ф} = [Ni Nj] Фj= [Ni Nj] [Фi Фj]Т(9.6)Функция Ni = 1 в узле с номером i и равна нулю в j-м узле. Аналогично функция Nj =1 в узле с номером j и равна нулю в i-м узле. Эти значения характерны для функций формы.Они равны 1 в одном определенном узле и обращаются в 0 в остальных узлах.Пример 9.1. Одномерный симплекс-элемент используется для аппроксимации температуры в стержне. Узлы 1 и 2 имеют координаты 1,5 и 6 см соответственно.
Известно, чтотемпература в узлах 1 и 2 равна 120 и 90 градусов соответственно. Требуется определитьтемпературу в точке х = 4 см и градиент температуры внутри элемента.Решение: Пользуясь выражением (9.5) для одномерного симплекс – элемента, можно записать закон изменения температуры внутри КЭ:(Xj-x)Lt =Ti +(x-Xi)LTjXi=1,5 см;Ti=120oC;X j=6,0 см;oTj=90 C;x=4 см;L = (Xj – Xj ) = 4,5 см.Подставляя данные в формулу для температуры получаем:(1,5 - 4)(4 – 1,5)t=120o +90o = 103,33 oC4,54,5Для градиента температуры имеем:Данные КЭ:dtdx= -TiL+TjL= -120o4,5+90o= -6,67 oC/см4,5Двумерный симплекс – элемент показан на рисунке 9.3 – это треугольник с прямолинейными сторонами и тремя узлами, по одному в каждой вершине. Примем последовательную логическую нумерацию узлов элемента против часовой стрелки, начиная от произвольно выбранного i-го узла.
Узловые значения скалярной величины обозначим через Фi,Фj, Фk, а координатные пары трех узлов - через (Xi, Yi), (Xj, Yj), (Xk, Yk).Рис. 5.3Интерполяционный полином в данном случае примет вид: = 1 +2 x +3 yРис. 5.4(5.7)33В узлах выполняются следующие условия: = Фi при x = Xi и y = Yi = Фj при x = Xj и y = Yj = Фk при x = Xk и y = YkПодстановка их в (9.7) приводят к системе трех уравнений:Фi = 1 + 2 Xi + 3 YiФ j = 1 + 2 Xj + 3 YjФk = 1 + 2 Xk + 3 Yk(9.8)Обозначим площадь симплекс – треугольника буквой А. Можно показать, что определительсистемы (9.8) связан с А (рис. 5.4) соотношением:[XjYk – XkYj + XiYj – XiYk + XkYi – XjYi] =2AРешая систему (9.8) с учетом (9.8) и вводя обозначения:Ai=(Xj Yk – Xk Yj);Aj = (Xk Yi – Xi Yk),Ak =(Xi Yj – Xj Yi),Bi=(Yj – Yk);Bj = (Yk – Yi),Bk = (Yi – Yj),Ci=(Xk – Xj),Cj =(Xi – Xk),Ck = (Xj – Xi),(9.9)получим значения искомых коэффициентов:1 = 0,5 А –1 [ Ai Фi + Aj Фj + Ak Фk ]2 = 0,5 А –1 [ Bi Фi + Bj Фj + Bk Фk ]3 = 0,5 А –1 [ Ci Фi + Cj Фj + Ck Фk ]Подставляя значения 1, 2, 3 в (9.7) и преобразуя получаемые выражения к виду, подобному (9.6), получим выражение для скалярной величины : = Ni Фi + Nj Фj + Nk Фk(9.10)где:Ni =Ai+Bix+Ciy2A;Nj =Aj+Bjx+Cjy2A;Nk =Ak+Bkx+Cky2A(9.11)Значение Ni в i-м узле составит: Ni = 0,5 А –1 [Ai + Bi x + Ci y] == 0,5 А –1 [Xj Yk – Xk Yj + (Yj – Yk) Xi + (Xk – Xj) Yi] == 0,5 А –1 [XjYk – XkYj + XiYj – XiYk + XkYi – XjYi] = 1Непосредственной проверкой можно показать, что в остальных узлах Ni = 0.Из (9.11) видно, что ФФ линейны по x и y, то есть, градиенты этой величины в направлениях Ox и Oy будут постоянны.
Заметим, что:дN= Вдx( = j, j, k)поэтому градиент в направлении оси Ох составит:дФдx=дNiФi +дxдNiФj +дxдNkФk = BiФi + BjФj + BkФkдx(9.12)Поскольку, переменные В и величины Ф начальных условий (при = i, j, k) фиксируются,как только задаются узловые координаты, то частная производная в (9.12) имеет постоянноезначение. Отсюда следует важный вывод: постоянство градиента внутри каждого элемента означает, что необходимо применять очень малые по величине элементы, чтобы аппроксимировать быстро меняющую функцию .34Пример 9.2. Требуется получить соотношение, определяющее элемент, и вычислитьзначение давления в точке В с координатами (2; 1,5), если заданы начальные значения: Pi =40 H/см2, Pj = 34 H/см2, Pk = 46 H/см2.Давление р внутри элемента определяется по формуле: р = Ni Рi + Nj Рj + Nk Рkгде ФФ Ni , Nj и Nk определяются по (9.11).Подставляя значения координат узлов в обозначения (9.9) для А , В, С (при = i, j,k), получим значения этих коэффициентов:Ai = (45)–(21,5) = 19;Bi = (0,5–5) = – 4,5 ;Ci = (2–4) = – 2 ;Aj = (20) –(05) = 0;Bj = (5 – 0) = 5;Cj = (0 – 2) = – 2;Рис.
9.4Вычисляем определитель:1 Xi Yi2A= 1 Xj Yj1 Xk YkAk =(00,5)–(40) = 0;Bk =(0 – 0,5) = – 0,5;Ck =(4– 0) = 4;Рис. 9.5=11104200,55=20-1=19После подстановки констант в ФФ выражение для р примет вид:p=[(19–4,5x–2y)Pi + (5x – 2y)Pj + (– 0,5x + 4y) Pk19Значение давления в точке В с координатами (2; 1,5) равно:p=740 +734 +546= 39,37 Н/см219Отметим два полезных свойства треугольного элемента. Во-первых, функция изменяется линейно между двумя любыми узлами. Так как узлы определяют границы элемента, меняется линейно вдоль каждой из трех его сторон.
Отсюда следует второе полезное свойство: любая линия, вдоль которой принимает одинаковые значения, есть прямая линия, пересекающая две стороны элемента. Исключением будет случай, когда во всех узлах значения 35одинаковые. Приведенные два свойства позволяют легко определять линии уровня скалярной величины. Обратимся к предыдущему примеру, чтобы проиллюстрировать эти свойства.Пример 9.3 (продолжение примера 9.2). Требуется определить линию уровня, соответствующую величине давления 42 Н/см2, для примера 9.2.Решение. Искомая линия пересекает стороны ik и kj.
Поскольку давление меняется линейно вдоль каждой из сторон треугольника, можно составить простые соотношения, позволяющие получить координаты точек на указанных сторонах, через которые проходит искомая линия. Для стороны jk имеем:(46 – 42)(46 – 34)=(2 – x)(2 – 4)=(2 – y) x = 2,67 см; y = 3,5 см(5 – 0,5)Аналогично вычислим координаты точки на стороне ik: x = 0,67 см, y = 1,67 см.Трехмерный симплекс – элемент показан на рисунке 9.5 – это тетраэдр, четыре узла которого обозначены индексами i, j, k, q, причем обход узлов i, j, k, q проведен, как и ранее, против часовой стрелки. Запишем интерполяционный полином для тетраэдра: = 1 + 2 x + 3 y + 4 z(5.13)Коэффициенты можно определить, используя следующие 4 условия в узлах:Фi= 1 + 2 Xi + 3 Yi + 4 ZiФj= 1 + 2 Xj + 3 Yj + 4 Zj(9.14)Фk = 1 + 2 Xk + 3 Yk+ 4 Zk Фq= 1 + 2 Xq + 3 Yq+ 4 ZqЭта система может быть решена с помощью правил Крамера и связана с вычислением 5-ти определителей.