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

78 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(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

### Community Treasure Hunt

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

Start Hunting!