CFD|多相流水平集算法原理(上)

   多相流指的是物质的状态，在现实工程应用中，物质通常具有三相：气相、液相和固相。多相流通常指的是在流动区域内存在两种或两种以上的相。其可以是包含有气体与液体的流动、气体与固体的流动，或者固体与液体的流动，也可以是包含有气液固三相物质的流动。

流体体积模型（ VOF ）

    VOF 模型主要用于跟踪两种或多种不相容流体的界面位置。在VOF模型中，界面跟踪是通过求解相连续方程完成，通过求出体积分量中急剧变化的点来确定分界面的位置。VOF 模型主要应用于分层流、自由液面流动、晃动、液体中存在大气泡的流动、溃坝等现象的仿真计算。其可以计算流动过程中分界面的时空分布。

    DPM模型用于分析颗粒在流体中的运移情况，颗粒可以是固体颗粒，也可以是液滴。DPM模型是在拉格朗日观点下进行的，即在计算过程中是以单个粒子为对象进行计算的，而不像连续相计算那样是在欧拉观点下以空间点为对象。

DPM模型仅适用于颗粒相体积分数小于10%的情况，且不能考虑颗粒体积。不考虑颗粒和颗粒之间的相互作用力，但可以考虑颗粒和流体之间的相互作用。

    水平集方法（Level Set）是一种广泛应用于具有复杂分界面的两相流动问题界面追踪的数值方法。在水平集方法中，分界面提供水平集函数进行捕捉及跟踪。由于水平集函数具有光滑及连续性的特性，其空间梯度能够精确的进行计算，因此可以精确地估算界面曲率及表面张力引起的弯曲曲率。然而水平集方法在保持体积守恒方面存在天然缺陷。水平集方法仅可用于两相流动的区域，且两种流体互不渗透，因此，水平集方法特别适合用于计算开放环境下的流动问题。

    NFX-CFD水平集算法应用ODDLS技术，本文重点介绍水平集模型ODDLS技术，它能以最大程度减少由于流体之间的界面处的物理特性不连续而导致的分析结果错误，ODDLS技术的基本概念是将共享接口的分析区域分解和分析为包含接口的分析区域。

SUMMARY

    This paper introduces a new stabilized finite element method based on the finite calculus (Comput. Methods Appl. Mech. Eng. 1998; 151:233–267) and arbitrary Lagrangian–Eulerian techniques (Comput. Methods Appl. Mech. Eng. 1998; 155:235–249) for the solution to free surface problems. The main innovation of this method is the application of an overlapping domain decomposition concept in the statement of the problem. The aim is to increase the accuracy in the capture of the free surface as well as in the resolution of the governing equations in the interface between the two fluids. Free surface capturing is based on the solution to a level set equation. The Navier–Stokes equations are solved using an iterative monolithic predictor–corrector algorithm (Encyclopedia of Computational Mechanics. Wiley: New York,2004), where the correction step is based on imposing the divergence-free condition in the velocity field by means of the solution to a scalar equation for the pressure. Examples of application of the ODDLS

formulation (for overlapping domain decomposition level set) to the analysis of different free surface flow problems are presented. Copyright ©2008 John Wiley & Sons, Ltd.

INTRODUCTION

    The prediction of the free surface motion of liquids is a topic of big relevance in many engineering fields. Despite recent advances in computational fluid dynamics, the development of an efficient,accurate and robust numerical algorithm for the analysis of problems with large free surface

deformation is still a challenging issue.

     The free surface flow problem can be considered as a particular case of the more general problem of predicting the interface between two immiscible fluids: the flowing liquid (typically water) and air. Computing the interface between two immiscible fluids is difficult because neither the shape

