# Improving the Accuracy of Fictitious Domain Method Using Indicator Function from Volume Intersection

Authors: Guoliang Chai; Junwei Su; Le Wang; Chunlei Yu; Yigen Zhang; Zhaolin Gu
Journal: Advances in Mathematical Physics
 (2019)

---

## Abstract

Fictitious domain method (FDM) is a commonly accepted direct numerical simulation technique for moving boundary problems. Indicator function used to distinguish the solid zone and the fluid zone is an essential part concerning the whole prediction accuracy of FDM. In this work, a new indicator function through volume intersection is developed for FDM. In this method, the arbitrarily polyhedral cells across the interface between fluid and solid are located and subdivided into tetrahedrons. The fraction of the solid volume in each cell is accurately computed to achieve high precision of integration calculation in the particle domain, improving the accuracy of the whole method. The quadrature over the solid domain shows that the newly developed indicator function can provide results with high accuracy for variable integration in both stationary and moving boundary problems. Several numerical tests, including flow around a circular cylinder, a single sphere in a creeping shear flow, settlement of a circular particle in a closed container, and in-line oscillation of a circular cylinder, have been performed. The results show good accuracy and feasibility in dealing with the stationary boundary problem as well as the moving boundary problem. This method is accurate and conservative, which can be a feasible tool for studying problems with moving boundaries.


---

## Body

## 1. Introduction

Moving boundary problems appear widely in nature (freely falling leaves) and engineering (particle flow, rotating blades of rotating machinery). In such problems, the boundary movement disturbs the flow of the surrounding fluid, resulting in complicated flow behaviour, which affects the dynamics of the moving object in turn and further changes its trajectory. There is a complicated nonlinear coupling relationship between flow behaviour and its trajectory. It is essential to track the influence of boundary movement on the flow and predict the complex behaviour of the moving boundary system accurately by the feedback of the flow field to the moving object.With the development of the computer, computational fluid dynamics (CFD) has become an important method to study the moving boundary problems. At present, the numerical simulation methods for moving boundary problems can be divided into the point-source method [1–3] and direct numerical simulation method [4–15]. In the point-source method, fluid-solid interactions are evaluated through point source models (drag model, lift model, etc.) and the solid (or particle) is much smaller than the mesh size. This situation is beyond the scope of this work. Direct numerical simulation techniques include arbitrary Lagrange-Euler (ALE) [4–6] and fixed mesh method [7–15]. For ALE, body-fitted mesh is usually applied, and the force and moment of fluid acting on the solid are integrated by the force on the boundary; then, the solid is moved and rotated by tracking the dynamic process and kinematics equation of the solid; finally, the fluid mesh is adjusted to calculate the flow for the next moment. This method can track the flow of fluid precisely. In this method, the fluid mesh needs to be moved when the solid is moving. If the solid domain changes greatly (e.g., the problem of rotating machinery), the mesh will greatly deform, resulting in a great decrease of mesh quality without change of mesh topology, affecting the prediction accuracy. If the mesh topology is changed, more complicated mesh reconstruction will be needed. This method applies only if the solid domain barely changes. For the fixed mesh method, the fluid mesh does not change during the solution procedure, and the behaviour of fluid is changed in the solid domain or solid boundary-covered mesh to indicate the effects of the solid on the flow. The immersed boundary method [7, 8] and fictitious domain method [9–15] belong to such methods. These methods avoid complicated mesh operations, and thus it has become an important method of dealing with moving boundary problems in recent years [16].Fictitious domain method is an important method for solving the moving boundary problems. This method divides the solution domain into fluid covering domain and solid covering domain. In the solid covering domain, a pseudo body force is exerted to make the fluid behaviour become consistent with the behaviour of the solid. This method was first proposed by Glowinski et al. [9]. They introduced the body force as a Lagrange multiplier to the momentum equation and solved it by iteration. Because of the saddle point problem, the computational load is increased. In order to overcome this disadvantage, Sharma and Patankar [11] proposed a non-Lagrangian multiplier/fictitious domain method, which simplifies the computation by calculating the body force explicitly, thereby significantly reducing the computational load. However, this method defines the body force at some Lagrangian points on the solid domain, while the particle moves in the Eulerian background mesh. Thus, the scheme of the body force mapping between the Eulerian mesh and the Lagrangian point determines the numerical accuracy. The nonsmoothness of the field obtained by the mapping incurs numerical instability. To overcome this problem, Gallier et al. [14] proposed a fictitious domain method in a fully Eulerian framework. The body force is defined in the centre of the mesh cell, rather than Lagrangian points. All the solution procedures are performed on the Eulerian mesh, and the information exchange between the Eulerian framework and the Lagrangian frame is avoided. However, the above process is achieved by defining a particle indicator function, which is 0 outside the particle and 1 inside the particle and smoothly transits from 0 to 1 within a grid near the particle boundary. The smooth transition process is specified by a Heaviside function, rather than the volume ratio of the particle domain to the arbitrarily polyhedral mesh cells, so the integral calculated by this particle indicator function does not possess the feature of conservation.Based on the fully Eulerian fictitious domain method of Gallier et al., a fictitious domain method based on the particle indicator function using intersection volume is proposed in this paper. This method marks the Eulerian arbitrarily polyhedral mesh cells at the interface between fluid and particles, firstly. The mesh cells are divided into a serial of tetrahedrons by a cell decomposing technique, and thus, the intersection volume between the particle and the mesh cell can be evaluated accurately. Adopting the interaction volume as the indicator function can guarantee the accuracy and conservation of the integral of the particle domain, thereby enhancing the overall numerical accuracy.


## 2. Fictitious Domain Method

In this paper, Gallier et al.’s fictitious domain method [14] is employed as an example to test the volume-intersection based indicator function. In order to make the method adapt to arbitrarily polyhedral mesh, the collocated finite volume method is employed in this work rather than the finite-difference method adopted in original work.


### 2.1. Governing Equations

The computational domainD contains fluid domain Df and particle domain DP, i.e., D = Df∪DP. The particle domain contains N rigid particles Pi with the density of ρP, and the interface between particle and fluid is ∂P. The fluid is incompressible viscous Newtonian fluid with the density of ρf, the dynamic viscous coefficient of μ, and the kinematic viscous coefficient of ν=μ/ρf. The governing equation of the fictitious domain model is expressed as follows(1)∂u∂t+u∇·u=-∇pρf+ν∇2u+λinD(2)∇·u=0inD(3)u=U+Ω×rinDP(4)ρf∫Pλdx=ρf-ρPρPMdUdt-g+Fint(5)ρf∫Pr×λdx=ρf-ρPρPJ·dΩdt-Ω×J·Ω+Tintwhere u and p are the velocity and pressures of fluid, respectively; λ is a pseudo body force, which is nonzero in the particle domain and zero elsewhere, and therefore the Navier-Stokes equation is also valid in the particle domain; M, J, U, and Ω are the mass, moment of inertia, translational velocity, and rotational angular velocity of a single particle, respectively; r is the position vector at the mass centre of a particle; g is gravitational acceleration; Fint and Tint are the forces and moments acting on the particle except for the fluid force.In the above governing equations, (1) and (2) are the momentum equation and continuity equation satisfied in the whole solution domain, respectively. The body force, λ, is enforced in the solid domain to make momentum equation valid in the whole solution domain. Equation (3) enforces rigid body constraint in the particle domain. Equations (4) and (5) are particle dynamic equations.


### 2.2. Numerical Procedure

Two subproblems can be divided based on the above equation descriptions, i.e., the fluid subproblems of (1) and (2) and the particle subproblems of (3)-(5).


#### 2.2.1. Fluid Subproblem

Equation (1) is the equation describing the fluid dynamics, which Gallier et al. [14] used the finite-difference method to discretize. In this paper, the finite volume method is employed. The semidiscrete form of (1) can be expressed as(6)apu=AH-∇pρfwhere ap is the diagonal value of the coefficient matrix obtained by discretization of (1), and AH is expressed as(7)AH=∑NaNu+bwhere aN is the implicit contribution coefficient of the neighbour cells to the current cell during discretization of (1) and b is the explicit discretization contribution except pressure. Since λ is unknown at this time, the explicit discretization is used here.Both sides of (6) are divided by ap yielding(8)u=AHap-∇pρfapThe velocity field obtained by (8) should satisfy the continuity equation. Substituting (8) into (2) yields(9)∇·1ρfap∇p=∇·AHapEquation (9) is a pressure equation derived by PISO algorithm [19]. A new pressure can be obtained by solving the equation. Equation (8) is used to update the velocity. In this paper, the collocated finite volume method is used for discretization. While the velocity on the mesh cell centre is updated, it is necessary to update the flow rate at the cell face using(10)φf=AHapf·Sf-∇pρfapf·Sfwhere mf is the interpolation of value of the body centre, m, on the face centre, Sf is the area vector of the cell face, and φf is the volume flow rate on the cell face.


