Дворецкий С.И. и др. - Комп. моделир. и оптимизация технологич. процессов и оборуд. [2003]

процесса теплообмена

cpρG(T(1)(τ)−T(1вх) )+ kтS(1)(T(1)(τ)− tk )= cpρV(1) dTτ(1) ; d

cpρG(T(2)(τ)−T(1))+ kтS2 (T(2)(τ)−tk )= cpρV(2) dTτ(2) ; d

…………………………………………

cpρG(T(m)(τ)−T(m−1))+ kтS(m)(T(m)(τ)−tk )= cpρV(m) dT(τm) . d

Если режим движения потока жидкости в аппарате описывается ячеечной моделью с обратными потоками (рис. 3.11), то математическая модель процесса теплообмена принимает вид:

cpρ[GT(1вх) (τ)+ gT(2)(τ)− (G − g)T(1)(τ)]± kтS(1)(T(1)(τ)−tk )= cpρV(1) dTdτ(1) ;

cpρ[(G + g)(T(1) −T(2))+ g(T(3) −T(2))]+ kтS(2)(T(2)(τ)−tk )= cpρV(2) dTτ(2) ; d

………………………………………………………………..

cpρ[(G + g)(T(i −1) −T(i))+ g(T(i +1) −T(i))]+ kтS(i)(T(i)(τ)−tk )= cpρV(i) dTdτ(i) ;

cpρ[((G + g)T(N −1) − (G + g)T(N ))]+ kтS(N )(T(N )(τ)−tk )= cpρV(N ) dT(τN ) . d

Зададим начальные условия для записанных выше уравнений динамики процесса теплообмена

T(1)(0)= T10 , T(i)(0)= Ti0 , T(N )(0)= TN 0 .

Для получения уравнений модели статики процесса теплообмена необходимо

dTdτ(i) = 0, i =1, 2, ..., N :

cpρ[GT(1вх) + gT(2) − (G + g)T(1)] + kтS(1)(T(1) −tk )= 0 , cpρ[(G + g)(T(1) −T(2))+ g(T(3) −T(2))]+ kтS(2)(T(2) −tk )= 0 ,

………………………………………………………………..

cpρ[(G + g)(T(i −1) −T(i))+ g(T(i +1) −T(i))]+ kтS(i)(T(i) −tk )= 0 ,

………………………………………………………………..

cpρ[((G + g)T(N −1) − (G + g)T(N ))]+ kтS(N )(T(N ) −tk )= 0 .

Пример. Оценим профиль температуры нагреваемого потока жидкости, исходя из различных гидродинамических моделей движения этого потока. Зададим условия осуществления теплообмена:

		G =1000 кг/ч ;	сp = 2520 Дж/(кгК) ;		ρ =1200 кг/м3 .
Обогрев осуществляется насыщенным водяным паром,					имеющим температуру t =120 oC . Диаметр
цилиндрической	поверхности теплообмена			равен d = 0,5 м .	Коэффициент	теплопередачи составляет
kт = 600 Вт/м2K , длина		теплообменника	–	1,5 м, параметры ячеечной		и диффузионной	модели:
n = 3, Dт = 3,54 10−4	м2/с,	на рис. Приведены		результаты расчета температурного профиля			по длине
теплообменника.
Они					свидетельствуют о значительном разбросе

температур

Более

температуры

ячеечная и

конечные

практически

температур

Т,0°CС
110	1			2
100
90	3			4
80
70
60	0,2	0,4	0,6	0,8	1,0	1,2	1,4 l,мм
0

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

1 – идеальноесмешение; 2 – идеальноевытеснение; 3 – ячеечная модель; 4 – диффузионная модель

для различных моделей гидродинамики.

реальный характер изменения

по длине теплообменника отражают

диффузионная модели. При этом

температуры для данных моделей

совпадают, но, тем не менее, профили

различаются существенно.

