Previous Up Next

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 ( T + β
∂ T
∂ x
) = k
∂2 T
∂ x2
− 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 ( T + β
∂ T
∂ x
) h ] dx = ∫[ k
∂2 T
∂ x2
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 ( T + β
∂ T
∂ x
) h − k
∂ T
∂ x
∂ h
∂ 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 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:

T = T(xj) Tj(x)


and:

h = h(xi) hi(x)


where xj 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 = ∫[ − ρ C ( T + β
∂ T
∂ x
) h(xi) hi(x) − h(xi) k
∂ T
∂ x
∂ hi(x)
∂ x
− 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 ( T + β
∂ T
∂ x
) hi(x) − k
∂ T
∂ x
∂ hi(x)
∂ x
− a T hi(x) + f hi(x) ] dx = 0


The time derivative T 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.
Previous Up Next