Mathematical Problems in Engineering

Volume 2018, Article ID 6154251, 9 pages

https://doi.org/10.1155/2018/6154251

## Modeling Free Surface Flows Using Stabilized Finite Element Method

Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Pisa, Via Uguccione della Faggiola 32, I-56126 Pisa, Italy

Correspondence should be addressed to Deepak Garg; ti.vgni@grag.kapeed

Received 3 January 2018; Revised 14 May 2018; Accepted 20 May 2018; Published 11 June 2018

Academic Editor: Luis J. Yebra

Copyright © 2018 Deepak Garg et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 [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 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.