Finite Element Method (Simplest case)

Motivation: Approximating solutions to ODE/PDE in a complex geometry

Steps:

  1. Domain Discretisation
  2. Derivation of Equations
  3. Assemble Element Equations
  4. Boundary Conditions
  5. Solution

Case:

Assume a cylinder with radius R and length L.

Changes in Temperature: \hat{\theta} = (T - T_{env})
Position of node: \eta = \frac{x}{L}
Constant: \mu^2 = L^2 \frac{2h}{kR}

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
fem_cylinder

Definition of the Boundary Conditions of the Cylinder

Step 1: Domain Discretisation

fem_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  \begin{cases}  \theta_i = a_1 + a_2 x_i\\  \theta_j = a_1 + a_2 x_j  \end{cases}

This then becomes \theta(x) = N_i (x) \theta_i + N_j (x) \theta_j,  with  \begin{cases} N_i (x) = \frac{x_j-x}{l}\\  N_j (x) = \frac{x-x_i}{l}  \end{cases}

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.

One thought on “Finite Element Method (Simplest case)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

w

Connecting to %s