Голдаев С.В., Ляликов Б.А. Основы математического моделирования в теплотехнике

Подождите немного. Документ загружается.

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

боковых и нижняя стороны ее теплоизолированы, а на верхней поддержива-

ется постоянная температура, равная 10 безразмерным единицам.

Математически рассматриваемый процесс описывается задачей Дирихле

для уравнения Лапласа:

∂

.:,0:0;0:,0:0

ttbyyytaxtx

(6.50)

Аналитическое решение этой задачи имеет вид [43]:

()

[]

(

)

[

]

()( )

[]

∑

∞

+⋅+

/1212

/12/12sin

abkshk

aykshaxk

yxt

ππ

(6.51)

Для повышения точности вычисления отношения гиперболических сину-

сов производилось по формуле

()

[]

(

)

[

]

()

(

)

2exp1

expexp

byrbyr

rbsh

rysh

−−

−

Используя аппроксимационные формулы для производных второго по-

рядка, получаем разностный аналог исходной задачи:

,10,0,0,0

,0,,1,0

1,,1,

,1,,1

====

−

+−

−+−+

niijmj

jijijijijiji

tttt

ttt

(6.52)

где h = a/m, l = b/n - шаги по координатам.

Замечание. При записи дискретных значений температур использованы

нижние двойные индексы, чтобы легче было представлять их в программе в

виде двухмерных массивов.

Для решения системы линейных уравнений воспользуемся методом Гаус-

са – Зейделя, который является одним из самых распространенных итераци-

онных алгоритмов, отличающийся простотой и легкостью программирования

[ 18], [20], [22], [23], [25].

Проиллюстрируем его применение на решении сист

емы трех уравнений.

Выражаем неизвестные x

, x

и x

соответственно из первого, второго и

третьего уравнений системы:

()

⎪

⎩

⎪

⎨

⎧

−−=

33232131313

22321312122

1131321211

axaxabx

Задавая начальные приближения неизвестных, находим их уточненные

значения следующим образом:

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

()

(

)

(

)

()

( ) () ()

()

() () ()

()

⎪

⎩

⎪

⎨

⎧

−−=

3323213131

2232131212

113132121

axaxabx

kkk

Для рассматриваемой задачи, в соответствии с данным алгоритмом, ре-

шение системы линейных уравнений представляется в виде

(

)

()

1,1,,1,1

⎟

⎠

⎞

⎜

⎝

⎛

+++

−+−+

ttktt

jijihljiji

(6.53)

Итерационный процесс по этой формуле прекращается при выполнении

в каждом узле сетки условия

(

)

(

)

<−

+ k

Ниже представлен листинг программы, реализующий описанный алго-

ритм. В теле программы имеются две процедуры, первая из них предназна-

чена для вычисления температур по формуле (6.51), а вторая – для численно-

го решения системы линейных уравнений.

program dirixle;

uses crt;

type mat = array[1..33,1..33] of real;

var a,b,e,u,h,l,t:real;l1,m,n,i,j:integer; f1,f2:mat;

procedure an(a,b,e,u,h,l:real;m,n:integer; var f:mat);

const pi=3.14159265; var p,q,x,y,s,r,z:real; i,j,k:integer;

begin

p:=pi/a; q:=-2*b; z:=4*u/pi;

for i:=2 to m div 2+1 do begin

x:=(i-1)*h; for j:=2 to n do begin

y:=(j-1)*l; s:=0; k:=1; repeat r:=k*p;

r:=sin(r*x)*(exp(r*(y-b))-exp(-r*(y+b)))/(k*(1-exp(r*q)));

s:=s+r; k:=k+2; until abs(r/s)<=e; f[i,j]:=z*s end end end;

procedure fd(e,u,t:real;m,n:integer; var f:mat);

var i,j,k:integer; r:real;

begin

for i:=1 to m+1 do

for j:=1 to n+1 do

if j=n+1 then f[i,j]:=u else f[i,j]:=0;

repeat k:=0;

for i:=2 to m do

for j:=2 to n do

begin

r:=(f[i-1,j]+f[i+1,j]+(f[i,j-1]+f[i,j+1])*t)/(2+2*t);

