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.

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 [13] 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 [1518] and space-time finite element method (STFEM) [1921] 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 [2224].

Numerical wave tank has become a popular choice for studying the natural disasters involving free surface flows such as tsunami waves [2527]. 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 Section 2 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 and (0, ) be the spatial and temporal domains, respectively, and let denote the boundary of . The Navier-Stokes equations in conservation variables arewhere , , and are the Euler flux, viscous flux, and source vectors, respectively, and are defined aswhere, , , , , , , , and are the density, pressure, component of velocity, component of viscous stress tensor, total energy per unit mass, component of heat flux vector, component of body force, and heat source, respectively. is the Kronecker delta function. The viscous stress tensor components are defined aswhere is the viscosity. The total energy per unit mass is defined as , where is the specific heat at constant volume, is the temperature, and is velocity. The heat flux is defined as , where is the thermal conductivity.

By change of variables, (1) can be converted from conservation variables into pressure primitive variables :where ; is the th Euler Jacobian matrix; is the diffusivity matrix with . 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, , where 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 is decomposed as , where , , and are slip, no-slip, and free surface parts of the boundary, respectively. For the boundary conditions readwhere is the given velocity vector; and 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 at time .(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 variables 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 be the space-time domain with boundary . Let the time interval be subdivided into an ordered series of time levels . Denoting the time interval as , a space-time slab is defined as with boundary , where is the spatial domain at time . Each space-time slab is divided into elements having boundary 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 is obtained by deforming the mesh nodes of the lower surface at .

The space-time mapping [30] from a physical space-time element to the master space-time element is defined aswhere subscript denotes the node index and superscript denotes the time level of a space-time element (=1 for and =2 for ); and are the spatial parts; and are the temporal parts; and are, respectively, the bilinear spatial and linear temporal shape functions. In case of free surface problems, the position of the mesh nodes at time is not known and is obtained by moving the vertices to their new position with deformation of the mesh nodes.where is the mesh velocity of node and . The Jacobian matrix for the space-time mapping over and its inverse transpose are computed asThe derivatives of the shape functions in physical and reference coordinate systems satisfy the following relationship:The fluid mesh velocity 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 and weighting function space , within each , find such that the following relation is satisfiedIn 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 and . is the residual of the equations and is the stabilization matrix for pressure primitive variables [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 space-time slab , the finite element trial solution and the weighting function are defined aswhere and are the vectors of unknowns at node at times and , respectively. Analogously, and are the vectors of nodal weighting function at times and , respectively. is the number of nodal points in th 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 of should be set to in (10).

3.2. Mesh Update

In free surface problems, the location of the free surface is calculated from the fluid velocity 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 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 space-time slab, the partial differential equations governing the displacement of the mesh arewhere is the Cauchy stress tensor which is dependent upon the strain tensor by the Hooke’s law for isotropic homogeneous medium.where is the mesh displacement in space-time slab; and are the Lamé coefficients which can be expressed in terms of modulus of elasticity and Poisson ratio :For 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:where 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., ). 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 byThe domain parameters are = 10 m and = 2 m. The simulation is performed over a structured mesh made of 6400 quadrangle elements. Fluid density is set to 1000 kg/. The gravitational acceleration is set to 9.81 . 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 asAt 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 s.

Figure 2 shows the velocity magnitude at various times. The wave reaches the right vertical wall at time = 7.6 . The simulation is carried out until the water wave returns to its initial position at = 0. The reflected wave takes approximately the same time to travel back to its initial position. The total run-up height against the right wall is 14.28 m, where is the maximum run-up height and can be calculated by the following relation [41]:The 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 norm error results.

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/ and 9.81 m/, respectively. The length of the rectangular domain is 32 with = 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 s.

Figure 5 shows the space-time characteristic curves for free surface elevation. The maximum amplitude (run-up) during the collision can be evaluated as [43]where , 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 = 7.5 s, which is close to the calculated value of 15.500 m by (26).

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 bywhere 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 is set to 0.01 m. The dynamic viscosity and the density of the fluid are set to 0.01 Pas and 1 kg/, 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 s.

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 with satisfying , the solution for one fluid is given bywhere is the kinematic viscosity of the fluid; k is the wave number; is the inviscid natural frequency; each is the root of the following algebraic equation:and with , , obtained by circular permutation of the indices and is the complementary error function defined asFigure 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.

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.

Data Availability

The numerical simulations data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The research leading to these results has received funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013) under the project NEMOH, REA, Grant Agreement no. 289976.