PROGRESSIVE WAVE SIMULATION USING STABILIZED EDGE-BASED FINITE ELEMENT METHODS

Flows involving waves and free-surfaces occur in several problems in hydrodynamics, such as sloshing in tanks, waves breaking in ships and motions of offshore platforms. The computation of such wave problems is challenging since the water/air interface (or free-surface) commonly present merging, fragmentation and cusps, leading to the use of interface capturing Arbitrary Lagrangian-Eulerian (ALE) approaches. In such methods the interface between the two fluids is captured by the use of a marking function which is transported in a flow field. In this work we simulate these problems with a 3D incompressible SUPG/PSPG parallel edgebased finite element flow solver associated to the Volume-ofFluid (VOF) method [1]. The hyperbolic equation for the transport of the marking function is also solved by a fully implicit parallel edge-based SUPG finite element formulation. Global mass conservation is enforced adding or removing mass proportionally to the absolute value of the normal velocity at the interface. The performance and accuracy of the proposed solution method is tested in the simulation of progressive waves and the interaction of a fixed cylinder with a progressive wave. INTRODUCTION Flows involving waves occur in several problems in hydrodynamics. Sloshing of liquids in tanks, wave breaking in ships, offshore platforms motions and green water on decks are important examples of this class of problems. The main computational challenge when solving such highly nonlinear problem is determining the evolution of the interface location. There are a large number of numerical methods devoted to the computation of free-surface problems. These methods are frequently classified as interface tracking and interface capturing methods. Interface tracking methods are based on a Lagrangian framework where the moving interface or boundary is explicitly tracked by the computational grid or by the particles of meshless methods which must be deformed or moved in order to follow the fluid flow. The deforming-spatialdomain/stabilized space time finite element formulation (DSD/SST) proposed by Tezduyar [1] and Tezduyar et al [2, 3] is a mesh based example of interface tracking method. Particle methods, such as those of Koshizuka et al [4] and Violeau and Issa [5] are examples of smoothed particle hydrodynamics (SPH) methods to the simulation of free-surface problems. However, these methods still present a high computational cost since they need to compute the interaction between the particles using search algorithms. As a compromise between the advantages offered by mesh based and meshless methods, Del Pin et al present in [6] the particle finite element method (PFEM) applied to free-surface flows. In this method the critical parts of the continuum are discretized with particles while the remaining parts are treated by a Lagrangian finite element formulation. Another technique mixing Lagrangian and Eulerian flavors was proposed in [7] by Takizawa et al. In this work the authors enhanced the Constrained Interpolation Profile method (CIP) for solving hyperbolic equations with a meshless Soroban grid. The resulting formulation was used to treat fluidobject and fluid-structure interaction in the presence of freesurfaces. As a cost effective alternative to interface tracking methods, interface capturing methods have emerged. Interface capturing methods are Eulerian in their concept, thus they rely on a unique and fixed computational grid to capture the interface evolution. In this class of methods the interface is represented by a scalar function which marks the regions filled with the fluids involved. In other words, the interface position is implicitly captured in a scalar marking function value and the interface evolution is determined by the additional cost of solving an advection equation for the marker. As opposed to interface tracking methods, interface capturing methods require little effort to represent all complicated features of moving interfaces. Additionally, the parallel implementation and postprocessing of interface capturing methods are straightforward. The main drawback of interface capturing methods is the need to average the fluid properties at the interface cells (elements) due to the discontinuity of the Eulerian representation of the interface. Moreover, the accuracy and computational cost of interface capturing methods are typically associated to grid resolution, properties of the marking function chosen to represent the interface and numerical methods for solving the fluid flow and marking function advection. The well known volume-of-fluid scheme (VOF), firstly proposed by Hirt and Nichols [8] for Cartesian grids, is an interface capturing technique which employs a step function ranging from 0 to 1 to represent the fraction of fluid within the grid cells. In this sense, the interface is implicitly represented by the partially filled cells. The main issues associated to VOF methods include the difficulty in advecting a discontinuous step function and the accurate modeling of surface tension effects. Level set methods [9] implement free-surface flows in a different manner than VOF by changing the marking function employed to represent the interface. Therefore, the fluids are associated to the range of the distance function signs while the interface is implicitly represented by the zero level set. However, the level set method suffers when the distance function looses its properties and must be rebuilt. In fact, the success of level set method lies in its ability of building and keeping a signed distance function without losing its properties. The enhanced-discretization interface-capturing method, firstly proposed by Tezduyar in [10] and the work of Lohner et al presented in [11], are both examples of unstructured grid formulations based on the finite element method to solve free-surface flows using a VOF marking function. In this work we use our VOF edge-based stabilized finite element solver [12] and [13], to deal with complex wave phenomena. The main characteristics of our solver are: Streamline-Upwind/Petrov-Galerkin (SUPG), PressureStabilizing/Petrov-Galerkin (PSPG) and Least Squares Incompressibility Constraint (LSIC) stabilized finite element formulation; implicit time marching scheme with adaptive time stepping control; advanced Inexact Newton solvers; edge-based data structures to save memory and improve performance; support to message passing and shared memory parallel programming models; and large eddy simulation extensions using a classical Smagorinsky model.


