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

Учитывая, что в данном случае S1= S2 = A и в силу сказанного, выражение (11.5) принимает вид:qT1S1dS +S2h(T3 – TOC)22dS = qT1А +h2(T3 – TOC)2 А(11.6)Таким образом, исходное уравнение для определения температуры в каждой точкестержня методом МКЭ примет вид:=V2[дTдx]2dV + qT1А +h2(T3 – TOC)2 А(11.7)Реализация метода МКЭ включает этапы:1. Определение подобластей (конечных элементов) и их узловых точек. В данномслучае, стержень может быть разбит на два одномерных симплекс – элемента, как это показано на рисунке (10.1-б) с узловыми значениями Т 1, Т2 и Т3.

Температура внутри элементовнаходится из формул:T[1] = N1[1] T1 + N2[1] T2 ;T[2] = N2[2] T2 + N3[2] T3 ;ФФ здесь согласно (9.5) равны:N1[1]=N2[2]=(X2 – x)L[1](X3 – x); N2[1]=; N3[2]=(x – X1);L[1](x – X2)(11.8)42L[2]L[2]2. Вычисление частных производных, входящих в выражение (11.7):дT[1]дx=дT[2]дx1(-T1+T2);L[1]=1L[2](-T2+T3)(11.9)3. Разделение интеграла в выражении (11.7) на два (по числу подобластей – конечных элементов, выделенных в пункте 1). Необходимость разбиения интеграла продиктованатем, что производная температуры по переменной х (градиент температуры по оси ОХ),входящая под знак интеграла, не является непрерывной в точке Т 3.

Учитывая, что dV=Adx,где А – площадь сечения стержня (А1 = А2 = А3 =А), после разделения и подстановки пределов интегрирования получаем выражение:x22дTдx[ A[1]]2dV =[[1]2Vx3дT 2] dx +  Aдx2[2][[2]x1дT 2] dx (11.10)дxx24., Проведение подстановки (11.9) в (11.10) и интегрирование:V2[дTдx]2dV =[1]A[1]2L[1][-T1+T2]2+[2]A[2][-T2+T3]22L[2](11.11)5. Выражаем функционал через узловые значения температуры, для чего объединяем выражения (11.7) и (11.11):[3][1][2](11.12)(-T1+T2)2 + C(-T2+T3)2 +qA[1]T1 + hA(-T3+TOC)2= C222Здесь приняты следующие обозначения:С(1) = (А(1)(1)/L(1));С(2) = (А(2)(2)/L(2))6. Получение системы алгебраических уравнений.

Правильными значениями Т1,Т2 и Т3 являются те, при которых величина функционала  достигает минимума.Приравнивая нулю первую производную функционала (11.12) по Т 1, получаем первоеуравнение системы:д(11.13)= C[1] T1 - C[1] T2 + qA[1] = 0дT1Аналогично получаем еще два уравнения:ддT2= -C[1] T1 + [C[1] +C[2] ]T2 -C[2] T3 = 0ддT3= -C[2] T2 + [C[3] +hA3 ]T3 - hA3TOC = 0(11.14)Запишем полученную систему в матричной форме:С(1)-С(1)0Т1-qA143-С(1) С(1)+С(2)-С(2)0-С(2)С(2)+hA3Т2Т3=0hA3TOC(11.15)В более общей матричной форме система примет вид:C  T = F(11.16)Матрица C в формуле (11.16) называется «глобальной матрицей жесткости».

В контексте задачи переноса тепла –это – «глобальная матрица теплопроводности». Векторстолбец F называется «глобальным вектором нагрузки». Искомый вектор [T] будем называть вектором решения.Пример 11.1. Рассчитать температурное поле в круглом стержне с площадью поперечного сечения A=1 см2 и длиной L=7,5 см с теплоизолированными стенками. К левому концустержня подводится тепловой поток q = 150 Вт/см2. Коэффициент теплопроводности материала стержня и коэффициент конвективного теплообмена на правом конце стержня соответственно равны: =75 Вт/(см  ОС), h = 10 Вт/(см2  ОС).