#### 2.2.2. Particle Subproblem

For the particle subproblem, the following equations are solved first(11)MUn+1=ρP-ρfρPMUn+gΔt+ρf∫Pun+1-λnΔtdx+Fint,nΔt(12)J·Ωn+1=ρP-ρfρPJ·Ωn-Ωn×J·ΩnΔt+ρf∫Pr×un+1-λnΔtdx+Tint,nΔtto obtain the translational velocity, Un+1, and rotational angular velocity, Ωn+1, at the new time step. Then, use the equation below(13)u~n+1=Un+1+Ωn+1×rinDPto update the fluid velocity in the particle covering domain and use this velocity to update the pseudo body force, λn+1.(14)λn+1=λn+u~n+1-un+1ΔtFinally, the second-order Adams-Bashforth difference scheme is used to discretize the displacement equation and update the particle position.(15)xn+1=xn+32ΔtUn+1-12ΔtUnFor nonspherical particle, it is necessary to rotate the particle by the angular velocity.


#### 2.2.3. Calculation of Integral of Particle Domain

Equations (11) and (12) contain the integrals of physical quantities on the particle domain. As a DNS (direct numerical simulation) method, the faces of background mesh cells do not coincide with the particle surface during the solving procedure, so it is necessary to use an approximation method to deal with the interface area. Gallier et al. [14] used the traditional approximate Heaviside function to calculate the integral term. The integral of any physical quantity, Φ, in the particle covering domain is approximated by (16)∫PΦdx=∫DPIxΦdxwhere I(x) is the indicator function, which is expressed by the following expression(17)Ix=1-Hψxwhere ψ(x) is the level-set function and is expressed as(18)ψx=x-xc-rwhere xc and r are the mass centre and radius of the particle, respectively; ψ(x) is a negative value inside the particle and a positive value outside the particle. The approximate Heaviside function can be expressed as(19)Hs=121+tanh⁡bsΔwhere ∆ is the grid spacing and the constant b=5 guarantees the smooth width within a grid. Eventually, the particle indicator function I(x) is 1 inside the particle and 0 outside the particle and smoothly transits from 1 to 0 at the interface of the two phases.The approximate Heaviside function guarantees a smooth transition within a grid around the interface and enhances the stability of the algorithm. However, it cannot guarantee the conservation. If the integrand is a constant, 1, the volume of the particle should be obtained by integration of (16). However, the use of approximate Heaviside function leads to a big difference between the calculated volume and actual volume of the particle, i.e., nonconservation of the volume. In the particle subproblem, (11)-(12) involve the integration of velocity and pseudo body force in the particle domain. The integration with approximate Heaviside function sacrifices the numerical accuracy.


## 2.1. Governing Equations

The computational domainD contains fluid domain Df and particle domain DP, i.e., D = Df∪DP. The particle domain contains N rigid particles Pi with the density of ρP, and the interface between particle and fluid is ∂P. The fluid is incompressible viscous Newtonian fluid with the density of ρf, the dynamic viscous coefficient of μ, and the kinematic viscous coefficient of ν=μ/ρf. The governing equation of the fictitious domain model is expressed as follows(1)∂u∂t+u∇·u=-∇pρf+ν∇2u+λinD(2)∇·u=0inD(3)u=U+Ω×rinDP(4)ρf∫Pλdx=ρf-ρPρPMdUdt-g+Fint(5)ρf∫Pr×λdx=ρf-ρPρPJ·dΩdt-Ω×J·Ω+Tintwhere u and p are the velocity and pressures of fluid, respectively; λ is a pseudo body force, which is nonzero in the particle domain and zero elsewhere, and therefore the Navier-Stokes equation is also valid in the particle domain; M, J, U, and Ω are the mass, moment of inertia, translational velocity, and rotational angular velocity of a single particle, respectively; r is the position vector at the mass centre of a particle; g is gravitational acceleration; Fint and Tint are the forces and moments acting on the particle except for the fluid force.In the above governing equations, (1) and (2) are the momentum equation and continuity equation satisfied in the whole solution domain, respectively. The body force, λ, is enforced in the solid domain to make momentum equation valid in the whole solution domain. Equation (3) enforces rigid body constraint in the particle domain. Equations (4) and (5) are particle dynamic equations.


## 2.2. Numerical Procedure

Two subproblems can be divided based on the above equation descriptions, i.e., the fluid subproblems of (1) and (2) and the particle subproblems of (3)-(5).


### 2.2.1. Fluid Subproblem

Equation (1) is the equation describing the fluid dynamics, which Gallier et al. [14] used the finite-difference method to discretize. In this paper, the finite volume method is employed. The semidiscrete form of (1) can be expressed as(6)apu=AH-∇pρfwhere ap is the diagonal value of the coefficient matrix obtained by discretization of (1), and AH is expressed as(7)AH=∑NaNu+bwhere aN is the implicit contribution coefficient of the neighbour cells to the current cell during discretization of (1) and b is the explicit discretization contribution except pressure. Since λ is unknown at this time, the explicit discretization is used here.Both sides of (6) are divided by ap yielding(8)u=AHap-∇pρfapThe velocity field obtained by (8) should satisfy the continuity equation. Substituting (8) into (2) yields(9)∇·1ρfap∇p=∇·AHapEquation (9) is a pressure equation derived by PISO algorithm [19]. A new pressure can be obtained by solving the equation. Equation (8) is used to update the velocity. In this paper, the collocated finite volume method is used for discretization. While the velocity on the mesh cell centre is updated, it is necessary to update the flow rate at the cell face using(10)φf=AHapf·Sf-∇pρfapf·Sfwhere mf is the interpolation of value of the body centre, m, on the face centre, Sf is the area vector of the cell face, and φf is the volume flow rate on the cell face.


### 2.2.2. Particle Subproblem

For the particle subproblem, the following equations are solved first(11)MUn+1=ρP-ρfρPMUn+gΔt+ρf∫Pun+1-λnΔtdx+Fint,nΔt(12)J·Ωn+1=ρP-ρfρPJ·Ωn-Ωn×J·ΩnΔt+ρf∫Pr×un+1-λnΔtdx+Tint,nΔtto obtain the translational velocity, Un+1, and rotational angular velocity, Ωn+1, at the new time step. Then, use the equation below(13)u~n+1=Un+1+Ωn+1×rinDPto update the fluid velocity in the particle covering domain and use this velocity to update the pseudo body force, λn+1.(14)λn+1=λn+u~n+1-un+1ΔtFinally, the second-order Adams-Bashforth difference scheme is used to discretize the displacement equation and update the particle position.(15)xn+1=xn+32ΔtUn+1-12ΔtUnFor nonspherical particle, it is necessary to rotate the particle by the angular velocity.


### 2.2.3. Calculation of Integral of Particle Domain

Equations (11) and (12) contain the integrals of physical quantities on the particle domain. As a DNS (direct numerical simulation) method, the faces of background mesh cells do not coincide with the particle surface during the solving procedure, so it is necessary to use an approximation method to deal with the interface area. Gallier et al. [14] used the traditional approximate Heaviside function to calculate the integral term. The integral of any physical quantity, Φ, in the particle covering domain is approximated by (16)∫PΦdx=∫DPIxΦdxwhere I(x) is the indicator function, which is expressed by the following expression(17)Ix=1-Hψxwhere ψ(x) is the level-set function and is expressed as(18)ψx=x-xc-rwhere xc and r are the mass centre and radius of the particle, respectively; ψ(x) is a negative value inside the particle and a positive value outside the particle. The approximate Heaviside function can be expressed as(19)Hs=121+tanh⁡bsΔwhere ∆ is the grid spacing and the constant b=5 guarantees the smooth width within a grid. Eventually, the particle indicator function I(x) is 1 inside the particle and 0 outside the particle and smoothly transits from 1 to 0 at the interface of the two phases.The approximate Heaviside function guarantees a smooth transition within a grid around the interface and enhances the stability of the algorithm. However, it cannot guarantee the conservation. If the integrand is a constant, 1, the volume of the particle should be obtained by integration of (16). However, the use of approximate Heaviside function leads to a big difference between the calculated volume and actual volume of the particle, i.e., nonconservation of the volume. In the particle subproblem, (11)-(12) involve the integration of velocity and pseudo body force in the particle domain. The integration with approximate Heaviside function sacrifices the numerical accuracy.


