Feb 20, 2004 ... Moving-mesh unstructured Finite Volume Method (FVM) is a good ... vertex-
based unstructured mesh motion solver designed to work with the ...
Automatic Mesh Motion for the Unstructured Finite Volume Method
Hrvoje Jasak
a,∗
ˇ Zeljko Tukovi´c
b
a
Nabla Ltd. The Mews, Picketts Lodge, Picketts Lane, Salfords, Surrey, RH1 5RG England
b
Faculty of Mechanical Engineering and Naval Architecture, University of Zagreb, Ivana Luˇci´ca 5, 10 000 Zagreb, Croatia
Abstract Moving-mesh unstructured Finite Volume Method (FVM) is a good candidate for tackling flow simulations where the shape of the domain changes during the simulation or represents a part of the solution. For efficient and user-friendly approach to the problem, it is necessary to automatically determine the point positions in the mesh, based on the prescribed boundary motion. In this paper, we present a vertex-based unstructured mesh motion solver designed to work with the movingmesh FVM. Motion is determined by solving the Laplace equation with variable diffusion on mesh points, using a tetrahedral decomposition of polyhedral cells. Cell decomposition and discretisation guarantees that an initially valid mesh remains geometrically valid for arbitrary boundary motion. Efficiency of the method is preserved by careful discretisation and the choice of iterative solvers, allowing the motion solver to efficiently couple with the FVM flow solver. This combination is tested on two free surface tracking flow simulations, including the simulation of free-rising air bubbles in water. Key words: Moving mesh, vertex motion, motion solver, unstructured, finite volume, free surface AMS: 74S10, 65M99, 76T10, 74S05
∗ Corresponding author Email address:
[email protected] ( Hrvoje Jasak ). URL: http://www.nabla.co.uk ( Hrvoje Jasak ).
Preprint submitted to Elsevier Science
20 February 2004
1
Introduction
There exists a number of physical phenomena in which the continuum solution couples with additional equations which influence the shape of the domain or the position of an internal interface. Examples of such cases include prescribed boundary motion simulations in pumps and internal combustion engines, free surface flows, where the interface between the phases is captured by the mesh, solidification and solid-fluid interaction, where the deposition or deformation of a solid changes the shape of the fluid domain etc. Two most popular approaches are based on tracking the front of interest, either by marker particles or an indicator variable (e.g. [1–3]), or by deforming the computational mesh to accommodate the interface motion. In the deforming mesh method, the computational mesh is adjusted to the shape of the boundary which is updated in every step of the transient simulation. Motion of all points internal to the mesh is based on the prescribed boundary motion. The main difficulty in tackling cases with variable geometry is maintaining the mesh quality. Several deforming mesh algorithms have been presented in the past, with various approaches to defining mesh motion. Behr and Tezduyar [4,5] use explicit algebraic expressions in the horizontal and vertical direction with a Finite Element (FE) Arbitrary Lagrangian-Eulerian (ALE) solver to simulate free-surface flows with mesh deformation. The most popular method to date is the spring analogy [6,7]. Here, all point-to-point connections within the mesh are replaced by linear springs and point motion is obtained as a response to the boundary loading. Several Finite Volume (FV) variants exist, e.g. [6,8,9], covering mostly unstructured triangular meshes in 2- and 3-D. However, this approach proved to lack robustness, particularly for arbitrarily unstructured meshes common in FV simulations. A review of merits and limitations of spring analogy and its variants is given by Blom [7]. In an effort to improve the robustness of the method, Farhat et al. [10,11] propose the addition of torsional springs to control all mechanisms of invalidating a tetrahedral cell. Numerous other approaches to creating a robust mesh motion solver include the use of Laplacian smoothing [12–15] with constant and variable diffusivity and the pseudo-solid approach [16–22] in the ALE FEM codes. Notably, in an effort to simultaneously control the position of the free boundary and mesh spacing next to it, Helenbrook [23] proposes the use of a biharmonic equation to govern mesh motion. From the FV viewpoint, research into dynamic mesh deformation seems to be limited to various forms of spring analogy and sometimes limited to triangular/tetrahedral cells. In order to remedy the lack of robustness, spring analogy 2
is sometimes used in conjunction with re-meshing techniques, providing a rescue path when the motion algorithm fails. In this paper we will present a general-purpose moving mesh algorithm developed to simulate deforming mesh cases compatible with arbitrarily unstructured FV solvers. A new second-order polyhedral “motion element” consistent with the FV mesh handling has been developed and used with a vertex-based solution method. A crucial part of the algorithm is that its efficiency matches the segregated FV flow solvers, both is terms of storage and CPU time requirements. The deforming mesh solver will be validated in isolation and as a part of a FV free surface flow solver based on the surface tracking approach. Robustness and efficiency of the motion solver will be examined on several 2- and 3-D test cases. The objective of the study is to assemble a FV surface tracking solver capable of performing Direct Numerical Simulation (DNS) of rising gas bubbles in liquids. For this purpose, the two solvers will be closely integrated, with particular attention to their joint efficiency and data sharing. Here, the flow equations are solved using a standard FVM approach and the motion of the free surface is obtained as a part of the solution. The rest of this paper will be organised as follows. In Section 2 the FV method for arbitrary moving volumes will be summarised. We will present the requirements on the automatic mesh motion system and review the deficiencies of past efforts in this direction in Section 3. A notable part of this effort is a review of mesh handling in an unstructured FVM code, together with the typical errors in the mesh structure (both topological and geometrical). Section 4 lays the foundation for a novel automatic mesh motion method, starting from the requirements on a robust motion system, choice of motion equation, solution variable, appropriate polyhedral cell decomposition and control of mesh quality through variable diffusion in the motion equation. A crucial part of development is the support for motion of arbitrary polyhedra. The new method is tested on two sample problems in Section 5. The paper is completed with two examples of free surface flows, including a simulation of a free-rising air bubble in water in 2- and 3-D and a closed with a short summary.
2
Finite Volume Method on Moving Meshes
A “static mesh” FVM is based on the integral form of the governing (conservation) equation over a Control Volume (CV) fixed in space. More generally, the integral form of the conservation equation for a tensorial property φ defined per unit mass in an arbitrary moving volume V bounded by a closed surface 3
S states: I I Z d Z ρφ dV + ds•ρ(u − ub )φ = − ds•ρqφ + Sφ dV, dt V
S
S
(1)
V
where ρ is the density, u is the fluid velocity, ub is the boundary velocity and qφ and Sφ are the surface and volume sources/sinks of φ, respectively. As the volume V is no longer fixed in space, its motion is captured by the motion of its bounding surface S by the velocity ub . Unstructured FVM discretises the computational space by splitting it into a finite number of convex polyhedral cells bounded by convex polygons which do not overlap and completely cover the domain. The temporal dimension is split into a finite number of time-steps and the equations are solved in a timemarching manner. A sample cell around the computational point P located in its centroid, a face f , its area vector sf and the neighbouring computational point N are shown in Fig. 1. z
N s f
x y
f P
Fig. 1. Finite volume cell.
Second-order FV discretisation of Eqn. (1) transforms the surface integrals into sums of face integrals and approximates them to second order using the mid-point rule: X (ρP φP VP )n − (ρP φP VP )o X + ρf (F − Fs )φf = − sf •ρqφ + Sφ VP , (2) ∆t f f
where the subscript P represents the cell values, f the face values and superscripts n and o the “new” and “old” time level, VP is the cell volume, F = sf •uf is the fluid flux and Fs is the mesh motion flux. The fluid flux F is usually obtained as a part of the solution algorithm and satisfies the conservation requirements (if any). 4
Compared with the FVM on a static mesh (e.g. [24]), Eqn. (2) shows only two differences: the temporal derivative introduces the rate of change of the cell volume and the mesh motion flux accounts for the grid convection. The relationship between the two is governed by the space conservation law [25]: I d Z dV − ds•ub = 0. dt V
(3)
S
While Eqn. (3) is always satisfied in the integral form, it also needs to be preserved in the discrete form: VPn − VPo X Fs = 0. − ∆t f
(4)
For this reason, the mesh motion flux Fs is calculated as the volume swept by the face f in motion during the current time-step rather than from the grid velocity ub , making it consistent with the cell volume calculation.
3
Deforming Mesh Simulation
From the case setup point of view, it is the fact that the shape of the domain changes in time that influences the solution; in fact, the boundary shape itself may in some cases be a part of the solution. Thus, one can distinguish between boundary motion and internal point motion. Boundary motion can be considered as given: it is either prescribed by external factors, e.g. piston and valve motion for in-cylinder flow simulations in internal combustion engines, or a part of the solution as in free surface tracking simulations. The role of internal point motion is to accommodate changes in the domain shape (boundary motion) and preserve the validity and quality of the mesh. Internal point motion influences the solution only through mesh-induced discretisation errors [26] and is therefore detached from the remainder of the problem. Consequently, internal point motion can be specified in a number of ways, ideally without user interaction. In the past, point motion in the FVM has been provided in various ways, ranging from pre-calculated point positioning, interaction with a pre-processor or a mesh generator, to the more flexible and user-friendly automatic methods. The latter provide great advantage, both in terms of easier and faster case setup, allowing for dynamically changing topology (e.g. adaptive refinement) or automatic improvements of mesh quality. 5
In summary, the objective of automatic mesh motion is to determine internal point motion (not involving topological changes) to conform with the given boundary motion while preserving mesh validity and quality.
3.1 Mesh Definition and Validity A valid mesh is a pre-requisite for a good numerical solution and a critical ingredient of automatic mesh motion. For this reason, we will briefly summarise the validity and quality measures from the FVM standpoint. The investigation of mesh validity can be separated into topological and geometrical tests. The first group contains the tests that can be performed without knowing the actual point positions, while the second deals with the shape of cells and the boundary. Note that it is the job of mesh generation to produce a mesh satisfying these requirements; here, we shall concentrate on typical errors found in real-life meshes and on methods of preserving the validity of an initially valid computational mesh.
3.1.1 Face-Based Mesh Definition Traditional points-and-cells mesh definition consists of a list of points and a list of cells defined in term of point labels. Additionally, the vertex ordering pattern is pre-defined for every available cell shape and allows the mesh faces to be calculated. This approach limits the number of available cell shapes which, while acceptable in the FEM (due to the fact that a shape function needs to be defined a-priori for every possible cell shape), it is unnecessarily limiting for the face-addressed FVM [27]. In the face-addressed mesh definition, a polyhedral mesh for the FVM is defined by the following components: • A list of points. For every point, its space co-ordinates are given; the point label is implied from its location in the list. Every point must be used in at least one face; • A list of polygonal faces, where a face is defined as an ordered list of point labels. Faces can be separated into internal (between two cells) and boundary faces. Every face must be used by at least one cell; • A list of cells defined in terms of face labels. Note that the cell shape is unknown and irrelevant for discretisation; • Boundary faces are grouped into patches, according to the boundary condition. A patch is defined as a list of boundary face labels. Additionally, face ordering is also enforced. For internal faces it is possible to 6
define the owner and neighbour cell such that owner appears first in the cell list. Face orientation is determined using the right-hand rule and it is such that the area vector points outwards from the owner cell. The face list will first collect all internal faces and then all boundary faces patch by patch in the order of patch definition. Internal faces are ordered to contain all faces from the first cell with the increasing neighbour label, followed by the faces owned by the second cell etc. This approach has proven to be robust and easy to handle as it enforces strict and unique face ordering 1 .
3.1.2 Topological Tests Face-addressed mesh definition is exposed to some mesh topology errors which are “impossible” in the points-and-cells approach; interestingly, the errors actually occur even there but typically remain undetected. Topological validity tests consist of the following criteria: • A point can appear in a face only once; • A face can appear in a cell only once. A face cannot belong to more than two cells. A boundary face can belong to only one patch; • Two cells can share no more than one face; • Collecting all faces from one cell and decomposing faces into edges, every edge must appear in exactly two cell faces; • Collecting all faces from the boundary and decomposing faces into edges, every edge must appear in exactly two boundary faces. The first four conditions control the validity of the mesh definition while the last two check that all cells and the boundary hull are topological closed. Additionally, mesh ordering rules are checked and enforced.
3.1.3 Geometrical Tests Geometrical tests deal with the positivity of face areas and cell volumes, as well as convexity and orientation requirements. In the context of second-order FVM, it is sufficient to use the weak definition of a convex shape. Here, the geometrical measures (face area and normal vector, face and cell centroid, volume swept by the face in motion etc.) for a polygonal face are calculated by decomposing the face into triangles. Two possible triangular decompositions of a polygon (with different severity of the convexness criteria) are shown in Fig. 2. 1
An alternative definition, giving the owner and neighbour cell index for mesh faces instead of the faces of a cell is also possible. The difference between the two is only in terms of enforcing the ordering of the face list.
7
(a) Using face centroid.
(b) Using internal edges.
Fig. 2. Decomposing a face into triangles.
A face is considered convex if all triangle normals point in the same direction. For a cell, where the metrics are calculated on a tetrahedral decomposition, an equivalent convexness definition is used. The tetrahedra are constructed using the (approximate) cell centroid and the triangles of the face decomposition. Geometrical validity criteria can be summarised as follows: • All faces and cells must be weakly convex; • All cells must be geometrically closed: the sum of outward-pointing face area vectors for a cell faces must be zero to machine tolerance; • The boundary must be geometrically closed (see above); • For all internal faces, the dot-product of the face area vector sf and the df = P N , Fig. 1, must be positive; this is usually termed the orthogonality test: df •sf > 0.
(5)
While some of the tests are clearly redundant (e.g. if a cell is closed topologically it is also closed geometrically), experience shows that various combinations of failures indicate common errors in mesh generation. We shall assume the existence of a topologically and geometrically valid mesh as a starting point for automatic mesh motion.
3.1.4 Preserving Mesh Quality in Motion During mesh motion, mesh topology remains unaffected and only the point positions change. Thus, preserving the mesh quality only relates to the geometrical tests. Moreover, once the convexness and orthogonality tests are satisfied, an initially valid mesh remains valid if no triangles or tetrahedra are inverted. We will use this criterion to examine existing automatic mesh 8
motion algorithms and devise a new approach.
3.2 Failure Modes of Simple Motion Algorithms
Let us consider the limitations and failure modes of some existing mesh motion techniques as a basis for the new approach.
3.2.1 Cell-Based Motion Equation The simplest suggestion for automatic mesh motion in the FV framework would be to re-use the available numerical machinery and solve a Laplace (or a linear elastic solid) equation to provide vertex motion. However, as the FVM provides the solution in cell centres and motion is required on the points, this necessarily leads to interpolation. Experience shows that is extremely difficult to construct an interpolation practice which stops the cells from flipping and degenerating even if the cell-centred motion field is bounded. Moreover, motion of corner points (belonging to only one cell) and intersections of free-moving boundaries cannot be reconstructed reliably. Finally, while the FVM is unconditionally bounded for the convection operator, on badly distorted meshes one needs to sacrifice either the second-order accuracy or boundedness in the Laplacian, due to the explicit (and unbounded) nature of the non-orthogonal correction [24]. A combination of the two has forced us to quickly abandon this approach with the lesson that a point-based solution for the motion equation is essential.
3.2.2 Spring Analogy The first obvious point-based solution strategy is the spring analogy [7]. Here, all edges in the mesh are replaced by elastic springs and point motion is obtained by “loading” the spring system with the prescribed boundary motion. The linear system of motion equations is obtained from the force balance in all points. While this seems intuitively right, a number of failure modes has been observed, especially on polyhedral cells. This has occured in spite of the fact that the system of linear equations describing point motion has been solved to machine tolerance before moving the mesh. In other words, the failure is associated with the final solution rather than intermediate mesh states or the relaxation procedure. Let us start by imagining a linear spring analogy solver and examining various failure modes using the rules from Section 3.1.3. After each failure, an improvement on the method will be suggested. 9
c’ a’
b’
a
b
c (a) Coincident points.
(b) Triangle flip.
Fig. 3. Failure modes for triangle motion in spring analogy.
The simplest failure mode is a situation where two points a and b of a triangle (separate at the beginning of the motion) end up on top of each other, Fig. 3(a), and the triangle degenerates into a line. The cause of the failure is the fact that with linear spring analogy the force resulting from an edge being collapsed to a point does not tend to infinity. This can be easily remedied, by introducing non-linearity, i.e. by making the stiffness coefficient length-dependent, using e.g. the exponential spring law. The price, however, is substantial: a linear elastic problem has been converted into a non-linear problem and an iterative solution procedure is required. Additional problems exist with relation to unloading (long non-linear springs tend to relax faster than the short ones), but for the case of presentation, we may consider this problem as solved. The second mode of failure is a triangle flip, Fig. 3(b). Here, a triangle becomes degenerate by pushing point c through the opposite edge without ever degenerating any of its edges. No edges are reduced to zero length and the non-linear spring fix mentioned above will not remedy the problem. Farhat et al. [10,11] propose a solution by introducing torsional springs in points which control internal edges of triangles. Additional equations are now obtained from the torque balance in all points. One can easily infer that a torsional spring will need to be replaced by its non-linear equivalent to produce a torque moment tending to infinity as the angle between two adjacent edges approaches zero, with additional corrections on unloading. Additionally, the problem of torsional springs becomes substantially more complex on 3-D polyhedra in relation to a tetrahedron flip. Thus, a reliable spring analogy setup would consist of linear and torsional springs with length/angle-dependent coefficient sensitised on loading/unloading, Fig. 10
Fig. 4. Combined linear-torsional spring analogy.
4. Taking a step back, it seems that a reliable spring-based mesh motion system has come at a very high price: a non-linear component-coupled system which is difficult and costly to solve. At the same time, the FEM community uses a linear Laplace equation with impunity, producing reliable and bounded motion fields at a fraction of a price, indicating that a fundamental re-think may be required.
4
Polyhedral Vertex-Based Motion Solver
Looking at the previous section, it seems that spring analogy continually attempts to overcome numerical deficiencies by introducing non-linearity. At the same time, one could clearly claim that the Laplacian operator is a perfect choice to govern mesh motion 2 , as it is unconditionally bounded; it only remains to devise a solution strategy that will preserve boundedness in the discrete form. One could claim that spring analogy is no more than a numerically clumsy approach to solving a Laplacian-based mesh motion equation on the points of the mesh. Taking a lead from the FEM practice, we can state the following requirements on the mesh motion solver: • A vertex-based solution method, avoiding the need for interpolation; • Use of iterative solvers for efficiency, implying diagonally dominant matrices resulting from discretisation; • No triangles or tetrahedra in the cell decomposition should be inverted – this should be guaranteed by the discretisation practice. 2
We shall examine the pseudo-solid approach as an alternative choice later in the text.
11
4.1 Tetrahedral Elements and Discretisation of a Laplacian The “obvious” solution of using a classical FEM solver to solve for mesh motion is rejected: to the authors’ knowledge, automatic definition of shape functions for arbitrary polyhedra does not exist and it is impractical to categorise all “allowed” cell shapes in the face-addressed mesh. Also, it is unclear whether such a shape function would produce a diagonally dominant matrix we are seeking for efficiency reasons. At the same time, tetrahedral finite elements for a Laplacian produce a diagonally equal matrix and second-order discretisation. Also, matrix coefficients tend to infinity for both degenerate situations in Fig. 3. We can prove this by considering a trick for calculating the off-diagonal matrix coefficient for a tetrahedron in real space, Fig. 5.
l k
si
i
j
Fig. 5. Calculating matrix coefficients for a tetrahedron in real space.
The coefficient contribution for a point pair (i, j) and the Laplacian operator can be calculated as: aij =
Z
VT
∇Ni •∇Nj dV =
si • s j . 9 VT
(6)
Here, N is the element shape function, si is the surface-normal vector on the triangular face opposite point i and VT is the volume of the tetrahedron. Consider a case where a tetrahedron approaches a degenerate state, either by a point approaching another point or the opposite face. In this case, the face area vectors si and sj will still be finite (and come close to being parallel), whereas the volume will approach zero. Thus, as the denominator in Eqn. (6) approaches zero and the nominator remains finite, aij tends to infinity. As a result, the tetrahedral FEM discretisation will remain bounded irrespective of mesh quality and, if used to solve for point motion, it guarantees no tetrahedra will be collapsed or flipped. 12
The final step in the analysis is the recognition that the cell measures we are trying to preserve are obtained by using a tetrahedral decomposition (which can be performed irrespective of the shape): it is precisely these tetrahedra whose quality we need to preserve. The matrix discretising the motion equation governed by a Laplacian as described above will be sourceless and diagonally dominant, allowing the use of iterative solvers. The motion variable will be bounded irrespective of mesh quality, thus answering all the motion solver requirements stated above.
4.2 Choice of Cell Decomposition
It remains to choose an appropriate decomposition of a polyhedron into tetrahedra; two methods investigated in this study are shown in Fig. 6.
(a) Cell decomposition.
*) )*! ) ) ) ) ) # ) # ) # ) # ! $ ! $ ! $ $ *! ! * ! * ! * ! * ! * ! * ! * ! * )*! ) ) ) ) ) ) ) ) # # # # ! * ! * ! * ! * ! * ! $ ! * ! $ ! * ! $ ! * $ )*)! )*! )*! )*! )*! )*! )*! )*! )*! #$! #$! #$! $! $! $! $#$#$# *)*)*)) *! *! *! *! *! *! *! *! ) ) ) ) ) # ) # ) # ) ) ) ) ) ) # ) # ) # ) ! $ ! $ ! $ *)! ! * ! * ! * ! * ! * ! * ! * ! * ))! ))! ))! ))! ))! ##! ))! ##! ))! ##! ))! $! $! $! $##$#$ **))**) *))! *! *! *! *! *! *! *! *! ! **! **! **! **! **! $$! **! $$! **! $$! **! )*! )**! ) ) ) ) ) ) ) # # # ! $ ! $ ! $ *! ! * ! * ! * ! * ! * ! * ! * )*)! ) ) ) ) ) # ) # ) # ) $#$#$# *)*)) )*! )*! )*! )*! )*! #$! )*! #$! )*! #$! )*! $! $! $! *! *! *! *! *! *! *! *! ) ) ) ) ) # ) # ) # ) *)! ) ) ) ) ) ) ) ) # # # *)! ! * ! * ! * ! * ! * ! $ ! * ! $ ! * ! $ ! * $#$#$# **)*)*) ) ) ) ) ) ) ) ) # # # *)! ! * ! * ! * ! * ! * ! $ ! * ! $ ! * ! $ ! * )*! ) ) ) ) ) # ) # ) # ) ! $ ! $ ! $ *! ! * ! * ! * ! * ! * ! * ! * ! * )*)! )*! )*! )*! )*! )*! #$! )*! #$! )*! #$! )*! $! $! $! $#$#$# *)*)*) *! *! *! *! *! *! *! *! )*! )*! )*! )*! )*! #$! )*! #$! )*! #$! )*! ) ) ) ) ) ) ) ) # # # *)! + + ,++*)! , , ! " ! " ! " ! " " ) ) ) ) ) ) ) ) # # # * ! * ! * ! * ! * ! $ ! * ! $ ! * ! $ ! * $#$#$# *)21*) ! ! ! ! ! + + ""! 2! ""!2! ""!2! ""!2! #$! #$! #$! $! $! 10/.-*,)2! 10/*)2! 10/*)2! 10/*)2! 1*)2! 1*)2! 1*)2! 210/.-,*)2! 210/.-*,)2! 2! 2! 2! "! "! "! "! """(' $! #&! #&! #&! / / / / / 10/.-210/.-! 1 1 1 1 1 1 1 1 . 0 . 0 0 0 0 # # # $ $ $ / / / / / 1 1 1 1 1 1 1 1 . 0 ! 2 . 0 ! 2 0 ! 2 0 ! 2 0 ! 2 ! 2 ! 2 ! 2 ' ' ' ' % % % ( ( ( ( &%$#&%&% 21212121 1 1 1 1 1 1 1 1 21.-! . ! 2 . ! 2 ! 2 ! 2 ! 2 ! 2 ! 2 ! 2 ' ' ' ' ' % % % ! & ! & ! & ! ( ! ( ! ( ! ( ( 1 1 1 1 1 1 1 1 21.-! ! 2 ! 2 ! 2 ! 2 ! 2 ! 2 ! 2 ! 2 . . ' ' ' ' ' % % % ! & ! & ! & ! ( ! ( ! ( ! ( ( 1.-2! 1.-2! 12! 12! 12! 12! 12! 12! 21.-! 2! 2! 2! 2! 2! '(! '(! '(! '(! % 2! % 2! % 2! (! 1.-.! 1.-.! 1.-.-- (! 1 (! 1 (! 1 ('('(' &! 1 &! 1 &! 1 &% 2121 2! '(! '(! '(! '(! -.-! ' ' ' ' ! . . ' ' ' ' ' ! ( ! ( ! ( ! ( ( --! .--! .! .--..- (! '(! '(! '(! '(! (! (! (! ('('(' ! ..-! ..! '(! '(! '(! '(! -.! '(! ' ' ' -.! -.! .! . '.! ' ' ' ! ( ! ( ! ( ('('(' -.-! . '(! '(! '(! '(! (! (! (! -.! .-.-.- (! ''! ''! ''! ''! .-! '(('' ( ( ( ( .-! ! . ' ' ' ' ! ( ! ( ! ( ! ( .-! ! . . ' ' ' ' ! ( ! ( ! ( ! ( (('(' -.! .! ! . . ' ' ' ' ! ( ! ( ! ( ! ( -.-! ! . . ' ' ' ' ! ( ! ( ! ( ! ( ! . . ' ' ' ' ! ( ! ( ! ( ! ( -.! .-! '(! '(! '(! '(! - (! .-!.! ' .-.- (! ' (! ' (! ' ('('('
(b) Cell-and-face decomposition.
Fig. 6. Decomposing a polyhedral cell into tetrahedra.
Consistent with face decomposition, Fig. 2, a cell is decomposed by introducing a point in its centroid and building tetrahedra above the triangular decomposition of a face. The two methods proposed here are the cell decomposition, Fig. 6(a), where additional points are introduced only in cell centres; and cell-and-face decomposition, Fig. 6(b), where points are introduced in both face and cell centres. The choice of decomposition is a balance between the convexness requirement, quality of resulting tetrahedra and the computational cost: cell decomposition 13
introduces an additional point in all cell centres, whereas cell-and-face decomposition adds a point in all face centres as well, giving a considerable increase in the number of unknowns. The issue of solution cost will be investigated further in Section 4.4. 4.3 Choice of Motion Equation In the FE mesh motion framework, we can choose several candidate equations to govern mesh motion. The most obvious choices are the Laplace equation with constant and variable diffusivity [12] and the small-strain formulation of the linear elastic model [16]. The difference between the two is that the former only allows direction-decoupled transfinite mapping, while the latter also allows rotation which may be preferable under some circumstances. However, this comes at a relatively high price: the pseudo-solid equation also couples the motion vector components (due to rotation), [28]. The choice here is either a further increase in storage associated with the coupled solution or an iterative segregated solution method. Numerical experiments in this direction show that the substantial increase in cost associated with the pseudo-solid equation is not justified by improved mesh quality; for this reason, the Laplace equation will be used in the rest of this study. Some evidence exists that introduction of variable diffusivity in the Laplacian can give substantial control over mesh quality. However, it is not immediately clear how to formulate the variable diffusivity field and this will be further investigated in numerical studies. As a guidance, increased diffusivity corresponds to increased stiffness of the “equivalent spring system”: a choice for governing the diffusivity can thus be a measure of mesh quality or distance to some chosen boundary. 4.4 Efficiency Concerns A critical requirement in this work is the motion solver which co-exists in the code with the FVM flow solver and matches it for efficiency of matrix assembly, execution speed and storage requirements. Unlike the classical face-addressed FV solver [26], where the matrix is assembled by looping over all mesh faces for all the operators and calculating the complete coefficient, the FEM assembles the matrix by looping over all elements and distributing coefficient contributions into the diagonal, off-diagonal and source vectors. A concern is the fact that the number of tetrahedra in the mesh will be considerably higher than the number of cells in the mesh. For example, in an all-hexahedral mesh of 100 000 CVs, the number of tetrahedra will 14
be 1.2 million for the cell decomposition or 2.4 million for the cell-and-face decomposition. Fortunately, the decomposition can be done “on-the-fly”, without storing the actual tetrahedra and, combined with an efficient matrix assembly algorithm, poses limited storage and matrix assembly overhead. Furthermore, additional efficiency is obtained by operating on a group of tetrahedra at a time (resulting from a single polyhedral cell) with no storage overhead. This is sometimes called the “mini-element” technique. Taking the efficiency concern a step further, we can compare the size of the pressure matrix in the FV solver (1 equation per cell) with the size of the motion matrix. In spite of the large number of tetrahedra, the number of equations is not excessive: a 100 000 cell mesh in 3-D produces a motion matrix of approx. 200 000 (cell decomposition) to 500 000 (cell-and-face decomposition); it remains to be seen whether the solution cost is acceptable. Looking at the polyhedral decomposition, Fig. 6, it can be seen that the cell centre point in the decomposition is connected only to the points of the current cell. We could thus examine the possibility of eliminating the equation for the centre point before assembling the matrix as point motion is needed only in the actual points of the mesh. However, it turns out that the elimination of additional variables would severely deteriorate the matrix condition number and preclude the use of iterative solvers. This option would therefore result in a substantial increase in the solution cost and is abandoned. In terms of storage requirements, one should note that in segregated FVM fluid flow solvers, the memory peak occurs during the pressure-velocity solution (using SIMPLE [29] or PISO [30]). Here, it is necessary to simultaneously store the momentum and pressure matrices. The mesh motion solver operates either before of after the pressure-velocity module and some storage will be reused (subject to dynamic memory handling within the code). This somewhat decreases the perceived storage peak of the motion solver relative to the FVM part of the algorithm. Simulations involving dynamic (or solution-dependent) mesh motion are typically done in the transient mode – this fact can be used for further optimisation of the motion solver. Apart from the choice of the motion equation, we are free to choose the motion variable to be either point position or motion velocity. With the use of iterative solvers and in transient simulations one can assume that the motion velocity field changes slower than point position and thus a better initial guess is available. For stationary meshes the velocity solution equals to zero everywhere and is less polluted by round-off errors than the (linear) point position field. For constant-velocity deformation the cost of solving the motion equation in terms of velocity becomes trivial (a good initial guess is available), which is not the case if point position is chosen as the primitive variable. 15
For better precision, the motion velocity on the boundary is calculated from the current and desired point position and the time-step. This approach avoids the accumulation of round-off errors associated with solving for motion velocity and using point position.
4.5 Final Form of the Motion Solver
In summary, the polyhedral mesh motion solver is constructed as follows: (1) Every polyhedral cell is split into tetrahedra by splitting its faces into triangles and introducing a point in cell centroid. Consistency in tetrahedral connectivity is obtained by using identical face decomposition for both cells sharing an internal face. (2) The Laplace operator: ∇•(γ∇u) = 0
(7)
with constant or variable diffusion field γ is chosen to govern mesh motion. Here, u is the point velocity field used to modify point positions: xnew = xold + u∆t,
(8)
where xold and xnew are the point positions before and after mesh motion and ∆t is the time-step. Eqn. (7) is discretised on the tetrahedral decomposition using standard second-order finite element method and produces a diagonally equal matrix. For efficiency reasons, the matrix coefficients are calculated in real space using Eqn. (6). (3) Boundary conditions for the motion equation are enforced from the known boundary motion; this may include free boundaries, symmetry planes, prescribed motion boundary etc. (4) The matrix is solved using an iterative linear equation solver; here the choice falls on the Incomplete Cholesky preconditioned Conjugate Gradient (ICCG) solver [31], also used by the FVM solver.
4.6 Implementation of the Mesh Motion Solver
Several code organisation issues need to be examined, dealing with component interaction and the actual implementation of the proposed motion solver. One could consider it an overkill to implement a fully-fledged FEM solver in order to move the mesh in an existing FVM code. In this study, the motion solver is implemented in FOAM [32,33], an object-oriented C++ continuum 16
mechanics library. The software is constructed to allow extensive code reuse, typically impossible in more traditional designs. FOAM currently implements a second- and fourth-order collocated FVM on arbitrarily unstructured meshes. It is written in operator form and has a class hierarchy designed to be shared between various discretisation practices. Examining the code structure, it emerges that the classes representing the computational mesh, fields and field algebra, boundary conditions, linear equation matrix and solver technology as well as problem setup and post-processing can be shared. All the components listed above are intrinsically discretisationindependent, totalling over 100 000 lines of code and provide an excellent headstart. What remains to be implemented are the FEM-specific parts: derived mesh handling adapted for the FE discretisation (and using raw mesh information from the base class), FE calculus (e.g. divergence and gradient operators etc.) and the actual FE discretisation operators (e.g. convection, diffusion, time derivative etc.) with boundary condition handling. This totals just under 21 000 lines of code, including both cell decomposition methods, and the handling of solver parallelisation (to be described in future publications). The code architecture allows us to keep the FEM implementation on the tetrahedral decomposition separate from the rest of the code, but still relying on the common components mentioned above. It is encapsulated in its own library and loaded on demand. The FE solver is also implemented in the operator form, currently providing the temporal derivative ∂u and various form ∂t of the div-grad operators: ∇•(γ∇u), ∇•(γ∇uT ) and ∇•[γI(∇•u)]. The solver has been validated in isolation by solving linear elasticity model problems. The actual motion solver is implemented by using the discretisation operators in the FEM library and packed for ease of use in a separate module, together with the necessary mesh checking and setup tools.
5
Examples of Mesh Motion
Numerical experiments show that the cell decomposition is sufficiently robust for 2-D and “trivial” 3-D meshes; for simulations on industrial-standard meshes with considerable geometrical complexity and cell distortion, the celland-face decomposition is preferred. Note that for the mesh motion solver, the original mesh will be decomposed into tetrahedra (as shown in Section 4.2) even if it is already composed exclusively of tetrahedra: in the face-based mesh format, this fact is not immediately obvious. 17
We shall now apply the novel motion algorithm on two test problems and examine the mesh quality for various definitions of non-constant diffusion fields.
5.1 Motion of a Cylinder The first test case consists of a triangular mesh in 2-D and a circle moving within it 3 . Constant motion velocity is prescribed for the points on the circle, the left and right boundary are fixed and a slip condition is prescribed at the top and bottom. Fig. 7 shows the mesh at several positions. It can be clearly seen that while the mesh remains valid, its quality is not satisfactory. The main objection is excessive mesh density in front of the moving cylinder and the overall reduction in quality due to stretching. It is still encouraging that even for such a large mesh deformation, the final mesh remains valid.
Fig. 7. Motion of a cylinder in a duct.
The simplest way of improving mesh quality is by introducing variable diffusivity. The rationale here is that the mesh quality is more important close to boundaries, whereas the “core” mesh has more freedom to deform. One could also postulate various “cell quality” indicators: as a cell becomes more 3
In reality, the mesh is 3-D and consists of prismatic elements, as the software only operates on 3-D meshes. Although this is not strictly necessary, motion equation is solved for all points in the mesh.
18
distorted, its further deformation should be avoided. In both cases, the way to diminish local distortion is by increasing local diffusivity in the motion operator. Note that the Laplacian operator remains bounded for arbitrary variation of the positive diffusion coefficient. Out of a variety of diffusion laws one could devise, the following have been tested in this study: • Distance-based methods, where a number of boundary patches are selected by the user and the diffusion field γ is a function of cell centre distance l to the nearest selected boundary: · Linear (D1): γ = 1l ; · Quadratic (D2): γ = l12 ; · Exponential (D3): γ = e−l . • Quality-based methods, where the diffusion field γ is a function of a cell quality measure: · Mean cell non-orthogonality (Q1); · Mean cell skewness (Q2); · Mixed (Q3) – a combination of Q1 and Q2.
(a) Constant.
(b) Linear distance-based (D1).
(c) Quadratic distance-based (D2). Fig. 8. Influence of diffusivity on mesh quality.
Fig. 8 shows the mesh comparison for various choices of diffusivity where the 19
moving circle has been chosen as the reference boundary. It can be clearly seen that distance-based diffusivity improves the overall mesh quality and especially mesh spacing near the circle. The result for quality-based method (Q1) was almost identical to Fig. 8(a), i.e. no significant improvement has been observed.
5.2 Pitching Airfoil The second test case consists of a pitching airfoil with a 2-D hybrid mesh. The airfoil moves according to the sinusoidal law, including both translation and pitching. Fig. 9 shows the mesh around the airfoil during three phases of motion using the D2 distance-based diffusivity. Of particular interest is the mesh around the trailing edge. The selected method preserves the mesh quality very well.
Fig. 9. Motion of a NACA airfoil, D2 diffusivity method.
A comparison of mesh quality for four diffusivity methods is shown in Fig. 10. In this case, all distance-based methods seem to perform very well, giving notable improvement over constant diffusivity. A quantitative confirmation of improved mesh quality is given in Fig. 11. Distance-based methods produce meshes with lowest non-orthogonality and skewness, while the constant diffusivity method considerably deteriorates the mesh. However, in all cases the resulting mesh remains valid. It is interesting to 20
(a) Constant.
(c) Quadratic (D2).
(b) Linear distance-based (D1).
distance-based
(d) Quality-based (Q1).
Fig. 10. Mesh quality around the trailing edge.
see how the mesh quality impacts on the computational cost. Fig. 11(c) plots the number of iterations in the ICCG solver for the motion equation solved to the same relative tolerance (1 × 10−9 ). For the methods with adverse effect on mesh quality, the solution cost increases considerably. The D2 solution method consistently requires around 130 solver iterations, comparable with the effort of the FV pressure solver for the same case. In realistic simulations, the transient nature of the flow introduces a time-step limitation on the FVM solver used for the flow equation, somewhat mitigating the cost of the motion solver. Also, in the iterative solution framework, for most cases it is enough to reduce the residual by approximately 2 orders of magnitude instead of converging the solution to machine tolerance for every time step. Partial convergence only becomes an issue on extremely distorted meshes (which can no longer be used for the FV solution) or in cases where one boundary comes to close proximity of another, e.g. Fig. 7. In conjunction with partial convergence, one can resort to an a-priori indicator of motion validity. The simplest solution is to manually check the mesh before executing the motion. However, this is expensive and unnecessary and 21
80
220 constant D1 D2 D3 Q1 Q2 Q3
75 70
180
65
160
60
Maximal skewness
Maximal non-orthogonality
constant D1 D2 D3 Q1 Q2 Q3
200
55 50 45
140 120 100 80
40
60
35
40
30 25
20 0
0.5
1
1.5
2
0
0.5
1
t, s
1.5
2
t, s
(a) Mesh non-orthogonality.
(b) Mesh skewness.
700 constant D1 D2 D3 Q1 Q2 Q3
600
Number of iterations
500
400
300
200
100
0 0
0.5
1 t, s
1.5
2
(c) Number of iterations. Fig. 11. Mesh quality indicators and computational effort.
a simpler method can be proposed. Consider two points A and B along an edge e = AB in space. The requirement of two points “crossing over” reduces to the condition of e•(∇u)•e < 1. If desired, this can be checked for all mesh edges before executing motion. Experience shows that the most robust way of dealing with problematic motion solutions is to halve the time-step size and/or tighten the solution tolerance before recalculating the motion solution for the current time-step. Comparing the performance of the two polyhedral decomposition methods, a curious fact has been noted. On reasonably complex meshes in 3-D it turns out that the cost of solution of the cell-and-face decomposition in total execution time becomes substantially lower than that for the cell decomposition. This is counter-intuitive, as the number of unknowns for the cell-and-face decomposition is considerably greater (in both cases, it is only the point positions that are used; the solution in cell and face centres is discarded). The crossover in solution cost is due to an increased number of iterations in the ICCG solver for the cell decomposition method. It seems that the higher quality of decomposition tetrahedra for the cell-and-face decomposition creates a better22
conditioned matrix which, for the same mesh, sometimes requires only a tenth of the number of iterations used with the cell decomposition.
6
Free Surface Flow Simulations
This study will be concluded by presenting two simulations of free surface flows using a surface tracking algorithm. Here, the fluid equations are solved in both phases and coupled across the free surface. The free surface is represented as a mesh interface whose motion depends on the flow solution. A schematic representation of the free surface condition is given in Fig. 12. Liquid
free surface(fs) pf s = pf ss − σKf s (∇u)f s =
Gas
µgas (∇u)f ss µliquid
free surface shadow (fss) (∇p)f ss = 0 uf ss = uf s
Fig. 12. Motion of a NACA airfoil, D2 diffusivity method.
Second-order FVM is used for the fluid flow and the automatic mesh motion solver described above adjusts the mesh. On the free surface, a double boundary condition is imposed: fixed pressure and zero flux condition need to be satisfied simultaneously. The fluid flow equations are solved using a segregated SIMPLE procedure, taking into account the kinematic and dynamic condition on the free surface, as well as surface tension. The “no-flux” condition is satisfied in an iterative sequence [34,35], providing the boundary condition for mesh motion on the free surface. The solution procedure enforces the fixed pressure boundary condition on the free surface and consequently a non-zero-flux is obtained. Position of the faces in the free surface patch is adjusted such that the area swept in motion equals the flow flux for the face. Clearly, the change in domain shape influences the pressure and flux solution and the procedure is repeated in an iterative manner for every time-step until the fixed surface pressure and zero flux condition are satisfied simultaneously. 23
6.1 Hydrofoil Under a Free Surface The first free surface flow simulation consists of an inclined NACA 0012 hydrofoil placed 0.203 m below a free surface in water at an angle of attack of 5◦ [36]. The Froude number for the flow is F r = 0.567. In this simulation only the water component is modelled. The mesh consists of 14 508 CVs. Fig. 13 shows the mesh and pressure distribution for an instant in time. As a result of the pressure distribution around the profile, waves are induced on the free surface and convected downstream. Motion of the free surface close to the inlet on the left is constrained to represent constant inlet depth. The wave motion only occurs in the vertical direction; in the horizontal direction the mesh remains undisturbed.
Fig. 13. Hydrofoil under a free surface.
6.2 Free-Rising Air Bubbles in Water The driving force behind this study is a desire to assemble a tool for Direct Numerical Simulation (DNS) of air bubbles in water, with the aim of providing lift and drag data needed for two-phase Eulerian modelling. Similar calculations, using the surface capturing method have been reported by Rusche [37], indicating some deficiencies in the handling of surface tension. A surface-tracking approach has the potential of treating surface tension to much higher accuracy and, provided the method is reliable and efficient, can be compared with [37]. Here, we shall report some initial results for free rising air bubbles in water in 2- and 3-D. The bubble is located in a large box and the flow boundary conditions are 24
adjusted such that it remains centred in the domain. The actual bubble trajectory and rising velocity transients can be followed through the changes in the boundary condition. In this simulation, the bubble rises through a quiescent fluid and the material properties of air and water are used including the surface tension effects. Note that, unlike in the surface capturing methods, the strength of surface tension or the jump is the density and viscosity do not pose a problem. Fig. 14 shows the mesh deformation and the pressure field around a 2-D air bubble of 1 mm diameter freely rising in water. After the initial transient, the bubble reaches terminal velocity and shape. The mesh in this simulation consists of 12 840 CVs in two disconnected regions and captures the interface coupled through the free surface condition.
Fig. 14. Free-rising air bubble in water in 2-D: pressure iso-lines and surface deformation.
A timing breakdown for a single iteration are given in Table 1. The simulation is performed on a Linux computer with a 2 GHz Intel Pentium IV processor with 1 GB of memory. Finally, Fig. 15 shows the flow field around the free rising air bubble of the same diameter in 3-D. The mesh consists of 561 920 cells with the near-surface resolution sufficient to resolve surfactant transport effects. A detailed break25
Table 1 Timing breakdown for a 2-D bubble simulation. Operation
Time [s]
Cumulative [s]
Building momentum matrix
0.12
0.12
Solving momentum equation
0.06
0.18
Building pressure matrix
0.06
0.24
Solving pressure equation
0.33
0.57
Building motion matrix
0.27
0.84
Solving motion equation
0.70
1.54
down of the timing for a single iteration of the simulation is given in Table 2 for the same platform as before.
Fig. 15. Free-rising air bubble in water in 3-D: pressure iso-lines and surface pressure. Table 2 Timing breakdown for a 3-D bubble simulation. Operation
Time [s]
Cumulative [s]
Building momentum matrix
6.34
6.34
Solving momentum equation
3.16
9.50
Building pressure matrix
3.09
12.59
Solving pressure equation
29.02
41.61
Building motion matrix
12.80
54.41
Solving motion equation
39.95
94.36
26
In both simulations, the cost associated with the mesh motion solver is 50 − 60% of the complete cost of simulation, which is high but acceptable. It has been noted that the cost balance is more favourable in 3-D and on large meshes, due to the change in the balance of the number of cells and points. However, the selected mesh motion algorithm is inherently parallel, both in terms of selected discretisation and the choice of linear equation solvers. A combination of a massively parallel FVM flow solver already available in FOAM and a parallel motion solver working on the identical mesh decomposition offers considerable scope in terms of reduced execution time per time-step. Good parallel efficiency seems to be the way to afford the cost of long transient runs needed to accumulate sufficient DNS statistics.
7
Summary and Future Work
In this study we have examined the requirements of deforming mesh simulations in a Finite Volume framework, with the objective of developing a robust and reliable automatic mesh motion tool. The purpose of a mesh motion solver is to determine the point positions for the mesh based on the prescribed boundary motion – this can be prescribed by external events or calculated as a part of the solution. While performing automatic mesh motion, it is crucial to preserve the validity and quality of the mesh. Having analysed several popular automatic mesh motion approaches and their advantages and drawbacks, we have settled on a second-order quasi-tetrahedral Finite Element method and the Laplace operator to govern the motion. Support for polyhedral cells is provided using the “mini-element” technique, where each polyhedron is, for the purposes of motion discretisation, split into tetrahedra and a second-order shape function is used. Analysis shows that it is precisely those tetrahedra that control the quality and validity of the faceaddressed unstructured FV mesh and that the chosen method of discretisation guarantees to preserve their quality. Furthermore, the chosen method of discretisation produces a symmetric positive definite matrix ideal for iterative linear equation solvers. The quality of the mesh in motion is preserved by prescribing non-constant diffusion field in the Laplace operator. Several techniques have been tested, most notably the distance-based diffusion, where the coefficient depends on the distance between the cell centre and the nearest boundary of interest. A combination of the above components with a second-order FV flow solver creates a robust and efficient dynamic mesh motion solver capable of handling free surface flows using a surface-tracking algorithm. The solver has been tested on two free surface flows, including a simulation of free rising air 27
bubbles in water. Overall, the cost of the automatic motion solver is about 50 − 60% of the overall cost of simulation. In future work, the flow solver will be used in DNS simulations of gas bubbles in liquids, establishing a base for phase interaction modelling in Eulerian multiphase simulations.
References [1] F. Harlow, J. Welsh, Numerical computation of time dependent viscous flows with free surface, Phys. Fluds 8 (1965) 2182–2189. [2] C. Hirt, B. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comp. Physics 39 (1981) 201–225. [3] O. Ubbink, R. Issa, ”a method for capturing sharp fluid interfaces on arbitrary meshes, J. Comp. Physics 153 (1999) 26–50. [4] M. Behr, T. Tezduyar, Finite element solution strategies for large-scale flow simulations, Comput. Methods Appl. Mech. Engrg. 113 (1994) 3–24. [5] I. G¨ uler, M. Behr, T. Tezduyar, Parallel finite element computation of freesurface flows, Computational Mechanics 23 (1999) 117–123. [6] J. Batina, Unsteady Euler airfoil solutions using unstructured dynamic meshes, AIAA Journal 28 (8) (1990) 1381–1388. [7] F. Blm, Considerations on the spring analogy, Int. J. Num. Meth. Fluids 32 (2000) 647–668. [8] B. Perot, R. Nallapati, A moving unstructured staggered mesh method for the simulation of incompressible free-surface flows, J. Comp. Physics 184 (2003) 192–214. [9] Y. Zhao, J. Tai, F. Ahmed, Simulation of micro flows with moving boundaries using high-order upwind method on unstructured grids, Computational Mechanics 28 (2002) 66–75. [10] C. Farhat, C. Degand, B. Koobus, M. Lesoinne, Torsional springs for twodimensional dynamic unstructured meshes, Comput. Methods Appl. Mech. Engrg. 163 (1998) 231–245. [11] C. Degand, C. Farhat, A three-dimensional torsional spring analogy method for unstructured dynamic meshes, Computers and Structures 80 (2002) 305–316. [12] R. L¨ohner, C. Yang, Improved ALE mesh velocities for moving bodies, Communications in Numerical Methods in Engineering 12 (1996) 599–608. [13] D. Littlefield, The use of r-adaptivity with local, intermittent remesh for modeling hypervelocity impact and penetration, Int. J. Impact Engrg. 26 (2001) 433–442.
28
[14] A. Masud, T. Hughes, A space-time Galerkin/least-squares finite element formulation of the navier-stokes equations for moving domain problems, Comput. Methods Appl. Mech. Engrg. 146 (1997) 91–126. [15] I. Robertson, S. Sherwin, Free-surface flow simulation using hp/spectral elements, J. Comp. Physics 155 (1999) 26–53. [16] M. Behr, F. Abraham, Free-surface flow simulations in the presence of inclined walls, Comput. Methods Appl. Mech. Engrg. 191 (2002) 5467–5483. [17] R. Cairncross, P. Schunk, T. Baer, R. Rao, P. Sackinger, A finite element method for free surface flows of incompressible fluids in three dimensions. Part I. Boundary fitted mesh motion, Int. J. Numer. Meth. Fluids 33 (2000) 375–403. [18] G. Chiandussi, G. Bugeda, E. O˜ nate, A simple method for automatic update of finite element meshes, Commun. Numer. meth. Engrg. 16 (2000) 1–19. [19] A. Johnson, T. Tezduyar, Mesh update strategies in parallel finite element computations of flow problems with moving boundaries and interfaces, Comput. Methods Appl. Mech. Engrg. 119 (1994) 73–94. [20] P. Sackinger, P. Schunk, R. Rao, A Newton-Raphson pseudo-solid domain mapping technique for free and moving boundary problems: A finite element implementation, J. Comp. Physics 125 (1996) 83–103. [21] M. Souli, J. Zolesion, Arbitrary Lagrangian-Eulerian and free surface methods in fluid mechanics, Comput. Methods Appl. Mech. Engrg. 191 (2001) 451–466. [22] H. Zhou, J. Derby, An assessment of a parallel, finite element method for three-dimensional, moving-boundary flows driven by capillarity for simulation of viscous sintering, Int. J. Numer. Meth. Fluids 36 (2001) 841–865. [23] B. Helenbrook, Mesh deformation using the biharmonic operator, Int. J. Numer. Meth. Engrg. 56 (2003) 1007–1021. [24] H. Jasak, A. Gosman, Automatic resolution control for the finite volume method. Part 1: A-posteriori error estimates, Numerical Heat Transfer, Part B 38 (3) (2000) 237–256. [25] I. Demirdˇzi´c, M. Peri´c, Space conservation law in finite volume calculations of fluid flow, Int. J. Num. Meth. Fluids 8 (9) (1988) 1037–1050. [26] H. Jasak, Error analysis and estimation in the finite volume method with applications to fluid flows, Ph.D. thesis, Imperial College, University of London (1996). [27] H. Jasak, H. Weller, A. Gosman, High resolution NVD differencing scheme for arbitrarily unstructured meshes, Int. J. Numer. Meth. Fluids 31 (1999) 431–449. [28] H. Jasak, H. Weller, Application of the finite volume method and unstructured meshes to linear elasticity, Int. J. Num. Meth. Engineering 48 (2) (2000) 267– 287.
29
[29] S. Patankar, Numerical heat transfer and fluid flow, Hemisphere Publishing Corporation, 1981. [30] R. Issa, Solution of the implicitly discretized fluid flow equations by operatorsplitting, J. Comp. Physics 62 (1986) 40–65. [31] M. Hestens, E. Steifel, Method of conjugate gradients for solving linear systems, Journal of Research 29 (1952) 409–436. [32] H. Weller, G. Tabor, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using object orientated techniques, Computers in Physics 12 (6) (1998) 620 – 631. [33] H. Jasak, H. Weller, N. Nordin, In-cylinder cfd simulation using a c++ objectoriented toolkit, SAE Technical Paper 2004-01-0110 (2004). [34] S. Muzaferija, M. Peri´c, Computation of free surface flows using interface-tracking and interface-capturing methods, Computational Mechanics Publications, Southampton, 1998, Ch. 3. ˇ Tukovi´c, H. Jasak, Unstructured finite volume free surface tracking algorithm [35] Z. with automatic mesh motion, to be published. [36] J. Ferziger, M. Peri´c, Computational methods for fluid dynamics, Springer Verlag, Berlin-New York, 1995. [37] H. Rusche, Computational fluid dynamics of dispersed two-phase flows at high phase fractions, Ph.D. thesis, Imperial College, University of London (2003).
30
