Курс лекций - Математическое моделирование технических объектов, страница 9

В матричной форме система имеет (9.14) вид:{Ф} = [C] {}(9.15)где:{Ф}T = [Фi Фj Фk Фq];{}T = [i j k q];(9.16)1XiYiZi1XjYjZj = [C](19.7)1 Xk Yk Zk1 Xq Yq ZqСтрока коэффициентов  в (9.16) может быть получена обращением матрицы [C] [C]–1 с последующим умножением (9.15) на [C]–1.{} = [C]–1 [ Ф ](9.18)Интерполяционный полином (9.13) в матричной форме имеет вид:1 = 1 + 2 x + 3 y + 4 z = [ 1 x y z] 224Поэтому с учетом (9.18) имеем: = [ 1 x y z ] [C]–1 [ Ф ](9.19)Определитель матрицы [C] равен шести объемам тетраэдра.Пример 9.4.

Определить ФФ, используя процедуру обращения матрицы длясимплекс – элемента на рисунке 9.5.36Решение. По значениям координат узлов составим матрицу [C]ветствующую ей обратную матрицу [C] –1 :1 1 2 10 601 0 0 0 =[C]0 -3 3=[C]-1 = 11 2 0 063 -1 -11 1 0 30 -1 -1(слева) и соот00-12Для определения ФФ воспользуемся матричным представлением интерполяционного полинома (9.6), согласно которому  = [N] {Ф}, откуда, учитывая выражение (9.19), имеем:[N] = [ 1 x y z ] [C]–1то есть:[N] =6-3-1-103-1-1(6 – 3x – y – z );1616 [1 x y z]003000-12или:[N] =[y2;16(6 – 3x – y – z );16(– y + 2z )]Таким образом, ФФ рассматриваемого элемента имеют вид:y1N1 =;N2 =(6 – 3x – y – z );26N3 =1(3x – y – z );6N4 =16(– y + 2z )5.3.

Интерполяционные полиномы.При обсуждении в предыдущем разделе интерполяционных соотношений для отдельных конечных элементов (симплекс – элементов) мы не фиксировали числовые значения узловых координат и ориентацию КЭ, выбирая их так, как было удобно для изложения сутипроблемы. Подобный «произвол» весьма является важным достоинством метода КЭ, поскольку свобода выбора размеров (при выборе узловых координат) и ориентации КЭ позволяет составлять самые общие вычислительные алгоритмы и подпрограммы, моделирующиеповедение отдельных КЭ. Обычно эти подпрограммы составляют основу библиотек конечных элементов в САПР теплового, прочностного и других видов анализа конструкции.

Указанные подпрограммы могут быть далее использованы без изменения при рассмотрении областей с самыми разнообразными границами. Поскольку при решении задачи анализа поведения ИТО в заданной области производится ее дискретизации (разбиение на КЭ), то будемназывать указанную область – дискретизированной областью (D-область).Перейдем к выводу системы уравнений для области в целом.

Другими словами, будемрешать задачу включения каждого элемента в заданную область. Эта задача требует решитьсначала проблему выбора и преобразования систем координат, в которых заданы конечныеэлементы и сама D-область.37Систему координат, связанную с элементом, будем называть местной, а систему координат, в которой задана D-область – глобальной.Местная система координат. Получение системы уравнений для узловых значений неизвестных величин включает интегрирование по площади элемента функций формы или ихчастных производных. Введение местной системы координат может существенно упроститьпроцесс интегрирования.Рассмотрим механизм преобразования интерполяционных соотношений, записанных вглобальной системе координат, к виду, представляющему эти соотношения в глобальной системе координат.

С этой целью рассмотрим треугольный элемент, в котором скалярная величина представляется согласно (9.10) как:  = Ni Фi + Nj Фj + Nk Фk, а ФФ определяются по(9.11). Поместим местную систему координат в центре элемента, имеющего координаты (Xc,Yc), где:XC =(Xi + Xj + Xk)3иYC =(Yi + Yj + Yk)3(9.20)Обозначив через s (t) – абсциссу (ординату) местной системы координат, запишемформулы преобразования координат (рисунок 9.6):x = Xc + s;y= Xc + t(9.21)ФФ Ni в глобальной системе координат, как было установлено ранее, имеет вид:Ni = 0,5 А –1 [Ai + Bi x + Ci y]Рис. 9.6Подставляя сюда вместо x и y их выражения через s и t, получим:Ni = 0,5 А –1 [Ai + Bi (Xc+s) + Ci (Yc+t)]или:Ni = 0,5 А –1 [(Ai + Bi Xc+ Ci Yc) + Bi s + Ci t](9.22)В результате преобразования Bi и Ci остаются неизменными и по-прежнему умножаются на независимые переменные. Константа же Ai – изменяется.