## 2.2.1. Fluid Subproblem

Equation (1) is the equation describing the fluid dynamics, which Gallier et al. [14] used the finite-difference method to discretize. In this paper, the finite volume method is employed. The semidiscrete form of (1) can be expressed as(6)apu=AH-∇pρfwhere ap is the diagonal value of the coefficient matrix obtained by discretization of (1), and AH is expressed as(7)AH=∑NaNu+bwhere aN is the implicit contribution coefficient of the neighbour cells to the current cell during discretization of (1) and b is the explicit discretization contribution except pressure. Since λ is unknown at this time, the explicit discretization is used here.Both sides of (6) are divided by ap yielding(8)u=AHap-∇pρfapThe velocity field obtained by (8) should satisfy the continuity equation. Substituting (8) into (2) yields(9)∇·1ρfap∇p=∇·AHapEquation (9) is a pressure equation derived by PISO algorithm [19]. A new pressure can be obtained by solving the equation. Equation (8) is used to update the velocity. In this paper, the collocated finite volume method is used for discretization. While the velocity on the mesh cell centre is updated, it is necessary to update the flow rate at the cell face using(10)φf=AHapf·Sf-∇pρfapf·Sfwhere mf is the interpolation of value of the body centre, m, on the face centre, Sf is the area vector of the cell face, and φf is the volume flow rate on the cell face.


## 2.2.2. Particle Subproblem

For the particle subproblem, the following equations are solved first(11)MUn+1=ρP-ρfρPMUn+gΔt+ρf∫Pun+1-λnΔtdx+Fint,nΔt(12)J·Ωn+1=ρP-ρfρPJ·Ωn-Ωn×J·ΩnΔt+ρf∫Pr×un+1-λnΔtdx+Tint,nΔtto obtain the translational velocity, Un+1, and rotational angular velocity, Ωn+1, at the new time step. Then, use the equation below(13)u~n+1=Un+1+Ωn+1×rinDPto update the fluid velocity in the particle covering domain and use this velocity to update the pseudo body force, λn+1.(14)λn+1=λn+u~n+1-un+1ΔtFinally, the second-order Adams-Bashforth difference scheme is used to discretize the displacement equation and update the particle position.(15)xn+1=xn+32ΔtUn+1-12ΔtUnFor nonspherical particle, it is necessary to rotate the particle by the angular velocity.


## 2.2.3. Calculation of Integral of Particle Domain

Equations (11) and (12) contain the integrals of physical quantities on the particle domain. As a DNS (direct numerical simulation) method, the faces of background mesh cells do not coincide with the particle surface during the solving procedure, so it is necessary to use an approximation method to deal with the interface area. Gallier et al. [14] used the traditional approximate Heaviside function to calculate the integral term. The integral of any physical quantity, Φ, in the particle covering domain is approximated by (16)∫PΦdx=∫DPIxΦdxwhere I(x) is the indicator function, which is expressed by the following expression(17)Ix=1-Hψxwhere ψ(x) is the level-set function and is expressed as(18)ψx=x-xc-rwhere xc and r are the mass centre and radius of the particle, respectively; ψ(x) is a negative value inside the particle and a positive value outside the particle. The approximate Heaviside function can be expressed as(19)Hs=121+tanh⁡bsΔwhere ∆ is the grid spacing and the constant b=5 guarantees the smooth width within a grid. Eventually, the particle indicator function I(x) is 1 inside the particle and 0 outside the particle and smoothly transits from 1 to 0 at the interface of the two phases.The approximate Heaviside function guarantees a smooth transition within a grid around the interface and enhances the stability of the algorithm. However, it cannot guarantee the conservation. If the integrand is a constant, 1, the volume of the particle should be obtained by integration of (16). However, the use of approximate Heaviside function leads to a big difference between the calculated volume and actual volume of the particle, i.e., nonconservation of the volume. In the particle subproblem, (11)-(12) involve the integration of velocity and pseudo body force in the particle domain. The integration with approximate Heaviside function sacrifices the numerical accuracy.


## 3. Intersecting Volume-Based Particle Indicator Function

### 3.1. Locating the Cells near the Solid Boundary

Based on the level-set function from (18), the mesh cells of ψ(x)> 0 are set to I(x) = 0 first, and the rest are set to I(x) = 1. The cells near the solid boundary are located using the level set functions ψ(xowner) and ψ(xneighbour). The subscripts owner and neighbour represent the current cell and the neighbour cell, respectively. If ψ(xowner)ψ(xneighbour) < 0, two cells are located on both sides of the interface, respectively. Thus, they are marked as the cells near the solid interface. All the cells are checked to locate the mesh cells near the interfaces. This is illustrated in Figure 1, where the circular solid line is the actual boundary of the circled particle, and the shaded cells are the marked ones at the particle interface.Figure 1
Mesh cells near the particle boundary.


### 3.2. Subdivision of Mesh Cells near the Solid Boundary

In order to calculate the intersection volumes between the mesh cells and the particles accurately, the cells near the interface are decomposed into tetrahedrons with different levels. The arbitrarily polyhedral cell is taken as an example. The process is as follows.Firstly, an arbitrary polyhedron is decomposed into several tetrahedrons.Figure2(a) is a polyhedral cell centred at Xi. The method given by Su et al. [20] is applied to decompose the polyhedron. In this method, we connect the centre to each vertex of the polyhedron and twelve pentagonal pyramids are obtained. One of the pyramids, ν0ν1ν2ν3ν4Xi, is taken separately, as shown in Figure 2(b), and Xf and Sf are the face centre and area vector of the pentagon ν0ν1ν2ν3ν4, respectively. Three tetrahedrons, ν0ν1ν2Xi, ν0ν2ν3Xi, and ν0ν3ν4Xi, can be obtained by connecting ν0, ν2, and ν3. The above operation is repeated for the rest of the pyramids so that the polyhedron can be decomposed into 36 tetrahedrons as shown in Figure 2(c).Figure 2
Decomposition of a polyhedral cell into tetrahedrons: (a) original polyhedron cut into pyramids; (b) a pyramid cut into tetrahedrons; (c) decomposition result.
(a) (b) (c)Secondly, a tetrahedron can be decomposed into eight small tetrahedrons further to improve the accuracy.The tetrahedron shown in Figure3(a) can be subdivided into eight subtetrahedrons. Four small tetrahedrons and one octahedron can be obtained by connecting any two of the midpoints of six edges of the tetrahedron along four surfaces as shown in Figure 3(b). The central octahedron given in Figure 3(c) can be subdivided into four tetrahedrons by connecting a pair of nonadjacent vertices, and the final tetrahedrons are shown in Figure 3(d). Finally, eight small tetrahedrons are obtained. According to this strategy, the tetrahedrons obtained can be decomposed into several smaller tetrahedrons further. The process of decomposing a polyhedron into several tetrahedrons is called Level 0 subdivision. The decomposed tetrahedron can be decomposed further if needed. Each subsequent decomposition process is called one higher level subdivision. The number of tetrahedrons after decomposition is eight times the number before the decomposition.Figure 3
Subdivision of a tetrahedron: (a) original tetrahedron; (b) initial subdivision; (c) central octahedron; (d) subdivision of central octahedron [17].
(a) (b) (c) (d)


### 3.3. Calculation of Particle Indicator FunctionI(x) of the Cells near the Solid Boundary

