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

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

(Программа написана под руководством профессора кафедры микро- и наноэлектроники СПбГЭТУ "ЛЭТИ" Рындина Е.А.)

Моделируемые данные:

Нестационарное уравнение теплопроводности:

где t — время; x,y — координаты; T(x,y) — искомая функция распределения абсолютной температуры по координатам; ρ(x,y) — плотность вещества; C(x,y) — удельная теплоемкость вещества; k(x,y) — коэффициент теплопроводности вещества; f(x,y) — плотность мощности источников тепла на прямоугольной области с граничными условиями Дирихле или Неймана на границах x=x_min, x=x_max, y=y_min, y=y_max и начальными условиями первого или второго рода на отрезке времени [t_min,t_max].

Координатная сетка G_x,y={(x_i,y_j)| i=1…I; j=1…J} задается в программе массивами L и W и числами шагов Sx и Sy; аналогично задается временная сетка G_t={t_m| m=1…M}. Начальные и граничные условия:

Дискретизируем уравнение, начальные и граничные условия; получаем следующее:

Рис. 1 Двумерное тело

Рис. 2,3 Временной и координатный шаблоны

Выражая температуру в следующем временном слое, зная предыдущий; обозначаем T_i,j,m+1 как T_i,j^*, а T_i,j,m+1 как T_i,j, получаем:

T_i,j = g_t, m = 1;

,

что позволяет рассчитывать температуру в следующем временном слое, зная её в настоящем, в чем и заключается суть явной схемы.

Результаты расчета:

Рис. 4 Изначальное распределение температуры в теле

Рис. 5 Распределение температуры в следующий момент времени dt

Рис. 6 Распределение температуры в теле в следующий момент времени

Рис. 7 Распределение температуры в теле в следующий момент времени

Рис. 8 Распределение температуры в теле в следующий момент времени

Рис. 9 Распределение температуры в теле в следующий момент времени

Рис. 10 Распределение температуры в теле в следующий момент времени

Код программы:

clear all
close all
clc

L=[30 150 230];
L=L.*1e-3;
L(3)=L(3)-L(1)-L(2);

W=[10 10 20 10 50 10 20 30 30 130];
W=W.*1e-3;
%задаем толщины слоев по x и y, переводим в мм

kM=[1000 100 200 100 300 10 10000];
rM=[10000 2000 1000 2000 1000 2000 300];
cM=[10000 100 1000 100 1000 100 300];
%значения коэффициента теплопроводности, плотности и теплоемкости на
%участках; M — массив

gt=300;
gxmin=-50;
gxmax=-50;
gymin=300;
gymax=300;
%значения для начального и граничных условий

tmax=25;
tzero=0.5*tmax;
F=1e7;
GR=1e4;
%максимальное время исследования, время обнуления f, функция плотности мощности в заданном
%участке и частота выведения графиков

St=10; % число отрезков по времени и координатам
Sx=5;
Sy=7;

x(1)=0;
kL(1)=kM(1);
rL(1)=rM(1);
cL(1)=cM(1);
for i=1:length(kM)
    x=[x max(x)+W(i)/Sx:W(i)/Sx:max(x)+W(i)]; % формируем сетки — каждый слой разбиваем на заданное число шагов, присоединяем к матрице x
    kL=[kL ones(1,Sx).*kM(i)]; %единичные матрицы, так как постоянные внутри слоя; ones(число строк, число столбцов)
    rL=[rL ones(1,Sx).*rM(i)]; 
    cL=[cL ones(1,Sx).*cM(i)];
end
I=length(x);
dx=diff(x);
kL=kL'; %транспонируем в столбцы
rL=rL';
cL=cL';

y(1)=0; %сетка по y
for i=1:length(L)
    y=[y max(y)+L(i)/Sy:L(i)/Sy:max(y)+L(i)];
end
J=length(y);
dy=diff(y);

k=kL;
r=rL;
c=cL;
for j=2:J
    k=[k kL];
    r=[r rL];
    c=[c cL];
end

dt = 1e-7;
t(1) = 0;
% Tpast

for j=1:J
    for i=1:I
        Tp(i,j) = gt; % начальное условие
    end
end

mesh(y.*1e3,x.*1e3,Tp-273) % первый график
xlabel('y, mm','FontSize',19)
ylabel('x, mm','FontSize',19)
zlabel('T, ^oC','FontSize',19)
grid on
xlim([min(y.*1e3) max(y.*1e3)])
ylim([min(x.*1e3) max(x.*1e3)])
colormap([0 0 0])
set(gcf,'Position',[680 42 1241 954]);
print(gcf,'-djpeg','Figure_1')
pause(1e-3)

while max(t)<tmax
    t = [t max(t)+dt];
    M=length(t);
    f=zeros(I,J);
    for j=Sy+1:2*Sy+1 %ненулевая функция внутри заданного прямоугольника, а не везде
        for i=1:I
            if x(i)>=W(9) && x(i)<=W(9)+W(8) && max(t)<tzero
                f(i,j)=F;
            end
        end
    end
    
    for j=1:J % граничные условия, согласно уравнениям
        T(1,j) = Tp(2,j) + dx(1)*gxmin;
        T(I,j) = Tp(I-1,j) + dx(I-1)*gxmax;
    end
    
    for i=2:I-1
        T(i,1) = gymin;
        T(i,J) = gymax;
    end
        
    for j=2:J-1
        for i=2:I-1 % температура в следующем временном слое, согласно уравнению
            T(i,j) = Tp(i,j) + dt/r(i,j)/c(i,j)*...
            (2/(dx(i)+dx(i-1))*((k(i+1,j)+k(i,j))/2/dx(i)*(Tp(i+1,j)-Tp(i,j)) - (k(i,j)+k(i-1,j))/2/dx(i-1)*(Tp(i,j)-Tp(i-1,j)))+...
            2/(dy(j)+dy(j-1))*((k(i,j+1)+k(i,j))/2/dy(j)*(Tp(i,j+1)-Tp(i,j)) - (k(i,j)+k(i,j-1))/2/dy(j-1)*(Tp(i,j)-Tp(i,j-1)))+...
            f(i,j));
        
        end
    end
    
    NN='Figure_';
    if M/GR-fix(M/GR) == 0
        NNN=[NN num2str(M)];
        figure
        mesh(y.*1e3,x.*1e3,T-273)
        xlabel('y, mm','FontSize',19)
        ylabel('x, mm','FontSize',19)
        zlabel('T, ^oC','FontSize',19)
        grid on
        xlim([min(y.*1e3) max(y.*1e3)])
        ylim([min(x.*1e3) max(x.*1e3)])
        colormap([0 0 0])
        set(gcf,'Position',[680 42 1241 954]);
        print(gcf,'-djpeg',NNN)
        pause(1e-3)
    end
    
    Tp = T; % делаем настоящий момент времени предыдущим, переходя к следующему слою
end

В заключение, отмечена необходимость достаточно малого шага по времени для устойчивости схемы, что ощутимо увеличивает время расчетов, по сравнению, например, с расчетами прямым методом решения СЛАУ. Программа проста в организации, подходит для работы на системах с малым количеством ресурсов, сразу выводит графики в ходе расчета.