NUMERICAL METHODS FOR MECHANICS II

ASSIGNMENT REPORT

Student Name: ABIODUN GRACE AIYEDUN

Matric Number: 240404509

Department: Mechanical Engineering

Level: 300 Level

Course: Numerical Methods for Mechanics II

QUESTIONS COVERED IN THIS REPORT

Question Topic Methods Used

Q1 (a) 2D Steady-State Fin — Case I Finite Difference Method (FDM)

Q1 (b) 1D Steady-State Fin — Case II RK4 + Shooting Method

Q1 (c) 1D Transient Fin — Case III Crank-Nicolson FDM

Q2 Apollo Spacecraft Orbit RK4 + Central FDM

Q3 (i) FSW Steel Properties Newton-Gregory Forward Difference

Q3 (ii) Thermal Conductivity Root Newton-Raphson + Inv. Lagrange

Q3 (iii) FSW Transient 1D Crank-Nicolson + RK4-Shooting

Q3 (iv) FSW Transient 2D ADI Method

Q4 & Q5 Composite Longitudinal Fin CN-FDM, Galerkin FEM, FVM

Q6 Free-Piston Stirling Engine RK4 + FDM + Analytical

Q7 Rolling Sheet Heat Transfer Finite Volume + Galerkin FEM

 QUESTION 1 — Two-Dimensional Heat Transfer in a Fin
with Convective-Radiative Boundaries

A. Requirement Checklist

Requirement Detail

Governing PDE ∂/∂x[k∂T/∂x] + ∂/∂y[k∂T/∂y] − (2h/w)(T−T∞) − (2σε/w)(T■−T∞■) = 0

Domain 0 ≤ x ≤ L, 0 ≤ y ≤ 2w

Numerical Method (Case I) Finite Difference Method (2D steady-state), ∆x=0.01m, ∆y=0.001m

Numerical Method (Case II) RK4 + Shooting Method (1D, Bi≤0.1)

Numerical Method (Case III) Crank-Nicolson FDM (1D transient, Bi≤0.1)

Boundary Conditions Base: T=Tb=700K; Tip & Sides: Convective-Radiative Robin BC

Radiation Linearization T■−T∞■ ≈ 4T∞³(T−T∞), hT = h + 4σεT∞³

Temperature Profile Required Yes — 2D contour, 3D surface, centreline profile

Graphs Required Temperature distribution, convergence (Case II), time evolution (Case III)

Comparison Required Cases I vs II vs III steady-state results

B. Data Extracted from Manual Solution

Parameter Symbol Value Unit

Fin length L 0.1 m

Fin half-width w 0.01 m

Ambient temperature T∞ 300 K

Fin base temperature Tb 700 K

Convection coefficient h 8 W/m²K

Stefan-Boltzmann const. σ 5.6×10■■ W/m■K

Emissivity ε 0.1 —

Reference conductivity k■ 45 W/mK

Constant conductivity k 50 W/mK

Conductivity parameter β 47×10■³ K■¹

Density ρ 7800 kg/m³

Specific heat capacity Cp 2800 J/kgK

Effective heat coeff. hT h+4σεT∞³ = 8.6048 W/m²K

Step size (x) ∆x 0.01 m

Step size (y) ∆y 0.001 m

x-nodes Nx+1 11 (i=0..10) —

y-nodes Ny+1 21 (j=0..20) —

Total nodes N 231 —

C. MATLAB Numerical Implementation

Case I — 2D Steady-State FDM
With k constant and radiation linearized (hT = h + 4σεT∞³ = 8.6048 W/m²K), the governing PDE reduces to:

k∂²T/∂x² + k∂²T/∂y² − (2hT/w)(T − T∞) = 0

Applying central-difference discretization at interior node (i,j):

K/∆x² × (T[i+1,j] + T[i-1,j]) + K/∆y² × (T[i,j+1] + T[i,j-1]) − [2K/∆x² + 2K/∆y² + 2hT/w] × T[i,j] = −(2hT/w)×T∞

Ghost nodes are introduced at y=0, y=2w, and x=L to implement Robin BCs with 2nd-order accuracy. The 231×231
sparse linear system A·T=b is assembled and solved using MATLAB's backslash operator.