Приведенный пример подчеркивает важность учета реальной структуры потоков в аппарате и его адекватного описания гидродинамическими моделями.

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

Для простоты будем рассматривать одномерные процессы теплопроводности. Они имеют место, например, в длинном тонком металлическом стержне, нагреваемом с одного из торцов при условии, что стержень изотропен. Его начальная температура в любом поперечном сечении не зависит от y, z (это условие должно выполняться и на торцах стержня), а потерями тепла с боковой поверхности можно пренебречь.

Рассмотрим произвольное сечение стержня с координатой x. Пусть ρ(x), cp(x), k(x) – соответст-

венно плотность, удельная теплоемкость и коэффициент теплопроводности в точках этого сечения. Запишем уравнение распространения этого тепла в стержне (уравнение теплопроводности) на некотором отрезке (x1, x2) за некоторый промежуток времени (t1, t2), применяя закон сохранения энергии (в интегральной форме)

t2		∂T				∂T				x2t2
∫	k		(x,τ)		−k		(x,τ)		dτ+ ∫∫F(ξ,τ)dξdτ=
		∂x		x=x		∂x		x=x
t1				2				1		x1 t1	.
	x
= ∫2cpρ[T(ξ,t2 ) −T(ξ,t1)]dξ
	x1

Предположим, что функция T (x,t) имеет непрерывные производные Txx и Tt .

Пользуясь теоремой о среднем, получаем равенство

ла, м3.

4 Внутри потока может возникать или поглощаться теплота. Выделение быть описано плотностью тепловых потоков F(x, y, z, τ) в точке (x, y, z) в момент зультате действия этих источников за промежуток времени (τ, τ + ∆τ) выделится

лоты

где cp – удельная теплоемкость,

Коэффициент теплоотдачи α выражает количество тепла, отданного единицей поверхности (S = 1 м2) в единицу времени (τ = 1 с) при разности температур (θ−T ). Заметим, что α не является

постоянной величиной, а зависит от многих параметров.

3 Количество теплоты, которое необходимо сообщить однородному телу, чтобы повысить его температуру на величину ∆T равно

Q = cpm∆T = cpρV∆T ,

Дж ; m – масса тела, кг; ρ – плотность тела, кг/м3; V – объем те-

кг К

	∂T					∂T
k		(x, τ)			− k		(x, τ)			∆t + F(x4 ,t4 )∆x∆t =
	∂x					∂x
			x=x	2				x=x
							1		τ=t3
					= {cpρ[T (ξ,t2 ) −T (ξ,t1)]}ξ=ξ						∆x,
											3

которое при помощи теоремы о конечных приращениях можно преобразовать к виду:

∂		∂T				∂T
	k		(x, τ)	∆x∆t + F (x4	, t4 )∆x∆t = cpρ			∆x∆t,
		∂x
∂x			x=x5			∂t x=x3
			t =t3				t =t5

где t3, t4, t5 и x3, x4, x5 – промежуточные точки интервалов (x1, x2) и (t1, t2).

dQdτ = α(θ−T )S ,

где θ, T – температура поверхности твердого тела и потока, соответственно; α – коэффициент теп-

лоотдачи,		Вт
			.
		м2град

теплоты может времени τ. В ре-

количество теп-

dQ = F(x, y, z, τ)dxdydzdτ ,

или в интегральной форме:

τ2 x2 y2 z2

Q = ∫ ∫ ∫ ∫F(x, y, z, τ)dxdydzdτ .

τ1 x1 y1 z1

Отсюда после сокращения на произведение ∆x ∆t

получим:

∂T

∂

(k

∂T )

+ F(x,t)

x=x4

= c

ρ

.

p

∂t

∂x

∂x

x=x5

t =t4

t =t5

t =t3

x=x3

Все эти рассуждения относятся к произвольным промежуткам

(x1, x2) и (t1, t2). Переходя к

