Example 4.3

%%HTProjEX4point3

clear

clc

% given values

k = 25; % W/m-K

dx = 0.001; % m

h_outer = 1000; % W/m^2-K

h_inner = 200; % W/m^2-K

T_outer = 1700; % K

T_inner = 400; % K

% blank matrix A and constant vector C

A = zeros(21,21);

C = zeros(21,1);

% row then column

% node 1

A(1,1) = -(2+(h_outer*dx/k));

A(1,2) = 1;

A(1,7) = 1;

C(1) = -(h_outer*dx/k)*T_outer;

% node 2

A(2,1) = 1;

A(2,2) = -2*((h_outer*dx/k) + 2);

A(2,3) = 1;

A(2,8) = 2;

C(2) = -2*(h_outer*dx/k)*T_outer;

% node 3

A(3,2) = 1;

A(3,3) = -2*((h_outer*dx/k) + 2);

A(3,4) = 1;

A(3,9) = 2;

C(3) = -2*(h_outer*dx/k)*T_outer;

% node 4

A(4,3) = 1;

A(4,4) = -2*((h_outer*dx/k) + 2);

A(4,5) = 1;

A(4,10) = 2;

C(4) = -2*(h_outer*dx/k)*T_outer;

% node 5

A(5,4) = 1;

A(5,5) = -2*((h_outer*dx/k) + 2);

A(5,6) = 1;

A(5,11) = 2;

C(5) = -2*(h_outer*dx/k)*T_outer;

% node 6

A(6,5) = 1;

A(6,6) = -(2+(h_outer*dx/k));

A(6,12) = 1;

C(6) = -2*(h_outer*dx/k)*T_outer;

% node 7

A(7,1) = 1;

A(7,7) = -4;

A(7,8) = 2;

A(7,13) = 1;

C(7) = 0;

% node 8

A(8,2) = 1;

A(8,7) = 1;

A(8,8) = -4;

A(8,9) = 1;

A(8,14) = 1;

C(8) = 0;

% node 9

A(9,3) = 1;

A(9,8) = 1;

A(9,9) = -4;

A(9,10) = 1;

A(9,15) = 1;

C(9) = 0;

% node 10

A(10,4) = 1;

A(10,9) = 1;

A(10,10) = -4;

A(10,11) = 1;

A(10,16) = 1;

C(10) = 0;

% Node 11

A(11,5) = 1;

A(11,10) = 1;

A(11,11) = -4;

A(11,12) = 1;

A(11,17) = 1;

% node 12

A(12,6) = 1;

A(12,11) = 2;

A(12,12) = -4;

A(12,18) = 1;

C(12) = 0;

% node 13

A(13,7) = 1;

A(13,13) = -4;

A(13,14) = 2;

A(13,19) = 1;

C(13) = 0;

% node 14

A(14,8) = 1;

A(14,13) = 1;

A(14,14) = -4;

A(14,15) = 1;

A(14,20) = 1;

C(14) = 0;

% node 15

A(15,9) = 2;

A(15,14) = 2;

A(15,15) = -2*(3+(h_inner*dx/k));

A(15,16) = 1;

A(15,21) = 1;

C(15) = -2*(h_inner*dx/k)*T_inner;

% node 16

A(16,10) = 2;

A(16,15) = 1;

A(16,16) = -2*((h_inner*dx/k) + 2);

A(16,17) = 1;

C(16) = -2*(h_inner*dx/k)*T_inner;

% node 17

A(17,11) = 2;

A(17,16) = 1;

A(17,17) = -2*((h_inner*dx/k) + 2);

A(17,18) = 1;

C(17) = -2*(h_inner*dx/k)*T_inner;

% node 18

A(18,12) = 1;

A(18,17) = 1;

A(18,18) = -(2+(h_inner*dx/k));

C(18) = -(h_inner*dx/k)*T_inner;

% node 19

A(19,13) = 1;

A(19,19) = -2;

A(19,20) = 1;

% node 20

A(20,14) = 1;

A(20,19) = 1;

A(20,20) = -4;

A(20,21) = 2;

% node 21

A(21,15) = 1;

A(21,20) = 1;

A(21,21) = -(2+(h_inner*dx/k));

C(21) = -2*(h_inner*dx/k)*T_inner;

% solve for nodal temperatures

T = A \ C;

% results

disp('Nodal temperatures for Example 4.3 (K):')

for n = 1:21

fprintf('T%-2.0f = %8.2f K\n', n, T(n))

end

Example 5.11

%% Group Project

% Example Problem 5.11 Explicit and Implicit Finite Difference Techniques

clear; clc; close all;

% Given

k = 30; % W/m*K

alpha = 5e-6; % m²/s

rho_c = k / alpha; % J/m³*K

h = 1100; % W/m²*K

T_inf = 250; % Celsius

q1 = 1e7; % W/m³

q2 = 2e7; % W/m³

L = 0.01; % m

dx = 0.002; % m (L/5)

N = 6; % Number of nodes (0 to 5)

t_end = 1.5; % Simulation end time [s]

% Biot and Stability

Bi = h * dx / k;

% Stability criterion for explicit

Fo_max = 0.5 / (1 + Bi);

dt_max = Fo_max * dx^2 / alpha;

% Initial Temperature Distribution

T_s0 = T_inf + q1 * L / h;

x_nodes = (0:5) * dx; % x positions: 0, 2, 4, 6, 8, 10 mm

