implicit numerical method: 1-D unsteady state heat transfer using finite difference method
显示 更早的评论
% 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 个评论
Joydeep Bagchi
2024-2-24
编辑：Joydeep Bagchi
2024-2-24
Further, you have a two-layer problem. Where do you discretize the transmission condition
T1 = T2
lambda1*dT1/dx = lambda2*dT2/dx
at the interface ?
I can only advise you to use "pdepe" (at least at the beginning) to have a valid reference solution.
The point at the interface should be part of the mesh for this solver.
I'm confused about the thermal diffusivity settings in your code. First you write
if y(i) <= layer2_thickness
k = layer2_k;
...
else
k = layer1_k;
...
end
but when you incorporate the boundary conditions, you work with layer1_k at y(1) and layer2_k at y(end).
This doesn't make sense in my opinion.
Joydeep Bagchi
2024-2-24
编辑：Joydeep Bagchi
2024-2-24
if you consider a paper and across its thickness of length say 1mm so 0.5 mm is of maybe paper 1 type and 0.5 mm may be of paper 2 type so the properties in those 2 layers will be different
I know. That's why you have to establish that temperature and heat flux at the interface point are continuous:
T1 = T2
lambda1*dT1/dx = lambda2*dT2/dx
where T1 and T2 are the limiting temperatures from the left and from the right, respectively.
Writing T1 = T2 and lambda1*dT1/dx = lambda2*dT2/dx does not mean that you impose temperature at the interface, but that you need these two conditions to compute the "correct" temperature that satisfies these two conditions.
Joydeep Bagchi
2024-2-25
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,:)
Joydeep Bagchi
2024-2-25
采纳的回答
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 个评论
Joydeep Bagchi
2024-2-25
I will apply the full problem statement to this now as it was only the part of the code, let's see if it can handle the big change in thermal conductivity value accurately or not
Note that pdepe should only be used for parabolic-elliptic pdes. Thus if your problem statement has equations without a second-order term, you should be careful.
Here, the condition at the interface is not defined as you have mentioned earlier, right?
According to the documentation of "pdepe", continuity of temperature and heat flux are established by the "pdepe" discretization if the interface point is part of the xmesh.
Joydeep Bagchi
2024-2-25
Joydeep Bagchi
2024-2-25
Have I defined the source term wrongly or defined the boundary condition for phase=3 wrongly because if I run the program then it is giving error that spatial discretization is not possible.
At the moment, you don't define a boundary condition at yr for phase 3:
pr = 0;
%pr = h_air*(ur-T_ambient)
qr = 0;
If radiation is present, isn't it a comination of convective and radiative heat transfer at the surface ?
lambda*dT/dy = h*(T-Tambient) + e*(T^4-Tambient^4)
or something of that kind ? Are you sure a source term should be set within the layers ?
Joydeep Bagchi
2024-2-25
Most probably you read about a model that is not discretized in spatial direction, but determines a single temperature only. Here, the energy input is directly set in the equation. You have a partial differential equation where usually (apart from chemical reaction heat) energy input is set across the boundaries.
Anyhow: Mathematically, a partial differential equation of second order needs two boundary conditions. And you have set only one.
Joydeep Bagchi
2024-2-25
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.
更多回答（1 个）
William Rose
2024-2-24
编辑：William Rose
2024-2-24
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
类别
在 帮助中心 和 File Exchange 中查找有关 General Physics 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
