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

Для значений во внутренних узлах согласно (15) составляем систему из 4-х уравнений с четырьмя неизвестными:F23 – 2F22 + F21(1/3)2+F32 – 2F22 + F12(1/3)2+ e – 1/3F23– F21- F22 = e – 1/3 (e 1/3 + 1)(2/3)F24 – 2F23 + F23(1/3)2+F33 – 2F23 + F13(1/3)2+ e – 2/3F24– F22- F23 = e – 1/3 (e 2/3 + 1)(2/3)F34 – 2F33 + F32(1/3)2+F43 – 2F33 + F23(1/3)2+ e – 2/3F34– F32- F33 = e – 2/3 (e 2/3 + 1)(2/3)F33 – 2F32 + F31(1/3)2+F42 – 2F32 + F22(1/3)2+ e – 1/3F33– F31- F32 = e – 2/3 (e 2/3 + 1)(2/3)Введем для краткости обозначения: F2,2  x, F2,3  y, F3,2  z, F3,3  u. Подставляяих в полученную систему и проводя предварительные вычисления, получим систему из 4-хследующих линейных уравнений: –37x + 10,1y +9z = –16,5+8,2x –37y +9u = –35,45+9y +8,2z – 37u = –18,7 +9x – 37z+10,1u = –7,42Программа вычисления корней данной системы приведена ниже:Система имеет следующее решение:F2,2  x=1.0119F2,3  y=1,4288 F3,2  z=0,7233 F3,3  u=1,0133.uses crt; const n=4;type qw=array[1..n,1..n] of real; Linia=array[1..n]of real;const MotL:qw=((-37,10.1,9,0),(8.2,-37,0,9),(0,9,8.2,-37),(9,0,-37,10.1));BotL:Linia= (-16.5,-35.45,-18.7,-7.42);var m1,m2:qw; x,b1,b2:Linia; aa,ss,zz:real;i,j,k,q,tt:integer;Procedure CoeFA(i,j,k:byte);begin m2[i,j]:=m1[i,j] - m1[i,k]*m1[k,j]/m1[k,k] End;Procedure FreeB(i,k:byte);begin b2[i]:=b1[i]-b1[k]*m1[i,k]/m1[k,k] End;BEGIN clrscr; For i:=1 To n Do x[i]:=0; For i:=1 to n Do b1[i]:=BotL[i];For i:=1 to n Do For j:=1 To n Do m1[i,j]:=MotL[i,j];For i:=1 to n Do For j:=1 to n Do m2[i,j]:=0;For j:=1 to n Do m2[1,j]:=m1[1,j]; b2[1]:=b1[1];For tt:=2 To n DoBegin For i:=tt to n Do m2[i,tt]:=0;For i:=tt To n Do For j:=2 to n Do CoeFA(i,j,tt-1);For i:=tt to n Do FreeB(i,tt-1);For i:=1 to n Do For j:=1 To n Do m1[i,j]:=m2[i,j];For i:=1 to n Do b1[i]:=b2[i] End;zz:=0; X[n]:=b2[n]/m2[n,n];For i:=n-1 DownTo 1 DoBegin zz:=b2[i]; q:=n;For j:=i To n-1 Do Begin zz:=zz-m2[i,q]*x[q]; dec(q); End;x[i]:=zz/m2[i,i] End;For i:=1 to n Do WriteLn('X',i,'=',X[i]:6:4); Repeat Until KeyPressed END.Система имеет следующее решение:F2,2  x=1.0119F2,3  y=1,4288 F3,2  z=0,7233 F3,3  u=1,0133.192.3.

Оценка погрешности дискретной моделиПри разностном решении ДУ в частных производных основным источником ошибокявляются погрешности от замены производных конечными разностями. Эти погрешности называются погрешностями дискретизации. Таким образом, в теории разностных схем основной является проблема наилучшего приближения к ДУ с помощью разностных соотношений,или наилучшей аппроксимации дифференциальных операторов – разностными.Погрешности дискретизации зависят от следующих факторов: способа замены дифференциальных уравнений разностными; от конфигурации элементов конструкции (формы рассматриваемой области); внешних воздействий (граничных условий); длительности рассчитываемого процесса.Определим порядок погрешности дискретизации, который определяется способом замены дифференциальных операторов в задаче – разностными, то есть порядком аппроксимации.

Порядок аппроксимации показывает, каким образом снижаются погрешности с уменьшением шага сетки. Если порядок аппроксимации – первый, то погрешности пропорциональны шагу, если – второй, то – квадрату шага и так далее.Покажем, как определить порядок аппроксимации на примере замены производных конечными разностями. Допустим, что мы хотим заменить первую производную в точке 0 (рисунок 3) идля этого наметим два узла сетки в точках x=-a и x=h-a. Будем считать функцию F и ее производную в точке 0 известными. Воспользовавшись разложением в ряд Тейлора, находим значенияфункции на концах отрезка при x=-a и x=h-a:дFaд2 Fa2д3 Fa2F(–a) = F –++…дX  1!дX2  2!дX3  3!дF(h-a)д2 F(h-a)2F(h–a) = F +++…дX 1!дX2  2!Далее определим значение конечной разности:F(h–a) – F(–a)h=дFдX+д2 FдX2(h-2a)2!+д3 FдX3h2 – 3ah +3a23!+…Погрешность от замены первой производной конечной разностью будет равна:F(h–a) – F(–a)дFд2 F(h-2a)д3 Fh2 – 3ah +3a2+=+…hдXдX2 2!дX3 3!При а=0 разность будет правой, при a=h - левой (9-б и 9-а соответственно).