INTRODUCTION
Flows involving waves occur in several problems in hydrodynamics. Sloshing of liquids in tanks, wave breaking in ships, offshore platforms motions and green water on decks are important examples of this class of problems. The main computational challenge when solving such highly nonlinear problem is determining the evolution of the interface location. There are a large number of numerical methods devoted to the computation of free-surface problems. These methods are frequently classified as interface tracking and interface capturing methods. Interface tracking methods are based on a Lagrangian framework where the moving interface or boundary is explicitly tracked by the computational grid or by the particles of meshless methods which must be deformed or moved in order to follow the fluid flow. The deforming-spatialdomain/stabilized space time finite element formulation (DSD/SST) proposed by Tezduyar [1] and Tezduyar et al [2,3] is a mesh based example of interface tracking method. Particle methods, such as those of Koshizuka et al [4] and Violeau and Issa [5] are examples of smoothed particle hydrodynamics (SPH) methods to the simulation of free-surface problems. However, these methods still present a high computational cost since they need to compute the interaction between the particles using search algorithms. As a compromise between the advantages offered by mesh based and meshless methods, Del Pin et al present in [6] the particle finite element method (PFEM) applied to free-surface flows. In this method the critical parts of the continuum are discretized with particles while the remaining parts are treated by a Lagrangian finite element formulation. Another technique mixing Lagrangian and Eulerian flavors was proposed in [7] by Takizawa et al. In this work the authors enhanced the Constrained Interpolation Profile method (CIP) for solving hyperbolic equations with a meshless Soroban grid. The resulting formulation was used to treat fluidobject and fluid-structure interaction in the presence of freesurfaces.
As a cost effective alternative to interface tracking methods, interface capturing methods have emerged. Interface capturing methods are Eulerian in their concept, thus they rely on a unique and fixed computational grid to capture the interface evolution. In this class of methods the interface is represented by a scalar function which marks the regions filled with the fluids involved. In other words, the interface position is implicitly captured in a scalar marking function value and the interface evolution is determined by the additional cost of solving an advection equation for the marker. As opposed to interface tracking methods, interface capturing methods require little effort to represent all complicated features of moving interfaces. Additionally, the parallel implementation and postprocessing of interface capturing methods are straightforward. The main drawback of interface capturing methods is the need to average the fluid properties at the interface cells (elements) due to the discontinuity of the Eulerian representation of the interface. Moreover, the accuracy and computational cost of interface capturing methods are typically associated to grid resolution, properties of the marking function chosen to represent the interface and numerical methods for solving the fluid flow and marking function advection. The well known volume-of-fluid scheme (VOF), firstly proposed by Hirt and Nichols [8] for Cartesian grids, is an interface capturing technique which employs a step function ranging from 0 to 1 to represent the fraction of fluid within the grid cells. In this sense, the interface is implicitly represented by the partially filled cells. The main issues associated to VOF methods include the difficulty in advecting a discontinuous step function and the accurate modeling of surface tension effects. Level set methods [9] implement free-surface flows in a different manner than VOF by changing the marking function employed to represent the interface. Therefore, the fluids are associated to the range of the distance function signs while the interface is implicitly represented by the zero level set. However, the level set method suffers when the distance function looses its properties and must be rebuilt. In fact, the success of level set method lies in its ability of building and keeping a signed distance function without losing its properties. The enhanced-discretization interface-capturing method, firstly proposed by Tezduyar in [10] and the work of Lohner et al presented in [11], are both examples of unstructured grid formulations based on the finite element method to solve free-surface flows using a VOF marking function.
In this work we use our VOF edge-based stabilized finite element solver [12] and [13], to deal with complex wave phenomena. The main characteristics of our solver are: Streamline-Upwind/Petrov-Galerkin (SUPG), Pressure-Stabilizing/Petrov-Galerkin (PSPG) and Least Squares Incompressibility Constraint (LSIC) stabilized finite element formulation; implicit time marching scheme with adaptive time stepping control; advanced Inexact Newton solvers; edge-based data structures to save memory and improve performance; support to message passing and shared memory parallel programming models; and large eddy simulation extensions using a classical Smagorinsky model.