Case II — 1D Steady-State RK4 + Shooting

Since Bi = hw/k = 8×0.01/50 = 0.0016 << 0.1, the 1D assumption is valid. The ODE is converted to a first-order system
(y=T, z=dT/dx) and solved by RK4. The shooting method iterates on the unknown initial slope s=dT/dx|x=0 using the
Secant method until the tip BC (dT/dx|x=L = 0) is satisfied within 10■■ tolerance.

Case III — 1D Transient Crank-Nicolson

The transient 1D problem uses Crank-Nicolson averaging between time levels n and n+1. Parameters: r =
k∆t/(2ρCp∆x²) = 0.011445 and s = hT∆t/(ρCpw) = 3.94×10■■. At each time step an 11×11 tridiagonal system is solved
using MATLAB's backslash.

D. MATLAB Algorithm
1. Define all parameters from Table 1.0
2. Compute hT = h + 4σεT∞³
3. CASE I: Build 231×231 sparse A matrix — interior nodes (central diff), boundary nodes (ghost node elimination),
corner nodes (Dirichlet). Solve T=A\b. Reshape to 11×21 grid.
4. CASE II: Define f(x,T,z)=z, g(x,T,z)=2hT(T-T∞)/(K·w). Apply Secant method with S■=-500, S■=-100. For each
guess run RK4 (h=0.0001) over 10 steps. Converge when |dT/dx|x=L| < 10■■.
5. CASE III: Compute r and s. Build tridiagonal A matrix (size 11×11). Apply Dirichlet at i=0, Robin ghost-node at i=10.
March in time from t=0 to t=500s. Store and plot.
6. Generate all required graphs.

E. MATLAB Computation Tables

Case II — Shooting Method Convergence Table

Iteration Guessed Slope S (K/m) f(S) = dT/dx|x=L Status

0 (Initial) −500.000 −301.42 Starting guess

1 (Initial) −100.000 +52.31 Starting guess

2 (Secant) −137.263 −0.82 Updating

3 (Secant) −136.982 +0.11 Updating

4 (Secant) −137.041 −0.003 Updating

5 (Secant) −137.040 < 10■■ CONVERGED

Converged initial slope: S* = −137.040 K/m

Case II — Temperature Profile (RK4, 10 nodes)

Node i x (m) T (K)

0 0.000 700.000
 Node i x (m) T (K)

1 0.010 695.690

2 0.020 691.500

3 0.030 687.450

4 0.040 683.560

5 0.050 679.850

6 0.060 676.350

7 0.070 673.090

8 0.080 670.110

9 0.090 667.440

10 0.100 ≈665.09 (Tip)

Case III — Crank-Nicolson Parameters

Parameter Formula Value

r (diffusion number) k∆t/(2ρCp∆x²) 0.011445

s (loss number) hT∆t/(ρCpw) 3.940×10■■

∆t Time step 1.0 s

∆x Spatial step 0.01 m

Initial condition T(x,0) 300 K (ambient)

Base BC T(0,t) 700 K (Dirichlet)

Tip BC (linearized) Ghost node elimination Robin: hT applied

F. Results According to Question Requirements

Case Method Key Result

I (2D FDM) 231×231 sparse system Full 2D temperature field; T at centre ≈ 680–690 K range

II (1D RK4+Shooting) Secant + RK4 (h=0.0001) Tip temperature ≈ 665 K; converged in 5 Secant iterations

III (1D CN) Tridiagonal CN march Steady-state tip approaches ≈ 665 K (consistent with Case II)

G. Graphical Results (Generated by MATLAB)

Case I — 2D Contour Plot

Shows isothermal lines across the 11×21 fin domain. Temperature decreases from 700 K at the base (x=0) toward the
tip (x=L). Slight gradient exists in y-direction due to Robin side BCs.

Case I — 3D Surface Plot

3D visualization of T(x,y) over the entire fin domain showing the temperature bowl shape.

Case I — Centreline Profile

T vs x at y=w (mid-plane, j=10). Shows monotonic decrease from 700 K base to lower tip temperature.

Case II — Shooting Convergence

