Модель явлений переноса (теплопроводность, диффузия)

1. Задача: Имеется прямоугольная пластина, коэффициент температуропроводности которого зависит от координаты a=a(x,y). Задано начальное распределение температуры T=T(x,y), а также мощность q_n и координаты (x_n, y_n) источников тепла. Необходимо изучить изменение температуры различных точек пластины с течением времени.

2. Теория. Уравнение теплопроводности для пластины:

где q(x,y) -- скорость увеличения температуры, обусловленная наличием в точке с координатами (x,y) источников (поглотителей тепла) тепла.

Приращение температуры равно:

Разобъем прямоугольную пластину на элементы, образующие вертикальные столбцы с номером i и горизонтальные строки с номером j. Пусть в некоторый дискретный момент времени t температура элемента i-ого столбца j-ой строки была равна T_i,j(t). Чтобы найти ее значение в следующий момент времени t+dt необходимо численно решить представленное выше дифференциальное уравнение отдельно по столбцам и отдельно по строкам.

Решая одномерную задачу теплопроводности по столбцам, получаем:

T_i,j(t+dt)=T_i,j(t)+a[(T_i,j+1 -2T_i,j+T_i,j-1)/h²]dt+q_i,jdt,

где q_i,j -- быстрота увеличения температуры элемента вследствие наличия в нем источника тепла, h=1. Аналогичным образом решается уравнение теплопроводности по строкам:

T_i,j(t+dt)=T_i,j(t)+a[(T_i+1,j -2T_i,j+T_i-1,j)/h²]dt+q_i,jdt,

Организовав два цикла по i и j, в которых перебираются все элементы пластины, вложенные в цикл по времени, можно расчитать распределение температуры пластины в последующие моменты времени. Чтобы исключить влияние направления перебора элементов на результаты моделирования, следует организовать 4 цикла, в которых элементы пластины перебираются слева направо, справа налево, снизу вверх, сверху вниз.

3. Алгоритм.

1. Задают коэффициенты температуропроводности a и b вдоль осей координат x и y, начальное распределение температуры T_i,j(0) различных элементов пластины, координаты и мощности источников тепла q_i,j. Считают, что t=0.

2. Путем последовательного перебора элементов i--ого столбца от j=1 до M с помощью конечно-разностного уравнения (1) пересчитывают их температуры в следующий дискретный момент времени. Повторяют расчеты для всех столбцов.

3. Повторяют ту же самую процедуру в противоположном направлении, перебирая элементы от j=M до 1 и используя уравнение (1).

4. Путем последовательного перебора элементов j--ой строки от i=1 до N с помощью конечно-разностного уравнения (2) пересчитывают их температуры в следующий дискретный момент времени. Повторяют расчеты для всех строк.

5. Повторяют ту же самую процедуру в противоположном направлении, перебирая элементы от i=N до 1 и используя уравнение (2).

6. Вводят текущее распределение температуры на экран, закрашивая элементы с различной температурой различными цветами.

7. Возвращение к операции 2. Если цикл по t закончился, -- выход из цикла.

4. Компьютерная реализация. Предлагаемая компьютерная программа позволяет исходя из начального распределения температуры и наличия источников тепла, расчитать температуру различных точек прямоугольной пластины в дискретные моменты времни.

program PROGRAMMA7;uses crt, graph;const n=100; m=100; h=1; dt=0.2;var ii,jj,kk,i,j,DriverVar, ModeVar, ErrorCode : integer;t: array[1..N, 1..M] of real; q,a,b,bb :real; naprav,uslovie: boolean; procedure Init; {---- Инициализация графики ----}begin DriverVar:=Detect; InitGraph(DriverVar,ModeVar,'c:\bp\bgi'); ErrorCode:=GraphResult; if ErrorCode <> grOK then Halt(1); end; procedure Param_sred; {---Коэффициент температуропроводности---}begin if j<70 then a:=2 else a:=1; end; procedure Istoch; {--- Источники тепла ---}begin if ((i>45)and(i<55))and((j>70)and(j<75)) then q:=50 else q:=0; end; procedure Nach_uslov; {-- Начальное распределение температуры --}begin For i:=1 to N do For j:=1 to M do begin uslovie:=((j<65)and(j>45)and(i>20)and(i<30)) or((j<45)and(j>35)and(i>50)and(i<60)); if uslovie=true then t[i,j]:=450 else t[i,j]:=1; end; end; procedure Raschet; {---- Расчет температуры ----}begin Istoch; Param_sred; t[i,j]:=t[i,j]+a*(t[i,j+1]-2*t[i,j]+t[i,j-1])*dt/(h*h)+q; if naprav=true then t[i,j]:= t[i,j]+a*(t[i+1,j]-2*t[i,j]+t[i-1,j])*dt/(h*h); end; procedure Draw;{---- Вывод на экран ----}begin if t[i,j]>50 then setcolor(2); if t[i,j]>300 then setcolor(12); if (t[i,j]<300)and(t[i,j]>120) then setcolor(10); if (t[i,j]<120)and(t[i,j]>70) then setcolor(3); if (t[i,j]<70)and(t[i,j]>30) then setcolor(4); if (t[i,j]<30)and(t[i,j]>20) then setcolor(5); if (t[i,j]<20)and(t[i,j]>10) then setcolor(7); if t[i,j]<10 then setcolor(15); rectangle(i*5+50,j*5,i*5+54,j*5+4); end; BEGIN {---- Основная программа ----}Init; Nach_uslov;Repeat kk:=kk+1;For i:=2 to N-1 do For j:=2 to M-1 do begin naprav:=false; Raschet; end;For j:=2 to M-1 do For i:=2 to N-1 do begin naprav:= true; Raschet; end;For i:=2 to N-1 do For jj:=2 to M-1 do begin j:=M+1-jj; naprav:=true; Raschet; end;For j:=2 to M-1 do For ii:=2 to N-1 do begin i:=N+1-ii; naprav:=false; Raschet; end;if kk/2=round(kk/2) then For i:=2 to N-1 do For j:=2 to M-1 do Draw;until KeyPressed; CloseGraph;END.

Тема 6,7. Динамические системы. Реализация алгоритма для экологических систем