пределу при x1, x2 → x и t1, t2 → t , получим уравнение

∂

∂T

∂T

,

k

+ F (x,t) = cpρ

∂t

∂x

∂x

называемое уравнением теплопроводности.

Частные случаи.

1. Если стержень однороден, то ρ, cp, k = const, и мы получаем линейное уравнение теплопро-

водности

	k	Tt = a2Txx + f (x,t) ,
где a2 =		– коэффициент температуропроводности; f (x,t) =	F(x,t)	.
	cpρ		cpρ

Если источники отсутствуют, т.е. F(x,t) = 0 , то уравнение теплопроводности примет вид

Tt = a2Txx .

2. В случае теплообмена с окружающей средой, подчиняющегося закону Ньютона, количество тепла, теряемого стержнем, рассчитываемого на единицу длины и времени, равноF0 = α(T −θ) ,где

θ(x, t) – температура окружающей среды, α – коэффициент теплообмена.

Поскольку в нашем приближении не учитывается распределение температуры по сечению, то действие поверхностных источников эквивалентно действию объемных источников тепла. Таким образом, плотность тепловых источников в точке x в момент времени t равна F = F1(x,t) − α(T − θ) ,

где F1(x, t) – плотность других источников тепла.

Если стержень однороден, то уравнение теплопроводности с боковым теплообменом имеет

2		α		F1 (x,t)
следующий вид: Tt = a Txx − h1T + f (x, t) , где h1	=		; f (x, t) = h1θ(x,t) +		– известная функция.
		cpρ		cpρ

3. Коэффициенты cp и k, как правило, являются медленно меняющимися функциями температуры. Поэтому сделанное выше предположение о постоянстве этих коэффициентов возможно

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

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

∂		∂T			∂T
	k (T, x)			+ F(x, t) = cp (T, x) ρ(T, x)		.
					∂t
∂x		∂x

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

Начальное условие состоит в задании значений функции T (x, t) в начальный момент t0 , т.е.

T (x, t0 ) = T0 (x), 0 ≤ x ≤ l .

Граничные условия могут быть различными в зависимости от температурного режима на торцах стержня. Рассматривают три основных типа граничных условий.

1 Наторцахстержнявлюбоймоментвременизадаетсятемпература:

T (0, t) = T1(t), T (l, t) = T2 (t), t > 0 .

	2 На торцах стержня задаются потоки теплоты как функции			времени:
∂T	(0, t)= ν1(t),	∂T	(l, t)= ν2 (t) .
∂x		∂x

К этим условиям мы приходим, если задана величина теплового потока Q(l,t) , протекающего

через торцевое сечение стержня,

Q(l, t) = −k ∂∂Tx (l, t),

откуда ∂∂Tx (l,t)= ν(t) , где ν(t) – известная функция, выражающаяся через заданный поток Q(l,t) по

формуле

ν(t) = − Q(l, t) . k

3 На торцах стержня задаются линейные соотношения между производной и функцией, например, для x = l:

∂∂Tx (l, t) = −h2 [T (l, t) − θ(t) ].

Это граничное условие соответствует теплообмену по закону Ньютона на поверхности тела с окружающей средой, температура которой θ известна.

Пользуясь двумя выражениями для теплового потока, вытекающего через сечение x = l:

Q = α(T − θ) и Q = −k			∂T		получаем математическую формулировку третьего граничного условия в
			∂x
виде
	∂T	(l, t) = −h2 [T (l, t) − θ(t) ]		,
	∂x

где h2 = α/ k – коэффициент теплообмена; θ(t) – некоторая заданная функция.

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

Более сложный (нелинейный) вариант условий на торцах отвечает сильно нагретому и поэтому излучающему энергию стержню, не контактирующему с какими-либо телами. Тогда в еди-

ницу времени стержень теряет на своих границах (торцах) энергию, равную σT 4 (0, t) и σT 4 (l, t) со-

