# Modeling Free Surface Flows Using Stabilized Finite Element Method

Authors: Deepak Garg; Antonella Longo; Paolo Papale
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This work aims to develop a numerical wave tank for viscous and inviscid flows. The Navier-Stokes equations are solved by time-discontinuous stabilized space-time finite element method. The numerical scheme tracks the free surface location using fluid velocity. A segregated algorithm is proposed to iteratively couple the fluid flow and mesh deformation problems. The numerical scheme and the developed computer code are validated over three free surface problems: solitary wave propagation, the collision between two counter moving waves, and wave damping in a viscous fluid. The benchmark tests demonstrate that the numerical approach is effective and an attractive tool for simulating viscous and inviscid free surface flows.


---

## Body

## 1. Introduction

Numerical simulations of flow with moving boundaries have received a growing interest in last few years. Many engineering applications such as ship hydrodynamics, dam break, sloshing in tanks, shallow water, and mold filling [1–3] are modeled with free surface flows. Moreover, many natural risk problems such as tsunami propagation and flood flows require the numerical methods to tackle free surface flows [4, 5]. The continuous variation of the domain shape with time during the numerical solution of such problems is a challenging issue.Moving domain problems are modeled with interface-capturing and interface-tracking techniques. In interface-capturing methods such as the level set [6, 7] and the volume of fluid [8, 9] a fixed grid is used. The location of the moving interface is captured by an additional equation for an artificial scalar field. A big disadvantage of these methods is that the level set method captures the free surface accurately but does not conserve the mass whereas the volume of fluid conserves the mass but does not capture the free surface accurately. Contrarily, in interface-tracking methods the location of moving boundary is directly tracked by the computational mesh. Generic-front-tracking [10], volume-tracking [11, 12], immersed boundary method [13], and smoothed particle hydrodynamics (SPH) method [14] are few examples of interface-tracking methods. Interface-tracking approaches are split in Lagrangian and Arbitrary Lagrangian-Eulerian (ALE). In Lagrangian, all mesh nodes follow the fluid motion, while in ALE method only the boundary nodes move with the fluid velocity and the interior nodes move arbitrarily. Interface-tracking techniques are known to provide greater accuracy over interface-capturing techniques because they define the interface more accurately.In the framework of finite element methods (FEM), ALE method [15–18] and space-time finite element method (STFEM) [19–21] have been successfully used to solve moving domain problems. In traditional ALE method, the governing equations are written in ALE form by introducing the mesh velocity and solved by a semidiscrete approach; i.e., first the spatial discretization is done by standard Galerkin method leaving the time untouched and then time is discretization by some finite difference scheme. In contrast, in STFEM there is no need to include the mesh velocity explicitly in the Navier-Stokes equations as it is mandatory for ALE. In STFEM both space and time are discretized simultaneously generating space-time slabs. In moving domain problems the computation of mesh velocity is treated as an additional subproblem and can be calculated by any mesh deformation approach; e.g., see [22–24].Numerical wave tank has become a popular choice for studying the natural disasters involving free surface flows such as tsunami waves [25–27]. In most cases, the wave tank considers an inviscid flow. Our aim in this work is to develop a robust interface-tracking numerical code for simulating viscous and inviscid flows to be used for the natural hazard problems such as the tsunami wave amplification analysis. We use the stabilized STFEM [28] for solving flow equations. The location of the free surface is tracked by the elastic deformation method [19]. The fluid flow and the mesh motion problem are coupled by an iterative segregated algorithm.The contents of this paper are organized as follows: In Section2 we present flow governing equations and the boundary conditions. In Section 3 the partitioned algorithm is introduced followed by the detailed numerical approaches for the solution of fluid flow and mesh updating problems. The validation of the numerical method is done in Section 4 over selected free surface flow problems. In Section 5 the concluding remarks are drawn.


## 2. Flow Governing Equations

