Case:

Assume a cylinder with radius $ R $ and length $ L $.

where $ h $ is the heat transfer coefficient and $ k $ is the thermal conductivity of the cylinder.

Governing equation of the system: $ \frac{d^{2} \hat{\theta}}{d \eta^2} – \mu^2 \hat{\theta} = 0 $

with Boundary conditions:

• $ \frac{d \hat{\theta}}{d \eta} = 0 $ at $ \eta = 0 $ at the insulated end
• $ \hat{\theta}_{tip} = (T_{tip}-T_{env}) $ at $ \eta = 1$

Definition of the Boundary Conditions of the Cylinder

Step 1: Domain Discretisation

We break down the domain into (5) smaller sub-sections.

Step 2: Derivation of Equations

• Shape functions (or interpolation functions or basis functions)

$ \theta(x) = a_1 + a_2 x $  giving $ \theta_i = a_1 + a_2 x_i $ and $\theta_j = a_1 + a_2 x_j $.

This then becomes $ \theta(x) = N_i (x) \theta_i + N_j (x) \theta_j $,  with  $ N_i (x) = \frac{x_j-x}{l}$ and $N_j (x) = \frac{x-x_i}{l}$.

Step 3: Assemble Element Equations

$ \frac{d\theta}{dx} = \frac{-1}{l} \theta_i + \frac{1}{l} \theta_j $

To minimise the errors between our estimate solution $ \theta $ within the element and the exact solution, Galerkin method is applied:

$ \int_{rod} N_k (\frac{d^2 \theta}{d\eta^2} – \mu^2 \theta)d\eta = 0 $

Here we have 2 shape functions $ N_i $ and $ N_j $ per element.

We then define the position along the rod as $ \eta_e = \frac{l}{L} $, and have:

$ \frac{d\theta}{d\eta} = L\frac{d\theta}{dx} = \frac{-1}{\eta_e} \theta_i + \frac{1}{\eta_e} \theta_j $

$ N_i (\eta) = 1-\frac{\eta}{\eta_e} $

$ N_j (\eta) = \frac{\eta}{\eta_e} $

Integrating $ \int_{rod} N_k (\frac{d^2 \theta}{d\eta^2} – \mu^2 \theta)d\eta = 0 $:

$ 0=\frac{1}{\eta_e}(\theta_i – \theta_j) + \frac{d\theta}{d\eta}+\frac{\mu^2 \eta_e}{6}(2\theta_i+\theta_j) $

$ 0=\frac{1}{\eta_e}(-\theta_i+\theta_j)-\frac{d\theta}{d\eta}+\frac{\mu^2 \eta_e}{6} (\theta_i+2\theta_j) $

Turning these into matrix representation:

$ \frac{1}{\eta_e}  \begin{pmatrix} 1 & -1 \\ -1 & 1 \end{pmatrix}  \begin{Bmatrix} \theta_i \\ \theta_j \end{Bmatrix} + \frac{\mu^2 \eta_e}{6} \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \begin{Bmatrix} \theta_i \\ \theta_j \end{Bmatrix} + \begin{Bmatrix} \frac{d\theta}{d\eta}\\ -\frac{d\theta}{d\eta} \end{Bmatrix} $

Let $ \mu^2 = 3 $ and $ \eta_e = 0.2 $:

$ \begin{pmatrix} 5.2 & -4.9 \\ -4.9 & 5.2 \end{pmatrix} \begin{Bmatrix} \theta_i \\ \theta_j \end{Bmatrix} + \begin{Bmatrix} \frac{d\theta}{d\eta}\\ -\frac{d\theta}{d\eta} \end{Bmatrix} $

Likewise we have the element characteristics of all other elements. A matrix is assembled to monitor the element charactertics of the elements of interest.

For instance if we look at the first two elements (between nodes 0, 1 and 2):

$ \begin{pmatrix} 5.2 & -4.9 & 0\\ -4.9 & 5.2+5.2 & -4.9\\ 0 & -4.9 & 5.2 \end{pmatrix} \begin{Bmatrix} \theta_0 \\ \theta_1\\ \theta_2 \end{Bmatrix} = \begin{Bmatrix} \frac{-d\theta}{d\eta}\\ \frac{d\theta}{d\eta}-\frac{d\theta}{d\eta} \\ \frac{d\theta}{d\eta} \end{Bmatrix} = \begin{Bmatrix} \frac{-d\theta}{d\eta}\\ 0 \\ \frac{d\theta}{d\eta} \end{Bmatrix} $

Or, if we look at the entire domain:

$ \begin{pmatrix} 5.2 & -4.9 & 0 & 0 & 0 & 0\\ -4.9 & 5.2+5.2 & -4.9 & 0 & 0 & 0\\ 0 & -4.9 & 5.2+5.2 & -4.9 & 0 & 0\\ 0 & 0 & -4.9 & 5.2+5.2 & 4.9 & 0\\ 0 & 0 & 0 & -4.9 & 5.2+5.2 & -4.9\\ 0 & 0 & 0 & 0 & -4.9 & 5.2 \end{pmatrix} \begin{Bmatrix} \theta_0 \\ \theta_1\\ \theta_2\\ \theta_3\\ \theta_4\\ \theta_5 \end{Bmatrix} = \begin{Bmatrix} \frac{-d\theta}{d\eta}\\ \frac{d\theta}{d\eta}-\frac{d\theta}{d\eta}\\ \frac{d\theta}{d\eta}-\frac{d\theta}{d\eta}\\ \frac{d\theta}{d\eta}-\frac{d\theta}{d\eta}\\ \frac{d\theta}{d\eta}-\frac{d\theta}{d\eta}\\ \frac{d\theta}{d\eta} \end{Bmatrix} = \begin{Bmatrix} \frac{-d\theta}{d\eta}\\ 0 \\ 0 \\ 0 \\ 0\\  \frac{d\theta}{d\eta} \end{Bmatrix} $

But – we have not considered the boundary conditions.

Step 4: Boundary Conditions

Earlier in the problem definition, we learnt that:

• $\frac{d \hat{\theta}}{d \eta} = 0$ at $ \eta = 0$ at the insulated end
• $ \hat{\theta}_{tip} = (T_{tip}-T_{env})$ at $ \eta = 1$

So if we substitute in these values into our element matrix, we get:

$ \begin{pmatrix} 5.2 & -4.9 & 0 & 0 & 0 & 0\\ -4.9 & 5.2+5.2 & -4.9 & 0 & 0 & 0\\ 0 & -4.9 & 5.2+5.2 & -4.9 & 0 & 0\\ 0 & 0 & -4.9 & 5.2+5.2 & 4.9 & 0\\ 0 & 0 & 0 & -4.9 & 5.2+5.2 & -4.9\\ 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix} \begin{Bmatrix} \theta_0 \\ \theta_1\\ \theta_2\\ \theta_3\\ \theta_4\\ \theta_5 \end{Bmatrix} = \begin{Bmatrix} 0\\ 0 \\ 0 \\ 0 \\ 0 \\ \hat{\theta}_{tip} \end{Bmatrix} $ 

Step 5: Solution

After formulating the system of simultaneous equations into a matrix, we simply invert the matrix to get the solution to the problem, ie. the temperature at each of the nodes.

Closing:

This is a very simple overview of the finite element method. However, in complex geometries where elements are not define as “uniformly” as we have seen in this case, several other mathematical techniques such as Interpolation methods and Gaussian Integration. This will be covered in the more advanced article.

Adapted from: Joan J. Cerda’s Simulationsmethoden 1WS 09 10 Tutorial.