After subdivision, the mesh cells near the solid boundary are divided into several tetrahedrons, which are in three states, outside the particle, inside the particle, and across the solid interface. For the last one, an algorithm similar to Rabbitz’s algorithm [17] is used to determine the intersection points between the interface and the tetrahedron edges. Rabbitz determines the intersection points by planes intersection. In this paper, the level-set function obtained by (18) is used to determine the intersection points. For the level-set value of any two vertices, v1 and v2, of a tetrahedron across the interface, if ψ(xv1)ψ(xv2)<0, the edge containing v1 and v2 intersects with the interface, and the intersection point xi can be determined by computing ψ(xi) = 0 on that edge.Then a tetrahedron is divided into two parts by connecting any two of these intersection points along the tetrahedron surfaces (Figure4(a), the shaded plane), and one shaded part is inside the particle (Figures 4(b) and 4(c)), which is either a tetrahedron or not. If not, it can be further decomposed into several tetrahedrons. And with all these tetrahedrons, the volume of the solid part is easy to obtain.Figure 4
All the circumstances of (a) a tetrahedron intersected by a surface and resulting in (b) one or (c) another shaded geometry within the surface, which is already a tetrahedron or can be decomposed into several tetrahedrons by the red lines.
(a) (b) (c)The solid volume in a cell near the solid boundary is the sum of the volumes,∑Vp, of all geometries inside the particle (including the tetrahedrons inside and the shaded part of the tetrahedron across the interface that is inside). The solid volume fraction evaluated by ∑Vp / Vcell is adopted as the particle indicator function I(x) of the cell.Notice that we assume the interface within the tetrahedron is flat (Figure4(a)), which is not true. Numerical errors are inevitable. However, theoretically, a higher subdivision level leads to a smaller error. Thus, with sufficient levels of the subdivision, the required numerical accuracy can be obtained.


## 3.1. Locating the Cells near the Solid Boundary

Based on the level-set function from (18), the mesh cells of ψ(x)> 0 are set to I(x) = 0 first, and the rest are set to I(x) = 1. The cells near the solid boundary are located using the level set functions ψ(xowner) and ψ(xneighbour). The subscripts owner and neighbour represent the current cell and the neighbour cell, respectively. If ψ(xowner)ψ(xneighbour) < 0, two cells are located on both sides of the interface, respectively. Thus, they are marked as the cells near the solid interface. All the cells are checked to locate the mesh cells near the interfaces. This is illustrated in Figure 1, where the circular solid line is the actual boundary of the circled particle, and the shaded cells are the marked ones at the particle interface.Figure 1
Mesh cells near the particle boundary.


## 3.2. Subdivision of Mesh Cells near the Solid Boundary

In order to calculate the intersection volumes between the mesh cells and the particles accurately, the cells near the interface are decomposed into tetrahedrons with different levels. The arbitrarily polyhedral cell is taken as an example. The process is as follows.Firstly, an arbitrary polyhedron is decomposed into several tetrahedrons.Figure2(a) is a polyhedral cell centred at Xi. The method given by Su et al. [20] is applied to decompose the polyhedron. In this method, we connect the centre to each vertex of the polyhedron and twelve pentagonal pyramids are obtained. One of the pyramids, ν0ν1ν2ν3ν4Xi, is taken separately, as shown in Figure 2(b), and Xf and Sf are the face centre and area vector of the pentagon ν0ν1ν2ν3ν4, respectively. Three tetrahedrons, ν0ν1ν2Xi, ν0ν2ν3Xi, and ν0ν3ν4Xi, can be obtained by connecting ν0, ν2, and ν3. The above operation is repeated for the rest of the pyramids so that the polyhedron can be decomposed into 36 tetrahedrons as shown in Figure 2(c).Figure 2
Decomposition of a polyhedral cell into tetrahedrons: (a) original polyhedron cut into pyramids; (b) a pyramid cut into tetrahedrons; (c) decomposition result.
(a) (b) (c)Secondly, a tetrahedron can be decomposed into eight small tetrahedrons further to improve the accuracy.The tetrahedron shown in Figure3(a) can be subdivided into eight subtetrahedrons. Four small tetrahedrons and one octahedron can be obtained by connecting any two of the midpoints of six edges of the tetrahedron along four surfaces as shown in Figure 3(b). The central octahedron given in Figure 3(c) can be subdivided into four tetrahedrons by connecting a pair of nonadjacent vertices, and the final tetrahedrons are shown in Figure 3(d). Finally, eight small tetrahedrons are obtained. According to this strategy, the tetrahedrons obtained can be decomposed into several smaller tetrahedrons further. The process of decomposing a polyhedron into several tetrahedrons is called Level 0 subdivision. The decomposed tetrahedron can be decomposed further if needed. Each subsequent decomposition process is called one higher level subdivision. The number of tetrahedrons after decomposition is eight times the number before the decomposition.Figure 3
Subdivision of a tetrahedron: (a) original tetrahedron; (b) initial subdivision; (c) central octahedron; (d) subdivision of central octahedron [17].
(a) (b) (c) (d)


## 3.3. Calculation of Particle Indicator FunctionI(x) of the Cells near the Solid Boundary

After subdivision, the mesh cells near the solid boundary are divided into several tetrahedrons, which are in three states, outside the particle, inside the particle, and across the solid interface. For the last one, an algorithm similar to Rabbitz’s algorithm [17] is used to determine the intersection points between the interface and the tetrahedron edges. Rabbitz determines the intersection points by planes intersection. In this paper, the level-set function obtained by (18) is used to determine the intersection points. For the level-set value of any two vertices, v1 and v2, of a tetrahedron across the interface, if ψ(xv1)ψ(xv2)<0, the edge containing v1 and v2 intersects with the interface, and the intersection point xi can be determined by computing ψ(xi) = 0 on that edge.Then a tetrahedron is divided into two parts by connecting any two of these intersection points along the tetrahedron surfaces (Figure4(a), the shaded plane), and one shaded part is inside the particle (Figures 4(b) and 4(c)), which is either a tetrahedron or not. If not, it can be further decomposed into several tetrahedrons. And with all these tetrahedrons, the volume of the solid part is easy to obtain.Figure 4
All the circumstances of (a) a tetrahedron intersected by a surface and resulting in (b) one or (c) another shaded geometry within the surface, which is already a tetrahedron or can be decomposed into several tetrahedrons by the red lines.
(a) (b) (c)The solid volume in a cell near the solid boundary is the sum of the volumes,∑Vp, of all geometries inside the particle (including the tetrahedrons inside and the shaded part of the tetrahedron across the interface that is inside). The solid volume fraction evaluated by ∑Vp / Vcell is adopted as the particle indicator function I(x) of the cell.Notice that we assume the interface within the tetrahedron is flat (Figure4(a)), which is not true. Numerical errors are inevitable. However, theoretically, a higher subdivision level leads to a smaller error. Thus, with sufficient levels of the subdivision, the required numerical accuracy can be obtained.


## 4. Results and Discussion

### 4.1. Particle Volume Evaluation

In this example, the accuracy of the volume calculation of circular and spherical particles on 2D and 3D structured background mesh at different subdivision levels. The circular and spherical particles have the diameter,d = 1 m, the grid spacing is set to d/10, d/20, d/30, and d/40, respectively, and the subdivision level is set to 1, 2, 3, and 4, totally 16 cases. Tables 1 and 2 give the particle volumes ∑IxVcell and relative errors by two indictor functions in 2D and 3D, respectively. The exact volumes for 2D circle and 3D volume are 0.785398 and 0.523599, respectively. It can be observed from the tables that, with a grid fining, our method produces much smaller relative errors than the approximate Heaviside function method, showing interaction-volume based particle indicator function I(x) is more accurate.Table 1
Comparison of results and errors of the particle volume computation in 2D.
Grid Spacing Volume and Errors Gallier [14] Subdivision level 1 2 3 4 d/10 Volume 0.793864 0.785212 0.785347 0.785384 0.785394 Error 1.08e-2 2.37e-4 6.50e-5 1.80e-5 5.02e-6 d/20 Volume 0.788934 0.785353 0.785386 0.785395 0.785397 Error 4.50e-3 5.92e-5 1.56e-5 4.33e-6 1.21e-6 d/30 Volume 0.787999 0.785379 0.785393 0.785397 0.785398 Error 3.31e-3 2.50e-5 6.90e-6 1.91e-6 5.31e-7 d/40 Volume 0.787515 0.785387 0.785395 0.785397 0.785398 Error 2.70e-3 1.41e-5 3.88e-6 1.07e-6 2.96e-7Table 2
Comparison of results and errors of the particle volume computation in 3D.
Grid spacing Volume and Errors Gallier [14] Subdivision level 1 2 3 4 d/10 Volume 0.529842 0.523251 0.523503 0.523572 0.523591 Error 1.19e-2 6.64e-4 1.83e-4 5.09e-5 1.42e-5 d/20 Volume 0.527607 0.523512 0.523575 0.523592 0.523597 Error 7.66e-3 1.66e-4 4.60e-5 1.28e-5 3.56e-6 d/30 Volume 0.525825 0.52356 0.523588 0.523596 0.523598 Error 4.25e-3 7.39e-5 2.04e-5 5.66e-6 1.58e-6 d/40 Volume 0.524104 0.523577 0.523593 0.523597 0.523598 Error 9.65e-4 4.15e-5 1.15e-5 3.19e-6 8.94e-7It can be found from Tables1 and 2 that the relative error is less than 5e-5 if the subdivision level is greater than or equal to 2 for the grid spacing with d/20. The subdivision level is set to 2 in the later tests, which is a good compromise between computational load and accuracy.Figure5 shows the evolution of particle volume with time obtained by two methods when the particle moves on the background mesh in 2D and 3D, respectively. It is not difficult to find that approximated Heaviside function method is not conservative because the results significantly oscillate with the particle motion, while our method in this paper is conservative and the results are stable.Figure 5
Evolution of predicted volume of a moving particle with time: (a) 2D; (b) 3D. solid line: the present numerical results; dash line: numerical results with Gallier’s method [14].
(a) (b)


