1 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:
ρ C ( Ṫ + β | | ) =
k | | − 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:
| ∫ | [ ρ C ( Ṫ + β | | ) h ] dx =
| ∫ | [ k | | 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 = | ∫ | [ − ρ C ( Ṫ + β | | ) h
− k | | | | − 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 Tj(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 hi(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:
and:
where xj is the position of node j, and summation convention is used.
For a two-noded linear element this would be:
and:
The weighting functions are substituted in the weak form of the equilibrium equation:
R = | ∫ | [ − ρ C ( Ṫ + β | | ) h(xi) hi(x) −
h(xi) k | | | | −
a T h(xi) hi(x) + f h(xi) hi(x) ] dx = 0
|
Now demanding that this is true for arbitrary nodal weighting function values h(xi) gives
for each i:
R = | ∫ | [ − ρ C ( Ṫ + β | | ) hi(x) −
k | | | | −
a T hi(x) + f hi(x) ] dx = 0
|
The time derivative Ṫ will be approximated with
the Euler backward scheme 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(xj).
Thus the last expression gives a set of equations for the nodal temperatures T(xj).
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 |
ρ | dens |
C | condif_capacity |
β | 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 |
hi(x) | h[inol] |
∂ hi(x)/∂ 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_rhside
which contains the residue R for the nodes.
-
The term − ρ C Ṫ hi(x) is implemented in general.cc.
The time derivative Ṫ is evaluated lumped.
The implementation of this term is marked by
the comment // inertia (lumped).
- The term − ρ C β ∂ T/∂ x hi(x) is
implemented in general.cc.
The implementation of this term is marked
by the comment // convection with respect to mesh.
- The term − k ∂ T/∂ x ∂ hi(x)/∂ 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 hi(x) is implemented in condif.cc.
The implementation of this term is marked
by the comment // absorption.
- The term + f hi(x) is implemented in elem.cc.
The implementation of this term is marked
by the comment // element force.