Implementation of the convection-diffusion equation

In this section we will show how the convection-diffusion equation is implemented. We do this for the one-dimensional equation. See the users manual for the convection-diffusion (heat) equation:

$$\rho C ( \dot{T} + \beta \frac{\partial T}{\partial x} ) = k \frac{\partial^2 T}{\partial {x}^2} - a T + f$$

where we take here for simplicity a prescribed temperature $T$ at $x=0$ and $x=1$. Multiplication with a weighting function (=test function) $h$ and integration over the volume gives:

$$\int [ \rho C ( \dot{T} + \beta \frac{\partial T}{\partial ... ...t [ k \frac{\partial^2 T}{\partial {x}^2} h - a T h + f h ] dx$$

where $h$ is $0$ at $x=0$ and $x=1$ but arbitrary otherwise. Integration by parts of the conductive term, and moving everything to the right gives for the residue $R$ in the equation:

$$R = \int [ - \rho C ( \dot{T} + \beta \frac{\partial T}{\pa... ...artial x} \frac{\partial h}{\partial x} - a T h + f h ] dx = 0$$

where the boundary terms resulting from the integration by parts disappear because $h$ is 0 at $x=0$ and $x=1$. This form of the equilibrium equation is called the 'weak form'.

The weak form above is still perfectly equal to the heat equation. Now, however, we start approximating with the finite element method. The temperature field $T$ is approximated with finite element interpolation functions $T_j(x)$ where $j$ denotes the node number. In the standard galerkin approach, the same functions for the weighting field is used, so $h$ is approximated with $h_i(x)$. Remark 1: we use indices $i$ for the weighting function to prevent conflicting indices of the temperature and weighting functions. Remark 2:$h$ can now be called a virtual temperature field, since it consists of the same set of functions as the temperature field itself. So:

$$T = T(x_j) T_j(x)$$

and:

$$h = h(x_i) h_i(x)$$

where $x_j$ is the position of node $j$, and summation convention is used. For a two-noded linear element this would be:

$$T = T(0) (1-x) + T(1) x$$

and:

$$h = h(0) (1-x) + h(1) x$$

The weighting functions are substituted in the weak form of the equilibrium equation:

$$R = \int [ - \rho C ( \dot{T} + \beta \frac{\partial T}{\pa... ...}{\partial x} - a T h(x_i) h_i(x) + f h(x_i) h_i(x) ] dx = 0$$

Now demanding that this is true for arbitrary nodal weighting function values $h(x_i)$ gives for each $i$:

$$R = \int [ - \rho C ( \dot{T} + \beta \frac{\partial T}{\pa... ...partial h_i(x)}{\partial x} - a T h_i(x) + f h_i(x) ] dx = 0$$

The time derivative $\dot{T}$ will be approximated with the Euler backward scheme $\frac{T(t+dt)-T(t)}{dt}$. Integrals are approximated by numerically integrating over a set of integration points. The temperature $T$ in the last expression is understood to be interpolated by the finite element interpolation functions, and so depends on the nodal temperatures $T(x_j)$. Thus the last expression gives a set of equations for the nodal temperatures $T(x_j)$.

We now will discuss where the seperate terms in this expression are implemented in the finite element program. First we link terms from the expression for $R$ with the program source:

mathematical term	program
$R$	element_rhside
$\rho$	dens
$C$	condif_capacity
$\beta$	condif_flow
$k$	condif_conductivity
$a$	condif_absorption
$f$	element_volume_force
$T(t+dt)$	new_dof
$T(t)$	old_dof
$dt$	dtime
$i$	inol
$h_i(x)$	h[inol]
$\frac{\partial h_i(x)}{\partial x}$	new_d

Different parts of $R$ are implemented in different routines. First, it will be collected in the variable element_rhside. Later, element_rhside will be used to fill node_rhsidewhich contains the residue $R$ for the nodes.

• The term $- \rho C \dot{T} h_i(x)$ is implemented in general.cc. The time derivative $\dot{T}$ is evaluated lumped. The implementation of this term is marked by the comment // inertia (lumped).
• The term $- \rho C \beta \frac{\partial T}{\partial x} h_i(x)$ is implemented in general.cc. The implementation of this term is marked by the comment // convection with respect to mesh.
• The term $- k \frac{\partial T}{\partial x} \frac{\partial h_i(x)}{\partial x}$ is implemented in general.cc. The implementation of this term is marked by the comment // diffusive terms (with partial integration).
• The term $- a T h_i(x)$ is implemented in condif.cc. The implementation of this term is marked by the comment // absorption.
• The term $+ f h_i(x)$ is implemented in elem.cc. The implementation of this term is marked by the comment // element force.