### 4.2. Flow around a Cylinder

The flow around a cylinder is a classical hydrodynamic problem. A large number of physical and numerical experiments have been carried out previously. The accuracy and validity of this method in the stationary boundary are verified by the flow around a circular cylinder with different Reynolds numbers. In this example, the computational domain is a40d×20d rectangular domain (diameter of the cylinder: d = 0.025 m). The centre of the cylinder is located at (10d, 10d), as shown in Figure 6. The fluid flows into the left inlet and out of the right outlet, with the inflow velocity, u (different for different Re), fluid density, ρf = 1.293 kg·m−3, and kinetic viscosity coefficient, ν = 1.5×10−5 m2·s−1. The inlet flow velocity is specified, and the pressure is zero-gradient, while at the outlet, Neumann boundary condition is specified for velocity and Dirichlet boundary condition for pressure. The walls at top and bottom are set as the symmetric boundary condition. The Reynolds numbers,Re(= ρfd|u| / ν), of 10, 20, 40, 100, and 200, are selected as test conditions in this example, with the corresponding inlet velocities of 6×10−3 m·s−1, 1.2×10−2 m·s−1, 2.4×10−2 m·s−1, 6×10−2 m·s−1, and 0.12 m·s−1, respectively. In order to reduce the computational cost, a refinement mesh is used within a 5d×5d square, whose centre is the same as the cylinder’s and edges are parallel to the computational domain’s, and d/5×d/5 meshes are used at the rest of the area in this example. Three levels of grid spacings are employed: d/10 for Case 1, d/20 for Case 2, and d/30 for Case 3, respectively. The corresponding time steps are set so that the Courant number Co (= U0d/Δt) maintains around 0.2.Figure 6
Sketch of the computational domain for flow around a circular cylinder.Table3 gives the drag coefficients Cd, at different Re and grid spacings, at the ultimately stable or cyclically stable flow around a circular cylinder. The drag coefficients are obtained from the following equation.(20)Cd=2Fhydρfdu2The computed Cd are convergent with the decrease of grid spacing. For the purpose of comparison, this table also gives the numerical results from the method in [14] (Gallier) and the experimental and numerical results in [21–23]. It can be seen from the table that our results are better almost in all cases than Gallier’s method and closer to the experimental and calculated results in the references.Table 3
Comparison of the present computation of drag coefficient (Cd) with other results for flow past a cylinder.
R e C d Present Gallier et al. [14] Nishioka et al. [21] Rajani et al. [22] Takami et al. [23] Case1 Case2 Case3 (Numerical) (Experimental) (Numerical) (Numerical) 10 2.991 3.037 3.049 3.116 2.74 2.839 2.754 20 2.126 2.164 2.174 2.236 1.95 2.059 2.003 40 1.573 1.604 1.611 1.668 1.49 1.539 1.536 100 1.302 1.321 1.328 1.376 — 1.335 — 200 1.246 1.335 1.381 1.402 — 1.337 —Figure7 shows the streamlines at Re of 10, 20, and 40 of Case 3 when the flow becomes eventually stable. It can be seen that stable recirculation zones (whirlpool vortex) are generated behind the cylinders in all three cases, and the length, Lw, of the recirculation zone increases with the increasing Re. Table 4 shows the comparison of the length, Lw/d (d is the diameter of the cylinder), of recirculation zone calculated by our method at Re of 20 and 40 for different grid spacings with those in [23–25]. The calculated results are convergent with the refinement of the mesh and in good agreement with those in the references.Table 4
The non-dimensional length of recirculation zone (Lw/d) at Re = 20, 40.
R e L w / d Present Takami et al. [23] Silva et al. [24] Ye et al. [25] Case1 Case2 Case3 (Numerical) (Numerical) (Numerical) 20 0.95 1.05 1.07 0.932 1.05 0.92 40 2.03 2.25 2.31 2.325 2.54 2.27Figure 7
Streamlines for flow around a circular cylinder of Case 3 atRe = (a) 10; (b) 20; (c) 40.
(a) (b) (c)Figure8 shows the velocity contour of cyclically stable flow at Re of 100 and 200 of Case 3. In both cases, the whirlpool vortex behind the cylinder begins to regularly shed and forms a regular vortex street. Table 5 shows the comparison of Strauhal number (St=fd / u, where f is the vortex shedding frequency) calculated by this method in both cases with those in [7, 26, 27]. The calculated results are converted with the decrease of grid spacing and are in good agreement with those in the references.Table 5
Comparison of the numerical predicted Strouhal number (St) with different methods for flow around a circular cylinder.
R e S t Present Liu et al. [26] Uhlmann [7] Tseng and Ferziger [27] Case1 Case2 Case3 100 0.162 0.165 0.166 0.164 0.169 0.164 200 0.189 0.193 0.194 0.192 — —Figure 8
Velocity contour for flow around a circular cylinder of Case 3 atRe = (a) 100, (b) 200.
(a) (b)The simulation results in this example indicate that, for the static boundary problem, the method in this paper is of good accuracy in a wide range ofRe.


### 4.3. Single Spherical Particle in a Creeping Shear Flow

In this section, a single force-free torque-free sphere in a linear creeping shear flow is considered to verify the feasibility of this method in dealing with sphere rotation. A spherical particle of diameterd=0.1 m is suspended at the centre of a cubic Couette cell, as shown in Figure 9. The particle is located at x = 0 and the cell size is L=10d. A velocity of γL/2 and -γL/2 in the flow direction is assigned to the top and bottom boundaries (the dark ones), respectively, where γ=6.71×10−3s−1 is the shear rate. And for the rest of the boundaries, periodic conditions are used. The kinematic viscosity of the flow is ν = 10−3 m2s−1. Three levels of grid spacings are used: Δ=d/5 for Case 1, Δ=d/7.5 for Case 2, and Δ=d/10 for Case 3. And the corresponding time steps are 0.01s, 0.005s, and 0.004s, respectively.Figure 9
Sketch of the computational domain for a single spherical particle in creeping shear flow.In a creeping shear flow, the effect of convection term is so little that it can be neglected, and unsteady Stokes equations are solved with respect to time-stepping. Eventually, the whole system reaches a steady state, and the particle is rotating in a constant angular velocity, which is available of analytical solutions [14].Figure10 shows a comparison between the profile of the computed nondimensional velocity 2u2/γd on different grid spacings on the centreline (x2=0) of the shear plane (x3=0), where the subscript represents the direction of the flow, and the analytical solution [28]. For the coarsest mesh, Case 1, the computed results near the outside region of the particle surface are slightly lower than the analytical solution. With the refinement of the mesh, the results become more and more close to the analytical solution both inside and outside the particle. The linear profile inside the particle (2|x1|/d<1) is due to the rigid body motion constraint of the particle rotation and a velocity 2|u2|/γd=0.5 is found at the surface of the particle (2|x1|/d=1). The rotational velocities Ω3/γ of all three cases are -0.4990, -0.4995, and -0.4996, respectively, and gradually approach the theoretical value -0.5, The particle translational velocity is also solved and turns out to be of the order of 10−8, which is acceptable with respect to the expected value 0.Figure 10
Evolution of the centreline fluid velocity2u2/γd with position 2x1/d for a single sphere located at the origin (x=0).The simulation results in this case show that this method is feasible for the simulation of spherical particle rotation.


