# Computation of the van Keken Problem with the Volume-of-Fluid Interface Tracking Method#

*This section is a co-production of Jonathan Robey and E. G. Puckett.*

One can also model the van Keken problem with the volume-of-fluid (VOF) interface tracking algorithm in ASPECT. In fact, this problem is particularly well-suited to being computed with the VOF method, since it consists of two distinct, immiscible fluids and interface tracking algorithms are specifically designed not to allow the two fluids to mix. In particular, assuming the computation is sufficiently well-resolved, the fluids will not mix at sub-grid scales over the entire duration of the computation. However, note that this implies that all computations of the van Keken problem made with the VOF method must necessarily be with discontinuous initial conditions. Finally, since one is often interested in a high resolution image of the shape of the interface between the two fluids, one can use the output of the VOF method to examine the interface within individual cells or within regions consisting of groups of cells.

Another advantage the VOF method has over modeling moving interface problems with compositional fields is that for problems in which the interface occupies a relatively small part of the computational domain, all of the computational work in the VOF method is done in cells that lie in a neighborhood of the interface, rather than in the entire computational domain.

As noted above, when the interface is discontinuous, the van Keken problem is a version of the Rayleigh-Taylor problem, which is unstable to perturbations of all wavelengths whether the two fluids have the same viscosity or different viscosities (e.g. see Chandrasekhar [1961]). Therefore, it is extremely sensitive to the initial conditions. In order to address this sensitivity, we do not use the default approach of computing the initial material volume fractions using a composition quadrature. Instead we compute the initial volume fractions using a signed distance function \(\phi\) as follows [Robey, 2019, Robey and Puckett, 2019].

First we create the function \(\phi\), which has the following two properties: 1) it is positive in the region that contains one of the fluids, which we will refer to as fluid 1, and negative in the complement of this region, which we will refer to as fluid 2, and 2) at each point in the domain the magnitude \(| \phi |\) of \(\phi\) is the distance to the boundary between the two fluids or materials. In the computations shown here, we use an approximation \(\tilde{\phi}\) to \(\phi\) such that the difference between \(\tilde{\phi}\) and \(\phi\) is small enough for the purposes of making the computations high-quality. The use of an approximation as opposed to the function itself is due to obtaining an appropriate function for almost all nontrivial (not a line or a circle) boundaries is extremely difficult. In this case, because the boundary is close to horizontal, we use the vertical distance to the boundary based on the argument that the gradient will not differ sufficiently from \(1\) to induce errors in the initialization computation. The primary advantage of choosing this particular initialization algorithm is that it allows us to more accurately reproduce the initial condition on a sub-grid scale than would otherwise be possible on the coarser grid on which we compute the time evolution of the interface.

```
subsection Compositional fields
set Number of fields = 1
set Compositional field methods = volume of fluid
end
subsection Volume of Fluid
set Number initialization samples = 16
end
subsection Initial composition model
set Model name = function
set Volume of fluid initialization type = C_1:level set
subsection Function
set Variable names = x,z
set Function constants = pi=3.1415926
set Function expression = 0.2+0.02*cos(pi*x/0.9142)-z
end
end
```

The relevant sections of the parameter file for this type of initialization of
the VOF method appears immediately above. In particular, the combination of
`Number of initialization samples`

with the `level set`

initialization type
indicates that our initialization will consist of dividing each grid cell into
\(16 \times 16\) subcells and the distance to the given initial interface
\(f(\mathbf x)\), provided in `Function expression`

, is computed in each of the
256 subcells. We then use this information to compute a piecewise linear
interface approximation to \(f(\mathbf x)\). The volume fraction in each subcell
is then found in the manner described in Robey [2019], Robey and Puckett [2019]). This initialization procedure provides a much finer and thus, more
accurate, initial condition than the standard VOF initialization procedure
described above.

While the visualization configuration in a typical parameter file is sufficient for most purposes, when using the VOF method one has the ability to see the division between the fluids reconstructed by the VOF algorithm in each cell. This is accomplished by plotting the zero contour of a field \(\tilde\psi\) that is generated to be \(0\) on the reconstructed interface, positive in the region with fluid 1, and negative in the region with fluid 2. However \(\tilde{\psi}\) does not satisfy the requirement that the magnitude is equal to the distance to the interface as would be required for the signed distance function \(\phi\). The modifications to the parameter file that are necessary in order to draw the reconstructed boundary as a contour are shown immediately below. The full configuration file for this version of the benchmark problem can be found at cookbooks/van-keken-vof/van-keken-vof.prm.

```
subsection Postprocess
set List of postprocessors = visualization, velocity statistics
subsection Visualization
set List of output variables = volume of fluid values
set Output format = vtu
set Time between graphical output = 100
subsection Volume of Fluid
set Output interface reconstruction contour = true
end
end
end
```

