finite element analysis title slide

Adapted from: Wikimedia / CC BY-SA 4.0

UHI logo

ESF logo

SFC logo

Introduction to the Finite Element Analysis of the Energy Balance Equation

In the previous lessons, we established the Energy Balance Equation eq. (16) that governs heat transfer through a solid. This equation is termed the ‘Governing Equation’. This governing equation is a set of partial differential equations, which is difficult to solve mathematically.

The basic principle of Finite Element Analysis (FEA) is to consider a solid object through which this governing equation is applied as an assembly of individual small regions called ‘elements’. In FEA we estimate the temperature variable acting through each element in a predefined manner – this estimation adds a degree of error to our analysis, but this error is often acceptably low.

Note that in FEA, each element is assumed to be connected to its neighbouring elements and the location of these connections are termed the ‘nodes’. The complete assembly of elements is termed the ‘mesh’. The process of representing a solid object as an assemblage of elements is known as discretisation. To assess the behaviour of heat transfer through the object, our governing equation, the Energy Balance Equation, needs to be expressed in its discretised form.

example of nodes and elements in a finite element analysis

Image adapted from Wikimedia / CC BY-SA 3.0

Discretisation by Finite Difference Method

To discretise the governing equation, we first need to re-write the Energy Balance Equation in Finite Difference Form. As a reminder, this is our Energy Balance Equation:

$k (Δ v) (\frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} + \frac{\partial^{2} T}{\partial z^{2}}) + ρ (Δ v) C_{p} \frac{\partial T}{\partial t} + (Δ v) Q = 0$

Step 1: Consideration of the Second Differential of the Temperature wrt. Displacement

The first step is to evaluate the second derivative terms (∂²T/∂x², ∂²T/∂y² and ∂²T/∂z²) using Taylor Series Expansions.

Considering only the x-dimension; we shall do one expansion for (x - Δx) and a second for (x + Δx). Note that this is called the Central Difference Method as we are evaluating the temperature at a given node based upon the temperature at its neighbouring nodes.

Here are our two Taylor Series Expansions:

$f (x - Δ x) = f (x) + (- Δ x) f' (x) + \frac{1}{2} {(- Δ x)}^{2} f'' (x) + . . .$

$f (x + Δ x) = f (x) + (+ Δ x) f' (x) + \frac{1}{2} {(+ Δ x)}^{2} f'' (x) + . . .$

The ‘+ …’ indicates higher order terms in the above expansions that shall be ignored. Once again, this implies that the resulting solution approximates the true solution but shall not be exact, i.e., our solution shall contain some degree of error. This has important implications for stability conditions, which we shall consider later.

If we add these two Taylor Series Expansions together, we get:

$f (x - Δ x) + f (x + Δ x) =$ $[f (x) + (- Δ x) f' (x) + \frac{1}{2} {(- Δ x)}^{2} f'' (x)]$ $+$ $[f (x) + (+ Δ x) f' (x) + \frac{1}{2} {(+ Δ x)}^{2} f'' (x)]$

$f (x - Δ x) + f (x + Δ x) = 2 f (x) + {(Δ x)}^{2} f'' (x)$

Note that we are seeking to establish a solution for ∂²T/∂x², which is the term f’’(x)

So, if we rearrange the function for f’’(x), we get:

$f'' (x) = \frac{f (x - Δ x) - 2 f (x) + f (x + Δ x)}{Δ x^{2}}$

If we substitute our function (f’’(x)) with our second derivative term (∂²T/∂x²), this equation can be rewritten as:

${\frac{\partial^{2} T}{\partial x^{2}} |}_{x} = \frac{T_{(x - 1)} - 2 T_{(x)} + T_{(x + 1)}}{Δ x^{2}}$

This equation is effectively a finite difference representation of the second derivative term (∂²T/∂x²) in the Energy Balance Equation. This equation states a relationship for temperature at any specific node based on the temperature at its neighbouring nodes.

Just for clarity, it should be noted that this equation provides a governing relationship applicable for all internal nodes at the same point in time. Therefore, we shall add the time (t) term into the equation to retain consistent notation.

Note that this equation does not hold true for nodes at the boundary of the object, i.e., those subject to exposure to the surroundings.

${\frac{\partial^{2} T}{\partial x^{2}} |}_{x}^{t} = \frac{T_{(x - 1)}^{t} - 2 T_{(x)}^{t} + T_{(x + 1)}^{t}}{Δ x^{2}}$ eq.(20)

example diagram

