Enthalpy

clc
clear all
global t Tbf Cth Cf L Wtot Ww V tend dt lamda rho Cp T0 a TU TD TR TL H0
t=0;
Tbf=-0.47+273.15;
Cth=2845;
Cf=1899;
L=45000;
Wtot=0.25;
Ww=0.12;
V=Ww/Wtot;
tend=10000;
dt=0.1;
lamda=57;
rho=7900;
Cp=450;
T0=270.15;
Lx=20;
Nx=20;
Ly=20;
Ny=20;
dx=Lx/(Nx-1);
dy=Ly/(Ny-1);
a=lamda/(rho*Cp);
Kx=(dx)^2/(a*dt);
Ky=(dy)^2/(a*dt);
TU=-12+273.15; %температура сверху
TD=-3.5+273.15; %температура снизу
TL=-2.5+273.15; %температура слева
TR=60+273.15; %температура справа
T=zeros(Nx,Ny);
H=zeros(Nx,Ny);
H0=Cf*(T0-Tbf);
Hnew=H0;

% Задаем начальные условия для энтальпии

for i=1:1:Nx
    for j=1:1:Ny
        H(i,j)=H0;
    end
end

% Задаем начальные условия для температуры

for i=1:1:Nx
    for j=1:1:Ny
        T(i,j)=T0;
    end
end

    % Апроксимация граничных условий

    for j=1;
        for i=1;
            Hnew(i,j)=H(i,j)+Kx*(T(i+1,j)-2*T(i,j)+TL)+Ky*(T(i,j+1)-2*T(i,j)+TD);
        end
    end

    for j=Ny;
        for i=1;
            Hnew(i,j)=H(i,j)+Kx*(T(i+1,j)-2*T(i,j)+TL)+Ky*(TU-2*T(i,j)+T(i,j-1));
        end
    end

    for j=1;
        for i=Nx;
            Hnew(i,j)=H(i,j)+Kx*(TR-2*T(i,j)+T(i-1,j))+Ky*(T(i,j+1)-2*T(i,j)+TD);
        end
    end

    for j=Ny;
        for i=Nx;
            Hnew(i,j)=H(i,j)+Kx*(TR-2*T(i,j)+T(i-1,j))+Ky*(TU-2*T(i,j)+T(i,j-1));
        end
    end

    for i=1;
        for j=2:1:Ny-1
            Hnew(i,j)=H(i,j)+Kx*(T(i+1,j)-2*T(i,j)+TL)+Ky*(T(i,j+1)-2*T(i,j)+T(i,j-1));
        end
    end

    for i=Nx;
        for j=2:1:Ny-1
            Hnew(i,j)=H(i,j)+Kx*(TR-2*T(i,j)+T(i-1,j))+Ky*(T(i,j+1)-2*T(i,j)+T(i,j-1));
        end
    end

    for j=1;
        for i=2:1:Nx-1
            Hnew(i,j)=H(i,j)+Kx*(T(i+1,j)-2*T(i,j)+T(i-1,j))+Ky*(T(i,j+1)-2*T(i,j)+TD);
        end
    end

    for j=Ny;
        for i=2:1:Nx-1
            Hnew(i,j)=H(i,j)+Kx*(T(i+1,j)-2*T(i,j)+T(i-1,j))+Ky*(TU-2*T(i,j)+T(i,j-1));
        end
    end

    % Апроксимация уравнения теплопроводности

while t<tend
    for i=2:1:Nx-1
        for j=2:1:Ny-1
            Hnew(i,j)=H(i,j)+Kx*(T(i+1,j)-2*T(i,j)+T(i-1,j))+Ky*(T(i,j+1)-2*T(i,j)+T(i,j-1));
        end
    end
    H=Hnew;
    t=t+dt;
end