Question 3.8: A 2-in thick plate of 1060 alloy is uniformly heated to 300 ...

Figure 3.28 shows a 1D transient heat transfer element described in the
local coordinate system (LCS). It consists of node i and node j. The state variables at two nodes are the temperatures (T_{i},T_{j}) and temperature gradient (\dot{T} _{i},\dot{T} _{j}) . The elemental behavior is governed by

k\frac{∂^{2}T }{∂x^{2} }-c\frac{∂T}{∂t}+\dot{q} =0 (3.117)

where k is conductivity, ρ is the density, and c is the specific heat.
The element model is developed based on Eq. (3.117) as follows. Let the shape functions of an element be S_{i} (x)=1-x/L and S_{j} (x)=x/L, using the shape functions as the test function in the Galerkin method leads to the weak-form solution Eq. (3.117) as (Bi 2018),

[K_{m} ]{\begin{Bmatrix} T_{i} \\ T_{j} \end{Bmatrix}}^{(p+ 1)} =\left\{Q_{m} \right\} (3.118)

where p is the step index and [K_{m} ] is the modified conductive matrix as

[K_{M}]=\left(\frac{\rho cAL}{6\theta \Delta t}\begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} +\frac{kA}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\right) (3.119)

[Q_{M} ] is the modified load vector determined by the state of the precedent step (i.e., step i) as

\left\{Q_{m}\right\} =\begin{Bmatrix}Q_{i} \\Q_{j} \end{Bmatrix} +\frac{\rho cAL}{6\theta \Delta t} \begin{bmatrix}2 & 1 \\ 1 & 2 \end{bmatrix} {\begin{Bmatrix} T_{i} \\ T_{j} \end{Bmatrix}}^{(p)}+\left(\frac{1}{\theta }-1\right) \frac{\rho cAL}{6} \begin{bmatrix}2 & 1 \\ 1 & 2 \end{bmatrix} {\begin{Bmatrix} \dot{T} _{i} \\ \dot{T} _{j} \end{Bmatrix}}^{(p)} (3.120)

After the temperature at (p + 1) have been found, the temperature gradients can be derived by .

{\begin{Bmatrix} \dot{T} _{i} \\ \dot{T} _{j} \end{Bmatrix}}^{(p+1)}=\frac{1}{\theta \Delta t}\begin{Bmatrix} T^{(p+1)}_{i}  -T^{(p)}_{i} \\T^{(p+1)}_{j}  -T^{(i)}_{j} \end{Bmatrix} -\left(\frac{1}{\theta }-1\right) {\begin{Bmatrix} \dot{T} _{i} \\ \dot{T} _{j} \end{Bmatrix}}^{(p)} (3.121)

Since the plate is symmetric about its neutral axis, a 1D FEA model is developed for half of the plate thickness. Accordingly, Fig. 3.29 shows a decomposition which leads to an FEA model with 5 elements and 6 nodes. Every element has the length of 0.2-in.
Let the Fourier number be F_{0} =0.5, and the step time Δt is found from Eq. (3.107) as

\Delta t=\frac{(\Delta x)^{2}F_{0} }{\alpha } \quad \quad (3.107) \ \\ =\frac{(0.2)^{2}(0.5) }{({0.002675}/{(0.0975437\ast 0.214961)})} =0.1568(s)

If the history of the temperature change in 5 s is concerned, the total number N of the time steps is

N=\frac{(5)}{\Delta t} =32(steps)

Let the Euler parameter θ = 0.5 in Eq. (3.121) where the Crank-Nicolson difference is used in the approximation. At each time step p from 1 to N, the system model can be assembled from 5 element models expressed in Eq. (3.118) as

[K]^{(G)}\begin{Bmatrix} T_{1} \\ T_{2} \\ T_{3} \\ T_{4} \\ T_{5} \\ T_{6} \end{Bmatrix} ^{(p+1)} =\left\{Q\right\} ^{(p)} =\begin{Bmatrix} 2c_{1}T_{1}^{(p)} +c_{1}T_{2}^{(p)}+2c_{2}\dot{T} _{1}^{(p)} +c_{2}\dot{T} _{2}^{(p)} \\ c_{1}T_{1}^{(p)} +4c_{1}T_{2}^{(p)}+ c_{2}\dot{T} _{1}^{(p)}+4c_{2}\dot{T} _{2}^{(p)}+c_{1}T_{3}^{(p)}+c_{2}\dot{T} _{3}^{(p)} \\ c_{1}T_{2}^{(p)}+4c_{1}T_{3}^{(p)}+c_{2}\dot{T} _{2}^{(p)}+4c_{2}\dot{T} _{3}^{(p)}+c_{1}T_{4}^{(p)} +c_{2}\dot{T} _{4}^{(p)} \\ c_{1}T_{3}^{(p)} +4c_{1}T_{4}^{(p)}+c_{2}\dot{T} _{3}^{(p)} +4c_{2}\dot{T} _{4}^{(p)} +c_{1}T_{5}^{(p)}+c_{2}\dot{T} _{5}^{(p)} \\ c_{1}T_{4}^{(p)}+4c_{1}T_{5}^{(p)} +c_{2}\dot{T} _{4}^{(p)} +4c_{2}\dot{T} _{5}^{(p)}+c_{1}T_{6}^{(p)} +c_{2}\dot{T} _{6}^{(p)} \\ c_{1}T_{5}^{(p)} +2c_{1}T_{6}^{(p)} +c_{2}\dot{T} _{5}^{(p)} +2c_{2}\dot{T} _{6}^{(p)} \end{Bmatrix}     (3.122)

where

[K]^{(G)}=\begin{bmatrix} 0.0312 & -0.0045 & 0 & 0 & 0 & 0 \\ -0.0045 & 0.0624 & -0.0045 & 0 & 0 & 0 \\ 0 & -0.0045 & 0.0624 & -0.0045 & 0 & 0 \\0 & 0 & -0.0045 & 0.0624 &-0.0045 & 0 \\ 0 & 0 & 0 & -0.0045 & 0.0624 & -0.0045 \\ 0 & 0 & 0 & 0 & -0.0045 & 0.0312 \end{bmatrix} (3.123)

c_{1} and c_{2} are two constants as below

\left.\begin{matrix}c_{1}=\frac{\rho cAL}{6\theta \Delta t}=0.0089 \\\\\ c_{2}=\left(1-\frac{1}{\theta } \right) \frac{\rho cAL}{6}= -0.00069893 \end{matrix} \right\} (3.124)

At each step p, after \left\{T\right\} ^{(p+1)} is found from Eq. (3.118), the temperature gradients on nodes can be updated based on Eq. (3.121) as

\begin{Bmatrix} \dot{T} _{1} \\ \dot{T} _{2} \\ \dot{T} _{3} \\\dot{T} _{4} \\\dot{T} _{5}\\\dot{T} _{6} \end{Bmatrix} ^{(p+1)}=\begin{Bmatrix} \frac{1}{\theta \Delta t}\left(T^{(p+1)}_{1}-T^{(p)}_{1} \right)-\left(\frac{1}{\theta }-1 \right)\dot{T} ^{(p)}_{1} \\ \frac{1}{\theta \Delta t}\left(T^{(p+1)}_{2}-T^{(p)}_{2} \right)-\left(\frac{1}{\theta }-1 \right)\dot{T} ^{(p)}_{2} \\ \frac{1}{\theta \Delta t}\left(T^{(p+1)}_{3}-T^{(p)}_{3} \right)-\left(\frac{1}{\theta }-1 \right)\dot{T} ^{(p)}_{3}\\ \frac{1}{\theta \Delta t}\left(T^{(p+1)}_{4}-T^{(p)}_{4} \right)-\left(\frac{1}{\theta }-1 \right)\dot{T} ^{(p)}_{4} \\ \frac{1}{\theta \Delta t}\left(T^{(p+1)}_{5}-T^{(p)}_{5} \right)-\left(\frac{1}{\theta }-1 \right)\dot{T} ^{(p)}_{5} \\ \frac{1}{\theta \Delta t}\left(25-T^{(p+1)}_{6} \right)-\left(\frac{1}{\theta }-1 \right)\dot{T} ^{(p)}_{6} \end{Bmatrix} (3.125)

Note that the calculation for the temperature gradient of the boundary node 6 is different since the temperature at this node is fixed as 25 °C. Accordingly, when the calculation result at step p is moved to next step p + 1, the temperature at node 6 must be set back to 25 °C as the  boundary constraint. The temperature changes on nodes with respect to time are obtained in Table. 3.9.

The plots of temperature on node 1 (i.e., the center of plate) are illustrated in Figs. 3.30 and 3.31 respectively.

Table 3.9 The History of nodal temperatures in 5 s