T_init = (q1 / (2*k)) * (L^2 - x_nodes.^2) + T_s0;

%% Explicit Method

dt_exp = 0.3; % Time step [s] (within stability limit)

Fo_exp = alpha * dt_exp / dx^2; % Fourier number

nsteps = round(t_end / dt_exp);

T_exp = T_init;

T_exp_history = zeros(nsteps+1, N);

T_exp_history(1,:) = T_exp;

% Source term coefficient

S = q2 * dx^2 / k;

for p = 1:nsteps

T_new = T_exp;

% Node 0 (midplane, symmetry: T_{-1} = T_1)

T_new(1) = Fo_exp (2T_exp(2) + S) + (1 - 2*Fo_exp) * T_exp(1);

% Interior nodes 1–4

for m = 2:5

T_new(m) = Fo_exp (T_exp(m-1) + T_exp(m+1) + S) + (1 - 2Fo_exp) * T_exp(m);

end

% Node 5 (surface, convection)

T_new(6) = 2*Fo_exp (T_exp(5) + BiT_inf + S/2) + (1 - 2*Fo_exp - 2*Fo_exp*Bi) * T_exp(6);

T_exp = T_new;

T_exp_history(p+1,:) = T_exp;

end

%% Implicit Method

dt_imp = 0.3;

Fo_imp = alpha * dt_imp / dx^2;

nsteps_imp = round(t_end / dt_imp);

T_imp = T_init'; % Column vector

T_imp_history = zeros(nsteps_imp+1, N);

T_imp_history(1,:) = T_imp';

A = zeros(N, N);

% Node 0 (midplane symmetry)

A(1,1) = 1 + 2*Fo_imp;

A(1,2) = -2*Fo_imp;

% Interior nodes 1-4

for m = 2:5

A(m, m-1) = -Fo_imp;

A(m, m) = 1 + 2*Fo_imp;

A(m, m+1) = -Fo_imp;

end

% Node 5 (surface)

A(6,5) = -2*Fo_imp;

A(6,6) = 1 + 2*Fo_imp + 2*Fo_imp*Bi;

for p = 1:nsteps_imp

b = T_imp;

% Add source terms

b(1) = b(1) + Fo_imp q2 dx^2 / k;

for m = 2:5

b(m) = b(m) + Fo_imp q2 dx^2 / k;

end

b(6) = b(6) + Fo_imp q2 dx^2 / k + 2*Fo_imp Bi T_inf; % surface

T_imp = A \ b;

T_imp_history(p+1,:) = T_imp';

end

%% Spatial temperature profile at t = 1.5 s

figure('Name','Spatial Profile at t=1.5s','Position',[100 100 700 450]);

x_mm = x_nodes * 1000; % convert to mm

plot(x_mm, T_exp_history(end,:), 'b-o', 'LineWidth', 2, 'MarkerSize', 8, 'DisplayName', 'Explicit');

hold on;

plot(x_mm, T_imp_history(end,:), 'r-s', 'LineWidth', 2, 'MarkerSize', 8, 'DisplayName', 'Implicit');

xlabel('x (mm)', 'FontSize', 12);

ylabel('Temperature (Celsius)', 'FontSize', 12);

title('Temperature Distribution at t = 1.5 s', 'FontSize', 13);

legend('Location','best', 'FontSize', 11);

grid on;

%% Long-run transient (t = 0 to 400 s)

dt_long = 0.3;

Fo_long = alpha * dt_long / dx^2;

nsteps_L = round(400 / dt_long);

t_vec = (0:nsteps_L) * dt_long;

T_long = T_init;

T0_hist = zeros(1, nsteps_L+1); T0_hist(1) = T_long(1);

T5_hist = zeros(1, nsteps_L+1); T5_hist(1) = T_long(6);

for p = 1:nsteps_L

T_new = T_long;

T_new(1) = Fo_long*(2*T_long(2) + S) + (1-2*Fo_long)*T_long(1);

for m = 2:5

T_new(m) = Fo_long*(T_long(m-1)+T_long(m+1)+S) + (1-2*Fo_long)*T_long(m);

end

T_new(6) = 2*Fo_long*(T_long(5)+Bi*T_inf+S/2) + (1-2*Fo_long-2*Fo_long*Bi)*T_long(6);

T_long = T_new;

T0_hist(p+1) = T_long(1);

T5_hist(p+1) = T_long(6);

end

% Analytical new steady-state

T_s_new = T_inf + q2*L/h;

T0_ss_new = q2/(2*k)*L^2 + T_s_new;

figure('Name','Transient History 0-400s','Position',[820 100 700 450]);

plot(t_vec, T0_hist, 'r-', 'LineWidth', 2, 'DisplayName', 'T_0 (midplane)');

hold on;

plot(t_vec, T5_hist, 'b-', 'LineWidth', 2, 'DisplayName', 'T_5 (surface)');

yline(T0_ss_new, 'r--', 'LineWidth', 1.2, 'DisplayName', sprintf('T_0 SS = %.1fCelsius', T0_ss_new));

yline(T_s_new, 'b--', 'LineWidth', 1.2, 'DisplayName', sprintf('T_5 SS = %.1fCelsius', T_s_new));

xlabel('t (s)', 'FontSize', 12);

ylabel('Temperature (°C)', 'FontSize', 12);

title('Transient Temperature History (Explicit, \Deltat = 0.3 s)', 'FontSize', 13);

legend('Location','best', 'FontSize', 11);

xlim([0 400]); grid on;