ответственно. В результате получаются условия:

σT 4 (0, t) = − k(T)	∂T			, σT 4 (l, t) = −k(T)	∂T		, t > 0 .
	∂x		x=0		∂x		x=l

Моделирование процесса диффузии газа в полой трубке.

Если среда неравномерно заполнена газом, то имеет место диффузия его из мест с более высокой концентрацией в места с меньшей концентрацией. Это же явление имеет место и в растворах, если концентрация растворенного вещества в объеме не постоянна.

Рассмотрим процесс диффузии в полой трубке или в трубке, заполненной пористой средой, предполагая, что во всякий момент времени концентрация газа (раствора) по сечению трубки одинакова. Тогда процесс диффузии может быть описан функцией c(x, t) , представляющей

концентрацию в сечении x в момент времени t .

Согласно закону Нернста, масса газа, протекающая через сечение x за промежуток времени (t, t + ∆t) , равна

dG = −D ∂∂cx (x,t)Sdt =WSdt ;

W = −D ∂∂cx ,

где D – коэффициент диффузии; S – площадь сечения трубки; W (x,t) – плотность диффузионного

потока, равная массе газа, протекающего в единицу времени через единицу площадки. По определению концентрации, количество газа в объеме V равно,

G = cV .

Отсюда получаем, что изменение массы газа на участке трубки (x1, x2 ) при изменении концентрации на ∆C равно

x2

∆G = ∫ε(x)∆cSdx ,

x1

где ε(x) – коэффициент пористости.

При выводе уравнения диффузии будем считать, что в трубке нет источников вещества, и диффузия через стенки трубки отсутствует.

Составим уравнение баланса массы газа на участке (x1, x2 ) за промежуток времени (t1,t2 ) :

t1				∂u	(x2		, τ) − D(x1)	∂u	(x2		x2
S ∫		D(x2 )		∂x				∂x		, τ)	dτ = S ∫ε(ξ)[u (ξ,t2 )−u (ξ,t1 )]dξ .
t											x
1											1
		Отсюда, подобно выводу уравнения теплопроводности, получим уравнение
∂			∂c			∂c	,
		D		= ε
∂x						∂t
			∂x

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

Если коэффициент диффузии постоянен D = const , то уравнение диффузии принимает вид

ct = a2cxx , где a2 = Dε .

Если коэффициент пористости ε =1 , а коэффициент диффузии постоянен, то уравнение диффузии имеет вид ct = Dcxx .

Для получения единственного решения уравнения диффузии необходимо к уравнению присоединить начальные и граничные условия.

Как бы глубоки и разнообразны ни были методы качественного анализа математических моделей, область их применимости всегда ограничена. Это – либо простые, главным образом, линейные, либо отдельные фрагменты сложных, в том числе нелинейных моделей. Единственным универсальным способом исследования моделей является применение численных методов для нахождения приближенного решения поставленной задачи. Для решения нелинейных задач теплопроводности и диффузии широко применяется метод конечных разностей [12], который состоит из двух этапов: на первом строятся дискретные аналоги исходных моделей и изучаются их свойства, на втором дискретные уравнения (как правило, нелинейные алгебраические уравнения высокой размерности) решаются численно.

О переходек дискретныммоделям теплопроводности идиффузии.

Метод конечных разностей состоит в следующем. Область непрерывного изменения аргументов ( x и t ) заменяется конечным (дискретным) множеством точек (узлов), называемым сеткой.

Вместо функции непрерывного аргумента рассматриваются функции дискретного аргумента, определенные в узлах сетки и называемые сеточными функциями. Производные, входящие в дифференциальные уравнения, заменяются (аппроксимируются) при помощи соответствующих разностных соотношений. В результате такой замены дифференциальное уравнение заменяется системой алгебраических (разностных) уравнений. Начальные и краевые условия то же заменяются разностными начальными и краевыми условиями. Естественно требовать, чтобы полученная таким образом разностная краевая задача была разрешима, и ее решение при увеличении числа N узловсетки приближалось (сходилось) крешению исходной задачи.

