Level Set Method
The level set method is a powerful tool for capturing the interface between two immiscible fluids. It is based on the idea of a signed distance function, which is a function that assigns a signed distance to each point in the domain. The level set function is defined as follows:
where \(\phi\) is the level set function, \(\mathbf{x}\) is the position vector, and \(t\) is the time.
Material Properties
The density and viscosity are defined as functions of the level set field:
where \(H(\phi)\) is the Heaviside function:
Time Discretization
For a single level, the momentum equation (4) is advanced by a fractional step method with the approximate projection [11, 12] to enforce the incompressibility condition (equation (3)). The LS advection equation (5) is updated using the Godunov scheme.
At the beginning of each time advancement of level \(l\), the velocity \(\boldsymbol{u}^{n,l}\) and the LS function \(\phi^{n,l}\) at time \(t^{n,l}\) are given. Owing to the fractional step method we used, the pressure is staggered in time and thus the pressure \(p^{n-1/2,l}\) at time \(t^{n-1/2,l}\) is known. Because the single-level advancement concerns only level \(l\), we omit the superscript \({}^l\) in this section. To obtain the updated velocity \(\boldsymbol{u}^{n+1}\), pressure \(p^{n+1/2}\), and LS function \(\phi^{n+1}\) on level \(l\), the solver performs the following steps:
Advance the LS function as
(8)\[\phi^{n+1} = \phi^{n} - \Delta t \left[\nabla \cdot (\boldsymbol{u}\phi) \right]^{n+1/2}\]The advection term in the above equation is calculated using the Godunov scheme.
Solve the intermediate velocity \(\boldsymbol{u}^{*}\) semi-implicitly
(9)\[\begin{split}\begin{aligned} &\boldsymbol{u^{*,n+1}}-\frac{\Delta t}{2\rho(\phi^{n+1/2})Re}{\nabla} \cdot {\mu(\phi^{n+1})}{\nabla}\boldsymbol{u^{*,n+1}} = \boldsymbol{u^n}-\Delta t \left[\nabla \cdot (\boldsymbol{uu})\right]^{n+1/2} + \\ &\frac{\Delta t}{\rho(\phi^{n+1/2})}\bigg[-\nabla p^{n-1/2}+ \frac{1}{2Re}{\nabla}\cdot{\mu(\phi^{n})}{\nabla}\boldsymbol{u^n} + \rho(\phi^{n+1/2}) \frac{z}{Fr^2} - \frac{1}{We}\kappa(\phi^{n+1/2})\delta(x^{n+1/2})\boldsymbol{n}\bigg]. \end{aligned}\end{split}\]In equation (9), the detailed discretization of the advection term \(\nabla \cdot (\boldsymbol{uu})\), viscous term \({\nabla}\cdot({\mu(\phi)}{\nabla}\boldsymbol{u})\), and surface tension term \(\kappa(\phi)\delta(x)\boldsymbol{n}/We\) are implemented using standard finite difference schemes. The LS function at \(t^{n+1/2}\) is calculated by
(10)\[\phi^{n+1/2} = \frac{1}{2}(\phi^{n}+\phi^{n+1})\]where \(\phi^{n+1}\) is obtained from step 1 (equation (8)). The \(\rho(\phi^{n+1/2})\), \(\mu(\phi^{n})\), and \(\mu(\phi^{n+1})\) are then obtained from equations (6) and (7).
Apply the projection method to obtain the pressure and a solenoidal velocity field. To conduct the level projection, a temporary variable \(\boldsymbol{V}\) is defined as
(11)\[\boldsymbol{V} = \frac{\boldsymbol{{u^{*,n+1}}}}{\Delta t} + \frac{1}{\rho(\phi^{n+1/2})} \nabla p^{n-1/2}\]Then the updated pressure \(p^{n+1/2}\) is calculated by
(12)\[L^{cc,\mathrm{level}}_{\rho^{n+1/2}} p^{n+1/2} = \nabla \cdot \boldsymbol{V}\]where \(L^{cc,\mathrm{level}}_{\rho^{n+1/2}}p^{n+1/2}\) is a density-weighted approximation to \(\nabla \cdot (1/\rho^{n+1/2} \nabla p^{n+1/2})\). Finally, the velocity can be calculated as
(13)\[\boldsymbol{{u^{n+1}}} = \Delta t \left(\boldsymbol{V} - \frac{1}{\rho^{n+1/2}} \nabla p^{n+1/2}\right)\]As defined in the AMReX framework, \(\nabla \cdot\) and \(\nabla\) are the cell-centered level divergence operator \(D^{cc,\mathrm{level}}\) and level gradient operator \(G^{cc,\mathrm{level}}\), respectively. The level gradient operator \(G^{cc,\mathrm{level}}\) is not the minus transpose of the level divergence operator \(D^{cc,\mathrm{level}}\), i.e., \(G^{cc,\mathrm{level}} \neq -(D^{cc,\mathrm{level}})^T\). As a result, the idempotency of the approximate projection \(\boldsymbol{P} = I - G^{cc,\mathrm{level}}(L^{cc,\mathrm{level}})^{-1}D^{cc,\mathrm{level}}\) is not ensured, i.e., \(\boldsymbol{P}^{2} \neq \boldsymbol{P}\). Yet, this nonidempotent approximate projection is stable and appears to be well-behaved in various numerical tests and practical applications [13]. Notably, for a uniform single grid with periodic boundary conditions, Lai theoretically proved that this approximate projection method is stable, in that \(\|\boldsymbol{P}\| \leq 1\). It should be noted that the approximate projection is applied to the intermediate velocity \(\boldsymbol{{u^{*,n+1}}}\) (equation (11)). Compared with the form that projects the increment velocity \(\boldsymbol{u^{*,n+1}}-\boldsymbol{u^n}\), e.g. as that used in other methods, the projection method used here can reduce the accumulation of pressure errors and lead to a more stable algorithm. The effectiveness and stability of this approximate projection has been validated through various numerical tests.
Reinitialize the LS function \(\phi\) to maintain \(\phi\) as a signed distance function of the interface and guarantee the conservation of the mass of the two phases. In this step, a temporary LS function \(d(\boldsymbol{x},\tau)\) is updated iteratively using the following pseudo evolution equation:
(14)\[\frac{\partial d}{\partial \tau}=S(\phi)(1-|\nabla d|)\]with the initial condition
(15)\[d(\boldsymbol{x},\tau = 0)=\phi^{n+1}(\boldsymbol{x})\]where
(16)\[S(\phi) = 2\left(H(\phi)-1 / 2\right)\]Here, \(\tau\) is the pseudo time for iterations. A second-order essentially non-oscillatory (ENO) scheme is used to discretize the distance function and a second-order Runge–Kutta (RK) method is applied for the pseudo time advancing. To ensure the mass conservation, \(d(\boldsymbol{x},\tau)\) is further corrected by minimizing the differences of the volume of each fluid between \(\tau=0\) and the final iteration. Finally, the LS function \(\phi\) is re-initialized by the volume corrected \(d\) [10, 14, 15].
At last, we give a summary of the single-level advancement algorithm as follows [10].