Central and Forward difference schemes heat conduction for two-layer materials

Hi all,

I am trying to implement a central and forward difference schemes heat conduction for two layer materials. My problem is I don't how to insert the continuity equation at the interface layer into the for loop.

I tried this code assuming it is one layer material it worked fine. The only problem is how to set the continuity condition at the interface between the two layers.

I would very much appriciate if someone could help me fix this problem.

Thank you!

clc

clear

T = 4.1563; %final time second

Nt = 999; % number of time step

dt = T/Nt; % time step

D = 94.32*0.0254; %total thickness[m]

Nx = 100; %number of space step

dx = D/Nx; %space step

%material 1 thickness and thermal properties

D1 = 22*0.0254; % thickness [m]

K1 = 0.96; % thermal conductivity [W/(m*K)]

c1 = 840/84600; % specific heat capacity [J/(kg*K)]

rho1 = 2000; % density [kg/m^3]

%material 2 thickness and thermal properties

D2 = 72.32*0.0254; % thickness [m]

K1 = 1.21; %thermal conductivity [W/(m*K)]

c1 = 850; %specific heat capacity [J/(kg*K)]

rho1 = 2360; % density [kg/m^3]

% stability check, beta parameter in FTCS must be less than 1

b = 2.*K*dt/(rho*c*dx^2);

disp(b);

z_mesh = linspace(z_sensor(1), z_sensor(end), 100);

% the initial condition

T0 =interp1(z_sensor, temp(1,:), z_mesh);

for i=1:Nx

u(i,1)= T0(i);

end

% implementation of the explicit FTCS method for two materials

for t =1:Nt

q(t) = q(t);

Ta(t) =Ta(t);

hc(t)=hc(t);

%boundary condition at x = 0

u(1,t+1)= u(1,t) + 2*K1*dt*(u(2,t)-u(1,t))/(rho1*c1*dx^2) + 2*(hc(t)*(Ta(t) -u(1,t))+q(t))*dt/(rho1*c1*dx);

%compute temperature in the interior nodes

for i=2:Nx-1

u(i,t+1) = u(i,t) + (K1*dt)/(rho1*c1*dx^2)*(u(i+1,t) - 2.*u(i,t) + u(i-1,t));

%condition continuity at the interface layer, x = 22

% u(i,t+1) = u(i,t) + 2*K1*dt/((rho1*c1*dx1 +rho2*c2*dx2)*dx1) * (u(i-1,t)-u(i,t)) + 2*K2*dt/((rho1*c1*dx1 +rho2*c2*dx2)*dx2) * (u(i+1,t)-u(i,t));

end

u(Nx,t+1) = u(Nx,t); % at x=94.32 -K*du/dx =0

end

t_mesh = convertTo(time, 'datenum') - convertTo(time(1), 'datenum');

% Make a small suite of plots showing the modelling results.

plots(z_mesh, time, u, z_sensor, temp)

% Make a movie of the evolving temperature profile.

movie(z_mesh, t_mesh, u, z_sensor, temp, 'San.avi')

[Torsten]

Use "pdepe" and use the interface point as a point of the mesh. This will automatically ensure continuity of temperature and heat flux at the interface.

[Sanley Guerrier]

Hi Torsten,

why we connot just use the continuity at the interface with this method to do it?

Is there a pdepe example that uses the interface point as a point of the mesh available? That could be helpful.

Thank you

[Torsten]

xmesh uses xinterface as mesh point:

xstart = 0;
xend = 1;
xinterface = 0.25;
xmesh1 = linspace(xstart,xinterface,25);
xmesh2 = linspace(xinterface,xend,75);
xmesh = [xmesh1,xmesh2(2:end)]

[Sanley Guerrier]

Torsten-

This is how I try to sep up pdepe using xinterface as mesh point you provided. Is this set up right? Thank you

xstart = 0;

xend = 1;

xinterface = 0.25;

xmesh1 = linspace(xstart,xinterface,25);

xmesh2 = linspace(xinterface,xend,75);

xmesh = [xmesh1,xmesh2(2:end)]

% pde function

function [c,f,s] = heatpde(xmesh, t, u, dudx)

c = 1;