LetΩt∈Rd(d=1,2,3) and (0, ttol) be the spatial and temporal domains, respectively, and let Γt denote the boundary of Ωt. The Navier-Stokes equations in conservation variables U=[ρ,ρvi,ρe]′ are(1)U,t+Fi,i=Ei,i+Swhere Fi, Ei, and S are the Euler flux, viscous flux, and source vectors, respectively, and are defined as(2)Fi=ρviρvjvi+pδjiρvie+pviEi=0τjiτijvj-qiS=0ρbjρbiui+rwhere, ρ, p, vi, τij, e, qi, bi, and r are the density, pressure, ith component of velocity, ijth component of viscous stress tensor, total energy per unit mass, ith component of heat flux vector, ith component of body force, and heat source, respectively. δij is the Kronecker delta function. The viscous stress tensor components are defined as(3)τij=μvi,j+vj,i-23μvi,iδijwhere μ is the viscosity. The total energy per unit mass is defined as e=cvT+v2/2, where cv is the specific heat at constant volume, T is the temperature, and v={vi} is velocity. The heat flux is defined as qi=-κT,i, where κ is the thermal conductivity.By change of variables, (1) can be converted from conservation variables U into pressure primitive variables Y=[p,v,T]′:(4)A0Y,t+AiY,i=KijY,j,i+Swhere A0=U,Y; Ai=Fi,Y is the ith Euler Jacobian matrix; K=[Kij] is the diffusivity matrix with KijY,j=Ei. The above equations are supplemented with an appropriate equation of state relating pressure to other variables. Most compressible flows can be modeled with the ideal gas law, p=ρRT, where R is the specific gas constant. For incompressible flows the density is set to a constant which implies that the isothermal compressibility and isobaric expansiveness are zero. For isothermal flows the energy equation can be excluded from the computations by setting a constant value of temperature. The equations are completed by the following boundary conditions. The boundary Γt is decomposed as Γt=Γn∪Γd∪Γf, where Γn, Γd, and Γf are slip, no-slip, and free surface parts of the boundary, respectively. For (0,ttot) the boundary conditions read(5)v·n=0andt·τn=0onΓn(6)v=vdonΓd(7)-pI+τn=0onΓfwhere vd is the given velocity vector; t and n are the tangential and outward normal unit vectors on the boundary. At free surface, no restriction is imposed on the velocity and the boundary conditions are similar to those normally imposed at an outlet.


## 3. Numerical Scheme

In free surface problems, the location of the moving boundary is unknown within the fluid mechanics problem and must be determined together with the solution of flow equations. The geometry and flow field evolve in time which makes the problem time-dependent and nonlinear. We decouple the fluid flow and the mesh evolution problems and solve iteratively through a partitioned approach. For each time step, the algorithm is composed of the following steps:(1)
Solve the Navier-Stokes equations (4) on a given domain Ω(tn) at time tn.(2)
Solve mesh deformation problem by assigning fluid velocity at the free surface.(3)
Update the interior mesh and go to next time step.


### 3.1. Solution of Navier-Stokes Equations