### 4.4. Sedimentation of a Circular Particle in a Closed Container

The sedimentation process of a circular particle in a closed container is simulated in this case to verify the feasibility of this method in dealing with moving boundary problems. The computational model is shown in Figure11. The computational domain is 8d×24d (diameter of a particle: d = 2.5 × 10−3 m), and the initial position of the particle is (4d,8d). The fluid is stationary at the initial time. The circular particle settles by gravity, whose density is ρp = 1250 kg·m−3. The dynamic viscosity and density of the fluid are μ = 0.01 kg·m−1 and ρf = 1000 kg·m−3, respectively. For a closed container, nonslip solid walls conditions are prescribed at the surrounded boundaries.Figure 11
Sketch of the computational domain for a circular particle settling in a closed container.In this example, the sedimentation process of the particle is simulated. At the beginning, the particle begins to settle by gravity, which is more than buoyant force. With the increase of the settling velocity, the drag force on the particle increases gradually and the settling acceleration decreases gradually. Finally, the settling velocity reaches a constant value due to the equilibrium of three forces, gravity, buoyancy, and drag force. For the purpose of convergence test, a series of grid spacings are considered, as shown in Table6, with a series of properly chosen time steps to maintain the Courant number around 0.2.Table 6
A series of grid spacings with time steps for different cases.
Case 1 2 3 4 5 6 7 Grid spacing/m d/12.5 d/18.75 d/25 d/31.25 d/37.5 d/43.75 d/50 Time step/s 7.5e-4 5e-4 3.75e-4 3e-4 2.5e-4 2.14e-4 1.975e-4Figure12(a) gives the evolution of particle vertical velocity with time, and Figure 12(b) is the enlarged view of area A in Figure 12(a). It is clear that the computed results are convergent with the decrease of grid spacing. For the purpose of comparison, Kang’s method [29] and Glowinski’s method [30] are also covered. It can be found that the method presented in this paper has closer results in all cases to Kang’s and especially Glowinski’s methods than Gallier’s (configuration is the same as case 2). Figure 13(a) shows the evolution of particle vertical position with time, and Figure 13(b) is the enlarged view of area A in Figure 13(a). Similar to Figure 12, the results are convergent with the decrease of grid spacing, and all of them are closer to Kang’s method than Gallier’s method. Glowinski did not present the evolution of particle vertical position with time in his paper, but the vertical velocity obtained by this method is consistent with that in his paper, so the vertical position predicted in this paper should not be very different from that calculated by Glowinski’s method.Figure 12
Evolution of particle vertical velocity with time (a). (b) is the enlarged view of area A in (a).
(a) (b)Figure 13
Evolution of particle vertical position with time (a). (b) is the enlarged view of area A in (a).
(a) (b)The simulation results illustrate that the method presented in this paper can deal with the moving boundary problems very well.


### 4.5. In-Line Oscillation of a Circular Cylinder

To further verify the feasibility of this method in dealing with moving boundary problems, an in-line oscillation of a circular cylinder is selected as an example for study. As shown in Figure14, the computational domain is 30d×20d (d=1 cm is the diameter of the cylinder), the circumference is set to be the second boundary condition with zero velocity gradient, and the initial position of the cylinder is (0, 0). The fluid is stationary at the initial time, with the density, ρf=1000 kg·m−3. The cylinder is periodically reciprocating according to the following expression(21)Xct=-Asin⁡2πtTwhere Xc(t) is the abscissa of the centre of the cylinder; A and T are amplitude and period, respectively. Reynolds number, Re, and Keulegan-Carpenter number, KC, which determine the characteristics of the flow, are calculated as follows(22)Re=Uc,maxdν(23)KC=Uc,maxTdwhere Uc,max is the maximum velocity of the cylinder during oscillation and ν is the coefficient of kinematic viscosity of the fluid.Re = 100,KC= 5 are selected in this example, similar to the experimental conditions selected by Dutsch et al. [18]. The grid spacing is d /20 in size.Figure 14
Sketch of the computational domain for a circular cylinder oscillating in-line in a quiescent fluid.Figure15 shows the stream force coefficient obtained by different methods within a period after stabilization. Cx is calculated according to the following formula(24)Cx=2FxρfUc,max2dwhere Fx is the horizontal resultant force on the cylinder. Fx is calculated according to the following formula(25)Fx=-ρf∫Pλxdx=-ρf∫DPIxλxdxFor the purpose of comparison, the results from Gallier’s method and experimental data from Dutsch are also given in the figure. It can be seen from the figure that although Gallier’s method provides smoother calculation results, the absolute values at the crest and trough are larger than the experimental values; on the contrary, in spite of small oscillation of the results from this method at 1.25T and 1.75T, the results are closer to the experimental values mainly due to the fact that the particle indicator function I(x) of Gallier’s method is absolutely smooth at the interface of two phases but cannot guarantee the conservation of the integral term in Formula (25), resulting in the smooth change of Cx with time, but inaccurate. However, I(x) from this method is absolutely conservative, but there are some discontinuous changes at the interface with the motion of the cylinder, resulting in small oscillations at 1.25T and 1.75T when the cylinder switches its direction. Yet the results obtained by this method agree with Dutsch’s experimental results.Figure 15
Comparison of the streamwise force coefficient.Figure16 shows the vorticity around the cylinder at different phase times within a period after stabilization. It can be spotted that the results in this paper are in good agreement with those from Dutsch et al.Figure 16
Vorticity contour (negative values dashed) at different phase time: (a) numerical results of Dutsch [18]; (b) the present numerical results.
(a) (b)Figures17, 18, and 19 show the profiles of the component u in the x direction and v in the y direction at x = -6, 0, 6, and 12 mm at the phase time t=nT+T/2, t=nT+7T/12, and t=nT+11T/12, respectively. For the purpose of comparison, the experimental data and numerical results from Dutsch are also presented. It can be seen that the results in this paper are in good agreement with the data from Dutsch.Figure 17
Comparison of the velocity components at four cross sectionsx =(a) -6, (b) 0, (c) 6, (d) 12mm at t=nT+T/2. Solid line: the present numerical results; ●: experimental data [18]; ▲: numerical results [18].
(a) (b) (c) (d)Figure 18
Comparison of the velocity components at four cross sectionsx =(a) -6 mm, (b) 0 mm, (c) 6 mm, (d) 12 mm at t=nT+7T/12. Solid line: the present numerical results; ●: experimental data [18]; ▲: numerical results [18].
(a) (b) (c) (d)Figure 19
Comparison of the velocity components at four cross sectionsx =(a) -6 mm, (b) 0 mm, (c) 6 mm, (d) 12 mm at t=nT+11T/12. Solid line: the present numerical results; ●: experimental data [18]; ▲: numerical results [18].
(a) (b) (c) (d)The simulation results from this example further indicate that the method in this paper is feasible in dealing with the moving boundary problems.


## 4.1. Particle Volume Evaluation

