sci.math.num-analysis SUNA49, 3-D Theorem Preliminaries ================================= This poster is a generalization of the 2D theory in SUNA11 to three dimensions. Consider the following partial differential equation: dQx/dx + dQy/dy +dQz/dz = 0 ("d" denotes partial differentiation) Here (Qx,Qy,Qz) can be interpreted as the heat-flow in a 3-dimensional medium, or as the velocities (u,v,w) in a 3-D incompressible flow field. According to the Finite Domain technique, a kind of Finite Difference Method, this equation is first integrated over a certain control volume: /// ||| (dQx/dx + dQy/dy + dQz/dz) dx.dy.dz = 0 /// Partial integration (Gauss' theorem) results in an equivalent formulation, which contains a surface-integral: // OO (Qx.dAx + Qy.dAy + Qz.dAz) where: (dAx,dAy,dAz) = infinitesimal area // * normal vector According to the (Galerkin) Finite Element Method, the conservation equation must first be multiplied by a weighting function F , and then integrated over a finite element volume, according to: /// ||| F.(dQx/dx + dQy/dy + dQ/dz) dx.dy.dz = 0 /// Partial integration results in an expression with surface-integrals, which all become part of the boundary conditions, and a term for the bulk material: /// - ||| (dF/dx.Qx + dF/dy.Qy + dF/dz) dx.dy.dz (: mind the minus sign) /// For the sake of generality, both the Finite Volume and the Galerkin formulation are assumed to be applicable at a curvilinear grid, consisting of brick shaped cells or elements. For the sake of simplicity, we wish to restrict ourselves to regular molecules and elements, such as the classical nine point star, and cubic bricks. The general curvilinear grid is mapped upon a regular grid by a transformation of coordinates x(h,k.l) , y(h,k,l) , z(h,k,l) in such a way that the _whole_ problem may be assumed to be regular in the (h,k,l) plane. This is the Finite Volume point of view, essentially a _global_ transformation. From a Finite Element point of view, however, such coordinate transformations are always carried out _locally_/element-wise: by isoparametrics. Restrictions of orthogonality do not apply in a Finite Element context. In either case, integrands of the Finite Volume and the Galerkin formulations must be transformed first into such regular (h,k,l) coordinates. The vector (Qx,Qy,Qz) is expressed in components at the (h,k,l) coordinate system as follows: | Qx | | dx/dh | | dx/dk | | dx/dl | | xh xk xl | | Qh | | Qy | = Qh.| dy/dh | + Qk.| dy/dk | + Ql.| dy/dl | = | yh yk yl | | Qk | | Qz | | dz/dh | | dz/dk | | dz/dl | | zh zk zl | | Ql | Here we made use of the following substitutions: xh = dx/dh ; xk = dx/dk ; xl = dx/dl yh = dy/dh ; yk = dy/dk ; yl = dy/dl zh = dz/dh ; zk = dz/dk ; zl = dz/dl Differentials of the global coordinates are transformed as follows: | dx | | dx/dh dx/dk dx/dl | | dh | (with proper interpretation of | dy | = | dy/dh dy/dk dy/dl | | dk | partial and non-partial "d's") | dz | | dz/dh dz/dk dz/dl | | dl | Last but not least, the chain rules for dT(h,k,l) are just the inverse of: | dT/dx | | dx/dh dy/dh dz/dh | -1 | dT/dh | | dT/dy | = | dx/dk dy/dk dz/dk | | dT/dk | | dT/dz | | dx/dl dy/dl dz/dl | | dT/dl | Parametrize the top and bottom surfaces of the finite volume by h and k . Then the two local base vectors at these surfaces are given by: (dx/dh, dy/dh, dz/dh) and (dx/dk, dy/dk, dz/dk) The normal at the surfaces is given by +/- the outer product of these vectors, and the components of the infinitesimal area's (dAx,dAy,dAz) are dh.dk times the normal: | dAx | | dy/dh.dz/dk - dz/dh.dy/dk | | dAy | = | dz/dh.dx/dk - dz/dk.dx/dh | .dh.dk | dAz | | dx/dh.dy/dk - dx/dk.dy/dh | The integrand of the Finite Volume formulation can now be written as follows, for one of the face surfaces: | Qx Qy Qz | | dAx | | dAy | = | dAz | | Qh Qk Ql | | dx/dh dy/dh dz/dh | | dy/dh.dz/dk - dz/dh.dy/dk | | dx/dk dy/dk dz/dk | | dz/dh.dx/dk - dz/dk.dx/dh | .dh.dk | dx/dl dy/dl dz/dl | | dx/dh.dy/dk - dx/dk.dy/dh | With MAPLE, we will calculate the inverse of the above 3 x 3 matrix: > with(linalg); > A := array([[ xh, yh, zh ], > [ xk, yk, zk ], > [ xl, yl, zl ]]); > inverse(A); [ yk zl - zk yl yh zl - zh yl yh zk - zh yk ] [ ------------- - ------------- ------------- ] [ %1 %1 %1 ] [ ] [ xk zl - zk xl xh zl - zh xl xh zk - zh xk ] [ - ------------- ------------- - ------------- ] [ %1 %1 %1 ] [ ] [ xk yl - yk xl xh yl - yh xl xh yk - yh xk ] [ ------------- - ------------- ------------- ] [ %1 %1 %1 ] %1 := xh yk zl - xh zk yl - xk yh zl + xk zh yl + xl yh zk - xl zh yk # The determinant of the original matrix is recognized in expression %1 . # It's precisely the Jacobian determinant J , of course. > quit; The last column in MAPLE's result gives us the thing we need. Because It proves that the product of the Jacobian matrix and the outproduct vector must be equal to the last column of the unity matrix times the Jacobian determinant: | Qx Qy Qz | | dAx | | Qh Qk Ql | | 0 | | dAy | = | 0 | . dh.dk = Ql.J.dh.dk | dAz | | J | Which gives a remarkebly simple result: || // || OO (Qx.dAx + Qy.dAy + Qz.dAz) = || // || // // // || || Qh.J.dk.dl + || Qk.J.dh.dl + || Ql.J.dh.dk + | surface integral || // // // | over all six faces || // // // | of the unit cube, || - || Qh.J.dk.dl - || Qk.J.dh.dl - || Ql.J.dh.dk . | parent of the brick || // // // Same kind of procedure for the Galerkin integral: { dF/dx.Qx + dF/dy.Qy + dF/dz.Qz}.dx.dy.dz = [ dF/dh dF/dk dF/dl ] | xh xk xl | -1 | xh xk xl | | Qh | | yh yk yl | . | yh yk yl | | Qk | .J.dh.dk.dl | zh zk zl | | zh zk zl | | Ql | Same kind of remarkebly simple result: || /// || - ||| (dF/dx.Qx + dF/dy.Qy + dF/dz.Qz) dx.dy.dz = || /// || /// || - ||| J. (dF/dh.Qh + dF/dk.Qk + dF/dl.Ql) dh.dk.dl || /// In order to shorten notation, we shall make everywhere the replacement: Qh := J.Qh ; Qk := J.Qk ; Ql := J.Ql - * 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