The Navier-Stokes equations are solved for pressure primitive variablesY which are applicable on both compressible and incompressible flows [29]. We use time-discontinuous stabilized space-time finite element method to solve the equations. The numerical technique allows using the same order interpolation functions for all solution variables. In space-time finite element method both space and time are discretized simultaneously by taking a tensor product of basis functions for the spatial domain and one-dimensional basis function in time direction. We consider piecewise continuous approximation in space and discontinuous approximation in time. Let Q=Ωt×(0,ttot) be the space-time domain with boundary P=Γt×(0,ttot). Let the time interval (0,ttot) be subdivided into an ordered series of time levels 0=t0<t1<⋯<tn<tn+1<⋯<tN=ttot. Denoting the nth time interval as In=(tn,tn+1), a space-time slab is defined as Qn=Ωn×In with boundary Pn=Γn×In, where Ωn is the spatial domain at time t=tn. Each space-time slab is divided into (nel)n elements Qne=Ωne×In having boundary Pne=Γne×In with discontinuous basis functions across the slab boundaries in time direction. The space-time mesh is considered structured in time and the space-time slabs are spatially deformed. The mesh for the upper surface of the slab at tn+1 is obtained by deforming the mesh nodes of the lower surface at tn.The space-time mapping [30] from a physical space-time element Qne to the master space-time element Q~ne=Ω~ne×In is defined as(8)x=∑i=12∑a=1nelnodesNaξTiθxait=∑i=12Tiθtiwhere subscript a denotes the node index and superscript i denotes the time level of a space-time element (i=1 for t=tn and i=2 for t=tn+1); x and ξ are the spatial parts; t and θ are the temporal parts; N(ξ) and T(θ) are, respectively, the bilinear spatial and linear temporal shape functions. In case of free surface problems, the position xa2 of the mesh nodes at time tn+1 is not known and is obtained by moving the vertices xa1 to their new position with deformation of the mesh nodes.(9)xa2=xa1+vamΔtwhere vam is the mesh velocity of node a and Δt=tn+1-tn. The Jacobian matrix for the space-time mapping ϕ over Qne and its inverse transpose are computed as(10)JST=∂x,t∂ξ,θ=x,ξx,θ′t,ξt,θ(11)JST-1′=x,ξ-1′θ,x′ξ,tθ,t=x,ξ-1′0′-vmx,ξ-1′2ΔtThe derivatives of the shape functions in physical and reference coordinate systems satisfy the following relationship:(12)NaTi,xNaTi,t=JST-1′Na,ξTiNaTi,θThe fluid mesh velocity vm determines the working frame of reference. If it is zero then we are in Eulerian frame; if it is equal to the fluid particle velocity then we are in Lagrangian frame; if it is equal to the fluid velocity for the moving boundaries and arbitrary for the internal domain then we recover the ALE frame.Remark 1.
The space-time mapping is similar to the ALE mapping except for the temporal derivative terms. In ALE method the time is kept continuous in course of the spatial discretization; therefore the space-time Jacobian coincides with the spatial Jacobian. However, in space-time FEM the time is discretized simultaneously with space which results in a product of spatial Jacobian and a time difference term.The space-time finite element weighted residual formulation of (4) is written as follows [28]: given a trial function space Tnh and weighting function space Wnh, within each Qn,n=0,1,⋯,N-1, find Yh∈Tnh such that ∀Wh∈Wnh the following relation is satisfied(13)∫Qn-W,th·UYh-W,ih·FiYh+W,ih·EiYh-Wh·SdQn+∫Ωtn+1-Whtn+1-·UYhtn+1-dΩ-∫Ωtn+Whtn+·UYhtn-dΩ+∑e=1neln∫QneL′Wh·τLYh-SdQn-∫PnWh·UYhn^·vmdPn-∫PnWh·EiYh-FiYhn^dPn=0In the above formulation, the first and last two integrals are the Galerkin terms obtained from integration by parts. The second and third integrals are the jump terms to enforce continuity of the solution between consecutive slabs. The jump terms facilitate the projection of the solution in case of remeshing. The fourth integral is the least-square stabilization term with differential operator L=A0(∂/∂t)+Ai(∂/∂xi) and L′=A0′(∂/∂t)+Ai′(∂/∂xi). LYh-S is the residual of the equations and τ is the stabilization matrix for pressure primitive variables Yh [29, 31].The space-time discretization can be done as either constant in time or linear in time [28]. In constant in time, the interpolation functions are considered piecewise linear in space and constant in time, while in linear-in-time interpolation functions are linear both in space and in time. For linear-in-time functions within the nth space-time slab Qn, the finite element trial solution Yh and the weighting function Wh are defined as(14)Yh=∑a=1nnpnNaT2Ya2+NaT1Ya1(15)Wh=∑a=1nnpnNaT2Wa2+NaT1Wa1where Ya2 and Ya1 are the vectors of unknowns at node a at times tn+1- and tn+, respectively. Analogously, Wa2 and Wa1 are the vectors of nodal weighting function at times tn+1- and tn+, respectively. (nnp)(n) is the number of nodal points in nth space-time slab. The substitution of trial and weighting functions in (13) results in two systems of nonlinear equations which are solved by implicit third-order predictor-corrector method [28].Remark 2.
The computation cost of the algorithm can be reduced by implementing the constant-in-time discretization approximation which is first-order time accurate and consists of a linear system having half number of unknowns in comparison to linear-in-time approximation. For constant-in-time the values oft,θ should be set to Δt in (10).


### 3.2. Mesh Update

In free surface problems, the location of the free surface is calculated from the fluid velocityv at each time step of the computation. Once the displacement of the moving surface is known, the internal nodes can be moved to reduce the distortion of the mesh elements. In order to evaluate the new location of the interior mesh nodes, a mesh deformation technique is used. Mesh deformation can be handled arbitrarily by matching the normal velocity of the mesh with the normal velocity of the moving interface. This requires the computation of a unit normal vector at each node of the moving surface. Since the finite element representation is only C0 continuous there is no unique normal at a vertex. Often, a separate normal can be defined for each element that meets at that point and an arithmetic average of the surrounding normals is taken as the normal vector at the mesh node [32]. To avoid this problem as an alternative we match both normal and tangential components of the velocity at the free surface. For mesh deformation, several methods have been proposed in the literature such as Laplacian [22], biharmonic [23], spring [33], transfinite interpolation [34], and the elastic deformation methods [35, 36]. We use the elastic deformation method which has been effectively used in [19, 24]. In elastic deformation method, the computational grid is considered as an elastic body and is deformed subject to applied boundary conditions. For nth space-time slab, the partial differential equations governing the displacement of the mesh are(16)∇·σ=0onΩtnwhere σ is the Cauchy stress tensor which is dependent upon the strain tensor ϵ by the Hooke’s law for isotropic homogeneous medium.(17)σ=λtr⁡ϵI+2μϵϵ=∇um+∇umT2where um is the mesh displacement in nth space-time slab; λ and μ are the Lamé coefficients which can be expressed in terms of modulus of elasticity E and Poisson ratio ν:(18)λ=νE1+ν1-2νμ=E21+νFor ν→0.5 the material is nearly incompressible and the elastic problem becomes ill-conditioned. For incompressible materials a neo-Hookean elasticity model can be used [37]. The mesh equation (16) is solved by applying the following boundary conditions:(19)um=0onΓd(20)um·n^=0onΓn(21)um=vΔtonΓfwhere v is the known fluid velocity. In order to avoid excessive distortion and to maintain a better mesh quality near the free surface, the Jacobian based stiffening approach is employed [24, 36].


