Robust control volume finite element methods for numerical wave tanks using extreme adaptive anisotropic meshes

Multiphase inertia‐dominated flow simulations, and free surface flow models in particular, continue to this day to present many challenges in terms of accuracy and computational cost to industry and research communities. Numerical wave tanks and their use for studying wave‐structure interactions are a good example. Finite element method (FEM) with anisotropic meshes combined with dynamic mesh algorithms has already shown the potential to significantly reduce the number of elements and simulation time with no accuracy loss. However, mesh anisotropy can lead to mesh quality‐related instabilities. This article presents a very robust FEM approach based on a control volume discretization of the pressure field for inertia dominated flows, which can overcome the typically encountered mesh quality limitations associated with extremely anisotropic elements. Highly compressive methods for the water‐air interface are used here. The combination of these methods is validated with multiphase free surface flow benchmark cases, showing very good agreement with experiments even for extremely anisotropic meshes, reducing by up to two orders of magnitude the required number of elements to obtain accurate solutions.

options, 1,2 making them very attractive and versatile for very different problems. The classical approach for modeling multiphase flow problems with FEM (CVFEM) is using finite element (FE) discretization for the pressure (continuous Galerkin, CG) and velocity (discontinuous Galerkin, DG) fields, typically with polynomials of order 0 or 1 and with control volumes (CVs) for the phase field. [3][4][5] A common way to reduce the computational cost is the optimization of the mesh, reducing the number of elements in regions of the domain where resolution is not needed. Using anisotropic meshes (high ratio side-to-side elements) can help by reducing the number of nodes even further. Dynamic mesh optimization (DMO) schemes can also be used in order to significantly reduce computational cost without losing accuracy by adapting the mesh at different time steps. [6][7][8][9][10][11][12] This is particularly attractive for multiphase problems where the air-water interface can be resolved with a moving high-resolution anisotropic mesh, such NWTs.
Using highly optimized anisotropic meshes can often lead to elements with large angles, which can make the model take significantly longer (creating matrices that are difficult to solve) or even fail to converge. 13,14 Moreover, DMO could be deemed impractical when highly anisotropic meshes are required if these are seen to be associated with model instability. These limitations have restricted the use of anisotropic meshes and DMO, making good quality isotropic (similar side-to-side element) meshes preferable by most due to their advantages in terms of model stability (creating matrices that are easy to solve), which leads to more elements and a higher computational cost.
To tackle elements with large angles for Darcy flow applications, Salinas et al 15 proposed a simple and robust CV discretization for the pressure field specifically to overcome mesh quality limitations. This work was based on the classic spatial discretizations P 0 DGP 1 and P 1 DGP 2 , that is, using polynomials of order 0 and 1 for the velocity field (with DG) and polynomials of order 1 and 2 for the pressure field (with CG). Although this approach is very attractive for its simplicity and robustness against challenging meshes, the authors believe that the spatial discretization with polynomials of order 0 for the velocity leads to lower accuracy for resolving the momentum equation, which is critical for inertia dominated flow problems. Alternatively, the P 1 DGP 2 discretization gives very high accuracy, but this accuracy comes at a much higher cost due to discretizing the pressure field with second-order polynomials.
This article proposes using a novel P 1 DGP 1 (CV), discretizing with CVs, the pressure field, and with DG of order 1, the critical velocity field. This has the potential to add stability and robustness to the classical P 1 DGP 1 when using highly anisotropic meshes and DMO for computational cost optimization.
Due to the simplicity of the P 1 DGP 1 (CV) approach, other schemes can be added to the model for improving the accuracy of the solution. Different schemes are considered to better represent the physics of the free surface of the water in NWT (see Section 2.4.1), such as a compressive advection scheme 16 to simulate very sharp water-air interfaces. Two benchmark cases are used for validation, focusing on NWTs and their main challenges of wave breaking and wave loading. A collapsing water column is modeled, following experiments carried out by Kleefsman et al 17 for validation of water impact loading, and a solitary wave breaking on a slope is compared to experiments by Li,18 where robust wave breaking simulations with very challenging anisotropic meshes are validated with the novel CV discretization for pressure.
The remainder of this article is structured as follows. Section 2 describes the governing equations of the model, their discretization, and how they are solved. Section 3 shows the results of the validation cases and compares them with experiments. Finally, this article concludes with a discussion chapter and remarks on the scope for future adoption of such methods for NWTs.