С учетом (9.10 и 9.20) имеем:(Ai + Bi Xc+ Ci Yc) = 2А/3Таким образом, ФФ в местной системе координат равна:Ni = 0,5 А –1 [(2А/3) + ( Yj –Yk ) s + ( Xk – Xj ) t ]Аналогично получим выражения для других функций формы:Nj = 0,5 А –1 [(2А/3) + ( Yi –Yk ) s + ( Xk – Xi ) t ]Nk = 0,5 А –1 [(2А/3) + ( Yi –Yj ) s + ( Xj – Xi ) t ]38Известно, что интеграл от функции, заданной в глобальной системе, может быть вычислен в местной системе координат с помощью соотношения:Rf(x,y)dxdy=R^f{x[s,t],y[s,t]J}dsdtгде: R и R^ – старая и новая области интегрирования, J - отношение площадей в двухсистемах (J=Аxy/Ast).

Форму элементов при переходе от местной системы координат к глобальной оставляют без изменений, поэтому R = R^, кроме того, местная система и глобальная система координат являются декартовыми, поэтому J=1. Следовательно, в нашем случае:R^f(x,y)dxdy=R^f{x[s,t],y[s,t]}dsdt(9.23)Рассмотрим теперь задачу включения каждого конечного элемента, заданного в местнойсистеме координат, в рассматриваемую D – область. Для этого необходимо выразить интерполяционные уравнения для каждого КЭ, используемого в ансамбле, через глобальные координаты и глобальные узловые значения.С этой целью рассмотрим показанную на рисунке 9.7 пятиэлементную конфигурациюконечных двумерных симплекс – элементов, покрывающую некоторую D – область.

Координаты всех узлов считаются известными. Узлы перенумерованы от 1 до 6, а порядковый номер конечного элемента указан на рисунке в скобках внутри соответствующего конечногоэлемента.Условия Ф, выполняющиеся в узлах ( = 1,2, … , 6), представляют собой глобальныестепени свободы. Для составления интерполяционных полиномов, действующих внутри каждого конечного элемента, отметим символом «звездочка» один из узлов внутри каждого конечного элемента. Тогда, по аналогии с полученным ранее выражением (9.10), для элемента[1] можно записать:[1] = N2[1]Ф2 + N3[1]Ф3 + N1[1]Ф1(9.24)Здесь выражения Nqпредставляют ФФ конечного элемента [1] в q–м узле (q=1,2,3). Дляих правильного вычисления в формулы (9.11) и ( 9.9) следует подставлять значения глобальных координат. Только в этом случаевыражение (9.24) действительно позволит учесть соответствие индексов элемента [ i, j, k]глобальным номерам узлов.

Согласно принятой на рисунке (9.7) нумерации узлов Э q, соответствие [i j k] примет вид:[1]Э1: [231], Э2: [324], Э3: [534], Э4: [635], Э5: [136](9.25)Учитывая принятую нумерацию индексов, приходим к следующей совокупности уравнений для всех конечных элементов D – области:39Рис. 9.7[1][2][3][4][5]=====N2[1]Ф2N3[2]Ф3N5[3]Ф5N6[4]Ф6N1[5]Ф1+++++N3[1]Ф3N2[2]Ф2N3[3]Ф3N3[4]Ф3N3[5]Ф3+++++N1[1]Ф1N4[2]Ф4N4[3]Ф4N5[4]Ф5N6[5]Ф6Пример 10.1.Получить ФФ N6[4] и N6[5] в заданной на рисунке 9.7 D – области.Решение.1.

