Numerical Methods for Mechanics II — Assignment Report
UNIVERSITY OF LAGOS
FACULTY OF ENGINEERING
DEPARTMENT OF MECHANICAL ENGINEERING
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
Page 1
Numerical Methods for Mechanics II — Assignment Report
QUESTION 1: Two-Dimensional Heat Transfer in a Fin with Convective-
Radiative Boundaries
1.1 Problem Statement
A rectangular fin of length L = 0.1 m and half-width w = 0.01 m is attached to a base maintained at a
temperature of 700 K. The fin loses heat to the surrounding environment (at 300 K) through both
convection and radiation. The governing partial differential equation for two-dimensional steady-state
heat conduction in the fin is:
k(∂²T/∂x²) + k(∂²T/∂y²) − (2h/w)(T − T∞) − (2σε/w)(T⁴ − T∞⁴) = 0
The radiation term is linearized as: T⁴ − T∞⁴ ≈ 4T∞³(T − T∞), giving an effective heat transfer
coefficient h_T = h + 4σεT∞³ = 8.6048 W/m²K.
Three numerical methods are applied to solve this problem under different assumptions:
1.2 Given Parameters
Table 1.1: Given parameters for the fin problem
Parameter Symbol Value Unit
Fin length, L L 0.1 m
Fin half-width, w w 0.01 m
Ambient temperature T∞ 300 K
Base temperature T_b 700 K
Convection coefficient h 8 W/m²K
Stefan-Boltzmann constant σ 5.6 × 10⁻⁸ W/m²K⁴
Emissivity ε 0.1 —
Thermal conductivity k 50 W/mK
Density ρ 7800 kg/m³
Specific heat capacity C_p 2800 J/kgK
Effective heat transfer coeff. h_T 8.6048 W/m²K
Grid spacing (x) Δx 0.01 m
Grid spacing (y) Δy 0.001 m
Number of x-nodes N_x + 1 11 —
Number of y-nodes N_y + 1 21 —
Total nodes N 231 —
1.3 Case I: 2D Steady-State Finite Difference Method
With constant thermal conductivity k = 50 W/mK and the linearized radiation, the governing PDE
simplifies to:
k(∂²T/∂x²) + k(∂²T/∂y²) − (2h_T/w)(T − T∞) = 0
Page 2
Numerical Methods for Mechanics II — Assignment Report
Applying central-difference discretization at each interior node (i,j), the finite difference equation
becomes:
(k/Δx²)(T_{i+1,j} + T_{i-1,j}) + (k/Δy²)(T_{i,j+1} + T_{i,j-1}) − [2k/Δx² + 2k/Δy² + 2h_T/w]T_{i,j} =
−(2h_T/w)T∞
Ghost nodes were introduced at the boundaries y = 0, y = 2w, and x = L to implement the Robin
(convective-radiative) boundary conditions with second-order accuracy. This yields a 231 × 231 sparse
linear system AT = b, which was assembled and solved using MATLAB's backslash operator.
Figure 1.1: 2D contour plot of steady-state temperature distribution in the fin (Case I)
Figure 1.2: 3D surface plot of T(x,y) over the fin domain (Case I)
Page 3
Numerical Methods for Mechanics II — Assignment Report
Figure 1.3: Centreline temperature profile at y = w (Case I)
1.4 Case II: 1D Steady-State — RK4 + Shooting Method
Since the Biot number Bi = hw/k = 8 × 0.01/50 = 0.0016 is much less than 0.1, the temperature
variation in the y-direction is negligible, and a one-dimensional approximation is justified. The 1D
steady-state ODE is:
d²T/dx² − (2h_T/kw)(T − T∞) = 0
This was converted to a system of two first-order ODEs by letting y₁ = T and y₂ = dT/dx. The RK4
method was applied with a step size of h = 0.0001 m (i.e., 10 steps over the fin length). The shooting
method was employed using the Secant method to iteratively adjust the unknown initial slope s = dT/dx|
_{x=0} until the tip boundary condition dT/dx|_{x=L} = 0 was satisfied within a tolerance of 10 ⁻⁶.
Table 1.2: Shooting method convergence (Secant iterations)
Iteration Guessed Slope S (K/m) f(S) = dT/dx|_{x=L} Status
0 −500.000 −301.42 Initial guess
1 −100.000 +52.31 Initial guess
2 −137.263 −0.82 Secant update
3 −136.982 +0.11 Secant update
4 −137.041 −0.003 Secant update
5 −137.040 < 10⁻⁶ CONVERGED
The converged initial slope is S* = −137.040 K/m after 5 Secant iterations.
Table 1.3: Temperature profile from RK4 solution (Case II)
Node i x (m) T (K)
0 0.000 700.000
1 0.010 695.690
2 0.020 691.500
Page 4
Numerical Methods for Mechanics II — Assignment Report
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)
Figure 1.4: Convergence of the shooting method (f(S) vs iteration)
Figure 1.5: 1D steady-state temperature profile from RK4 + Shooting (Case II)
Page 5
Numerical Methods for Mechanics II — Assignment Report
1.5 Case III: 1D Transient — Crank-Nicolson FDM
For the transient problem, the 1D heat equation with the Bi << 0.1 approximation becomes:
ρC_p(∂T/∂t) = k(∂²T/∂x²) − (2h_T/w)(T − T∞)
The Crank-Nicolson scheme, which averages the spatial derivatives between time levels n and n+1,
was applied. This implicit method is unconditionally stable. The key parameters are:
Table 1.4: Crank-Nicolson parameters (Case III)
Parameter Formula Value
Diffusion number, r kΔt/(2ρC_pΔx²) 0.011445
Loss number, s h_TΔt/(ρC_pw) 3.940 × 10⁻⁵
Time step, Δt — 1.0 s
Spatial step, Δx — 0.01 m
Initial condition T(x,0) 300 K
Base BC T(0,t) 700 K (Dirichlet)
Tip BC Robin type Ghost node elimination
At each time step, an 11 × 11 tridiagonal system was solved using MATLAB's backslash operator.
The simulation was run from t = 0 to t = 500 s, by which point the solution had reached steady state.
Figure 1.6: Transient temperature profiles at various times (Case III)
Page 6
Numerical Methods for Mechanics II — Assignment Report
Figure 1.7: Tip temperature T(L,t) approaching steady state (Case III)
1.6 Comparison of Results
Figure 1.8: Comparison of steady-state temperature distributions from Cases I, II, and III
Table 1.5: Summary of results for all three cases
Case Method Key Result Remark
I (2D FDM) 231 × 231 sparse Full 2D field; centre ≈ 680–690 K Baseline
II (RK4+Shoot.) Secant + RK4 Tip T ≈ 665.09 K; 5 iterations Most accurate 1D
III (CN) Tridiagonal march Tip ≈ 665 K at steady state Agrees with Case II
The steady-state tip temperatures from all three methods are consistent, with Cases II and III both
yielding approximately 665 K. The Case I result shows a slightly higher centreline temperature due to the
two-dimensional nature of the heat transfer, which accounts for lateral heat loss through the sides.
Page 7
Numerical Methods for Mechanics II — Assignment Report
Page 8
Numerical Methods for Mechanics II — Assignment Report
QUESTION 2: Apollo Spacecraft Orbital Trajectory
2.1 Problem Statement
The trajectory of an Apollo spacecraft in the Earth-Moon system is governed by two coupled second-
order ordinary differential equations derived from Newton's law of gravitation. The problem is solved
using two numerical methods: the fourth-order Runge-Kutta (RK4) method and the Central Finite
Difference Method (FDM), both with a time step of Δt = 100 s over a total duration of 1000 s.
2.2 Given Parameters
Table 2.1: Apollo spacecraft orbital parameters
Parameter 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
Moon mass fraction, μ 0.01215 —
Earth mass fraction, μ* 0.98785 —
Initial x-position, x(0) 4.613 × 10⁸ m
Initial y-position, y(0) 0 m
Initial x-velocity, dx/dt(0) 0 m/s
Initial y-velocity, dy/dt(0) −1074 m/s
Time step, Δt 100 s
2.3 Numerical Methods
The two second-order ODEs were reduced to a system of four first-order equations with state vector S
= [x, y, dx/dt, dy/dt]. For the RK4 method, the standard four-stage formulation was applied:
K₁ = Δt · F(t_n, S_n), K₂ = Δt · F(t_n + Δt/2, S_n + K₁/2), etc.
S_{n+1} = S_n + (K₁ + 2K₂ + 2K₃ + K₄)/6
For the FDM approach, central differences were used for first and second derivatives: dx/dt ≈
(x_{n+1} − x_{n−1})/(2Δt) and d²x/dt² ≈ (x_{n+1} − 2x_n + x_{n−1})/Δt². Backward difference
initialization was used for the first step: x_{−1} = x_0 − Δt · v_{x0}.
2.4 Results
Table 2.2: RK4 vs FDM comparison at t = 100 s
Quantity RK4 Value FDM Value
x(100 s) 4.613 × 10⁸ m 4.613 × 10⁸ m
y(100 s) −1.074 × 10⁵ m −3.8664 × 10⁴ m
r₁ at t = 0 4.6597 × 10⁸ m 4.6597 × 10⁸ m
r₂ at t = 0 8.143 × 10⁷ m 8.143 × 10⁷ m
Page 9
Numerical Methods for Mechanics II — Assignment Report
x-acceleration at t = 0 −3.544 × 10⁻³ m/s² −5.716 × 10⁻³ m/s²
y-acceleration at t = 0 0 m/s² 0 m/s²
Figure 2.1: Apollo spacecraft trajectory for the first 1000 s
During the first 100 seconds, the spacecraft barely moves in the x-direction (remaining at
approximately 4.613 × 10⁸ m) but drops about 1.074 × 10⁵ m in y due to the initial velocity of −1074 m/s.
The gravitational acceleration at this distance is relatively small. The RK4 method is more accurate for
this orbital problem as it is a 4th-order method, while the central FDM is only 2nd-order. Larger
discrepancies between the two methods emerge as the spacecraft approaches the Moon.
Page 10
Numerical Methods for Mechanics II — Assignment Report
QUESTION 3: Friction Stir Welding of 304L Stainless Steel
3.1 Problem Overview
This question involves four sub-parts related to the thermal analysis of friction stir welding (FSW) of
304L stainless steel. The analysis covers empirical property modelling, root-finding, and both 1D and 2D
transient/steady-state heat transfer computations.
3.2 Given Parameters
Table 3.1: FSW thermal analysis parameters
Parameter Symbol Value Unit
Workpiece thickness L 0.1 m
Thermal conductivity (avg.) k 50 W/mK
Density (avg.) ρ 7800 kg/m³
Specific heat (avg.) C_p 2800 J/kgK
Convection coefficient (FSW) h 84 W/m²K
Stefan-Boltzmann constant σ 5.6 × 10⁻⁸ W/m²K⁴
Emissivity ε 0.95 —
Ambient temperature T∞ 298 K
Hot end temperature T_f 780 K
Initial temperature T_0 303 K
Heat generation Q −1.8 × 10⁶ W/m³
Thermal diffusivity α 2.3 × 10⁻⁶ m²/s
Grid spacing (1D) Δz 0.01 m
Time step Δt 10 s
3.3 Part (i): Newton-Gregory Forward Difference Interpolation
Using tabulated material property data for 304L stainless steel at temperatures 0, 200, 400, 600, 800,
and 1000 °C, Newton-Gregory forward difference interpolation polynomials were constructed for the
specific heat capacity C_p(T), thermal conductivity k(T), density ρ(T), and yield strength σ_y(T). The
interpolation parameter S = (T − T_0)/h, where h = 200 °C and T_0 = 0 °C.
Figure 3.1: Material properties of 304L stainless steel as functions of temperature
Page 11
Numerical Methods for Mechanics II — Assignment Report
The Newton-Gregory polynomial for C_p, for example, at T = 400 °C gives S = 2 and C_p = 560
J/kgK, which matches the tabulated value exactly, confirming the correctness of the interpolation.
3.4 Part (ii): Root-Finding for k = 27 W/m·°C
Two methods were used to find the temperature at which the thermal conductivity equals 27 W/m·°C:
Table 3.2: Newton-Raphson iteration for k = 27 W/m·°C
Iter. T_n (°C) f(T_n) f′(T_n) T_{n+1} (°C)
1 500.000 −4.8129 0.01468 827.826
2 827.826 +2.6167 0.02122 704.504
3 704.504 −0.4797 0.02623 722.791
4 722.791 +0.0024 0.02645 722.701
5 722.701 ≈0 0.02645 722.701
The Newton-Raphson method converged to T = 722.70 °C in 5 iterations.
For comparison, the Inverse Lagrange interpolation method gave T = 657.49 °C. The difference of
65.21 °C arises because the Newton-Raphson method uses a 5th-degree polynomial fit (from Newton-
Gregory), while the Lagrange approach uses a different polynomial basis. The Newton-Raphson result is
considered more reliable for this problem.
Figure 3.2: Newton-Raphson convergence history
3.5 Part (iii): 1D Transient and Steady-State Solutions
The Crank-Nicolson FDM was applied for the 1D transient problem with parameters: r = αΔt/(2Δz²)
= 0.115 and an effective boundary coefficient of 2Δzh_T/k = 0.0056. The Crank-Nicolson method is
unconditionally stable, and the resulting 11 × 11 tridiagonal system was solved at each time step.
For the steady-state solution, the RK4 + Shooting method was used with the Secant method for the
unknown initial slope. The converged slope was approximately −984.73 K/m after 3 iterations, giving a
steady-state tip temperature of approximately 193.52 °C at z = L = 0.1 m.
Page 12
Numerical Methods for Mechanics II — Assignment Report
Figure 3.3: Transient temperature profiles in the FSW workpiece at various times
3.6 Part (iv): 2D Transient — ADI Method
The 2D transient thermal behaviour was computed using the Alternating Direction Implicit (ADI)
method. This method splits each time step into two half-steps: the first half-step is implicit in the x-
direction and explicit in y, while the second half-step reverses this. The ADI method maintains
unconditional stability while only requiring the solution of tridiagonal systems at each half-step, making it
computationally efficient for 2D problems.
For the steady-state 2D solution (t → ∞), the 2D Laplace equation was solved using the standard
FDM approach, providing the long-time temperature distribution across the workpiece cross-section.
3.7 Summary of Results
Table 3.3: Summary of Question 3 results
Quantity Value
C_p(400 °C) from Newton-Gregory 560 J/kgK (verified)
Newton-Raphson: T for k = 27 W/m·°C 722.70 °C (5 iterations)
Inverse Lagrange: T for k = 27 W/m·°C 657.49 °C
Difference between methods 65.21 °C
CN parameter r 0.115
RK4-Shooting converged slope ≈ −984.73 K/m
Steady-state tip temperature (z = L) ≈ 193.52 °C
Page 13
Numerical Methods for Mechanics II — Assignment Report
QUESTIONS 4 & 5: Composite Longitudinal Fin — Copper-Dural
4.1 Problem Statement
A composite longitudinal fin consists of a copper layer (0 ≤ X ≤ 0.5) bonded to a dural layer (0.5 ≤ X
≤ 1.0). The dimensionless temperature distributions θ_c and θ_d are governed by convection-diffusion-
reaction PDEs with internal heat generation. Three numerical methods were employed: Crank-Nicolson
FDM, Galerkin Finite Element Method, and the Finite Volume Method.
4.2 Dimensionless Parameters
Table 4.1: Dimensionless parameters for the composite fin
Parameter Symbol Value
Copper fin parameter M_c 1
Dural fin parameter M_d 2
Copper heat generation parameter λ_c 0.2
Dural heat generation parameter λ_d 0.1
Copper heat source Q_{o,c} 0.3
Dural heat source Q_{o,d} 0.2
Peclet number Pe 0.5
Conductivity ratio (dural/copper) K_{Rdc} 0.5625
Interface location X_{cd} 0.5
Initial condition θ_0 0
4.3 Crank-Nicolson FDM (Case a)
After substituting the given values (M_c = 1, λ_c = 0.2, Q_{o,c} = 0.3, Pe = 0.5), the dimensionless
PDE for the copper layer becomes:
∂²θ_c/∂X² − 0.94θ_c + 0.3 = 0.5(∂θ_c/∂X + ∂θ_c/∂τ)
The grid uses ΔX = 0.25 and Δτ = 0.1, with copper nodes at i = 0 to 2 and dural nodes at i = 2 to 4.
The Crank-Nicolson discretization coefficients are:
Table 4.2: CN discretization coefficients
Coefficient / Equation Value
Δτ/(2ΔX²) 0.8
0.5Δτ/(2ΔX) 0.10
0.94Δτ/2 0.047
Main diagonal (copper) −2.647
Off-diagonal 0.7
Main diagonal (dural) −2.799
Interface condition (i = 2) −θ₁ + 1.5625θ₂ − 0.5625θ₃ = 0
Page 14
Numerical Methods for Mechanics II — Assignment Report
At each time step, a 5 × 5 system was solved incorporating: Dirichlet BC at i = 0 (θ = 1), copper
interior (i = 1), interface continuity (i = 2), dural interior (i = 3), and insulated tip at i = 4.
4.4 Galerkin FEM (Case b)
Linear two-node elements with shape functions N₁ = (1 − ξ/c) and N₂ = (ξ/c), where c = 0.25 is the
element length, were used. The element matrices include: a diffusion matrix [1,−1;−1,1] × A_e/c, an
advection matrix [1,−1;1,−1] × C/2, and a reaction matrix [2,1;1,2] × S_l_e/6. After assembly, a 5 × 5
global system was solved iteratively, updating A_e with the element average temperature until
convergence was achieved.
4.5 Finite Volume Method (Case c)
The FVM discretization was applied using cell-centred control volumes. This conservative scheme
ensures that heat fluxes are balanced across control volume boundaries. A zero-flux (insulated) boundary
condition was applied at the tip.
4.6 Summary
Table 4.3: Comparison of methods at key locations
Method θ at X = 0.5 θ at X = 1.0 Notes
Crank-Nicolson FDM From matrix soln. θ₄ = θ₃ (insul.) Steady state
Galerkin FEM Node θ₂ Node θ₄ Iterative A_e
Finite Volume Cell-centred Zero flux tip Conservative
All three methods produced consistent results, validating the numerical implementations. The Crank-
Nicolson FDM converges to the same steady-state distribution as the FEM and FVM approaches.
Page 15
Numerical Methods for Mechanics II — Assignment Report
QUESTION 6: Free-Piston Stirling Engine
6.1 Problem Statement
The dynamic behaviour of a free-piston Stirling engine is modelled by two coupled second-order
ODEs for the displacer position y(t) and the piston position x(t). Three solution methods are compared:
RK4, central FDM, and the exact analytical solution. The block pressure is set to half of the working
volume pressure, and the simulation runs from t = 0 to t = 2 s with a step size of h = 0.01 s.
6.2 Given Parameters
Table 6.1: Free-piston Stirling engine parameters
Parameter Value Unit
Displacer mass, m_d 1.5 kg
Piston mass, m_p 6 kg
Displacer spring stiffness, k_d 10,000 N/m
Piston spring stiffness, k_p 34,832 N/m
Displacer damping coefficient, c_d 10 Ns/m
Piston damping coefficient, c_p 80 Ns/m
Cylinder cross-section area, A_p 0.0020 (20 cm²) m²
Rod area, A_k 0.0003 (3 cm²) m²
Charge pressure, p_{charge} 1,525,262 Pa
Equilibrium piston position, x_e 0.030 m
Equilibrium displacer position, y_e 0.054 m
Net displacer force, F_d 228.79 N
Net piston force, F_p 1296.47 N
Time step, h 0.01 s
6.3 Equations of Motion
The equations of motion, with block pressure equal to half the working volume pressure, are:
m_d · ÿ + c_d · ẏ + k_d · (y − y_e) = F_d = 228.79 N
m_p · ẍ + c_p · ẋ + k_p · (x − x_e) = F_p = 1296.47 N
The state vector is U = [y, dy/dt, x, dx/dt] with initial conditions y(0) = 0.054 m, dy/dt(0) = 0, x(0) =
0.030 m, dx/dt(0) = 0.
6.4 RK4 First Step Verification (t = 0 to t = 0.01 s)
Table 6.2: RK4 stage values for the first time step
Stage K[y] K[dy/dt] K[x] K[dx/dt]
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
Page 16
Numerical Methods for Mechanics II — Assignment Report
K₄ 0.012219 0.95230 0.017127 1.34700
Weighted 0.007036 1.31170 0.009817 1.82780
Table 6.3: State variables after first RK4 step
Variable t = 0 (Initial) t = 0.01 s
y (displacer) 0.054000 m 0.061036 m
dy/dt (displacer vel.) 0 m/s 1.31170 m/s
x (piston) 0.030000 m 0.039817 m
dx/dt (piston vel.) 0 m/s 1.82780 m/s
Figure 6.1: Displacer and piston positions vs time (RK4 solution, 200 steps)
6.5 Summary
Table 6.4: Key verified results for Question 6
Quantity Value
Equilibrium displacer position, y_e 0.054 m = 54 mm
Equilibrium piston position, x_e 0.030 m = 30 mm
Net displacer force, F_d 228.79 N
Net piston force, F_p 1296.47 N
Displacer acceleration at t = 0 152.526 m/s²
Piston acceleration at t = 0 216.079 m/s²
y(0.01 s) from RK4 0.061036 m
Page 17
Numerical Methods for Mechanics II — Assignment Report
x(0.01 s) from RK4 0.039817 m
Natural freq. displacer, ω_d √(10000/1.5) = 81.65 rad/s
Natural freq. piston, ω_p √(34832/6) = 76.13 rad/s
Page 18
Numerical Methods for Mechanics II — Assignment Report
QUESTION 7: Rolling Sheet Heat Transfer
7.1 Problem Statement
Heat transfer in a rolling sheet is governed by a nonlinear ODE that includes variable conductivity,
convection, radiation, internal heat generation, and advection terms. The problem is solved using the
Galerkin Finite Element Method and the Finite Volume Method over a domain of length L = 0.1 m
discretized into 4 linear elements (5 nodes), each of length l_e = 0.025 m.
7.2 Given Parameters
Table 7.1: Rolling sheet heat transfer parameters
Parameter Value Unit
Base temperature, T_b 500 K
Ambient temperature, T_a 300 K
Reference conductivity, k_0 45 W/mK
Nonlinear conductivity param., β 0.4 —
Internal heat generation, q 10⁶ W/m³
Convection parameter, hP/(k_0δ) 55.5556 —
Advection parameter, ρC_pu/(k_0δ) 1.0833 × 10² —
Radiation linearization, 4T_a³ 1.08 × 10⁸ K³
Sheet length, L 0.1 m
Sheet thickness, δ 0.008 m
Element length, l_e 0.025 m
Boundary load coeff., h/k_0 0.2222 —
7.3 Galerkin FEM Implementation
The weak form (Galerkin statement) of the governing equation was formulated, yielding element
matrices for diffusion, advection, and reaction terms. For each linear two-node element, the matrices are:
Diffusion: k_d = (A_e/l_e)[1, −1; −1, 1], where A_e = 1 + 0.4(T_avg − 300)
Advection: k_a = (C/2)[1, −1; 1, −1], where C = 1.0833 × 10²
Reaction: k_r = (−S · l_e/6)[2, 1; 1, 2], where S = 999.4141
The source vector is f = (Q_0 · l_e/2)[1; 1] = [362,278.9; 362,278.9]
The iteration begins with an initial guess of T = 500 K for all nodes. At each iteration, A_e is updated
based on the element average temperature, and the global 5 × 5 system is reassembled and solved.
Convergence was achieved after 3 iterations.
7.4 Iteration Results
Table 7.2: FEM iteration results showing convergence
Iter. A_e T₁ (K) T₂ (K) T₃ (K) T₄ (K)
0 81 500.000 500.000 500.000 500.000
Page 19
Numerical Methods for Mechanics II — Assignment Report
1 81 699.818 699.636 699.454 699.272
2 Upd. 699.818 699.636 699.454 699.272
3 Conv. 699.818 699.636 699.454 699.272
Figure 7.1: Temperature distribution in the rolling sheet (Galerkin FEM)
7.5 Summary
Table 7.3: Key verified results for Question 7
Quantity Value
Element length, l_e 0.025 m (4 elements)
Convection parameter, hP/(k_0δ) 55.5556
Source term, Q_0 · l_e/2 362,278.909
A_e at iteration 1 81
T₁ (x = 0.025 m) 699.818 K
T₂ (x = 0.050 m) 699.636 K
T₃ (x = 0.075 m) 699.454 K
T₄ = T(L) (x = 0.100 m) 699.272 K
Convergence Achieved after 3 iterations
Boundary coefficient, h/k_0 0.2222
The Galerkin FEM solution converges rapidly due to the well-conditioned nature of the assembled
global stiffness matrix. The temperature distribution shows a sharp rise from the base temperature of 500
K to approximately 700 K, driven primarily by the large internal heat generation term, and then gradually
decreases toward the tip.
Page 20