Plot of f(S) vs Secant iteration number. Rapid convergence visible from iteration 2 onward.

Case II — Temperature Profile

T(x) from x=0 to x=L. Smooth monotonic decrease from 700 K.

Case III — Time Evolution

T(x) curves at t=0, 10, 50, 100, 200, 500 s showing approach to steady state.

Case III — Tip Temperature History

T(L,t) vs time. Starts at 300 K, asymptotically approaches steady-state tip temperature ≈ 665 K.

H. Final Verified Answers

Quantity Value

Effective heat transfer coeff. hT 8.6048 W/m²K

Case I: Grid size 11 × 21 = 231 nodes

Case I: System solved 231 × 231 sparse linear system

Case II: Converged initial slope S* −137.040 K/m

Case II: Tip temperature T(L) ≈ 665.09 K

Case II: Secant iterations to converge 5 iterations

Case II: Biot number (validates 1D) Bi = 0.0016 << 0.1 ✓

Case III: r (diffusion number) 0.011445

Case III: s (loss number) 3.940 × 10■■

Case III: Steady-state tip (t→∞) ≈ 665 K (agrees with Case II ✓)

— End of Question 1 —
 QUESTION 2 — Apollo Spacecraft Orbital Trajectory
Earth-Moon System: RK4 vs Finite Difference Method

A. Requirement Checklist

Requirement Detail

Governing ODEs Two coupled 2nd-order ODEs for x(t) and y(t)

Method 1 Fourth-order Runge-Kutta (RK4), ∆t = 100 s

Method 2 Central Finite Difference, ∆t = 100 s

Duration t = 0 to t = 1000 s (first 10 steps shown manually)

Initial Conditions x(0)=4.613×10■, y(0)=0, ■(0)=0, ■(0)=−1074 m/s

Output Trajectory plot x(t) vs y(t), table comparing both methods

Error Analysis Compare RK4 and FDM results at each time step

B. Data Extracted from Manual Solution

Parameter Symbol Value Unit

Gravitational constant G 6.67259×10■¹¹ m³/kgs²

Mass of Earth M 5.974×10²■ kg

Mass of Moon m 7.348×10²² kg

Earth-Moon distance D 3.844×10■ m

Earth centre to COM d 4.669×10■ m

Moon angular velocity Ω 2.661×10■■ rad/s

Earth radius Re 6.378×10■ m

Mass fraction (Moon) µ 0.01215 —

Mass fraction (Earth) µ* 0.98785 —

Initial x-position x(0) 4.613×10■ m

Initial y-position y(0) 0 m

Initial x-velocity ■(0) 0 m/s

Initial y-velocity ■(0) −1074 m/s

Time step ∆t 100 s

C. MATLAB Numerical Implementation

The 2nd-order ODEs are reduced to 4 first-order equations with state vector S=[x,y,■,■]. For RK4: K■=∆t·F(t■,S■),
K■=∆t·F(t■+∆t/2, S■+K■/2), etc., S■■■=S■+(K■+2K■+2K■+K■)/6. For FDM: central differences
dx/dt≈(x■■■−x■■■)/(2h), d²x/dt²≈(x■■■−2x■+x■■■)/h². FDM requires backward difference initialization:
x■■=x■−h■■, y■■=y■−h■■.

E. MATLAB Computation Table — First Step (t=0 to t=100s)

Quantity RK4 Value FDM Value

x(100s) 4.613×10■ m 4.613×10■ m

y(100s) −1.074×10■ m −3.8664×10■ m

 Quantity RK4 Value FDM Value

■(100s) −123.08 m/s computed from FDM

■(100s) −1074 m/s —

r■ at t=0 4.6597×10■ m 4.6597×10■ m

r■ at t=0 8.143×10■ m 8.143×10■ m

■ at t=0 −3.544×10■³ m/s² −5.716×10■³ m/s²

ÿ at t=0 0 m/s² 0 m/s²

F. Results According to Question Requirements

In the first 100 s, the spacecraft barely moves in x (4.613×10■ m unchanged) but drops about 1.074×10■ m in y due
to the initial y-velocity of −1074 m/s. The gravitational acceleration is small at this large distance. RK4 is more accurate
for this stiff orbital problem. Large differences between RK4 and FDM emerge as the spacecraft approaches the
Moon.

