two layer transient conduction heat transfer with pde

I write a code to caculate 1D transient heat conduction with pde. There are 2 layers, the left is concrete contacting with room air and the right is insulation material contacting with atomsphere air. The left and the right boundaries are both natural convection. The initial temperature is 24 degC. the room air temperature is 32 degC and the atomsphere temperature is 35 deg C. The program works but the results looks like wrong. Anybody can help? Thanks in advance! Naihua

%% caculate 1D transient conduction of two layer

clc
clear 
global td tk1 tk2 urMax uAtm

%% parameters

A   = 420;             % heat transfer area, m2
td  = 72 * 3600;       % time duaration (72h), s

Tr0   = 24;            % initial room / wall temperature, deg C 
TrMax = 32;            % room temperature, deg C
Tatm  = 35;            % atomsphere temperature, deg C
urMax = TrMax - Tr0;
uAtm  = Tatm  - Tr0;

tk1    = 0.61;          % concrete wall thickness, m
tk2    = 0.1;           % insulation wall thickness, m

%% define domain
x1 = linspace(0, tk1, 1000);
x2 = linspace (tk1+tk1/1000, tk1+tk2, 1000);
x  = [x1 x2];
t = linspace (0, td, 2000);

%% solve pde
m   = 0;
sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);

%% figure
%----------------------------------------------------------------------------
% plot sol, ie. (x, t)
colormap hot
imagesc(x,t/3600,sol)
colorbar;
ylim([0 72])
yticks([0:6:72])
xlabel('Distance, m', 'interpreter', 'latex')
ylabel('Time, h', 'interpreter', 'latex')
title('Temperatue increase (deg C) for $0\le x\le 0.61m$ and $0 \le t\le 72h$','interpreter', 'latex')

%% define pde problem
%----------------------------------------------------------------------------
% pde function
function [c,f,s] = heatpde(x, t, u, dudx)
    global tk1 k1 k2

    c = 1;
    if x <= tk1

    %% define concrete 
        k1     = 1.54;          % concrete thermal conductivity (W/mK)
        rho1   = 2400;          % concrete density (kg/m^3)
        Cp1    = 840;           % concrete specific heat (J/kgK)
        a1     = k1/(rho1*Cp1); % concrete thermal diffusivity (m^2/s)
        f      = a1*dudx;
    
    else

        %% define insulation
        k2     = 0.024;         % insulation thermal conductivity (W/mK)
        rho2   = 40;            % insulation density (kg/m^3)
        Cp2    = 1400;          % insulation specific heat (J/kgK)
        a2     = k2/(rho2*Cp2); % insulation thermal diffusivity (m^2/s)
        f      = a2*dudx;
    end

    s      = 0;
end
%% initial condtion
function u0 = heatic(x)
    u0 = 0;
end

%% boundary condition

function [pl,ql,pr,qr] = heatbc(xl, ul, xr, ur, t)
    
    global td k1 k2 urMax uAtm hr_int hr_ext
    
    % left boundary - natural convection
    hr_int    = 3.804;   % natural convection heat transfer coefficient, W/(m2-K)
    pl        = hr_int * (urMax - ul);
    ql        = k1;  
    
    % right boundary - natural convection
    hr_ext    = 3.804;    % natural convection heat transfer coefficient, W/(m2-K)
    pr        = hr_ext * (ur - uAtm);
    qr        = k2; 
end

Answer 1

[Torsten]

The boundary condition for pdepe reads

p + q*f = 0.

In your case this leads to

hr_int * (urMax - ul) + k1*k1/(rho1*Cp1)*dT/dx = 0 at x = xl
hr_ext * (ur - uAtm)  + k2*k2/(rho2*Cp2)*dT/dx = 0 at x = xr

I doubt that this is what you want.

[Naihua WANG]

Thanks, dearTorsten, I will try it! Naihua

[Torsten]

I will try it!

?

[Naihua WANG]

Dear Torsten, actually the function pdepe has a lot limits for the above mentioned function with different thermal diffusion coefficient and also for different boundary condition. I finally solved it with the PDE Toolbox.

[Torsten]

Your boundary condition setting is wrong, as I showed you. The problem can easily be solved with pdepe (see below).

%% caculate 1D transient conduction of two layer

clc
clear 
global td tk1 tk2 urMax uAtm

%% parameters

A   = 420;             % heat transfer area, m2
td  = 72 * 3600;       % time duaration (72h), s

Tr0   = 24;            % initial room / wall temperature, deg C 
TrMax = 32;            % room temperature, deg C
Tatm  = 35;            % atomsphere temperature, deg C
urMax = TrMax - Tr0;
uAtm  = Tatm  - Tr0;

tk1    = 0.61;          % concrete wall thickness, m
tk2    = 0.1;           % insulation wall thickness, m

%% define domain
x1 = linspace(0, tk1, 1000);
x2 = linspace (tk1+tk1/1000, tk1+tk2, 1000);
x  = [x1 x2];
t = linspace (0, td, 2000);

%% solve pde
m   = 0;
sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);

%% figure
%----------------------------------------------------------------------------
% plot sol, ie. (x, t)
colormap hot
imagesc(x,t/3600,sol)
colorbar;
ylim([0 72])
yticks([0:6:72])
xlabel('Distance, m', 'interpreter', 'latex')
ylabel('Time, h', 'interpreter', 'latex')
title('Temperatue increase (deg C) for $0\le x\le 0.61m$ and $0 \le t\le 72h$','interpreter', 'latex')

plot(x,[sol(1,:);sol(200,:);sol(400,:);sol(600,:);sol(800,:);sol(1000,:);sol(end,:)])

%% define pde problem
%----------------------------------------------------------------------------
% pde function
function [c,f,s] = heatpde(x, t, u, dudx)
    global tk1 

    if x <= tk1

    %% define concrete 
        k1     = 1.54;          % concrete thermal conductivity (W/mK)
        rho1   = 2400;          % concrete density (kg/m^3)
        Cp1    = 840;           % concrete specific heat (J/kgK)
        a1     = k1/(rho1*Cp1); % concrete thermal diffusivity (m^2/s)
        c      = rho1*Cp1;
        f      = k1*dudx;
        s      = 0;
    
    else

        %% define insulation
        k2     = 0.024;         % insulation thermal conductivity (W/mK)
        rho2   = 40;            % insulation density (kg/m^3)
        Cp2    = 1400;          % insulation specific heat (J/kgK)
        a2     = k2/(rho2*Cp2); % insulation thermal diffusivity (m^2/s)
        c      = rho2*Cp2;
        f      = k2*dudx;
        s      = 0;
    end

end
%% initial condtion
function u0 = heatic(x)
    u0 = 0;
end

%% boundary condition

function [pl,ql,pr,qr] = heatbc(xl, ul, xr, ur, t)
    
    global urMax uAtm 
    
    
    % left boundary - natural convection
    hr_int    = 3.804;   % natural convection heat transfer coefficient, W/(m2-K)
    pl        = hr_int * (urMax - ul);
    ql        = 1;  
    
    % right boundary - natural convection
    hr_ext    = 3.804;    % natural convection heat transfer coefficient, W/(m2-K)
    pr        = hr_ext * (ur - uAtm);
    qr        = 1; 
end