A three-dimensional robust volume-of-fluid solver based on the adaptive mesh refinement

Xin Zhao

Department of Mechanics, Beijing Institute of Technology, Beijing 10 0 081, China

Keywords: Adaptive mesh refinement Volume of fluid method Surface tension Interfacial flow Projection method

ABSTRACT The present study provides a three-dimensional volume-of-fluid method based on the adaptive mesh refinement technique.The projection method on the adaptive mesh is introduced for solving the incom- pressible Navier-Stokes equations.The octree structure mesh is employed to solve the flow velocities and the pressure.The developed solver is applied to simulate the deformation of the cubic droplet driven by the surface tension without the effect of the gravity.The numerical results well predict the shape evolution of the droplet.

Many engineering problems have multiscale properties.The presence of the interfacial flow makes the numerical solution for such problem even harder [1] .To obtain the efficient and robust numerical method for solving the multiscale flow with the sharp phase interface, the adaptive mesh refinement technique is intro- duced.With the dynamic adaptive mesh method, one can sig- nificantly reduce the computational cost and improve the com- putational accuracy.Up to now, more and more commercial and open-sourced computational fluid dynamics codes [ 2 , 3 ] involve the adaptive mesh refinement technique to speed up the simulation, especially for interfacial flows.The present study provides a simpli- fied adaptive volume-of-fluid solver to achieve both the efficiency and the accuracy.

The governing equations for the present study are the incom- pressible Navier-Stokes equations:

Here,uis the fluid velocity vector,ρis the fluid density,pis the pressure,μisthedynamic viscosity, andisthe

deformation tensor.For surface tension term,σis the surface ten- sion coefficient,κis the curvature of the interface,δis the Dirac delta function of the interface, andnis the unit normal vector of the interface.

For the gas-liquid two phase flow, we use the volume fraction functionαof the liquid to identify the gas domain (α= 0), the liquid domain (α= 1), and the phase interface (0<α<1).With the volume fraction function we could represent the fluid density and dynamic viscosity as the following:

Here,ρgandρlare the gas density and the liquid density, andμgandμlare the gas dynamic viscosity and the liquid dynamic viscosity.The governing equation for the volume fraction function of the liquid is shown as the followings:

To solve the gas-liquid two phase flow on the octree structure mesh, we need to generate the multiscale adaptive mesh.An ex- ample to represent the quadtree mesh in two dimensions is shown in Fig.1 , which is reproduced for the Ref.[2] .In the mesh genera- tion process, the adaptive mesh refinement criterion is set to see if the volume fraction and the vorticity reach the threadhold values.This means in a grid cell if the norm of the gradient of the liquid volume fraction | ∇α|>εαor the norm of the vorticity |ω|>εω, we need to generate finer mesh in this grid cell.εαandεωare the threadhold values given in advance.To achieve such meshing pro- cess, we need first to define different cells used in the algorithm [4] :

Root cells :The base of the cell tree.The cells on the level 0 as shown in Fig.1 .

Children cells :The direct descendants of a cell.One cell has eight children cells in the octree mesh of three dimensions and four children cells in the quadtree mesh of two dimensions as shown in Fig.1 .

Parent cell :The direct ancestor of a cell.

Leaf cells :The cells do not have children.

Neighbour cells :The neighbours of a given cell.There are six neighbour cells in three dimensions and four neighbour cells in two dimensions.

Cell level :Record the level of a cell as shown in Fig.1 .

For convenience, the cell level difference of a given cell and its neighbour cells is 0 or 1.

With the adaptive mesh, we could construct the numerical schemes to solve the governing equations.Firstly, we introduce the numerical schemes to solve Navier-Stokes equations.For the advec- tion term, the second-order upwind explicit scheme is employed, and the discretization for the spatial derivative is constructed with the cell value and its neighbour cells as the Refs.[ 5 , 6 ].For exam- ple, the spatial derivative alongxaxis as shown in Fig.2:

Here,his the grid size of the given cellC,uCis the value of the cellC,uR is the value of the right neighbour of the cellC, anduL is the value of the left neighbour of the cellC.Also,uL is the value of the