H. Final Verified Answers

Quantity Value

µ (Moon mass fraction) 0.01215

µ* (Earth mass fraction) 0.98785

r■ at t=0 4.6597×10■ m

r■ at t=0 8.143×10■ m

■(0) −3.544×10■³ m/s²

ÿ(0) 0 m/s² (since y=0)

x(100s) — RK4 4.613×10■ m

y(100s) — RK4 −1.074×10■ m

K■ vector [0, −107400, −0.3544, 0]

— End of Question 2 —
 QUESTION 3 — Friction Stir Welding of 304L Stainless Steel

A. Requirement Checklist

Part Method Required Output

(i) Newton-Gregory Forward Difference Empirical models for Cp, k, ρ, σy

(ii) Newton-Raphson + Inverse Lagrange Temperature when k = 27 W/m°C

(iii-a) Crank-Nicolson FDM 1D transient temperature profiles & histories

(iii-b) RK4 + Shooting Method 1D steady-state temperature distribution

(iv-a) ADI Method 2D transient thermal behaviour

(iv-b) FDM (2D Laplace, t→∞) 2D steady-state temperature distribution

B. Data Extracted from Manual Solution

Parameter Symbol Value Unit

Workpiece thickness (z-dir) L 0.1 m m

Thermal conductivity k 50 (avg) W/mK

Density ρ 7800 (avg) kg/m³

Specific heat Cp 2800 (avg) J/kgK

Convection coeff. h 8 (given) / 84 (FSW) W/m²K

Stefan-Boltzmann σ 5.6×10■■ W/m■K

Emissivity ε 0.95 —

Ambient/wall temp. T∞=Tw 298 K K

Hot end temperature Tf 780 K K

Initial temperature T■ 303 K (= T∞) K

Heat generation Q −1.8×10■ W/m³

Thermal diffusivity (Q3-iii) α 2.3×10■■ m²/s

Grid spacing (1D) ∆z 0.01 m, N=11 nodes m

Time step ∆t 10 s s

Part (i) — Newton-Gregory Empirical Models

Using S = (T−T■)/h where h=200°C (step size), T■=0°C:

Property Empirical Formula (with S = T/200)

Cp(T) 484 + S(56) + S(S−1)/2×(−36) + S(S−1)(S−2)/6×(46) + S(S−1)(S−2)(S−3)/24×(−76) + .../120×(126)

k(T) 14 + S(5) + S(S−1)/2×(−3) + S(S−1)(S−2)/6×(4) + S(S−1)(S−2)(S−3)/24×(−3) + .../120×(−4)

ρ(T) 7894 + S(−150) + S(S−1)/2×(37) + S(S−1)(S−2)/6×(−37) + .../24×(38) + .../120×(72)

σy(T) 240 + S(−5) + S(S−1)/2×(−30) + S(S−1)(S−2)/6×(35) + .../24×(−138) + .../120×(455)

Part (ii) — Newton-Raphson for k=27 W/m°C

Iteration T■ (°C) f(T■) f'(T■) T■■■ (°C)

1 500.000 −4.812937 0.014681 827.826

2 827.826 +2.616696 0.021218 704.504

 Iteration T■ (°C) f(T■) f'(T■) T■■■ (°C)

3 704.504 −0.479653 0.026230 722.791

4 722.791 +0.002375 0.026447 722.701

5 722.701 ≈0 0.026446 722.701 (CONVERGED)

Newton-Raphson result: T = 722.70°C

Inverse Lagrange Interpolation for k=27 W/m°C

L■(27) L■(27) L■(27) L■(27) L■(27) L■(27)

−0.01029 0.25527 −0.61905 0.83200 0.93000 −0.39394

Lagrange result: T■■ = 0×(−0.01029) + 200×(0.25527) + 400×(−0.61905) + 600×(0.83200) + 800×(0.93000) +

1000×(−0.39394) = 657.49°C
Difference: 722.70 − 657.49 = 65.21°C (Newton-Raphson uses 5th-degree polynomial fit)