if abs(r-f[i,j])>e*abs(r) then k:=1; f[i,j]:=r; end;

until k=0; end;

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

begin clrscr;

repeat write('a,b,e,u,m,n '); readln(a,b,e,u,m,n);

h:=a/m; l:=b/n; t:=sqr(h/l);

writeln(' a=', a:5,' b=', b:5, ' e=',e:5,' u=',u:5,' h=',h:5,' l=', l:5);

an(a,b,e,u,h,l,m,n,f1); fd(e,u,t,m,n,f2);

writeln(' - «analit. chislen.');

for i:=2 to m div 2+1 do

for j:=2 to n {div 2+1} do

if (i mod 2 = 1) and (j mod 2 =1) then

writeln('f[',i:2,',',j:2,']=',f1[i,j]:8:5,' f[',i:2,',',j:2,']=',f2[i,j]:8:5);

until readkey =#27 end.

В диалоговом режиме осуществляется ввод:

l1 – условное число, при нулевом значении которого срабатывает первая

процедура;

a =4,b = 4 – размеры системы;

ε = 0,001- погрешность решения системы разностных уравнений;

tb = 10 – безразмерная температура верхней стенки;

m = 8, n = 8 – число разбиений области интегрирования вдоль координат

x и y.

Затем вычисляются шаги сетки h и l по

коорди

натам x и y, а также значе-

ние переменной k

, необходимой в подпрограмме.

После обращения к подпрограммам аналитического или разностного ме-

тодов осуществляется вывод результатов на дисплей.

Ниже приведена выборка результатов расчетов температуры по данной

программе. Вследствие симметрии задачи на печать выведены значения

температуры в узлах левой половины прямоугольной области.

a,b,e,u,m,n 4 4 0.005 10 8 8

a= 4.0E+00 b= 4.0E+00 e= 5.0E-03 u= 1.0E+01

t= 1.0E+00 h= 5.0E-01 l= 5.0E-01

analit chiclen

x=1.500 y=1.500 ua= 0.67973 u= 0.67535

x=1.500 y=2.500 ua= 1.82030 u= 1.82119

x=1.500 y=3.500 ua= 4.31932 u= 4.30024

x=2.500 y=1.500 ua= 0.95412 u= 0.94569

x=2.500 y=2.500 ua= 2.50003 u= 2.47942

x=2.500 y=3.500 ua= 5.40426 u= 5.34837

Как видно из этих данных, результаты численного решения хорошо со-

ответствуют аналитическим, что свидетельствует о достоверности вычисли-

тельного эксперимента.

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

7. МЕТОДЫ РЕШЕНИЯ НЕСТАЦИОНАРНЫХ ЗАДАЧ

ТЕПЛОПРОВОДНОСТИ

7.1. Аналитические методы решения обыкновенных

дифференциальных уравнений

Дифференциальным называется уравнение, содержащее неизвестные

функции у, независимую переменную х и производные от неизвестных

