CFD|多相流之水平集ODDLS原理(中)

HOW TO CALCULATE A SIGNED DISTANCE TO THE INTERFACE

   Based on the flow properties, such as density or viscosity, we can calculate an initial guess foras follows:

   Clearly, the level set function obtained from Equation (17) is not a signed distance to. We need to reinitiate in order to guarantee that fulfils definition (5). We use a technique based on the following property to reinitiate the level set function as a signed distance to [9]:

   If is a signed distance function to , then

where denotes the Euclidian norm in 

    We use Equation (18) to recalculate. In case Equation (18) is not satisfied, then

    The residual r in Equation (19) can be understood as the time derivate of with respect to a pseudo-time τ, i.e.

   Then, the stationary solutions to Equation (20) (i.e. when =0) satisfy Equation (18). In summary, for a given time t ∈[0,T], we calculate as the stationary solution to the following hyperbolic problem:

Problem (21) can be rewritten in a more convenient manner as

   As the level set values identify the free surface between the two fluids, the following relations can be obtained

    where n is the normal vector to the interface , oriented from fluid 2 to fluid 1, and k is the local curvature of the interface.

    From Equation (22) it can be easily shown that the reinitialization process begins from the interface and it propagates to the whole domain. This property is attractive from the computational point of view, as we are interested only in the accurate values of in the region close to the interface. Then we only need to reach the steady-state solution to Equation (22) in a narrow band close to the interface as proposed in [7]. The reinitialization algorithm used in this work has been optimized using the concept of nodal levels described next.

   Given a level set distribution , for every time step t one can assign a level  to the node  based on the following rules:

   A nodal point i of coordinates  belongs to the level =1, if (

, t )≥0, and it is connected to at least one node j with (, t)<0.

   A nodal point i  of coordinates  belongs to the level =−1, if ( , t)<0, and it is connected to at least one node j with( , t )≥0.

   A nodal point i of coordinates belongs to the level >1, if ( , t)>0, and it is connected to at least one node of level −1.

   A nodal point i of coordinates belongs to the level <−1, if ( , t)<0, and it is connected to at least one node of level +1.

   Once the levels have been assigned to the mesh points, the reinitialization is performed in a narrow band of n levels around the interface. Equation (22) is solved using a forward Euler scheme until the steady state in the n levels is achieved. k iterations with k>n are performed in such a way that for every iteration i=1, . . .,k. The iterations are carried out only for those nodes j of the mesh with .

   Owing to the evolution of the interface defined by Equation (12), the level set function does not fulfil Equation (18) after some time. It is therefore recommended to apply the reinitialization scheme described above at every time step.

  At a moving interface the jump condition applies. For immiscible fluids we can express the jump condition at the free surface as

where  is the coefficient of surface tension (a constant of the problem), s is any unit vector tangent to the interface and [·] defines the jump across the interface

Taking into account that，i =1, 2, and using that s·n=0, the second equation in (24) yields

  From Equation (26) we can conclude that，where  denotes the orthogonal subspace to v. Then it is clear that

   Integrating over any closed surface  containing , in both sides of Equation (27) yields

   Applying the Gauss theorem to the right-hand side of Equation (28), we can deduce that λ=0 or equivalently

   Introducing Equation (29) into the first condition in Equation (24), the following equivalent condition is obtained

  This form of the jump condition is very useful in the development of the method presented in the following sections.

   In the classical level set formulation for multiphase flow problems [9], jump condition (24) can be introduced into the momentum equation by adding an additional source term, i.e.

   where  is Dirac’s delta function on . As the properties of the fluids change discontinuously across the interface, the direct solution to Equation (31) has important numerical drawbacks: inaccurate resolution of the pressure at the interface, fluid mass losses, inaccurate interface location, etc.

   A common approach to overcome these difficulties is to smooth the discontinuous transition of the flow properties and Dirac’s delta functions across the interface using of a smoothed Heaviside function. This function is usually expressed as

   where 2ε is the thickness of the transition area. Then a property  is approximated as