Governing equations
In this two phase modeling, a multicomponent approach 19,20 is used. The component advection equations are formulated within the pressure and continuity equations, which results in local mass balance. A number of components ( C ) exist in one or more phases (one is assumed here), where i is the volume fraction of component i, i = 1, 2, … ,  C , and the system is constrained to: For each component i, the equation of conservation of mass is defined as: where t is time, i is the component density, and u is the velocity vector. The equation of motion of an incompressible viscous fluid is: where p is the pressure field, = ∑  c i=1 i i is the bulk dynamic viscosity, = ∑  c i=1 i i is the bulk density, and g is the gravity acceleration vector.

Temporal discretization
The adaptive time step and discretization schemes described by Pavlidis et al 16 are used here. Time step size is adapted based on the Courant-Friedrichs-Levy (CFL) number. An adaptive discretization scheme is used for solving the component volume fraction field based on a parameter , which varies between 0 and 1 following equation (42) from Pavlidis et al. 16 For = 0, an explicit forward Euler scheme is used, which results in a compressive interface, whereas for = 1, a numerically robust implicit backward Euler scheme is used. The advantage of an adaptive discretization scheme is that it is using forward Euler whenever possible (ie, the time step is below the order of the grid Courant number) and implicit Euler (or a combination of the two) where required (ie, the time step is of the order or above the grid Courant number), leading to a combined very robust method and a more compressive interface. The velocity field is solved with an implicit backward Euler scheme in order to guarantee stability.

Computational grids and spatial discretization
This article proposes a novel P 1 DGP 1 (CV) based on the classic P 1 DGP 1 , this is linear discontinuous velocity and CV pressure spatial discretizations. When compared to classical P 1 DGP 1 , the P 1 DGP 1 (CV) approach is more robust when applied to very anisotropic elements, thanks to discretizing the pressure field with CV, and faster because of its higher convergence speed see Salinas et al 15 . In order to keep higher accuracy for the pressure field, a balanced pressure decomposition method is used, 21 which solves the hydrostatic pressure separately with a quadratic element discretization scheme. This method is found to lead to very significant improvements for the P 1 DGP 1 FE pair, at a negligible increase in computational cost, see Maddison et al. 21 The pressure decomposition applies only to the momentum equation and is omitted in the spatial discretization description section of this article for simplicity. Assuming FE and CV representations for u and p, respectively, for the discretization of the governing equations, these can be expressed in terms of their FE basis functions Q k and CV basis function P j : being  u and  p the velocity and pressure number of degrees of freedom, respectively. Q k are the Lagrangian shape functions of order n (here we use n equal to 1) and P j the CV shape functions of order 1 (ie, they take the form of (1 0 0) (0 1 0) (0 0 1) for 2D cases). Discretization of the mass conservation equation (Equation (2)) is well established and can be found in Pavlidis et al. 20 The discretization of the momentum equation (Equation (3)) is, however, different due to the pressure term (resolved here using CV) and is explained here. Multiplying the equation (Equation (3)) with a linear DG basis function Q k , applying integration by parts twice over each element E, with the time discretization method, and Green's theorem for the pressure term, the discrete form of the momentum equation can be formulated as: where V E is the domain, Γ CV j is the surface of the CV, and Γ Ω p is the domain surface. S k ADV and S k VIS are the advection and viscous terms, respectively, described by Xie et al. 22 p CV is the pressure field resolved with CV and p bc the pressure value at the boundary. Note that in this case, because of P j , we cannot explicitly evaluate the gradient of the pressure field as in the standard FE approach, so we therefore need to apply the Green theorem. We have seen that this combination improves the quality of the pressure matrix, helping the convergence and solution of the pressure field.

The transport equation
The transport or flux of each component's volume fraction is modeled with the transport equation: where i is the transported constituent, in this case i = i i , and is the diffusion term. Here a Petrov-Galerkin method 23 is used to obtain as an anti-diffusion term, as described by Pavlidis et al. 16 This method is equivalent to applying a negative diffusion on the interface in order to compensate for the model induced artificial diffusion and make the interface F I G U R E 2 Flow chart of the pressure projection method used here for solving the system of equations more compressive. The flux or transport limiter, which guarantees mass continuity and the method to be conservative and here also controls the compressiveness of the interface, is based on the normalized variable diagram (NVD) and given by Pavlidis et al. 16

Solving the system of equations
A pressure projection method 22 is used to solve the system of equations. GMRES Krylov subspace solver 24 from the PETSc framework 25 with the Boomer algebraic multigrid preconditioner from the HYPRE library 26 is used to solve for pressure, velocity, and transport. The full algorithm for solving the system of equations is shown in Figure 2.