GOVERNING EQUATIONS
Incompressible Fluid Flow Let sd n Ω ⊂ » be the spatial domain, where n sd is the number of space dimensions. Let Γ denote the boundary of Ω . We consider the following velocity-pressure formulation of the Navier-Stokes equations governing the incompressible flow of two immiscible fluids within an Arbitrary Lagrangian-Eulerian frame [14]: where ρ and u are the density and velocity, mesh u is the mesh velocity, f is the body force vector carrying the gravity acceleration effect and σ σ σ σ is the stress tensor given as ( ) where p is the pressure, I is the identity tensor, T is the deviatoric stress tensor and ( ) u ε ε ε ε is the strain rate tensor defined as In the present work a large eddy simulation (LES) approach to turbulence is considered by the use of a classic Smagorinsky turbulence model [15]. In this model, the viscosity µ is augmented by a subgrid-scale viscosity SGS µ proportional to the norm of the local strain rate tensor and to a filter width h defined here as the cubic root of the element volume, where S C is the Smagorinsky constant, taken as 0.1.
The essential and natural boundary conditions associated with equations (1) and (2) can be imposed at different portions of the boundary Γ and represented by, In the above equation the first four integrals on the left hand side represent terms that appear in the Galerkin formulation of the problem (1)-(8), while the remaining integral expressions represent the additional terms which arise in the stabilized finite element formulation. Note that the stabilization terms are evaluated as the sum of element-wise integral expressions, where n el is the number of elements in the mesh. The first summation corresponds to the SUPG term and the second to the PSPG term. We have evaluated the SUPG and PSPG stabilization parameters according to Tezduyar et al. [2], as follows: here u h is the local velocity vector, and ν is the kinematic viscosity.
In Equation (9), the last summation is the Least Squares Incompressibility Constraint (LSIC) term [3], added to prevent spurious oscillations at high Reynolds number flows. The LSIC stabilization parameter is, The discretization of Equation (9) leads us to a nonlinear system of equations to be solved at each time step.

Interface Capturing
In volume-of-fluid (VOF) method [8], a scalar marking function can be employed to capture the position of the interface between the fluids by simply using the fluids fraction relationship.
The volume-of-fluid can be stated as: assuming the value 1 in regions filled with fluid A, e.g., water, and the value 0 in regions filled with fluid B, e.g., air, the position of the fluid interface will be defined by the isovalue contour Finally, the function ( ) φ x is driven by a velocity field u satisfying the following transport equation, given in conservative form as, In the volume-of-fluid (VOF) formulation the fluid density and viscosity, employed in the fluid flow solution, are interpolated across the interface as follows: where subscripts A and B denote the values corresponding to each fluid.
The finite element formulation of Equation (12) can be written as follows: where h S φ and h V φ are standard test and weight finite element spaces. The first two integrals represent the Galerkin formulation of Equation (12), while the first element-wise summation represents the SUPG and the second summation term is the nonlinear discontinuity-capturing term, useful when sharp gradients and/or boundary layers are present in directions other than the streamlines [16]. The evaluation of SUPG τ and δ stabilization terms follow the definitions described in [17] and [16] respectively. The discretization of Equation (15) leads us to a non-linear ordinary differential equation system, due to the discontinuity-capturing term. In this work, we have adopted the YZ β discontinuity capturing term as proposed in [18] where δ is computed as: where ( ) e R φ is the element residual of Equation (15) Note that if 1 β = and the reference value 1 φ = , the YZ β discontinuity capturing term renders to the Consistent Approximated Upwind (CAU) method [16]. In order to avoid non physical results, values lying outside the range [0,1] were truncated with the following function:

