Next Article in Journal
Variability in the Chemical Composition of Spring Waters in the Postomia River Catchment (Northwest Poland)
Previous Article in Journal
Effects of Warming on Aquatic Snails and Periphyton in Freshwater Ecosystems with and without Predation by Common Carp
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Development of a Three-Dimensional Hydrodynamic Model Based on the Discontinuous Galerkin Method

State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University, Tianjin 300072, China
*
Author to whom correspondence should be addressed.
Water 2023, 15(1), 156; https://doi.org/10.3390/w15010156
Submission received: 30 November 2022 / Revised: 24 December 2022 / Accepted: 26 December 2022 / Published: 30 December 2022
(This article belongs to the Topic Computational Fluid Dynamics (CFD) and Its Applications)

Abstract

:
Though the discontinuous Galerkin method is attracting more and more applications in many fields due to its local conservation, high-order accuracy, and flexibility for resolving complex geometries, only a few three-dimensional hydrodynamic models based on the discontinuous Galerkin method are present. In this study, a three-dimensional hydrodynamic model with a σ-coordinate system in the vertical direction is developed. This model is discretized in space using the discontinuous Galerkin method and advanced in time with the implicit-explicit Runge–Kutta method. Numerical tests indicate that the developed model is convergent and can obtain better results with a smaller computational time when a higher approximation order is adopted. Other tests with exact solutions also indicate that the developed model can well simulate the vertical circulation under the effect of surface wind stress and the flow under the combined effect of wind stress and Coriolis acceleration terms. The simulation results of tidal flow in part of Bohai Bay, China, indicate that the model can be used for the simulation of tidal wave motion in realistic situations.

1. Introduction

With the assumption of hydrostatic pressure distribution, the three-dimensional shallow water equations (3D SWEs) can be derived from the Reynolds-averaged three-dimensional Navier–Stokes equations. They can predict the vertical distribution of primitive variables and can be used together with advection-diffusion-type equations to simulate the transport and fate of substances, such as temperature, salt, and sediment.
In the last few decades, 3D SWEs have been intensively studied [1], and many hydrodynamic models have been developed and widely used. A brief categorization of these models, based on the numerical methods and the type of grids used in the horizontal direction, is listed in Table 1. For the models at an early stage, the conservative finite difference method on structured grids was usually used. They can be implemented straightforwardly and are computationally efficient, while they approximate the coastlines as staircases and prevent the flexible implementation of variable resolution. Later, unstructured, grid-based hydrodynamic models, such as the finite element and finite volume method, were brought into the main stream, since the coastlines could be accurately represented to give a reasonable simulation at the regional scale [2], and the area of most concern could be simulated using local refinement grids.
In addition to the numerical methods mentioned above, the discontinuous Galerkin (DG) method has received much attention in recent years. It can be viewed as a hybrid between the finite element method and the finite volume method and enjoys most of the strengths of both methods, such as local conservation, compactness, a high order of accuracy, and good application to unstructured meshes. It has been successfully applied in the numerical discretization of two-dimensional shallow water equations [3], two-dimensional non-hydrostatic shallow water equations [4], and Boussinesq-type equations [5].
Table 1. Categorization of the commonly used hydrodynamic model.
Table 1. Categorization of the commonly used hydrodynamic model.
Grid TypeNumerical Method
Finite Difference MethodFinite Volume MethodFinite Element Method
Structured gridROMS [6], POM [7], MoM [8], GETM [9], TRIM [10], DELFT3D [11], EFDC [12], ECOMSED [13], COHERENS [14]MITgcm [15], MOHID [16]
Unstructured grid FVCOM [17], UnTRIM [18], HydroInfo [19]FESOM [20], ICOM, SELFE [21], ADCIRC [22], SCHISM [23], TELEMAC [24]
In the realm of 3D SWEs with DG solutions, Dawson and Aizinger [25] first developed a three-dimensional barotropic model with a z-coordinate system adopted in the vertical direction and presented the corresponding stability analysis [26]. The vertical eddy viscosity terms in the model are explicitly advanced in time. Due to the stiffness of these terms, the time step must be small enough, which limits the applicability of the model. Later, Aizinger et al. [27] made improvements by introducing the baroclinic terms and treating both the vertical eddy viscosity and the vertical diffusion terms implicitly. This model is called UTBEST3D and has been validated with realistic applications. In the same period, Blaise [28] and Comblen [29] developed the three-dimensional baroclinic model named SLIM and adopted it for the simulation of various baroclinic phenomena. Due to the discontinuous nature of the DG method, the water depth between the adjacent elements is generally not continuous such that the lateral boundaries are mismatched, and the post-process measures are needed to smooth the water depth such that the calculation of numerical flux will not be affected. To circumvent the problem of boundary mismatch between adjacent elements, Gandham [30] mapped the time-varying z-coordinate system to a fixed domain by invoking the σ-coordinate transformation and developed a three-dimensional barotropic model. However, only one layer is adopted in the vertical direction, which may lead to serious numerical errors when the water depth is large. Conroy and Kubatko [31] developed a three-dimensional barotropic model with a σ-coordinate system adopted in the vertical direction and developed a fast and accurate method for the calculation of depth-averaged velocity. An arbitrary number of layers can be used in the vertical direction. However, non-linear advective terms are not included, and the vertical eddy viscosity terms are explicitly advanced in time. To the authors’ knowledge, there is currently no DG method-based three-dimensional hydrodynamic model that includes non-linear advective terms, treats the vertical eddy viscosity terms implicitly, and can use an arbitrary number of vertical layers in the σ-coordinate system simultaneously.
In this paper, a three-dimensional hydrodynamic model that includes nonlinear advective terms and can use an arbitrary number of vertical layers in the σ-coordinate system is developed based on the DG method. The paper is organized as follows. In Section 2, we provide the governing equations of the 3D SWEs in the σ-coordinate system. Then, we present the space discretization and time stepping of the model in Section 3. Numerical experiments are given in Section 4. Finally, conclusions are drawn in Section 5.

2. Governing Equations

Using the notations given in Table 2, the governing equations of the three-dimensional hydrodynamic model in the physical domain are as follows [13]:
u x + v y + w z = 0 ,
u t + u u x + u v y + u w z f v = g η x + z K m u z + x K h u x + y K h u y ,
v t + u v x + v v y + v w z + f u = g η y + z K m v z + x K h v x + y K h v y ,
The surface boundary conditions for u , v , and w are
K m u z , K m v z = 1 ρ 0 τ s x , τ s y ,
w = η t + u η η x + v η η y .
Likewise, the bottom boundary conditions for u , v , and w are
w = u b b x + v b b y ,
K m u z , K m v z = C f u c 2 + v c 2 u c , v c ,
where C f is determined by matching a logarithmic bottom layer to the model at height z 0 and calculated as [32]
C f = κ ln z 0 + z 0 b / z 0 b 2 .
To accurately represent the bottom and surface geometry, the σ transformation developed by Philips [33] is adapted
t * = t , x * = x , y * = y , σ = z η x , y , t D x , y , t .
According to the principle of chain differentiation, partial derivatives of the variable ζ = ζ ( x , y , z , t ) in the physical domain can be written in the following form in the computational domain:
ζ t = ζ t * + ζ σ σ t , ζ x = ζ x * + ζ σ σ x , ζ y = ζ y * + ζ σ σ y , ζ z = ζ σ σ z .
Substituting Equation (10) back into Equations (1)–(3), and ignoring the superscript “*”, we can obtain the governing equations in the computational domain as
D t + D u x + D v y + ω σ = 0 ,
U t + E x + G y + H σ = S 0 + S f + S d , v + S d , h .
with
U = D u D v ,   E = D u 2 + 1 / 2 g D 2 z b 2 D u v ,   G = D u v D v 2 + 1 / 2 g D 2 z b 2 ,   H = ω u ω v .
S 0 = g η z b / x g η z b / y ,   S f = f D v f D u .
S d , v = σ K m D 2 D u σ σ K m D 2 D v σ ,   S d , h = x K h D u x + y K h D u y x K h D v x + y K h D v y .
where ω is related to w by
ω = D d σ d t = w z s t + σ D t u z s x + σ D x v z s y + σ D y .
The surface boundary condition in the computational domain is given by
ω = 0 , K m D 2 D u σ , K m D 2 D v σ = 1 ρ 0 τ s x , τ s y ,
and the bottom boundary condition is given by
ω = 0 , K m D 2 D u σ , K m D 2 D v σ = C f u b 2 + v b 2 u b , v b ,
Integrating Equation (11) from bottom to surface and applying the boundary condition about ω leads to the primitive continuity equation as
D t + D U ¯ x + D V ¯ y = 0

