Central and Forward difference schemes heat conduction for two-layer materials
7 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Sanley Guerrier
am 16 Mär. 2024
Kommentiert: Sanley Guerrier
am 16 Mär. 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')
Akzeptierte Antwort
Torsten
am 16 Mär. 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
8 Kommentare
Torsten
am 16 Mär. 2024
Bearbeitet: Torsten
am 16 Mär. 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.
Weitere Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!