Central and Forward difference schemes heat conduction for two-layer materials
8 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Sanley Guerrier
le 16 Mar 2024
Commenté : Sanley Guerrier
le 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')
Réponse acceptée
Torsten
le 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
8 commentaires
Torsten
le 16 Mar 2024
Modifié(e) : Torsten
le 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.
Plus de réponses (0)
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!