sci.math.numanalysis
SUNA36, Nonuniform 3D convection
==================================
In [06] general, but inverse, Skew Upwind Finite Element Difference schemes
were derived for pure steady, and _uniform_, convection in three dimensions.
Analytically: u.dT/dx + v.dT/dy + w.dT/dz = 0
Discretization: U.(T1  T0) + V.(T2  T0) + W.(T3  T0) = 0
Upwinding: U <= 0 ; V <= 0 ; W <= 0
u = (x1  x0).U + (x2  x0).V + (x3  x0).W
Transformation: v = (y1  y0).U + (y2  y0).V + (y3  y0).W
w = (z1  z0).U + (z2  z0).V + (z3  z0).W
It is necessary to determine the general inverse of the equations, in order
to arrive at a SUFED formulation for pure, steady but nonuniform convection.
For the sake of simplicity, select (x0,y0,z0) as the origin (0,0,0). Then:
 U   + y2.z3  y3.z2  x2.z3 + x3.z2 + x2.y3  x3.y2   u 
J .  V  =   y1.z3 + y3.z1 + x1.z3  x3.z1  x1.y3 + x3.y1   v 
 W   + y1.z2  y2.z1  x1.z2 + x2.z1 + x1.y2  x2.y1   w 
Here: J = x1.y2.z3 + x2.y3.z1 + x3.y1.z2  x1.y3.z2  x2.y1.z3  x3.y2.z1
We are looking for a physically acceptable interpretation of this result. First
the Jacobian determinant J is just (6 times) the _volume_ of the tetrahedron,
which is spanned by the vectors (x1,y1,z1) , (x2,y2,z2) , (x3,y3,z3) . Consider
a boundary triangle, which is spanned by the vectors (x2,y2,z2) and (x3,y3,z3).
The normal at the surface of this triangle times the area of it is equal to the
outer product of these two vectors (apart from a factor 2 ):
 + y2.z3  z2.y3 
2 x triangle area x normal = (x2,y2,z2) x (x3,y3,z3) =   x2.z3 + z2.x3 
 + x2.y3  y2.x3 
It is seen that J.U equals two times the inner product (.) of the area vector
and the velocity vector (u,v,w).
J.U = [ (x2,y2,z2) x (x3,y3,z3) ] . [ u , v , w ] = 2.dF1
Similarly: J.V = [ (x3,y3,z3) x (x1,y1,z1) ] . [ u , v , w ] = 2.dF2
J.W = [ (x1,y1,z1) x (x2,y2,z2) ] . [ u , v , w ] = 2.dF3
All three expressions indicate a flux dFk through an area which is opposite
to the vertex k with temperature Tk in the tetrahedron. An expression for
the node with temperature T0 seem to be missing, however. Well, not really,
since the missing component is recovered by J.(U + V + W) , according to:
J.U.T1 + J.V.T2 + J.W.T3  J.(U + V + W).T0 = 0
Which is easily derived using the first equation of this article. Doing so:
_ _ _ _ _ _ _ _
J.(U + V + W) = (r2 x r3 + r3 x r1 + r1 x r2) . u where denotes "vector".
_ _ _ _ _
= [ (r2  r1) x (r3  r1) ] . u =  2.dF0 .
_ _ _ _
Indeed, the vectors (r2  r1) and (r3  r1) span the triangle opposite to T0 .
Hence dF0 indicates the flux through that surface (watch out for the sign).
For _uniform_ 3D convection, we finally conclude from the above: 2 x
dF0.T0 + dF1.T1 + dF2.T2 + dF3.T3 = 0 and dF0 + dF1 + dF2 + dF3 = 0
Nonuniform velocities [ uk , vk , wk ] are introduced, instead of the hitherto
uniform fields [u,v,w]. Then, a method for 3D becomes a rather straightforward
generalization of the twodimensional case.
1. Adopt linear tetrahedra (three dimensional simplices) as the basic elements.
Determine the four fluxes dF , according to the recipe:
_ _
2.dF1 = (u1,v1,w1) . (r2 x r3) where u1,v1,w1 are midplane velocities.
^^^^^^^^
Notes: Apart from (sensible) permutations among the numbers (0,1,2,3);
dF1 is negative for flow coming in, positive for flow going out.
Calculated velocities can be at vertices, midside or midplane positions. SUFED,
or the modified "FLOTRAN" scheme, will work in all cases, being dependent only
upon fluxes, which are defined by _midplane_ velocities. The latter can always
be expressed in velocity values which are defined elsewhere. The reverse is not
true. For example, relationships between midplane and vertexvelocities are:
u1 := (u0 + u1 + u2)/3 ; v1 := (v0 + v1 + v2)/3 ; w1 := (w0 + w1 + w2)/3
Substitute these into 2.dF0 + 2.dF1 + 2.dF2 + 2.dF3 =
_ _
 1/3.(u1 + u2 + u3, v1 + v2 + v3, w1 + w2 + w3) . ( r21 x r31 ) +  No room
