VOF-PLIC test case of sheared bubble (2D)
This test case solves the advection of a disc of fluid in a sheared analytical velocity field. It is particularly adapted to measure mass-conservation. The objective is to compare two VOF-PLIC (Piecewise Linear Interface Construction) advection methods, the conservative method presented in [1] and the non conservative method presented in [2].
The simulation takes place in a 2D box of coordinate \( (0,0) \) and \( (1,1) \). This domain is discretized uniformly with N cells in each directions. A bubble of fluid A is initialized in a domain full of fluid B. The coordinates of the center of the bubble at initial time is \( (0.5,0.75) \) with a \(0.15m\) radius. The simulation runs during \(T=1.5s\).
Four velocity fields are proposed in this test case, from case A to case D. Velocity fields of cases C and D are respectively the symetric of cases A and B with respect to the line \( x=0.5 \).
CASE | \(u(x,y)\) | \(v(x,y)\) |
---|---|---|
A | \( x \) | \(-y \) |
B | \( x \) | \(1-y \) |
C | \(-1+x \) | \(-y \) |
D | \(-1+x \) | \(1-y \) |
![Figure 1: Velocity field for the CASE A] (velocity_field_bubble_distorsion.png)
To ensure that CFL of the method is always respected, we set a time step \( \Delta t \) such as
\begin{align} \Delta t = \frac{0.2}{N} \end{align}
We adopt the developpement of the exact solution proposed in [2] to our test case A.
\begin{align} u(x,y) = x \\ v(x,y) = -y \end{align}
We also know that :
\begin{align} u(x,y) = \frac{\partial x(t)}{\partial t} \Rightarrow u'-u=0 \qquad (1) \\ v(x,y) = \frac{\partial y(t)}{\partial t} \Rightarrow v'+v=0 \qquad (2) \end{align}
(1) and (2) are ordinary differencial equations with initial conditions defined such as :
\begin{align} x(t=0) = x_0 \\ y(t=0) = y_0\\ u(t=0) = x_0 \\ v(t=0) = -y_0 \end{align}
then solving (1) and (2) we obtain the following paraametric equation system :
\begin{align} x(t) = x_0e^{t} \\ y(t) = y_0e^{-t} \end{align}
which describes hyperbolic network curves. Let consider fluid particles which at initial time take part in the \( (x_{c0},y_{c0}) \) centered circle of radius \( R_0 \), we obtain :
\begin{align} (x_0-x_{c0})^2+(y_0-y_{c0})^2 = R_0^2 \end{align}
Let \( x_c(t) \), \( y_c(t) \) and \( x(t) \), \( y(t) \) the coordinate of the center and fluid particle of the bubble at the initial time. Then we get :
\begin{align} (x(t)-x_{c}(t))^2 = (x_0-x_{c0})^2e^{2t} \\ (y(t)-y_{c}(t))^2 = (y_0-y_{c0})^2e^{-2t} \end{align}
so that leads to the following ellipse cartesian equation of the form :
\begin{align} \left( \frac{x(t)-x_{c0}}{R_0e^t}\right)^2 + \left( \frac{y(t)-y_{c0}}{R_0e^{-t}}\right)^2 = 1 \end{align}
with for semi major axis \( a = R_0e^t \) and semi minor axis \( b = R_0e^{-t} \).
Finally, to initialize an ellipse in Notus we need to define a polygon with a given number of points \(n\). The coordinates of each point \(i\) are given by:
\begin{align} x(i) = acos(i) = R_0e^{0.8}cos(i\pi/n) \\ y(i) = bsin(i) = R_0e^{-0.8}sin(i\pi/n) \end{align}
Next figures show isocontour of the volume fraction value 0.5 for cases A and B at initial time, intermediate and final time. Since isocontour for conservative and non conservative methods are surimposed, only one is shown.
![Figure 2: Isocontour of the volume fraction value 0.5 for the CASE A, from t=0s, t=0.4s to t=0.8s] (CASEB.png)
![Figure 3: Isocontour of the volume fraction value 0.5 for the CASE B, from t=0s, t=0.4s to t=0.8s] (CASEA.png)
In order to evaluate mass loss \( \varepsilon \) during the simulation, we compute the difference between the volume fraction of fluid at initial time and at time T.
\begin{align} \varepsilon = \sum \limits^{N^2}_{i=1}f_i - \sum \limits^{N^2}_{i=1} f_i^0 \end{align}
To evaluate mesh convergence, we use the \( L_1 \) norm of the difference between the initial and the current volume fraction:
\begin{align} L_1 = \sum \limits^{N^3}_{i=1} |f_i - f^0_i| \\ \end{align}
We also compute local convergence order \( p \) of the method evaluating the evolution of the error with the mesh refinement:
\begin{align} p = \frac{ln\left( \frac{E_2}{E_1} \right)}{ln\left (\frac{N_1}{N_2} \right)} \end{align}
where \( E_1 \) and \(E_2 \) are error for mesh 1 and mesh 2, of size \( N_1\) and \(N_2\) respectively.
The following tables show, for different spatial discretizations the material loss for the two diffents advection methods.
\begin{align} \textbf{Non-conservative method} \end{align}
Grid size N | \( \Delta x \) | \( \varepsilon_{air}\) (non conservative method) | \( L_1 \) | p |
---|---|---|---|---|
8 | 1/8 | \( 3.75 \times 10^{-2} \) | \( 3.75 \times 10^{-2}\) | |
16 | 1/16 | \( 6.97 \times 10^{-4} \) | \( 9.36 \times 10^{-3}\) | \(2.00 \) |
32 | 1/32 | \( -3.53 \times 10^{-4} \) | \( 2.39 \times 10^{-3}\) | \(1.97 \) |
64 | 1/64 | \( 1.76 \times 10^{-4} \) | \( 1.05 \times 10^{-3}\) | \(1.18 \) |
128 | 1/128 | \( 8.82 \times 10^{-5} \) | \( 4.69 \times 10^{-4}\) | \(1.16 \) |
256 | 1/256 | \( 4.42 \times 10^{-5} \) | \( 2.30 \times 10^{-4}\) | \(1.02 \) |
512 | 1/512 | \( 2.21 \times 10^{-5} \) | \( 1.08 \times 10^{-4}\) | \(1.09 \) |
\begin{align} \textbf{Conservative method} \end{align}
Grid size N | \( \Delta x \) | \( \varepsilon_{air}\) (conservative method) | \( L_1 \) | p |
---|---|---|---|---|
8 | 1/8 | \( -1.11 \times 10^{-16} \) | \( 3.76 \times 10^{-2}\) | |
16 | 1/16 | \( 5.55 \times 10^{-16} \) | \( 9.25 \times 10^{-3}\) | \( 2.02 \) |
32 | 1/32 | \( -1.11 \times 10^{-16} \) | \( 2.37 \times 10^{-3}\) | \(1.96 \) |
64 | 1/64 | \( 0.0 \) | \( 1.01 \times 10^{-3}\) | \(1.23 \) |
128 | 1/128 | \( 1.11 \times 10^{-16} \) | \( 4.75 \times 10^{-4}\) | \(1.08 \) |
256 | 1/256 | \( 1.11 \times 10^{-16} \) | \( 2.34 \times 10^{-4}\) | \(1.02 \) |
512 | 1/512 | \( 1.11 \times 10^{-16} \) | \( 1.10 \times 10^{-4}\) | \(1.09 \) |
Results confirm that even though the method presented in [2] shows pretty good results concerning mass conservation, the one presented in [1] is fully conservative (up to machine precision). Moreover, the convergence order is near to one for both methods and nearly same errors are noticed.
[1] G.D. Weymouth, Dick K.-P. Yue, Conservative Volume-of-Fluid method for free-surface simulations on Cartesian-grids, Journal of Computational Physics 229 (2010) 2853–2865
[2] Jérôme Breil. Modélisation du remplissage en propergol de moteur a propulsion solide. Mécanique des fluides [physics.class-ph]. Université de Bordeaux 1, 2001. Français. tel-0147869. https://tel.archives-ouvertes.fr/tel-01478691