SOLUTION PROCEDURE
The computational solution kernels consist of predictor multicorrector time integration schemes as described in [17], [19] for both incompressible fluid flow and interface transport equations. The generalized trapezoidal rule is employed in the time discretization. The nonlinearities due to the convective term on the Navier-Stokes equation are treated by an Inexact Newton-GMRES algorithm as described in Elias et al [20]. In this solution algorithm, at the beginning of the nonlinear iterations in each time step, the algorithm computes large linear tolerances, producing fast nonlinear steps. As the iterations progress towards the solution, the inexact nonlinear method adapts the GMRES tolerances to reach the desired accuracy. A nodal-block diagonal preconditioner is employed for the flow and a simple diagonal preconditioner for the marker. Moreover, for both the fluid flow and the marking function advection we use an adaptive time stepping procedure based on a Proportional-Integral-Derivative (PID) controller (see [21]) for further details). Most of the computational effort spent in this solution procedure is due to the matrix-vector products within the GMRES driver for both flow and marker. To improve the computational efficiency with respect to standard element-byelement and sparse matrix-vector storage schemes, we adopt an edge-based data structure in order to minimize indirect memory addressing, diminish floating point operation counts (flops) and memory requirements, as described in Elias et al [22] and Coutinho et al [23] for both the Navier-Stokes equations and the marking function advection. Further computational gains are obtained from data preprocessing performed by the EdgePack library -a package to improve cache reutilization based on reordering and grouping techniques [24]. All computational kernels are parallelized covering most of distributed and shared memory systems as well as vectorized and pipelined processors [25].

Enforcing Mass Conservation
The most challenging feature for a good interface capturing method resides in its ability to preserve the mass of the species involved. According to Lohner [11], in VOF methods the mass loss can be associated to reasons such as: interface smearing due to numerical diffusion of the step function, inexact divergence free velocity field and boundary conditions, undulations in the solution of the marking function advection. Level set methods suffer when the marking function loose its signed distance property. These problems have been reported by many authors and are still the subject of several researches. In this work we have followed the procedure proposed in [11] and [12] to overcome mass losses. In this procedure the mass lost/gained are found by comparing the expected value, composed by the initial mass plus the inlet and outlet fluxes, at the end of each time step. Therefore, the values to be added or removed are made proportional to the absolute value of the normal velocity of the interface given by the amount computed from Equation (19) guarantees that the mass correction will act mainly in regions where the interface is moving faster, while keeping the stationary regions untouched. Therefore, the portion of mass correction corresponding to each element is computed and projected onto the global nodes by L 2 projection. In [12], Elias and Coutinho evaluate different test problems with the help of the mass preservation algorithm showing its efficiency.
In order to guarantee the effectiveness of the mass correction method we must have an accurate procedure for computing the volume of the fluid phases. In this work we have adopted a method, described in [12], which makes use of the interface explicitly determined for those elements lying on the interface. To introduce the method, firstly we need to identify the elements as filled, partially filled or empty. This can be done by computing an element key according to the nodal value of the marking function φ. Thus, if φ ≥ 0.5 we assign a value 1, otherwise, it is marked as 0. Finally, the element key is found by summing its nodal values. Note that the empty elements can be skipped from the computation since the void region will be the difference between the domain and the filled region. Moreover, for elements partially filled, the filling volume can be found from the sub-tetrahedra and wedges formed by the interface plane cutting the original element.