parent cell of cellsuL1 ,uL2 ,uL3 , anduL4 withuR1 ,uR2 ,uR3 , anduR4 are the virtual sub-cells of the celluR , and the value ofuR1 is obtained by the linear interpolation method.With this finite difference method on the octree mesh, we could further construct the spatial discretization scheme for the viscous terms by using the central difference scheme.With the projection method, we could construct the multilevel method to solve the pressure Poisson equation as the Ref.[2] .

With the spatial discretization on the adaptive mesh, we could solve the fluid velocities and the pressure.Based on the distribu- tion of the fluid velocities and the liquid volume fraction, we could calculate the new distribution of the liquid volume fraction for the next time step.The advection of the liquid volume fraction is sim- ulated by the geometric volume of fluid method as the Refs.[ 1 , 3 ].Figure 3 shows the sketch on the calculation for the advection ofthe liquid volume fraction on the two dimensional quadtree mesh.Since we have set the cell level difference of a given cell and its neighbour cells is 0 or 1, we only have two cases for calculating the advection of the liquid volume fraction as shown in Fig.3 .In Fig.3 , we could find the change of the liquid volume on the given cell is the sum on the area of the dark blue parts.Based on the geometric volume of fluid method, the area of the dark blue parts can be calculated exactly.

With the obtained distribution of the liquid volume fraction, we could calculate the surface tension term in the momentum equa- tions.Firstly, the unit normal vector of the phase interface can be calculated by the following formula:

The Dirac delta function of the phase interface is defined by

The curvature of the phase interface can be calculated as

Therefore, the final expression for the surface tension term is

In Eq.(11) , the gradient and divergence operators can be cal- culated by the central finite difference method.The procedure is similar to the solving method for momentum equations.

With the above numerical methods, we could develop the solver to simulate the interfacial flow driven by the surface ten- sion.To test the solver, the deformation of a cubic water droplet driven by the surface tension in the air is simulated.The proper- ties of the water and the air are given in Table 1 .The initial volume of the cubic water droplet is 0.125 mm3, and the theoretical radius of the droplet at the finial stable state is 0.31 mm.

Table 1.The physical parameters of the water and the air.

In the simulation, the highest level is set as 4.The grid size for the level 0 is 0.125 mm, and the grid size for the level 4 is 0.0078125 mm.The grid number for the level 0 is 16 ×16 ×16 .The estimated total grid number for the uniform grid size 0.0078125 mm is about 17 million.The maximum total grid num- ber for the adaptive simulation is 0.3 million, and it is about 0.18% of the estimated total grid number for the uniform grid size.

Figure 4 shows the numerical results on the evolution of the iso-surface of the droplet and the adaptive mesh mapped on the two dimensional slice.Figure 5 shows the evolution of the con- tours of the droplet on the two dimensional slice.For the results, we could find the shape of the cubic droplet change to sphere fast at the starting stage.However, in the following stage, the droplet tends to the cubic shape due to the inertial effect, and the droplet is not stable.After longer time, the droplet starts to change to the spheric shape gradually.Finally, the droplet obtains the stable spheric shape att= 1 ms .The whole process satisfies the law of physics..

To further evaluate the accuracy of the solver, we plot the com- parison between the theory and the simulation on the radius of the droplet (or the inverse of the curvature of the interface) for different positions at the final stable state(Fig.6).We can find that the simulation result is very close to the theoretical value, but the radius of the droplet for the simulation result is slightly smaller than the theoretical value.This means the total volume of the droplet slightly decreases during the deformation process.However, the error is acceptable.

In sum, we introduce a simplified three-dimensional volume- of-fluid solver based on the adaptive mesh refinement technique.The phase interface tracking method and the fluid solving numeri- cal method are constructed based on the dynamic octree adaptive mesh.The adaptive solver is employed to simulate the deformation of a cubic droplet driven by the surface tension.The numerical re- sults on the radius of the droplet at the final stable state are close to the theoretical value.This comparison demonstrates the accu- racy of the adaptive solver.

Declaration of Competing Interest

The authors declare that they have no known competing finan- cial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work was supported by the National Natural Science Foun- dation of China (No.41776194).