implicit numerical method: 1-D unsteady state heat transfer using finite difference method

53 次查看（过去 30 天）
Joydeep Bagchi
Joydeep Bagchi 2024-2-23
编辑： Joydeep Bagchi 2024-2-27
% 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 个评论
显示 6更早的评论
Torsten
Torsten 2024-2-25
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,:)

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2024-2-24
编辑：Torsten 2024-2-24
Ran in:
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 个评论
显示 7更早的评论
Joydeep Bagchi
Joydeep Bagchi 2024-2-25
ok, thankyou for the information. But I only have the data of intensity of IR lamp which is 138 kW, radiation length and width, no other data is available. What else could I set up?
Torsten
Torsten 2024-2-26
编辑：Torsten 2024-2-27
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
https://en.wikipedia.org/wiki/Thermal_radiation
how the radiative heat flux between surfaces is usually modelled.

请先登录，再进行评论。

更多回答（1 个）

William Rose
William Rose 2024-2-24
编辑：William Rose 2024-2-24
@Joydeep Bagchi,
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
  1 个评论
Joydeep Bagchi
Joydeep Bagchi 2024-2-24
Thankyou @William Rose for the reply.
I am running three different matlab files so the constants are same at the beginning, just the time stepping loop is different. all three methods should give about same results and implicit methods should be more robust and unconditionally stable. because with explicit method, i am getting the solution but it heavily depends on parameter 'r' and it depends on density,thermal conductivity and specific heat capacity. So if I change any one of the parameters by a large magnitude then my solution differs a lot which should not be the case.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Heat Transfer 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by