and Dirac’s delta function is approximated as

   This approach also presents several deficiencies: the accuracy of the results depends on the ad hoc thickness of the transition area 2ε ; the jump condition is not properly captured; the interface propagation velocity is not correctly calculated; etc.

    Other schemes, such as the ghost fluid method [10] or Lagrangian methods such as the PFEM[4],attempt to overcome the problem of smoothing of properties at the interface region. However, these methods have several drawbacks, such as limitations of conservative and convergence properties

in some cases or an increase in the computational effort. As a conclusion, a method able to deal with a discontinuous transition of properties across the interface, the jump interface condition and computationally affordable for large models is required.

    The overlapping domain decomposition method described in the following sections has been developed to overcome these difficulties.

STABILIZED FIC-FEM FORMULATION

  It is well known that the standard Galerkin approximation of the incompressible Navier–Stokes equations may suffer from numerical instabilities from two main sources. The first is due to the advective character of the equations, which induces oscillations for high values of the velocity.

The second source has to do with the mixed character of the equations, which limits the stability of the solution to the satisfaction of the well-known inf–sup condition. The stabilization technique chosen in this work to overcome both problems is based on the FIC method presented in [12].

   The stabilized FIC form of the governing differential equations (11) and (12) is expressed as

The boundary conditions for the stabilized FIC problem are expressed as

 The underlined terms in Equations (35) and (36) introduce the necessary stabilization for the numerical solution to the Navier–Stokes problem [13, 15, 17].

   Note that the terms and  denote the residuals of Equations (11) and (12), i.e.

   The characteristic length vectors  and contain the dimensions of the finite space domain where the balance laws of mechanics are enforced. At the discrete level these dimensions are of the order of the element or grid dimension used for the computation. Details on obtaining the

FIC-stabilized equations and the recommendation for computing the characteristic length vectors can be found in [15–17].

OVERLAPPING DOMAIN DECOMPOSITION METHOD

   To overcome the drawbacks related to the discontinuity of fluid properties across the interface and to impose the interfacial boundary condition we split the domain Ω into two overlapping subdomains. Based on this subdivision of the domain, we apply a Dirichlet–Neumann overlapping

domain decomposition technique with appropriate boundary conditions in order to satisfy the interfacial condition.

  Let  be a finite element partition of domain Ω, where

 are the element nodes, is the element spatial domain,are the element shape functions and K is the total number of elements in the partition.

    We assume that K satisfies the following approximation property:

for a fixed instant t, t ∈(0,T]. Let us consider a domain decomposition of domain Ω into three disjoint subdomains andin such a way thatand  where are the elements of the finite element partition K, such as, and  are the elements of the finite element partition, such asThe geometrical domain decomposition is completed with

From this partition, we define the following two overlapping domains:

Some comments about the partition are given next.

   In order to simplify the notation, we will omit the time dependency of domains in the remainder of this section.

   Let, i =1, 2, be the boundary of , then

   It has to be noted that contains an additional termthat is not included into . This term comes from the presence of an interface inside Ω. As we are using a capturing technique, the position of the interface will not lie on mesh nodes, as usually happens in a tracking technique.

   Therefore, some elements can be intersected by the interface. Then is defined as follows:

Note that is not coincident with , but the following condition is satisfied:

Summarizing, we have

Figure 1 is a sketch of the domain partition described above.

   We have two restrictions for the definition of the boundary condition on : fluid velocities must be compatible at the interface and the jump boundary condition (24) must be satisfied on .

   The domain decomposition technique chosen is of Dirichlet–Neumann type and it allows us to fulfil both restrictions.

   Let us introduce a standard notation. We denote by a variable related to . For example, is the velocity field solution to problem (36), where the domain Ω has been replaced by .

    We apply the Dirichlet conditions on  and the Neumann conditions on . For the boundary  , we use the compatibility of velocities at the interface:

   For the boundary , we use the jump boundary condition. We are looking for a traction vector  such that

