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

5 views (last 30 days)
Sanley Guerrier on 16 Mar 2024
Commented: Sanley Guerrier on 16 Mar 2024
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 on 16 Mar 2024
Edited: Torsten on 16 Mar 2024
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 on 16 Mar 2024
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

Torsten on 16 Mar 2024
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
Torsten on 16 Mar 2024
Edited: Torsten on 16 Mar 2024
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 on 16 Mar 2024
it works now with the new script. Thanks much!