Температура окружающей среды равна ТОС=40 ОС.Решение.1. Тепло подводится к стержню, поэтому тепловой поток q следует записывать со знаком«минус»: q = - 150 Вт/см2.2. Рассчитываем значение термов, входящих в коэффициенты матриц C и F:С(1) =(А(1)(1)/L(1))=(175/3,75)=20Вт/(смОС),С(2) =(А(2)(2)/L(2))=(175/3,75)=20Вт/(смОС),hA3=10Вт/(смОС), -qA1= -(-150)1 = 150Вт/см,hA3TOC=10140 = 400Вт/см.3. Окончательная система уравнений примет вид:20-200-2040-200-2030Т1Т2Т3=15004004. Решением полученной системы являются следующие узловые значения температуры:Т1=70 оС, Т2=62,5 оС ; Т3=55 оС.Проблема реализации МКЭ на ЭВМ.

Процедура минимизации  приводит к системеуравнений, которые решаются относительно узловых значений температур. Однако, с точкизрения реализации процедуры минимизации  на ЭВМ, целесообразно функционал (11.4)представить в виде суммы вида:m = 1 + 2 +…+m =  1i =1(11.17)где: m – количество конечных элементов, на которые разбивается ИТО.Дело в том, что в библиотеке САПР, реализующей минимизацию функционала  наЭВМ, содержаться модели не всего ИТО, а именно конечных элементов (например, симплекс– элементов), причем мощность указанной библиотеки КЭ и определяет функциональные44возможности САПР ИТО.

В процессе решения задачи ЭВМ (в соответствии с заданием напроектирование) автоматически объединяет модели конечных элементов в единую модельИТО. В этой связи, представляется целесообразным описать последовательность шагов получения системы линейных уравнений (11.16), используя в качестве исходного шага разбиение(11.17).

Тем более, что эта процедура и является центральной в работе инженера, моделирующего поведение ИТО на ЭВМ.Из примера (11.1) ясно, что функционалы по отдельным конечным элементам, выраженные через узловые значения, имеют вид:2(L[1])21 =V[1]2(L[2])22 =V(-T1+T2)2dV +S(11.18)[1](-T2+T3)2dV +[2]qT1 dSSh(T3+TOC)2dS2[2]Проведем дифференцирование (1) системы (11.18) по всем узловым значениям:д(1)дT1(L[1])2=(-T1+T2) (-1)dV +V[1]д(1)дT2=д(1)дT3=0S[1](L[1])2 q dS(-T1+T2) dVV[2]Вычисляя в этой системе интегралы, и применяя обозначения, принятые в формуле (11.12),получим следующую систему уравнений в обычной и матричной форме:д[1]= + C[1] T1 - C[1] T2 + qA [1]дT1д[1]= - C[1] T1 + C[1] T2 + 0дT2д[1]=0+0+0дT3д[1]дT1д[1]дT2д[1]дT3=C[1]-C[1]0=-C[1]C[1]0=000qA[1]T1T2T3+0045Для краткости изложения будем далее обозначать ее так:д[1]д[T]Запишем систему уравнений (11.19) в матричной форме для первого КЭ:д(1)д[T]=[ C (1) ]  [ T ] +[ F ](11.19)В отличие от системы уравнений (11.16) в системе (11.19) матрица коэффициентов Cназывается «матрицей жесткости элемента».

Ее название в контекстезадачи переноса тепла - «матрица теплопроводности элемента». Вектор-столбец Fкак и ранее является «глобальным вектором нагрузки».Проведем теперь дифференцирование второй компоненты (2) системы (11.18)по всем узловым значениям:(1)д(2)дT1=0д(2)дT2=д(2)дT3=(L[2])2(-T2+T3)( -1) dV(L[2])2(-T2+T3) dV +V[2]V[2]h (T3-TOC) dSS[2]После вычисления интегралов получим систему уравнений:д (2)дТ1д (2)дТ2д (2)дТ3=0+0+0=0+С(2)Т2-С(2)Т3=0-С(2)Т2+(С(2)+hA3)Т3Или в матричной форме: (2) Т1 (2) Т2 (2) Т3=000=0С(2)-С(2)=0-С(2)(С(2)+hA3)Т1Т2Т30+0+hA3Учитывая аддитивный характер функционала  можно утверждать, что для его минимизации поузловым значениям необходимо, чтобы выполнялось равенство:46дд[T]=д[1]д[T]+д[2]д[T]=0(11.20)Поэтому, складывая выражения для обеих компонент функционала в матричном виде,получим окончательную систему уравнений:С(1)-С(1)0Т1qA10(1)(1)(2)(2)-СС +С-СТ2+00=(2)(2)0-СС +hA3Т3-hA3TOC0Данная система идентична системе (11.15).

