implicit numerical method: 1-D unsteady state heat transfer using finite difference method
64 views (last 30 days)
Show older comments
% Constants
Ly = 0.0015; % Total thickness of the system (meters)
T = 40; % Total simulation time (seconds)
Ny = 51; % Number of spatial grid points
Nt = 50000; % Number of time steps
velocity_x = 0.2; % line speed in m/s
% Thermal properties of different layers
layer1_thickness = 0.001; % Thickness of layer 1 (meters)
layer1_k = 0.15; % Thermal conductivity of layer 1 (W/m-K)
layer2_thickness = 0.5e-3; % Thickness of layer 2 (meters)
layer2_k = 0.16; % Thermal conductivity of layer 2 (W/m-K)
rho_1 = 1250; %density of layer 1 (kg/m³)
rho_2 = 1520; %density of layer 2 (kg/m³)
cp_1 = 1350; %specific heat capacity of layer 1 (J/kg-K)
cp_2 = 1460; %specific heat capacity of layer 2 (J/kg-K)
T_roller = 150; %Temperature of hot roller (°C)
h_roller = 500; %Convective heat transfer coefficient (W/m2-K)
T_cooling_roller = 22; %Temperature of cooling temperature(°C)
% layer3_thickness = 0.002; % Thickness of layer 3 (meters)
% layer3_k = 0.3; % Thermal conductivity of layer 3 (W/m-K)
% Ambient temperature and convection properties
T_ambient = 25; % Ambient temperature (°C)
h_air = 12; % Convective heat transfer coefficient of air (W/m²-K)
radiation_length = 1.280; %radiation length in m
radiation_width = 2.300; % radiation width in m
radiation_power = 138e3; %W
radiative_flux = 50e3; %W/m2
absorption = 45; % in %
performance = 95; % in %
emissivity = 0.91;
%if radiative flux is given
%Net_radiative_intensity = radiative_flux*(absorption/100)*(performance/100);
%if radiation power is given
Net_radiative_intensity = (radiation_power/(radiation_length*radiation_width))*(absorption/100)*(performance/100)*emissivity;
% Discretization
dy = Ly / (Ny - 1);
dt = T / Nt;
y = linspace(0, Ly, Ny);
t = linspace(0, T, Nt);
% Initialize temperature matrix
u = zeros(Ny, Nt);
initial_temperature1 = 25;
initial_temperature2 = 25;
%Assigning initial values to the temperature matrix
for i=1:Ny
if y(i) <= layer2_thickness
u(i,1) = initial_temperature1;
else
u(i,1) = initial_temperature2;
end
end
distance_travelled = zeros(1,Nt); %initial distance travelled
dia_roller = 0.8; %diameter of hot roller (m)
dia_cooling_roller = 0.57; %diameter of cooling roller(m)
contact_angle = 180; %contact angle/wrap angle for hot_roller in °
contact_angle_cool1 = 120; %contact angle/wrap angle for cooling_roller1 in °
contact_angle_cool2 = 220; %contact angle/wrap angle for cooling_roller2 in °
contact_angle_cool3 = 190; %contact angle/wrap angle for cooling_roller3 in °
y_1 = pi*dia_roller*contact_angle/360; %in contact with hot roller
y_2 = 0.3; %after heating roller
p0 = y_1;
p1 = p0+y_2;
% Time-stepping loop (explicit method)
for n = 1:Nt - 1
% Update position in x-direction
distance_travelled(n+1) = distance_travelled(n) + velocity_x * dt;
% Check if the total distance is within the contact length
if distance_travelled(n) <= p0
for i = 2:Ny - 1
% Determine the thermal conductivity based on the layer
if y(i) <= layer2_thickness
k = layer2_k;
rho = rho_2;
cp = cp_2;
alpha = k/(rho*cp);
r = alpha*dt/(dy^2);
else
k = layer1_k;
rho = rho_1;
cp = cp_1;
alpha = k/(rho*cp);
r = alpha*dt/(dy^2);
end
u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
end
% Apply convective boundary condition
u(Ny, n + 1) = ((h_roller * T_roller * dy) + (layer2_k * u(Ny - 1, n + 1))) / (layer2_k + h_roller * dy);
u(1, n + 1) = ((h_air * T_ambient * dy) + (layer1_k * u(2, n + 1))) / (layer1_k + h_air * dy);
else
for i = 2:Ny - 1
% Determine the thermal conductivity based on the layer
if y(i) <= layer2_thickness
k = layer2_k;
rho = rho_2;
cp = cp_2;
alpha = k/(rho*cp);
r = alpha*dt/(dy^2);
else
k = layer1_k;
rho = rho_1;
cp = cp_1;
alpha = k/(rho*cp);
r = alpha*dt/(dy^2);
end
u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
end
% Apply convective boundary condition
u(1, n + 1) = ((h_air * T_ambient * dy) + (layer2_k * u(2, n + 1))) / (layer2_k + h_air * dy);
u(Ny, n + 1) = ((h_air * T_ambient * dy) + (layer1_k * u(Ny - 1, n + 1))) / (layer1_k + h_air * dy);
end
% Break the loop when the total distance is reached
if distance_travelled(n + 1) >= p9
break;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time-stepping loop (crank nicholson method)
for n = 1:Nt - 1
% Update position in x-direction
distance_travelled(n+1) = distance_travelled(n) + velocity_x * dt;
% Check if the total distance is within the contact length
if distance_travelled(n) <= p0
for i = 2:Ny - 1
% Determine the thermal conductivity based on the layer
if y(i) <= layer2_thickness
k = layer2_k;
rho = rho_2;
cp = cp_2;
alpha = k/(rho*cp);
r = alpha*dt/(2*dy^2);
else
k = layer1_k;
rho = rho_1;
cp = cp_1;
alpha = k/(rho*cp);
r = alpha*dt/(2*dy^2);
end
u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n) + u(i + 1, n + 1) - 2 * u(i, n + 1) + u(i - 1, n + 1));
end
% Apply convective boundary condition
u(Ny, n + 1) = ((h_roller * T_roller * dy) + (layer2_k * u(Ny - 1, n + 1))) / (layer2_k + h_roller * dy);
u(1, n + 1) = ((h_air * T_ambient * dy) + (layer1_k * u(2, n + 1))) / (layer1_k + h_air * dy);
else
for i = 2:Ny - 1
% Determine the thermal conductivity based on the layer
if y(i) <= layer2_thickness
k = layer2_k;
rho = rho_2;
cp = cp_2;
alpha = k/(rho*cp);
r = alpha*dt/(dy^2);
else
k = layer1_k;
rho = rho_1;
cp = cp_1;
alpha = k/(rho*cp);
r = alpha*dt/(dy^2);
end
u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n) + u(i + 1, n + 1) - 2 * u(i, n + 1) + u(i - 1, n + 1));
end
% Apply convective boundary condition
u(1, n + 1) = ((h_air * T_ambient * dy) + (layer2_k * u(2, n + 1))) / (layer2_k + h_air * dy);
u(Ny, n + 1) = ((h_air * T_ambient * dy) + (layer1_k * u(Ny - 1, n + 1))) / (layer1_k + h_air * dy);
end
% Break the loop when the total distance is reached
if distance_travelled(n + 1) >= p1
break;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time-stepping loop (backward euler method)
for n = 1:Nt - 1
% Update position in x-direction
distance_travelled(n+1) = distance_travelled(n) + velocity_x * dt;
% Check if the total distance is within the contact length
if distance_travelled(n) <= p0
for i = 2:Ny - 1
% Determine the thermal conductivity based on the layer
if y(i) <= layer2_thickness
k = layer2_k;
rho = rho_2;
cp = cp_2;
%alpha = k/(rho*cp);
%r = alpha*dt/(2*dy^2);
else
k = layer1_k;
rho = rho_1;
cp = cp_1;
%alpha = k/(rho*cp);
%r = alpha*dt/(dy^2);
end
% Calculate coefficient matrix for implicit equation
alpha = k / (rho * cp);
r = alpha * dt / dy^2;
A = diag(1 + 2 * r * ones(Ny-2, 1)) + diag(-r * ones(Ny-3, 1), 1) + diag(-r * ones(Ny-3, 1), -1);
%u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n) + u(i + 1, n + 1) - 2 * u(i, n + 1) + u(i - 1, n + 1));
%end
% Apply convective boundary condition
u(Ny, n + 1) = ((h_roller * T_roller * dy) + (layer2_k * u(Ny - 1, n + 1))) / (layer2_k + h_roller * dy);
u(1, n + 1) = ((h_air * T_ambient * dy) + (layer1_k * u(2, n + 1))) / (layer1_k + h_air * dy);
% Right-hand side of the implicit equation
b = u(2:end-1, n) + r * (u(3:end, n) - 2 * u(2:end-1, n) + u(1:end-2, n));
% Solve the implicit equation
u(2:end-1, n + 1) = A \ b;
end
else
for i = 2:Ny - 1
% Determine the thermal conductivity based on the layer
if y(i) <= layer2_thickness
k = layer2_k;
rho = rho_2;
cp = cp_2;
%alpha = k/(rho*cp);
%r = alpha*dt/(2*dy^2);
else
k = layer1_k;
rho = rho_1;
cp = cp_1;
%alpha = k/(rho*cp);
%r = alpha*dt/(dy^2);
end
% Calculate coefficient matrix for implicit equation
alpha = k / (rho * cp);
r = alpha * dt / dy^2;
A = diag(1 + 2 * r * ones(Ny-2, 1)) + diag(-r * ones(Ny-3, 1), 1) + diag(-r * ones(Ny-3, 1), -1);
%u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n) + u(i + 1, n + 1) - 2 * u(i, n + 1) + u(i - 1, n + 1));
%end
% Apply convective boundary condition
u(Ny, n + 1) = ((h_roller * T_roller * dy) + (layer2_k * u(Ny - 1, n + 1))) / (layer2_k + h_roller * dy);
u(1, n + 1) = ((h_air * T_ambient * dy) + (layer1_k * u(2, n + 1))) / (layer1_k + h_air * dy);
% Right-hand side of the implicit equation
b = u(2:end-1, n) + r * (u(3:end, n) - 2 * u(2:end-1, n) + u(1:end-2, n));
% Solve the implicit equation
u(2:end-1, n + 1) = A \ b;
end
end
% Break the loop when the total distance is reached
if distance_travelled(n + 1) >= p1
break;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the temperature distribution over time
% Find the last valid time step where u(:, Nt) is non-zero
last_valid_time_step = find(u(1, :) ~= 0, 1, 'last');
figure();
plot(t(1:last_valid_time_step), u(1, 1:last_valid_time_step), ...
'b', 'DisplayName', 'Unterseite');
hold on
plot(t(1:last_valid_time_step), u(50, 1:last_valid_time_step), ...
'r', 'DisplayName', 'Oberseite');
hold on
plot(t(1:last_valid_time_step), u(25, 1:last_valid_time_step), ...
'g', 'DisplayName', 'Interface');
xlabel('Time (seconds)');
ylabel('Temperature (°C)');
yticks(20:10:200);
xticks(0 :1 :t(1,last_valid_time_step));
legend('Location', 'Best');
title('Temperature Distribution in composite');
grid on;
axis([0 t(1,last_valid_time_step) 20 200]);
hold off
I want to apply implicit method to the 1-D unsteady state heat transfer problem to diminsh the effect of large thermal conductivity or very small densities or specific heat capacities. I was writing the code according to the classic example given in book and youtube videos but I am unable to match upto the solution provided by Explicit numerical method. Also there is no example of using crank nicholson or backward euler using convective boundary conditions or conductive boundary conditions.
Could you explain why is this so?
8 Comments
Torsten
on 25 Feb 2024
I usually mesh from a to the interface point and from the interface point to b and concatenate the grids (like in the pdepe code below).
I don't understand what you mean by
u(Ny/2,:)=u(Ny+1/2,:)
Accepted Answer
Torsten
on 24 Feb 2024
Edited: Torsten
on 24 Feb 2024
pde()
function pde()
% Constants
Ly = 0.0015; % Total thickness of the system (meters)
T = 40; % Total simulation time (seconds)
Ny = 51; % Number of spatial grid points
Nt = 50000; % Number of time steps
velocity_x = 0.2; % line speed in m/s
% Thermal properties of different layers
layer1_thickness = 0.001; % Thickness of layer 1 (meters)
layer1_k = 0.15; % Thermal conductivity of layer 1 (W/m-K)
layer2_thickness = 0.5e-3; % Thickness of layer 2 (meters)
layer2_k = 0.16; % Thermal conductivity of layer 2 (W/m-K)
rho_1 = 1250; %density of layer 1 (kg/m³)
rho_2 = 1520; %density of layer 2 (kg/m³)
cp_1 = 1350; %specific heat capacity of layer 1 (J/kg-K)
cp_2 = 1460; %specific heat capacity of layer 2 (J/kg-K)
T_roller = 150; %Temperature of hot roller (°C)
h_roller = 500; %Convective heat transfer coefficient (W/m2-K)
T_cooling_roller = 22; %Temperature of cooling temperature(°C)
% layer3_thickness = 0.002; % Thickness of layer 3 (meters)
% layer3_k = 0.3; % Thermal conductivity of layer 3 (W/m-K)
% Ambient temperature and convection properties
T_ambient = 25; % Ambient temperature (°C)
h_air = 12; % Convective heat transfer coefficient of air (W/m²-K)
radiation_length = 1.280; %radiation length in m
radiation_width = 2.300; % radiation width in m
radiation_power = 138e3; %W
radiative_flux = 50e3; %W/m2
absorption = 45; % in %
performance = 95; % in %
emissivity = 0.91;
%if radiative flux is given
%Net_radiative_intensity = radiative_flux*(absorption/100)*(performance/100);
%if radiation power is given
Net_radiative_intensity = (radiation_power/(radiation_length*radiation_width))*(absorption/100)*(performance/100)*emissivity;
%Discretization in space
ny1 = ceil(layer2_thickness/Ly*Ny);
ny2 = Ny-ny1;
y1 = linspace(0,layer2_thickness,ny1);
y2 = linspace(layer2_thickness,Ly,ny2);
y = [y1,y2(2:end)];
% Initialize temperature matrix
initial_temperature1 = 25;
initial_temperature2 = 25;
dia_roller = 0.8; %diameter of hot roller (m)
dia_cooling_roller = 0.57; %diameter of cooling roller(m)
contact_angle = 180; %contact angle/wrap angle for hot_roller in °
contact_angle_cool1 = 120; %contact angle/wrap angle for cooling_roller1 in °
contact_angle_cool2 = 220; %contact angle/wrap angle for cooling_roller2 in °
contact_angle_cool3 = 190; %contact angle/wrap angle for cooling_roller3 in °
y_1 = pi*dia_roller*contact_angle/360; %in contact with hot roller
y_2 = 0.3; %after heating roller
p0 = y_1;
p1 = p0+y_2;
t1 = p0/velocity_x;
t2 = p1/velocity_x-t1;
nt1 = ceil(t1/(t1+t2)*Nt);
nt2 = Nt-nt1;
t11 = linspace(0,t1,nt1);
t12 = linspace(t1,t1+t2,nt2);
%pdepe settings
m = 0;
phase = 1;
sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t11);
u1 = sol(:,:,1);
phase = 2;
sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t12);
u2 = sol(:,:,1);
plot([t11,t12],[[u1(:,1);u2(:,1)],[u1(:,25);u2(:,25)],[u1(:,50);u2(:,50)]])
grid on
function [c f s] = pdefun(y,t,u,dudy)
if y <= layer2_thickness
c = rho_2*cp_2;
f = layer2_k*dudy;
s = 0;
else
c = rho_1*cp_1;
f = layer1_k*dudy;
s = 0;
end
end
function u = icfun(yq)
if phase==1
if yq <= layer2_thickness
u = initial_temperature1;
else
u = initial_temperature2;
end
else
u = interp1(y,u1(end,:),yq);
end
end
function [pl ql pr qr] = bcfun(yl,ul,yr,ur,t)
if phase==1
pl = -h_air*(ul-T_ambient);
ql = 1;
pr = h_roller*(ur-T_roller);
qr = 1;
else
pl = -h_air*(ul-T_ambient);
ql = 1;
pr = h_air*(ur-T_ambient);
qr = 1;
end
end
end
9 Comments
Torsten
on 26 Feb 2024
Edited: Torsten
on 27 Feb 2024
As I said: I'd compute the radiation heat flux at the upper surface, add it to the convective heat flux and remove the source term. But I'm not really a process engineer.
See the chapter "Heat Transfer between surfaces" under
how the radiative heat flux between surfaces is usually modelled.
More Answers (1)
William Rose
on 24 Feb 2024
Edited: William Rose
on 24 Feb 2024
I think you need to divide the long term (u(i+1,n)-2*u(i,n)+...) by 2, on the right side of your Crank-Nicolson update equation. But still, I think it is wrong, because Crank-Nicolson is an implicit method which one solves with a triadiagonal matrix. It is just a different tridiagonal matrix than the one in the backward Euler implicit method. I have not evvaluatred the handling of the boundary conditions.
This is a complicated model. I would test the 3 algorithms on a simpler model first, to see if they give the same results - or close-enough-to-the-same.
How do you know the algoritms give different results? The Crank-Nicolson results appear to overwrite the explicit method results. Then the backward Euler appears to overwrite the Crank-Nicolson. Am I missing something?
You can compute r1 and r2 (for layers 1 and 2), and matrices A1 (using r1) and A2 (using r2) at the start of your code, before all the loops.
r1 = layer1_k*dt/(rho_1*cp_1*2*dy^2); % don't need alpha
r2 = layer2_k*dt/(rho_2*cp_2*2*dy^2);
% equation for A is different for backward Euler and C-N
A1be=[triagonal matix using r1];
A2be=[triagonal matix using r2];
A1cn=[triagonal matix using r1];
A2cn=[triagonal matix using r2];
Then, when looping over over position, you do
for i = 2:Ny - 1
% Determine the thermal conductivity based on the layer
if y(i)<=layer2_thickness
r=r2; A=A2be; % or A2cn
else
r=r1; A=A1be; % or A1cn
end
u(i,n+1) = u(i,n) + r*(stuff);
end
and so on. This will make your code shorter and simpler.
Explicit method update equation (with some spaces removed)
u(i,n+1) = u(i,n) + r*(u(i+1,n)-2*u(i,n)+u(i-1,n));
Your Crank-Nicholson update equation
u(i,n+1) = u(i,n) + r*(u(i+1,n)-2*u(i,n)+u(i-1,n)+u(i+1,n+1)-2*u(i,n+1)+u(i-1,n+1));
Yor Implicit method (backward Euler) update equation
% Coefficient matrix for implicit equation
A = diag(1+2*r*ones(Ny-2,1)) + diag(-r*ones(Ny-3,1),1) + diag(-r*ones(Ny-3,1),-1);
% Right-hand side of the implicit equation
b = u(2:end-1,n) + r * (u(3:end,n)-2*u(2:end-1,n)+u(1:end-2,n));
u(2:end-1,n+1) = A\b; % solve
See Also
Categories
Find more on Heat Transfer in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!