При этомпогрешности соответственно составят:22д3 F h }{ д F2  h } + { для правой разности:дX2!дX33!23дF  hдFh2 }{} + { для левой разности:дX22!дX33!В том и в другом случае погрешность пропорциональна шагу сетки, то есть имеет место первый порядок аппроксимации производной конечной разностью. Условно это можнозаписать в виде:дFдFFm+1n- FmnFm,n- Fm-1,n=+ O(h)и=+ O(h)hдYhдY20Для вторых разностей ошибка замены второй производной может быть определена аналогично.

Используя разложение функции F в ряд Тейлора вблизи точки X = mh, можно показать, что здесь имеет место второй порядок аппроксимации.3. Задача расчета теплового процесса на дискретной моделиВ электронно-вычислительной аппаратуре имеют место следующие процессы передачи тепла: конвекция, кондукция и лучеиспускание.

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

Блок залит наполнителем с коэффициентом теплопроводности К и удельной теплоемкостью С. Разобьеммысленно конструкцию на части прямоугольной формы, каждую из которых назовем элементомРис. 3Рис. 4Рис. 5Для более высокой точности расчета выберем элементы одинаковых размеров, причемсами размеры элементов примем минимально возможными.

В центре элемента выделим особую точку – узел сетки. Далее, попытаемся определить температуру в каждом узле сетки вкаждый момент времени. Для простоты будем считать, что блок однороден, то есть входящиев него материалы имеют одинаковую теплоемкость и коэффициент теплопроводности. Температуру, определяемую в узле сетки с координатами x, y, z, в момент времени t, обозначимкак tX,Y,Z, а в следующий момент времени как t+1X,Y,Z. Размеры блока, координаты и мощность тепловыделяющих радиоэлементов будем считать заданными.

Кроме того, для решениязадачи должны быть заданы начальные и граничные условия.В начальных условиях задачи необходимо указать температуру во всех узлах сеткиблока в начальный момент времени. Обычно при рассмотрении переходных процессов за начальный момент времени выбирается момент включения электрических цепей под нагрузку.До этого момента температура во всех узлах считается одинаковой и равной наружной температуре, например, комнатной (20ОС или 293ОК).В граничных узлах блока могут быть заданы различные граничные условия. Когда награнице задается значение самой функции, то есть температура – это граничное условие 1-города и решение получается наиболее простым.

Однако, к сожалению, только при грубом упрощении нестационарной задачи (то есть задачи с изменением температуры во времени)можно считать температуру на поверхности бока заданной, например, равной наружной температуре. Наиболее близкими к реальным условиям являются граничные условия 2-го рода,когда задаются плотности теплового потока по всей наружной поверхности блока. Далее мыбудем решать задачу с граничными условиями 2-го рода.3.1. Уравнение передачи тепла через элемент дискретной модели.21Запишем типовые уравнения движения теплоты. Для этого воспользуемся законом сохранения энергии: количество притекающей к данному элементу тепловой энергии равно количеству утекающей энергии плюс количество накапливающейся энергии. В рассматриваемом случае тепловая энергия не превращается в другие виды энергии, однако, другие видыэнергии могут превращаться в тепло.

Например, электрическая энергия целиком превращается в тепло, поэтому в уравнении теплового баланса нужно учесть количество энергии, выделяемой за счет электрических потерь.Рассмотрим прямоугольный элемент объема блока (рисунок 4). Количество энергии,притекающей и утекающей через боковые поверхности этого элемента, выражается через величину плотности тепловых потоков. Удельная плотность теплового потока J [ Дж/м2сек ]определяется количеством теплоты, проходящей через единичную площадь в единицу времени. Чтобы определить количество теплоты, проходящей через боковую грань элемента за некоторое время, необходимо соответствующую плотность теплового потока умножить наплощадь грани и на интервал времени:(JX+ – JX–) hYhZ + (JY+ – JY–) hXhZ + (JZ+ – JZ–) hXhY = C (14)где: J – удельная плотность тепловых потоков,  - время,  - приращение температуры.В правой части уравнения (14) записано количество теплоты, накапливаемой внутриэлемента за время .

Выполним в уравнении (14) следующие преобразования:1. Приведем количество теплоты в левой и правой части уравнения к единичному объему и кединице времени, для этого разделим все члены. Для этого разделим все члены на объемэлемента hXhZhY и на интервал времени .2. Представим приращение температуры  в узле с координатами i, j, k за интервал времени  ввиде разности температур в начале и в конце этого интервала: =  t+1 i, j, k –  t i, j, kВ результате получим уравнение:(JX+ – JX–)(JY+ – JY–)(JZ+ – JZ–)ijkt+1–ijkt++= CУД(15)hXhYhZТеперь в правой части уравнения (15) стоит не теплоемкость элемента, а удельная теплоемкость вещества (наполнителя), составляющего элемент. В целом правая часть определяет количество теплоты, которое накапливается в единичном объеме в единицу времени в томместе теплового поля, где расположен рассматриваемый элемент.

Теперь можно учесть тотепло G, которое выделяется в радиоэлементах за счет превращения электрической энергии втепловую. Поскольку удельное тепловыделение определяется через количество теплоты, выделяемой в единичном объеме за единицу времени, то можно прибавить соответствующийчлен к левой части уравнения (15). Приходим к выражению:(JX+ – JX–)(JY+ – JY–)(JZ+ – JZ–)ijkt+1–ijkt+++G = CУД(16)hXhYhZУдельное тепловыделение G стоит в левой части уравнения потому, что оно вносит теплоту в рассматриваемый объем.3.2. Уравнение теплопроводности для дискретной модели блокаВыразим плотности потоков J [Дж/м2с] через температуру в узлах сетки.