функций

)('''

,...,,

yyy :

(

)

(

)

0,...,,,, =

′′′

yyyyxF

(7.1)

Если неизвестные функции зависят от одной независимой переменной х,

то дифференциальное уравнение является обыкновенным. Если они зависят

от нескольких независимых переменных, дифференциальное уравнение на-

зывается уравнением с частными производными. Порядком дифференциаль-

ного уравнения называется наивысший из порядков производных, входящих

в уравнение. Дифференциальное уравнение считается линейным, если оно

линейно относительно искомой функции и ее производных. Совокупность

уравнений, содерж

ащих несколько неизвестных функций и их производных,

образует систему дифференциальных уравнений [37].

Решить дифференциальное уравнение, значит, найти его общий интеграл.

Под общим интегралом понимается соотношение между независимой пере-

менной, зависимой переменной и произвольными постоянными, число кото-

рых равно порядку дифференциального уравнения. Это соотношение при

всех допустимых з

начениях независимой переменной должно удовлетворять

данному дифференциальному уравнению. Общий интеграл геометрически

выражается семейством кривых. Для выделения одной кривой из этого се-

мейства, т. е. для получения частного решения, необходимо определить про-

извольную постоянную, задав начальные условия. В уравнениях n-го поряд-

ка эти условия налагаются на переменную х и ее производные, а именно - за-

даются зн

ачения:

(

)

....,,,:

0000

yyyxx =

(7.2)

Задача определения решения этого уравнения при заданных начальных

условиях называется задачей Коши для обыкновенного дифференциального

уравнения.

Так, дифференциальное уравнение первого порядка dy/dx = y имеет ре-

шение у = Сe

, где С - произвольная постоянная. При различных значениях

постоянной С получается семейство кривых. Выбор начального значения у =

при х = х

позволяет выделить из этого семейства одну определенную кри-

вую, т. е. найти частное решение.

Точные методы решения дифференциальных уравнений (разделение пе-

ременных, использование интегрирующего множителя, вариации постоян-

ной, операционный и др.) позволяют выразить его через элементарные или

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

специальные функции [36], [46]. Однако классы уравнений, для которых су-

ществуют точные методы решения, довольно узки и охватывают только ма-

лую часть возникающих на практике задач.

Пример 7.1. Плоская стальная стенка с одной стороны омывается газом с

температурой T

= 1073 K, коэффициент теплоотдачи от газов к стенке равен

= 150 Вт/(м

.К). Толщина стенки 2 мм. Другая сторона стенки, имеющая в

начальный момент температуру T

= 273 К, теплоизолирована. Теплофизиче-

ские характеристики материала стенки: λ = 17 Вт/м.К, с = 500 Дж/кг К, ρ

=7900 кг/м

. Влияние лучистого теплообмена не учитывать.

Определить время, в течение которого плоская стальная стенка нагреется

до температуры 990 К.

Решение. Как установлено в теории теплопроводности [19], [27], [30],

[43], при малых числах Био температуры на поверхности пластины или ци-

линдра мало отличаются от температур на оси, поскольку процесс выравни-

вания температуры в теле происходит существенно интенсивнее, чем отвод

теплоты с поверхности. В этом случая распределением температур

ы по тол-

щине пренебрегают и вместо уравнения с частными производными исполь-

зуют обыкновенное дифференциальное уравнение для осредненной темпера-

туры.

Вычисляя критерий Био, имеем: Bi = α

δ/λ = 0,0176. Таким образом,

можно использовать данное приближение. Количественный анализ, как от-

мечалось в гл.1, удобнее осуществлять в безразмерных переменных, которые

определяются следующим образом:

−

====

−

= ,,,,

(7.3)

Математически задача Коши в рассматриваемом примере записывается

так:

()

uuBi

dFo

−=

.0:0

uFo

(7.4)

Переменные в данном уравнении разделяются, и решение, удовлетво-

ряющее нулевому начальному условию, имеет вид:

(

)

[

]

BiFouu

−

= exp1

(7.5)

Безразмерные температуры среды и пластинки равны: u

= 2,93, u =2,63.

Подставляя эти значения в (6.5), находим критерий Фурье:

(

)

5,128//1ln

−

= BiuuFo

Отсюда

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

cFoc 5,119/

λρδτ

Известно, что любое уравнение n-го порядка можно заменить системой из

n уравнений первого порядка [37]. Для этого вводят новые переменные:

(

)

,...,,,

−

′

yyyyyy

в результате чего получают эквивалентную систему:

()

⎪

⎩

⎪

⎨

⎧

.,...,,,

,.....

yyyxf

(7.6)

Указанный прием позволяет свести решение дифференциального уравне-

ния n-го порядка к решению системы n - уравнений первого порядка.

Порядок системы зависит от сложности физико-математической модели

исследуемого процесса. В промышленной теплотехнике, как правило, одно-

временно протекает несколько процессов [1-6], [13], [14], [17], [47]. Напри-

мер, при мокрой очистке дымовых газов, наряду с процессом абсорбции,

осаждением пыли на каплях, т.е. массоперед

ачей, происходит изменение

температур обеих фаз, т.е. теплообмен. Так, предложенная в [47] модель теп-

ломассообмена и пылеулавливания в прямоточных распылительных аппара-

тах содержала двадцать четыре нелинейных обыкновенных дифференциаль-

ных уравнения, которые можно решить только численным методом.

В свою очередь алгоритмы решения одного уравнения первого порядка

распространяются на си

стемы таких уравнений. На основании этого в даль-

нейшем рассмотрение численных методов будем проводить на примерах ре-

шения задачи Коши:

у’ = f (x, y), при

yyxx

(7.7)

Разработаны и аналитические методы (последовательного дифференци-

рования, неопределенных коэффициентов, Пикара), дающие приближенное

решение дифференциального уравнения в виде аналитического выражения

[32], [36], [46].

Опишем алгоритм применения метода последовательных приближений

(или метод Пикара) для решения задачи Коши (6.7). Искомая зависимость

y(x) является пределом последовательности функций y

(x), которые находят-

ся по рекуррентной формуле

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

() ( )

dxyxfyxy

∫

−

(7.8)

Получим приближенное решение задачи Коши из рассмотренного выше

примера 7.1. Соотношение (7.8) для (7.4) запишется следующим образом:

()

∫

−

−=

nln

dFouuBiu

Первое приближение равно

∫

≅≅

FoBiudFouBiu

Подставляя его в предыдущее и выполняя интегрирование, получаем:

() ()

.2/1

BiFoFoBiudFoFoBiuuBiu

−≅−≅

∫

Сравним структуру второго приближения с аналитическим решением

(7.5). Разлагая в нем экспоненту в ряд и оставляя три первых члена, получаем

такую же зависимость, как и u

Поскольку при реализации метода Пикара возникает необходимость на-

хождения интегралов, то его эффективность снижается в случае сложных

выражений функции f(x,y).

Этот метод применим и для решения систем дифференциальных уравне-

ний невысокого порядка.

7.2. Численные методы решения обыкновенных

дифференциальных уравнений

Среди численных методов широко распространенными являются разно-

стные методы. Они основаны на введении некоторой разностной сетки в рас-

сматриваемой области. Значения производных, начальные и граничные усло-

вия выражаются через значения функций в узлах сетки, в результате чего по-

лучается система алгебраических уравнений, называемая разностной схемой.

Решая эту систему уравнений, можно найти в узлах сетк

и значения сеточных

функций, которые приближенно считаются равными значениям искомых

функций.

Сущность данных методов состоит в следующем [19], [21-24], [33], [37].

На отрезке решения [x

, x

] выбирается некоторое множество точек, назы-

ваемое сеткой:

< х

< ... < х

В полученных точках вычисляются приближенные значения у

, у

, ... ,

уn решения задачи Коши. Разность h = х

к+1

- x

называют шагом сетки. В

большинстве случаев величину h полагают постоянной. Тогда

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

= х

+ k h , k = 0,1,2, ... .

Метод, дающий формулу для вычисления значения у

l+1

по k предыдущим

значениям у

, у

l-1

, ... , у

l-k+1

называется k - шаговым методом. Если k =1, ме-

тод является одношаговым. Среди методов этого типа широко применяются

методы Рунге - Кутта. При k > 1 метод называется многошаговым. Примера-

ми таких методов могут служить методы прогноза и коррекции.

Рассмотрим наиболее употребительные методы решения задачи Коши

(7.7). При явных (одношаговых) методах решения каждое новое значение y

k+1

находится по известному предыдущему значению y

. Этим методам свойст-

венно «самостартование», а также возможность изменения шага h = x

k+1

– x

в процессе вычислений.

Простой метод Эйлера имеет первый порядок точности и записывается

так:

()

1 kkkk

xyhfyy

(7.9)

Эту схему интегрирования используют, когда функция f(x, y) имеет дос-

таточно сложный вид и затрачивается много времени на ее вычисление, а ко-

личество точек x

большое.

В соответствии с усовершенствованным методом Эйлера – Коши вначале

вычисляется приближенное значение функции простым методом Эйлера:

(

)

kkk

xyhfyy +=

(7.10)

а затем оно уточняется с помощью формулы

()

(

)

[]

++=

kkkk

yxfyxf

(7.11)

Для повышения точности интегрирования шаг h целесообразно делать

малым, что приводит к увеличению объема вычислений.

На рис. 7.1 представлено сравнение результатов решения примера 7.1 ме-

тодом Эйлера с аналитическим (сплошная кривая), которое иллюстрирует

отмеченную выше особенность описанных алгоритмов.

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

Рис. 7.1. Влияние шага интегрирования на точность решения Коши

методом Эйлера

Недостатком описанных явных методов является возникновение чи-

словой неустойчивости решения при большом шаге h. Установлено, напри-

мер, что в методе Эйлера накладывается следующее ограничение на шаг:

(

)

.//2 yfh

∂

−

≤

(7.12)

Точное аналитическое выражение устойчивости в практических зада-

чах из-за их сложности получить не удается, поэтому требуемую погреш-

ность достигают путем постепенного уменьшения шага и сравнения полу-

ченных результатов. Когда отличие их не будет превышать, скажем, 0,5% и

менее, то дальнейшее дробление шага прекращают.

Метод Рунге-Кутта имеет четвертый порядок точности и заключается в

предварительно

м вычислении четырех значений функции. Геометрическая

интерпретация данной схемы интегрирования состоит в том, что смещение

из точки (х

, у

,) в точку (х

k+1

, у

k+1

) происходит через промежуточные точки,

в каждой из которых определяется направление касательной. В итоге смеще-

ние из точки (х

,у

.) в (х

k+1

, у

k+1

.) выполняется вдоль некоторого усреднен-

ного направления. Расчет выполняется по следующим формулам:

2,4

1,6

0,8

0,0

0 40 80 120 Fo

h=0.2

h=0.1

h=0.05

Голдаев С.В., Ляликов Б.А. Основы математического моделирования

в теплотехнике: Учебное пособие. - Томск: Изд. ТПУ, 1999. – 106 с.

()

(

)

()()

[]

.6/)(2

,,,2/,2/

,2/,2/,,

43211

3423

121

hmmmmyy

myhxfmmyhxfm

myhxfmyxfm

kkkk

++++=

++=++=

(7.13)

Метод Рунге – Кутта устойчив по отношению к ошибкам в начальных

условиях и к искажениям в правой части уравнения, обладает достаточно вы-

сокой точностью, удовлетворяющей многим практическим приложениям,

легко программируется. Величину шага можно изменить на любом этапе вы-

числений. Среди немногочисленных недостатков этого метода выделяют от-

сутствие легко определяемой оценки ошибки алгоритма. Выбор шага в этих

условиях осуществляют путем многократного повторен

ия расчетов с умень-

шающимися шагами. Добавление пятого коэффициента m

в модификации

Мерсона не повышает точности метода, но позволяет выразить погрешность

ε через m

, m

, …, m

. Расчетные формулы имеют вид:

()( ) ()

(

)

() ( )()

()( )

.5/2/42/9,2/4

,2/44/33,3/

,8/33,2/3/

,2/,3/3/

,,3/3/,,3/

54315411

4315

314

213

121

mmmmmmmyy

mmmyhxfhm

mmyhxfhm

myhxfhmyxfhm

kkkk

−+−=+++=

−−++=

+++=

(7.14)

Неявные методы, у которых y

k+1

находится по значению f(y

k+1

), ли-

шены этого недостатка. Так, наиболее простая абсолютно устойчивая схема

имеет вид:

(

)

(

)

(

)

[

]

kkkk

yfyfhyy

(7.15)

Однако она не всегда реализуема, например, если f(y) = - sin(y). Кроме

того, возникает необходимость на каждом шаге интегрирования решать

трансцендентное уравнение, что, как отмечалось в гл.3, требует выбора под-

ходящего метода.

Встречаются случаи, когда разностные методы не дают приемлемых ре-

зультатов. Это наблюдается в задачах с быстрорастущими решениями. Тогда

делают специальные математические преобразования, устраняющие или

сглажив

ающие эти особенности, после чего допускается использование рас-

сматриваемых численных методов [19], [23].

Описанные алгоритмы применимы и для систем обыкновенных диффе-

ренциальных уравнений. Например, расчетные формулы метода Рунге-Кутта

в случае системы двух уравнений