3. Numerical Implementation

In this section, we will detail the numerical implementation of our model based on the DG method. The domain partition and the discontinuous polynomial space are first described. Then, the numerical discretization of each part is presented, followed by the time-stepping scheme.

3.1. Domain Partition and Polynomial Space

For an arbitrary three-dimensional domain Ω 3 d , its horizontal projection is denoted as Ω 2 d and is partitioned by E l e 2 d non-overlapping triangular or quadrilateral elements; i.e., Ω 2 d , h = Ω e , e 1 , E l e 2 d . In the vertical direction, Ω 2 d , h is extruded by N L layers to obtain the three-dimensional computational domain Ω 3 d , h , and the computational domain consists of E l e 2 d × N L triangular prisms or quadrangular prisms. In this study, the three-dimensional element obtained by the extrusion of Ω e is denoted as Ω e , L , L 1 , N L , with the bottommost element denoted as Ω e , 1 and the upmost element denoted as Ω e , N L . Figure 1 gives a schematic view of two columns of three-dimensional computational elements with Figure 1a showing the quadrangular prisms and Figure 1b showing the triangular prisms. To illustrate the discontinuous nature of the DG method, the discontinuities between the adjacent computational elements have been exaggerated, and no gap exists in reality.
In addition, we introduce the following finite dimensional space of polynomials:
V h , 3 d = s L 2 Ω : Ω e , L Ω 3 d , h , s P N h , N v Ω e , L ,
where L 2 Ω is the space of the square-integrable functions and P N h , N v Ω e , L is the complete space of polynomials defined over Ω e , L , which is of order at most N h in the horizontal direction and at most N v in the vertical direction. In this study, V h , 3 d is spanned by N p Lagrangian functions, l i x , i 1 , N p , and the dimension, N p , depends on the approximation order and element type. For triangular prisms,
N p = N h + 1 N h + 2 2 × N v + 1 ,
while for quadrilateral prisms,
N p = N h + 1 2 × N v + 1 .
For any two adjacent elements Ω e 1 and Ω e 2 , let ε = Ω e 1 Ω e 2 be the unique interior interface between these two elements, and let n and n + be the outward unit normal for these two elements on ε . For scalar field c V h , 3 d , let c ε ± be the trace of c on interface ε from the interior of Ω e 1 2 . We further define the mean value c and the jump value c on interface ε as
c = c ε + + c ε / 2 , c = c ε + n + + c ε n .
On the boundary of the computational domain, they are defined as
c = c , c = c n ,
where c is the single value and n is the outward unit normal on the boundary edges. We note that the jump of scalar c is a vector, and we further denote its horizontal and vertical components by c h and c σ , respectively.

3.2. Numerical Discretization of Momentum Equations

For Equation (12), both the Coriolis acceleration terms and the vertical eddy viscosity terms are treated implicitly, while the others are treated explicitly.

3.2.1. Convective and Bottom Topography Terms

For the numerical discretization of Equation (12), we consider the following general formulation: find the local approximate solution vector U h V h , 3 d 2 such that for all test functions ϕ V h , 3 d and Ω e , L Ω 3 d , h we have
Ω   e , L U h t + E U h x + G U h y + H U h σ ϕ d x = Ω   e , L ϕ n F U h F * U h , U h + d x + Ω   e , L S 0 U h ϕ d x + D d , h D , K h , U h , ϕ + D d , f v , h D , K m , U h , f , ϕ .
where U h and U h + are the values of the local approximate solution vector U h on the interface of two adjacent elements, and n F U h = n x E U h + n y G U h + n σ H U h and n F * U h , U h + are the flux term and the numerical flux at the element boundary, respectively. For the three-dimensional computational mesh used in this study, the third component of the unit outward normal is zero, i.e., n σ = 0 , at the lateral boundary of the computational elements, the calculation of the numerical flux degenerates to that of the two-dimensional shallow water equations, and the Harten–Lax–van Leer with contact (HLLC) Riemann solver [34] is adopted. As n x = n y = 0 at the top and bottom boundaries of the computational elements, only the vertical convection needs to be considered, and the upwind flux is taken, i.e.,
n F * U h , U h + = 0.5 n σ ω U h / D 1 + sgn ( ω n σ ) + ω + U h + / D 1 + sgn ( ω + n σ + ) ,
where sgn is the sign function. At the top boundary of the computational domain σ = 0 and bottom boundary of the computational domain σ = 1 , the numerical flux is set to zero. Operators D d , h D , K h , U h , f , ϕ and D d , f v , h D , K m , U h , f , ϕ in Equation (25) are the discretization corresponding to the horizontal eddy viscosity terms and the vertical and Coriolis acceleration terms, which will be detailed in the forthcoming parts. The requirement that Equation (25) is valid for all ϕ V h , 3 d and can be imposed by considering Equation (25) for all the basis functions that span the functional space V h , 3 d , i.e., setting ϕ = l i x ,   i 1 , N p .
For the discretization of both the volume and the surface terms, the quadrature-free manner of Hesthaven and Warburton [35] is adopted. We assume that these terms also lie in the finite dimensional space V h , 3 d 2 and can be expressed as
E ( U h ) E h = i = 1 N P E i l i , G ( U h ) G h = i = 1 N P G i l i , H ( U h ) H h = i = 1 N P H i l i ,   S 0 ( U h ) S 0 , h = i = 1 N P S 0 , i l i , F ( U h ) n F h n = f = 1 N f p F h , f n f l f , F * U h , U h + n F h * n = f = 1 N f p F h , f * n f l f ,
where l f f = 1 N f p are the basis functions on the element boundaries and N f p is the total number of basis functions along all the boundaries. According to this approximation, we can obtain the following semi-discrete equation in matrix form:
e , L d U h , e , L d t = S x , e , L E h , e , L S y , e , L G h , e , L S σ , e , L H h , e , L + e , L f F h , e , L F h , e , L * + e , L S 0 , h , e , L ,
where U h , e , L = D u h , e , L , D v h , e , L T are the nodal values of the conservative variables at the N p interpolation points; E h , e , L , G h , e , L and H h , e , L are the non-linear flux terms at the N p interpolation points; F h , e , L = E h , e , L n x + G h , e , L n y + H h , e , L n σ and F h , e , L * are the norm flux term and the numerical flux term at the N f p facial interpolation points, respectively; S 0 , h , e , L is the bottom topography source term at the N p interpolation points; and e , L , S e , L = S x , e , L , S y , e , L , S σ , e , L , and e , L f are the mass matrix, stiff matrix, and edge-mass matrix, respectively, and can be expressed as
[ e , L ] i , j = Ω   e , L l i l j d x , [ e , L f ] i , f = Ω   e , L l i l f d x , [ S x , e , L ] i , j = Ω   e , L l i l j x d x ,   [ S y , e , L ] i , j = Ω   e , L l i l j y d x , [ S σ , e , L ] i , j = Ω   e , L l i l j σ d x .