## 3.1. Solution of Navier-Stokes Equations

The Navier-Stokes equations are solved for pressure primitive variablesY which are applicable on both compressible and incompressible flows [29]. We use time-discontinuous stabilized space-time finite element method to solve the equations. The numerical technique allows using the same order interpolation functions for all solution variables. In space-time finite element method both space and time are discretized simultaneously by taking a tensor product of basis functions for the spatial domain and one-dimensional basis function in time direction. We consider piecewise continuous approximation in space and discontinuous approximation in time. Let Q=Ωt×(0,ttot) be the space-time domain with boundary P=Γt×(0,ttot). Let the time interval (0,ttot) be subdivided into an ordered series of time levels 0=t0<t1<⋯<tn<tn+1<⋯<tN=ttot. Denoting the nth time interval as In=(tn,tn+1), a space-time slab is defined as Qn=Ωn×In with boundary Pn=Γn×In, where Ωn is the spatial domain at time t=tn. Each space-time slab is divided into (nel)n elements Qne=Ωne×In having boundary Pne=Γne×In with discontinuous basis functions across the slab boundaries in time direction. The space-time mesh is considered structured in time and the space-time slabs are spatially deformed. The mesh for the upper surface of the slab at tn+1 is obtained by deforming the mesh nodes of the lower surface at tn.The space-time mapping [30] from a physical space-time element Qne to the master space-time element Q~ne=Ω~ne×In is defined as(8)x=∑i=12∑a=1nelnodesNaξTiθxait=∑i=12Tiθtiwhere subscript a denotes the node index and superscript i denotes the time level of a space-time element (i=1 for t=tn and i=2 for t=tn+1); x and ξ are the spatial parts; t and θ are the temporal parts; N(ξ) and T(θ) are, respectively, the bilinear spatial and linear temporal shape functions. In case of free surface problems, the position xa2 of the mesh nodes at time tn+1 is not known and is obtained by moving the vertices xa1 to their new position with deformation of the mesh nodes.(9)xa2=xa1+vamΔtwhere vam is the mesh velocity of node a and Δt=tn+1-tn. The Jacobian matrix for the space-time mapping ϕ over Qne and its inverse transpose are computed as(10)JST=∂x,t∂ξ,θ=x,ξx,θ′t,ξt,θ(11)JST-1′=x,ξ-1′θ,x′ξ,tθ,t=x,ξ-1′0′-vmx,ξ-1′2ΔtThe derivatives of the shape functions in physical and reference coordinate systems satisfy the following relationship:(12)NaTi,xNaTi,t=JST-1′Na,ξTiNaTi,θThe fluid mesh velocity vm determines the working frame of reference. If it is zero then we are in Eulerian frame; if it is equal to the fluid particle velocity then we are in Lagrangian frame; if it is equal to the fluid velocity for the moving boundaries and arbitrary for the internal domain then we recover the ALE frame.Remark 1.
The space-time mapping is similar to the ALE mapping except for the temporal derivative terms. In ALE method the time is kept continuous in course of the spatial discretization; therefore the space-time Jacobian coincides with the spatial Jacobian. However, in space-time FEM the time is discretized simultaneously with space which results in a product of spatial Jacobian and a time difference term.The space-time finite element weighted residual formulation of (4) is written as follows [28]: given a trial function space Tnh and weighting function space Wnh, within each Qn,n=0,1,⋯,N-1, find Yh∈Tnh such that ∀Wh∈Wnh the following relation is satisfied(13)∫Qn-W,th·UYh-W,ih·FiYh+W,ih·EiYh-Wh·SdQn+∫Ωtn+1-Whtn+1-·UYhtn+1-dΩ-∫Ωtn+Whtn+·UYhtn-dΩ+∑e=1neln∫QneL′Wh·τLYh-SdQn-∫PnWh·UYhn^·vmdPn-∫PnWh·EiYh-FiYhn^dPn=0In the above formulation, the first and last two integrals are the Galerkin terms obtained from integration by parts. The second and third integrals are the jump terms to enforce continuity of the solution between consecutive slabs. The jump terms facilitate the projection of the solution in case of remeshing. The fourth integral is the least-square stabilization term with differential operator L=A0(∂/∂t)+Ai(∂/∂xi) and L′=A0′(∂/∂t)+Ai′(∂/∂xi). LYh-S is the residual of the equations and τ is the stabilization matrix for pressure primitive variables Yh [29, 31].The space-time discretization can be done as either constant in time or linear in time [28]. In constant in time, the interpolation functions are considered piecewise linear in space and constant in time, while in linear-in-time interpolation functions are linear both in space and in time. For linear-in-time functions within the nth space-time slab Qn, the finite element trial solution Yh and the weighting function Wh are defined as(14)Yh=∑a=1nnpnNaT2Ya2+NaT1Ya1(15)Wh=∑a=1nnpnNaT2Wa2+NaT1Wa1where Ya2 and Ya1 are the vectors of unknowns at node a at times tn+1- and tn+, respectively. Analogously, Wa2 and Wa1 are the vectors of nodal weighting function at times tn+1- and tn+, respectively. (nnp)(n) is the number of nodal points in nth space-time slab. The substitution of trial and weighting functions in (13) results in two systems of nonlinear equations which are solved by implicit third-order predictor-corrector method [28].Remark 2.
The computation cost of the algorithm can be reduced by implementing the constant-in-time discretization approximation which is first-order time accurate and consists of a linear system having half number of unknowns in comparison to linear-in-time approximation. For constant-in-time the values oft,θ should be set to Δt in (10).


