Численное решение модельного уравнения
Численное решение модельного уравнения диссипации,
конвекции и кинетики
СОДЕРЖАНИЕ
1.
Общая постановка задачи
2.
Постановка тестовых задач
3.
Методика решения тестовых задач
4.
Результаты вычислений
Список литературы
Приложения
Приложение 1: Описание программы
Приложение 2: Текст
программы
1. ОБЩАЯ ПОСТАНОВКА ЗАДАЧИ
Перенос
тепла (или вещества) теплопроводностью (для вещества соответственно диффузией)
и конвекцией описывается дифференциальным уравнением параболического типа:
( 1 )
где температура
(или концентрация). Пусть являются некоторыми
константами и . Уравнение (1) при указанных
выше предположениях называется модельным уравнением диссипации, конвекции и
кинетики. Слагаемые правой части имеют следующий физический смысл:
-
соответствует переносу тепла теплопроводностью (или вещества диффузией);
-
соответствует конвективному переносу;-
-
"кинетический член", соответствует источнику, пропорционально-
му
температуре или концентрации;
-
интенсивность внешних источников или стоков.
В дальнейшем будем
рассматривать только тепловую интерпретацию уравнения (1).
Численное решение
уравнения (1) будем искать в области :
( 2 )
при заданных начальных
значениях температуры: (
3 )
и граничных условиях.
Граничные условия
описывают режимы теплообмена с внешней средой:
при
;
при .
2. ПОСТАНОВКА ТЕСТОВЫХ ЗАДАЧ
В
качестве тестовых задач для температуры мною
были выбраны следующие пять функций:
(
9 )
( 10 )
( 11 )
( 12 )
( 13 )
Для функции (9) имеем:
Для функции (10):
Для функции (11):
Для функции (12):
Для функции (13):
Данные функции тестировались на отрезке по X: [0,
1], по времени: [0, 1], с количеством разбиений по этим отрезкам - 30.
3. МЕТОДИКА РЕШЕНИЯ ТЕСТОВЫХ ЗАДАЧ
Данная задача решается с
помощью двухслойной неявно конечно-разностной схемы.
Схема
реализуется в три этапа.
1 этап: находятся предварительные значения с помощью 4-х точечной неявной схемы:
( 5 )
2 этап: используется за два шага. Сначала находятся на полученном слое () с шагом , а
затем через . В
этом случае используется 4-х точечная неявная разностная схема:
( 6 )
( 7 )
3 этап: окончательные значения находятся
в виде линейной комбинации двух предварительных значений:
( 8 )
Для
решения (1) воспользуемся формулами (5) - (8). Данные уравнения представляют
трех диагональные матрицы, решаемые методом скалярной прогонки.
В начале нужно
преобразовать (5) – (7) к виду:
(
14 )
Тогда (5) примет вид:
;
;
.
Формула (6) преобразуется
в:
Т.е.
;
;
;
.
Формула (7) преобразуется в:
Т.е. ;
;
;
.
Далее решаем по формулам
скалярной прогонки:
( 15 )
( 16 )
Для определения , и воспользуемся данными граничными
условиями, т.е. формулой (4) и функцией . Так
если мы берём из формулы (9), то имеем:
Приведём это выражение к виду: .
Т.е. теперь мы имеем
и
:
Далее найдем конечное :
( 18 )
Проведя
аналогичные расчёты для заданных формулами
(10) – (13), мы получим соответствующие , и .
Далее мы можем решить системы методом прогонки и получить требуемый результат.
4. РЕЗУЛЬТАТЫ ВЫЧИСЛЕНИЙ
В
результате проведённых испытаний программа показала свою высокую надёжность.
Были получены следующие данные.
При расчёте с
использованием функции и входных данных ; ; ; ; ; ; на отрезке по X и по времени [0,1]
с шагом 0,033 был получен результат с ошибкой равной 0,0675.
Для функции при ; ; ; ; ; ; , на
том же промежутке, ошибка составляет 0,055.
С функцией и ; ; ; ; ; ; ошибка
примет значение 0,0435.
При и
условиях ; ; ; ; ; ; в результате возникает ошибка равная 0,0055.
И, наконец, если выбрана
функция и ; ; ; ; ; ; , то
ошибка составит 0,00255.
Т.е. можно сказать, что мы
имеем результат с первым порядком точности. Столь малую точность можно
объяснить тем, что производная, найденная при граничных условиях, так же имеет
первый порядок точности.
СПИСОК ЛИТЕРАТУРЫ
1.
А. Епанешников, В. Епанешников Программирование в среде Turbo-Pascal
7.0. - М.: Диалог - Мифи, 1996. - 288 с.
2.
Петухова Т. П., Сибирцев В.
В. Пакет прикладных программ для
численного моделирования процессов тепло- и массопереноса. – Караганда: Изд-во
КарГУ. 1993
3.
Фигурнов В. Э. IBM PC для пользователя. - М.: Инфра
- М, 1995. - 432 с.
Приложение 1
ОПИСАНИЕ ПРОГРАММЫ
Поставленная
задача была программно реализована на языке программирования Turbo-Pascal 7.0.
В состав
программы входят следующие файлы:
basis.pas - PAS-файл
основной части программы
(решение системы уравнений методом скалярной прогонки);
basis.v&v - EXE-файл
основной части программы (вызывается из START.PAS);
fun.bmp - BMP-фаил с
изображением функций;
inform.v&v - TXT-фаил с
информацией о программе (вызывается из START.PAS);
music.v&v - музыкальный
EXE-фаил (вызывается из START.PAS);
my_menu.pas - UNIT для
создания меню;
sea.exe - программа для
просмотра графических файлов;
start.pas - файл для
запуска всей программы;
u - файл с
результатами работы;
zastavka.v&v - EXE-фаил с
заставкой к основной программе
(вызывается из START.PAS).
Файл
START является, как
бы оболочкой программы, из которой вызываются другие файлы. Сам процесс решения
содержится в файле BASIS.
BASIS содержит
следующие процедуры и функции:
Function Fun_U (Xm,t:real):real;
Вход:
значение по X и значение по
времени t, а также
глобальная переменная выбранной
функции SelectFunction.
Действие:
вычисляет точное значение функции U при заданных X и t.
Выход:
Fun_U – значение
функции.
Function Fun_F (Xm,t,a,b,v:real):real;
Вход:
значение по X, по времени t, коэффициенты
, , и номер выбранной функции
Действие:
вычисляет значение функции F при заданных X, t, , , .
Выход:
Fun_F – значение
функции F.
Function Betta_Zero (time:real):
real;
Вход:
значение времени t и глобальные коэффициенты , , , номер выбранной
функции SelectFunction.
Действие:
вычисляет , используемое в методе скалярной
прогонки.
Выход:
Betta_Zero – значение .
Function U_End
(time,Alf,Bet:real): real;
Вход:
значение времени t, , и глобальные коэффициенты ,
, ,
номер выбран-
ной
функции
SelectFunction.
Действие:
вычисляет используемое в методе скалярной
прогонки.
Выход:
U_End – значение .
Procedure
PrintArray;
Вход:
использует глобальный массив данных U_m.
Действие:
выдает содержимое U_m на экран и в файл.
Выход:
вывод U_m.
Приложение 2
ТЕКСТ ПРОГРАММ Ы
Основная
часть программы выглядит так:
Program Basis;
Uses Crt; { Подключение библиотек }
Label
Metka1,Metka2; { Метки }
Var
a, b, v : real; { Коэффициенты,
задаются пользователем }
h, tau : real; { Шаг по X и по
времени соответственно }
X,x0 : real; { Конечное и
начальное значение X }
m,n,k : word; { Переменные
используемые в циклах для расчета }
T,t0 : real; { Конечное и
начальное значение времени }
Kol_voX, Kol_voT : word; {
Количество разбиений по X и по времени }
U_m,U_,_U_1_2,_U_1 : array [0..200] of
real; { Массивы результатов }
z : array [0..200] of real; {
Массив точных решений }
Xm : real; { Промежуточный X }
Alfa,Betta : array [0..200] of
real; { Массив коэффициентов используемых при скалярной прогонке }
a_progonka, b_progonka, c_progonka,
d_progonka : real; { Коэффициенты для скалярной прогонки }
Error : real; { Значение ошибки }
time : real; { Переменная
времени }
ch : char; { Код нажатой
клавиши }
SelectFunction:word; { Номер
выбранной функции }
U : text; { Переменная для
вывода результата в файл }
Alfa_1,Alfa_2,Betta_1,Betta_2 :
real; { Коэффициенты граничных условий }
Data : word; { Переменная режима
ввода начальных данных }
Function Fun_U (Xm,t:real):real; {
Функция U (точное решение) }
begin
If
SelectFunction=1 then Fun_U:=SQR(Xm)*Xm+SQR(t);
If
SelectFunction=2 then Fun_U:=SQR(Xm)*SQR(t)*t+10*Xm*t+SQR(SQR(t))*Xm;
If
SelectFunction=3 then Fun_U:=Xm*SIN(Xm*t)-4*SQR(Xm)*COS(t);
If
SelectFunction=4 then Fun_U:=t*EXP(Xm);
If
SelectFunction=5 then Fun_U:=SIN(Xm)+EXP(t);
end;
Function Fun_F
(Xm,t,a,b,v:real):real; { Функция F }
begin
if
SelectFunction=1 then Fun_F:=2*t-v*6*Xm+a*3*SQR(Xm)-b*(SQR(Xm)*Xm+SQR(t));
if
SelectFunction=2 then
Fun_F:=3*SQR(Xm)*SQR(t)+10*Xm+4*SQR(t)*t*Xm-v*2*SQR(t)*t+
a*(2*Xm*SQR(t)*t+10*t+SQR(SQR(t)))-b*(SQR(Xm)*SQR(t)*t+10*Xm*t+Xm*SQR(SQR(t)));
if
SelectFunction=3 then
Fun_F:=SQR(Xm)*COS(Xm*t)+4*SQR(Xm)*SIN(t)-v*(2*COS(Xm*t)*t-
Xm*SIN(Xm*t)*SQR(t)-8*COS(t))+a*(SIN(Xm*t)+Xm*t*COS(Xm*t)-8*COS(t)*Xm)-
b*(Xm*SIN(Xm*t)-4*SQR(Xm)*COS(t));
if
SelectFunction=4 then
Fun_F:=EXP(Xm)-v*(t*EXP(Xm))+a*(t*EXP(Xm))-b*(t*EXP(Xm));
if
SelectFunction=5 then
Fun_F:=EXP(t)-v*(-SIN(Xm))+a*(COS(Xm))-b*(SIN(Xm)+EXP(t));
end;
Function
Betta_Zero (time:real): real; { Функция Betta[0] для прогонки }
begin
If SelectFunction=1
then Betta_Zero:=(h/(Betta_1*h-Alfa_1))*(Alfa_1*3*SQR(x0)+
Betta_1*(SQR(x0)*x0+SQR(time)));
If
SelectFunction=2 then
Betta_Zero:=(h/(Betta_1*h-Alfa_1))*(Alfa_1*(2*x0*SQR(time)*time+
10*time+SQR(SQR(time)))+Betta_1*(SQR(x0)*SQR(time)*time+10*x0*time+SQR(SQR(time))*x0));
If
SelectFunction=3 then
Betta_Zero:=(h/(Betta_1*h-Alfa_1))*(Alfa_1*(SIN(x0*time)+
x0*time*COS(x0*time)-8*x0*COS(time))+Betta_1*(x0*SIN(x0*time)-4*SQR(x0)*COS(time)));
If
SelectFunction=4 then
Betta_Zero:=(h/(Betta_1*h-Alfa_1))*(Alfa_1*(time*EXP(x0))+
Betta_1*(time*EXP(x0)));
Betta_1*(SIN(x0)+EXP(time)));
end;
Function U_End
(time,Alf,Bet:real): real; { Функция Um для прогонки }
begin
If
SelectFunction=1 then U_End:=(Alfa_2*h*3*SQR(X)+Betta_2*h*(SQR(X)*X+SQR(time))
+
Bet*Alfa_2)/(Alfa_2-Alf*Alfa_2+h*Betta_2);
If
SelectFunction=2 then U_End:=(Alfa_2*h*(2*X*SQR(time)*time+10*time+SQR(SQR(time)))+
Betta_2*h*(SQR(X)*SQR(time)*time+10*X*time+SQR(SQR(time))*X)
+Bet*Alfa_2)/(Alfa_2-Alf*Alfa_2+h*Betta_2);
If
SelectFunction=3 then
U_End:=(Alfa_2*h*(SIN(X*time)+X*time*COS(X*time)-8*X*COS(time))+
Betta_2*h*(X*SIN(X*time)-4*SQR(X)*COS(time))+Bet*Alfa_2)/(Alfa_2-Alf*Alfa_2+h*Betta_2);
If
SelectFunction=4 then
U_End:=(Alfa_2*h*(time*EXP(X))+Betta_2*h*(time*EXP(X))+Bet*Alfa_2)/
(Alfa_2-Alf*Alfa_2+h*Betta_2);
If SelectFunction=5
then U_End:=(Alfa_2*h*(COS(X))+Betta_2*h*(SIN(X)+EXP(time))+Bet*Alfa_2)/
(Alfa_2-Alf*Alfa_2+h*Betta_2);
end;
Procedure
PrintArray; { Процедура печати массива U }
begin
WriteLn; For
m:=0 to Kol_voX do begin Write(U_m[m]:15:4); Write(U,U_m[m]:15:4); end;
WriteLn;
WriteLn(U);
end;
{ Основная программа }
Begin
Assign(U,'u'); { Файл для записи
значений функции }
Rewrite(U); { Открытие файла для
записи }
TextBackGround(0); { Выбор функции
для работы }
ClrScr; TextColor(10);
GoToXY(20,8); Write('Введите номер выбранной функции (1-5):');
Metka1:
ch:=ReadKey;
If ch='1'
then SelectFunction:=1
else If
ch='2' then SelectFunction:=2
else If ch='3' then SelectFunction:=3
else
If ch='4' then SelectFunction:=4
else If ch='5' then SelectFunction:=5
else
begin
Sound(400); Delay(100); NoSound; GoTo Metka1;
end;
GoToXY(59,8);TextColor(12);WriteLn(SelectFunction); TextColor(11);
GoToXY(11,12);
Write('Вы будете
работать со стандартными параметрами (цифра ~1~)');
GoToXY(22,13); Write('или введете свои
данные (цифра ~2~) ?');
Metka2:
ch:=ReadKey;
If ch='1'
then Data:=1
else If
ch='2' then Data:=2
else
begin
Sound(400); Delay(100); NoSound; GoTo Metka2;
end;
TextBackGround(9); TextColor(10); ClrScr;
{ Ввод начальных
данных }
WriteLn;
WriteLn('-------------------------------- Ввод данных
---------------------------------¬');
For k:=1 do 21
do WriteLn('¦ ¦');
WriteLn('L------------------------------------------------------------------------------');
TextColor(15);
Window(3,3,77,23); Write(' Введите область рассчета по X от: ');
If Data=1 then
begin
x0:=0; Write(x0:1:0); WriteLn;
end
else
ReadLn(x0);
Write(' до: ');
If Data=1 then
begin
X:=1; Write(X:1:0); WriteLn;
end
else
ReadLn(X);
WriteLn;
Write(' Введите количество разбиений по направлению X: ');
If Data=1 then
begin Kol_voX:=30; Write(Kol_voX:2); WriteLn; end else ReadLn(Kol_voX);
WriteLn;WriteLn;
Write(' Введите область рассчета по времени от: ');
If Data=1 then
begin t0:=0; Write(t0:1:0); WriteLn; end else ReadLn(t0);
Write(' до: ');
If Data=1 then
begin T:=1; Write(T:1:0); WriteLn; end else ReadLn(T);
If Data=1 then
begin Kol_voT:=30; Write(Kol_voT:2); WriteLn; end else ReadLn(Kol_voT);
WriteLn;WriteLn; WriteLn(' Введите коэффициенты'); Write('
a=');
If Data=1 then
begin a:=1; Write(a:1:0); WriteLn; end else ReadLn(a);
Write('
b=');
If Data=1 then
begin b:=1; Write(b:1:0); WriteLn; end else ReadLn(b);
Write('
v=');
If Data=1 then
begin v:=0.001; Write(v:1:3); WriteLn; end else ReadLn(v);
Write('
Alfa-1=');
If Data=1 then
begin Alfa_1:=1; Write(Alfa_1:1:0); WriteLn; end else ReadLn(Alfa_1);
Write('
Betta-1=');
If Data=1 then
begin Betta_1:=1; Write(Betta_1:1:0); WriteLn; end else ReadLn(Betta_1);
Write('
Alfa-2=');
If Data=1 then
begin Alfa_2:=1; Write(Alfa_2:1:0); WriteLn; end else ReadLn(Alfa_2);
Write('
Betta-2=');
If Data=1 then
begin Betta_2:=1; Write(Betta_2:1:0); WriteLn;TextColor(14);
Write(' Нажмите любую клавишу');
ReadKey; end else ReadLn(Betta_2);
{ Интерфейс
экрана при выдаче результата }
TextBackGround(3);
TextColor(1); Window(1,1,80,25); ClrScr; WriteLn;
WriteLn('г=====================
Результат ==========================¬');
For k:=1 to 21
do
WriteLn('¦ ¦');
WriteLn('===================================================================-');
TextColor(0);
TextBackGround(7); GoToXY(2,23);
WriteLn(' Для продолжения нажмите любую клавишу');
TextBackGround(3); Window(3,4,77,22);
TextColor(15); ClrScr;
{ Вычесление шага сетки }
tau:=(T-t0)/Kol_voT; h:=(X-x0)/Kol_voX;
{ Ввод данных
при time=t0 }
For m:=0 to
Kol_voX do
begin
Xm:=x0+h*m; U_m[m]:=Fun_U(Xm,t0);
end;
TextColor(14);
WriteLn('Время равно
',time:3:3); TextColor(15); WriteLn(U,'Время равно ',time:3:3);
PrintArray;
{ Начало рассчета }
time:=t0;
Repeat
time:=time+tau;
WriteLn;
TextColor(14); WriteLn('Время равно
',time:3:3); TextColor(15);
WriteLn(U,'Время равно
',time:3:3);
{ 1 этап. Решается методом
скалярной прогонки }
a_progonka:=(-2*v-a*h)/(2*SQR(h));
b_progonka:=(SQR(h)+2*v*tau-b*tau*SQR(h))/(SQR(h)*tau);
c_progonka:=(a*h-2*v)/(2*SQR(h));
Alfa[0]:=Alfa_1/(Alfa_1-Betta_1*h); Betta[0]:=Betta_Zero(time);
For m:=1 to
Kol_voX-1 do
begin
Alfa[m]:=-c_progonka/(a_progonka*Alfa[m-1]+b_progonka);
Betta[m]:=(Fun_F(x0+m*h,time+tau,a,b,v)+U_m[m]/tau-a_progonka*Betta[m-1])/
(a_progonka*Alfa[m-1]+b_progonka);
end;
U_[Kol_voX]:=U_End(time,Alfa[Kol_voX-1],Betta[Kol_voX-1]);
For
m:=Kol_voX-1 downto 1 do
U_[m]:=Alfa[m]*U_[m+1]+Betta[m];U_[0]:=Alfa[0]*U_[1]+Betta[0];
{ 2 этап, часть
первая. Решается методом скалярной прогонки }
a_progonka:=(-2*v-a*h)/(2*SQR(h));
b_progonka:=(2*SQR(h)+2*v*tau-b*tau*SQR(h))/(SQR(h)*tau);
c_progonka:=(a*h-2*v)/(2*SQR(h));
Alfa[0]:=Alfa_1/(Alfa_1-Betta_1*h); Betta[0]:=Betta_Zero(time);
For m:=1 to
Kol_voX-1 do
begin
Alfa[m]:=-c_progonka/(a_progonka*Alfa[m-1]+b_progonka);
Betta[m]:=(Fun_F(x0+m*h,time+tau/2,a,b,v)+2*U_m[m]/tau-a_progonka*Betta[m-1])/
(a_progonka*Alfa[m-1]+b_progonka);
end;
_U_1_2[Kol_voX]:=U_End(time,Alfa[Kol_voX-1],Betta[Kol_voX-1]);
For m:=Kol_voX-1
downto 1 do _U_1_2[m]:=Alfa[m]*_U_1_2[m+1]+Betta[m];
_U_1_2[0]:=Alfa[0]*_U_1_2[1]+Betta[0];
{ 2 этап, часть
вторая. Решается методом скалярной прогонки }
a_progonka:=(-2*v-a*h)/(2*SQR(h));
b_progonka:=(2*SQR(h)+2*v*tau-b*tau*SQR(h))/(SQR(h)*tau);
c_progonka:=(a*h-2*v)/(2*SQR(h));
Alfa[0]:=Alfa_1/(Alfa_1-Betta_1*h); Betta[0]:=Betta_Zero(time);
For m:=1 to
Kol_voX-1 do
Alfa[m]:=-c_progonka/(a_progonka*Alfa[m-1]+b_progonka);
Betta[m]:=(Fun_F(x0+m*h,time+tau,a,b,v)+2*_U_1_2[m]/tau-a_progonka*Betta[m-1])/
(a_progonka*Alfa[m-1]+b_progonka);
end;
_U_1[Kol_voX]:=U_End(time,Alfa[Kol_voX-1],Betta[Kol_voX-1]);
For
m:=Kol_voX-1 downto 1 do _U_1[m]:=Alfa[m]*_U_1[m+1]+Betta[m];
_U_1[0]:=Alfa[0]*_U_1[1]+Betta[0];
{ 3 этап.
Окончательное значение }
For m:=0 to
Kol_voX do
U_m[m]:=2*_U_1[m]-U_[m];
PrintArray; {
Вывод результата на экран и его запись в файл }
For m:=0 to Kol_voX do { Рассчет
точного значения функции }
begin
z[m]:=Fun_U(x0+m*h,time); end;
{ Вывод ошибки
расчета на экран и в файл }
Error:=0;
For m:=0 to
Kol_voX do
begin
a:=Abs(U_m[m]-z[m]); If Error<a then Error:=a;
end;
WriteLn; TextColor(4);
WriteLn('Максимальная ошибка для этого времени равна ',Error:10:7);
TextColor(15); WriteLn(U,'Максимальная
ошибка для этого времени равна ',Error:15:13);
WriteLn(U);
ReadKey;
Until
time>T;
Close(U); {
Закрытие файла с результатами }
End.