1/3.(u0 + u2 + u3, v0 + v2 + v3, w0 + w2 + w3) . ( r2 x r3 ) +  for the
1/3.(u0 + u1 + u3, v0 + v1 + v3, w0 + w1 + w3) . ( r3 x r1 ) +  vector
1/3.(u0 + u1 + u2, v0 + v1 + v2, w0 + w1 + w2) . ( r1 x r2 ) =  under
1/3. (u0,v0,w0) . [ (r2 x r3) + (r3 x r1) + (r1 x r2) ] + .... =  scores.
1/3. (u0,v0,w0) . [ (r2  r1) x (r3  r1) ] ... =  1/3 times  Sorry !
an expression having exactly the same form !!! ^^^^^^^^^^^
Like in the 2D case, this is consistent with the Gauss theorem:
// _ _ ///
 (u . dA) =  (du/dx + dv/dy + dw/dz) .dx.dy.dz
// midplane /// vertex velocities 1/2.1/3.J = volume tetrahedron
2. Next, a downstream point has to be determined. Several possibilities exist,
but in all cases we are finished, exept when there is ONE downstream point,
say at vertex (0):
dF1 >= 0 and dF2 >= 0 and dF3 >= 0 .
Again apart from permutations.
3. With help of the fluxes, an upstream position and an upstream temperature
has to be determined. Take (0) as the downstream point to be considered.
_
Call the upstream position r' . The isoparametric formulation for triangles is
used: _ _ _ _ _ _
r' = r1 + h.(r2  r1) + k.(r3  r1) = P
Joining P with (1,2,3), we see that the flux going through the plane (0,2,3)
also goes through the plane (P,2,3). Similarly, flux (0,3,1) = flux (P,3,1) and
flux (0,1,2) = flux (P,1,2). The 3 fluxes add up to: dF1 + dF2 + dF3 =  dF0 .
Now remember the _area_ coordinates of a triangle [02]. It becomes clear then,
that these areacoordinates must match precisely with the fluxes:
1  h  k =  dF1/dF0 ; h =  dF2/dF0 ; k =  dF3/dF0 .
_ _ _ _
Therefore: r' = (1  h  k).r1 + h.r2 + k.r3 which gives the
_ _ _ _ _
Recipe: r' =  dF1/dF0.r1  dF2/dF0.r2  dF3/dF0.r3 ; r' = (x',y',z')
Also: T' =  dF1/dF0.T1  dF2/dF0.T2  dF3/dF0.T3 : isoparametrics .
4. This is the deviant part. A streamlinevelocity (u*,v*,w*) is _defined_ by:
J/2. u* =  dF1.(x1  x0)  dF2.(x2  x0)  dF3.(x3  x0) _ _
J/2. v* =  dF1.(y1  y0)  dF2.(y2  y0)  dF3.(y3  y0) = + dF0.(r'  r0)
J/2. w* =  dF1.(z1  z0)  dF2.(z2  z0)  dF3.(z3  z0)
The factor /2 comes from the fact that terms 2.dFk have to be taken into
account. J = 6 * volume tetrahedron (0,1,2,3). The new vector (u*,v*,w*) is
directed _along_ the streamline, and _towards_ the downwind point (0). This,
in my not so humble opinion, is better than in the original FLOTRAN schemes,
where (u*,v*,w*) := (u0,v0,w0) . I wonder how that choice was motivated.
5. The length of the streamline and the mean velocity at the streamline are now
defined by:
ds = sqrt[(x'  x0)^2 + (y'  y0)^2 + (z'  z0)^2]
Us = sqrt[(u*)^2 + (v*)^2 + (w*)^2]
Consequently: J. Us/ds =  2.dF0 
^^^^^^^^^^^^^^^^^^^^
6. The Monotone Streamline Upwind scheme, _unified_ with SUFED (: Skew Upwind
Finite Element Differences) in three dimensions is then, at last, given by:
J. Us/ds. (T0  T') = 2.dF0.T0 + 2.dF1.T1 + 2.dF2.T2 + 2.dF3.T3 = 0

* 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
* Email: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 78 37 87 * Blood