## 3.2. Mesh Update

In free surface problems, the location of the free surface is calculated from the fluid velocityv at each time step of the computation. Once the displacement of the moving surface is known, the internal nodes can be moved to reduce the distortion of the mesh elements. In order to evaluate the new location of the interior mesh nodes, a mesh deformation technique is used. Mesh deformation can be handled arbitrarily by matching the normal velocity of the mesh with the normal velocity of the moving interface. This requires the computation of a unit normal vector at each node of the moving surface. Since the finite element representation is only C0 continuous there is no unique normal at a vertex. Often, a separate normal can be defined for each element that meets at that point and an arithmetic average of the surrounding normals is taken as the normal vector at the mesh node [32]. To avoid this problem as an alternative we match both normal and tangential components of the velocity at the free surface. For mesh deformation, several methods have been proposed in the literature such as Laplacian [22], biharmonic [23], spring [33], transfinite interpolation [34], and the elastic deformation methods [35, 36]. We use the elastic deformation method which has been effectively used in [19, 24]. In elastic deformation method, the computational grid is considered as an elastic body and is deformed subject to applied boundary conditions. For nth space-time slab, the partial differential equations governing the displacement of the mesh are(16)∇·σ=0onΩtnwhere σ is the Cauchy stress tensor which is dependent upon the strain tensor ϵ by the Hooke’s law for isotropic homogeneous medium.(17)σ=λtr⁡ϵI+2μϵϵ=∇um+∇umT2where um is the mesh displacement in nth space-time slab; λ and μ are the Lamé coefficients which can be expressed in terms of modulus of elasticity E and Poisson ratio ν:(18)λ=νE1+ν1-2νμ=E21+νFor ν→0.5 the material is nearly incompressible and the elastic problem becomes ill-conditioned. For incompressible materials a neo-Hookean elasticity model can be used [37]. The mesh equation (16) is solved by applying the following boundary conditions:(19)um=0onΓd(20)um·n^=0onΓn(21)um=vΔtonΓfwhere v is the known fluid velocity. In order to avoid excessive distortion and to maintain a better mesh quality near the free surface, the Jacobian based stiffening approach is employed [24, 36].


## 4. Numerical Examples

In order to demonstrate the effectiveness of the model and to validate the computational code, we present three numerical examples: solitary wave propagation, the collision of two waves, and wave damping. The first two examples are for inviscid flow and the third one is for viscous flow. In order to compute the inviscid flow the fluid viscosity is set to zero in (2). In case of no heat exchange this results in a zero diffusion flux vector (i.e., Ei=0). For all simulations, the spatial meshes are made of quadrilateral elements and the time slabs are formed simply by moving the spatial meshes in time. The mesh deformation is computed with the elastic deformation method (16) by setting Young’s modulus E = 1000 Pa and the Poisson ratio ν = 0.35.


### 4.1. Solitary Wave Propagation

