Finite elements/Heat equation time integration

Solving the Matrix Problem

The matrix equations for the Poisson problem do not involve time and can be solved directly using either direct or iterative methods for solving systems of equations. You have done that in your introductory course on finite elements.

For the time-dependent heat equation, a few extra steps are needed. This is because the equations we have developed so far still have continuous time derivatives which need to be approximated.

Recall equations (45)

${\text{(45)}}\qquad \mathbf {M} ~{\dot {\mathbf {T} }}+\mathbf {K} ~\mathbf {T} =\mathbf {f} ~.$

These equations are a coupled system of first-order ordinary differential equations rather than a system of algebraic equations.

One way of solving the system of differential equations (45) is to use the generalized midpoint rule.

Generalized midpoint rule

Consider the following initial value problem

${\begin{aligned}{\dot {T}}(t)&=g(T(t))\qquad {\text{for all}}\qquad t\in [0,\tau ]\\T(0)&=T_{n}~,\end{aligned}}$

where $g:\mathbb {R} ^{}\rightarrow \mathbb {R} ^{}$ is a smooth function.

The generalized midpoint rule can be used to solve this differential equation in an approximate manner.

Let us discretize the domain $[0,\tau ]$ into increments of $\Delta t$. Let $T_{n+1}:=T(t_{n+1})$ be the approximation to the exact value of $T(t)$ at time $t_{n+1}=t_{n}+\Delta t$.

The integration rule for the generalized midpoint rule is

${\begin{aligned}T_{n+1}&=T_{n}+\Delta t~g(T_{n+\theta })\\g(T_{n+\theta })&:=(1-\theta )~g(T_{n})+\theta ~g(T_{n+1})\qquad \theta \in [0,1]~.\end{aligned}}$

When we choose $\theta =0$, we get

${T_{n+1}=T_{n}+\Delta t~g(T_{n});\qquad \theta =0\implies {\text{Forward Euler - fully explicit.}}}$

When we choose $\theta =1$, we get

${T_{n+1}=T_{n}+\Delta t~g(T_{n+1});\qquad \theta =1\implies {\text{Backward Euler - fully implicit.}}}$

When we choose $\theta =1/2$, we get

${T_{n+1}=T_{n}+\Delta t~\left({\cfrac {g(T_{n})+g(T_{n+1})}{2}}\right);\qquad \theta ={\frac {1}{2}}\implies {\text{Midpoint rule, Crank-Nicolson.}}}$

Generalized midpoint rule for Heat equation

If we apply the generalized midpoint rule to the system of differential equations

$\mathbf {M} ~{\dot {\mathbf {T} }}+\mathbf {K} ~\mathbf {T} =\mathbf {f}$

we get

${\text{(46)}}\qquad {\begin{aligned}\mathbf {M} ~{\dot {\mathbf {T} }}_{n+1}+\mathbf {K} ~\mathbf {T} _{n+1}&=\mathbf {f} _{n+1}\\\mathbf {T} _{n+1}&=\mathbf {T} _{n}+\Delta t~{\dot {\mathbf {T} }}_{n+\theta }\\{\dot {\mathbf {T} }}_{n+\theta }&=(1-\theta )~{\dot {\mathbf {T} }}_{n}+\theta ~{\dot {\mathbf {T} }}_{n+1}\end{aligned}}$

The computational problem is to find $\mathbf {T} _{n+1}$ and ${\dot {\mathbf {T} }}_{n+1}$ given $\mathbf {T} _{n}$ and ${\dot {\mathbf {T} }}_{n}$.

We start at $t=0$ at which time the initial condition $\mathbf {T} _{0}$ is given. In that case, the value of ${\dot {\mathbf {T} }}_{0}$ can be calculated using

${\dot {\mathbf {T} }}_{0}=\mathbf {M} ^{-1}(\mathbf {f} _{0}-\mathbf {K} ~\mathbf {T} _{0})~.$

Implementation of the algorithm for subsequent times may take various forms.

Let us combine the second and third equations in (46) to get

${\text{(47)}}\qquad \mathbf {T} _{n+1}=\mathbf {T} _{n}+(1-\theta )~\Delta t~{\dot {\mathbf {T} }}_{n}+\theta ~\Delta t~{\dot {\mathbf {T} }}_{n+1}~.$

If we collect the known quantities at time $t_{n}$, we can write equation (47) as

${\text{(48)}}\qquad {\mathbf {T} _{n+1}={\tilde {\mathbf {T} }}_{n+1}+\theta ~\Delta t~{\dot {\mathbf {T} }}_{n+1}}$

where

${\tilde {\mathbf {T} }}_{n+1}:=\mathbf {T} _{n}+(1-\theta )~\Delta t~{\dot {\mathbf {T} }}_{n}~.$