if x <= xtinterface

%% define material 1

K1 = 0.96; % thermal conductivity [W/(m*K)]

Cp1 = 840/84600; % specific heat capacity [J/(kg*K)]

rho1 = 2000; % density [kg/m^3]

c1= rho1*Cp1

f = K1*dudx;

s =0;

else

%% define material 2

K2 = 1.21; %thermal conductivity [W/(m*K)]

Cp2 = 850; %specific heat capacity [J/(kg*K)]

rho2 = 2360; % density [kg/m^3]

c2= rho2*Cp2;

f = K2*dudx;

s = 0;

end

end

%% initial condtion

function u0 = heatic(x)

u0 = T0;

end

%% boundary condition

function [pl,ql,pr,qr] = heatbc(xl, ul, xr, ur, t)

% left boundary -K*du/dx = q

pl = q; %heat flux W/m^2

ql = 1;

% right boundary -K*du/dx =0

pr = 1;

qr = 0;

end

Answer 1

[Torsten]

heattransfer()

function heattransfer
xstart = 0;
xend = 1;
xinterface = 0.25;
xmesh1 = linspace(xstart,xinterface,25);
xmesh2 = linspace(xinterface,xend,75);
xmesh = [xmesh1,xmesh2(2:end)];
T0 = 10;
q = 10;
m=0;
tmesh=linspace(0,0.1,10);
sol = pdepe(m,@heatpde,@heatic,@heatbc,xmesh,tmesh);
plot(xmesh,[sol(1,:);sol(2,:);sol(3,:);sol(4,:);sol(5,:);sol(6,:);sol(10,:)])
% pde function
function [c,f,s] = heatpde(x, t, u, dudx)

if x <= xinterface
    %% define material 1 
    K1   = 0.96;                        % thermal conductivity [W/(m*K)]
    Cp1  = 840/84600;                     % specific heat capacity [J/(kg*K)]
    rho1 = 2000;                        % density [kg/m^3]
    c= rho1*Cp1;
    f      = K1*dudx;
    s =0;
    
else
    %% define material 2
    K2 = 1.21;                      %thermal conductivity [W/(m*K)]
    Cp2 = 850;                       %specific heat capacity [J/(kg*K)]
    rho2 = 2360;                    % density [kg/m^3]
    c= rho2*Cp2;
    f      = K2*dudx;
    s      = 0;
end
end

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

%% boundary condition
function [pl,ql,pr,qr] = heatbc(xl, ul, xr, ur, t)
% left boundary K1*du/dx = -q
pl        = q; %heat flux W/m^2
ql        = 1;  

% right boundary -K*du/dx =0
pr        = 0;
qr        = 1; 
end
end

[Sanley Guerrier]

Nice!

Thank you, Torsten

[Sanley Guerrier]

How do I extract the solutions?

I tried

u = sol(:,:,1);

But it didn't give anything

[Sanley Guerrier]

I tried to add another layer following your process the code didn't work and said "The entries of XMESH must be strictly increasing."

xstart = 0;

xend = 1;

xinterface = 0.25;

xinterface1 = 0.25;

xinterface2 = 0.5

xmesh1 = linspace(xstart,xinterface1,25);

xmesh2 = linspace(xinterface1,xinterface2,20);

xmesh3 = linspace(xinterface2, xend, 50)