This numerical problem is a classical check for the free surface flow models and has been studied in many works [16, 21, 38, 39]. The problem consists in the motion of an inviscid fluid in a tank with three fixed walls (bottom and both vertical walls). The upper surface is free to deform with gravitational force controlling the motion of the wave. The model problem is presented in Figure 1. The initial profile of the wave at the free surface η is given by(22)η=d+Hsech2⁡3H4d3xThe domain parameters are d = 10 m and H = 2 m. The simulation is performed over a structured mesh made of 6400 quadrangle elements. Fluid density is set to 1000 kg/m3. The gravitational acceleration is set to 9.81 m/s2. The initial conditions for this problem are set according to the approximations given by Laitone [40] for an infinite domain. The initial velocity is set as(23)v1=gdHdsech2⁡3H4d3x(24)v2=3gdHd3/2ysech2⁡3H4d3xtanh⁡3H4d3xAt the free surface, the boundary condition (7) is imposed. Both top corner points belong to the free surface. The fluid mesh nodes on the vertical walls of the container move in response to the motion of the top corner points. The time step is set to 10-3 s.Figure 1
Wave propagation: problem description.Figure2 shows the velocity magnitude at various times. The wave reaches the right vertical wall at time t = 7.6 s. The simulation is carried out until the water wave returns to its initial position at x = 0. The reflected wave takes approximately the same time to travel back to its initial position. The total run-up height d+Rmax against the right wall is 14.28 m, where Rmax is the maximum run-up height and can be calculated by the following relation [41]:(25)Rmaxd=2Hd+12Hd2+OHd3The numerical results are in close agreement with the analytical one proposed in [41]. The total run-up height at the right wall computed by our model is compared with [16] in Figure 3. The spatial and temporal convergence has been studied by considered two grids and several time step sizes. Table 1 lists the L2 norm error results.Table 1
L 2 norm solution error.
Grid Time step L 2 error 240 × 15 0.001 1.3583e-3 240 × 15 0.005 9.7483e-2 240 × 15 0.01 7.9872e-1 320 × 20 0.001 5.8368e-4 320 × 20 0.005 2.9428e-2 320 × 20 0.01 3.8183e-1Figure 2
Wave propagation: absolute velocity (m/s) from top to bottom att = 0, 1, 3, 7.6, 9, 12, and 15 s.Figure 3
Wave propagation: comparison of numerical results with [16].


### 4.2. Collision between Two Waves

This problem investigates a head-on collision between two solitary waves and was studied by Su and Mirie [42] for an inviscid, incompressible, homogeneous fluid by a perturbation method up to the third order of approximation. The problem was reformulated in [43] for an overtaking collision between two solitary waves and the results on phase shifts and the time history of interaction were verified with the theoretical solution of the Korteweg-de Vries equation. In [44] the problem was numerically investigated by using a velocity-vorticity based Arbitrary Lagrangian-Eulerian method. We simulate the collision of two solitary waves with two different initial amplitudes: 3 m on the left-hand side and 2 m on the right-hand side. The density and gravity are set to 1000 kg/m3 and 9.81 m/s2, respectively. The length of the rectangular domain is 32 d with d = 10 m. Initial velocities are set using Laitone’s relations (23) and (24). The initial profiles for the velocity are shown in Figure 4. The computational grid is made of 10000 elements. The boundary conditions for the fluid flow and the mesh problem are same as those for the wave damping of viscous fluid problem. The time step is 10-3 s.Figure 4
Wave collision: the initial profiles for horizontal (v1) and vertical velocity (v2) (m/s).Figure5 shows the space-time characteristic curves for free surface elevation. The maximum amplitude (run-up) during the collision can be evaluated as [43](26)Maximum  run-up=AL+BR+12ALBR+14ALBRAL+BRwhere AL, BR are the initial amplitudes of the solitary waves on the left and the right-hand sides, respectively. The simulation results give a maximum run-up height of 15.491 m at time t = 7.5 s, which is close to the calculated value of 15.500 m by (26).Figure 5
Wave collision: space-time characteristic curves for free surface elevation.


### 4.3. Wave Damping of a Viscous Fluid