DMO and mesh to mesh interpolation
The mesh adaptivity scheme used here is based on the anisotropic mesh optimization library of Pain et al. 7 A metric advection method 27 is used here in order to take into account the expected movement and direction of the high resolution area between mesh adapts. Mesh to mesh interpolation is used once the new mesh has been generated. 28 A standard consistent interpolation is used for the velocity field. For pressure and component fraction fields, a CG projection 29,30 is used because of its high order interpolation and mass conservation, which also leads to better interface capture as explained in Section 2.4.1.

Interface capture schemes
Running simulations with first-order spatial discretization schemes, and CG pressure field in particular, can lead to a more dissipative interface and therefore less accurate result. Different options are proposed here for keeping a sharp interface without increasing the computational cost: 1. A more compressive interface can be obtained with the adaptive time discretization scheme described in Section 2.2. 2. A CG projection 29,30 is used for the component volume fraction (see Section 2.4). This guarantees a conservative, nondissipative, and nonbounded projection of the water volume fraction field. 3. A high-order accurate Petrov-Galerkin-based compressive advection scheme 16 with a gradient-based scaling is applied.
This method described in Section 2.3.1 is equivalent to adding a negative diffusion into the advection equation. 4. A high-order FEM down-winding scheme is used to obtain fluxes on the CV boundaries, which are subjected to a flux-limiting NVD approach, [31][32][33] to obtain compressive and bounded solutions.

VALIDATION OF THE MODEL
The aim of this validation is to show the potential of the P 1 DGP 1 (CV) spatial discretization scheme for the optimization of accuracy and computational cost for FEM simulations. The two cases chosen for this article are believed to show such potential for NWTs. The results shown here were obtained using the open source Imperial College Finite Element Reservoir Simulator (IC-FERST) with the P 1 DGP 1 (CV) scheme.

2D breaking of a solitary wave
A solitary wave breaking on a slope is simulated in order to prove that the presented approach is capable of modeling wave breaking accurately. This case is run in 2D for sensitivity tests and full validation of the P 1 DGP 1 (CV) and interface capture schemes. The benchmark case used here is based on the experiments run by Li. 18 The dimensions and slope of the NWT are shown in Figure 3. Water depth and wave height at the deepest point are h=0.3048m and H=0.45*h, respectively. Half of the solitary wave is added to the domain by imposing the free surface and horizontal and vertical velocities as initial conditions. The sloped (nonslip), right, and side boundaries are closed (imposing velocities of 0). The top boundary is considered open with an imposed constant pressure of 0. Initial and inlet boundary conditions (left boundary) follow the analytical solution of a solitary wave. 34 Density and viscosity of the water and air components are detailed in Table 2. Mesh is adapted at t=0 with minimum and maximum element sizes in (x, y) of (h, h) and (h∕100, h∕100). Wave breaking, and this case in particular, is very sensitive to numerical wave energy dissipation. Different spatial discretizations are considered in order to evaluate accuracy of the results against experiments. 18 Figure 4 shows the wave profile at different breaking stages for simulations using a classical P 1 DGP 1 and P 1 DGP 1 (CV) discretizations with no compressive schemes for the interface, a P 1 DGP 1 (CV) with compressive interface, and the latter using an LES model. P 1 DGP 1 and P 1 DGP 1 (CV) show lower wave height and late breaking when compared to experiments, but no significant numerical energy dissipation is observed between them. This suggests that CV discretization for pressure does not add significant extra dissipation to the P 1 DGP 1 scheme for this case. The numerical energy dissipation can be significantly reduced by using compressive schemes (described in section 2.4.1). As shown in Figure 4 (green cross), the compressive interface scheme increases wave height and presents a similar wave breaking pattern to experiments, although it also adds an upward curving head to the plunging jet of the wave before impacting the forward free surface, which is not observed in experiments. This effect can be reduced by using a turbulence model such as LES, as shown in Figure 4 (blue x), which suggests that this effect may have been artificially added by the interface compressive scheme (see Section 5 for further discussion). Although the LES model is expected to add extra local dissipation on the breaking wave, this has no significant impact in overall energy dissipation (due to lower or nonturbulent flow before breaking) when comparing simulations with and without LES, as observed in Figure 4.