The quantity ${\tilde {\mathbf {T} }}_{n+1}$ is often called the predictor value of $\mathbf {T} _{n+1}$.

At this stage we can proceed in two ways (or more).

The v-form of the Generalized Midpoint Rule.

Substitute equation (48) into the first equation in (46) to get

$\mathbf {M} ~{\dot {\mathbf {T} }}_{n+1}+\mathbf {K} ~\left({\tilde {\mathbf {T} }}_{n+1}+\theta ~\Delta t~{\dot {\mathbf {T} }}_{n+1}\right)=\mathbf {f} _{n+1}~.$

Collect terms containing ${\dot {\mathbf {T} }}_{n+1}$ and rearrange to get

${\text{(49)}}\qquad {\left(\mathbf {M} +\theta ~\Delta t~\mathbf {K} \right)~{\dot {\mathbf {T} }}_{n+1}=\mathbf {f} _{n+1}-\mathbf {K} ~{\tilde {\mathbf {T} }}_{n+1}~.}$

We can solve equation (49) for ${\dot {\mathbf {T} }}_{n+1}$. Substitute this solution into equation (48) to get $\mathbf {T} _{n+1}$.

This approach is called the $\mathbf {v}$-form because the "velocity" or rate of change of the unknown quantity is calculated first followed by the actual quantity.

The d-form of the Generalized Midpoint Rule.

Substitute equation (48) into the first equation in (46) to get

$\mathbf {M} ~\left({\cfrac {\mathbf {T} _{n+1}-{\tilde {\mathbf {T} }}_{n+1}}{\theta ~\Delta t}}\right)+\mathbf {K} ~\mathbf {T} _{n+1}=\mathbf {f} _{n+1}~.$

Collect terms containing $\mathbf {T} _{n+1}$ and rearrange to get

${\text{(50)}}\qquad {\left(\mathbf {M} +\theta ~\Delta t~\mathbf {K} \right)\mathbf {T} _{n+1}=\theta ~\Delta t~\mathbf {f} _{n+1}+\mathbf {M} ~{\tilde {\mathbf {T} }}_{n+1}~.}$

The system of equations (50) can be solved for $\mathbf {T} _{n+1}$. Once we know $\mathbf {T} _{n+1}$, we can calculate ${\dot {\mathbf {T} }}_{n+1}$ using equation (48).

This approach is called the $\mathbf {d}$-form because the "displacement" or the actual unknown quantity is calculated first followed by its rate.

Explicit foward Euler method

FE System of ODEs:

$\mathbf {M} ~{\dot {\mathbf {T} }}+\mathbf {K} ~\mathbf {T} =\mathbf {f}$

We know:

$\Delta t~;\qquad \mathbf {T} _{n}~;\qquad {\dot {\mathbf {T} }}_{n}~;\qquad \mathbf {f} _{n+1}$

Apply generalized midpoint rule ($\theta =0$):

${\begin{aligned}\mathbf {M} ~{\dot {\mathbf {T} }}_{n+1}+\mathbf {K} ~\mathbf {T} _{n+1}&=\mathbf {f} _{n+1}\\\mathbf {T} _{n+1}&=\mathbf {T} _{n}+\Delta t~{\dot {\mathbf {T} }}_{n}\\\end{aligned}}$

Combine:

$\mathbf {M} ~{\dot {\mathbf {T} }}_{n+1}+\mathbf {K} ~\left(\mathbf {T} _{n}+\Delta t~{\dot {\mathbf {T} }}_{n}\right)=\mathbf {f} _{n+1}$

Solve for ${\dot {\mathbf {T} }}_{n+1}$.

Many other techniques are available for solving systems of time-dependent ODEs. We will look into some of them later.

Quality of Approximate Solutions

The usual engineer's approach is to stop after a solution has been obtained and assume that this solution is adequate. However, it is quite important to have some information about the quality of the approximation. Unless such information is available, the finite element solution is essentially useless because it could have little resemblance to the actual solution.

Verification is the process of determining if the numerical approximation is an accurate representation of the mathematical model. The first step in the process is to obtain a qualitative estimate of the error in the approximation. Functional analysis plays a vital role in determining these estimates of error.

The next step in the verification process is to obtain some information about whether an approximate solution converges to the exact solution as the mesh is refined. We can also determine what the rate of convergence or order of accuracy for a particular approach is. We will not get into the details of error estimation in this course except for a few specific cases.

The final step in the verification process involves comparisons of numerical results with known exact solutions and experimental results of well-characterized benchmark problems.

We also need to validate our models. Validation is the process of determining the degree to which our mathematical model represents physical reality (as far as the intended use of the model is concerned).

Later, we will discuss some aspects of verification and validation in the context of multi-physics problems.

← Solution of heat equation

Finite elements

Homework 1 →