Abstract
This paper discusses the numerical approximation of flow problems, in particular the large eddy simulation of turbulent flow and the simulation of laminar immiscible two-phase flow. The computationsfor both applications are performed with a coupled solution approach of the Navier-Stokes equations discretized with the finite element method .Firstly, a new implementation strategy for large eddy simulation of turbulent flow is discussed. The approach is based on the vibrational multiscale method, where scale ranges are separated by vibrational projection. The method uses a standard Navier–Stokes model for representing the coarser of the resolved scales, and adds a sub grid viscosity model to the smaller of the resolved scales.The scale separation within the space of resolved scales is implemented in a purely algebraic way based on a plain aggregation algebraic multigrid restriction operator. A Fourier analysis underlines the importance of projective scale separations and that the proposed model does not affect consistency of the numerical scheme. Numerical examples show that the method provides better results than other state-of-the-art methods for computations at low resolutions.Secondly, a method for modeling laminar two-phase flow problems in the vicinity of contact lines is proposed. This formulation combines the advantages of a level set model and of a phase field model: Motion of contact lines and imposition of contact angles are handled like for a phase field method, but the computational costs are similar to the ones of a level set implementation. The model is realized by formulating the Cahn–Hilliard equation as an extension of a level set model.The phase-field specific terms are only active in the vicinity of contact lines. Moreover, various aspects of a conservative level set method discretized with finite elements regarding the accuracy of force balance and prediction in jump of pressure between the inside and outside of a circular bubble are tested systematically. It is observed that the errors in velocity are mostly due to inaccuracies in curvature evaluation, whereas the errors in pressure prediction mainly come from the finite width of the interface. The error in both velocity and pressure decreases with increasing number of mesh points.
1. INTRODUCTION
The numerical simulation of fluid dynamics is one of the main fields in computational mathematics. Simulation is today both an alternative and a complement to experiments in many engineering disciplines as it helps in predicting the behavior of fluids. Moreover, simulation permits to systematically improve the design and material properties of products by means of optimization algorithms. Such a procedure would be tremendously expensive if it had to be done by prototyping and other experimental means. Simulation is thus a tool that renders technical processes more effective. One interesting area of research in fluid dynamics is the behavior of turbulent incompressible flow. Turbulent flows often occur in engineering applications, for example, in the air flow around a vehicle, the water flow in turbines of hydroelectric power plants as well as in water networks and oil pipelines with high throughput. For instance, relevant issues are the determination and the control of aerodynamic drag as well as efficiency of engines. Another area of active research is multiphase flow, a setting that involves two or more different fluids. There is a wide range of applications that involve multiple fluids and where numerical simulation is extensively used. An exceptionally challenging problem is computational combustion, where different fluids and flame fronts need to be represented. Often, the flows occurring in combustion
Are turbulent. But there are also a lot of multi-phase settings where the flow field is laminar and thus more regular. One such application is intravenous therapy in medicine. The objective of numerical simulations is to understand how bubbles consisting of oleaginous substances are transported and resolved in the circulation of blood. These processes are often strongly driven by surface tension, a force that strives after keeping different fluids separated. A similar application is the modeling of the absorption of air in the blood circulation in human lungs. Multiphase flow is a phenomenon that is particularly relevant also in subsurface flow. In this area, it is often assumed that the individual fluids are well separated, that is, the fluids do not mix and there are no additional particles resolved in the fluids. This thesis contributes to the simulation of such laminar well-separated two-phase flows. For subsurface flows, an additional difficulty arises, namely the interaction of systems consisting of different fluids with solids, called capillarity. This interaction takes place on so-called contact lines, locations where the interface (separating two fluids) meets a solid. The properties of the solid and the composition of the fluid give rise to a phenomenon called wetting, indicating which fluid preferably is in contact with the solid. Applications of this kind are the modeling of groundwater basins and oil extraction from subsurface reservoirs driven by water injection. The pores of rock where the fluid moves are often of the order of a few millimeters or less, whereas the actual simulation region can extend over thousands of square kilometers. In order to tackle problems with such an enormously wide range of scales, simplifications need to be applied. One example is to perform small-scale flow simulations with detailed rock structures in order to determine the permeability of various stone configurations. Such effective parameters are then used as input to a large scale simulation, usually based on Darcy’s law [17]. Other applications of multiple fluids in contact with solids are found in industrial processes like sintering, where mixtures of different fluids and their solidification are tracked; lab-on-a chip processes, which form the basic component for microsystem design in the biomedical industry; and inkjet printing.
This paper discusses the equations that describe the flow of individual fluids, the Navier–Stokes equations. Next, an introduction to turbulence is given, which is followed by a presentation of numerical models for the simulation of laminar two-phase flow.Thepapers included in this identify possible directions of future research.
2. NUMERICAL METHODS FOR THE NAVIER–STOKES EQUATIONS
The motion of fluids at low and medium velocities in a connected computational domain
(spatial dimension d = 2, 3) is given by the incompressibleNavier–Stokes equations [35], [62], [73],
The vector u denotes the d components of fluid velocity and p denotes the fluid’s (dynamic) pressure. The rank-2 tensor represents the viscous stress tensor for a Newtonian fluid. The fluid density is
denoted by ρ, and the fluid’s (dynamic) viscosity by µ. In single-phase fluids, density is constant due to the incompressibility assumption. The term f specifies external forces acting on the fluid.
The system of Navier–Stokes equations (2.1)–(2.2) is completed by a divergence- free initial velocity field,
With
And application-specific boundary conditions (cf. Gresho and Sani [32] for a discussion of various possibilities). In this paper, we consider flow subject to no-slip boundary conditions
And natural (inflow/outflow) boundary conditions
On
Wheren denotes the unit outer normal on the domain boundary ¡h and h is some external traction boundary force. The boundary subdomains ¡g and ¡h are assumed to be non-overlapping and to span the whole boundary. For the numerical approximation of the equations of fluid flows, there are various methods available. The most popular ones are the finite volume method, the finite difference method, and the finite element method. A comparative introduction of these methods can be found in Quartapelle [74]. The algorithms developed in this thesis are based on the finite element method (FEM). There is a vast literature discussing various aspects of the finite element discretization of system (2.1)–(2.2). Mathematical aspects are discussed, e.g., in Glowinski [28] and Gunzburger [35]. The book by Gresho and Sani [32] focuses on practical issues for implementation such as boundary conditions, discretization schemes, and stabilization. The work by Elman, Silvester&Wathen [23] and Turek [88] put special emphasis on algorithmic matters like efficient numerical linear algebra for the flow equations.
2.1 WEAK FORM AND SPATIAL DISCRETIZATION
The Navier–Stokes equations (2.1)–(2.2) depend on both space and time. Most discretization approaches split the approximation into separate approximations for spatial and temporal part. For setting up the spatial FE approximation, the first step is to rewrite the system of equations (2.1)–(2.2) in variational form. To this end, we define the space of admissible velocity solutions at a certain time instant by
, and the space of admissible pressure solutions by
thevariational problem corresponding to (2.1)–(2.2) is to find, at each instant in time, a pair
Suchthat
t (2.6)
Holds for all test function
the form
denotes the standard
Inner product.
The spatial discretization of system (2.6) replaces the (infinite-dimensional) solution function spaces Vu and Vpby finite-dimensional subspaces V hu and V h p.
Then, only the projection of theNavier–Stokes equations onto finite dimensional test function spaces V h u £V h p is considered, resulting in numerical approximations uh and ph. The discretization considered in this thesis use a decomposition of the computational domain into small subsets of characteristic size h, called elements. We mainly focus on quadrilateral elements in 2D and hexahedral (brick) elements in 3D. The basic functions that span Vu and Vpare chosen to be piecewise polynomials within the elements, and continuous over the whole domain . For the representation of the velocity u, equally shaped basis functions are used for each spatial component. For instance, we will denote a finite element approximation with quadratic basis functions for velocity and linear basis functions for pressure on the elements by Qd2 Q1. On each element, the basic functions for u are defined as a tensor product of Lagrangian interpolation polynomials of degree two in each coordinate direction [43], and the pressure basis functions as tensor product of linear Lagrangianinterpolants. Globally, the basis functions are nonzero in exactly one node, and zero on all the others. The support of the basis functions is confined to a patch of elements around the respective nodes.
2.2 STABILITY OF SPATIAL DISCRETIZATION
For the spatial discretization to be stable, a compatibility condition between the approximation of velocity and pressure must be satisfied, namely
Where the constant ¯ does not depend on the mesh size parameter h. This condition was independently derived by Ladyzhenskaya [58], Babuška [4], and Brezzi [12] and is called LBB condition or inf–sup condition. The condition ensures that the divergence matrix discretizing the is of full rank, and hence the discretized system has a unique solution for velocity and pressure.
The element so-called Taylor–Hood element pair, satisfiesthe LBB condition [82], whereas the element combination does not.We refer to Donea& Huerta [18] and Gresho&Sani [32] for an extensive discussionof LBB-stable elements. However, an element pair that does not satisfythe LBB condition (2.7) can yet be modified so that stability of discretizationis obtained. Papers I and II investigate flow problems with the element pair. For stabilization, an artificial diffusion on the pressure space is introducedin a systematic way, so that consistency is preserved. This approach [38], [45], [84] modifies the test functions in the discrete variational form (2.6), resulting in a Petrov–Galerkin type scheme. Similar stabilization approaches have been developed to avoid unphysical oscillations in velocity that otherwise appear in convection-dominated regimes for the standard central difference approximations introduced by the Galerkin FEM [18], [31], [32]. We refer to Braack et al. [10] for a review of various commonly used stabilization techniques.
2.3 TIME DISCRETIZATION
For discretization in time, it is usual to use some standard procedure developed for ordinary differential equations and differential–algebraic equations, see Donea& Huerta [18] and Hairer et al. [36], [37] for derivation and analysis of numerical schemes. For schemes that treat both convection and diffusion explicitly, linear stability theory sets a restrictive limit on the time step ¢t in terms of the spatial mesh size, namely
See, e.g., the discussion in [52]. For applications where the viscosity is not too small compared to the velocity, (semi-)implicit time stepping schemes are often preferred in order to avoid the diffusive stability constraint. Popular schemes are variants of the one-step theta method [51] and the BDF-2 method. For a coupled solution approach, stability requires implicit time discretization of the (differential–algebraic) terms . These terms enforce the divergence-free condition on the velocity in Lagrangian multiplier form, see, e.g., Gunzburger [35] and Lions [62].
Harlow and Welch [39] proposed to implement time propagation of The Navier–Stokes equations by fractional time stepping strategy, and which was then recast as a projection scheme by Chorin [15] and Temam [83]. This concept has been mathematically justified for the finite element method by Guermond and Quartapelle [33], [34]. Projection schemes first advance velocities subject to the momentum equation with pressure extrapolated from the old time step. The resulting velocity is in general not divergence-free, so a pressure Poisson equation is solved with forcing given by the divergence of the medium-step velocity in a second step, in order to correct the velocity. Projection schemes are to date the most popular variant for the numerical approximation of the Navier–Stokes equations [74]. The subsystems resulting from fractional time stepping are of easier structure than the original saddle-point system. However, each of the subproblems needs to be equipped with suitable boundary conditions for well-posedness, which are usually not the physical ones. This can yield a non-consistent approximation close to the boundary. We refer to Gresho&Sani [32] and Quartapelle [74] for discussion, and to the articles [13] and [52] for recent developments.
The work presented in this thesis applies an alternative approach where the system remains coupled. Coupled schemes solve for pressure and velocity as one algebraic system, with boundary conditions set according to physical reasoning. However, the algebraic system is in general more challenging because of its saddle point nature, and puts high demands on numerical linear algebra
2.4 NUMERICAL LINEAR ALGEBRA FOR COUPLED DISCRETIZATIONS
The introduction of spatial and temporal discretization schemes for the Navier–Stokes system (2.6) results in a nonlinear system to be solved. For resolving the non-linearity, fixed point methods or Newton-type iterations are usually applied, finally yielding a system of linear equations. The Lagrangian basis functions are localized, so the final coefficient matrix is sparse. For the solution of the saddle-point Navier–Stokes system, different techniques have been proposed. Efficient schemes employ multigrid methods with suitably defined smoothing schemes (see Trottenberg, Oosterlee&Schüller [85] and John [48]), or use preconditioners based on the Schur complement of the block matrix system defined by velocity and pressure matrices (cf. [55], [53], [22], see also the books by Turek [88] and Elman, Silvester&Wathen [23]).
3. A MULTISCALE APPROACH TO LARGE EDDY SIMULATION
The definition of a LES needs to distinguish between resolved and unresolved scales in order to formulate the subgrid viscosity term. Traditional methods use low-pass filters to extract the large-scale information from the Navier–Stokes equations and introduce the subgrid viscosity [77]. However, this spatial filtering is problematic close to domain boundaries. These problems can be traced back to so-called commutation errors [7], [20]. Moreover, finding appropriate boundary conditions for the averaged large scales is to date an open problem. An alternative approach to LES is based on the variation multiscale method [44], introduced by Hughes, Mazzei, and Jansen [46]. The ideas are closely related to the concept of dynamic multilevel methods, described in Dubois, Jauberteau&Temam [19], and the multiscale turbulence analysis in Collis [16]. In this framework, filters are replaced by variation projections onto appropriate subspaces. The space for admissible velocity solutions is decomposed into
u
is a space associated with mesh size 3h, representing large resolvedu, representing small resolved scales, and finally, the space of unresolved scales . We choose a resolution of 3h, however, different definitions of large resolved scales are common as well, see the review article by Gravemeier [30]. Applying an equivalent decomposition for pressure, we get (following Harten’s notation [40])
the form that contains all the left-hand-side terms of the weak form of the Navier–Stokes equations (2.6). The aim of any numerical method is to find numerical approximations uh, phin the resolved subspace V h u Vp according to the direct sum decomposition (3.3). A systematic way of doing so is to insert the decomposition (3.4) into the weak form (2.6), and to restrict the test functions to the space of resolved variables. Then, the nonlinear convective term is linearized and resorted in terms of resolved and unresolved quantities [46]. Hence, we search for Contains linear terms in (uˆ ,pˆ) and the linearized convection around uh, and the term B2 Contains the quadratic convection term. For the numerical approximation to be computable, these two terms need to be substituted by some modeling that eliminates the unresolved variables.
Note that the sub grid viscosity is only applied to the smaller of the resolved scales, as discussed, e.g., in the review articles by Grave Meier [30] and John [49]. A physical explanation for application of eddy viscosity only on small resolved scales is that energy transport occurs mainly on neighboring scale groups, so that the effect of a not adequately resolved viscosity is best modeled on the fine resolved scales [66].The implementation of a method according to the framework (3.5) with approximation (3.6) needs to distinguish between coarse and fine resolved scales. Standard methods for doing this are based on different polynomial orders for the FE basis functions, or two different grid levels as used, e.g., in geometric multigrid. Such methods have been presented by John & Kaya [50] and Grave Meier [29]. Papers discuss an implementation where this separation is performed on the algebraic equation system by using routines from algebraic multigrid.
4. SIMULATION OF TWO-PHASE FLOW
We consider the dynamics of two immiscible incompressible fluids separated by an interface ¡ as displayed in Fig. 4.1. By 1 we denote the region occupied by the first fluid, and by 2 the region occupied by the second. We assume that the flow is laminar, i.e., the Reynolds number is low for the considerations in this chapter.
Figure 4.1:
Schematic view of a two-phase flow problem. Fluid 1 occupies region 1, and fluid 2 occupies 2. The interface separates the two fluids.
The motion of each fluid is described by the Navier–Stokes equations (2.1) – (2.2) using the variables and 1, 2. Extra conditions are necessary to describe the behavior on the interface. A widespread formulation is to pose the two-phase flow problem on the whole domain with variables u and p that take the respective velocity and pressure values in the individual domains. Generally, the fluid densities ½i and viscosities ¹i are different for the two fluids. Hence, the material parameters ½ and ¹ in the global momentum equation have a discontinuity on the interfaceHere, ¾ is a constant specifying the relative strength of the surface tension, · denotes the interface curvature, nthe direction normal to the interface (pointing into 2), and ±¡ is a delta function that localizes the force to the interface. The interface ¡ is transported by the local fluid motion The variables in a single-domain formulation comprise kinks and discontinuities due to changes in material parameters and forces: The discontinuous viscosity introduces a discontinuity in the symmetric gradient "(u). Similarly, due to the distributional character of the delta function ±¡ there is a jump in pressure proportional to ¾·. Any change in curvature tangential to the interface gives rise to a discontinuity in the pressure derivative normal to the interface as well as the viscous stress tensor. The conditions describing the irregularities in velocity and pressure over the interface are called jump conditions, see Edwards, Brenner &Wasan [21] and Lee &LeVeque [59]. Because of different fluid densities, gravitational forces are often essential for the dynamics of two-phase flow
4.1 OVERVIEW OF NUMERICAL MODELS
The two main challenges in the numerical simulation of two-phase flow are the representation of the interface as it evolves in time, and the evaluation of the distributional force term (4.1). Several discrete approaches have been proposed during the last decades. There are two main strategies to couple the interface evolution problem to the Navier–Stokes equations discretized on a fixed grid
The first strategy is to explicitly track the interface position by introducing an additional discrete set of marker points or marker elements that describe the interface. These methods are called front tracking and were introduced in different variants by Peskin [71] and Glimm et al. [27], and reviewed in the articles [78], [86], and [89]. A special class of front tracking schemes are the immersed boundary method by Peskin (see Peskin’s [72] review paper) and a more accurate variant, the immersed interface method, by LeVeque and co-workers [59], [60], [61]. The evolution of the interface in front tracking schemes is accomplished by Lagrangian advection of the marker points, and the surface tension force is generally evaluated based on a discrete approximation of the delta function around the interface. Normal vectors and curvatures can be calculated using the interface points by straight-forward geometric identities. Explicit descriptions are very powerful because they allow to include general models of the interface, like elastic forces in immersed boundary and interface implementation [72]. However, the methods are quite difficult to implement because the marker points need to interact with the fluid grid in order to describe advection and interface forces. Straight-forward implementations are also prone to unphysical changes of the area/volume of the respective fluid. To retain an accurate interface representation, the marker points need to be redistributed when the interface is deformed. The second strategy is front-capturing, where the interface is implicitly defined. Historically, the first scheme of this kind was the marker-and-cell method proposed by Harlow and Welch [39], and can be considered to be a volume tracking method. Here, one fluid is colored by marker particles whose location is adverted with the flow field. The position of the interface can then be reconstructed from this particle field. An extension of this approach is to replace the marker particles by a marker function. This approach has been pursued by Noh & Woodward [65] and Hirt& Nichols [42], who introduced the volume-of-fluid (VOF) method. VOF includes an additional variable that stores the fraction of the first fluid on the total fluid for each cell. From the volume fractions, the position of the interface can be reconstructed, e.g., by defining line segments. For details of the reconstruction, including higher order schemes, we refer to Scardovelli&Zaleski [78] and Sethian [79], and references therein. The advection of the interface is usually implemented by increasing or decreasing the volume fraction depending on the composition of neighboring cells and the velocity field. To put surface tension force in VOF into practice, several different approaches are common. One is to use discrete delta functions in combination with the reconstructed interface, and another employs the continuous surface tension approach by Brackbill, Kothe, and Zemach [11], see also the review by Tryggvason et al. [86]. The advantage of VOFmethods is the volume conservation and a natural mechanism for breakup and fusion of bubbles. However, the reconstruction of the interface and the interface motion are considerably more complicated to implement here than for front tracking methods. While VOF methods operate on the discrete level, a continuous front capturing framework is provided by the level set method, introduced by Osher&Sethian [70] and first applied to two-phase incompressible flow by Sussman, Smereka&Osher [81]. Yet another method is provided by the phase field method, formulated by Jacqmin [47] for the simulation of two-phase flow. The methods discussed are based on the latter two approaches, which are discussed in detail in the following two sections.
4.2 LEVEL SET MODELS
The basic idea of level set methods is to define the interface implicitly as the zero contour line of a function Á that is defined in . Level set methods allow for a straight-forward evaluation of normal vectors by setting and curvature. Both these quantities can be computed locally from on the same computational mesh as the equations of fluid flow. The mechanism for moving the interface in level set methods is advection of the function with local fluid velocity.
This is an additional partial differential equation coupled to the Navier–Stokes equations. The advection can be solved with standard tools for hyperbolic transport equations, like upwind finite difference/volume schemes or stabilized finite element methods. The advection can be implemented so that the volume is conserved, since the velocity is divergence-free
Figure 4.2:
Representing an interface by level set functions Á.
The most popular choice for a level set function is the signed distance function Á, depicted in Fig. 4.2(a). We refer to the books by Sethian [79] as well as Osher and Fedkiw [69] for an extensive discussion of thesemethods. A signed distance function encodes the distance of a particular point to the interface, and different signs identify the two fluids. A signed distance function has the property almost everywhere. A benefit of this choice is that the integral in the surface tension force can be simplified to a one-dimensional discrete delta function as Z
wherex denotes an arbitrary point in the computational domain and y(s) is the image of a parametrization of the interface, see, e.g., [79]. However, it has been shown by Engquist, Tornberg, and Tsai [24] that such an approximation can lead errors in the force unless the width h of the delta function is varied. Recent implementations therefore reconstruct the interface position from Á and apply a d-dimensional discrete delta function that does not show these consistency problems. Since the information from Á is only needed in a region around the interface, so-called narrow-band level set implementations are widely used. These methods calculate Á only on a few grid cells around the interface, often in combination with high-resolution grids. We refer to Adalsteinsson and Sethian [1] for details. To start a level set calculation, an initial profile needs to generated given a certain interface position. Moreover, the flow field as well as inaccuracies in the numerical scheme deform the signed distance function in the course of the simulation. Therefore, algorithms have been proposed to (re-)initialize the signed distance function. These algorithms enforce the property for example, by solving the partial differential equation
to steady state, where S(Á0) is a sign function with value 1 in the first fluid and value ¡1 in the other (see, e.g. [69]). However, discretization of the reinitialization scheme do not preserve mass, so that the overall level set implementation based on signed distance functionsmay suffer from unphysical volume changes in the fluid phases.
Another choice is to use a smoothed color/indicator function as depicted in Fig. 4.2(b). Such a function is (almost) constantly 1 away from the interfaces, with positive value in one region and negative value in the other. Around the interface, there is a transition region where the function smoothly switches from one value to the other. This function Á can be calculated from the reinitialization equation (cf. Olsson and co-workers [67], [68])
The factor two comes fromthe fact that the function switches from¡1 toÅ1. This scheme for evaluating the interface force corresponds to the continuous surface tension model proposed by Brackbill, Kothe, and Zemach [11]. A disadvantage of the formulation with a smoothed color function is that the evaluation of curvature is potentially less accurate because it is calculated from a steep profile. This generally increases resolution requirements compared to the signed distance function approach. Another drawback of the smoothed indicator function is that the function needs to be defined on the whole domain in order to conserve volume.
4.3 PHASE FIELD MODELS
All the methods discussed above are based on the mathematical model for interfacial tension, assuming a continuous PDE model for fluid dynamics that does not take processes on the level of atoms and molecules into account. The numerical approximations aim at discretely reproducing the effect of a delta function, and evaluating the curvature to determine the strength of the surface tension. A different approach takes chemical considerations of the interface into account. The so-called van der Waals hypothesis [90] states that immiscible fluids actually mix on molecular level, forming a diffuse interface. The profile of the interface is given by a balance of random molecular motion and molecular attraction, as devised by Cahn and Hilliard [14]. The considerations are based on the free energy of a molecule, given.
WhereC denotes the concentration and ª(C) Æ (C Å1)2(C ¡1)2. Similar to the smoothed color function (4.6), the function seeks to attain the two equilibrium values C Æ §1, which is driven by the double-well potential ª(C). The gradient term on the other hand enforces a smooth transition. The system seeks to minimize the energy.
Which is the same form as for the continuous surface tension model, that is, equation multiplied by curvature ·. Boyer [9] extended the phase field model to describe motion of fluids with different densities and viscosities. The phase field model can be augmented in order to include contact forces of the interface with solid walls [47], which results in flux-boundary conditions on the chemical potential Ã. This approach allows for setting static contact angles that correspond to physical wetting properties. The direct inclusion of contact line dynamics is different from the other approaches in computational multiphase dynamics, where additional mechanisms as slip-velocities and near-wall interface diffusion need to be added, cf. [76], [92].
A considerable challenge in the simulation of phase field models is the high resolution requirements. Even though the interface width can be chosen considerably larger than molecular dynamics would suggest (see Villanueva and Amberg [91]), accurate results require sophisticated numerical approximations and fine computational meshes. Furthermore, the equations are nonlinear and contain fourth order derivatives. Several attempts have been made to make this approach more competitive, see for example the work by Kay and Welford [54].
5. CONCLUSIONS
The method proposed in this paper aims at providing phase-field-like features without having to use all its terms everywhere. This is done so that the phase field method is reduced to a level set formulation in regions far away from contact lines. It will be necessary to perform an in-depth convergence study of level set models and phase field models on several test cases in two-phase flow dynamics. In particular, the resolution requirements of the two methods should be quantified in order to identify the precise benefits of the hybrid method. In our hybrid model, the equations have been coupled by a smoothly varying switch function. The influence of the switch function needs to be considered in an energy estimate in order to ensure stability of the coupling of the equations, which has up to date only been verified numerically. Since the numerical implementation of the equations is based on a system of two equations to avoid the direct appearance of fourth derivatives, the analysis needs to be performed in a way to correctly reflect the discretized state of the equations.
We also aim at using our models for larger problems, especially three-dimensional simulations. For this purpose, parallel implementations of the flow solver and the coupled level-set/Navier–Stokes system are necessary. Within the software that was used, it is possible to use parallel assembly and linear algebra routines. For a Stokes problem coupled to an advection equation, a parallel version for up to about 50 processors has already been implemented by the author of this paper. The objective is to apply this implementation framework to the coupled Navier–Stokes/level-set system and to eventually extend it to a massively parallel implementation for even larger systems. The hybrid level-set/phase-field method is a potential tool for improving the simulation of multi-phase flow in porous media [17]. A first step in that direction is to apply the new model to more complicated geometries, like, e.g., channels with oscillating walls or flow around small obstacles. Small scale results are a helpful tool in the simulation of subsurface flow like the prediction of the flow in oil and gas reservoirs. The imposition of a pressure gradient on a small-scale material configuration induces a flow field, which can be used to determine values for the permeability of that configuration. The permeability tensor is one of the main input parameter in coarse-scale simulations based on Darcy’s law. In this, a study regarding error sources in the discrete approximation with the conservative level set implementation is presented. These results point out directions where the method needs to be improved. One such direction concerns the accuracy in the curvature calculations. A particularly easy approach is to use higher order finite elements for representing the level set function. When keeping the interface width constant, increasing the order of approximation considerably increases resolution of the interface, and thus, accuracy of normal and curvature that rely on the interface representation. This improvement should be advantageous for the overall accuracy, even though a small imbalance in the force representation is introduced. Such a strategy increases the costs for the level set portion of the algorithm if the same mesh is used as for the Navier– Stokes equations, though. An alternative would be to coarsen the level set mesh in regions away from the interface. This is reasonable since in these remote regions, only the momentum and continuity balances need to be represented accurately, while the level set equation does not contain any additional information. Similarly, the interface region can be resolved by more mesh points than there are used for the Navier–Stokes equations. However, this requires a special implementation for the evaluation of the surface tension force in the finite element algorithms because the variables are related to different meshes. An alternative hybrid approach that might be taken into account is to not mix the level set formulation with the phase field formulation as has been done in this paper, but to use the phase field model as a micro-scale input to a macro-scale formulation based on a level set description of two-phase flow. A heterogeneous multiscale approach has been proposed by Ren and E [75] for modelling the interaction of macro-scale fluid flow by micro-scale information from molecular dynamics. The response on molecular level provides information about the expected motion of contact lines during one time step. This coupled model behaves differently than large-scale level set model, which would assume zero velocity. It is conceivable that a similar interaction could be applied for a micro scale simulation done with the phase field model. The difficulty of such a micro macro coupling is to find appropriate transfer operations from the micro scales to the macro scales. Suitable boundary conditions for the Cahn–Hilliard sub problem given the solution of a macro-scale level set model are needed, as well as a shear velocity for the macro-scale Navier–Stokes equations, based on the results from a micro model.
REFERENCES:
[1] D. Adalsteinsson and J. A. Sethian. A Fast Level Set Method for Propagating Interfaces.J. Comput. Phys., 118(2):269–277, 1995.
[2] D. M. Anderson, G. B. McFadden, and A. A. Wheeler. Diffuse-Interface Methods in FluidMechanics.Annu. Rev. FluidMech., 30:139–165, 1998.
[3] E. Aulisa, S. Manservisi, and R. Scardovelli.A novel representation of the surface tension force for two-phase flow with reduced spurious currents.Comput.MethodsAppl.Mech. Engrg., 195:6239–6257, 2006.
[4] I. Babuška. Error-bounds for finite element method.Numer. Math., 16(4):322–333, 1971.
[5] W. Bangerth, R. Hartmann, and Kanschat G. deal.II — a General Purpose Object Oriented Finite Element Library.ACMTrans.Math. Softw., 33(4):article no. 24, 2007.
[6] Y. Bazilevs, V. M. Calo, J. A. Cottrell, T. J. R. Hughes, A. Reali, and G. Scovazzi.Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows.Comput.MethodsAppl.Mech. Engrg., 197:173–201, 2007.
[7] L. Berselli, C. Grisanti, and V. John.On Communtation Errors in the Derivation of Space Averaged Navier–Stokes Equations.Technical Report 12, Dipartimento di MatematicaApplicata, Università di Pisa, 2004.
[8] R. P. Beyer and R. J. Leveque.Analysis of a one-dimensional model for the immersed boundary method.SIAM J. Numer. Anal., 29(2):332–364, 1992.
[9] F. Boyer. A theoretical and numerical model for the study of incompressible mixture flows. Comput.& Fluids, 31:41–68, 2002.
[10] M. Braack, E. Burman, V. John, and G. Lube.Stabilized finite element methods for the generalizedOseen problem.Comput.MethodsAppl.Mech. Engrg., 196:853–866, 2007.
[11] J. U. Brackbill, D. B. Kothe, and C. Zemach. A Continuum Method for Modeling Surface Tension.J. Comput. Phys., 100:335–354, 1992.
[12] F. Brezzi. On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers.Rev. Française Automat. Informat.RechercheOpérationelleSér. Rouge, 8(R-2):129–151, 1974.
[13] D. L. Brown, R. Cortez, and M. L. Minion. Accurate Projection Methods for the Incompressible Navier–Stokes Equations.J. Comput. Phys., 168(2):464–499, 2001.
[14] J.W. Cahn and J. E. Hilliard.Free Energy of a NonuniformSystem. I. Interfacial Free Energy. J. Chem. Phys., 28(2):258–267, 1958.
[15] A. J. Chorin. Numerical Solution of the Navier–Stokes Equations.Math.Comput., 22:745–762, 1968.
[16] S. S. Collis. Monitoring unresolved scales in multiscale turbulencemodeling.Phys. Fluids, 13:1800–1806, 2001.
[17] C. Crowe. Multiphase Flow Handbook. CRS Press, Boca Raton, 2006.
[18] J. Donea and A. Huerta.Finite ElementMethods for Flow Problems.J.Wiley& Sons, Chichester, 2003.
[19] T. Dubois, F. Jauberteau, and R. Temam.Dynamic Multilevel Methods and the Numerical Simulation of Turbulence. Cambridge University Press, Cambridge, 1999.
[20] A. Dunca, V. John, andW. J. Layton. The Commutation Error of the Space Averaged Navier–Stokes Equations on a Bounded Domain. In G. P. Galdi, J. G. Heywood, and R. Rannacher, editors, Contributions to Current Challenges in Mathematical FluidMechanics, volume 3, pages 53–78. BirkhäuserVerlag, Basel, 2004.
[21] D. A. Edwards, H. Brenner, and D. T. Wasan. Interfacial Transport Processes and Rheology.Butterworth–Heinemann, Stoneham, 1991.
[22] H. Elman, D. Silvester, and A. Wathen.Performance And Analysis of Saddle Point Preconditioners for The Discrete Steady-State Navier–Stokes Equations. Numer.Math., 90(4):665–688, 2002.
[23] H. Elman, D. Silvester, and A. Wathen.Finite Elements and Fast Iterative Solvers with Applications in Incompressible Fluid Dynamics.Oxford Science Publications,Oxford, 2005.
[24] B. Engquist, A.-K.Tornberg, and R. Tsai. Discretization of Dirac delta functions in level set methods. J. Comput. Phys., 207:28–51, 2005.
[25] M. M. Francois, S. J. Cummins, E. D. Dendy, D. B. Kothe, J. M. Sicilian, and M. W. Williams. A balanced-force algorithm for continuous and sharp interface surface tension models within a volume tracking framework. J. Comput. Phys., 213:141– 173, 2006.
[26] M. W. Gee, C. M. Siefert, J. J. Hu, R. S. Tuminaro, and M. G. Sala. ML 5.0 Smoothed Aggregation User’s Guide.Technical Report 2006-2649, Sandia National Laboratories, 2006.
[27] J. Glimm, J. Grove, B. Lendquist, O. A.McBryan, and G. Tryggvason.The bifurcation of tracked scalar waves.SIAM J. Sci. Stat. Comput., 9(1):61–79, 1988.
[28] R. Glowinski. Finite element methods for incompressible viscous flow. In J. L. Lions and P. G. Ciarlet, editors, NumericalMethods for Fluids (3), Handbook of NumericalAnalysis, volume 9, pages 3–1176. North-Holland, Amsterdam, 2003.
[29] V. Gravemeier. Scale-separating operators for variationalmultiscale large eddy simulation of turbulent flows. J. Comput. Phys., 212:400–435, 2006.
[30] V. Gravemeier. The VariationalMultiscale Method for Laminar and Turbulent Flow.Arch. Comput.Meth. Engrg., 13(2):249–324, 2006.
[31] P. M. Gresho and R. L. Sani.Incompressible Flow and the Finite Element Method, volume one: Advection–Diffusion. JohnWiley& Sons, Chichester, 2000.
[32] P. M. Gresho and R. L. Sani.Incompressible Flow and the Finite Element Method, volume two: Isothermal Laminar Flow. JohnWiley& Sons, Chichester, 2000.
[33] J. L. Guermond and L. Quartapelle.Calculation of Incompressible Viscous Flow by an Unconditionally Stable Projection FEM. J. Comput. Phys., 132:12–33, 1997.
[34] J. L. Guermond and L. Quartapelle. A projection FEM for variable density incompressible flows. J. Comput. Phys., 165(1):167–188, 2000.
[35] M. D. Gunzburger. Finite Element Methods for Viscous Incompressible Flows. Academic Press, Boston, 1989.
[36] E. Hairer, S. P. Nørsett, and G. Wanner.Solving Ordinary Differential Equations I. Nonstiff Problems. Springer-Verlag, Berlin, 2nd edition, 1993.
[37] E. Hairer and G. Wanner.Solving Ordinary Differential Equations II.Stiff and Differential-Algebraic Problems. Springer-Verlag, Berlin, 1991.
[38] P. Hansbo and A. Szepessy. A velocity–pressure streamline diffusion FEM. Comput. Methods Appl.Mech.Engrg., 84(2):175–192, 1990.
[39] F. H. Harlow and J. E. Welch. Numerical Calculation of Time-Dependent Viscous Incompressible FlowofFluidwith Free Surface. Phys. Fluids, 8(12):2182–2189, 1965. [40] A. Harten. Multiresolution Representation of Data: A General Framework. SIAM J.Numer. Anal., 33(3):1205–1256, 1996.
[41] M. A. Heroux, R. A. Bartlett, V. E. Howle, R. J. Hoekstra, J. J. Hu, T. G. Kolda, R. B. Lehoucq, K. R. Long, R. P. Pawlowski, E. T. Phipps, A. G. Salinger, H. K. Thornquist, R. S. Tuminaro, J.M.Willenbring,Williams A., and K. S. Stanley. An overview of the Trilinos project.ACM Trans.Math.Softw., 31(3):397–423, 2005.
[42] C.W. Hirt and B. D. Nichols.Volume of fluid (VOF) method for the dynamics of free boundaries.J. Comput. Phys., 39(1):201–225, 1981.
[43] T. J. R. Hughes. The finite element method: linear static and dynamic finite element analysis. Dover,Mineola, NY, 2000.
[44] T. J. R. Hughes, G. R. Feijóo, L. Mazzei, and J.-B. Quincy. The variational multiscale method – a paradigm for computational mechanics.Comput.MethodsAppl.Mech.Engrg., 166:3–24, 1998.
[45] T. J. R. Hughes, L. P. Franca, and M. Balestra. A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuska–Brezzi condition: A stable Petrov–Galerkin formulation of the Stokes problem accomodatingequalorder interpolations. Comput.MethodsAppl.Mech. Engrg., 59:85–99, 1986.
[46] T. J. R. Hughes, L. Mazzei, and K. E. Jansen. Large Eddy Simulation and the variationalmultiscale
method. Comput.Visual. Sci., 3:47–59, 2000.
[47] D. Jacqmin. Calculation of two-phase Navier–Stokes flows using phase-field modeling. J. Comput. Phys., 155(1):96–127, 1999.
[48] V. John. On the efficiency of linearization schemes and coupled multigrid methods in the simulation of a 3D flow around a cylinder.Int. J. Numer.Meth. Fluids, 50:845– 862, 2006.
[49] V. John. On large eddy simulation and variational multiscale methods in the numerical simulation of turbulent incompressible flows.Appl. Math., 50(4):321–353, 2008.
[50] V. John and S. Kaya. A finite element variational multiscale method for the Navier– Stokes equations.SIAM J. Sci. Comput., 26(5):1485–1503, 2005.
[51] V. John, G. Matthies, and J. Rang. A comparison of time-discretization/linearization approaches for the incompressibleNavier–Stokes equations.Comput.MethodsAppl.Mech. Engrg., 195(44-47):5995–6010, 2006.
[52] H. Johnston and J.-G.Liu. Accurate, stable and efficientNavier–Stokes solvers based on explicit treatment of the pressure term. J. Comput. Phys., 199(1):221–259, 2004.
[53] D. A. Kay, D. Loghin, and A. Wathen. A preconditioner for the steady-state Navier– Stokes equations.SIAM J. Sci. Comput., 24(1):237–256, 2002.
[54] D. A. Kay and R. R.Welford. Efficient Numerical Solution of Cahn–Hilliard–Navier–Stokes Fluids in 2D.SIAM J. Sci. Comput., 29(6):2241–2257, 2007.
[55] A. Klawonn. An Optimal Preconditioner for a Class of Saddle Point Problems with a Penalty Term.SIAM J. Sci. Comput., 19(2):540–552, 1998.
[56] A. N. Kolmogorov. The Local Structure of Turbulence in Incompressible Flows (in Russian).Proceedings of the USSR Academy of Sciences, 30:299–303, 1941.
[57] D. B. Kothe,W. J. Rider, S. J.Mosso, J. S. Brock, and J. I. Hochstein. Volume tracking of interfaces having surface tension in two and three dimensions. In Proceedingsof the 34th Aerospace SciencesMeeting and Exhibit, Technical Report AIAA 96-0859, Reno, NV, 1996.
[58] O. A. Ladyzhenskaya. The mathematical theory of viscous incompressible flow. Gordon
and Breach Science, New York, 1969.
[59] L. Lee and R. LeVeque. An immersed interface method for incompressible Navier– Stokes equations.SIAM J. Sci. Comput., 25(3):832–856, 2003.
[60] R. J. LeVeque and Z. Li.The immersed interface method for elliptic equations with discontinuous coefficients and singular sources.SIAM J. Numer. Anal., 31(4):1019– 1044, 1994.
[61] R. J. LeVeque and Z. Li. Immersed interface method for Stokes flow with elastic boundaries or surface tension. SIAM J. Sci. Comput., 18(3):709–735, 1997.
[62] P.-L. Lions.Mathematical Topics in Fluid Mechanics, volume 1: Incompressible Models. Oxford University Press, Oxford, 1996.
[63] P.Moin and K.Mahesh. Direct numerical simulation: a tool in turbulence research. Annu. Rev. FluidMech., 30:539–578, 1998.
[64] R. D. Moser, J. Kim, and N. N. Mansour. Direct numerical simulation of turbulent channel flow up to Re¿ Æ 590.Phys. Fluids, 11(4):943–945, 1999.
[65] W. Noh and P. Woodwards. SLIC (simple line interface calculation). In Proceedings of the 5th International Conference onNumericalMethods in FluidDynamics, pages330–340, Enschede, 1976. Springer-Verlag.
[66] A. A. Oberai, V. Gravemeier, and G. C. Burton.Transfer of Energy in the Variational Multiscale Formulation of LES.Technical report, Center for Turbulence Research, Stanford, Summer Program 2004.
[67] E. Olsson and G. Kreiss. A conservative level set method for two phase flow. J. Comput. Phys., 210:225–246, 2005.
[68] E. Olsson, G. Kreiss, and S. Zahedi. A conservative level set method for two phase flow II. J. Comput. Phys., 225:785–807, 2007.
[69] S.Osher and R. P. Fedkiw. Level SetMethodsandDynamic Implicit Surfaces. Applied Mathematical Sciences.Springer, New York, 2003.
[70] S. Osher and J. A. Sethian. Fronts Propagating with Curvature-Dependent Speed: Algorithms Based on Hamilton–Jacobi Formulations. J. Comput. Phys., 79(1):12– 49, 1988.
[71] C. S. Peskin. Numerical Analysis of Blood Flow in the Heart.J. Comput. Phys., 25(3):220–252, 1977.
[72] C. S. Peskin. The immersed boundary method.ActaNumerica, 11:479–517, 2002.
[73] S. B. Pope. Turbulent Flows. Cambridge University Press, Cambridge, 2000.
[74] L. Quartapelle. Numerical solution of the incompressible Navier–Stokes equations, volume 113 of International Series of Numerical Mathematics.Birkhäuser-Verlag, Basel, 1993.
[75] W. Ren and W. E. Heterogeneous multiscale method for the modeling of complex fluids and micro-fluidics.J. Comput. Phys., 204(1):1–26, 2005.
[76] M. Renardy, Y. Renardy, and J. Li. Numerical Simulation of Moving Contact Line ProblemsUsing a Volume-of-FluidMethod.J. Comput. Phys., 171(1):243–263, 2001.
[77] P. Sagaut. Large Eddy Simulation for Incompressible Flows. Springer-Verlag, Berlin, 3rd edition, 2006.
[78] R. Scardovelli and S. Zaleski. Direct Numerical Simulation of Free-Surface and Interfacial Flow. Ann. Rev. FluidMech., 31:567–603, 1999.
[79] J. A. Sethian. Level Set Methods and Fast Marching Methods. Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science.Cambridge University Press, 2000.
[80] J. Smagorinksy. General Circulation Experiments with the Primitive Equations. I. The Basic Experiment. MonthlyWeather Rev., 91(3):99–164, 1963.
[81] M. Sussman, P. Smereka, and S. Osher. A level set method for computing solutions to incompressible two-phase flow. J. Comput. Phys., 114(1):146–159, 1994.
[82] C. Taylor and P. Hood. A numerical solution of the Navier–Stokes equations using the finite element technique.Comput. Fluids, 1(1):73–100, 1973.
[83] R. Teman. Uneméthoded’approximation de la solution des équations de Navier–
Stokes.Bull. Soc.Math. France, 96:115–152, 1968.
[84] T. E. Tezduyar. Stabilized Finite Element Formulations for Incompressible Flow Computations.Adv. Appl.Mech., 28:1–44, 1992.
[85] U. Trottenberg, C. Oosterlee, and A. Schüller. Multigrid. Elsevier Academic Press, London, 2001.
[86] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, and Y.-J. Jan. A Front Tracking Method for the Computations of Multiphase Flow.J. Comput. Phys., 169(2):708–759, 2001.
[87] R. Tuminaro and C. Tong. Parallel smoothed aggregation multigrid: aggregation
strategiesonmassively parallel machines. In J. Donnelley, editor, Super Computing
2000 Proceedings
, 2000.
[88] S. Turek. Efficient Solvers for Incompressible Flow Problems: An Algorithmic and Computational Approach. Springer, Berlin, 1999.
[89] S.O.UnverdiandG. Tryggvason.A Front Tracking Method for Viscous, Incompressible, Multi-fluid Flows.J. Comput. Phys., 100(1):25–37, 1992.
[90] J. D. van derWaals. The thermodynamic theory of capillarity under the hypothesis of continuous variation of density.J. Stat. Phys. (originally published in Dutch inVerhandel. Konink.Akad.Weten. Amsterdam, 1893), 20(2):200–244, 1979.
[91] W. Villanueva and G. Amberg.Some generic capillary-driven flows.Int. J. Multiphase Flow, 32:1072–1086, 2006.
[92] S. Zahedi, K. Gustavsson, and G. Kreiss. A conservative level set method for contact line dynamics. J. Comput. Phys., 228(17):6361–6375, 2009.