3.2.2. Horizontal Eddy Viscosity Terms

For the discretization of the horizontal eddy viscosity terms in Equation (12), the local discontinuous Galerkin method by Cockburn and Shu [36] is adopted. To discretize the horizontal eddy viscosity terms in the momentum equation about D u , we first express it in the following form:
r = x y K h D x y u ,
where r is an auxiliary function. Following the philosophy of the local discontinuous Galerkin method, we introduce the auxiliary variables and cast Equation (30) into the following form:
q ^ = q ^ x , q ^ y = x y u , q = υ h q ^ , r = x y q ,
with υ h = K h D and x y = / x , / y . For the numerical discretization of Equation (31), we assume a local approximate vector solution u h , q ^ h = u h , q ^ h , x , q ^ h , y V h , 3 d × V h , 3 d 2 and q h = q h , x , q h , y V h , 3 d 2 such that for all ϕ , π V h , 3 d × V h , 3 d 2 , we have
Ω   e , L π q ^ h x y u h d x = Ω   e , L π n h u h u h * d x , Ω   e , L π q h d x = Ω   e , L υ h q ^ h π d x , Ω   e , L ϕ r h d x = Ω   e , L ϕ x y q h d x Ω   e , L ϕ n h q h q h * d x ,
where the terms marked with superscript “*” are the flux terms and are given as
u h * = u h , q h * = υ h x y u h τ u h h ,
where τ is the penalty parameter defined as
τ = N + 1 N + 3 3 n 0 2 A V max υ h , υ h + ,
where N = max N h , N v , n 0 is the number of neighbors of the element, i.e., five for triangular prisms and six for quadrilateral prisms, A is the area of the interface, and V is the average volume of the two adjacent elements.
We further suppose q ^ h , q h , υ h , r h , u h V h , 3 d 2 × V h , 3 d 2 × V h , 3 d × V h , 3 d × V h , 3 d and express these terms as
q ^ h = i = 1 N p q ^ h , e , L , i l i x ,   q h = i = 1 N p q h , e , L , i l i x ,   υ h = i = 1 N p υ h , e , L , i l i x ,   r h = i = 1 N p r h , e , L , i l i x , u h = i = 1 N p u h , e , L , i l i x ,
where q ^ h , e , L , i = q ^ h , e , L , x , i , q ^ h , e , L , y , i and q h , i = q h , e , L , x , i , q h , e , L , y , i are the auxiliary variables at the ith interpolation point. Substituting Equation (35) back into Equation (32), we have
e , L q ^ h , e , L , x = S x , e , L u h , e , L e , L , f u h , e , L u h , e , L * n x , e , L q ^ h , e , L , y = S y , e , L u h , e , L e , L , f u h , e , L u h , e , L * n y ,
q h , e , L , x = υ h , e , L q ^ h , e , L , x ,   q h , e , L , y = υ h , e , L   q   ^ h , e , L , y ,
e , L r h = S x , e , L q h , e , L , x + S y , e , L q h , e , L , y e , L , f q h , e , L q h , e , L *   n h ,
where   n h = n x , n y is the vector of the horizontal components of the outward unit normal outward at the facial interpolation points.
The right-hand side of Equation (38) corresponds to the numerical discretization of the second-order term for the horizontal momentum equation about D u , and the same procedure applies to that about D v .

3.2.3. Vertical Eddy Viscosity and Coriolis Acceleration Terms

For the discretization of the vertical eddy viscosity and Coriolis acceleration terms in Equation (12), we first express it in the following form:
r = σ K m D 2 D u σ f D v ,   i n   Ω K m D 2 D u σ n σ = g N D u ,   o n   Ω N D u = g D D u ,   o n   Ω D
where g N D u is the Neumann boundary condition about D u on Ω N and is given in Equations (17) and (18); g D D u is the Dirichlet boundary condition about D u on Ω D . The symmetric interior penalty discontinuous Galerkin (SIPG) method [37] is adopted, and the corresponding primal form is given as
Ω   h , 3 d υ v ϕ σ D u h σ d x ε I D u h σ υ v ϕ σ d x ε I υ v D u h σ τ D u h σ ϕ σ d x Ω D ϕ υ v D u h σ n σ d x Ω D D u h υ v ϕ σ n σ d x + Ω D τ ϕ n σ 2 D u h d x + Ω   h , 3 d ϕ f D v h d x = Ω   h , 3 d ϕ r d x + Ω D τ ϕ n σ 2 g D D u d x Ω D g D D u υ v ϕ σ n σ d x + Ω N ϕ n σ 2 g N D u d x .
where υ v = K m / D 2 is the vertical eddy viscosity coefficient in the computational domain. The primal form about D v is similar to that about D u . In this study, the Neumann boundary conditions for D u and D v given in Equation (18) are treated implicitly, and the bottom velocity module u c 2 + v c 2 is linearized by using the old values to simplify the calculation. As Equation (40) is only solved in the vertical direction, this equation is independent for each column of prisms and can be solved independently.

3.3. Numerical Discretization of the Primitive Continuity Equation

The water depth is calculated according to Equation (19), and the semi-discrete form is given as
e d D h , e d t = S x , e DU h , e S y , e DV h , e + e f DU h , e DU h , e * n x + e f DV h , e DV h , e * n y ,
where e , S e = S x , e , S y , e and e f are the mass matrix, stiff matrix and edge-mass matrix defined over the two-dimensional computational element Ω e , respectively. Their definitions are analogous to the three-dimensional counterpart defined in Equation (29). D h , e , DU h , e , and DV h , e are the nodal values of the water depth and the depth-averaged velocities in the x and y directions. DU h , e n x + DV h , e n y and DU h , e * n x + DV h , e * n y are the normal flux term and the numerical flux term defined over the boundary of the two-dimensional computational element Ω e and are obtained through the depth integration of their three-dimensional counterpart to preserve the local mass conservation properties.

3.4. Calculation of Vertical Velocity

According to Equations (11) and (19), the governing equation for vertical velocity ω is given as
ω σ = D U D u x + D V D v y .
For the space discretization of Equation (42), suppose that for all ϕ V h , 3 d and Ω e , L Ω h , 3 d , we have
Ω   e , L ϕ ω h σ d x + Ω   e , L ϕ D u h D U h x + D v h D V h y d x = Ω   e , L L a t ϕ n x D u h D U h D u h D U h * d x + Ω   e , L L a t ϕ n y D v h D V h D v h D V h * d x + Ω   e , L T o p ϕ n σ ω h ω h * d x + Ω   e , L B o t ϕ n σ ω h ω h * d x ,
where a h is the approximate solution to a and belongs to V h , 3 d . In Equation (43), the boundary of Ω e , L is split into three parts, i.e., the lateral boundary Ω e , L L a t , the top boundary Ω e , L T o p , and the bottom boundary Ω e , L B o t , and terms marked with an asterisk are the numerical flux terms. Here, the numerical flux n x D u h * + n y D v h * defined on Ω e , L L a t is calculated according to the HLLC Riemann solver, and n x D U h * + n y D V h * is the depth-averaged counterpart. At the bottom boundary, Ω e , L B o t , ω h * is taken as follows:
ω h * = ω h , e , L 1 T o p ,   L > 1 , 0 ,   L = 1 .
The numerical flux is taken to be the boundary condition at the bottom in the case of the bottommost element or using the solution from the element below. Thus, the integral defined over Ω e , L T o p diminishes. With such a definition of numerical flux about ω h * given in Equation (44), we can obtain the vertical velocity ω h layer by layer in each water column starting at the bottom and using the solution from the element below to provide an initial condition.