Пусть область изменения аргументов (x, t) есть прямоугольник Π = (0 ≤ x ≤1, 0 ≤ t ≤ T ). Построим на отрезке 0 ≤ x ≤1 сетку {xi = ih, i = 0,1,..N1} с шагом h =1N1 и сетку {t j = jτ, j = 0,1,..., N2} с шагом τ = TN2 на отрезке 0 ≤ t ≤T .

Множество узлов (xi , t j ) с координатами xi = ih и t j = jτ назовем сеткой в прямоугольнике Π и обозначим через ωhτ сетку{(xi = ih, t j = jτ), i = 0, N1; j = 0, N2 }. Эта сетка равномерна по каждой из переменных x и t .

Пусть у – сеточная функция, заданная на ωhτ . Будем обозначать yij = y(xi ,t j ) значение сеточной функции у в узле (xi , t j ) сетки ωhτ . Непрерывной функции T (x, t) или c(x, t), где (x, t) Π , будем ставить в соответствие сеточную функцию yii = T (xi , t j ) c (xi , t j ).

Рассмотрим теперь производную v′x функции v(x) . Заменить ее разностным выражением можно

бесчисленным множеством способов. Простейшими являются замены вида v

′

vi − vi −1

−

– левая

x

~

= L v

i

h

h

разностная производная (левое разностное отношение); v

′

vi+1 − vi

+

– правая разностная про-

x

~

=

L v

i

h

h

vi +1 − vi −1

изводная; v

′

0

– центральная разностная производная, где знак ~ означает соответ-

x

~

= L v

i

2h

h

ствие или аппроксимацию.

Обращаясь к формулам для L±h , видим, что L−hvi и L+hvi

аппроксимируют Lv = v′

с первым по-

рядком. Выражения для L−v

i

содержат значения v в двух узлах x = x

i

и x = x

i−1

сетки. Говорят, что

h

оператор L−h является двухточечным или оператором первого порядка.

Множество узлов, значения сеточной функции в которых входят в выражение Lh v , называют

шаблоном оператора Lh в точке xi . Очевидно, что шаблон оператора L−h состоит из узлов xi и xi−1 ,

а шаблон L+h						– из узлов xi				и xi+1 .
	Возьмем теперь трехточечный оператор, определенный на шаблоне xi−1 , xi , xi+1 :
(σ)	v		+			−		=	σvi +1 + (1− 2σ)vi − (1− σ)vi −1		,
L		i	= σL v	i	+ (1− σ)L v		i
h			h			h				h

где σ – произвольное число. В частности при σ = 12 получаем центральную разностную производ-

ную

0

, которая аппроксимирует v (x) со вторым порядком.

Lhvi

′

Рассмотрим теперь вторую производную Lϑ = ϑ . Выберем трехточечный шаблон, состоящий

′′

из узлов xi−1 ,

xi , xi+1

и рассмотрим разностный оператор

vi +1 − 2vi + vi −1

vi +1 − vi

− vi − vi −1

L

v

i

=

=

h

h

,

h

h2

h

1

(L+v

)=

vi +1 − vi − vi − vi −1

=

vi+1 − 2vi + vi −1

.

L

v

i

=

i

− L−v

h

h

h

h

h

h i

h

h2

На практике аппроксимация производных на многоточечных шаблонах используется редко, так как при увеличении шаблона обычно увеличивается объем вычислительной работы, и ухудшаются качества получающихся разностных операторов (в смысле устойчивости).

Рассмотримболее сложный оператор Lu = ∂u

−

∂2u

, где u = u(x,t) – функция двух аргументов x и t , ме-

∂t

∂x2

няющаяся в области Π=(0≤x≤1,0≤t ≤T). Введем сетку ωhτ ={(xi =ih, t j = jτ), i =

; j =

} с шагами h =1 N1 ;

0, N1

0, N2

τ = T N2 . Произведем замену

∂u

u j+1

− u j

∂2u

u j

−2u j +u j

∂t

~

i

i

= utj,+i1 ; с

∂t2

~

i−1

i i+1

= uxxj , i .

τ

h2

В результате получим разностный оператор

u j+ 1

−u j

u j

−2u j

+u j

Lhτ uij+1 =

i

i

−

i−1

i

i

+1

.

τ

h2

Этот оператор определен на шаблоне, состоящем из четырех точек (xi ,t j+1 ),		(xi ,t j ), (xi−1,t j ),
(xi+1,t j ) (рис. 3.13, а).
Оператор Lhτ определен не во всех узлах ωhτ , а только при 0 < i < N и	j > 0 , т.е. во внутренних
узлах. В остальных узлах, называемых граничными, должны быть	заданы	начальные и
краевые

		(i, j+1)				(i-1, j+1)		(i, j+1)		(i+1, j+1)
σ=0
(i-1, j)		(i, j)		(i+1, j)				(i, j)
		а)						б)

Рис. 3.13 Четырехточечный шаблон

условия. Оператор Lhτ имеет первый порядок аппроксимации по τ и второй по h

max Lhτuij − (Lu)ij = 0 (h2 + τ).

ωhτ

Аппроксимируемэтотже оператор Lu нашаблоне вида (рис. 3.13, б)

В результате получим оператор

	Lhτuij +1 =	uij +1 −uij	−	uij−+11 − 2uij +1 + uij++11	,
		τ		h2

аппроксимирующий Lu с тем же порядком точности, что и предыдущий оператор.

Рассмотрим пример постановки разностной задачи для уравнения теплопроводности

∂T

= ∂2T + f (x, t), 0 < x <1, 0 < t < tk ;

∂t

∂x2

T (x, 0)= T0 (x), 0 ≤ x ≤1 ; T (0, t)= µ1(t); T (1, t)= µ2 (t), 0 ≤ t ≤ tk .

Введемравномернуюсетку ωhτ = {(xi = ih, t j

= jτ), i =

; j =

} и запишем соответствующую разно-

0, N1

0, N2

стную краевую задачу

yij +1 − yij

=

yij−1 − 2yij + yij+1

+ ϕij +1,

0 < i < N, j > 0

τ

h2

y0

= T (x );

y j

= µ

(t

j

);

y j

= µ

2

(t

j

),

где ϕj +1

= f (x , t

j

).

i

0 i

0

1

N

i

i

Определим yij +1 : yij +1 = (1− 2γ)yij + γ(yij−1 + yij+1 )+ τϕij +1 , где γ =τ h2 .

Если yij известно,

то по этой формуле можно определить yij +1 во всех узлах i =1, 2,..., N1 −1 (на

слое j +1 ). Так как при

j = 0 задано начальное условие yi0 = T0 (xi ), то последняя формула позволяет

определить от слоя к слою значения yij +1 во всех внутренних узлах сетки ωhτ , используя при этом

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

	yij +1 − yij	=	yij−+11 − 2 yij +1 + yij++11			+ ϕij +1 .
	τ
				h2
	В этом случае для определения								yij +1 на новом слое j +1 получаем систему алгебраических
уравнений
	γy j +1 − (1 + 2γ)y j +1			+ γy j +1	= −y j − τϕj +1			, 0 < i < N .
	i −1		i	i+1	i		i		1

Такая схема называется неявной или схемой с опережением.

Понятие об устойчивости разностных сил.

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

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

Ошибки вычисления можно рассматривать как возмущения начальных данных или правой части уравнения. Отсюда следует, что от схемы надо требовать, чтобы решение разностной задачи мало менялось при малом изменении входных данных задачи (правой части краевых и начальных условий) или, иными словами, чтобы решение непрерывно зависело от входных данных при измельчении сетки. Если это требование выполняется, то такая схема называется устойчивой, в противном случае схема неустойчива.

