sci.math.num-analysis SUNA56, 3-D Potential Flow ========================== Potential flow in three dimensions is described analytically by the equation of Laplace: [ (d/dx)^2 + (d/dy)^2 + (d/dz)^2 ] f = 0 Here x , y , z = space coordinates ; f = unknown function ; d = partial . In order to arrive at a Finite Element formulation, the equation is weighted weighted with a test function F (Galerkin method). Partial integration over a space volume then results in: /// - ||| [ dF/dx.df/dx + dF/dy.df/dy + dF/dz.df/dz ].dx.dy.dz + : bulk /// // + || F.[ df/dx.nx + df/dy.ny + df/dz.nz ].dA = 0 : boundary // _ n = (nx, ny, xz) = surface normal ; dA = surface area element // The boundary integral is also written as: || F.df/dn.dA = 0 // We could refer to Zienkiewicz "The Finite Element Method" (3th edition) chapter 3.3 . The derivation there is a bit clumsy, however (especially formula 3.22). But the following conclusion is right: a) "natural" boundary conditions with df/dn = 0 are automatically satisfied. Many ingredients for our own finite element formulation have been developed in previous SUNA posters. Especially the theorem that overlapping tetrahedra in a brick are equivalent with numerical integration at the vertices of that brick. We have seen that the latter is an integration method which makes F.E.M. fully equivalent with a Finite Volume Method (: SUNA50). Let the unit cube be the parent element of the brick, and numbering according to our well known convention: number - 1 = vertex coordinates (z,y,x) selected from {0,1} and conceived as bits. Then the tetrahedra under consideration are enumerated as follows: 1,2,3,5 , 2,4,1,6 , 3,1,4,7 , 4,3,2,8 , 5,7,6,1 , 6,5,8,2 , 7,8,5,3 , 8,6,7,4 | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | z 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 y 0 0 1 0 0 1 0 0 1 0 1 1 1 1 0 1 0 1 0 0 0 0 1 0 1 1 0 1 1 0 1 1 x 0 1 0 0 1 1 0 1 0 0 1 0 1 0 1 1 0 0 1 0 1 0 1 1 0 1 0 0 1 1 0 1 Given the number of the brick's corner (1,...8), the other vertex-numbers of a tetrahedron can be found by making them _one_ bit different from the first one. For reasons of comparison, also the classical partitions of an eight cornered brick into tetrahedra were taken into account (see Zienkiewicz' book par 6.3). There exist two alternatives for a subdivision into five tetrahedra; and there exist four alternatives for a subdivision into six. Care must be taken that in the classical F.E. case everything should be according to ideology: the volume of tetrahedra is 1/6 times a determinant, and there exist no diagonals crossed at the brick faces. The latter means that the two types of subdivisions into 5 tetrahedra should be used in a mixed fashion: one type for uneven and the other type for even numbered brick elements. This works for systematical numberings, and if the total number of bricks in every mesh direction is uneven. For the 6 tetrahedra subdivisions, a mere enumeration of all vertex quadruples would have been sufficient. Instead, this enumeration is calculated by a routine "coding", which could do the same job for more than 3 dimensions (iff such ever occurs). Last, but not least, also the classical trilinear interpolation with numerical integration at (Gauss & other) points was implemented. The SUNA case above is found back here for integration points positioned at the brick's corner nodes. The SUNA case, together with the above classical tetrahedron partitions, are implemented in the two programs, which are associated with this article, by a variable called "MODE". As follows: MODE Kind of brick element ---- --------------------- 0 known analytical solution 1 tetrahedra at the corners (SUNA) 2 composite of 5 tetrahedra (2 types) 3,4,5,6 different composites of 6 tetrahedra 7,8,9 numerical integration at Gauss points, at vertices (SUNA) & close to midpoint. The composite elements are implemented as "elem002.for", bricks with numerical integration are implemented as "elem004.for". The Fortran statement function kube(i,n) surely works also for hypercubes. The theory behind the elements is handled by several SUNA posters, but is repeated here shortly for convenience. For the case of a tetrahedron, the basic formula's are (with "d" partial): df/dh = dx/dh.df/dx + dy/dh.df/dy + dz/dh.df/dz df/dk = dx/dk.df/dx + dy/dk.df/dy + dz/dk.df/dz df/dl = dx/dl.df/dx + dy/dl.df/dy + dz/dl.df/dz here h,k,l = local element coordinates, x,y,z = global coordinates. Discretizations for d(x,y,z,f)/d(h,k,l) can be found by elementary means. In order to arrive at expressions for df/d(x,y,z) , the system of equations must be inverted, resulting in well known formulas for d(h,k,l)/d(x,y,z) . And: df/dx = dh/dx.df/dh + dk/dx.df/dk + dl/dx.df/dl df/dy = dh/dy.df/dh + dk/dy.df/dk + dl/dy.df/dl df/dz = dh/dz.df/dh + dk/dz.df/dk + dl/dz.df/dl Here df/d(h,k,l) = f2 - f1 , f3 - f1 , f4 - f1 . Thus: | df/dx | | -dh/dx-dk/dx-dl/dx dh/dx dk/dx dl/dx | | f1 | | df/dy | = | -dh/dy-dk/dy-dl/dy dh/dy dk/dy dl/dy | | f2 | | df/dz | | -dh/dz-dk/dz-dl/dz dh/dz dk/dz dl/dz | | f3 | ------------------------------------------------ | f4 | Approximation matrix Multiplying the transpose of the Approximation matrix with itself results in the 4 x 4 element matrix for a tetrahedron, describing a potential problem. For numerically integrated brick elements, the problem is splitted up into 8 sub-problems, which in fact means that an element-matrix is constructed for each of the integration points. Expressions d(x,y,z)/d(h,k,l) are specified at an integration point before the Jacobian matrix is inverted. Derivatives df/d(h,k,l) fall apart into coefficients for (f1, f2, f3, f4), resulting in a slightly different procedure for calculating the 4 x 3 Approximation matrix of trilinear bricks for the 3-D potential problem. The 9 elements for potential flow were tested for two different cases, where also the analytical solution is well known (the first one thanks to Usenet): a. Potential flow around a (unit) sphere b. Potential flow around a (unit) cylinder The analytical potential function for case (a) is given by: (r + 1/2r^2).z/r The analytical potential function for case (b) is given by: (r + 1/r).x/r Case (a) is implementented as the Fortran program "suna56a.for". Case (b) is implementented as the Fortran program "suna56b.for". The layout of the finite element mesh is calculated in place. Attention is paid to a "fancy statement function", which actually relates Finite Difference grids to Finite Element meshes. It is defined in "suna56b" as: nr(i,j,k)=NDZ*NDT*(i-1)+NDZ*(j-1)+k Meshing in case of the sphere, instead of the cylinder, is much more tricky, especially at the North and South Poles (:-). Also the accompanying statement function is somewhat more complicated. Basically, spherical coordinates are used for problem (a): x = r.cos(theta).cos(phi) y = r.cos(theta).sin(phi) with discrete values for (r, phi, theta). z = r.sin(theta) And, of course, cylinder (or polar) coordinates for (b): x = r.cos(psi) y = r.sin(psi) with discrete values for (r, psi, z). z = z Degenerate bricks are used where appropriate: at the pole region of the sphere problem. Take a look at the "suna56a.inp" AVS UCD output to see what I mean. Explicit boundary conditions are: - Potential = 0 for one of the boundaries. Due to this B.C. the numerical solution would be very boring if we used more than just _one_ element in the direction parallel to that boundary. - Velocity = 1 for one of it's components, at "infinity" (outer boundary). Due to this B.C. the potential along the outer boundary can be calculated as a function of one of the coordinates. Implicit boundary conditions are given by the fact that the normal derivative of the potential is zero, if nothing else is done. A comparison is made between the numerical solutions and the exact one. Well, don't shoot at me for that! Look at the outcomes yourself, at "suna56*.out". Especially in case of the sphere potential problem, for certain elements, the result is ultimately dissapointing. I think degeneration of bricks in certain regions of the mesh is responsible for it. As expected, the results for MODE= 1 are identical to those for MODE=8 . The SUNA element is among most reliable elements. It seems to behave well even in the degenerate case (: where some of it's tetrahedra are simply eliminated for "det=0". OK, so that's the trick ;-). A piece of postprocessing is done also, in order to complete the AVS UCD file (which is meant for a course I'm preparing) with "cell data". Three equations are set up, according to the theory of "Stream Function & Flow Potential" (14), when extended to three dimensions. All three are of the form: P2 - P1 = (x2 - x1).u + (y2 - y1).v + (z2 - z1).w In order to have a reasonable condition of the 3 equations, opposite midpoints of the brick's faces are used for evaluating the potential function values. This kind of "numerical differentiation" seems to work very well. - * 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