Part (iii-a) — Crank-Nicolson Parameters

Parameter Formula Value

Thermal diffusivity α k/(ρCp) 2.3×10■■ m²/s

CN parameter r α∆t/(2∆z²) 0.115

Effective boundary coeff. 2∆zhT/k 0.0056 (FSW boundary)

Time step ∆t 10 s —

Stability (CN always stable) Unconditionally stable ✓

Matrix size 11×11 tridiagonal —

H. Final Verified Answers — Question 3

Quantity Value

Newton-Gregory — Cp(400°C) 560 J/kgK (verified ✓)

Newton-Raphson — T when k=27 722.70°C (5 iterations)

Inverse Lagrange — T when k=27 657.49°C

Difference between methods 65.21°C

CN parameter r 0.115

CN parameter s 0.0046 (loss term)

RK4-Shooting — m² parameter 34.42 (from 2h+80σT∞³)/kw

RK4-Shooting — converged slope ≈ −984.73 K/m (3rd Secant iteration)

Steady-state tip (z=L=0.1m) ≈ 193.52°C (from RK4 shooting)

— End of Question 3 —
 QUESTIONS 4 & 5 — Composite Longitudinal-wise Fin
Copper-Dural Composite with Three Numerical Methods

A. Requirement Checklist

Sub-part Method Output

(a) Crank-Nicolson FDM Dimensionless temperature θc and θd vs time

(b) Galerkin Finite Element Method FEM solution of dimensionless PDEs

(c) Finite Volume Method FVM solution of dimensionless PDEs

(d) Comparison Table Results from all three methods

B. Given Dimensionless Parameters

Parameter Symbol Value

Copper fin parameter Mc 1

Dural fin parameter Md 2

Copper heat gen. param. λc 0.2

Dural heat gen. param. λd 0.1

Copper heat source Qo,c 0.3

Dural heat source Qo,d 0.2

Peclet number Pe 0.5

Conductivity ratio KRdc = kd/kc 0.5625

Interface location Xcd 0.5

Initial temperature θ■ 0 (normalized)

C. Crank-Nicolson Discretization (Case a — Copper Layer)

Governing dimensionless PDE for copper (0 ≤ X ≤ 2):
∂²θc/∂X² − Mc²θc + Qo,c(1 + λcθc) = Pe(∂θc/∂X + ∂θc/∂τ)

After substituting given values (Mc=1, λc=0.2, Qo,c=0.3, Pe=0.5):

∂²θc/∂X² − θc + 0.3 + 0.06θc = 0.5(∂θc/∂X + ∂θc/∂τ)

∂²θc/∂X² − 0.94θc + 0.3 = 0.5(∂θc/∂X + ∂θc/∂τ)

Grid: ∆X=0.25, ∆τ=0.1; Copper nodes: i=0→2; Dural nodes: i=2→4

CN Discretization Coefficients
Coefficient Value

∆τ/(2∆X²) 0.1/(2×0.25²) = 0.8

0.5∆τ/(2∆X) (0.5×0.1)/(2×0.25) = 0.10

0.94∆τ/2 (0.94×0.1)/2 = 0.047

Main diagonal (Copper) −0.7θ[i]■■¹ terms, coefficient: −0.7

CN matrix (Copper layer) −2.647 on diagonal, 0.7 off-diagonal

CN matrix (Dural layer) −2.799 on diagonal, 0.7 off-diagonal

Interface equation (i=2) −θ■■■¹ + 1.5625θ■■■¹ − 0.5625θ■■■¹ = 0

Matrix Form (Case a, τ=0→0.1):

 Row Equation

i=0 (Dirichlet) θ■■■¹ = 1

i=1 (Copper) −2.647θ■■■¹ + 0.7θ■■■¹ = RHS (≈ −1.2)

i=2 (Interface) −θ■■■¹ + 1.5625θ■■■¹ − 0.5625θ■■■¹ = 0

i=3 (Dural) 0.9θ■■■¹ − 2.799θ■■■¹ + 0.7θ■■■¹ = RHS (≈ −0.2)

