sci.math.numanalysis
SUNA15, Nonuniform 2D convection
==================================
In the sequel, references [*] correspond with SUNA* articles.
Skew Upwind Finite Element Difference schemes were developed for a _uniform_
velocity field [05,06,07]. This obviously imposes severe limitations upon the
applicability of the theory. The present article overcomes this limitations.
We will restrict ourselves in what follows, however, to incompressible fluids.
According to [05], equation (5), the general SUFED scheme for 2D pure steady
uniform convection is (apart from a different offset for numbering):
[ + (y3  y1).u  (x3  x1).v ].(T1  T2) +
[  (y2  y1).u + (x2  x1).v ].(T1  T3) = 0
According to [14] Finite Differences schemes for the Stream Function F are:
F2  F1 = (y2  y1).u  (x2  x1).v 3
F3  F1 = (y3  y1).u  (x3  x1).v o
/ \
Consequently: / \
1 o  o 2
(F3  F1).(T1  T2)  (F2  F1).(T1  T3) = 0
This expression can be written in several other ways:
= F1.(T2  T3) + F2.(T3  T1) + F3.(T1  T2)
= T1.(F3  F2) + T2.(F1  F3) + T3.(F2  F1)
Note that the temperature at a node is associated with streamfunction values
that are _opposite_ to this node. Which seems a sensible result, because the
streamfunction has something to do with flow that comes in at a triangle side,
and such a flow can only influence the opposite node.
The above scheme was found independently, and earlier, in a more specialized
form, by my excollegue Jan Tulp. We worked at Neratoom at that time, and he
called it a "dupwind" scheme. (I don't know of any references. Neratoom does
not exist anymore.)
The form [ (F3  F1).(T1  T2)  (F2  F1).(T1  T3) ] has the following
analytical equivalent:
 dF/dk.dT/dh + dF/dh.dT/dk = 0 > (dT/dh)/(dT/dk) = (dF/dh)/(dF/dk)
With other words: the gradients of the isotherms are the same as the gradients
of the streamlines, or the isotherms follow the streamlines. Which is of course
true for pure steady convection.
Another analytical analogue follows with the continuity equation (conservation
of mass) and Green's theorem:
u.dT/dx + v.dT/dy = d(u.T)/dx + d(v.T)/dy . Integrating:
// / /
 [d(T.u)/dx + d(T.v)/dy] dx.dy = O T.u.dy  T.v.dx = O T.dF
// / /
which resembles the form: T1.(F3  F2) + T2.(F1  F3) + T3.(F2  F1)
Concerning the other form: F1.(T2  T3) + F2.(T3  T1) + F3.(T1  T2)
reverse engineering is applied:
/ / //
O F.dT = O F.[dT/dx.dx + dT/dy.dy] =  [d/dx(F.dT/dy)  d/dy(F.dT/dx)] dx.dy
/ / //
// //
=  [dF/dx.dT/dy  dF/dy.dT/dx] dx.dy =   [ u.dT/dx + v.dT/dy ] dx.dy
// //
/ / //
Conclusion: O T.dF =  O F.dT =  [ u.dT/dx + v.dT/dy ] dx.dy , analytically.
/ / //
The scheme will be applied now to a patch of triangles around a mesh point (0),
according to the figure below:
5 o  o 4 Suppose for a moment there is NO upwinding. Then:
/ \ / \
6 o  0  o 3 T0.(F2  F1) + T1.(F0  F2) + T2.(F1  F0) +
\ / \ / T0.(F3  F2) + T2.(F0  F3) + T3.(F2  F0) +
o  o T0.(F4  F3) + T3.(F0  F4) + T4.(F3  F0) + ... +
1 2 T0.(F1  F6) + T6.(F0  F1) + T1.(F6  F0) = 0
Simplified: T0.( 0 ) + T1.(F6  F2) + T2.(F1  F3) + T3.(F2  F4) + .... = 0
 NO upwinding always results in a zero coefficient for the temperature ( T0 )
 to be calculated. Which is commonly recognized as weird.
Rules for upwinding the nonuniform 2D scheme will be derived in the sequel.
For convenience, we first introduce FLUXES instead of streamfunction values.
The relationship between the former and the latter is established by means of:
3 dF1 = F3  F2
o dF2 = F1  F3 ; dF1 + dF2 + dF3 = 0
/ \ dF3 = F2  F1 or  dF2  dF3 = dF1
2 * * 1
/ \ ! Take a good look at
o  *  o this way of numbering.
1 3 2 At this point the velocity field becomes nonuniform !
According to [14], the streamfunction is discretized along a meshline (1,2) as:
F2  F1 = (y2  y1).u3  (x2  x1).v3 & same for (3,2) and (1,3) .
^ ^
It is emphasized that velocities are located at staggered, midside positions.
Furthermore, it is advantageous to establish any numbering conforming the idea
that temperatures Tk are influenced by opposite (uk,vk), according to [05], or
the idea that temperatures and fluxes are kind of "conjugate" quantities.
Upwinding is introduced next. It is remembered from above that:
 dF2.(T1  T2)  dF3.(T1  T3) = 0  ______________
  
At the other hand, equation (3) in [05] reads: ==>  dF2 >= 0 
  dF3 >= 0 
(U <= 0 ).(T1  T2) + (V <= 0 ).(T1  T3) = 0  ______________
Upwinding
Thus: T1.dF1 + T2.dF2 + T3.dF3 = 0
^^ <0 >=0 >=0 : rule of coefficients [Patankar] OK
Theorem: there is only ONE or NO upwind point in a triangle.
Lemma: there can not be two or three upwind points.
Proof: suppose there are two upwind points, say (1) and (2). Then:
dF2 >= 0 & dF3 >= 0 for (1) > F1 >= F3 & F2 >= F1
dF1 >= 0 & dF3 >= 0 for (2) > F3 >= F2 & F2 >= F1
Conclusion: F1 = F2 = F3 < F1 >= F3 >= F2 & F1 <= F2
This means that 0.T1 + 0.T2 + 0.T3 = 0 so temperatures can be anything.
The above is contradictory to the assumption that one of the <= must be < .
Beware of the special cases! There are essentially _two_ of them:
First case: F1 = F2 , F2 = F3 therefore F3 = F1 .
2______ 1 For example a CORNER in a heat exchanger. Then F1 = F2 , F2 = F3 ,
 < / therefore 0.T1 + 0.T2 + 0.T3 = 0 > temperatures _undefined_.
 / Because of other triangles in the mesh, this in fact only affects
 / the corner temperature: T2 = undefined , while T3 = T1 . A more
3 severe situation occurs when a vortex is present. Then equations
like T1 = T2 = T3 = T4 = T5 = ... = T1 occur, which render many
temperatures _independent_ of any boundary conditions!
Second case: dF1 = 0 , dF2 > 0 , dF3 < 0 (: apart from 1,2,3 permutations)
For example if a meshline coincides with a streamline: F2  F3 .
Then dF1.T1 + dF2.T2 + dF3.T3 = 0 & dF2 + dF3 = 0 > dF3.(T3  T2) = 0 .
Consequently: T2 = T3 , as it should be. This case especially occurs with flow
BOUNDARIES. It may be considered as good practice to set up separate equations
for the streamline(s) along a fixed boundary, because numerical inaccuracy with
an upwinding procedure for the bulkflow might too often "miss" points there.

* 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