We made a number of computations of the van Keken problem with the VOF method in order to compare the wall clock times with computations using a DG compositional field. We ran both on the same cluster at global refinements 5–8 using one node with four CPUs and refinements 9 and 10 using two nodes with 16 CPUs. Our results are shown in Table 7. In all of the computations shown in Table 7 we used a CFL number of \(\sigma=0.5\). Due to the change in the CFL number from \(\sigma = 1.0\) in Table 6 to \(\sigma = 0.5\) in Table 7 and the difference between HPC clusters on which the computational results shown in the two tables were made, we can’t make a direct quantitative comparison between the data in Table 6 and Table 7.

However, we can compare the required run time for a VOF computation to that for a DG computation. We note that the use of the VOF advection algorithm significantly reduces the required computation time in all cases, frequently requiring less than half the time required by the DG compositional field.

We now examine the RMS velocity data shown in Fig. 160. Other than for the case of 5 levels of uniform global refinement, the curves for the RMS velocities for \(6\), \(7\), \(8\), \(9\) and \(10\) levels of refinement in Fig. 160 are nearly indistinguishable.

Upon examining the solution at the final time, we note that the general structure of the solution shown in Fig. 161 matches the form and the general structure found in other versions of this benchmark such as in Fig. 150. We also note that the differences in the shape of the interface based on a single refinement as shown in Fig. 161 are minor, although still slightly visible. This is to be expected as refinement is a perturbation of the initial condition at a smaller wave length.

TODO

Table title should contain references, however, Sphinx does not appear to support references in table titles at the moment. Original caption:
Comparison of runtimes for the van Keken problem with VOF
and a DG compositional field, in which the initial conditions for
DG smoothed are as described in The van Keken thermochemical composition benchmark. The times shown are for the full
computation, ending at \(t_\text{end} = 2000\) with a CFL number of
\(\sigma=0.5\) in both cases.
All of these computations were made with ASPECT version 2.2.0-pre
(master, commit `ef542ecc2`

) in release mode on the Peloton2 cluster at
U.C. Davis. We note that the change in the CFL number \(\sigma\) and the
differing choice of cluster makes a direct quantitative comparison
between this table and Table 6 invalid due to
too many confounding factors.

Global Refinement |
Number of Processors |
VOF |
DG |
---|---|---|---|

5 |
4 |
1.33 minutes |
2.57 minutes |

6 |
4 |
8.51 minutes |
19.5 minutes |

7 |
4 |
1.15 hours |
2.49 hours |

8 |
4 |
8.53 hours |
19.6 hours |

9 |
16 |
16.30 hours |
2.72 days |

10 |
16 |
5.17 days |
>6.00 days |

The consistency of the results shown here differs noticeably from the behavior of the problem with discontinuous initial conditions when computed with the FEM and DG advection algorithms. One possible reason for these differences is the specialized initialization procedure used for the volume of fluid method, which permits a much more consistent initialization by reducing the variation in the initial condition when the initial mesh is refined.

To study this feature of our algorithm and the sensitivity of the problem to the precise initial condition, we vary the size of the initial interface perturbation and examine the sensitivity of the final results to a small change in the initial conditions. Specifically, we vary the amplitude \(a\) of the cosine function in the initial conditions, as shown below.

```
subsection Initial composition model
set Model name = function
set Volume of fluid initialization type = C_1:level set
subsection Function
set Variable names = x,z
set Function constants = pi=3.1415926, a=0.02
set Function expression = 0.2+a*cos(pi*x/0.9142)-z
end
end
```

In these computations we vary the value of \(a\) from its usual value of \(a = 0.02\) to \(5\% = 0.001\) below its usual value to \(5\%\) above its usual value in increments of \(0.01\). In other words, we compare the values for \(a =0.019\), \(0.020\), and \(0.021\). Upon examination of Fig. 162, we see a visible variation in the location of the second peak, although the overall shape of the curve remains consistent with the curves in Fig. 160. The size of this variation in the initial conditions cannot be expected to be reproduced using the standard compositional quadrature initialization procedure for VOF unless the cell size is on the scale of the change in the value of \(a\); i.e., \(h \, \lessapprox \, \Delta a = 0.001\). We also note that the smoothing parameter which would produce a \(10^{-3}\leq C \leq 1 - 10^{-3}\) band on the order of the same size as the amplitude variation shown here, would be approximately \(2.8 \cdot 10^{-4}\). This perturbation is much smaller than any of the changes in width of the smoothed regions in the computations shown in Fig. 156. In summary, these results demonstrate the sensitivity of the discontinuous version of the van Keken problem to even extremely small variations in the initial conditions.