xmesh = [xmesh1,xmesh2, xmesh3(2:end];

Error using pdepe

The entries of XMESH must be strictly increasing.

Error in heattransfer (line 58)

sol = pdepe(m,@heatpde,@heatic,@heatbc,xmesh,tmesh);

function heattransfer

xstart = 0;

xend = 1;

xinterface = 0.25;

xinterface1 = 0.25;

xinterface2 = 0.5

xmesh1 = linspace(xstart,xinterface1,25);

xmesh2 = linspace(xinterface1,xinterface2,20);

xmesh3 = linspace(xinterface2, xend, 50)

xmesh = [xmesh1,xmesh2, xmesh3(2:end];

T0 = 10;

q = 10;

m=0;

tmesh=linspace(0,0.1,10);

sol = pdepe(m,@heatpde,@heatic,@heatbc,xmesh,tmesh);

plot(xmesh,[sol(1,:);sol(2,:);sol(3,:);sol(4,:);sol(5,:);sol(6,:);sol(10,:)])

% pde function

function [c,f,s] = heatpde(x, t, u, dudx)

if x <= xinterface1

%% define material 1

K1 = 0.96; % thermal conductivity [W/(m*K)]

Cp1 = 840/84600; % specific heat capacity [J/(kg*K)]

rho1 = 2000; % density [kg/m^3]

c= rho1*Cp1;

f = K1*dudx;

s =0;

else if x<= xinterface2

%% define material 2

K2 = 1.21; %thermal conductivity [W/(m*K)]

Cp2 = 850; %specific heat capacity [J/(kg*K)]

rho2 = 2360; % density [kg/m^3]

c= rho2*Cp2;

f = K2*dudx;

s = 0;

else

k3 =1;

Cp3 = 900;

rho3 = 1500;

c =rho3*Cp3;

f = K3*dudx;

s = 0;

end

end

%% initial condtion

function u0 = heatic(x)

u0 = T0;

end

%% boundary condition

function [pl,ql,pr,qr] = heatbc(xl, ul, xr, ur, t)

% left boundary K1*du/dx = -q

pl = q; %heat flux W/m^2

ql = 1;

% right boundary -K*du/dx =0

pr = 0;

qr = 1;

end

end

[Torsten]

heattransfer()

function heattransfer
xstart = 0;
xend = 1;
xinterface1 = 0.25;
xinterface2 = 0.5;
xmesh1 = linspace(xstart,xinterface1,25);
xmesh2 = linspace(xinterface1,xinterface2,20);
xmesh3 = linspace(xinterface2, xend, 50);
xmesh = [xmesh1,xmesh2(2:end), xmesh3(2:end)];
    T0 = 10;
    q = 10;
    m=0;
    tmesh=linspace(0,0.1,10);
    sol = pdepe(m,@heatpde,@heatic,@heatbc,xmesh,tmesh);
    plot(xmesh,[sol(1,:);sol(2,:);sol(3,:);sol(4,:);sol(5,:);sol(6,:);sol(10,:)])
    
    % pde function
    function [c,f,s] = heatpde(x, t, u, dudx)
    if x <= xinterface1
        %% define material 1 
        K1   = 0.96;                        % thermal conductivity [W/(m*K)]
        Cp1  = 840/84600;                     % specific heat capacity [J/(kg*K)]
        rho1 = 2000;                        % density [kg/m^3]
        c= rho1*Cp1;
        f      = K1*dudx;
        s =0;
    elseif x<= xinterface2
            %% define material 2
            K2 = 1.21;                      %thermal conductivity [W/(m*K)]
            Cp2 = 850;                       %specific heat capacity [J/(kg*K)]
            rho2 = 2360;                    % density [kg/m^3]
            c= rho2*Cp2;
            f      = K2*dudx;
            s      = 0;
    else
        K3 =1;
        Cp3 = 900;
        rho3 = 1500;
        c =rho3*Cp3;
        f      = K3*dudx;
        s      = 0;
    end
    end
    %% initial condtion
        function u0 = heatic(x)
            u0 = T0;
        end
        %% boundary condition
        function [pl,ql,pr,qr] = heatbc(xl, ul, xr, ur, t)
            % left boundary K1*du/dx = -q
            pl        = q; %heat flux W/m^2
            ql        = 1;  
            % right boundary -K*du/dx =0
            pr        = 0;
            qr        = 1; 
        end
    end

[Torsten]

How do I extract the solutions?

I tried

u = sol(:,:,1);

But it didn't give anything

sol(i,j) is the temperature in xmesh(j) at time tmesh(i).

[Sanley Guerrier]

I used disp(sol) I could see the solutions but what I want is to save it in my workspace.

[Torsten]

Then define "sol" as output argument from function "heattransfer" (preferably with "xmesh" and "tmesh").

[xmesh,tmesh,sol] = heattransfer()

and do the plotting in the script part.

Maybe defining parameters in the script part and using them is inputs to "heattransfer" is also an option.

[Sanley Guerrier]

it works now with the new script. Thanks much!