A second point to note for clarity; this equation provides a governing relationship applicable for all internal locations in the x-direction only. The same process of discretisation should be carried out for ∂²T/∂y² and ∂²T/∂z² terms, resulting in:

${\frac{\partial^{2} T}{\partial y^{2}} |}_{y}^{t} = \frac{T_{(y - 1)}^{t} - 2 T_{(y)}^{t} + T_{(y + 1)}^{t}}{Δ y^{2}}$ eq.(21)

${\frac{\partial^{2} T}{\partial z^{2}} |}_{z}^{t} = \frac{T_{(z - 1)}^{t} - 2 T_{(z)}^{t} + T_{(z + 1)}^{t}}{Δ z^{2}}$ eq.(22)

Step 2: Consideration of the First Differential of the Temperature wrt. Time

The second step is to evaluate the first derivative term (i.e.: ∂T/∂t). This is a considerably simpler process as we do not need to undertake any Taylor Series Expansions to discretise this term. If we consider an infinitesimally small time step (Δt) defining the interval between an initial time (t) and the next point in time (t+1), we can say that the rate of change of temperature is the difference in temperature divided by the time interval, i.e.:

${\frac{\partial T}{\partial t} |}_{n = x : y : z}^{t} = \frac{T_{n}^{(t + 1)} - T_{n}^{(t)}}{Δ t}$ eq.(23)

Note that this is called the Forward Difference Method as we are seeking to define the temperature at any given point based only through the forward direction of time. In this case the Forward Difference Method is used because the heat flow as a function of time is an irreversible process as governed by the 2nd Law of Thermodynamics.

It should also be noted that eq. (23) is stated in terms of n which represents any given node in a mesh, which could be express in terms of its x, y and z coordinates.

Step 3: Constructing the Discrete Form of the Governing Equation

The final step is to substitute the discrete functions derived in Step 1 and Step 2 back into the governing equation, thus:

$k (Δ V) (\frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} + \frac{\partial^{2} T}{\partial z^{2}}) + ρ (Δ V) C_{p} \frac{\partial T}{\partial t} + (Δ V) Q = 0$

becomes:

$k (Δ V) ([\frac{T_{(x - 1)}^{t} - 2 T_{(x)}^{t} + T_{(x + 1)}^{t}}{Δ x^{2}}] + [\frac{T_{(y - 1)}^{t} - 2 T_{(y)}^{t} + T_{(y + 1)}^{t}}{Δ y^{2}}] + [\frac{T_{(z - 1)}^{t} - 2 T_{(z)}^{t} + T_{(z + 1)}^{t}}{Δ z^{2}}]) + ρ (Δ V) C_{p} [\frac{T_{n}^{(t + 1)} - T_{n}^{(t)}}{Δ t}] + (Δ V) Q = 0$

We can re-arrange the above formula to isolate the future term for any internal node, aka: $T_{n}^{(t + 1)}$

$T_{n}^{(t + 1)} =$ $T_{n}^{(t)} - \frac{k}{ρ C_{p}} ((\frac{Δ t}{Δ x^{2}}) (T_{(x - 1)}^{t} - 2 T_{(x)}^{t} + T_{(x + 1)}^{t}) +$ $(\frac{Δ t}{Δ y^{2}}) (T_{(y - 1)}^{t} - 2 T_{(y)}^{t} + T_{(y + 1)}^{t}) +$ $(\frac{Δ t}{Δ z^{2}}) (T_{(z - 1)}^{t} - 2 T_{(z)}^{t} + T_{(z + 1)}^{t})) - Δ t Δ V Q$ eq.(24)

This is the fully discretised form of our governing equation.

Note: Von-Neumann Stability Analysis

As we noted in Step 1, the process of discretisation contains some degree of inherent numerical error. Von Neumann Stability Analysis considers how these errors are magnified with each step forward in time.

In a stable finite difference scheme, the errors inherent in each step do not result in magnification of errors as time progresses, i.e., errors either stay constant or naturally decay over time. In an unstable finite difference scheme, the errors are amplified on each step such that as time (t) increases the solution becomes chaotic.

Proof of the Von Neumann Stability Analysis applied to the heat equation is outside the scope of this material, however, students should note the following stability conditions, which must be adhered to in our finite difference scheme produced in Step 3. For the solution to be stable, our time step functions with regards to the square of mesh dimensions must be less than ½, i.e.:

$|- \frac{k}{ρ C_{p}} (\frac{Δ t}{Δ x^{2}})| \leq \frac{1}{2}$