3.5. Time Stepping

In the discretization presented above, the ordinary differential equations Equation (28) and Equation (41) are obtained for the momentum equations and the primitive continuity equation, and they can be formulated as
d y d t = f exp l ( y ( t ) , t ) + f i m p l ( y ( t ) , t ) ,
where y ( t ) denotes the vector of all discrete degrees of freedom of a step and f exp l y ( t ) , t and f i m p l y ( t ) , t denote the terms treated explicitly and implicitly. For the momentum equations, the explicit terms consist of nonlinear advection, horizontal eddy viscosity, and bottom topography source terms, while the implicit terms consist of the vertical eddy viscosity and Coriolis acceleration terms. For the primitive continuity equation, all terms are treated explicitly.
In this study, the implicit-explicit Runge–Kutta method is adapted to integrate system (45) in time. For an s-stage implicit Runge–Kutta method defined by coefficients a i , j , b i , and c i and a σ = s + 1 stage explicit Runge–Kutta method defined by coefficients a ^ i , j , b ^ i , and c ^ i , this time-stepping method is formulated as
y i = y n + Δ t j = 1 i a ^ i + 1 , j K ^ j + Δ t j = 1 s a i , j K j ,   i 1 , s , y n + 1 = y n + Δ t j = 1 σ b ^ j K ^ j + Δ t j = 1 s b j K j ,
where y n are the known values at time t n , Δ t is the time step, K ^ 1 = f e x p l y n , t n , K ^ i + 1 = f e x p l t n + c ^ i + 1 Δ t , y i , K i = f i m p l y i , t n + c i Δ t , and y n + 1 are the unknown values at t n + 1 = t n + Δ t . The time-stepping method adopted is derived by Ascher [38], and the corresponding coefficients are
c i a i , j b j = γ γ 0 1 1 γ γ 1 γ γ ,
c ^ i a ^ i , j b ^ j = 0 0 0 0 γ γ 0 0 1 δ 1 δ 0 δ 1 δ 0 ,
with γ = 2 2 / 2 and δ = 1 1 / 2 γ . The global time-stepping algorithm from t n to t n + 1 = t n + Δ t is as follows:
  • Calculate the explicit terms K 1 e x p l and K D , 1 e x p l for the three-dimensional momentum equations and primitive continuity equation according to the variables at t n .
  • Calculate the intermediate water depth as
    D i = D n + Δ t j = 1 i a ^ i + 1 , j K D , j e x p l , i 1 , 2
    Likewise, the intermediate conservative variables U i and K i i m p l are
    U i = U n + Δ t j = 1 i a ^ i + 1 , j K j e x p l + Δ t j = 1 i a i , j K j i m p l . i 1 , 2
    With U i available, the intermediate depth-averaged momenta D U i and D V i are calculated through the depth-integration of U i , followed by the calculation of the intermediate vertical velocity ω i . Later, the explicit terms K i + 1 e x p l and K D , i + 1 e x p l for the three-dimensional momentum equations and primitive continuity equation, are calculated according to the intermediate variables.
  • Calculate the water depth D n + 1 and the conservative variables U n + 1 at t n + 1 as
    U n + 1 = U n + Δ t j = 1 3 b ^ j K j e x p l + Δ t j = 1 2 b j K j i m p l ,   D n + 1 = D n + Δ t j = 1 3 b ^ j K D , j e x p l .
  • Integrate U n + 1 along the water depth to obtain the final depth-averaged momenta D U n + 1 and D V n + 1 , followed by the calculation of the final vertical velocity ω n + 1 .

4. Numerical Tests

In this section, several numerical experiments are conducted to verify the performance of the developed model. For all the experiments, the acceleration due to gravity, g , is set to 9.81 m/s2. In addition, all the experiments are run on an Intel Xeon E5-2620 processor with 16 GB of internal memory. Our program is parallelized using OpenMP to run on six cores.

4.1. Manufactured Solution

This test is used to verify the convergence property of the developed model. As the 3D SWEs are a non-linear system, designing a nontrivial function satisfying the equations and the boundary conditions at the same time is not a trivial task. Following the philosophy presented by Salari and Knupp [39], we manually give the solution that satisfies the continuity equation and solve a Dirichlet problem for a modified system that includes a forcing term.
For this test, the Coriolis acceleration terms are not included, the horizontal eddy viscosity coefficient K h is taken as 0.1 m2 s−1, the vertical eddy viscosity coefficient K m is taken as 0.01 m2 s−1, and the analytical solution is given as follows:
b x , y , t = 2 0.005 x + y , η x , y , t = 0.01 sin 0.01 x + t + sin 0.01 y + t , u x , y , σ , t = 0.1 D 1 + σ 3 sin 0.01 x + t , v x , y , σ , t = 0.1 D 1 + σ 3 sin 0.01 x + t .
The vertical velocity at an arbitrary location is obtained through the integration of the continuity equation, Equation (11), from the bottom to the studied location. The horizontal projection of the horizontal domain is Ω 2 d = 0 , 100 × 0 , 100 , and the whole simulation lasts for 86.4 s.
Tests are run with successively refined meshes. For the coarsest mesh, there is one element in the horizontal direction and one layer in the vertical direction. We refine the mesh up to four times by partitioning each quadrilateral prism into eight and compute the L 2 error of the numerical results against the analytical solutions as
L 2 = 1 e Ω e d x e Ω e Q e S x Q e e x a c t x 2 d x 1 / 2 ,
where Q e S x is the numerical result and Q e e x a c t x is the analytical solution. To measure the convergence properties of the model, we further define the convergence rate (CR) as follows
C R = log 2 L 2 m / L 2 m 1 / log 2 L m / L m 1 ,
where L 2 m and L 2 m 1 are the numerical errors on two successive refined meshes and L m and L m 1 are the characteristic lengths of the meshes, which are set to be the longest edge length of the quadrilateral prisms. Table 3 and Table 4 give the L 2 error for each variable and the corresponding CR for N h = N v = 1 and 2, respectively. N e denotes the number of elements in the horizontal direction and N L the number of vertical layers. We note that for water depth D , horizontal momenta D u and D v converge at the optimal rate, i.e., CR = 2 for the former case and 3 for the latter, whereas the vertical velocity ω converges at the suboptimal rate for both cases. We owe this to the fact that the vertical velocity is computed from the numerical results for D u , D v , and their depth-averaged counterparts using the continuity Equation (42), and the CR for vertical velocity would be influenced as the horizontal momenta are not exact, as pointed out by Dawson and Aizinger [25].
Figure 2 shows the relation between the L 2 error and DOFs and the normalized CPU time. Figure 2a indicates that the developed model is convergent, as the L 2 error for each variable diminishes in both cases as the computational mesh is refined. Figure 2b indicates that we can achieve a given error tolerance with fewer DOFs on a coarser mesh when a higher approximation order is adapted. For example, when the error tolerance about water depth D is taken as 1.0 × 10 6 , the normalized CPU time is about 0.033 when the approximation order is N h = N v = 2 , whereas it is about 0.18 when the approximation order is N h = N v = 1 .