This problem consists in solving the small-amplitude motion of the free surface of a viscous fluid in a rectangular tank. The initial profile of the free surface is given by(27)hx=1.5+a0sin⁡π0.5-xwhere a0 is the amplitude of the initial sinusoidal perturbation of the movement. The problem has been investigated in [45, 46]. The simulation domain for the problem is shown in Figure 6. The initial perturbation a0 is set to 0.01 m. The dynamic viscosity and the density of the fluid are set to 0.01 Pa·s and 1 kg/m3, respectively. The gravitational acceleration is set to unity. The computational mesh is composed of 50 × 50 quadrangle elements. The initial conditions for the fluid are zero velocity and a hydrostatic pressure profile. The fluid is allowed to slip at all walls of the container. The time step is set to 10-3 s.Figure 6
Viscous sloshing: schematic description of the problem.Under the effect of gravity, the left part of fluid tends to move downwards and enforces the right part of the fluid to move upwards. The first phase of this movement lasts for the first 1.7 s and the fluid reaches its minimum amplitude of 1.4992 m at the left wall. After that, the motion is reversed and continues periodically. The flow is dominated by the viscous forces which are responsible for the damping of the movement. For small-amplitude waves, the initial-value problem has been solved analytically in [47]. For free surface y=a(t)g(x,y) with g satisfying (∂x2+∂y2+k2)g=0, the solution for one fluid is given by(28)at=4ν2k4a0fνk2t8ν2k4+ω02+∑i=14ziZiω02a0zi2-νk2ezi2-νk2tfzitwhere ν is the kinematic viscosity of the fluid; k is the wave number; ω02=gk is the inviscid natural frequency; each zi is the root of the following algebraic equation:(29)z4+2k2νz2+4k6ν3z+ν2k4+ω02=0and Z1=(z2-z1)(z3-z1)(z4-z1) with Z2, Z3, Z4 obtained by circular permutation of the indices and f(.) is the complementary error function defined as(30)fx=2π∫x∞e-t2dtFigure 7 compares the simulation results with the analytical solution [47] for the amplitude of the top left node over time. Numerical results match remarkably well with the analytic solution.Figure 7
Viscous sloshing: time evolution of the amplitude at the top left corner point.


## 4.1. Solitary Wave Propagation

This numerical problem is a classical check for the free surface flow models and has been studied in many works [16, 21, 38, 39]. The problem consists in the motion of an inviscid fluid in a tank with three fixed walls (bottom and both vertical walls). The upper surface is free to deform with gravitational force controlling the motion of the wave. The model problem is presented in Figure 1. The initial profile of the wave at the free surface η is given by(22)η=d+Hsech2⁡3H4d3xThe domain parameters are d = 10 m and H = 2 m. The simulation is performed over a structured mesh made of 6400 quadrangle elements. Fluid density is set to 1000 kg/m3. The gravitational acceleration is set to 9.81 m/s2. The initial conditions for this problem are set according to the approximations given by Laitone [40] for an infinite domain. The initial velocity is set as(23)v1=gdHdsech2⁡3H4d3x(24)v2=3gdHd3/2ysech2⁡3H4d3xtanh⁡3H4d3xAt the free surface, the boundary condition (7) is imposed. Both top corner points belong to the free surface. The fluid mesh nodes on the vertical walls of the container move in response to the motion of the top corner points. The time step is set to 10-3 s.Figure 1
Wave propagation: problem description.Figure2 shows the velocity magnitude at various times. The wave reaches the right vertical wall at time t = 7.6 s. The simulation is carried out until the water wave returns to its initial position at x = 0. The reflected wave takes approximately the same time to travel back to its initial position. The total run-up height d+Rmax against the right wall is 14.28 m, where Rmax is the maximum run-up height and can be calculated by the following relation [41]:(25)Rmaxd=2Hd+12Hd2+OHd3The numerical results are in close agreement with the analytical one proposed in [41]. The total run-up height at the right wall computed by our model is compared with [16] in Figure 3. The spatial and temporal convergence has been studied by considered two grids and several time step sizes. Table 1 lists the L2 norm error results.Table 1
L 2 norm solution error.
Grid Time step L 2 error 240 × 15 0.001 1.3583e-3 240 × 15 0.005 9.7483e-2 240 × 15 0.01 7.9872e-1 320 × 20 0.001 5.8368e-4 320 × 20 0.005 2.9428e-2 320 × 20 0.01 3.8183e-1Figure 2
Wave propagation: absolute velocity (m/s) from top to bottom att = 0, 1, 3, 7.6, 9, 12, and 15 s.Figure 3
Wave propagation: comparison of numerical results with [16].


## 4.2. Collision between Two Waves