In this example, the accuracy of the volume calculation of circular and spherical particles on 2D and 3D structured background mesh at different subdivision levels. The circular and spherical particles have the diameter,d = 1 m, the grid spacing is set to d/10, d/20, d/30, and d/40, respectively, and the subdivision level is set to 1, 2, 3, and 4, totally 16 cases. Tables 1 and 2 give the particle volumes ∑IxVcell and relative errors by two indictor functions in 2D and 3D, respectively. The exact volumes for 2D circle and 3D volume are 0.785398 and 0.523599, respectively. It can be observed from the tables that, with a grid fining, our method produces much smaller relative errors than the approximate Heaviside function method, showing interaction-volume based particle indicator function I(x) is more accurate.Table 1
Comparison of results and errors of the particle volume computation in 2D.
Grid Spacing Volume and Errors Gallier [14] Subdivision level 1 2 3 4 d/10 Volume 0.793864 0.785212 0.785347 0.785384 0.785394 Error 1.08e-2 2.37e-4 6.50e-5 1.80e-5 5.02e-6 d/20 Volume 0.788934 0.785353 0.785386 0.785395 0.785397 Error 4.50e-3 5.92e-5 1.56e-5 4.33e-6 1.21e-6 d/30 Volume 0.787999 0.785379 0.785393 0.785397 0.785398 Error 3.31e-3 2.50e-5 6.90e-6 1.91e-6 5.31e-7 d/40 Volume 0.787515 0.785387 0.785395 0.785397 0.785398 Error 2.70e-3 1.41e-5 3.88e-6 1.07e-6 2.96e-7Table 2
Comparison of results and errors of the particle volume computation in 3D.
Grid spacing Volume and Errors Gallier [14] Subdivision level 1 2 3 4 d/10 Volume 0.529842 0.523251 0.523503 0.523572 0.523591 Error 1.19e-2 6.64e-4 1.83e-4 5.09e-5 1.42e-5 d/20 Volume 0.527607 0.523512 0.523575 0.523592 0.523597 Error 7.66e-3 1.66e-4 4.60e-5 1.28e-5 3.56e-6 d/30 Volume 0.525825 0.52356 0.523588 0.523596 0.523598 Error 4.25e-3 7.39e-5 2.04e-5 5.66e-6 1.58e-6 d/40 Volume 0.524104 0.523577 0.523593 0.523597 0.523598 Error 9.65e-4 4.15e-5 1.15e-5 3.19e-6 8.94e-7It can be found from Tables1 and 2 that the relative error is less than 5e-5 if the subdivision level is greater than or equal to 2 for the grid spacing with d/20. The subdivision level is set to 2 in the later tests, which is a good compromise between computational load and accuracy.Figure5 shows the evolution of particle volume with time obtained by two methods when the particle moves on the background mesh in 2D and 3D, respectively. It is not difficult to find that approximated Heaviside function method is not conservative because the results significantly oscillate with the particle motion, while our method in this paper is conservative and the results are stable.Figure 5
Evolution of predicted volume of a moving particle with time: (a) 2D; (b) 3D. solid line: the present numerical results; dash line: numerical results with Gallier’s method [14].
(a) (b)


## 4.2. Flow around a Cylinder

The flow around a cylinder is a classical hydrodynamic problem. A large number of physical and numerical experiments have been carried out previously. The accuracy and validity of this method in the stationary boundary are verified by the flow around a circular cylinder with different Reynolds numbers. In this example, the computational domain is a40d×20d rectangular domain (diameter of the cylinder: d = 0.025 m). The centre of the cylinder is located at (10d, 10d), as shown in Figure 6. The fluid flows into the left inlet and out of the right outlet, with the inflow velocity, u (different for different Re), fluid density, ρf = 1.293 kg·m−3, and kinetic viscosity coefficient, ν = 1.5×10−5 m2·s−1. The inlet flow velocity is specified, and the pressure is zero-gradient, while at the outlet, Neumann boundary condition is specified for velocity and Dirichlet boundary condition for pressure. The walls at top and bottom are set as the symmetric boundary condition. The Reynolds numbers,Re(= ρfd|u| / ν), of 10, 20, 40, 100, and 200, are selected as test conditions in this example, with the corresponding inlet velocities of 6×10−3 m·s−1, 1.2×10−2 m·s−1, 2.4×10−2 m·s−1, 6×10−2 m·s−1, and 0.12 m·s−1, respectively. In order to reduce the computational cost, a refinement mesh is used within a 5d×5d square, whose centre is the same as the cylinder’s and edges are parallel to the computational domain’s, and d/5×d/5 meshes are used at the rest of the area in this example. Three levels of grid spacings are employed: d/10 for Case 1, d/20 for Case 2, and d/30 for Case 3, respectively. The corresponding time steps are set so that the Courant number Co (= U0d/Δt) maintains around 0.2.Figure 6
Sketch of the computational domain for flow around a circular cylinder.Table3 gives the drag coefficients Cd, at different Re and grid spacings, at the ultimately stable or cyclically stable flow around a circular cylinder. The drag coefficients are obtained from the following equation.(20)Cd=2Fhydρfdu2The computed Cd are convergent with the decrease of grid spacing. For the purpose of comparison, this table also gives the numerical results from the method in [14] (Gallier) and the experimental and numerical results in [21–23]. It can be seen from the table that our results are better almost in all cases than Gallier’s method and closer to the experimental and calculated results in the references.Table 3
Comparison of the present computation of drag coefficient (Cd) with other results for flow past a cylinder.
R e C d Present Gallier et al. [14] Nishioka et al. [21] Rajani et al. [22] Takami et al. [23] Case1 Case2 Case3 (Numerical) (Experimental) (Numerical) (Numerical) 10 2.991 3.037 3.049 3.116 2.74 2.839 2.754 20 2.126 2.164 2.174 2.236 1.95 2.059 2.003 40 1.573 1.604 1.611 1.668 1.49 1.539 1.536 100 1.302 1.321 1.328 1.376 — 1.335 — 200 1.246 1.335 1.381 1.402 — 1.337 —Figure7 shows the streamlines at Re of 10, 20, and 40 of Case 3 when the flow becomes eventually stable. It can be seen that stable recirculation zones (whirlpool vortex) are generated behind the cylinders in all three cases, and the length, Lw, of the recirculation zone increases with the increasing Re. Table 4 shows the comparison of the length, Lw/d (d is the diameter of the cylinder), of recirculation zone calculated by our method at Re of 20 and 40 for different grid spacings with those in [23–25]. The calculated results are convergent with the refinement of the mesh and in good agreement with those in the references.Table 4
The non-dimensional length of recirculation zone (Lw/d) at Re = 20, 40.
R e L w / d Present Takami et al. [23] Silva et al. [24] Ye et al. [25] Case1 Case2 Case3 (Numerical) (Numerical) (Numerical) 20 0.95 1.05 1.07 0.932 1.05 0.92 40 2.03 2.25 2.31 2.325 2.54 2.27Figure 7
Streamlines for flow around a circular cylinder of Case 3 atRe = (a) 10; (b) 20; (c) 40.
(a) (b) (c)Figure8 shows the velocity contour of cyclically stable flow at Re of 100 and 200 of Case 3. In both cases, the whirlpool vortex behind the cylinder begins to regularly shed and forms a regular vortex street. Table 5 shows the comparison of Strauhal number (St=fd / u, where f is the vortex shedding frequency) calculated by this method in both cases with those in [7, 26, 27]. The calculated results are converted with the decrease of grid spacing and are in good agreement with those in the references.Table 5
Comparison of the numerical predicted Strouhal number (St) with different methods for flow around a circular cylinder.
R e S t Present Liu et al. [26] Uhlmann [7] Tseng and Ferziger [27] Case1 Case2 Case3 100 0.162 0.165 0.166 0.164 0.169 0.164 200 0.189 0.193 0.194 0.192 — —Figure 8
Velocity contour for flow around a circular cylinder of Case 3 atRe = (a) 100, (b) 200.
(a) (b)The simulation results in this example indicate that, for the static boundary problem, the method in this paper is of good accuracy in a wide range ofRe.


## 4.3. Single Spherical Particle in a Creeping Shear Flow

In this section, a single force-free torque-free sphere in a linear creeping shear flow is considered to verify the feasibility of this method in dealing with sphere rotation. A spherical particle of diameterd=0.1 m is suspended at the centre of a cubic Couette cell, as shown in Figure 9. The particle is located at x = 0 and the cell size is L=10d. A velocity of γL/2 and -γL/2 in the flow direction is assigned to the top and bottom boundaries (the dark ones), respectively, where γ=6.71×10−3s−1 is the shear rate. And for the rest of the boundaries, periodic conditions are used. The kinematic viscosity of the flow is ν = 10−3 m2s−1. Three levels of grid spacings are used: Δ=d/5 for Case 1, Δ=d/7.5 for Case 2, and Δ=d/10 for Case 3. And the corresponding time steps are 0.01s, 0.005s, and 0.004s, respectively.Figure 9
Sketch of the computational domain for a single spherical particle in creeping shear flow.In a creeping shear flow, the effect of convection term is so little that it can be neglected, and unsteady Stokes equations are solved with respect to time-stepping. Eventually, the whole system reaches a steady state, and the particle is rotating in a constant angular velocity, which is available of analytical solutions [14].Figure10 shows a comparison between the profile of the computed nondimensional velocity 2u2/γd on different grid spacings on the centreline (x2=0) of the shear plane (x3=0), where the subscript represents the direction of the flow, and the analytical solution [28]. For the coarsest mesh, Case 1, the computed results near the outside region of the particle surface are slightly lower than the analytical solution. With the refinement of the mesh, the results become more and more close to the analytical solution both inside and outside the particle. The linear profile inside the particle (2|x1|/d<1) is due to the rigid body motion constraint of the particle rotation and a velocity 2|u2|/γd=0.5 is found at the surface of the particle (2|x1|/d=1). The rotational velocities Ω3/γ of all three cases are -0.4990, -0.4995, and -0.4996, respectively, and gradually approach the theoretical value -0.5, The particle translational velocity is also solved and turns out to be of the order of 10−8, which is acceptable with respect to the expected value 0.Figure 10
Evolution of the centreline fluid velocity2u2/γd with position 2x1/d for a single sphere located at the origin (x=0).The simulation results in this case show that this method is feasible for the simulation of spherical particle rotation.


