sci.math.num-analysis SUNA36, Non-uniform 3-D 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 non-uniform 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_ 3-D convection, we finally conclude from the above: 2 x dF0.T0 + dF1.T1 + dF2.T2 + dF3.T3 = 0 and dF0 + dF1 + dF2 + dF3 = 0 Non-uniform velocities [ uk , vk , wk ] are introduced, instead of the hitherto uniform fields [u,v,w]. Then, a method for 3-D becomes a rather straightforward generalization of the two-dimensional 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 vertex-velocities 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 2-D 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 area-coordinates 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 streamline-velocity (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 * E-mail: Han.deBruijn@RC.TUDelft.NL --| Fax: +31 15 78 37 87 ----* Blood