4.2. Tide Wave Propagation in a Semi-Closed Bay

This test is used to indicate the accuracy of the developed model for the simulation of both the horizontal and vertical velocity fields. Figure 3 shows the horizontal locations of Gauge points A (2500, 5000), B (52,500, 5000), and C (92,500, 5000) and a sketch of the horizontal domain. This domain is 100,000 m long and 10,000 m wide and is partitioned by a uniform quadrilateral element with a size of 1000 m in both directions. The still water depth is h = 12 m. A cosine tidal wave with amplitude A = 0.25 m and period T = 12 h is imposed at the left boundary, i.e., x = 0, whereas the other boundaries are all treated as slip walls. Both the eddy viscosity terms and the Coriolis acceleration terms are not included.
According to Neuman and Pierson [40], the analytical solution of the surface elevation η , the horizontal velocities u and v , and the vertical velocity w at a point P(x, y, z) at time t are given as
η x , t = A cos ω g h L x / cos ω g h L cos ω t , u x , t = A g h h / cos ω g h L sin ω g h L x sin ω t , v x , t = 0 , w z , t = h + z η + h A ω sin ω t cos ω g h L x / cos ω g h L h + z η + h A 2 ω h 1 cos 2 ω L / g h cos ω t sin ω t .
where ω = 2 π / T is the frequency. For the initial condition, all the velocity components are set to zero, and the surface elevation is set according to the analytical solution. The number of vertical layers is set to 12, and we run this simulation for 24 h.
Figure 4 gives the comparison between the numerical result and the analytical solution for surface elevation η , horizontal velocity u , and vertical velocity w at Gauges A, B, and C, and the vertical coordinate for the three Gauges is z   = −2.5 m. It can be seen that the numerical solution is generally in good agreement with the analytical solution, indicating that the three-dimensional model developed in this study has good accuracy in the simulation of water level and velocity.

4.3. Wind-Induced Water Circulation

In this test, the wind-driven circulation in a rectangular closed basin is simulated. Following Liu et al. [41], we run this case to test the accuracy of the proposed numerical model in predicting the vertical stratified circulation.
With a homogeneous Dirichlet boundary condition adapted at the bottom boundary, Koutitas et al. [42] derived the analytical solution of the horizontal velocity in the center of the basin as
u σ = τ s x σ + 1 h 3 σ + 1 4 K v ρ 0 .   σ 1 , 0
where h is the still water depth of the basin. Likewise, Huang [43] also derived that when the Neumann boundary condition is adapted, the result reads
u σ = 1 6 K v g η x h 2 3 σ 2 1 + τ s x h 2 ρ 0 K v 2 σ + 1 ,   σ 1 , 0 η x = 3 2 τ s x ρ 0 g h 2 K v + C f u c h 3 K v + C f u c h ,
where u c is the velocity magnitude at the center of the bottommost element, and C f u c = 0.005   m / s in this study. In this part, these two situations are simulated.
This case is set as follows: the basin is 2000 m long and 800 m wide, h = 10   m , and τ s x = 1.5   P a . The Coriolis acceleration and horizontal eddy viscosity terms are ignored, and the vertical eddy viscosity coefficient K m is set to be 0.001 m2 s−1; the whole simulation lasts for 3600 s.
Figure 5a,b shows the comparison between the numerical results and the analytical solution for the vertical distribution of horizontal velocity at the center of the basin. The comparison indicates that the numerical results match well with the analytical solution for both cases.
We also run the simulations with different numbers of vertical layers ( N L ). To quantitatively measure the difference between different simulations, we define the root mean square error (RMSE) as
RMSE = j = 1 N T u a , j u j 2 N T ,
where N T is the number of simulated results, u j is the numerical result at point j , and u a , j is the analytical solution at the same point. Table 5 gives the RMSE for different simulations. We find that the more vertical layers there are, the better the result will be for each case.

4.4. Generation of the Ekman Profile

Following Trahan et al. [44], we run this test to check whether the developed model can generate the Ekman layer, where the Coriolis acceleration terms, the pressure gradient, and the turbulent drag are in balance. The analytic solution for the Ekman velocity profile was developed by Ekman [45] and is expressed as
u = u 0 e z / D f cos π 4 z D f , v = v 0 e z / D f sin π 4 z D f , v 0 = τ s x ρ K v f 0.5 , D f = 2 K v f 0.5 , f = 2 ϖ sin θ .
where D f is the folding depth, K v = 0.1 m2/s, τ s x = 0.1 Pa, and θ is the geophysical latitude ( 45 ° ). Theoretically, the flow direction at the surface would rotate counterclockwise 45 degrees to that of the wind stress.
The computational domain is 40,000 km long and 40,000 km wide and is partitioned by a uniform quadrilateral element with a size of 2000 km in both directions. A constant depth of 200 m is set, and 50 vertical layers are used. The bottom drag coefficient C f is set to 0.0025. We run this simulation for 200 days. Figure 6 gives the comparison of the analytic solution and the computed Ekman profile. The numerical results accurately match the analytic solution and individual velocity components with an average error of less than 0.001 m/s.

4.5. Tidal Flow in Bohai Bay

For the last case, we simulate tidal flow in the western region of Bohai Bay to test the ability of the developed model for the simulation of actual tidal wave motion problems. The Bohai Bay is one of the three major bays of the Bohai Sea in China. It is located in the west of the Bohai Sea and connected with Hebei and Shandong Province and Tianjin City. It’s a typical semi-closed sea area. The average tidal range in Bohai Bay is 2.5 m, and its tidal current is irregularly semi-diurnal. In recent years, a lot of port and reclamation projects in Bohai Bay have been carried out, which has led to a large number of large-scale engineering structures. The accurate simulation of the tidal current field is of great significance for the convection–diffusion process of temperature, salt, pollutants, etc., as well as the sediment transport in Bohai Bay.
The domain geometry for this case is shown in Figure 7a; the still water depth ranges from 38 m at the east open boundary to 2 m at the coast. There are five measuring gauges in the computational domain, and the distribution of these points is shown in Figure 7b. The computational domain is partitioned by 14,627 triangle grids, and the edge length of the finite element mesh ranges from 2 km to 100 m. Five equidistant layers are used in the vertical direction. Wetting/drying is not accounted in the present study.
The tidal elevations at the east open boundaries were provided by Chinatide [46], which is a tidal forecast system through nine harmonic constants of Q1, P1, O1, K1, N2, M2, S2, K2, and Sa. The Coriolis parameter was taken as 9.1557 × 10 5 s 1 , and the horizontal eddy viscosity coefficient was set to 70 m2 s−1. The vertical eddy viscosity coefficient was modelled using the k ϵ turbulence closures provided by GOTM [47], and the coupling between GOTM and the hydrodynamic model is established following the same philosophy presented by Tuomas et al. [32]. The bottom roughness length is set as 0.1 mm. Field measurements from 12/08/2018 to 12/09/2018 are used to verify the developed model.
Figure 8 shows the flow field at the surface and the middle and bottom layers at the time of flood tide (left) and ebb tide (right), respectively. As we can see, there is little difference between the velocity fields of the surface and the middle layers. Because of the bottom friction, the velocity of the bottom layer is smaller than that of the surface and the middle layers, which is about 60% of the magnitude of the latter two. The results of the three-dimensional velocity field show that the model can reflect the three-dimensional characteristics of tidal flow.
Figure 9 shows the comparison between the simulated surfaced elevation and the field measurements at Gauges P1 and P2. The simulated data are quantitatively very close to the field measurements, and the tidal phasing is virtually identical. The left column of Figure 10 gives the comparison between the simulated velocity and the field measurements at gauges P3, P4, and P5. Small discrepancies exist between the simulated velocity and the field measurements. This is especially the case for P5, where the maximum simulated velocity is about 15% smaller than that of the field measurements. It can be attributed to that the computational mesh used in this simulation is relatively coarse, and the wetting/drying is not considered. Although differences exist, the simulated data are in phase with the field measurements as the comparison between the simulated flow direction and the field measurements shows in the right column of Figure 10. Generally, the developed hydrodynamic model can be used for the simulation of tidal wave motion in realistic situations.

5. Conclusions

In this study, a three-dimensional hydrodynamic model based on the Discontinuous Galerkin method is developed. The σ -coordinate system is used in the vertical direction to circumvent the mismatch of lateral faces due to the discontinuous nature of the solution between different elements. Non-linear advective terms are included, and the model is marched in time with the implicit-explicit Runge–Kutta method. The vertical eddy viscosity and the Coriolis acceleration terms are discretized with the symmetric interior penalty discontinuous Galerkin method and are treated implicitly, while others are treated explicitly.
The developed model is validated with a series of numerical tests. The case with the manufactured solution indicates that the model is convergent, and the L 2 error about the variables diminishes as we increase the approximation order or refine the computational mesh. Calculating the CR indicates that water depth D and horizontal momenta D u and D v can all converge at the optimal order, whereas the vertical velocity converges in a suboptimal order. This case also indicates that we can obtain a better numerical result on a computational mesh with fewer DOFs when a higher approximation is adapted. Cases of tide-induced three-dimensional flow in a semi-closed bay, wind-induced water circulation, and the generation of the Ekman profile indicate that the developed model can well simulate the velocity field in both horizontal and vertical directions, can reasonably simulate the vertical stratified circulation, and can simulate water motion under the combined effect of wind stress and Coriolis acceleration terms. For the Bohai Bay case, the simulated surface elevation, velocity, and flow direction match reasonably well with the field measurements, indicating the model can be used for the simulation of tidal wave motion in realistic situations.
More applications of the developed model and development of the wetting/drying treatment and the baroclinic model will be presented in forthcoming work.

Author Contributions

Formal analysis, G.R.; Funding acquisition, Q.Z.; Methodology, G.R.; Software, G.R. and Z.C.; Supervision, Q.Z.; Validation, G.R.; Visualization, G.R.; Writing—original draft, G.R.; Writing—review and editing, Q.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This study was financially supported by the Joint Funds of the National Natural Science Foundation of China (Grant No. U1906231), the Natural Science Foundation of Tianjin (Grant No. 19JCZDJC40200).

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lee, J.; Lee, J.; Yun, S.-L.; Kim, S.-K. Three-Dimensional Unstructured Grid Finite-Volume Model for Coastal and Estuarine Circulation and Its Application. Water 2020, 12, 2752. [Google Scholar] [CrossRef]
  2. Blaise, S.; Comblen, R.; Legat, V.; Remacle, J.-F.; Deleersnijder, E.; Lambrechts, J. A discontinuous finite element baroclinic marine model on unstructured prismatic meshes. Ocean. Dyn. 2010, 60, 1371–1393. [Google Scholar] [CrossRef] [Green Version]
  3. Li, L.; Zhang, Q. Development of an efficient wetting and drying treatment for shallow-water modeling using the quadrature-free Runge-Kutta discontinuous Galerkin method. Int. J. Numer. Methods Fluids 2021, 93, 314–338. [Google Scholar] [CrossRef]
  4. Ran, G.; Zhang, Q.; Li, L. A depth-integrated non-hydrostatic model for nearshore wave modelling based on the discontinuous Galerkin method. Ocean. Eng. 2021, 232, 108661. [Google Scholar] [CrossRef]
  5. Di Pietro, D.A.; Marche, F. Weighted interior penalty discretization of fully nonlinear and weakly dispersive free surface shallow water flows. J. Comput. Phys. 2018, 355, 285–309. [Google Scholar] [CrossRef] [Green Version]
  6. Shchepetkin, A.F.; McWilliams, J.C. The regional oceanic modeling system (ROMS): A split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean. Model. 2005, 9, 347–404. [Google Scholar] [CrossRef]
  7. Mellor, G. Users Guide for a Three-Dimensional, Primitive Equation, Numerical Ocean Model; Princeton University: Princeton, NJ, USA, 2004. [Google Scholar]
  8. Griffies, S.M. Elements of the Modular Ocean Model (MOM); GFDL Ocean Group Technical Report; NOAA: Princeton, NJ, USA, 2012. [Google Scholar]
  9. Burchard, H.; Bolding, K. GETM: A General Estuarine Transport Model. Scientific Documentation; European Commission Technical Report EUR 20253 EN; European Commission: Brussels, Belgium, 2002; p. 157. [Google Scholar]
  10. Casulli, V.; Cheng, R.T. Semi-implicit finite difference methods for three-dimensional shallow water flow. Int. J. Numer. Methods Fluids 1992, 15, 629–648. [Google Scholar] [CrossRef]
  11. Deltares Systems. User Manual of Delft3D-FLOW—Simulation of Multi-Dimensional Hydrodynamic Flows and Transport Phenomena, Including Sediments; Report of Delft Hydraulics; Deltares: Delft, The Netherlands, 2003; p. 497. [Google Scholar]
  12. Hamrick, J.; Wu, T. Computational Design and Optimization of the EFDC/HEM3D Surface Water Hydrodynamic and Eutrophication. In Next Generation Environmental Models and Computational Methods; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 1997; pp. 143–161. [Google Scholar]
  13. Blumberg, A.F.; Mellor, G.L. A Description of a Three-Dimensional Coastal Ocean Circulation Model; AGU: Washington, DC, USA, 1987; pp. 1–16. [Google Scholar]
  14. Luyten, P.J.; Jones, J.E.; Proctor, R. A numerical study of the long-and short-term temperature variability and thermal circulation in the North Sea. J. Phys. Oceanogr. 2003, 33, 37–56. [Google Scholar] [CrossRef]
  15. Alistair, A.; Jean-Michel, C.; Stephanie, D.; Constantinos, E.; David, F.; Gael, F.; Baylor, F.-K.; Patrick, H.; Chris, H.; Ed, H. Mitgcm User Manual; MIT: Cambridge, MA, USA, 2018. [Google Scholar]
  16. Martins, F.; Leitão, P.; Silva, A.; Neves, R. 3D modelling in the Sado estuary using a new generic vertical discretization approach. Oceanol. Acta 2001, 24, 51–62. [Google Scholar] [CrossRef] [Green Version]
  17. Chen, C.; Beardsley, R.; Cowles, G.; Qi, J.; Lai, Z.; Gao, G.; Stuebe, D.; Xu, Q.; Xue, P.; Ge, J. An Unstructured-Grid, Finite-Volume Community Ocean Model: FVCOM User Manual; SMAST/UMASSD Technical Report-13-0701; MIT: Cambridge, MA, USA, 2012; p. 404. [Google Scholar]
  18. Casulli, V.; Walters, R.A. An unstructured grid, three-dimensional model based on the shallow water equations. Int. J. Numer. Methods Fluids 2000, 32, 331–348. [Google Scholar] [CrossRef]
  19. Nan, Z.; Sheng, J.; Congfang, A.; Weiye, D. An integrated water-conveyance system based on Web GIS. J. Hydroinform. 2018, 20, 668–686. [Google Scholar] [CrossRef] [Green Version]
  20. Wang, Q. The Finite Element Ocean Model and Its Aspect of Vertical Discretization. Ph.D. Thesis, Universität Bremen, Bremen, Germany, 2007. [Google Scholar]
  21. Zhang, Y.; Baptista, A.M. SELFE: A semi-implicit Eulerian–Lagrangian finite-element model for cross-scale ocean circulation. Ocean. Model. 2008, 21, 71–96. [Google Scholar] [CrossRef]
  22. Luettich, R.A.; Westerink, J.J.; Scheffner, N.W. ADCIRC: An Advanced Three-Dimensional Circulation Model for Shelves, Coasts, and Estuaries. Report 1, Theory and methodology of ADCIRC-2DD1 and ADCIRC-3DL; Department of the Army, US Army Corps of Engineers Waterways Experiment Station, Coastal Engineering Research Center: Vicksburg, MS, USA, 1992. [Google Scholar]
  23. Zhang, Y.J.; Ye, F.; Stanev, E.V.; Grashorn, S. Seamless cross-scale modeling with SCHISM. Ocean. Model. 2016, 102, 64–81. [Google Scholar] [CrossRef] [Green Version]
  24. Hervouet, J.-M. Hydrodynamics of Free Surface Flows: Modelling with the Finite Element Method; John Wiley & Sons: Hoboken, NJ, USA, 2007. [Google Scholar]
  25. Dawson, C.; Aizinger, V. A discontinuous Galerkin method for three-dimensional shallow water equations. J. Sci. Comput. 2005, 22, 245–267. [Google Scholar] [CrossRef]
  26. Aizinger, V.; Dawson, C. The local discontinuous Galerkin method for three-dimensional shallow water flow. Comput. Methods Appl. Mech. Eng. 2007, 196, 734–746. [Google Scholar] [CrossRef]
  27. Aizinger, V.; Proft, J.; Dawson, C.; Pothina, D.; Negusse, S. A three-dimensional discontinuous Galerkin model applied to the baroclinic simulation of Corpus Christi Bay. Ocean. Dyn. 2013, 63, 89–113. [Google Scholar] [CrossRef]
  28. Blaise, S. Development of a Finite Element Marine Model. Ph.D. Thesis, UCL-Université Catholique de Louvain, Ottignies-Louvain-la-Neuve, Belgium, 2009. [Google Scholar]
  29. Comblen, R. Discontinuous Finite-Element Methods for Two-and Three-Dimensional Marine Flows. Ph.D. Thesis, UCL-Université Catholique de Louvain, Ottignies-Louvain-la-Neuve, Belgium, 2010. [Google Scholar]
  30. Gandham, R. High Performance High-Order Numerical Methods: Applications in Ocean Modeling. Ph.D. Thesis, Rice University, Houston, TX, USA, 2015. [Google Scholar]
  31. Conroy, C.J.; Kubatko, E.J. Hp discontinuous Galerkin methods for the vertical extent of the water column in coastal settings part I: Barotropic forcing. J. Comput. Phys. 2016, 305, 1147–1171. [Google Scholar] [CrossRef]
  32. Kärnä, T.; Legat, V.; Deleersnijder, E.; Burchard, H. Coupling of a discontinuous Galerkin finite element marine model with a finite difference turbulence closure model. Ocean. Model. 2012, 47, 55–64. [Google Scholar] [CrossRef]
  33. Phillips, N.A. A coordinate system having some special advantages for numerical forecasting. J. Atmos. Sci. 1957, 14, 184–195. [Google Scholar] [CrossRef]
  34. Toro, E.F. Shock-Capturing Methods for Free-Surface Shallow Flows; John Wiley and Sons: Hoboken, NJ, USA, 2001. [Google Scholar]
  35. Hesthaven, J.S.; Warburton, T. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2007. [Google Scholar]
  36. Cockburn, B.; Shu, C.-W. The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM J. Numer. Anal. 1998, 35, 2440–2463. [Google Scholar] [CrossRef]
  37. Di Pietro, D.A.; Ern, A. Mathematical Aspects of Discontinuous Galerkin Methods; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2011; Volume 69. [Google Scholar]
  38. Ascher, U.M.; Ruuth, S.J.; Spiteri, R.J. Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Appl. Numer. Math. 1997, 25, 151–167. [Google Scholar] [CrossRef]
  39. Salari, K.; Knupp, P. Code Verification by the Method of Manufactured Solutions; SAND2000-1444; Sandia National Laboratories: Albuquerque, NM, USA, 2000. [Google Scholar]
  40. Neumann, G.; Pierson, W.J., Jr. Principles of Physical Oceanography; Prentice-Hall: Englewood Cliffs, NJ, USA, 1966. [Google Scholar]
  41. Liu, X.; Mohammadian, A.; Sedano, J.Á.I.; Kurganov, A. Three-dimensional shallow water system: A relaxation approach. J. Comput. Phys. 2017, 333, 160–179. [Google Scholar] [CrossRef] [Green Version]
  42. Koutitas, C.; O’Connor, B. Modeling three-dimensional wind-induced flows. J. Hydraul. Div. ASCE 1980, 106, 1843–1865. [Google Scholar] [CrossRef]
  43. Huang, W. Three Dimensional Numerical Modeling Of Circulation And Water Quality Induced By Combined Sewage Overflow Discharges. Ph.D. Thesis, Department of Ocean Engineering, University of Rhode Island, Kingston, RI, USA, 1993. [Google Scholar]
  44. Trahan, C.; Savant, G.; Berger, R.; Farthing, M.; McAlpin, T.; Pettey, L.; Choudhary, G.; Dawson, C.N. Formulation and application of the adaptive hydraulics three-dimensional shallow water and transport models. J. Comput. Phys. 2018, 374, 47–90. [Google Scholar] [CrossRef]
  45. Price, J.F.; Weller, R.A.; Schudlich, R.R. Wind-driven ocean currents and Ekman transport. Science 1987, 238, 1534–1538. [Google Scholar] [CrossRef]
  46. Li, M.; Zheng, J. Introduction to Chinatide software for tide prediction in China seas. J. Waterw. Harb. 2007, 28, 65–68. (In Chinese) [Google Scholar]
  47. Burchard, H.; Bolding, K.; Villarreal, M.R. GOTM, a General Ocean Turbulence Model: Theory, Implementation and Test Cases; Technical Report EUR 18745 EN; European Commission: Brussels, Belgium, 1999.
Figure 1. Schematic view of a column of three-dimensional computational elements. (a) quadrangular prisms; (b) triangular prisms.
Figure 1. Schematic view of a column of three-dimensional computational elements. (a) quadrangular prisms; (b) triangular prisms.
Water 15 00156 g001
Figure 2. L2 error versus DOFs and normalized CPU time for water depth D, vertical velocity ω, and horizontal momenta Du and Dv at different approximation orders (Nh = Nv = 1(black), Nh = Nv = 2(red)) on a sequence of meshes. (a) L2 error versus DOFs; (b) L2 error versus normalized CPU time.
Figure 2. L2 error versus DOFs and normalized CPU time for water depth D, vertical velocity ω, and horizontal momenta Du and Dv at different approximation orders (Nh = Nv = 1(black), Nh = Nv = 2(red)) on a sequence of meshes. (a) L2 error versus DOFs; (b) L2 error versus normalized CPU time.
Water 15 00156 g002
Figure 3. Location of the gauge points and sketch of the two-dimensional computational domain.
Figure 3. Location of the gauge points and sketch of the two-dimensional computational domain.
Water 15 00156 g003
Figure 4. Comparison between the numerical result (black line) and the analytical solution (red circle) for surface elevation η , horizontal velocity u , and vertical velocity w at Gauges (AC).
Figure 4. Comparison between the numerical result (black line) and the analytical solution (red circle) for surface elevation η , horizontal velocity u , and vertical velocity w at Gauges (AC).
Water 15 00156 g004aWater 15 00156 g004b
Figure 5. Comparison between the numerical results and the analytical solution for the vertical distribution of horizontal velocity at the center of the basin. (a) Homogeneous Dirichlet boundary condition; (b) Neumann boundary condition.
Figure 5. Comparison between the numerical results and the analytical solution for the vertical distribution of horizontal velocity at the center of the basin. (a) Homogeneous Dirichlet boundary condition; (b) Neumann boundary condition.
Water 15 00156 g005
Figure 6. Comparison between the numerical results and the analytical solution for individual velocity components and velocity magnitude for the Ekman profile.
Figure 6. Comparison between the numerical results and the analytical solution for individual velocity components and velocity magnitude for the Ekman profile.
Water 15 00156 g006
Figure 7. Computational domain geometry for tidal flow in the west region of Bohai Bay. (a) Computational domain geometry; (b) Gauge points.
Figure 7. Computational domain geometry for tidal flow in the west region of Bohai Bay. (a) Computational domain geometry; (b) Gauge points.
Water 15 00156 g007
Figure 8. Flow field at different levels at the flood tide (left column) and the ebb tide (right column). (a) surface ( σ = 0 ); (b) middle ( σ = 0.5 ); (c) bottom ( σ = 1.0 ).
Figure 8. Flow field at different levels at the flood tide (left column) and the ebb tide (right column). (a) surface ( σ = 0 ); (b) middle ( σ = 0.5 ); (c) bottom ( σ = 1.0 ).
Water 15 00156 g008
Figure 9. Comparison between the numerical results and the measured surface elevation at Gauges P1 and P2 for the Bohai Bay case. (a) P1; (b) P2.
Figure 9. Comparison between the numerical results and the measured surface elevation at Gauges P1 and P2 for the Bohai Bay case. (a) P1; (b) P2.
Water 15 00156 g009
Figure 10. Comparison between the simulated results and the measured data for velocity (left column) and flow direction (right column) at Gauges P3, P4, and P5 for the Bohai Bay case. (a) P3; (b) P4; (c) P5.
Figure 10. Comparison between the simulated results and the measured data for velocity (left column) and flow direction (right column) at Gauges P3, P4, and P5 for the Bohai Bay case. (a) P3; (b) P4; (c) P5.
Water 15 00156 g010
Table 2. Notations for the governing equations of the three-dimensional hydrodynamic model.
Table 2. Notations for the governing equations of the three-dimensional hydrodynamic model.
VariablesExplanation
u x , y , z , t velocity of water in x direction
v x , y , z , t velocity of water in y direction
w x , y , z , t velocity of water in z direction
gacceleration due to gravity
θ angle of geographical latitude
ϖ magnitude of the angular velocity of the Earth
f = 2 ϖ sin θ Coriolis parameter
K m vertical eddy viscosity coefficient
K h horizontal eddy viscosity coefficient
ρ 0 density
τ s x wind stress in x direction
τ s y wind stress in y direction
η x , y , t the surface elevation
u η x , y , t surface velocity of water in x direction
v η x , y , t surface velocity of water in y direction
b x , y bottom elevation
u b x , y , t bottom velocity of water in x direction
v b x , y , t bottom velocity of water in y direction
z 0 x , y , t half height of the bottommost element
u c x , y , t velocity   of   water   at   z 0 in x direction
v c x , y , t velocity   of   water   at   z 0 in y direction
C f x , y , t drag coefficient
κ von Karman constant
z 0 b x , y bottom roughness parameter
ω x , y , σ , t vertical velocity in the computational domain
U ¯ x , y , t depth-averaged velocity in x direction
V ¯ x , y , t depth-averaged velocity in y direction
D x , y , t = η x , y , t b x , y water depth
Table 3. N h = N v = 1 , the L2 error for each variable, and the corresponding CR.
Table 3. N h = N v = 1 , the L2 error for each variable, and the corresponding CR.
NeNLEr(D)CR(D)Er(Du)CR(Du)Er(Dv)CR(Dv)Er(ω) CR(ω)
110.00033 0.00169 0.00169 5.224 × 1 0 5
42 6.125 × 1 0 5 2.420.000312.470.000312.47 1.277 × 1 0 5 2.03
164 1.340 × 1 0 5 2.19 6.086 × 1 0 5 2.33 6.086 × 1 0 5 2.33 3.748 × 1 0 6 1.77
648 3.170 × 1 0 6 2.08 1.350 × 1 0 5 2.17 1.350 × 1 0 5 2.17 1.239 × 1 0 6 1.60
25616 7.617 × 1 0 7 2.06 3.083 × 1 0 6 2.13 3.083 × 1 0 6 2.13 4.527 × 1 0 7 1.45
Table 4. N h = N v = 2 , the L2 error of each variable, and the corresponding CR.
Table 4. N h = N v = 2 , the L2 error of each variable, and the corresponding CR.
NeNLEr(D)CR(D)Er(Du)CR(Du)Er(Dv)CR(Dv)Er(ω) CR(ω)
11 9.997 × 1 0 6 0.00015 0.00015 4.95 × 1 0 6
42 1.353 × 1 0 6 3.10 1.804 × 1 0 5 3.10 1.804 × 1 0 5 3.10 5.71 × 1 0 7 3.11
164 1.703 × 1 0 7 3.04 2.187 × 1 0 6 3.04 2.187 × 1 0 6 3.04 1.13 × 1 0 7 2.33
648 2.196 × 1 0 8 3.05 2.649 × 1 0 7 3.05 2.649 × 1 0 7 3.05 3.00 × 1 0 8 1.91
25616 3.276 × 1 0 9 3.02 3.264 × 1 0 8 3.02 3.264 × 1 0 8 3.02 6.99 × 1 0 9 2.10
Table 5. RMSE for simulations with different vertical layers and different boundary conditions.
Table 5. RMSE for simulations with different vertical layers and different boundary conditions.
NL5101520
Homogeneous Dirichlet boundary condition
RMSE0.01070.00510.00410.0037
Neumann boundary condition
RMSE0.02190.01380.01140.0103
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Ran, G.; Zhang, Q.; Chen, Z. Development of a Three-Dimensional Hydrodynamic Model Based on the Discontinuous Galerkin Method. Water 2023, 15, 156. https://doi.org/10.3390/w15010156

AMA Style

Ran G, Zhang Q, Chen Z. Development of a Three-Dimensional Hydrodynamic Model Based on the Discontinuous Galerkin Method. Water. 2023; 15(1):156. https://doi.org/10.3390/w15010156

Chicago/Turabian Style

Ran, Guoquan, Qinghe Zhang, and Zereng Chen. 2023. "Development of a Three-Dimensional Hydrodynamic Model Based on the Discontinuous Galerkin Method" Water 15, no. 1: 156. https://doi.org/10.3390/w15010156

APA Style

Ran, G., Zhang, Q., & Chen, Z. (2023). Development of a Three-Dimensional Hydrodynamic Model Based on the Discontinuous Galerkin Method. Water, 15(1), 156. https://doi.org/10.3390/w15010156

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop