Добавил:

Upload Опубликованный материал нарушает ваши авторские права? Сообщите нам.

Вуз:

Московский государственный областной университет

Предмет:

[НЕСОРТИРОВАННОЕ]

Файл:

Heat_conduction_problems

.pdf

Скачиваний:

Добавлен:

24.03.2015

Размер:

4.71 Mб

Скачать

☆

<<< < Предыдущая 1 2 3 4 5 6 7 8 9 1011 / 1811 12 13 14 15 16 17 18 > Следующая >>>

Rewrite(g);

for i:=1 to N do

writeln(g,' ',h*(i-1):6:3,' ',T[i]-273:8:5); close(g);

close(f1);

end.

Результаты, полученные на основе чисто неявной схемы, полностью совпадают с результатами, полученными по явно-неявной схеме, число итераций при этом не превышает 3.

3.2. ОДНОМЕРНОЕ УРАВНЕНИЕ ТЕПЛОПРОВОДНОСТИ С НЕЛИНЕЙНЫМИ ГРАНИЧНЫМИ УСЛОВИЯМИ (ИЗЛУЧЕНИЕ НА ГРАНИЦЕ)

В качестве примера рассмотрим теплоперенос в бесконечной пластине. Тем самым пренебрегаются два направления переноса тепла, и анализируется одномерное уравнение теплопроводности. На границах области решения будет моделироваться теплообмен за счет конвекции и излучения. Теплоперенос излучением будем рассматривать на основе закона Стефана-Больцмана. Таким образом, сформулированная физическая постановка математически будет выглядеть так:

				ρc ∂T			= λ∂2T , 0 < x < L ;	(44)
						∂t	∂x2
		t = 0 : T =T0 , 0 ≤ x ≤ L;
		x = 0 : −λ		∂T		= κ1 (T e1 −T )+ ε1σ((T e1 )4 −T 4 ), t > 0, κ1 > 0;		(45)
		x = 0 : −λ				= κ1 (T e1 −T )+ ε1σ((T e1 )4 −T 4 ), t > 0, κ1 > 0;		(45)
				∂x
		x = L : λ	∂T		= κ2 (T e2 −T )+ ε2σ((T e2 )4 −T 4 ), t > 0, κ2 > 0;
		x = L : λ			= κ2 (T e2 −T )+ ε2σ((T e2 )4 −T 4 ), t > 0, κ2 > 0;
			∂x
где ε , ε	2	– приведенная степень черноты, σ =5.669 10−8 Вт (м2 К4 ) –
1	2
постоянная Стефана-Больцмана.
Основной интерес в сформулированной краевой задаче
представляют нелинейные граничные условия.
Проведем дискретизацию нелинейных граничных условий III рода
с погрешностью O(h) .
Определим			первые				прогоночные коэффициенты α1 иβ1	из

соотношения T1 = α1 T2 + β1 .

Итак, из второго соотношения (45) следует:

101

− λ

∂T

= κ1 (T e1 −T

x=0 )+ ε1σ((T e1 )4 − (T

x=0 )4 );

∂x

x=0

− λ

T2 −T1

= κ1 (T e1 −T1 )+ ε1σ((T e1 )4 − (T1 )4 ).

κ1 h

Введем обозначение

≡ Bi , тогда

ε1σh

(T e1 )4

ε1δh

−T = Bi

T e1 − Bi

−

(T )4

;

);

Bi1

ε1σh

((T

T1 =

T2 +

)

− (T1 )

1 + Bi

λ(1 + Bi

)

α1

;

1+ Bi1

Bi1

ε1σh

((T e1 )4 −(T1 )4 )

β1 =

T e1 +

1+ Bi1

λ(1+ Bi1 )

или

α1

;

λ + h κ1

(46)

h κ1

ε1σh

((T e1 )4 −(T1 )4 ).

β1 =

T e1 +

λ + h κ1

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

не будет

выполняться условие, вида:

s+1

s – номер

−T1

≤ ε , где

итерации,

ε – точность вычислений.

Правое граничное условие используют для определения

температуры TN .

= κ2 (T e2 −T

x=L )+ ε2σ((T e2 )4 − (T

x=L )4 );

∂T

∂x

x=L

− (TN )4 );

TN −TN −1

= κ2 (T e2 −TN )+ ε2σ((T e2 )4

т.к. TN −1 = αN −1 TN +βN −1 , то

TN − αN −1

TN −βN −1 = Bi2 (T e2 −TN )+

ε2σh

((T e2 )4 − (TN )4 );

102

		βN −1	+ Bi2	T e2 +	ε2σh	((T e2 )4	− (TN )4 )
		βN −1	+ Bi2	T e2 +		((T e2 )4	− (TN )4 )
	TN =				λ
	TN =			1 + Bi2 − αN −1
				1 + Bi2 − αN −1
				или
TN =	λ βN −1 + h κ2 T e2 + ε2σh((T e2 )4 − (TN )4 ).							(47)
			h κ2 + λ (1 − αN −1 )

В результате получили нелинейное уравнение (47) для определения температуры на правой границе. Это уравнение также можно решить наиболее простым методом – методом простых итераций.

Проведем дискретизацию нелинейных граничных условий (45) с погрешностью O(h2 ) . Предположим, что на границе выполняется уравнение теплопроводности (44). Разложим функцию Т(x) в ряд Тейлора в окрестности точки x = 0 до членов второго порядка

относительно h: T n+1	=T n+1	+ h	∂T	n+1	+ h2	∂2T	n+1	. Используя



2	1		∂x	x=0	2	∂x2
							x=0

соотношение (44) получим:

T n+1 =T n+1 + h

∂T

n+1

ρсh2

∂T

∂x

x=0

2 λ

∂t

∂T

n+1

T n+1

−T n+1

ρсh

∂T

n+1

−

∂x

x=0

2 λ

∂t

x=0

n+1

;

x=0

	T n+1		−T n+1		ρсh	T n+1	−T n
=		2	1	−		1	1	.
			h		2 λ
							τ

Из второго соотношения (45):

((T e1 )4 − (T1n+1 )4 ).

∂T

n+1

= κ1 (T1n+1

−T e1 )−

ε1σ

∂x

x=0

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

T2n+1 −T1n+1

−

ρсh

T1n+1 −T1n

κ1

(T1n+1 −T e1 )−

ε1σ((T e1 )4 − (T1n+1 )4 );

2 λ

n+1

ρсh2

n+1

ρсh2

κ h

n+1

κ h

−T

−

2 λ τ

λ τ

−

ε1σh

((T e1 )4 −(T1n+1 )4 );

n+1

ρсh2

κ h

n+1

ρсh2

κ h

ε σh

e1 4

n+1 4

=T2

T1 +

T +

((T ) −(T1 ) );

2 λ τ

103

T n+1 =

2 λ τ

T n+1

ρсh2

T n

ρсh2 + 2 τ(λ + κ h)

ρсh2 + 2 τ(λ + κ h) 2

2τκ h

2τε σh

e1 4

n+1 4

((T

) −(T1

) ).

ρсh2 + 2 τ(λ + κ h)

Таким образом,

2 λ τ

α1 =

;

ρсh

+ 2

τ(λ + κ1h)

ρсh2

2τκ1h

β1 =

(48)

ρсh2 + 2 τ(λ + κ h)

2τε σh

n+1 4

+ ρсh2 + 2 τ(λ + κ h)((T

)

−(T1

)

Определим TN , используя правое граничное условие.

T n+1 =T n+1

− h

∂T

n+1

∂2T

n+1

=T n+1 − h ∂T

n+1

+ ρсh2

∂T

n+1 ;

N −1

∂x

x=L

∂x2

∂x

x=L

2 λ

∂t

x=L

TNn−+11 =TNn+1 − κ2h (T e2 −TNn+1 )− ε2σh ((T e2 )4 − (TNn+1 )4 )+ ρсh2

TNn+1 −TNn

2 λ

Таким образом, т.к. TN −1 = αN −1 TN +βN −1

T n+1 +β

=T n+1

−

κ2h T e2

κ2h T n+1

ρсh2

T n+1

−

ρсh2

T n

N −1

2 λ τ

2 λ τ N

−

ε2σh ((T e2 )4 −(TNn+1 )4 );

n+1

ρсh2

TN 1−αN −1 +

2 λ τ

=βN −1 +

2 λ τ

+ ε2σh

((T e2 )4 −(TNn+1 )4 );

2τκ2hTe2

n+1

2λτβN−1

2λτ(1−α

)+2τκ h+ρсh2

2λτ(1−α

)+2τκ h+ρсh2

N−1

(49)

ρсh2

2τε σh

n+1

+2λτ(1−α )+2τκ h+ρсh2 TN

+2λτ(1−α )+2τκ h+ρсh2 ((T

)

−(TN

)

N−1

Определим температурное поле в бетонной пластине через 600,

1800, 3600 и 7200 секунд. Толщина пластины

L = 0.3 м. Начальная

температура

=500 С.

Бетон имеет

следующие

теплофизические

104

характеристики λ = 0,9 Вт/(м ºC), ρ = 2000 кг/м3, с = 920 Дж/(кг ºC). На

границе	x = 0 и x = L пластина	контактирует	с		внешней средой
( κ =1000 Вт (м2 0 С), T e1 =30 0C ,		ε = 0.5 и	κ	2	=500 Вт (м2 0 С),
1		1		2
T e1 = 70	0C, ε2 = 0.2 ).

Блок-схема к рассматриваемой задаче имеет вид:

105

106

Ниже приведен листинг программы для решения рассматриваемой задачи (на языке программирования Pascal)

uses crt; const mf=500;

sigma=5.669e-8; eps=1e-5;

type

vector=array[1..mf] of real;

var {раздел описания переменных, которые мы будем использовать в

программе}
i, j, N	: integer;
T, Tn, alfa, beta	: vector;
ai, bi, ci, fi, d, d1	: real;
lamda, ro, c, h, tau	: real;
kapa1, kapa2, Te1, Te2	: real;
eps1, eps2	: real;
T0, L, t_end, time	: real;
f, g	: text;
begin
clrscr;

{с клавиатуры вводим все необходимые входные параметры}

Writeln('Введите количество пространственных узлов, N'); Readln(N);

Writeln('Введите окончание по времени, t_end'); Readln(t_end);

Writeln('Введите толщину пластины, L'); Readln(L);

Writeln('Введите коэффициент теплопроводности материала пластины, lamda');

Readln(lamda);

Writeln('Введите плотность материала пластины, ro'); Readln(ro);

Writeln('Введите теплоемкость материала пластины, c'); Readln(c);

Writeln('Введите коэффициент теплообмена на границе х = 0, kapa1'); Readln(kapa1);

Writeln('Введите коэффициент теплообмена на границе х = L, kapa2'); Readln(kapa2);

Writeln('Введите температуру внешней среды относительно границы х

= 0, Te1'); Readln(Te1);

107

Writeln('Введите температуру внешней среды относительно границы х

= L, Te2'); Readln(Te2);

Writeln('Введите приведенную степень черноты на границе х = 0, eps1'); Readln(eps1);

Writeln('Введите приведенную степень черноты на границе х = L, eps2'); Readln(eps2);

Writeln('Введите начальную температуру, T0'); Readln(T0);

{определяем расчетный шаг сетки по пространственной координате} h:=L/(N-1);

{определяем расчетный шаг сетки по времени} tau:=t_end/100.0;

{определяем поле температуры в начальный момент времени} for i:= 1 to N do

T[i]:=T0;

{проводим интегрирование нестационарного уравнения теплопроводности}

time:=0;

while time<t_end do {используем цикл с предусловием} begin

{увеличиваем переменную времени на шаг τ} time:=time+tau;

{определяем alfa начальный прогоночный коэффициент на основе левого граничного условия, используя соотношение (48)} alfa[1]:=2.0*tau*lamda/(2.0*tau*(lamda+kapa1*h)+ro*c*sqr(h));

{запоминаем поле температуры на предыдущем временном слое} for i:=1 to N do

Tn[i]:=T[i];

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

repeat

{определяем beta начальный прогоночный коэффициент на основе левого граничного условия, используя соотношение (48), при этом начинаем итерационный цикл по левому граничному условию}

d:=T[1];

beta[1]:=(ro*c*sqr(h)*Tn[1]+2.0*tau*kapa1*h*Te1+2.0*tau*eps1 *sigma*h*(sqr(sqr(Te1))-sqr(sqr(d))))/(2.0*tau*(lamda +kapa1*h)+ro*c*sqr(h));

108

{цикл с параметром для определения прогоночных коэффициентов по формуле (8)}

for i:= 2 to N-1 do begin

{ai, bi, ci, fi – коэффициенты канонического представления системы уравнений с трехдиагональной матрицей}

ai:=lamda/sqr(h);

bi:=2.0*lamda/sqr(h)+ro*c/tau;

ci:=lamda/sqr(h); fi:=-ro*c*Tn[i]/tau;

{alfa[i], beta[i] – прогоночные коэффициенты} alfa[i]:=ai/(bi-ci*alfa[i-1]); beta[i]:=(ci*beta[i-1]-fi)/(bi-ci*alfa[i-1]); end;

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

repeat d1:=T[N];

{определяем значение температуры на правой границе на основе правого граничного условия, используя соотношение (49)}

T[N]:=(ro*c*sqr(h)*Tn[N]+2.0*tau*(lamda*beta[N-1]+kapa2*h*Te2 +eps2*sigma*h*(sqr(sqr(Te2))-sqr(sqr(d1))))) /(ro*c*sqr(h)+2.0*tau*(lamda*(1-alfa[N-1])+kapa2*h));

until abs(d1-T[N])<=eps; {значение температуры справа определили} {используя соотношение (7) определяем неизвестное поле температуры}

for i:= N-1 downto 1 do T[i]:=alfa[i]*T[i+1]+beta[i];

until abs(d-T[1])<=eps; {значение температуры справа определили} end; {цикл с предусловием окончен}

{выводим результат в файл}

Assign(f,'res.txt');

Rewrite(f);

Writeln(f,'Толщина пластины L = ',L:6:4); Writeln(f,'Число узлов по координате N = ',N);

Writeln(f,'Коэффициент теплопроводности материала пластины lamda = ',lamda:6:4);

Writeln(f,'Плотность материала пластины ro = ',ro:6:4); Writeln(f,'Теплоемкость материала пластины с = ',c:6:4); Writeln(f,'Начальная температура T0 = ',T0:6:4);

109

Writeln(f,'Коэффициент теплообмена kapa1 = ',kapa1:6:4); Writeln(f,'Коэффициент теплообмена kapa2 = ',kapa2:6:4); Writeln(f,'Температура внешней среды Te1 = ',Te1:6:4); Writeln(f,'Температура внешней среды Te2 = ',Te2:6:4); Writeln(f,'Приведенная степень черноты eps1 = ',eps1:6:4); Writeln(f,'Приведенная степень черноты eps2 = ',eps2:6:4); Writeln(f,'Результат получен с шагом по координате h = ',h:6:4); Writeln(f,'Результат получен с шагом по времени tau = ',tau:6:4); Writeln(f,'Температурное поле в момент времени t = ',t_end:6:4); close(f);

Assign(g,'tempr.txt');

Rewrite(g);

for i:=1 to N do

writeln(g,' ',h*(i-1):10:8,' ',T[i]:8:5); close(g);

end.

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

Рис. 23. Распределения температуры по толщине пластины в различные моменты времени

110

<<< < Предыдущая 1 2 3 4 5 6 7 8 9 1011 / 1811 12 13 14 15 16 17 18 > Следующая >>>

Соседние файлы в предмете [НЕСОРТИРОВАННОЕ]

#
08.03.2016139.78 Кб10Frantsuzskaya_revolyutsia_Prichiny_etapy_itogi.doc
#
19.11.2019638.98 Кб3gl_6.doc
#
17.04.2019712.7 Кб18gotovye_otvety_po_filosofii.doc
#
24.09.2019373.76 Кб8GP_OTVET .doc
#
08.03.2016143.36 Кб32Grazhdanskaya_voyna_v_SShA.doc
#
24.03.20154.71 Mб21Heat_conduction_problems.pdf
#
17.08.201979.36 Кб4HELPFUL HINTS.doc
#
14.04.2019280.58 Кб13history_of_psychology_answers_2008 (1).doc
#
24.08.20192.48 Mб25Hrestomatiya_po_filosofii_A.R._Abdullin_Ufa-200...doc
#
21.09.2019118.27 Кб0idioms.doc
#
08.03.2016124.46 Кб44ie_prarodiny_chistovik.docx