$|- \frac{k}{ρ C_{p}} (\frac{Δ t}{Δ y^{2}})| \leq \frac{1}{2}$

$|- \frac{k}{ρ C_{p}} (\frac{Δ t}{Δ z^{2}})| \leq \frac{1}{2}$

These equations place some important constraints on our analysis. In particular it should be noted that if we are at the boundary of stability, i.e.: $v = \frac{1}{2}$ , and we wish to half the mesh size (reducing $Δ x$ , $Δ y$ or $Δ z$ by 50%) then we are required to decrease the time step size by a factor of 4 in order to remain stable.

Step 4: Process for Solving the Discretised Equation

The following setup conditions are required in order to solve the problem:

The initial temperature distribution at internal nodes through the object $T_{n}^{0}$ – this is termed the initial condition. Where n which represents any internal node in the mesh with x, y and z coordinates
The temperature at the boundaries of the object, e.g. $T_{0}^{t}$ – these are termed the boundary condition.
A law that governs how heat is transferred between points along the bar – this is termed the governing equation, which we derived in Steps 1 to 3.

Using eq. (24), and with known initial conditions and boundary conditions, we are able to project the temperature at any location one step forward in time, $T_{n}^{(t + 1)}$ . Undertaking this calculation for all internal nodes through the object, we can then continue to project our analysis further into the future, i.e., for $T_{n}^{(t + 2)}$ and then for $T_{n}^{(t + 3)}$ etc. This is termed a Time Marching Calculation.

Example 7

To illustrate how we may use the discretised governing equation in practice we shall approximate the heat transfer through a simplified physical problem: How does heat diffuse through a bar of length (L) over time (t)?

For simplicity we have placed one end of the bar at position x = 0 and the other end of the bar is at position x = L in the illustration above.

It shall be assumed that the bar is much longer in the x direction than it is thick (in the y or z directions) and therefore we are only concerned with the diffusion of heat as a function of x and t. By not considering diffusion in the y or z direction; this is a one-dimensional heat transfer problem. It shall also be assumed that the bar is perfectly insulated along its length such that no heat is entering or leaving the bar in the y or z direction. It shall also be assumed that no heat is generated from within the bar itself, i.e., Q = 0.

With these assumptions in mind, our mathematical model of the temperature distribution depends on the following:

The initial temperature distribution at internal nodes through the object $T_{x}^{0}$ – the initial condition.
The temperature at each end face of the bar. $T_{0}^{t}$ and $T_{L}^{t}$ – the boundary condition
A law that approximates how heat is transferred between points along the bar – the discretised governing equation.

Note that the discretised governing equation, eq. (24) stated in the x-direction only and with Q = 0 simplifies to:

$T_{n}^{(t + 1)} = T_{n}^{(t)} + v (\frac{Δ t}{Δ x^{2}}) (T_{(x - 1)}^{T} - 2 T_{(x)}^{T} + T_{(x + 1)}^{t})$

where: $v = (- \frac{k}{ρ C_{p}})$ and is termed the diffusivity of the material.

In the problem that we shall be solving each of the above requirements shall be known. This particular kind of problem is named an Initial Boundary Value Problem (IBVP).

In solving by discretisation, we are effectively reducing the problem statement as stated previously from:

How does heat diffuse through a bar of length (L) over time (t)?

How does the heat distribute across each ‘node’ through a bar of length (L) over time (t)?

In this case, instead of trying to solve the governing equation for T_x (the temperature at any point along the bar), we shall solve the governing equation for T₀, T₁, T₂ and T₃ (the temperature at each of these selected nodes).

We may then linearly interpolate between those points to estimate any other value of temperature along the bar as illustrated.

interpolated values

We may then linearly interpolate between those points to estimate any other value of temperature along the bar.

linear graphing of temperature points

In terms of terminology, the points along the bar for which we shall solve the governing equation are the ‘nodes’ and those areas in between nodes are the ‘elements’.

In our example we shall evaluate heat transfer through the bar subject to the following conditions:

The length (L) of the bar is 6 m.
When constructing the mesh, the bar is broken down into 3 elements (4 nodes).
The diffusivity of the material is 9.
The temperature at x = 0 is held constant at 100°C.
The temperature at x = L is held constant at 0°C.
The initial temperature for internal nodes is 20°C.

Solve the governing equation for the 2 internal nodes: T₁ and T₂.as we progress through time to t = 1 and then t = 2.

Example 7 – Answer

How does the heat distribute across each ‘node’ through a bar of length (L) over time (t)?

