Основы вычислительного теплообмена и гидродинамики - Аникеев А.А., Молчанов А.М., Янышев Д.С., страница 13

,=0∂ u1∂ uN(7.13)Координатные функции fn строятся на основе функций формы каждогоэлемента. Каждая из функций формы равна 1 в «своей» узловой точке и нулюво всех других узлах данного элемента. Вне элемента все его функции формытакже считаются равными нулю. Очевидно, что для каждого элементатребуется количество функций формы, соответствующее количеству узлов вэлементе. Таким образом, функция формы m-ого элемента, равная 1 впринадлежащей ему узловой точке n, является как бы «представителем»координатнойфункции.аппроксимируетсяПоэтомусуммойтемпературноепроизведенийегополефункцийвэлементеформынаприближённые значения температуры в узлах этого элемента.

Очевидно, чтодля каждого элемента аппроксимация своя, но границах при этом должнасохраняться непрерывность поля температуры.1137.3Выбор типа элементов и составление функцийформыВыбор типа элементов является первым этапом в решении задачи спомощью МКЭ. Существует несколько типов конечных элементов. Видыплоских элементов представлены на Рисунок 7.2.Рисунок 7.2 Виды конечных элементов для плоских задач.Как видно из рисунка, чаще всего применяют элементы треугольнойформы(ихещёназывают2-симплексами),четырёхугольникиикриволинейные треугольники.При решении трёхмерных зада конечные элементы могут представлятьсобой тетраэдры (иначе – 3-симлексы), шестигранники, призмы и т.д.

(см.Рисунок 7.3)Рисунок 7.3 Трёхмерные конечные элементы114Каждому из видов элементов соответствует свой особый вид функцийформы. Здесь мы ограничимся лишь рассмотрением линейных плоскихтреугольных элементов (Рисунок 7.2., a). О функциях формы более сложныхэлементов можно прочесть, например, в [22].Название «линейные» для рассматриваемых элементов обусловлено тем,что при наличии на границе элемента лишь трёх узлов можно использоватьфункции формы, зависящие от координат только линейно:F (x, y ) = a + bx + cy(7.14)Более сложные элементы (Рисунок 7.2., б, г) позволяют использоватьболее сложные законы распределения температуры в элементе.

Например длятреугольного элемента с 6 узлами (Рисунок 7.2., б) функция формы имеет видквадратного полинома:F (x, y ) = a + bx + cy + dxy + ex 2 + fy 2(7.15)Итак, рассчитаем коэффициенты в формуле (7.14).Рассмотрим единичный треугольный элемент, обозначив его узлы как i, j,k.Рисунок 7.4 Треугольный элемент115Расчет длин сторон треугольных линейных элементов производится последующим формулам:L jk =(xi − x j )2 + (yi − y j )2 (x j − xk )2 + (y j − yk )2 Lki =(xk − xi )2 + ( y k − yi )2Lij =(7.16)Площадь рассчитывается следующим образом:S=x j y k − x k y j + y j xi − y k xi + x k y i − x j y i2(7.17)Данные формулы выводятся в курсе аналитической геометрии (см.например [7]).Исходя из приведённых выше условий, по которым функция формыравна единице в «своём» узле и нулю во всех других узлах, можно рассчитатькоэффициенты в формуле (7.14):x j y k − xk y j ai =2Sy− ykj b =; i2Sx − xj ci = k2Sx k y i − xi y ka j =2Sy − yi bj = k;2Sxi − xk c j = 2Sxi y j − x j yiak =2Syi b = − yj k2Sx j − xi ck =2S(7.18)Необходимо отметить, что коэффициенты ai , a j , ak как правило, врасчете не участвуют.В п.

7.2 говорилось, что температура в элементе может бытьпредставлена как произведение функций формы на значение температуры вузлах:u (m ) ( x, y ) = ui Fi(m ) ( x, y ) + u j F j(m ) ( x, y ) + u k Fk(m ) (x, y )116(7.19)Исходя из этого выражения и формулы (7.14), запишем выражение дляпроизводных от температуры по координатам:∂u= bi ui + b j u j + bk u k∂x∂u= ci ui + c j u j + ck u k∂y(7.20)Как видно из этого выражения, значение производных по координатам, а,следовательно, и градиента температуры в каждом элементе постоянно.Теперь приведём некоторые свойства функций формы, которые будутиспользованы нами в дальнейшем:Интеграл от любой функции формы по площади элемента равен 1/3площади элемента:∫∫Fn(m ) (x, y )dxdy =SmS (m ); n = i, j , k3(7.21)Интегралы по стороне Lij вычисляются по следующим формулам:∫ F (x, y )dl =iLij7.4Lij2;∫Fi2(x, y )dl =LijLij3∫ F (x, y )F (x, y )dl =;ijLijLij6(7.22)Система уравнений МКЭИсходя из выбранного способа представления температуры (7.10),функционал (7.9) можно представить как сумму функционалов по всемэлементам:MI=∑I(m)(7.23)m =1I (m ) =2  ∂ u 2 ∂u λ x  dxdy + αu 2 − 2q wu dl+λ−2QuyV∂ y  ∂ x (m)S L( m )∫(∫∫117)(7.24)Второй интеграл в (7.24) вычисляется только для граничных элементов.Окончательное условие минимума полученного функционала с учётом(7.13) может быть представлено в виде:∂∂ un M (m )  M ∂ I (m )I == 0, n = 1, ...