## 4.4. Sedimentation of a Circular Particle in a Closed Container

The sedimentation process of a circular particle in a closed container is simulated in this case to verify the feasibility of this method in dealing with moving boundary problems. The computational model is shown in Figure11. The computational domain is 8d×24d (diameter of a particle: d = 2.5 × 10−3 m), and the initial position of the particle is (4d,8d). The fluid is stationary at the initial time. The circular particle settles by gravity, whose density is ρp = 1250 kg·m−3. The dynamic viscosity and density of the fluid are μ = 0.01 kg·m−1 and ρf = 1000 kg·m−3, respectively. For a closed container, nonslip solid walls conditions are prescribed at the surrounded boundaries.Figure 11
Sketch of the computational domain for a circular particle settling in a closed container.In this example, the sedimentation process of the particle is simulated. At the beginning, the particle begins to settle by gravity, which is more than buoyant force. With the increase of the settling velocity, the drag force on the particle increases gradually and the settling acceleration decreases gradually. Finally, the settling velocity reaches a constant value due to the equilibrium of three forces, gravity, buoyancy, and drag force. For the purpose of convergence test, a series of grid spacings are considered, as shown in Table6, with a series of properly chosen time steps to maintain the Courant number around 0.2.Table 6
A series of grid spacings with time steps for different cases.
Case 1 2 3 4 5 6 7 Grid spacing/m d/12.5 d/18.75 d/25 d/31.25 d/37.5 d/43.75 d/50 Time step/s 7.5e-4 5e-4 3.75e-4 3e-4 2.5e-4 2.14e-4 1.975e-4Figure12(a) gives the evolution of particle vertical velocity with time, and Figure 12(b) is the enlarged view of area A in Figure 12(a). It is clear that the computed results are convergent with the decrease of grid spacing. For the purpose of comparison, Kang’s method [29] and Glowinski’s method [30] are also covered. It can be found that the method presented in this paper has closer results in all cases to Kang’s and especially Glowinski’s methods than Gallier’s (configuration is the same as case 2). Figure 13(a) shows the evolution of particle vertical position with time, and Figure 13(b) is the enlarged view of area A in Figure 13(a). Similar to Figure 12, the results are convergent with the decrease of grid spacing, and all of them are closer to Kang’s method than Gallier’s method. Glowinski did not present the evolution of particle vertical position with time in his paper, but the vertical velocity obtained by this method is consistent with that in his paper, so the vertical position predicted in this paper should not be very different from that calculated by Glowinski’s method.Figure 12
Evolution of particle vertical velocity with time (a). (b) is the enlarged view of area A in (a).
(a) (b)Figure 13
Evolution of particle vertical position with time (a). (b) is the enlarged view of area A in (a).
(a) (b)The simulation results illustrate that the method presented in this paper can deal with the moving boundary problems very well.


## 4.5. In-Line Oscillation of a Circular Cylinder

To further verify the feasibility of this method in dealing with moving boundary problems, an in-line oscillation of a circular cylinder is selected as an example for study. As shown in Figure14, the computational domain is 30d×20d (d=1 cm is the diameter of the cylinder), the circumference is set to be the second boundary condition with zero velocity gradient, and the initial position of the cylinder is (0, 0). The fluid is stationary at the initial time, with the density, ρf=1000 kg·m−3. The cylinder is periodically reciprocating according to the following expression(21)Xct=-Asin⁡2πtTwhere Xc(t) is the abscissa of the centre of the cylinder; A and T are amplitude and period, respectively. Reynolds number, Re, and Keulegan-Carpenter number, KC, which determine the characteristics of the flow, are calculated as follows(22)Re=Uc,maxdν(23)KC=Uc,maxTdwhere Uc,max is the maximum velocity of the cylinder during oscillation and ν is the coefficient of kinematic viscosity of the fluid.Re = 100,KC= 5 are selected in this example, similar to the experimental conditions selected by Dutsch et al. [18]. The grid spacing is d /20 in size.Figure 14
Sketch of the computational domain for a circular cylinder oscillating in-line in a quiescent fluid.Figure15 shows the stream force coefficient obtained by different methods within a period after stabilization. Cx is calculated according to the following formula(24)Cx=2FxρfUc,max2dwhere Fx is the horizontal resultant force on the cylinder. Fx is calculated according to the following formula(25)Fx=-ρf∫Pλxdx=-ρf∫DPIxλxdxFor the purpose of comparison, the results from Gallier’s method and experimental data from Dutsch are also given in the figure. It can be seen from the figure that although Gallier’s method provides smoother calculation results, the absolute values at the crest and trough are larger than the experimental values; on the contrary, in spite of small oscillation of the results from this method at 1.25T and 1.75T, the results are closer to the experimental values mainly due to the fact that the particle indicator function I(x) of Gallier’s method is absolutely smooth at the interface of two phases but cannot guarantee the conservation of the integral term in Formula (25), resulting in the smooth change of Cx with time, but inaccurate. However, I(x) from this method is absolutely conservative, but there are some discontinuous changes at the interface with the motion of the cylinder, resulting in small oscillations at 1.25T and 1.75T when the cylinder switches its direction. Yet the results obtained by this method agree with Dutsch’s experimental results.Figure 15
Comparison of the streamwise force coefficient.Figure16 shows the vorticity around the cylinder at different phase times within a period after stabilization. It can be spotted that the results in this paper are in good agreement with those from Dutsch et al.Figure 16
Vorticity contour (negative values dashed) at different phase time: (a) numerical results of Dutsch [18]; (b) the present numerical results.
(a) (b)Figures17, 18, and 19 show the profiles of the component u in the x direction and v in the y direction at x = -6, 0, 6, and 12 mm at the phase time t=nT+T/2, t=nT+7T/12, and t=nT+11T/12, respectively. For the purpose of comparison, the experimental data and numerical results from Dutsch are also presented. It can be seen that the results in this paper are in good agreement with the data from Dutsch.Figure 17
Comparison of the velocity components at four cross sectionsx =(a) -6, (b) 0, (c) 6, (d) 12mm at t=nT+T/2. Solid line: the present numerical results; ●: experimental data [18]; ▲: numerical results [18].
(a) (b) (c) (d)Figure 18
Comparison of the velocity components at four cross sectionsx =(a) -6 mm, (b) 0 mm, (c) 6 mm, (d) 12 mm at t=nT+7T/12. Solid line: the present numerical results; ●: experimental data [18]; ▲: numerical results [18].
(a) (b) (c) (d)Figure 19
Comparison of the velocity components at four cross sectionsx =(a) -6 mm, (b) 0 mm, (c) 6 mm, (d) 12 mm at t=nT+11T/12. Solid line: the present numerical results; ●: experimental data [18]; ▲: numerical results [18].
(a) (b) (c) (d)The simulation results from this example further indicate that the method in this paper is feasible in dealing with the moving boundary problems.


## 5. Conclusions

A fictitious domain method based on the particle indicator function from intersection volume is developed in this paper for dealing with moving boundary problems. In this method, the background mesh cells across the interface between fluid and particles are located and subdivided into tetrahedrons. The volume fraction of the particle domain within each cell is accurately computed to achieve high precision of integration calculation in the particle domain, improving the accuracy of the whole method. The integral value of the particle domain calculated by this function is more accurate, which greatly improves the overall calculation accuracy of the fictitious domain method. Volume evaluation based on cells subdivision techniques shows that the technology can accurately calculate the volume of static particles and guarantee the conservation of the calculated integral on the moving particle. The test results of the flow around a cylinder, a single sphere in a creeping shear flow, the sedimentation of a circular particle in a closed container, and in-line oscillation of a circular cylinder show that this fictitious domain method has good accuracy in dealing with stationary boundary and moving boundary problems. This method has better performance in simulation of stationary and moving particle-fluid motion than the method using approximate Heaviside function. This method is more accurate and conservative. It can be adopted to deal with moving and stationary boundary problems.



---

*Source: 5450313-2019-06-26.xml*