This problem investigates a head-on collision between two solitary waves and was studied by Su and Mirie [42] for an inviscid, incompressible, homogeneous fluid by a perturbation method up to the third order of approximation. The problem was reformulated in [43] for an overtaking collision between two solitary waves and the results on phase shifts and the time history of interaction were verified with the theoretical solution of the Korteweg-de Vries equation. In [44] the problem was numerically investigated by using a velocity-vorticity based Arbitrary Lagrangian-Eulerian method. We simulate the collision of two solitary waves with two different initial amplitudes: 3 m on the left-hand side and 2 m on the right-hand side. The density and gravity are set to 1000 kg/m3 and 9.81 m/s2, respectively. The length of the rectangular domain is 32 d with d = 10 m. Initial velocities are set using Laitone’s relations (23) and (24). The initial profiles for the velocity are shown in Figure 4. The computational grid is made of 10000 elements. The boundary conditions for the fluid flow and the mesh problem are same as those for the wave damping of viscous fluid problem. The time step is 10-3 s.Figure 4
Wave collision: the initial profiles for horizontal (v1) and vertical velocity (v2) (m/s).Figure5 shows the space-time characteristic curves for free surface elevation. The maximum amplitude (run-up) during the collision can be evaluated as [43](26)Maximum  run-up=AL+BR+12ALBR+14ALBRAL+BRwhere AL, BR are the initial amplitudes of the solitary waves on the left and the right-hand sides, respectively. The simulation results give a maximum run-up height of 15.491 m at time t = 7.5 s, which is close to the calculated value of 15.500 m by (26).Figure 5
Wave collision: space-time characteristic curves for free surface elevation.


## 4.3. Wave Damping of a Viscous Fluid

This problem consists in solving the small-amplitude motion of the free surface of a viscous fluid in a rectangular tank. The initial profile of the free surface is given by(27)hx=1.5+a0sin⁡π0.5-xwhere a0 is the amplitude of the initial sinusoidal perturbation of the movement. The problem has been investigated in [45, 46]. The simulation domain for the problem is shown in Figure 6. The initial perturbation a0 is set to 0.01 m. The dynamic viscosity and the density of the fluid are set to 0.01 Pa·s and 1 kg/m3, respectively. The gravitational acceleration is set to unity. The computational mesh is composed of 50 × 50 quadrangle elements. The initial conditions for the fluid are zero velocity and a hydrostatic pressure profile. The fluid is allowed to slip at all walls of the container. The time step is set to 10-3 s.Figure 6
Viscous sloshing: schematic description of the problem.Under the effect of gravity, the left part of fluid tends to move downwards and enforces the right part of the fluid to move upwards. The first phase of this movement lasts for the first 1.7 s and the fluid reaches its minimum amplitude of 1.4992 m at the left wall. After that, the motion is reversed and continues periodically. The flow is dominated by the viscous forces which are responsible for the damping of the movement. For small-amplitude waves, the initial-value problem has been solved analytically in [47]. For free surface y=a(t)g(x,y) with g satisfying (∂x2+∂y2+k2)g=0, the solution for one fluid is given by(28)at=4ν2k4a0fνk2t8ν2k4+ω02+∑i=14ziZiω02a0zi2-νk2ezi2-νk2tfzitwhere ν is the kinematic viscosity of the fluid; k is the wave number; ω02=gk is the inviscid natural frequency; each zi is the root of the following algebraic equation:(29)z4+2k2νz2+4k6ν3z+ν2k4+ω02=0and Z1=(z2-z1)(z3-z1)(z4-z1) with Z2, Z3, Z4 obtained by circular permutation of the indices and f(.) is the complementary error function defined as(30)fx=2π∫x∞e-t2dtFigure 7 compares the simulation results with the analytical solution [47] for the amplitude of the top left node over time. Numerical results match remarkably well with the analytic solution.Figure 7
Viscous sloshing: time evolution of the amplitude at the top left corner point.


## 5. Conclusions

Our objective in this work was to develop a numerical scheme for simulating viscous and inviscid free surface flows for the numerical study of the wave tank. We used a partitioned algorithm to iteratively solve the fluid flow and the mesh deformation problems. The numerical technique introduced mesh velocity in an implicit manner and avoided the need to write the flow equations in the ALE form a priori. We presented the detailed formulation of the numerical scheme covering space-time mapping and discretization strategies. We successfully validated our numerical code on one viscous and two inviscid flow test cases. For the simulations, our choice of initial meshes preserved a good quality for the entire simulation times and allowed avoiding remeshing for any arbitrary time interval. The accuracy of the numerical results leads to concluding that the numerical approach is highly efficient and stable for a variety of free flow problems.The numerical scheme is clearly not applicable for free surface flows involving the wave breaking and the fluid separation. In such cases, free surface undergoes a large deformation and a portion of the fluid split from the main body leading to the breakdown of the numerical scheme. Such problems can be modeled with multicomponent/multiphase flows. Future developments will concern the extension of the developed algorithm to multicomponent flows and fully three-dimensional cases which are relevant and crucial for the natural risk problems.



---

*Source: 6154251-2018-06-11.xml*