В соответствии с выражениями (9.11) и (9.9) запишем общее выражение для ФФ в произвольной точке элемента 4:Ni =Ai+Bix+Ciy2A={ (XjYk – XkYj + (Yj – Yk) x + (Xk – Xj) y }2A2. Пользуясь выражением (9.25), запишем соответствие индексов для конечного элементаЭ4 имеем: [i=6j=3k=5]. Заменяя индексы их значениями, получим для N6[4] :N6[4] ={ (X3Y5 – X5Y3 + (Y5 – Y3) x + (X3– X5) y }2A3. Аналогичное соответствие индексов для конечного элемента Э5 имеем: [i=1j=3k=6].Заменяя индексы их значениями, получим для N6[5]:N6[5] ={ (X3Y5 – X5Y3 + (Y3 – Y5) x + (X5– X3) y }2AПоследние две формулы показывают, что ФФ N6[5] и N6[4] – совершенно разныефункции, аппроксимирующие заданный функционал соответственно в пределах конечныхэлементов 5 и 4.

Однако, в самом шестом узле обе эти функции принимают единичные значения, поскольку числители обоих формул в этой точке принимают значения определителя(9.8–а), равные 2А.406. Решение краевых задач методом конечных элементов.До настоящего времени мы рассмотрели: вопросы аппроксимации непрерывной функции на отдельном элементе и методику получения множества кусочно-непрерывных функций(КНФ), аппроксимирующих данную непрерывную функцию в D–области. Это множествоКНФ определяется числовыми значениями узловых величин. Однако остался открытым вопрос получения для узловых величин таких числовых значений, при которых множествоКНФ, определенных для конечных элементов, аппроксимирует с заданной точностью интересующий исследователя физический параметр ИТО. Рассмотрим порядок получения системы уравнений, решение которых позволит это сделать.6.1.

Задача переноса тепла в стержне.Постановка задачи. Выберем в качестве ИТО одномерный стержень с коэффициентом теплопроводности , показанный на рисунке 11.1-а. Стержень имеет теплоизолированную боковую поверхность. К левому концу стержня подводится тепловой поток заданной интенсивности q (Вт/см2). На правом конце стержня происходит конвективный обмен тепла скоэффициентом теплообмена – h (Вт/см2 оС).

Температура окружающей среды – Тос (оС). Поскольку стержень теплоизолирован, потерь тепла через боковую поверхность не происходит.Требуется определить температурное поле вдоль стержня в установившемся режиме.Известно, что для данной модели распределение температуры внутри стержня описывает следующее дифференциальное уравнение:д2 Tдx2=0(11.1)б)a)Рис. 10.1При этом, поскольку в установившемся режиме в точках приложения (при х=0) и отвода (х=L) тепла тепловая энергия не должна «задерживаться», должны быть соблюденыследующие граничные условия:– на левом конце стержня (х=0):–дTдx+q =0(11.2)на правом конце стержня (х=L):дTдx+ h (T – TОС) = 0(11.3)41Если тепло отводится от стержня, тепловой поток q должен быть положителен, в противном случае – отрицателен.Исследования методами вариационного исчисления показывают, что с математической точки зрения в интересующем нас установившемся режиме должен достигать минимума следующий функционал:=V22dV[ дT]дxS [ QT ++h(T – TOC)22] dS(11.4)Учитывая, что боковая поверхность стержня теплоизолирована, приведенный функционал можно представить в следующем виде:=V22dV + S1[ дT]дx(qT ) dS +S2h(T – TOC)2 dS2(11.5)С физической точки зрения функционал (11.5) моделирует непрерывность тепловогопотока в установившемся тепловом режиме.

Это означает, что в любой момент временисумма подводимой (через поверхность S1) к стержню и рассеиваемой им (через поверхностьS2) тепловой энергии равна энергии, сосредоточенной в объеме (V) стержня. В противномслучае, не отводимый от стержня избыток тепловой энергии будет продолжать нагреватьстержень, что противоречит условию установившегося режима.Поскольку, с одной стороны, установившийся режим описывается дифференциальным уравнением (11.1) с граничными условиями (11.2 и 11.3), а, с другой стороны, функционал (11.4) достигает минимума именно в установившемся режиме, то минимум функционала (11.4) и является решением ДУ (11.1) с граничными условиями (11.3).Температура стержня во всех точках сечения S1 (S2) одинакова и равна неизвестной пока(но постоянной в стационарном режиме) величине – Т1 (Т3).