sci.math.research

SUNA05, 2-D Skew Upwind F.E.D.
==============================

This article is part of the "Series on Unified Numerical Approximations".


          SKEW UPWIND FINITE ELEMENT DIFFERENCES (SUFED)

      FOR PURE STEADY UNIFORM CONVECTION IN TWO DIMENSIONS


Steady pure 2-D convection of heat is governed by the following equation:

    U.dT/dx + V.dT/dy = 0                                                (1)

Here (x,y) = plane coordinates; (U,V) = flow velocity; T = temperature.


1. Standard F.D. scheme

It is assumed for the moment that the coordinate system applies to a grid which
is typically Finite Difference / Volume: equidistant and regular (and unity).
Consider a piece of this grid around a mesh point P.
Four points of the compass are denoted as usual: South, North, East, West.
Let there be a wind, blowing from the North-East. With other words, there is
convective heat transfer in the quadrant PEN. Then, according to the rule of
Upwind Differencing, the following equation holds:

    U.(Tp - Te) + V.(Tp - Tn) = 0    where   U <= 0 , V <= 0             (2)

This scheme is generally accepted as a proper discretization for its analytical
idealization: equation (1).
                                2
              N                 o
              |                 | \
              |                 |   \
       W ---- P ---- E        0 o ---- o 1
              |
              |
              S      F.D.  ==>  F.E.

However, there are no five, but only three nodal values of  T  involved in this
upwind scheme. Why employ a five point star, if three nodes are sufficient?
Switching to the Finite Element viewpoint, we can abstract herefrom a linear
triangle, which is at the same time the simplest 2-D element we can think of.
Instead of the F.D. compass-nomenclature, a F.E. local numbering (starting with
0) should then be adopted. And the upwind scheme (2) simply translates into:

    U.(T0 - T1) + V.(T0 - T2) = 0    where   U <= 0 , V <= 0             (3)


2. Standard F.E. transformation

Now a unit triangle can be transformed, according to standard F.E. techniques.
Then the unit base-vectors are transformed into:

    (1,0)  ->  (x1 - x0,y1 - y0)   ;   (0,1)  ->  (x2 - x0,y2 - y0)

This in turn induces a transformation of the velocities:

    u = (x1 - x0) .U + (x2 - x0) .V
    v = (y1 - y0) .U + (y2 - y0) .V

The inverse of this transformation is the one that's needed:

    U = + (y2 - y0)/J .u - (x2 - x0)/J .v
    U = - (y1 - y0)/J .u + (x1 - x0)/J .v                                (4)

Where:    J = (x1 - x0).(y2 - y0) - (x2 - x0).(y1 - y0)  ; assume that  J > 0 .

Substitution of (4) into (3) already yields the general 2-D SUFED scheme:

                 { + (y2 - y0).u - (x2 - x0).v } .(T0 - T1) +
                 { - (y1 - y0).u + (x1 - x0).v } .(T0 - T2) = 0
provided that:                                                           (5)
                   + (y2 - y0).u - (x2 - x0).v  <=  0
                   - (y1 - y0).u + (x1 - x0).v  <=  0

In fact, the most important part of the 2-D theory is completed herewith. Since
it is known now how to perform upwind differencing on an arbitrary curvilinear
(Finite Element) mesh, for pure steady uniform convection.


3. Skewed F.D. scheme

In order to make the meaning of the above theory clearer, consider a patch of
triangles, like the following:
                                                          n o ----- o ne
                                                            | II  / |
triangle  I :  x0 = y0 = y1 = 0  ,  x1 = x2 = y2 = 1 ;      |   /   |
triangle II :  x0 = y0 = x1 = 0  ,  y1 = x2 = y2 = 1 .      | /  I  |
                                                          p o ----- o e

Substitute this into the general 2-D SUFED scheme (5), resulting in:

triangle  I :  (u - v).(T0 - T1) + v.(T0 - T2) = u.(T0 - T1) + v.(T1 - T2) = 0
triangle II :  (v - u).(T0 - T1) + u.(T0 - T2) = v.(T0 - T1) + u.(T1 - T2) = 0

The upwind conditions  U <= 0  and  V <= 0  are transformed to  u <= 0  and
v <= 0  and for triangle  I :  |u| >= |v|  ,  for triangle II :  |u| <= |v|  .

Switching again to the F.D. viewpoint:

   u.(Tp - Te) + v.(Te - Tne) = 0   for  |u| >= |v|
   v.(Tp - Tn) + u.(Tn - Tne) = 0   for  |u| <= |v|                      (6)

   where:   u <= 0  and  v <= 0

It is remarked that for  u = v  both parts of the scheme are identical,
resulting in:
                Tp = Tne                                                 (7)

This assures a smooth transition between the two parts of the scheme.

Note that equation (6) only represents a quarter of the differencing scheme.
The full scheme can be obtained by adding six more triangles.
Moreover, it can be written as one equation (thanks to David Paterson):

   - Max(|u|,|v|).Tp + Max(-u-|v|,0).Te + Max(-v-|u|,0).Tn + Max(u-|v|,0).Tw
   + Max(v-|u|,0).Ts + Max(0,Min(-u,-v)).Tne + Max(0,Min(u,-v)).Tnw
   + Max(0,Min(u,v)).Tsw + Max(0,Min(-u,v)).Tse = 0                      (8)

It is asserted that the above Skew Upwind Finite Element Difference schemes are
all both conservative and monotone. When assembling the matrix of coefficients,
this matrix is diagonally dominant. For (8) it has a larger bandwidth than for
the non-skewed (standard U.D.) scheme, but it is expected that the "artificial
diffusion" for scheme (8) will be less.

It is implicitly assumed in the above that the velocity field does not vary
across an element, hence that it is uniform. This is one of the restrictions
that has to be removed later.

To be continued ...
-
* Han de Bruijn; Applications&Graphics | "A little bit of Physics * No
* TUD Computing Centre; P.O. Box 354   | would be NO idleness in  * Oil
* 2600 AJ  Delft; The Netherlands.     | Mathematics" (HdB).      * for
* E-mail: Han.deBruijn@RC.TUDelft.NL --| Fax: +31 15 78 37 87 ----* Blood