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

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

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

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

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

где 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 как T_i,j^*, а T_i,j,m-1 как T_i,j, получаем:

то есть T_i,j^* неявно зависит от T_i,j.

На каждом временном шаге решается СЛАУ путем создания соответствующих матриц коэффициентов и свободных членов и решения матричного уравнения Ax=B => x=A^-1B.

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

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

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

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

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

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

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

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

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

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

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

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

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

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=10;
tzero = 0.5*tmax;
F=1e7;
GR = 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)];
    kL=[kL ones(1,Sx).*kM(i)];
    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;
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-1;
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
    
    %%% не удаляемˆ; циклы без изменения, вместо уравнений, записать
    %%% коэффициенты и свободные члены
    % I*(j-1)+i
    A=zeros(I*J, I*J);
    B=zeros(I*J, 1);
    
    for j=1:J
        A(I*(j-1)+1, I*(j-1)+1)=-1/dx(1);
        A(I*(j-1)+1, I*(j-1)+2)=1/dx(1);
        % первый индекс такой же; A и B — индексы при переменных в граничных условиях
        B(I*(j-1)+1)=gxmin;
        
        A(I*(j-1)+I, I*(j-1)+I)=1/dx(I-1);
        A(I*(j-1)+I, I*(j-1)+I-1)=-1/dx(I-1);
        B(I*(j-1)+I)=gxmax;
    end
    
    for i=2:I-1
        A(I*(1-1)+i, I*(1-1)+i) = 1;
        B(I*(1-1)+i) = gymin;
        A(I*(J-1)+i, I*(J-1)+i) = 1;
        B(I*(J-1)+i) = gymax;
    end
        
    for j=2:J-1
        for i=2:I-1
           A(I*(j-1)+i, I*(j-1)+i)=r(i, j)*c(i,j)/dt+...
               2/(dx(i)+dx(i-1))*((k(i+1, j)+k(i, j))/2/dx(i)+(k(i, j)+k(i-1, j))/2/dx(i-1))+...
               2/(dy(j)+dy(j-1))*((k(i, j+1)+k(i, j))/2/dy(j)+(k(i, j)+k(i, j-1))/2/dy(j-1));
           A(I*(j-1)+i, I*(j-1)+i+1)=-2/(dx(i)+dx(i-1))*(k(i+1, j)+k(i, j))/2/dx(i);
           A(I*(j-1)+i, I*(j-1)+i-1)=-2/(dx(i)+dx(i-1))*(k(i-1, j)+k(i, j))/2/dx(i-1);
           A(I*(j-1)+i, I*(j+1-1)+i)=-2/(dy(j)+dy(j-1))*(k(i, j+1)+k(i, j))/2/dy(j);
           A(I*(j-1)+i, I*(j-1-1)+i)=-2/(dy(j)+dy(j-1))*(k(i, j-1)+k(i, j))/2/dy(j-1);
           B(I*(j-1)+i)=f(i, j) + r(i, j)*c(i,j)/dt*Tp(i,j);
           % первый индекс — номер уравнения, не меняется
        
        end
    end
    
    TL=A^(-1)*B; % решаем СЛАУ
    for i=1:I
        for j=1:J
            T(i, j)=TL(I*(j-1)+i);
        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

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