nor the position of the interface between the two fluids is known a priori. There are basically two approaches for computing free surfaces in this kind of flows: interface-tracking and interfacecapturing methods. The former computes the motion of the flow particles based on a Lagrangian  approach, where the numerical domain adapts itself to the shape and position of the free surface.Different numerical techniques, such as the smoothed particle hydrodynamics (SPH) method [1, 2]and the particle finite element method (PFEM), belong to this kind [3–5]. In interface-tracking methods, the free surface is treated as a boundary of the computational domain where the kinematic and dynamic boundary conditions are applied. The main problems of this approach are the large computational effort required due to the need of updating the analysis domain every time step and the difficulty in imposing mass continuity in an accurate way.Standard interface-capturing methods consider both fluids as a single effective fluid with variable properties [6–11]. The interface is considered as a region of sudden change in the fluid properties.

    This approach requires an accurate modelling of the jump in the properties of the two fluids taking into account that the free surface can move, bend and reconnect in arbitrary ways. Furthermore,the imposition of the exact boundary conditions in the interface is usually simplified.

    The outline of this paper is as follows. In the following section, the statement of the governing equations of two incompressible and immiscible fluids is presented. Then, the treatment of the interfacial boundary condition is analysed and the FIC-stabilized problem is presented. The overlapping domain decomposition methodology is then applied to the stabilized problem. The discretization of the FIC-governing equations using equal-order linear finite element is described and finally, the arbitrary Lagrangian–Eulerian (ALE) version of the method is introduced. The accuracy of the new formulation for the analysis of free surface flows is verified in different validation cases.

STATEMENT OF THE PROBLEM

    The velocity and pressure fields of two incompressible and immiscible fluids moving in the domain during the time interval (0,T] can be described by the incompressible Navier–Stokes equations for multiphase flows, also known as the non-homogeneous incompressible Navier–Stokes equations [20]:

where μ is the dynamic viscosity and I is the identity matrix. Boldface characters are used to denote vector and tensor variables.

   As we consider problems with a moving interface inside the domain, all subdomains and functions defined hereafter are time dependent.

    Forbe the part of the domain Ω occupied by fluid 1 and letbe the part of the domain Ω occupied by fluid 2. Therefore,and are two disjoint subdomains of Ω. Then

where ‘int’ denotes the topological interior and the over-bar indicates the topological adherence of a given set, for more details see [21]. The system of equations (1) must be completed with the necessary initial and boundary conditions, as shown below.

    It is usual in the literature to consider that the first equation of system (1) is equivalent to impose a divergence-free velocity field (the third equation in (1)), as the density is taken as a constant.However, for multiphase incompressible flows, the density cannot be considered to be constant in

Ω×(0,T]. In fact, it is possible to define theρ,μ fields as follows:

Let→be a function named the level set function hereafter and defined as follows:

whereis the distance of point x to the interface between the two fluids, denoted by, at time t .

Therefore, it is possible to re-write definition (4) as

Let us express the density fields in terms of the level set function as

Then, the density derivatives can then be expressed as

Substituting relation (9) into the first equation of system (1) gives

The multiphase Navier–Stokes problem (1) is therefore equivalent to solving the following system of equations:

coupled with the equation

Equation (12) defines the transport of the level set function due to the velocity field obtained by solving Equations (11).

    As a conclusion, the free surface capturing problem can be described by Equations (11) and (12).The interface between the two fluids is defined by the level set 0 of .

   It is possible to prove, assuming the variables of the problem as sufficiently smooth, that Equations (1), or equivalently the system given by Equations (11) and (12), have a unique global solution [21].

   Denoting by an over-bar the prescribed values and by n the normal to the boundary, the boundary conditions of Equations (11) and (12) are

   where the boundary of the domain Ω has been split into three disjoint sets:,, where the Dirichlet and Neumann boundary conditions are imposed, and where mixed conditions are imposed. Mixed Dirichlet–Neumann boundary conditions are usually applied when wall functions

are used for modeling the behaviour of the flow close to solid walls. Vectors g, s in the above equation span the space tangent to . Owing to the fact that (12) is an hyperbolic equation, only boundary conditions must be imposed at the inlet boundary noted by and defined as

Then the boundary condition for Equation (12) is Finally, the initial conditions for the problem are Remark 1

  The initial condition for  can be directly derived from the initial position of the interface. Using the definition of given in Equation (5), can be computed as the signed distance to the initial interface position, . An efficient algorithm to calculate the signed distance to the interface is presented next.

From：