Progressive Wave
This case represents the simulation of generation and propagation of a progressive regular wave, from a piston-type wave maker. The computational domain is presented in Figure   1. It has 120 meters length, 20 meters height and 0.5 meters width. Figure 1 also presents the initial condition and initial mesh configuration. Figure 1. Initial condition and mesh -67,320 nodes and 254,052 tetrahedral elements.
The movement is prescribed by a periodic function at the left face, simulating a piston-type wave maker. To avoid mesh entangling we introduce a buffer region where the periodic function describing the nodes displacement suffers a damping effect, given by, where A is the piston displacement amplitude, x and 0 x are position coordinates, ω is the piston angular frequency, t refers to the time and ϕ is the phase angle. In this problem, the piston moves periodically, with amplitude 0.40 m from its original position and the oscillation period is 4 s. All nodes, from the piston up to 20 m, also oscillate gradually according to Equation (20).
The time step was fixed in 0.005 seconds and the simulation time is 40 seconds. Linear GMRES tolerance was set to 10 -3 and 10 -5 for the Navier-Stokes equations and scalar marking function equation, respectively. The Navier-Stokes Inexact-Newton flow solver converges in 3 iterations while the scalar marking function equation converges in 6 iterations. Free surface elevation at t = 0s, t = 8s, t = 16s, t = 24s, t = 32s and t = 40s.
According to linear wave theory, the wave celerity can be computed as the ratio between its length and period [26].  Three positions were chosen to track pressure and velocity fields. Figure 4 shows the monitored points. The analytical solution was calculated from equations given by the potential formulation for regular progressive wave propagation, known as linear theory [26]. The velocity components and pressure equations are shown below: where g denotes gravity. Figure 5 shows the pressure field for all monitored points. The numerical simulation is compared to the linear theory. The transient phase of the pressure time history is due to the initial condition. In all three points, after the initial transient, a good agreement between the numerical and analytical pressure solutions is observed.
Mean pressure is presented in  Velocity field for each one of three monitored points is shown in Figures 6 -11. As it happens in pressure solution, the simulation needs an adjustment time to converge to steady-state solution. Note that, there is a good agreement of the frequency, but the amplitudes are not so close. On the other hand, it is important to observe that linear theory can not cope with viscosity effects and non-linear free surface effects.

Progressive wave interaction with a fixed cylinder
This case refers to the interaction of a progressive wave with a fixed circular cylinder placed on the center of the domain, as shown in Figure 12.   Figure 13 shows the computational domain and the partitions for parallel processing. The computation was performed in parallel (MPI) in a SGI Altix XE 1200, employing 8 of its 88 Xeon cores (2.66 GHz). The maximum inexact-Newton linear tolerance was set to 10 -3 while the linear tolerance for the marking function equation was 10 -5 . The nonlinear loops were stopped after the relative residual as well as the relative step increment decreased 3 orders of magnitude. In addition, a fixed time step of 0.005 was employed. For 50 seconds of simulation it was spent about 70 hours of CPU time. Figure 14 presents five snapshots representing the time solutions at 10, 20, 30, 40 and 50 seconds.  Figure 14. Progressive wave over a cylinder. Free surface elevation at (a) t = 10s, (b) t = 20s, (c) t = 30s, (d) t = 40s, (e) t = 50s. Figure 15 presents the kinetic energy of the fluid. Since the fluid is initially at rest, the kinetic energy is introduced by the piston. The oscillatory movement of the piston, the wavecylinder and wave-opposite face interactions are captured in the kinetic energy oscillatory behavior. Figure 16 shows the water mass preservation with respect to the initial water mass, and reveals how efficiently the mass correction algorithm works. Note that most of the disturbances occur at the initial stage, up to 1 second, which coincides with the moment that the piston pushes the fluid, which is initially at rest. It is important to remind that the mass correction algorithm is based on the velocity field normal to the free surface, making the correction effects more pronounced when the interface is moving faster. The maximum computed water loss was around 0.005% which is one order of magnitude smaller than the adopted nonlinear stopping criteria. It means that material losses are negligible. The history of the horizontal hydrodynamic force exerted by the fluid over the cylinder is illustrated by Figure 17. The amplitude of the linear force evaluated using the MacCamy and Fuchs analytical formulation is equal to 612,299 N [29] and by the WAMIT® code is 591,257 N, [30]. These results, based on linear diffraction theory, are close to the numerical results obtained by the program. The causes of the changes in the time history of the horizontal force (between 37s and 45s for example) could at first be imputed to the wall reflections at the end of the channel, but deserves more investigation

CONCLUSIONS
In this work we have successfully extended our parallel incompressible flow solver to problems involving wave interactions with moving physical boundaries. This goal was reached by introducing Arbitrary Lagrangian-Eulerian formulation into the stabilized finite element formulation. We have validated the present formulation simulating a progressive wave generated from a piston-like wave maker and its interaction with a vertical circular cylinder. The main conclusions are that some important parameters in progressive wave behavior, as celerity, wavelength and height, are in good agreement with existing linear theoretical values. In the problem of the progressive wave interaction with a fixed circular cylinder, we have computed the horizontal force over the cylinder and the kinetic energy of the fluid, which is introduced by the piston movement. We noticed the linear force given by the MacCamy and Fuchs analytical formulation and the WAMIT® code are close to the numerical results obtained in the simulation. Furthermore, we verified in this case that the maximum loss computed is around 0.005%, which evidences the efficiency of global mass conservation procedure.