3D collapsing water column or dam break
The case of a collapsing water column is commonly used as a benchmark case for numerical model validation. [35][36][37][38] The experiments of Kleefsman et al 17 are believed to be the most relevant to compare with numerical simulation due to the presence of a rectangular obstacle in the collapsing water column path and several available pressure measurements from sensors on the obstacle's front and top faces. Simulation results shown here are based on the P 1 DGP 1 (CV) with the interface capture schemes presented in Section 2.4.1. The experiment is based on a classic dam break test (see Figure 5). Pressure sensors are located at the front and top faces of the obstacle, as shown in Figure 5B. Water level gauges are located on the water column, at 0.72 m from the tank wall, and 1 m from the front of the water column.  The top boundary of the numerical tank is considered as an open roof (P = 0) and side and bottom boundaries are closed. Viscosities and densities of air and water are shown in Table 2.
The numerical mesh is initially adapted by the mesh adaptivity algorithm for water component volume fraction with maximum and minimum element sizes of 1 m and 0.01 m, respectively. The aim is to exploit a highly anisotropic mesh with a maximum aspect ratio of 100. Figure 6 shows the adapted mesh at time t = 0 in preparation for highly resolved interface tracking with the previously described compressive interface scheme. Table 3 shows the mesh quality in terms of element angles at t = 0, where up to 8% of elements present a very high angle (usually associated with lower quality mesh but with no negative impacts in this approach). Simulations are run with adaptive time step for 0.1 and 1 CFL numbers for the velocity field in order to provide a full range of time step sizes and their impact on results. The simulation is run for 5 seconds from the dam break. Figure 7 shows the snapshots with a very good representation of water drops in the air after impact (left) and air trapping behind the obstacle (right). These results show the potential of resolving complex two-phase flow with droplets and trapped air bubbles, which can be significantly improved by adding surface tension schemes. 22,39 A complete sequence of snapshots of experiments and numerical model results for CFL of 0.1 are presented in Figure 8. Pressures measured at the obstacle front and top faces are shown in Figure 9. Pressure peaks are zoomed in Figure 10. Water level measurements from experiments and numerical simulations are shown in Figure 11. Experiments and COM-Flow results, which is a numerical method based on finite volumes method (FVM) presented in the original paper by Kleefsman et al,17 are also shown for comparison. Peak pressures at P1 (lower front face) show very good agreement of P 1 DGP 1 (CV) with experiments for both CFL cases. The peak pressure from the model simulation with the presented formulation is expected to be higher (25 kPa in this case) than in the experiment. This is due to pressure sensor missing the real pressure maximum value in experiments by measuring at lower frequency than the simulation time stepping, as shown in Figure 10.
P3 (higher up the front face) results show very good agreement for CFL=0.1 and a slightly smaller peak pressure for CFL=1, compared to the experiments. The authors believe this is induced by a higher energy dissipation due to the longer time step. In both cases, the presented method shows better agreement than COMFlow 17 simulation results. Pressures measured at the top face (P5 and P7) of the obstacle show very good agreement for all cases using P 1 DGP 1 (CV), with the different pressure spikes assumed to be related to water drops falling and turbulent flow on the surface of the obstacle before this is completely submerged. Good agreement with experiments is also found on water level measurements ( Figure 11). Small spikes are observed in pressure plots (for CFL=0.1), but these are very low level fluctuations when compared to other good models such as COMFlow.
Numerical energy dissipation is expected for P 1 DGP 1 (CV) due to solving the pressure field with CV (low order), which should become more evident at the later stages of the simulation. This is observed in Figure 11 as a delay on H2 (A) water level peak at 1.8 second, and as a smoothing effect on peaks at 4.8 second (H2), and 2.8 second and 3.8 second (H4). However, this does not affect the general good agreement with experiments. Due to DMO, the number of elements varies during the simulation (see Figure 12). For CFL=1, a maximum of 2× 10 5 elements is achieved at 1.3 second. For CFL=0.1, a maximum threshold is set to 5× 10 5 . This threshold is reached at 0.8 second and therefore does not affect simulated pressure peaks at sensors P1 and P3. Mesh resolution requirements, that is, maximum (1 m) and minimum (0.01 m) edge element lengths, are the same for all three cases. Time step remains between 10 −3 second and 5×10 −3 second for CFL=1 and between 5×10 −4 second and 3×10 −3 second for CFL=0.1. Overall, the global behavior of the P 1 DGP 1 (CV) scheme results using CFL=0.1 and CFL=1 agree well with the experiment. Impact peaks are in good agreement, especially for P3 with CFL=0.1, which was missed by COMFlow 17 simulation.