To answer this question, we shall first define the initial conditions, i.e., the temperature at each node at time t = 0.

The problem specifies initial conditions as:

The temperature at x = 0 is held constant at 100°C.

$T_{n (0)}^{t (0)} = 100 ° C$

The temperature at x = L is held constant at 0°C.

$T_{n (L)}^{t (0)} = 0 ° C$

The initial temperature for internal nodes is 20°C.

$T_{n (1)}^{t (0)} = T_{n (2)}^{t (0)} = 20 ° C$

Our physical problem therefore takes the initial form illustrated here.

solution diagram

Consideration of mesh size (Δx)

The mesh has been selected such that the bar (of length 6 m) is discretised into 4 nodes (i.e., 3 elements). The size of Δx is therefore given by:

$Δ x = (\frac{L}{# elements}) = \frac{6}{3} = 2 m$

Confirmation of Von-Neuman Stability

The Von-Neuman Stability criterion is stated as:

$|v (\frac{Δ t}{Δ x^{2}})| \leq \frac{1}{2}$

In our problem, the diffusivity (v) is stated as 0.9, therefore:

$0.9 (\frac{1}{2^{2}}) = 0.225 \leq \frac{1}{2}$

Consideration of first time step t(1)

The temperatures at each face are held constant therefore we can already infer the temperatures at the boundary nodes:

$T_{n (0)}^{t (1)} = 100 ° C$ and $T_{n (L)}^{t (1)} = 0 ° C$

The temperatures for each internal node are controlled by the governing equation and are therefore affected their current temperature and the current temperature of their neighbours.

$T_{n (1)}^{t (1)}$ is given by:

$T_{n}^{(t + 1)} = T_{n}^{(t)} + v (\frac{Δ t}{Δ x^{2}}) (T_{(x - 1)}^{t} - 2 T_{(x)}^{t} + T_{(x + 1)}^{t})$

$T_{n (1)}^{t (1)} = T_{n (1)}^{(0)} + v (\frac{Δ t}{Δ x^{2}}) (T_{n (0)}^{t (0)} - 2 T_{n (1)}^{t (0)} + T_{n (2)}^{t (0)})$

$T_{n (1)}^{t (1)} = 20 + 0.9 (\frac{1}{2^{2}}) (100 - 2 \times 20 + 20) = 38.0 ° C$

$T_{n (2)}^{t (1)}$ is given by:

$T_{n (2)}^{t (1)} = T_{n (2)}^{(0)} + v (\frac{Δ t}{Δ x^{2}}) (T_{n (1)}^{t (0)} - 2 T_{n (2)}^{t (0)} + T_{n (L)}^{t (0)})$

$T_{n (1)}^{t (1)} = 20 + 0.9 (\frac{1}{2^{2}}) (20 - 2 \times 20 + 0) = 15.5 ° C$

Consideration of first time step t(2)

Once again, the temperatures at each face are held constant therefore we can already infer the temperatures at the boundary nodes:

$T_{n (0)}^{t (2)} = 100 ° C$ and $T_{n (L)}^{t (2)} = 0 ° C$

The temperatures for each internal node are now dependent upon their updated temperature and that of their neighbours.

$T_{n (1)}^{t (2)}$ is given by:

$T_{n (1)}^{t (2)} = T_{n (1)}^{(1)} + v (\frac{Δ t}{Δ x^{2}}) (T_{n (0)}^{t (1)} - 2 T_{n (1)}^{t (1)} + T_{n (2)}^{t (1)})$

$T_{n (1)}^{t (2)} = 38.0 + 0.9 (\frac{1}{2^{2}}) (100 - 2 \times 38.0 + 15.5) = 46.9 ° C$

$T_{n (2)}^{t (2)}$ is given by:

$T_{n (2)}^{t (2)} = T_{n (1)}^{t (1)} + v (\frac{Δ t}{Δ x^{2}}) (T_{n (1)}^{t (1)} - 2 T_{n (2)}^{t (1)} + T_{n (L)}^{t (1)})$

$T_{n (1)}^{t (2)} = 15.5 + 0.9 (\frac{1}{2^{2}}) (38 - 2 \times 15.5 + 0) = 17.1 ° C$

In more complex systems (with greater mesh sizes, great time periods, complex geometry and analysis in multiple dimensions) the above analysis quickly becomes impractical. In these cases, we use dedicated mathematical software packages to undertake the analysis on our behalf. In the next lessons, we shall introduce the ANSYS software environment and learn how to use the steady-state thermal toolset to undertake heat transfer analyses of complex systems.