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.
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 domain contains fluid domain and particle domain , i.e., = . The particle domain contains rigid particles with the density of , and the interface between particle and fluid is . The fluid is incompressible viscous Newtonian fluid with the density of , the dynamic viscous coefficient of , and the kinematic viscous coefficient of . The governing equation of the fictitious domain model is expressed as followswhere and 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; , , , and are the mass, moment of inertia, translational velocity, and rotational angular velocity of a single particle, respectively; is the position vector at the mass centre of a particle; is gravitational acceleration; and 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 aswhere is the diagonal value of the coefficient matrix obtained by discretization of (1), and is expressed aswhere is the implicit contribution coefficient of the neighbour cells to the current cell during discretization of (1) and 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 yieldingThe velocity field obtained by (8) should satisfy the continuity equation. Substituting (8) into (2) yieldsEquation (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 usingwhere is the interpolation of value of the body centre, , on the face centre, is the area vector of the cell face, and is the volume flow rate on the cell face.
2.2.2. Particle Subproblem
For the particle subproblem, the following equations are solved firstto obtain the translational velocity, , and rotational angular velocity, , at the new time step. Then, use the equation belowto update the fluid velocity in the particle covering domain and use this velocity to update the pseudo body force, .Finally, the second-order Adams-Bashforth difference scheme is used to discretize the displacement equation and update the particle position.For 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 where is the indicator function, which is expressed by the following expressionwhere is the level-set function and is expressed aswhere and are the mass centre and radius of the particle, respectively; is a negative value inside the particle and a positive value outside the particle. The approximate Heaviside function can be expressed aswhere ∆ is the grid spacing and the constant =5 guarantees the smooth width within a grid. Eventually, the particle indicator function 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 > 0 are set to = 0 first, and the rest are set to = 1. The cells near the solid boundary are located using the level set functions and . The subscripts and represent the current cell and the neighbour cell, respectively. If < 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.
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.
Figure 2(a) is a polyhedral cell centred at . 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, , is taken separately, as shown in Figure 2(b), and and are the face centre and area vector of the pentagon , respectively. Three tetrahedrons, , , and , can be obtained by connecting , , and . 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).
(a)
(b)
(c)
Secondly, a tetrahedron can be decomposed into eight small tetrahedrons further to improve the accuracy.
The tetrahedron shown in Figure 3(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.
(a)
(b)
(c)
(d)
3.3. Calculation of Particle Indicator Function 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, and , of a tetrahedron across the interface, if , the edge containing and intersects with the interface, and the intersection point can be determined by computing = 0 on that edge.
Then a tetrahedron is divided into two parts by connecting any two of these intersection points along the tetrahedron surfaces (Figure 4(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.
(a)
(b)
(c)
The solid volume in a cell near the solid boundary is the sum of the volumes, ∑, 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 ∑ / is adopted as the particle indicator function of the cell.
Notice that we assume the interface within the tetrahedron is flat (Figure 4(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, = 1 m, the grid spacing is set to /10, /20, /30, and /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 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 is more accurate.
It can be found from Tables 1 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 /20. The subdivision level is set to 2 in the later tests, which is a good compromise between computational load and accuracy.
Figure 5 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.
(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 a rectangular domain (diameter of the cylinder: = 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, (different for different ), fluid density, = 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 (= / ), 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 square, whose centre is the same as the cylinder’s and edges are parallel to the computational domain’s, and meshes are used at the rest of the area in this example. Three levels of grid spacings are employed: for Case 1, for Case 2, and for Case 3, respectively. The corresponding time steps are set so that the Courant number (= ) maintains around 0.2.
Table 3 gives the drag coefficients , at different and grid spacings, at the ultimately stable or cyclically stable flow around a circular cylinder. The drag coefficients are obtained from the following equation.The computed 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.
Figure 7 shows the streamlines at 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, , of the recirculation zone increases with the increasing . Table 4 shows the comparison of the length, ( is the diameter of the cylinder), of recirculation zone calculated by our method at 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.
(a)
(b)
(c)
Figure 8 shows the velocity contour of cyclically stable flow at 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 ( / , where 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.
(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 of .
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 diameter 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 . A velocity of and - 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: for Case 1, for Case 2, and for Case 3. And the corresponding time steps are 0.01s, 0.005s, and 0.004s, respectively.
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].
Figure 10 shows a comparison between the profile of the computed nondimensional velocity on different grid spacings on the centreline () of the shear plane (), 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 () is due to the rigid body motion constraint of the particle rotation and a velocity is found at the surface of the particle (). The rotational velocities 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.
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 Figure 11. The computational domain is (diameter of a particle: = 2.5 × 10−3 m), and the initial position of the particle is . The fluid is stationary at the initial time. The circular particle settles by gravity, whose density is = 1250 kg·m−3. The dynamic viscosity and density of the fluid are = 0.01 kg·m−1 and = 1000 kg·m−3, respectively. For a closed container, nonslip solid walls conditions are prescribed at the surrounded boundaries.
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 Table 6, with a series of properly chosen time steps to maintain the Courant number around 0.2.
Figure 12(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.
(a)
(b)
(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 Figure 14, the computational domain is ( 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, =1000 kg·m−3. The cylinder is periodically reciprocating according to the following expressionwhere is the abscissa of the centre of the cylinder; and are amplitude and period, respectively. Reynolds number, , and Keulegan-Carpenter number, , which determine the characteristics of the flow, are calculated as followswhere 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 /20 in size.
Figure 15 shows the stream force coefficient obtained by different methods within a period after stabilization. is calculated according to the following formulawhere is the horizontal resultant force on the cylinder. is calculated according to the following formulaFor 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 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 with time, but inaccurate. However, 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 16 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.
(a)
(b)
Figures 17, 18, and 19 show the profiles of the component in the direction and in the direction at = -6, 0, 6, and 12 mm at the phase time , , and , 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The project supported by the National Science and Technology Major Project of the Ministry of Science and Technology of China (Grant No. 2016ZX05011001), the National Natural Science Foundation of China (Grant No. 21306145), and Nature Science Basic Research Plane in Shaanxi Province of China (Grant No. 2019JQ-335).