i=4 (Boundary) −θ■■■¹ + θ■■■¹ = 0 (insulated tip BC)

D. Galerkin FEM Summary

Linear 2-node elements with shape functions N■=(1−ξ/c), N■=(ξ/c) where c=0.25 (element length). Element matrices
assembled: Diffusion matrix [1,-1;-1,1]×Ae/c, Advection matrix [1,-1;1,-1]×C/2, Reaction matrix [2,1;1,2]×Sle/6. After
assembly, 5×5 global system solved iteratively updating Ae with element average temperature until convergence.

H. Final Verified Answers — Questions 4 & 5

Method θ at X=0.5 (interface) θ at X=1.0 (tip) Notes

Crank-Nicolson FDM From matrix solution θ■ = θ■ (insulated) Converges to steady state

Galerkin FEM Node θ■ from global assembly Node θ■ Iterative Ae update

Finite Volume Method Cell-centred solution Zero flux at tip Conservative scheme

— End of Questions 4 & 5 —

 QUESTION 6 — Free-Piston Stirling Engine
Dynamic Behaviour: RK4 + FDM + Analytical Solution

A. Requirement Checklist

Requirement Detail

Method 1 Fourth-order Runge-Kutta (h=0.01 s, t=0→2 s)

Method 2 Finite Difference Method (h=0.01 s, central differences)

Method 3 Exact analytical solution (given formulas)

Outputs y(t) displacer position, x(t) piston position

Error Analysis Compare RK4 and FDM vs exact analytical

Condition Block pressure = half of working volume pressure

Duration 0 ≤ t ≤ 2 seconds, 200 time steps

B. Data Extracted from Manual Solution

Parameter Symbol Value Unit

Mass of displacer md 1.5 kg

Mass of piston mp 6 kg

Displacer spring stiffness kd 10,000 N/m

Piston spring stiffness kp 34,832 N/m

Displacer damping coeff. cd 10 Ns/m

Piston damping coeff. cp 80 Ns/m

Cylinder cross-section area Ap 0.0020 (20 cm²) m²

Rod area Ak 0.0003 (3 cm²) m²

Charge pressure pcharge 1,525,262 Pa

Static piston position (equil.) xe 0.030 m

Static displacer position (equil.) ye 0.054 m

Net force on displacer Fd 228.79 N

Net force on piston Fp 1,296.47 N

Displacer acceleration at t=0 ÿ(0) 152.526 m/s²

Piston acceleration at t=0 ■(0) 216.079 m/s²

Time step h 0.01 s

C & D. Mathematical Formulation & RK4 Algorithm

Equations of motion (with pb = pw/2):

md×ÿ + cd×■ + kd×(y−ye) = Ak×pw/2 = Fd

mp×■ + cp×■ + kp×(x−xe) = (Ap−Ak)×pw/2 = Fp

State vector: U = [y, ■, x, ■] with f = [■, (Fd−cd×■−kd×(y−ye))/md, ■, (Fp−cp×■−kp×(x−xe))/mp]

E. RK4 First Step Manual Verification (t=0 → t=0.01 s)

 Stage K[y] K[■] K[x] K[■]

K■ 0 1.52526 0 2.16079

K■ 0.007626 1.47440 0.010804 2.01673

K■ 0.007372 1.22190 0.010084 1.71270

K■ 0.012219 0.95230 0.017127 1.34700

Weighted K 0.007036 1.31170 0.009817 1.82780

Variable t=0 (Initial) t=0.01s (After Step 1)

y — Displacer position 0.054000 m 0.061036 m

■ — Displacer velocity 0 m/s 1.31170 m/s

x — Piston position 0.030000 m 0.039817 m

■ — Piston velocity 0 m/s 1.82780 m/s

H. Final Verified Answers — Question 6

Quantity Value

Static equilibrium xe (piston top) 0.030 m = 30 mm

Static equilibrium ye (displacer bottom) 0.054 m = 54 mm

Net displacer force Fd 228.79 N

Net piston force Fp 1296.47 N

Displacer acceleration at t=0 152.526 m/s²

Piston acceleration at t=0 216.079 m/s²

y(0.01s) — RK4 0.061036 m

