Implementation of the Level Set Method into OpenFOAM for Capturing the Free Interface in Incompressible Fluid Flows Bitan Shu, Frank Dammel, Peter Stephan e-mail:
[email protected] Chair of Technical Thermodynamics, Technische Universit¨at Darmstadt Petersenstr. 30, 64287 Darmstadt, Germany
Abstract The topic of this work is the implementation of the level set method into OpenFOAM for capturing the free interface motion in incompressible fluid flows. The incompressible fluid flow is calculated with mass conservation equation and momentum transport equation. The free interface is indicated with the zero level set of a smooth distance function. The level set method has two important issues: the reinitialisation of the level set function and the mass conservation. During the reinitialisation of the level set function, an upwind scheme and a WENO scheme are applied to calculate the gradient. To avoid the loss of mass, the method which was introduced by Chang et al. [3] is used. As test case, the terminal velocities of rising air bubbles in water are calculated. Additionally, detachment of a bubble from a wall is simulated.
1 Introduction Capturing the motion of a free interface in incompressible multi-phase fluid flows is a challenging research field of computational fluid dynamics. The motion of a free interface is a important physical phenomenon in many physically interesting problems, e.g. combustion and nucleate boiling. The level set method is becoming more and more popular to capture the motion of a free surface since its introduction by Osher and Sethian [9]. In the level set method, the interface is captured implicitly by the level set function embedded in the fluid field. The simplest choice for the level set function is the signed distance from the interface. In this way the level set function is smooth and continuous, its spatial derivatives can be accurately determined for the calculation of the interface curvature. As in the volume of fluid (VOF) method, there is no difficulty to handle the topological change of the interface during its evolution, so that the detachment of a bubble from a wall and the merging of bubbles can be simulated without special efforts. In general, the implementation of the level set
1
method is straightforward, as well as the extension from 2D simulation to 3D simulation. Another advantage of the level set method is the maintenance of the sharp interface which in the VOF method can only be remained by the complicated reconstruction of the interface [16]. The maintenance of the sharp interface is the major motivation for this work. A sharp interface is required, e.g. if the heat transfer has to be calculated during simulation of fluid flows with phase change [15]. Sussman et al. [14] were the first that applied the level set method to incompressible fluid flows. In their work two-phase fluid flow with a large density ratio of about 1000 could be successfully simulated. In this work the level set method is combined with the PISO method [5] for solving the mass conservation equation and the momentum transport equation. This method is then applied to simulate rising air bubble in water, a problem with a large density ratio. The second example, the detachment of a bubble, shows the ability of this method to handle the topological change of the interface.
2 Numerical Formulation 2.1 Governing Equations The mass conservation equation for incompressible fluid flow is ∇ · u = 0.
(1)
Using the continuum surface force model [2], the momentum transport equation for incompressible flow is ∂u + ∇u · u = −∇p + ρg + ∇ · µ∇(u + uT ) + σκδ(φ)∇φ. (2) ρ ∂t The level set function φ is defined as negative in the gas phase, positive in the liquid phase and 0 at the interface. The equation for the advection of the level set function φ with the flow field reads: ∂φ + u · ∇φ = 0 (3) ∂t In equation (2) σ is the surface tension and κ is the interface curvature which can be calculated with the following equation: κ=∇·
∇φ |∇φ|
(4)
δ(φ) depends on the φ-field: δ(φ) =
(
0 1 2 (1
+ cos(πφ/ǫ))/ǫ
2
if |φ| ≥ ǫ if |φ| < ǫ.
(5)
The δ-function is the derivative of the Heaviside-function H: if φ ≤ −ǫ 0 H= 1 if φ ≥ ǫ (φ + ǫ)/(2ǫ) + sin(πφ/ǫ)/(2π) if |φ| < ǫ.
(6)
The H-function is used to smooth the density and the viscosity at the interface over a width of ǫ, which can be chosen to be 1.5 grid width (h). Thus, the density and the viscosity over the whole fluid field are: ρ = ρg + (ρl − ρg )H (7) and µ = µg + (µl − µg )H.
(8)
2.2 Solution of the Governing Equations Equations (1) and (2) in section 2.1 are discretised and solved with the PISO method [5] in the framework of the finite volume method, which is already available in OpenFOAM [1, 6, 10]. The equation for the level set function (3) is modified as follows: ∂φ + ∇ · (uφ) = 0, ∂t
(9)
which holds because of equation (1). This equation is more suitable for discretisation with the finite volume method in OpenFOAM. The level set function φ is defined as a volScalarField in OpenFOAM, and the term ∇ · (uφ) can be discretised with the first order upwind scheme.
2.3 Reinitialisation of the Level Set Function The level set function is defined as the distance function from the interface. However, it will not remain the distance function during the advection because the velocity field is not uniform in general. For this reason it must be reinitialised after some time steps with the following equation: ∂φ = sign(φ0 )(1 − |∇φ|) (10) ∂τ Here, τ is the artificial time. The level set field at the beginning of the iteration is: φ(x, 0) = φ0 (x) The steady state solution of the equation (10) satisfies the condition |∇φ| = 1, so that it is the distance function from the interface. To solve equation (10), the first order upwind scheme can be used to determine the gradient of the φ-field during the iterations. In the upwind scheme the gradient is calculated
3
on both sides of the cells [14]. Therefore, the following variables are defined: a = Dx− φij ≡ ∇x− ≡ (φij − φi−1,j )/h b = Dx+ φij ≡ ∇x+ ≡ (φi+1,j − φi,j )/h c = Dy− φij ≡ ∇y− ≡ (φij − φi,j−1 )/h d = Dy+ φij ≡ ∇y+ ≡ (φi,j+1 − φi,j )/h,
(11)
the smoothed sign function
and
Sǫ (φij,0 ) = q
φij,0 φ2ij,0 + ǫ2
,
p 0 + 2 − 2 + 2 − 2 pmax((a ) , (b ) ) + max((c ) , (d ) ) − 1 if φij > 0 G(φ)ij = max((a− )2 , (b+ )2 ) + max((c− )2 , (d+ )2 ) − 1 if φ0ij < 0 . 0 otherwise
(12)
Everywhere holds x+ = max(x, 0) and x− = min(x, 0). The equations defined above can be extended to the 3D case in similar manner. Equation (10) is updated with the following equation: +1 0 N φN = φN ij − ∆τ Sǫ (φij )G(φij ). ij
(13)
The artificial time step ∆τ must be smaller than one grid width and is chosen to be 0.1 grid width in the test cases. The stop criterion for the iterations is: P +1 |φN − φN ij | |φN ij ij |<α < ∆τ h2 , (14) E= M where h is the grid width as above and α is an artificial finite thickness adjacent to the interface, e.g. a = 5h. M is the number of grid points where |φN ij | < α. While the upwind scheme above is only of the first order, approximation of the gradients with the 5th order WENO scheme is more accurate. With the 5th order WENO scheme, the gradients in equations (11) and (12) are calculated with the level set values in 5 adjacent grids for each grid face. Details of the WENO scheme are given in [8, 12]. Note that the gradients calculated in this step can also be used to solve equation (3).
2.4 Conservation of the Mass It is found that the loss of mass occurs during the solution of equations (3) and (10). To avoid the loss of mass various studies have been carried out. Russo [11] suggested a method, which focuses on the reinitialisation of the level set function, and works well in case of a poorly predefined level set field. However, the simplest way to avoid the loss of mass during the solution of equation (10) is just not to update the values of the level set function in cells at the interface during the iterations.
4
A greater amount of the mass loss occurs during the solution of equation (3), if it is solved with the first order upwind scheme. With the 5th order WENO scheme the amount of the mass loss is considerable smaller. A possible reason for mass loss is the unphysical smoothing of density [7]. The coupled level set and volume of fluid method was introduced by [13] to overcome this type of mass loss. As the first attempt in this work the simple method introduced by Chang et al. [3] is implemented. After the iterations for the solution of equation (3), the following equation is solved to reset the level set function: ∂φ + (M0 − M (τ ))(−P + κ)|∇φ| = 0 ∂τ φ(x, 0) = φ0 (x).
(15)
In this equation τ is the artificial time as in equation (10) above and P is a small positive constant, which is chosen to be 1 in this work. M0 is the total mass of one of the two phases at the beginning of the iterations for solving equation (15). After solving equation (15), M (τ ) becomes M0 with a small tolerance. The mass is added into cells which are not determined physically, the mass is only globally conserved. This method can not be utilised for the simulation of the bubble detachment in section (3.2) because it is not accurate enough.
2.5 Summary of the Major Procedures of the Program The major procedures are summarised as follows: 1. create the fields before the time-loop 2. begin the time-loop 3. solve the level set function 4. solve equation (15) to conserve the mass 5. solve equation (10) to reinitialise the level set function 6. update the transport properties, e.g. density and viscosity, according to the new level set function 7. solve the mass conservation equation (1) and the momentum transport equation (2) with the PISO method 8. write out the solution if required 9. go to the next time step
3 Numerical Tests 3.1 Rising Bubble in Fluid To test the program, the terminal velocities of rising air bubbles in pure water are computed. The size of the computational domain is x×y, with x is about 3R (R : initial bubble radius)
5
and y is about 10R. The computational domain is 2D axisymmetric. The properties of the fluid are listed in Tab. 1. At the beginning the velocity in the whole computational domain
variable unit density [kg/m3 ] viscosity [m2 /s] surface tension [N/m]
air (phase 1) water (phase 2) 1.2045 998.2 15.11e-06 1.004e-06 0.07
Table 1: Properties used in the simulation of terminal velocities is set to 0. In the course of time the bubble rises upwards because of the buoyancy force. Fig. 1 shows the position of the weight centre of the bubble with the initial radius 0.8 mm as a function of time. At the end it reaches the constant terminal velocity because the buoyancy force is equal to the resistibility, and the terminal velocity is determined from the slope of the curve in this stage. 0.01 0.009
Bubble position [m]
0.008 0.007 0.006 0.005 0.004 0.003 0.002 0.001
0
0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 t [s]
Figure 1: Position of the bubble with initial radius 0.8 mm vs. time For the case with 0.8 mm as initial bubble radius, where the computational domain is 3 × 10 mm2 , a test about the influence of the grid width is carried out. At first, the computational domain is discretised with the grid width h = 1/10 mm as the reference case. Then it is discretised with the grid width h = 1/8 mm and h = 1/12 mm. With the coarse grid the calculated terminal velocity is about 10% smaller than in the reference case, while with the fine one it is 5% greater than in the reference case. With regard to the computation time, grid width h = 1/10 mm is used as reference for the other cases. In Fig. 2 the computed terminal velocities of bubbles with the initial radii 0.5 mm, 0.8 mm, 1.5
6
simulation:
Figure 2: Calculated terminal velocities in comparison to the data from experiments
mm, 2.0 mm and 3.0 mm are plotted in comparison to the data from the literature [4]. The results are all located in the area of the data which are obtained from the experiments. Fig. 3 presents the bubble with the initial radius R = 3.0 mm at different times. The strong deformation and oscillation of the bubble can be observed. In this simulation the bubble does not reach the final form because of the short simulation time and the small viscosity of the fluids. Note in this simulation the measure described in section 2.4 is not used, because the calculation of gradients with the 5th order WENO scheme is accurate enough for such a great bubble and leads only to negligible mass loss during solution of equation (3).
3.2 Detachment of a Bubble from a Wall In this case a bubble rests on the lower wall of the computational domain at the beginning. The properties of the fluids are the same as in Tab. 1. The boundary condition of the level set function at the lower is set to static contact angle. The way to deal with the contact angle at the wall is inherited from the implementation of the VOF method in OpenFOAM [1]. In this approach contact angle is considered as boundary condition, which indicates the direction of the surface tension at the wall. The direction of the surface tension is adjusted according to the contact angle after every calculation of the unit normal of the interface. Fig. 4 shows the position and the form of a bubble with an initial radius 1.5 mm and its
7
t = 0.0 s
t = 0.01s
t = 0.02s
t = 0.03s
t = 0.04s
t = 0.05 s
t = 0.06s
t = 0.07s
t = 0.08s
t = 0.09s
Figure 3: The rising air bubble with initial radius R = 3 mm (2D axisymmetric)
centre 1 mm above the wall at different times. In this simulation the static contact angle of the bubble at the lower wall is set to 60◦ , measured in the liquid phase. As mentioned in section 2.4, the measure to conserve mass could not be utilised in this simulation, otherwise the bubble does not detach from the wall.
8
t = 0s
t = 0.0025s
t = 0.0100s
t = 0.0150s
t = 0.0370s
t = 0.0371s
t = 0.0375s
interface, r = 1.5mm
axis
1mm air
liquid
t = 0.0300s
Figure 4: Bubble detaching from a wall at different times (2D axisymmetric)
4 Conclusions The major procedures and aspects of the implementation of the level set method into OpenFOAM are presented. The terminal velocities calculated with this implementation agree well with the experiments. The capability of the level set method to handle the topological change of the free interface is also presented with the simulation of the bubble detachment from a wall. Solving the equation for advection of the level set function is quite straightforward, while more effort is needed for the reinitialisation of the level set function and the conservation of the mass. These two points are also interesting research fields in general. The method for mass conservation shows its limitation in case of bubble detachment because of its ill distribution of masses. A method to conserve the mass locally can be the first extension of this work.
References [1] Documentation of OpenFOAM. www.openfoam.com. [2] J.U. Brackbill, D.B. Kothe, and C. Zemach. A continuum method for modeling surface tension. Journal of computational physics, 100:335–354, 1992.
9
[3] Y.C. Chang, Hou T.Y., B. Merriman, and S. Osher. A level set formulation of Eulerian interface capturing methods for incompressible fluid flows. Journal of computational physics, 124:449–464, 1996. [4] R. Clift, J.R. Grace, and M.E. Weber. Bubbles, Drops and Particles. Academic Press, Inc., London, 1978. [5] R.I. Issa. Solution of the implicitly discretised fluid flow equations by operatorsplitting. Journal of computational Physics, 62:40–65, 1986. [6] H. Jasak. Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows. PhD thesis, Department of Mechanical Engineering Imperial College of Science, Technology and Medicine University of London, June 1996. [7] S. Majumder and S. Chakraborty. New physically based approach of mass conservation correction in level set formulation for incompressible two-phase flows. Journal of Fluids Engineering, 127:554–563, 2007. [8] S. Osher and R. Fedkiw. Level Set Methods and Dynamics Implicit Surfaces. SpringerVerlag, New York Inc., 2000. [9] S. Osher and J.A. Sethian. Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 79:12–49, 1988. [10] H. Rusche. Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions. PhD thesis, Imperial College of Science, Technology and Medicine Department of Mechanical Engineering University of London, Exhibition Road, London SW7 2BX, December 2002. [11] G. Russo and P. Smereka. A remark on computing distance functions. Journal of computational Physics, 163:51–67, 2000. [12] C.W. Shu. Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws. NASA/CR-97-206253, ICASE Report No. 97-65, Nov. 1997. [13] M. Sussman and E.G. Puckett. A coupled level set and volume-of-fluid method for computing 3d and axisymmetric incompressible two-phase flows. Journal of Computational Physics, 162:301–337, 2000. [14] M. Sussman, P. Smereka, and S. Osher. A level set approach for computing solutions to incompressible two-phase flow. Journal of computational physics, 114:146–159, 1994. [15] S.W.J. Welch and J. Wilson. A volume of fluid based method for fluid flows with phase change. Journal of Computational Physics, 160:662–682, 2000. [16] D.L. Youngs. Time-dependent multi-material flow with large fluid distortion. In K.W. Morton and M.J. Baines, editors, Numerical Methods for Fluid Dynamics, volume 79, pages 273–285. Academic Press, 1982.
10