and  must guarantee that the jump boundary condition (31) holds, that is to say

   Let us write the variational problem considering the domain decomposition above described.

Domain 1

Domain 2

   An iterative strategy between the two domains is used to reach a converged global solution in the whole domain  Ω. It is expected that this global solution will satisfy both restrictions presented above. For this reason we define the following stop criteria:

The global velocity solution obtained from Equation (49) is used as the convective velocity for the transport of the level set function, i.e.

  Unfortunately there is not a theoretical result that can ensure the convergence of this method.However, our experience in applying this method to a wide range of problems involving moving interfaces has showed that the method is stable and robust.

  An expected property of the method is its dependency with the mesh size around the interface.This issue is directly related with how good is

approximated by , as has been pointed out in property (43). A fine mesh close to the interface significantly reduces the amount of iterations needed to reach a converged global solution, as well as the accuracy of the velocities and the pressure at the interface.

   This methodology can be viewed as a combination of domain decomposition and level set techniques. This justifies the name given to the new method: ODDLS.

DISCRETIZATION BY THE FEM

   In this section we present the final discrete system of equations associated with problem (48) and (49) using the FEM. Let us consider a uniform partition of the time interval of analysis [0,T],with time step size. We will denote by a superscript the time step at which the algorithmic solution is computed. We assume that  are known. If θ∈[0,1], the trapezoidal rule applied to Equations (48) and (49) consists in finding

 as the solution to the problem

where

   This problem is non-linear due to the convective term and the evolution of the free surface.Prior to the finite element discretization, we linearize it using a Picard method, which leads to an Oseen problem for each iteration step. Several strategies can be adopted to deal with the two iterative algorithms involved in the ODDLS method: the overlapping domain decomposition iterative scheme and the linearization scheme (Picard). In our case, for each iteration step of the linearization scheme, the domain decomposition scheme is also updated. This simultaneous iterative strategy requires to complete Equation (50) with a stop criterion for the Picard iteration.The same iteration index is used for both schemes.

    It is noteworthy that the adopted strategy reduces the computational effort compared with that of other nested iterative schemes. We denote by a superscript n the time step at which the algorithmic solution is computed and by a superscript j the iteration step of the domain decomposition scheme(which for a simultaneous iterative strategy will be the same as for the Picard iteration step). The form of the iterative scheme is as follows:

   The Galerkin approximation of Equation (53) is straightforward and stable, thanks to the FICstabilized formulation [13–18]. Equal polynomial order can be used for the approximation of the velocities and the pressure.Based on the finite element partition introduced in Equation (38),letand be two finite dimensional conforming finite element spaces, h being the maximum of the elements’ partition diameters.

Let be the finite element spaces for the test functions. The discretized form of Equations (55) is

Remark

Updating the level set equation is done once the convergence of the Navier–Stokes equations is obtained. Therefore, the complete iteration loop is as follows:

Step 1: Solve Equations (54) in .

Step 2: Solve Equations (54) in .

Step 3: Check convergence on pressure and velocity fields. If it is satisfied, go to Step 4, otherwise go to Step 1.

Step 4: Solve Equation (51).

Step 5: Check convergence on level set equation. If it is satisfied, advance to the next time step,otherwise go to Step 1.

ALE FORMULATION

   It is of interest in many applications to consider the movement of some parts of the analysis domain. In the mobile parts of the domain it is more convenient to use a Lagrangian formulation of the equations and update the spatial discretization every time step. Conversely, in the fixed areas

of the analysis domain, it is more efficient to use the standard Eulerian formulation.

    The ALE formulation of the residuals in (35) can be expressed as follows:

where  is the relative velocity between the local axes fixed to the fluid particle and the global reference axes of the problem [22, 23].