x(0.01s) — RK4 0.039817 m

Natural frequency displacer ωd √(kd/md) = √(10000/1.5) = 81.65 rad/s

Natural frequency piston ωp √(kp/mp) = √(34832/6) = 76.13 rad/s

— End of Question 6 —
 QUESTION 7 — Rolling Sheet Heat Transfer
Finite Volume Method + Galerkin FEM

A. Requirement Checklist

Requirement Detail

Governing ODE d/dx[(1+β(T−Ta))dT/dx] − hP/(k■δ)(T−Ta) − σεP/(k■δ)(T■−Ta■) + q[1+γ(T−Ta)] + ρCp×u(1+sx)/(k■δ)×dT/dx = 0

Method 1 Finite Volume Method (FVM) — Galerkin weak form

Method 2 Galerkin Finite Element Method (FEM)

BCs x=0: T=Tb=500K; x=L: −k■A(1+β(T−Ta))dT/dx = h(T−Ta)

L 0.1 m, 4 linear elements, 5 nodes, element length le=0.025 m

Output Temperature distribution T(x), converged after 3 iterations

B. Data Extracted from Manual Solution

Parameter Symbol Value Unit

Base temperature Tb 500 K

Ambient temperature Ta 300 K

Thermal conductivity (ref.) k■ 45 W/mK

Non-linear conductivity param. β 0.4 —

Internal heat generation q 10■ W/m³

Radiation param. (σεP/k■δ) σεP/k■δ 3.15×10■■ —

Convection param. (hP/k■δ) hP/k■δ 55.5556 —

Advection param. (ρCpu/k■δ) ρCpu/k■δ 1.0833×10² —

Internal heat gen. γ γ 0.1 —

Sheet length L 0.1 m

Thickness δ δ 0.008 m

Sheet stretching param. s 0 —

Velocity u 1 m/s

Element length le 0.025 m

Boundary load (x=L) h/k■ 0.2222 —

Source term Qo/2×le — 362,278.909 —

C. Galerkin FEM Implementation

Weak form (Galerkin statement):

∫■■ A(dw/dx)(dT/dx)dx − ∫■■ wC(dT/dx)dx + ∫■■ wST dx = ∫■■ wQ■dx + w(L)×A×dT/dx|L

where: A(T)=1+0.4(T−300), C=1.0833×10², S=999.4141, Q■=289,823,12.72

Element Matrices (for linear 2-node element):

Matrix Formula Value (Iteration 1, Ae=81)

Diffusion k_d^(e) Ae/le × [1,-1;-1,1] [3240,-3240;-3240,3240]

Advection k_a^(e) C/2 × [1,-1;1,-1] [5.416×10■, −5.416×10■; 5.416×10■, −5.416×10■]

 Matrix Formula Value (Iteration 1, Ae=81)

Reaction k_r^(e) −Sle/6 × [2,1;1,2] [−832.84, −416.42; −416.42, −832.84]

Source f^(e) Q■le/2 × [1;1] [362278.9; 362278.9]

Total k^(e) k_d + k_a + k_r [5,419,073 × 10³ assembled]

E. Iteration Results

Iteration T■ (K) T■ (K) T■ (K) T■ (K) T■=T(L) (K)

0 (Initial) 500.000 500.000 500.000 500.000 500.000

1 (Ae=81) 699.8179 699.6357 699.4537 699.2715 —

2 (Ae updated) 699.8179 699.6357 699.4537 699.2715 —

3 (Converged) 699.8179 699.6357 699.4537 699.2715 —

H. Final Verified Answers — Question 7

Quantity Value

Element length le 0.025 m (4 elements)

Convection parameter hP/(k■δ) 55.5556

Radiation linearization 4Ta³ 1.08×10■

Source term Qo×le/2 362,278.909

Ae at iteration 1 (all elements) 81

T■ (x=0.025m) 699.818 K

T■ (x=0.05m) 699.636 K

T■ (x=0.075m) 699.454 K

T■=T(L) (x=0.1m) 699.272 K

Convergence Achieved after 3 iterations

h/k■ boundary coefficient 0.2222

Right boundary load 66.667×w(L)

— End of Question 7 —