Разностные схемы для нелинейных уравнений теплопроводности (диффузии).

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

x2	t2		x2 t2
∫cv [T (x, t2 )−T (x, t1 )]dx = ∫		[W (x1, t)−W (x2 , t)]dt + ∫∫ f (x, t)dxdt ,
x1	t1		x1 t1
где T (x, t)	– температура, cv		– объемная теплоемкость, f (x, t) – плотность источников тепла;

W (x, t)= −k (x, t)∂∂Tx (x, t) – тепловой поток; k(x, t, T ) – коэффициент теплопроводности.

Если существуют непрерывные производные

∂T

и

∂

∂T

, то из уравнения баланса следует

∂t

k

∂x

∂x

дифференциальное уравнение теплопроводности

∂T

∂

∂T

c

∂t

=

k(x, t)

+ f (x, t).

∂x

∂x

Естественно при написании разностных уравнений, приближенно описывающих тот или иной процесс, исходить из уравнения баланса. Пусть дана сетка ( xi = ih , t j = jτ ). Для каждой элементар-

ной ячейки (прямоугольника) этой сетки пишется уравнение баланса, которое содержит интегралы от функции и ее производных вдоль границы ячейки. Для их вычисления необходимо предположение о профиле функций. В зависимости от выбора локальной интерполяции как по x , так и по t мы получим различные схемы. Вопрос о выборе интерполяции подчинен требованиям ус-

тойчивости, точности и простоты реализации Для примера рассмотрим стационарное уравнение теплопроводности

	∂		∂T
		k(x, t)		− q(x)T = − f (x, t), 0	< x <1;	k > 0,	q ≥ 0 ,
	∂x		∂x
где q(x)T			– мощность стоков тепла (при q ≤ 0 – источников), пропорциональная температуре T (x).

Выберем на отрезке 0 ≤ x ≤1 сетку ωh		= {xi = ih, i =			} с шагомh . Напишем уравнение баланса
			0, N1
тепла на отрезке xi −1 2 ≤ x ≤ xi+1 2 ,
	xi −1 2	= 1 (xi−1 + xi )= xi −1 +				h	;
xi+1 2	xi+1 2	2		2
Wi −1 2 −Wi +1 2 − ∫ q(x)T (x)dx + ∫ f (x)dx = 0 .
xi−1 2	xi−1 2
Возьмем простейшую аппроксимацию T = const = Ti				при xi −1 2 ≤ x ≤ xi+1 2 :

xi+1 2

∫ q (x)T (x)dx ≈ hdiTi ,

xi−1 2

xi+1 2

di = 1h ∫q(x)dx .

xi−1 2

Проинтегрируем равенство

dT

= −

W

на отрезке xi −1 ≤ x ≤ xi

dx

k

x

W

Ti−1

−Ti = ∫i

dx .

xi −1

k

Полагая W = const =Wi −1 2

при xi −1 ≤ x ≤ xi

, будем иметь

x

dx

Ti −1 −Ti =Wi −1 2

∫i

k(x)

xi−1

или

W

−1 2

= −a

Ti

−Ti−1

;

a

=

1

.

i

i

h

i

1

xi

dx

x∫

h

k(x)

i−1

x

dx

Отметим, что ∫i

есть тепловое сопротивление отрезка [xi−1, xi ]. Заменяя интеграл по одной

k(x)

x

i−1

1

xi

dx

1

1

xi

dx

1

1

2ki

−1ki

из формул

∫

≈

,

∫

≈

+

1

, получим ai = ki −1 2

,

ai =

и т.д.

h

k(x)

h

k(x)

2

ki−1

+ ki

ki −1 2

xi −1

ki−1

ki

xi−1

В результате получим разностную схему вида: