Question 3 of 7: One-dimensional finite element heat conduction in a rod

Question 3: One-dimensional finite element heat conduction in a rod (20 marks)

Question text not reproduced: the examination questions are © Engineers and Geoscientists BC. Open the official past paper (linked at the top of this page) to read the question, then follow the worked solution below.

Given. A circular rod of two unequal-conductivity elements with an internal source, a prescribed lateral loss and end convection at the right-hand face only.

Find. The nodal temperatures $T_2$ and $T_3$, the conduction heat flux carried by element 1, and the magnitude and direction of the heat flow rate entering or leaving at node 1.

Approach. Build the two element conductance matrices, add the three separate load contributions (volumetric source, lateral loss, end convection), assemble, impose the prescribed temperature at node 1, solve for $T_2$ and $T_3$, then recover the element flux and the nodal heat flow rate and close a global energy audit.

1. Work out the three geometric measures the loads scale with. The three effects do not share a geometry: the source scales with the cross-sectional area, the lateral loss with the perimeter, and the end convection with the end-face area. $$A = \frac{\pi\phi^{2}}{4} = \frac{\pi(0.15)^{2}}{4} = 0.0176715\ \text{m}^{2},\qquad P_{\text{er}} = \pi\phi = \pi(0.15) = 0.4712389\ \text{m}$$
2. Form the element conductance matrices. For a two-node bar element carrying one-dimensional conduction, $[\mathbf{k}^{(e)}] = \dfrac{K_{xx}A}{L_e}\begin{bmatrix}1&-1\\-1&1\end{bmatrix}$, so with $L_e = 3\ \text{m}$: $$k^{(1)} = \frac{(150)(0.0176715)}{3} = 0.8835729\ \text{W/}^\circ\text{C},\qquad k^{(2)} = \frac{(400)(0.0176715)}{3} = 2.3561945\ \text{W/}^\circ\text{C}$$
3. Build the element load vectors from the source and the lateral loss. A uniform volumetric source is lumped equally to the two nodes of a linear element, and the outgoing lateral flux is treated the same way but with the opposite sign because it removes energy: $$f_Q = \frac{QAL_e}{2} = \frac{(400)(0.0176715)(3)}{2} = 10.60288\ \text{W},\qquad f_{q} = -\frac{q_{\text{out}}P_{\text{er}}L_e}{2} = -\frac{(12)(0.4712389)(3)}{2} = -8.48230\ \text{W}$$ so each node of each element receives $f_e = 10.60288 - 8.48230 = +2.12058\ \text{W}$. The rod is a net gainer of energy per unit length — the internal generation slightly exceeds what the surface sheds.
4. Add end convection to node 3 only. Because the question states convection arises only at the right-hand end, the convective terms use the END-FACE area $A$, not a lateral surface matrix: a stiffness contribution $hA$ on the $K_{33}$ diagonal and a load $hT_\infty A$ in row 3. $$hA = (25)(0.0176715) = 0.4417865\ \text{W/}^\circ\text{C},\qquad hT_\infty A = (25)(120)(0.0176715) = 53.01464\ \text{W}$$ Using the lateral surface convection matrix $\frac{hP_{\text{er}}L}{6}\begin{bmatrix}2&1\\1&2\end{bmatrix}$ here would be wrong by roughly an order of magnitude.
5. Assemble the global system. Superimposing the two elements and the convective terms, $$\begin{bmatrix}0.8835729 & -0.8835729 & 0\\ -0.8835729 & 3.2397674 & -2.3561945\\ 0 & -2.3561945 & 2.7979810\end{bmatrix}\begin{Bmatrix}T_1\\T_2\\T_3\end{Bmatrix} = \begin{Bmatrix}2.12058 + Q_1\\ 4.24115\\ 55.13495\end{Bmatrix}$$ where $Q_1$ is the unknown nodal heat flow at the prescribed-temperature node.
6. Apply $T_1 = 300\ ^\circ\text{C}$ and solve rows 2 and 3. Moving the known column to the right-hand side, $$\begin{bmatrix}3.2397674 & -2.3561945\\ -2.3561945 & 2.7979810\end{bmatrix}\begin{Bmatrix}T_2\\T_3\end{Bmatrix} = \begin{Bmatrix}4.24115 + (0.8835729)(300)\\ 55.13495\end{Bmatrix} = \begin{Bmatrix}269.31302\\ 55.13495\end{Bmatrix}$$ Solving the $2\times 2$ system, $$\boxed{\;T_2 = 251.47\ ^\circ\text{C},\qquad T_3 = 231.47\ ^\circ\text{C}\;}$$ Both lie between the prescribed 300 °C and the ambient 120 °C and decrease monotonically to the right, as the physics demands.
7. Part (b) — heat flux over element 1. Within a linear element the temperature gradient is constant, so Fourier's law gives a single flux value, $$q_x^{(1)} = -K_{xx}^{(1)}\frac{T_2-T_1}{L_1} = -(150)\frac{251.4667-300}{3} = \boxed{\,+2426.7\ \text{W/m}^{2}\,}$$ The positive sign means the flux is in the $+x$ direction, i.e. from the hot prescribed end towards the convecting end. For comparison element 2 carries $q_x^{(2)} = -(400)(231.4667-251.4667)/3 = +2666.7\ \text{W/m}^{2}$; the flux grows along the rod because the internal generation outweighs the lateral loss.
8. Part (c) — heat flow rate and direction at node 1. The nodal heat flow is recovered from the first row of the assembled system, which was set aside when the temperature was prescribed: $$Q_1 = k^{(1)}\left(T_1 - T_2\right) - f_e = (0.8835729)(300 - 251.4667) - 2.12058 = 42.883 - 2.121$$ $$\boxed{\;Q_1 = +40.76\ \text{W}\ \text{entering the rod, flowing in the }+x\ \text{direction}\;}$$ The sign convention is that a positive nodal value is heat supplied to the mesh, so the constant-temperature reservoir at the left end is a source: it must push 40.76 W into the bar to hold that face at 300 °C.
9. Close the global energy audit. A finite element answer to a conduction problem should always be checked against an overall balance, input plus generation equals lateral loss plus end convection: $$QAL = (400)(0.0176715)(6) = 42.412\ \text{W},\qquad q_{\text{out}}P_{\text{er}}L = (12)(0.4712389)(6) = 33.929\ \text{W}$$ $$hA\left(T_3-T_\infty\right) = (0.4417865)(231.4667-120) = 49.245\ \text{W}$$ Then $40.762 + 42.412 = 83.174\ \text{W}$ in, and $33.929 + 49.245 = 83.174\ \text{W}$ out. The balance closes exactly, which confirms both the nodal temperatures and the recovered nodal flow.