3D breaking of a solitary wave
We now extend the 2D wave (see Section 3.1) to 3D in order to fully validate the P 1 DGP 1 (CV) scheme and show the potential of the presented numerical schemes in terms of mesh optimization. This is a particular case where one of the dimensions is initially of less importance compared with the other two dimensions. The goal is, therefore, to reduce the resolution in the y axis direction without compromising the results. For this case, a solitary wave with H∕h = 0.4 (h=0.3048 m) is simulated and compared to experiments. 18 The 3D domain consists of a layout as described in Figure 3 and a flume width of 0.5 m. The mesh is adapted at t = 0 with maximum element length in (x, y, z) of (h, h, h). Three different resolution settings for minimum element lengths are considered (see Table 4): Mesh A for a standard isotropic mesh; Mesh B for a highly anisotropic mesh; and Mesh C for an even higher anisotropy ratio mesh. Figure 13 shows the 3D mesh refinement around the interface for meshes A, B, and C. These three meshes have a starting number of elements of 472342 (A), 30915 (B), and 8742 (C), respectively. statistics for the three mesh settings considered here at t = 0. P 1 DGP 1 (CV) scheme results from mesh setups (B) and (C) can be seen in Figure 14, which shows a very close match up to the fourth breaking stage just before the plunging jet impacts the water surface. This indicates that the latest stages of the wave breaking process are more sensitive to mesh changes. However, the authors believe that, as seen in 2D simulations, this could be due to the high sensitivity from the compressive interface schemes used here.
It is important to note that no comparison between P 1 DGP 1 (CV) and the classical P 1 DGP 1 schemes could be shown here as the classical P 1 DGP 1 scheme failed to converge with mesh B and mesh C. A qualitative comparison with experiments of the wave breaking (plunging) is shown in Figure 15 with the mesh C setup. The mesh requirements of this simulation using P 1 DGP 1 (CV) reduces 2 and 1 orders of magnitude the number of nodes needed for an accurate result when compared to mesh A and mesh B setups (see Table 6), while it is in very good agreement with experiments. 18

DISCUSSION
The novel P 1 DGP 1 (CV) scheme is proven to be as accurate as P 1 DGP 1 standard CVFEM scheme, while being more robust and capable of resolving free surface flow problems with highly anisotropic meshes accurately. These code developments when applied to NWTs effectively accelerate FEM approaches with their typical higher accuracy, when compared to generally established faster NWTs based on FVM or SPH alternatives. [40][41][42] The FEM speed up can be attributed to a new level of mesh optimization with a significant reduction in the number of elements and computational cost, making P 1 DGP 1 (CV) more competitive and attractive. This model acceleration is also attributed to a faster convergence of the solution and to the optimization of the mesh. The validation cases used here have shown that the formulation is capable of providing an accurate solution with a reduction of up to 1 and 2 orders of magnitude compared with standard isotropic meshes or meshes with low anisotropy rates. Compressive interface schemes are really useful for tackling numerical energy dissipation. The interface capture scheme used here is proven to be conservative in terms of energy dissipation, but it generates artificial lower pressure near the crest of the wave when the wave breaking jet is created. Balancing uplift velocities and an artificially induced shape not seen in experiments are generated. This side effect of the compressive interface scheme can be smoothed with an LES model. Future work will aim to improve the performance of the interface capture scheme.
Traditionally more expensive unstructured grids can be combined with DMO helping optimize the mesh and reduce the computational cost of the simulations. However, the standard formulations would fail to solve the system with very obtuse elements. The presented formulation allows to push the boundaries of DMO by allowing to increase the aspect-ratio of elements and further reduce the computational costs while maintaining accuracy. Therefore, this is not only providing a new scheme to overcome a long standing problem for classical FEM models but also has a great potential for a wide range of applications. This is of particular interest for 2D and 3D cases where one or more dimensions have less importance in solving the physics of the problem during part of the simulation or parts of the domain, such as with NWT applications. Future work will look at NWTs with WSI, mainly working with rubblemound structures, where high accuracy and adaptive highly optimized meshes will have a significant advantage and will help contributing to the understanding the physics of the problem.

CONCLUDING REMARKS
The P 1 DGP 1 (CV) scheme is proven to be more robust than the P 1 DGP 1 scheme for bad quality meshes. The P 1 DGP 1 (CV) scheme is also proven capable of resolving free surface flow problems with highly anisotropic meshes accurately, facilitating a new level of mesh optimization with a significant reduction in the number of elements and computational cost. Energy dissipation under the P 1 DGP 1 (CV) scheme is of the same order of magnitude as for the standard P 1 DGP 1 scheme. The numerical energy dissipation in the validation cases is mostly attributed to the interface capture, which can be significantly reduced by applying compressive schemes for the interface. Compressive scheme side effects on the breaking jet can be smoothed with an LES model.
Overall, the novel and very robust P 1 DGP 1 (CV) scheme is able to provide very accurate results with a significant reduction of the number of elements. A speed up associated with reduction in element number of up to 2 orders of magnitude is demonstrated for the benchmark cases simulated for this article.