Таким образом показано, что система уравненийдля минимизации исходного функционала  может быть получена путем суммирования соответствующих матриц для элементов.В этой связи представляется актуальной проблема формального получения матриц теплопроводности и нагрузки для отдельных конечных элементов на основании анализа их физических и геометрических параметров, а также вопросы формального получения на их основе и решения глобальной системы уравнений, аппроксимирующей поведение ИТО. Совершенно очевидно, что такой подход позволяет выбирать характеристики элементов, наиболееприемлемые для каждой конкретной задачи.12. Уравнения метода конечных элементов: Задача переноса тепла.Вернемся к анализу функционала  (11.5), моделирующего непрерывность тепловогопотока через стержень в установившемся тепловом режиме. Пусть стержень разбивается наЕ симплекс–элементов. В пределах отдельного (e-го) элемента величины (е), q(е), h(е) считаем заданными и постоянными, а величины узловых температур Т i(е) и Тj(е) подлежат определению.

Для минимизации  по аналогии с выражением (11.20) потребуем выполнениясоотношения:Едд[T]д=Ед[T]де=1 д[T][e] =е=1[e]=0(12.1)где [e] – элементарный функционал, представляющий собой сумму интегралов для произвольного конечного элемента (например, для 1-го КЭ имеем: [e] =[1] и так далее). В связис этим, учитывая полученное выше выражение (11.5), имеем:E={ е=12[дT[e]дxV[e]]2 dV + (qT[e]) dS +S1[e]+h(T[e] – TOC)2 dS2}(12.2)S2[e]Для вычисления частных производных элементарных функционалов [e] в формуле(12.1) выразим интегралы в (12.2) через узловые значения температур.47Учитывая соотношение (9.6), запишем интерполяционную формулу для произвольногосимплекс – элемента в общем виде:Т (е) = Ni(е) Ti + Nj(е) Tj = [N(е)] {Т}(12.3)Вычислим далее значение частной производной Т(е) по координате х:дT[e]дxдNi [e]дx=Введем обозначение:[e]Bi=Ti +дNj [e]дxTj(12.4)дNi [e]дxдT[e] = [B[e]]{T}(12.5)дxЭто позволяет получить выражение для функционала (е) в матричной форме.

Согласно(12.2), (12.3) и (12.5) имеем:и запишем (12.4) в матричной форме:[e] = [e][B ]{T}[B[e]]{T}dV +2V[e]+q [N[e]]{T} dS +S1[e](h(TOC)2[e][e][e][N ] {T} [N ] {T} – hTOC[N ]{T} +22) dS(12.6)S2[e]Для вычисления искомых производных, в соответствии с исходной формулой (12.1),покажем предварительно, что в матричном виде:д( [B[e]] {T} [B[e]] {T} )дx= 2 B i [e] [B[e]] {T}(12.7)Действительно, левая часть приведенного тождества представляет собой искомую частную производную от квадрата частной производной Т(е) по Ti , представленную в матричной форме, которая по определению производной от сложной функции и с учетом (12.5) равна:дT[e] )2дT[e] )((()дддxдT[e]дx=2= 2 Bi[e] ( Bi[e]Ti + Bj[e]Tj )дTiдxдTiОткуда, переходя к матричной форме, получаем выражение (12.7).Итак, мы подготовили все необходимое для вычисления и представления в матричнойформе искомой системы уравнений (12.1).

Вычисления проведем в два этапа: на первом получим матрицы для конечных элементов, а на втором – объединим их в матрицы ИТО.Первый этап состоит в вычислении частных производных от элементарного функционала [e] (12.7) по всем узловым значениям температуры. Последовательно находим для конечного элемента 1.д[e]дT1=V[e]B1[e][B[e]]{T}dV+S1[e]q N1[e] dS +48+S2hN1[e][N[e]] {T} dS -[e]S2– hTOCN1[e] dS(12.6)[e]Вектор {T} не зависит от переменных интегрирования, поэтому, объединяя первое итретье слагаемое, вынося этот вектор за скобки и вычисляя производные для остальных узловых переменных конечного элемента 1, приходим к системе:В данной системе выделены элементы, представляющие собой транспонированные матрицы [В(е)]T и [N(e)]T.