, N m =1 m =1 ∂ u n∑∑(7.25)Заметим, что значение функционала I (m ) в любом элементе зависиттолько от температуры в узлах этого элемента. Соответственно в выражении(7.25) отличны от нуля будут производные от I (m ) по u n лишь техэлементов, которые включают в себя узел n.Это обстоятельство позволяет подойти к проблеме формированиясистемы уравнений МКЭ двумя способами.Первый способ заключается в переборе узлов. Поочерёдно выбираетсякаждый узел и определяется, в какие элементы он входит.

Затемзаписывается и приравнивается нулю сумма частных производных потемпературе в данном узле от функционалов соответствующих элементов.Так одно за другим формируются N уравнений системы МКЭ.Второй способ основан на переборе элементов. Рассматривается каждыйэлемент, определяется, какие узлы в него входят, после чего производные отфункционала этого элемента по температурам в данных узлах заносятся всоответствующие уравнения системы МКЭ.Второй способ получил более широкое распространение и легчепрограммируется, хотя при этом он менее нагляден. В дальнейшем мы нанесложном примере рассмотрим оба способа.Теперь получим явные выражения для производной от функционала потемпературе в узле i.11822 ∂u ∂ I (m )∂    ∂u λ x  + λ y  − 2QV u  dxdy + αu 2 − 2q wu dl  ==∂ ui∂ ui  (m )  ∂ x ∂yL( m )S 22 ∂ λ x  ∂ u   + ∂ λ y  ∂ u   − ∂ (2QV u ) dxdy += ∂ ui   ∂ x   ∂ ui   ∂ y   ∂ uiS ( m )∫(∫∫)∫∫(7.26) ∂∂∫  ∂ u (αu ) − ∂ u (2q u )dl2+wL( m )iiЗдесь мы воспользовались тем фактом, что производную от интеграламожно взять, продифференцировав подынтегральное выражение.Воспользовавшись формулами (7.20), получим значения для каждойпроизводной в выражении (7.26):∂∂ ui  ∂ u 2 λ x   = 2 b i ui + b j u j + bk u k bi  ∂ x  ()2∂   ∂u  λ y  = 2 c i ui + c j u j + ck u k ci∂ ui   ∂ y  ∂(2QV u ) = 2QV Fi (x, y )∂ ui()(7.27)∂(2qwu ) = 2qw Fi∂ ui( )∂αu 2 = 2α Fi ui + F j u j + Fk u k Fi∂ ui()Подставляя (7.27) в (7.26) и используя формулы (7.21) и (7.22), послепреобразований получим следующие выражения для интегралов, входящих внего:1192  ∂ u 2 ∂u λ x  dxdy =+−2QuλyV∂ y  ∂ x (m)S ∫∫(7.28)Q = 2S (n ) λ x bi bi ui + bi b j u j + bi bk u k + λ y ci ci ui + ci c j u j + ci ck u k − V 3 ()∫ (αu2L( m )() αu αu j q w − 2q wu dl = 2 Lij  i +−62  3)(7.29)Таким образом, формула (7.26) приобретает вид:Q∂ I (m )= 2 S (n ) λ x bi bi ui + bi b j u j + bi bk u k + λ y ci ci ui + ci c j u j + ci ck u k − V∂ ui3()()+ αu α u j q w + 2 Lij  i +−62  3(7.30)Для формирования системы уравнений удобно записать полученныевыше соотношения в матричной форме.

Для получения матричной записипринято использовать так называемую локальную нумерацию узлов исоответствующих им температур.Рассмотримтреугольныйэлементномерp,имеющийточкисглобальными номерами i, j, k, и будем условно считать узел i – первым, j –вторым, k – третьим:ui = u1( p ) ; u j = u 2( p ) ; u k = u3( p )(7.31)Для принятой системы нумерации выражение для частных производныхот функционала элемента примет вид:∂ I ( p)= g ( p )U ( p ) − φ ( p )∂uЗдесь под∂ I ( p)понимается вектор-столбец вида:∂u120(7.32) ∂ I ( p)   ∂ I ( p)   ∂ ui   ∂ u1( p )  ∂ I ( p)  ∂ I ( p)   ∂ I ( p) == ∂ u j   ∂ u 2( p ) ∂u ( p)   ( p) ∂I  ∂I  ∂ u   ∂ u( p) k  3 (7.33)g ( p ) – локальная матрица теплопроводности. g iip)(g =  g ji g ki( p)g ik   g11( p)g jk  =  g 21  ( p )g kk g 31g ijg jjg kj( p)g12( p)g 22( p)g 32( p) g13( p) g 23( p ) g 33(7.34)В общем виде, при отсутствии на гранях элемента какого-либо тепловогопотока её можно представить в виде: bi big ( p ) = λ x ⋅ S ( p ) ⋅  b j bi bk bibi b jb jb jbk b j ci cibi bk (p)b j bk  + λ y ⋅ S ⋅  c j cibk bk  ck cici c jc jc jck c jci ck c j ck ck ck (7.35)φ ( p ) – вектор тепловой нагрузки.( p) φi   φ1 φ ( p ) =  φ j  =  φ 2( p )  φ   φ ( p )  k  3 Вектор тепловой нагрузки включает в себя тепловые воздействия,оказываемыенатело(тепловыепотокинагранице,внутреннеетепловыделение и т.п.).Если бы тело состояло всего лишь из одного элемента, то вся системауравнений метода конечных элементов могла бы быть записана в виде:g (1)U (1) = φ (1)121(7.37)В реальных случаях такое крупное разбиение,естественно, неприемлемо.