A DISCONTINUOUS GALERKIN METHOD FOR THREE

16 downloads 0 Views 7MB Size Report
three-dimensional shallow water models. The rest of this paper is organized as follows. In the next section, we outline the governing equations and boundary ...
A DISCONTINUOUS GALERKIN METHOD FOR THREE-DIMENSIONAL SHALLOW WATER EQUATIONS CLINT DAWSON† AND VADYM AIZINGER† Abstract. We describe the application of a local discontinuous Galerkin method to the numerical solution of the three-dimensional shallow water equations. The shallow water equations are used to model surface water flows where the hydrostatic pressure assumption is valid. The authors recently developed a DG method for the depth-integrated shallow water equations. The method described here is an extension of these ideas to non depth-integrated models. The method and its implementation are discussed, followed by numerical examples on several test problems. Key words. Shallow water equations, discontinuous Galerkin method, free surface

1. Introduction. In this paper, we describe a local discontinuous Galerkin (LDG) method for the three-dimensional shallow water equations. The shallow water equations are derived from the incompressible Navier-Stokes equations defined on a domain with a moving free (top) surface. The shallow water assumption reduces the vertical momentum equation to the hydrostatic pressure relation ∂p = −ρg, ∂z where p is the hydrostatic pressure, ρ is density and g is gravitational acceleration. In many cases of practical importance, vertical effects are negligible and one can integrate the horizontal momentum equations and the continuity equation over the depth, applying appropriate boundary conditions at the bottom and free surface, to obtain the two-dimensional shallow water equations. Where eddy viscosity and/or density variations are important, for example, where salinity and temperature are spatially varying, then the full three-dimensional model should be employed. See [16, 15] for a discussion of shallow water models in both two and three dimensions. Shallow water equations are complicated by many nonlinear effects and are defined on complex domains involving varying bottom topography and coastal shorelines. Viscosity effects, especially horizontal viscosity, are usually relatively small, and algorithms which are stable and accurate for smooth to highly advective flows on general geometries are of interest for the numerical solution of these models. A substantial literature exists on the application of various finite difference and finite element methods to the three-dimensional shallow water equations; see for example [10, 12, 13]. The search is still on, however, for methods which are locally mass conservative, can handle very general types of elements, and are stable and accurate under highly varying flow regimes. Recently developed algorithms such as the discontinuous Galerkin (DG) method are therefore of great interest within the shallow water community. DG finite element methods are promising because of their ability to handle geometrically complex elements, use of shock-capturing numerical fluxes, adaptivity in polynomial order, and mass conservation properties; see [6] for a historical overview of DG methods. In [2, 3], we investigated DG and related finite volume methods for the solution of the two-dimensional shallow water equations. Viscosity (second-order † Center for Subsurface Modeling - C0200; Institute for Computational Engineering and Sciences (ICES); The University of Texas at Austin; Austin, TX 78712.

1

2

Dawson, Aizinger

derivative) terms are handled in this method through the so-called local discontinuous Galerkin (LDG) framework [4], which involves rewriting the model as a first order system of equations. In this paper, we describe the extension of the LDG method to three-dimensional shallow water models. The rest of this paper is organized as follows. In the next section, we outline the governing equations and boundary conditions for our model. In Section 3, we outline our solution strategy, which is described in detail in Sections 4 and 5. In Section 6, some preliminary numerical results are given for problems with smooth and rough analytical solutions, and a typical tidal flow problem. 2. Mathematical Formulation. The governing equations include momentum equations for the horizontal velocity components, a continuity equation, and boundary and initial conditions. The domain over which these equations are defined has a free surface described by a kinematic boundary condition. The conservative form of the momentum conservation equations can be written as follows: µ ¶ Z ξ ¡ ¢ g 1 ∂uxy v (1) ∇xy ρdz, +∇· uxy uT − D∇uxy +g∇xy ξ = f − ∇xy pa − −u ∂t ρ0 ρ0 z ³ ´ ∂ ∂ where ∇xy = ∂x , ∂y , ξ is the value of the z coordinate at the free surface, u =

(u, v, w)T is the velocity vector, uxy = (u, v)T is the vector of horizontal velocity components, f is the Coriolis coefficient, g is the free fall acceleration, ρ and ρ 0 are the actual and the reference water densities, pa is the atmospheric pressure, and D is a tensor of eddy viscosity coefficients defined as follows: ¶ µ Du 0 (2) D= 0 Dv with Du , Dv being 3x3 symmetric positive-definite matrices. The continuity equation can be written as ∂w ∂u ∂v + + = 0. ∂x ∂y ∂z

(3)

We augment the system with the following boundary conditions: • At the bottom, we have no normal flow (4)

u(zb ) · n = 0

and no slip for horizontal velocity components (5)

u(zb ) = v(zb ) = 0.

where zb is the value of the z coordinate at the sea bed and n = (nx , ny , nz )T is an exterior unit normal to the face. • The free surface boundary conditions have the form: (6)

∂ξ ∂ξ ∂ξ + u(ξ) + v(ξ) − w(ξ) = 0. ∂t ∂x ∂y

(7)

∂v ∂u = = 0. ∂n ∂n

A DG method for the 3D SWE

3

On the lateral boundaries, we can have several types of boundary conditions (note, that we assume all lateral boundaries to be strictly vertical; therefore, if n = (nx , ny , nz )T is an exterior normal to a lateral boundary face then nz = 0): • Land boundary: No normal flow (8)

un = u · n = 0,

and no drag on tangential flow ∂uτ = 0. ∂n

(9)

• Open sea boundary: Zero normal derivative of the horizontal velocity components ∂v ∂u = = 0, ∂n ∂n

(10)

and prescribed surface elevation ξs (x, y, t) (11)

ξ = ξs (x, y, t).

• River boundary: Prescribed normal un (x, y, t) and tangential uτ (x, y, t) velocities (12)

unx + vny = un (x, y, t),

−uny + vnx = uτ (x, y, t),

and prescribed surface elevation ξr (x, y, t) (13)

ξ = ξr (x, y, t).

• Radiation boundary: Zero normal derivative of the horizontal velocity components ∂u ∂v = = 0. ∂n ∂n

(14)

Analytically, the free surface elevation can be computed using (6), but, numerically, a more robust way to do it is to integrate the continuity equation (3) over the depth and, taking into account the boundary conditions at the top (6) and bottom (4), obtain the following 2D equation for the surface elevation commonly called the primitive continuity equation (PCE): (15)

∂ξ ∂ + ∂t ∂x

Z

ξ zb

∂ udz + ∂y

Z

ξ

vdz = 0. zb

The equation above or, alternatively, the wave equation derived by substituting the 2D momentum equations in (15) is a standard part of just about every 2- and 3D shallow-water model in existence today, though using this equation to compute the position of the free surface within a framework of a 3D solver based on the discontinuous Galerkin method is not trivial, unless we resort to solving an auxiliary 2D problem as is done in many 3D shallow water solvers. We will discuss this issue in some detail in section 4.

4

Dawson, Aizinger

Let us denote H = ξ − zb . Then, we can rewrite conservation equations (15), (1) in the following compact form: (16)

∂H + ∇xy · CH = 0, ∂t

∂uxy + ∇ · (CM − D∇uxy ) = F, ∂t  2 u + gH ³R ´T R ξ ξ uv where CH = zb u dz, zb v dz , CM = (Cu , Cv ) =  uw ! Ã R ξ ∂ρ g 1 ∂pa b + f v − − dz g ∂z ∂x ρ0 ∂x ρ0 z ∂x R ξ ∂ρ . F= g 1 ∂pa b − f u − − g ∂z ∂y ρ0 ∂y ρ0 z ∂y dz

(17)

 uv v 2 + gH , vw

3. General issues and solution strategy.

3.1. Solution strategy. The general solution strategy employed in our implementation is not substantially different from ones found in other 3D solvers. The main differences lie in the fact that all state variables are approximated in space with functions which are continuous on each element but can be discontinuous across inter-element boundaries. For time stepping, we use explicit Runge-Kutta methods, choosing the time step and scheme order appropriately to keep the error of the time stepping algorithm well within the error of space discretization. Within a time step (or a time substep in case of multistage Runge-Kutta methods), first, we solve the mass and momentum conservation equations. These equations are tightly coupled and must be dealt with simultaneously. Then, for given values of horizontal velocity components uxy = (u, v)T , we compute vertical velocity w using the discrete continuity equation to obtain a divergence-free velocity field. The latter equation is not time-dependent and is solved as an initial value problem layer by layer starting at the bottom and using the solution from the layer below (or the boundary conditions for u and v at the bottom in case of the bottommost layer) as an initial value. Every few time steps the mesh geometry is updated using computed values of the surface elevation. Frequency of mesh update can be chosen according to the type of problem we are solving. 3.2. Computational mesh. The most common type of mesh employed in 3D finite element simulations of shallow-water flow is a 2D grid projected vertically and subdivided into layers using a Cartesian or σ-stretched coordinate system. This approach agrees well with the physical anisotropy of the problem: The vertically directed gravity force usually is the main body force acting on the system. The grids used in our implementation also belong to this type. In our model, we use prismatic elements with triangular cross-section, strictly vertical lateral faces, and flat but not necessarily parallel top and bottom faces. An important feature of most modern 3D shallow water solvers, the dynamic free surface, turns out to be the main stumbling block on the way to a stable DG implementation. Since the surface elevation is one of our state variables living in a discontinuous approximation space, the grid determined by our discrete solution would also have a discontinuous free surface, thus, the lateral faces shared by neighboring elements in the uppermost layer of the mesh would, in general, not match. This causes the boundary integrals over those faces to be ill-defined. In order to avoid

5

A DG method for the 3D SWE

this problem, we smooth the free surface of our mesh (e.g. by least squares fit) and compute all 3D integrals on a grid with a continuous free surface. With regard to this procedure, it must be noted that the smoothing postprocessor does not degrade the local conservation properties of our scheme because it only slightly modifies the geometry of prisms in the surface layer of the mesh whereas the values of computed state variables remain unaffected. 3.3. Treatment of discontinuities. In order to better understand the issues discussed below we have to keep in mind the fact that we are solving a system of PDEs representing a conservation law using a numerical method in which the approximations to all state variables are discontinuous across inter-element boundaries. This gives rise to a Riemann shock-tube problem at every Gauss integration point on those boundaries. In addition, the mass and momentum conservation equations are tightly coupled; therefore, solutions to these Riemann problems must satisfy the entropy condition for every one of these equations at the same time. Which puts us in a difficult situation: (15) is a 2D equation and is not well-defined on the lateral boundaries of our 3D elements unlike the momentum equations (1). Uncoupling these equations or computing the surface elevation from (6) results in an unstable numerical model because, in this case, the boundary discontinuities are not resolved in an entropy conforming way. Since we don’t want to resort to solving an auxiliary 2D shallowwater problem for the surface elevation, we modify instead the boundary integrals in the discrete mass conservation equation to bring them in a 3D form, thus giving a well defined Riemann problem at every Gauss point on the lateral 3D faces. Here, it must also be noted that the Riemann problems we solve on horizontal faces are significantly less demanding in terms of stability than the ones on lateral faces due to the fact that the surface elevation is continuous across horizontal interelement interfaces. Our experience shows that using upwinding or central differences on those boundary faces does not cause any instability. The more complex Riemann problems on lateral interfaces are handled in our implementation with the help of Roe’s solver that will be discussed in detail in a separate section. 4. Space Discretization, the LDG method. To discretize our problem in space using the LDG method, we first introduce an auxiliary flux variable Q and rewrite the momentum conservation equations in the mixed form: √ ∂uxy + ∇ · (CM + DQ) = F, ∂t √ Q = − D∇uxy .

(18) (19)

Let Th be a general partition of our 3D domain Ω and let Ωe ∈ Th . To obtain a weak formulation, we multiply the equations above by smooth test functions φ, Ψ ; integrate them on each element Ωe ∈ Th ; and, finally, integrate by parts obtaining the following expressions: Z Z √ ∂uxy (CM + DQ) · n φ ds φ dxdydz + (20) ∂Ωe Ωe ∂t Z Z √ (CM + DQ) · ∇φ dxdydz = F φ dxdydz, − Ωe Ωe Z √ Z Z (21) uxy ∇ · Ψ dxdydz, uxy Ψ · n ds + D−1 Q Ψ dxdydz = − Ωe

∂Ωe

Ωe

6

Dawson, Aizinger

where n is a unit exterior normal to Ωe . This weak formulation is well defined for any uxy (t, x, y, z) ∈ H 1 (0, T ; V d−1 ); w(t, x, y, z) ∈ V, ∀t ∈ [0, T ]; φ(x, y, z) ∈ V d−1 ; Q(t, x, y, z) ∈ Y d−1 , ∀t ∈ [0, T ]; and Ψ (x, y, z) ∈ Y d−1 , where \ def (22) V = L2 (Ω) {u : u¯¯ ∈ H 1 (Ωe ), ∀Ωe ∈ Th }, Ωe

def

(23)

2

Y = L (Ω)

d

\

{q : q¯¯

Ωe

∈ H 1 (Ωe )d , ∀Ωe ∈ Th }.

Next, we seek to approximate (uxy (t, ·), w(t, ·), Q(t, ·)), the solution to the weak problem, with a function (uxy,h (t, ·), wh (t, ·), Qh (t, ·)) ∈ Uh × Wh × Zh , where Uh ⊂ V d−1 , Wh ⊂ V , and Zh ⊂ Y d−1 are some finite-dimensional subspaces. In order to do so, we can use our weak formulation with one important modification. Since our approximation space does not guarantee continuity across the inter-element boundaries, the integrands in the boundary integrals have to be replaced by suitably chosen numerical fluxes preserving consistency and stability of the method. The semi-discrete finite element solution (uxy,h (t, ·), Qh (t, ·)) is obtained by requiring that for any t ∈ [0, T ], all Ωe ∈ Th , and for all (φ, Ψ ) ∈ Uh × Zh the following holds: Z Z √ ∂uxy,h ˆ M,n + DQ ˆ h · n) φ ds (C φ dxdydz + (24) ∂t ∂Ωe Ωe Z Z √ F φ dxdydz, (CM + DQh ) · ∇φ dxdydz = − Ωe Ωe Z Z √ Z ˆ xy,h Ψ · n ds + (25) uxy,h ∇ · Ψ dxdydz, u D−1 Qh Ψ dxdydz = − Ωe

∂Ωe

Ωe

ˆ M,n is a solution to the Riemann problem for the nonlinear boundary flux where C ˆ h are set equal to the average of the values of the corresponding ˆ xy,h , Q CM · n; and u variables on both sides of the boundary. (Note, that there are other possible choices ˆ h .) ˆ xy,h and Q of u Discretization of the mass conservation equation is done in a similar way. Let Ωe,xy be the orthogonal projection of Ωe into the xy-plane. We multiply (16) by a smooth test function δ = δ(x, y), integrate it over Ωe,xy , and integrate by parts. Then, the mass balance in the water column corresponding to the 2D element Ωe,xy can be expressed as Z Z Z ∂H (26) CH · ∇xy δ dxdy = 0, CH · n δ ds − δ dxdy + Ωe,xy ∂t Ωe,xy ∂Ωe,xy ´T ³R Rξ ξ and using the facts that δ is independent of z Noting that CH = zb udz, zb vdz and that H > 0, we can transform the equation above as follows: Z Z X ∂H uxy H · n (27) δ dxdy + δ ds H Ωe ,xy ∂t Ωe ∈col(Ωe,xy ) ∂Ωe,lat Z X uxy · ∇xy δ dxdydz = 0, − Ωe ∈col(Ωe,xy )

Ωe

where ∂Ωe,lat denotes the lateral boundary faces of prism Ωe , and col(Ωe,xy ) is the set of 3D elements in the water column corresponding to Ωe,xy . Note, that the expression

7

A DG method for the 3D SWE def

above is well defined for any δ(x, y) ∈ H = L2 (Ωxy )

T

{h : h¯¯

Ωe,xy

∈ H 1 (Ωe,xy ), ∀Ωe ∈

Th } and H(t, x, y) ∈ H 1 (0, T ; H), where Ωxy is the orthogonal projection of the domain Ω into the xy-plane. Analogous to the momentum conservation equations, we seek an approximation Hh (t, ·) ∈ Hh to the solution of the weak problem H(t, ·), where Hh ⊂ H is some finite dimensional subspace. Using the weak formulation (28) and replacing integrands in the boundary integrals by a suitable numerical flux we obtain our semi-discrete finite element solution Hh (t, ·) ∈ Hh by requiring that for any t ∈ [0, T ], all Ωe ∈ Th , and for all δ ∈ Hh the following holds: Z

(28)

Ωe ,xy



Z X ∂Hh CˆH δ dxdy + δ ds ∂t ξs − z b Ωe ∈col(Ωe,xy ) ∂Ωe,lat Z X uxy,h · ∇xy δ dxdydz = 0,

Ωe ∈col(Ωe,xy )

Ωe

where CˆH is a solution to the Riemann problem for the normal boundary flux uxy H ·n, and ξs denotes the value of the z coordinate at the free surface of the smoothed mesh. This boundary flux formulation has several important advantages. It transforms integrals over 2D edges into integrals over lateral faces of 3D elements thus allowing us to solve the Riemann problem for elevation simultaneously with the corresponding problem for the momentum equations; it is consistent with the continuous formulation; and it takes into account the coupling between velocity and elevation which is crucial for stability of our numerical scheme. Finally, we turn to the space discretization for continuity equation (3). Unlike the mass and momentum conservation equations it is not time-dependent, its main role being computation of the vertical velocity component w to maintain a divergence-free velocity field. Regarding (3) and the boundary condition at the bottom as an initial value problem for w, we can solve it element by element in each water column starting at the bottom and using the solution from the element below to provide an initial condition. Multiplying (3) by a smooth test function σ, integrating it over Ωe , integrating by parts, and re-ordering terms we obtain a weak formulation: µ ¶ Z Z ∂σ (29) dxdydz = w w nz σ ds − ∂z Ωe ∂Ωe,top Z Z Z uxy · n σ ds, uxy · ∇xy σ dxdydz − w nz σ ds − Ωe

∂Ωe,bot

∂Ωe

where ∂Ωe,top , ∂Ωe,bot denote the top and bottom boundaries of element Ωe . We seek wh (t, ·) ∈ Wh , where Wh is some finite dimensional subspace of V , so that for given values of uxy,h (t, ·) ∈ Uh , for all Ωe ∈ Th , and for all σ ∈ Wh the following holds: µ ¶ Z Z ∂σ wh wh nz σ ds − (30) dxdydz = ∂z Ωe ∂Ω Z Z Z e,top − uxy,h · ∇xy σ dxdydz − wh nz σ ds − Cˆw σ ds, Ωe

∂Ωe,bot

∂Ωe

8

Dawson, Aizinger

where wh− is an initial value for wh taken from the element below (or a boundary condition at the bottom in the case of the bottommost element), and Cˆw is a numerical flux for the normal boundary flux function uxy,h · n. On the lateral boundaries Cˆw ˆH should be set equal to ξsC−z (comp. with (29)) in order to preserve the local mass b conservation properties of our numerical scheme. On horizontal faces it can be taken equal to the average or upwinded value of the corresponding variables. 5. Roe’s solver for the Riemann problem. In this section, we will see the main advantage of our new formulation for the surface elevation. Since we were able to modify the discrete mass conservation equation to transform the boundary integrals into 3D form (see (29)), we end up with a well-posed Riemann problem that can be solved to produce a numerical flux on the lateral boundaries satisfying the entropy condition (for a discussion of different entropy conditions see e.g. [11]). Let x be a point on Γ, where Γ is an interior lateral boundary face in our 3D mesh. Let n = (nx , ny , nz ) be a unit normal to Γ at x. We denote s = (H, u, v)T the vector of state variables (we don’t include w in s because it enters the normal fluxes only multiplied by nz , and nz = 0 for lateral boundary faces). Then, we define the left and right states sL , sR at x as follows: (31)

sL = lim s(x + εn), sR = lim s(x + εn). ε→0−

ε→0+

Our task is to compute an entropy solution to the Riemann problem for the normal flux function Cn = (uxy H · n, Cu · n, Cv · n) at x. Using Roe’s solver, we first linearize ˆ n (sL , sR ) s, and then, compute a solution for the linearized our problem Cn (s) ≈ R normal flux. ˆ n (sL , sR ) has to satisfy three conditions (see [14]): The matrix R ˆ n (sL , sR )(sR − sL ) = Cn (sR ) − Cn (sL ); (i) R ˆ n (sL , sR ) is diagonalizable with real eigenvalues; (ii) R ˆ n (sL , sR ) → Cn0 (s) smoothly as sL , sR → s, where (iii) R (32)



unx + vny gnx Cn0 (s) =  gny

Hnx 2unx + vny vnx

 Hny . uny unx + 2vny

ˆ i be the eigenvalues of R ˆ n (sL , sR ) and ˆri the corresponding eigenvectors . Then, Let λ an approximation to the normal flux Cˆn (sL , sR ) ≈ Cn (s) is computed as follows: (33)

Cˆn (sL , sR ) = Cn (sL ) +

3 X

ˆ − ˆri , αi λ i

i=1

ˆ − = min(0, λ ˆ i ) and αi are calculated from: where λ i (34)

3 X i=1

αi ˆri = sR − sL .

ˆ n (sL , sR ) equal to C 0 (s) with s = We claim that setting R n

1 2 (sL

+ sR ) satisfies

A DG method for the 3D SWE

9

ˆ n . Indeed, the eigenvalues and eigenvectors of C 0 (s) are found to be: conditions on R n λ1 = 32 Un − 21 a, (35)

λ2 = 32 Un , λ3 = 32 Un + 21 a,



 H r1 =  u − n2x (Un + a)  ; ny  v − 2(Un + a) 0 r2 =  −ny  ;  nx  H r3 =  u − n2x (Un − a)  ; n v − 2y (Un − a)

p where Un = unx + vny , and a = Un2 + 4gH. Therefore, condition (ii) is satisfied. Clearly, Cn0 (s) → Cn0 (s) smoothly as sL , sR → s. Finally, the first condition can be verified by simply substituting the appropriate values in (i). 6. Numerical Results. In this section, we show some preliminary numerical results obtained with our 3D solver. In the first test, we use a problem with a known smooth analytical solution to obtain convergence rates. Next, we compute a steady state solution to a supercritical flow problem to demonstrate stability of the method in a case involving discontinuities. Our last test example is a standard tidal flow simulation on a domain with complex boundary and bathymetry. In all runs, the time stepping was performed using explicit TVD Runge-Kutta schemes described in [5] of the order matching the order of the space discretization. 6.1. Analytical Test Problem with a Forcing Term. In this section, we will show some numerical results and convergence rates for a problem possessing a smooth analytical solution. Since designing a non-trivial function satisfying our system and boundary conditions is not a simple task, we chose a slightly different approach. We had our test function satisfy the boundary conditions and continuity equation and added a forcing term to the system to make this function a solution to a modified problem. We set vertical and horizontal eddy viscosities equal to zero and tested our numerical solver with the following analytical test example: (36) (37) (38)

(39)

ξ(t, x, y) = e(sin(ω(x + t)) + sin(ω(y + t))), u(t, x, y, z) = d(z − zb ) sin(ω(x + t)),

v(t, x, y, z) = d(z − zb ) sin(ω(y + t)), ¶ µ ∂zb ∂zb sin(ω(x + t)) + sin(ω(y + t)) − w(t, x, y, z) = d(z − zb ) ∂x ∂y 1 d ω(z − zb )2 (cos(ω(x + t)) + cos(ω(y + t))), 2

where e = 0.01, d = 0.1, ω = 0.01. Ωxy , the projection of our 3D domain Ω into the xy-plane, is a square [0, 100] × [0, 100], and the bathymetry varies linearly as zb (x, y) = −(2 − 0.005(x + y)). The time step was set to 1 second for the coarsest mesh and was reduced for finer meshes according to the CFL restriction. The problem was solved for t ∈ [0, T ] with T = 0.001 days. The elevation error was computed at time T using the L2 (Ωxy ) norm, and the error for velocity components was measured in the L2 ([0, T ], Ω) norm. In the tables below, Ne is the number of surface triangles,

10

Dawson, Aizinger Table 1 Piecewise constant approximation, Forward Euler

1 2 3 4 5 6

Ne 2 8 32 128 512 2048

Nl 1 2 4 8 16 32

Er(ξ) 1.088 0.422 0.212 0.104 0.048 0.023

C(ξ) 1.37 0.99 1.03 1.12 1.06

Er(u) 38.323 25.004 12.949 6.591 3.340 1.707

C(u) 0.62 0.95 0.97 0.98 0.97

Er(v) 38.323 25.004 12.949 6.591 3.340 1.707

C(v) 0.62 0.95 0.97 0.98 0.97

Er(w) 0.429 0.311 0.149 0.086 0.054 0.036

C(w) 0.46 1.06 0.79 0.67 0.58

Table 2 Piecewise linear approximation, 2-stage Runge-Kutta

1 2 3 4 5 6

Ne 2 8 32 128 512 2048

Nl 1 2 4 8 16 32

Er(ξ) 2.84e-1 7.25e-2 1.74e-2 4.24e-3 1.05e-3 2.63e-4

C(ξ) 1.97 2.06 2.04 2.01 2.00

Er(u) 8.36 2.68 0.82 0.23 0.062 0.017

C(u) 1.64 1.71 1.83 1.89 1.87

Er(v) 8.39 2.69 0.82 0.23 0.062 0.017

C(v) 1.64 1.71 1.83 1.89 1.87

Er(w) 2.28e-1 1.01e-1 5.03e-2 2.71e-2 1.48e-2 8.06e-3

C(w) 1.17 1.01 0.89 0.87 0.88

Table 3 Piecewise quadratic approximation, 3-stage Runge-Kutta

1 2 3 4 5

Ne 2 8 32 128 512

Nl 1 2 4 8 16

Er(ξ) 1.57e-2 3.42e-3 4.61e-4 6.96e-5 1.30e-5

C(ξ) 2.20 2.89 2.73 2.42

Er(u) 1.34e0 2.18e-1 2.80e-2 3.65e-3 4.90e-4

C(u) 2.62 2.96 2.94 2.90

Er(v) 1.34e0 2.18e-1 2.80e-2 3.65e-3 4.90e-4

C(v) 2.62 2.96 2.94 2.90

Er(w) 4.35e-2 1.18e-2 2.28e-3 5.57e-4 1.73e-4

C(w) 1.88 2.37 2.03 1.69

Nl the maximum number of equidistant vertical layers, Er(.) norm of the error, and C(.) convergence rate. Although the theory of DG methods does not contain any a priori error estimates for this class of problems, our numerical results (with an exception of those for w) indicate behavior similar to linear problems. In view of the fact that the vertical velocity component w is computed by solving the continuity equation (31) with the numerical approximations for u and v used as input data, unsurprisingly, it converges at a slower rate than u and v. 7. Supercritical Flow Problem. In this section, we test our numerical model on a 2D problem with a known analytical solution and compare results from 2D and 3D simulations. Supercritical channel flows subject to a change in cross-section can lead to formation of shock and rarefaction waves. Here, we take one particular configuration given in Zienkiewicz and Ortiz [17]. We are given a channel of uniform depth whose boundary wall is constricted on both sides with a constriction angle of 50 resulting in a cross-wave pattern. Flow is induced through the left boundary with no flow through the top and bottom boundaries. The inlet Froude number is defined by F r = √U , where U and H are the mean velocity and total depth of the fluid at gH

11

A DG method for the 3D SWE

the inlet, respectively. When F r > 1, the flow is said to be supercritical. For our test, we chose F r = 2.5. This problem can be solved analytically using methods presented by Ippen in [9]. The discontinuous true solution projected on our (smoothed) mesh is demonstrated in Figure 1. The coarse 2D mesh consists of 3155 triangular elements with no particular orientation. We project these elements vertically to obtain a 3D grid with prismatic elements. The 3D mesh can be refined horizontally by subdividing every triangle in four using edge bisection, or/and vertically by subdividing into layers. In our test runs, we set vertical and horizontal eddy viscosities equal to zero and used a no-drag boundary condition at the bottom.

Z 0.85 0.75 0.65 0.55 0.45 0.35 0.25 0.15 0.05 -0.05 -0.15 -0.25 -0.35 -0.45 -0.55 -0.65 -0.75 -0.85 -0.95

70 60 50

0.5 40

-0.5

20

Z

0 30

X

-1 0

10 0

10 20

-10 -20

30 40

Y

Fig. 1. Constricted channel flow problem: true solution, 1 layer

As we can see in Figures 2–3, the numerical solutions for this problem, certainly extremely demanding in terms of stability requirements, show no oscillatory behavior. This fact is especially interesting in the case of piecewise linear approximation because no slope limiters or any other stabilization techniques were employed in the course of the simulation. We also verified that increasing the number of layers does not affect the solution of this problems substantially (comp. table (4)) just as expected in case of an essentially 2D problem. In Table 4, we list L2 norms of the error in surface elevation for numerical solutions computed using different mesh refinements and approximation spaces. For comparison, we also show errors obtained using the 2D version of our simulator. 8. Tidal flow near Bahamas Islands. In this test run, we simulate tide-driven flow near the Bahamas Islands. The domain geometry and the coarse finite element mesh consisting of 1696 elements are shown in Figure 4. Finer meshes were obtained by subdividing each triangle into four using edge bisection (we refer to this mesh as “refined”) and by subdividing into vertical layers. The bathymetry is shown in Figure 5.

12

Dawson, Aizinger

Z 0.85 0.75 0.65 0.55 0.45 0.35 0.25 0.15 0.05 -0.05 -0.15 -0.25 -0.35 -0.45 -0.55 -0.65 -0.75 -0.85 -0.95

70 60 50

0.5 40

-0.5

20

Z

0 30

X

-1 0

10 0

10 20

-10 -20

30 40

Y

Fig. 2. Constricted channel flow problem: piecewise constants, 10 layers

Z 0.85 0.75 0.65 0.55 0.45 0.35 0.25 0.15 0.05 -0.05 -0.15 -0.25 -0.35 -0.45 -0.55 -0.65 -0.75 -0.85 -0.95

70 60 50

0.5 40

-0.5

20

-1 0

10 0

10 20

-10 -20

30 40

Y

Fig. 3. Constricted channel flow problem: piecewise linears, 10 layers Table 4 Error in surface elevation for supercritical flow problem

Model type 3D 2D

Coarse mesh, constants, 1 layer 2.85 2.75

Coarse mesh, constants, 10 layers 2.85 N/A

Refined twice mesh, constants, 1 layer 1.74 1.61

Coarse mesh, linears, 1 layer 1.53 1.42

Coarse mesh, linears, 10 layers 1.57 N/A

Z

0 30

X

13

A DG method for the 3D SWE

1.0x10

5

8.0x10

4

OPEN SEA

1.0x10

+05

8.0x10

+04

6.0x10

+04

DP 8.5 8 7.5 7 6.5 6 5.5 5 4.5 4 3.5 3 2.5 2 1.5

5-COMP. TIDE

4

4.0x10

4

2.0x10

4

Y

6.0x10

LAND Un = 0

1 4.0x10+04

4 3

2.0x10

+04

2 0.0x100

2.0x104

4.0x104

6.0x104

0.0x10+00

8.0x104

2.0x10+04

4.0x10+04

6.0x10+04

8.0x10+04

X

Fig. 4. Mesh for tide-driven flow past an island (Bahamas). Lengths are in meters.

Fig. 5. Bathymetry for tide-driven flow past an island (Bahamas). Lengths are in meters.

The following tidal forcing with time(t) in hours was imposed at the open sea boundary: ˆ ξ(t) (40)

t = 0.075 cos( 25.82 t + 0.095 cos( 23.94 t + 0.100 cos( 12.66 t + 0.395 cos( 12.42 t + 0.060 cos( 12.00

+ 3.40) + 3.60) + 5.93) + 0.00) + 0.75) (meters).

The simulations were cold-started and the tidal forcing was imposed gradually over a period of two days. The Coriolis parameter was set to 3.19 × 10−5 s−1 . The horizontal eddy viscosity was set equal to zero, the vertical eddy viscosity was modeled using a quadratic formula given by Davies in [7]: 1

(41)

vt = K t

(¯ u2 + v¯2 ) 2 , ωa

where u ¯, v¯ are depth averaged horizontal velocity components, ωa a typical long wave frequency taken as 10−4 s−1 , and Kt = 2 × 10−5 a dimensionless coefficient. The elevation and velocities were monitored at four different stations, whose coordinates in meters are (38666.66, 49333.32), (56097.79, 9612.94), (41262.60, 29775.73), and (59594.66, 41149.62), for the duration of days eleven and twelve of the simulation at an increment of six hundred seconds. A constant time step of 20 seconds was used in the runs performed on the coarse mesh, and a time step of 10 seconds was taken for the refined mesh. We conducted a number of numerical experiments using our 3D solver on this problem with different approximation spaces and meshes. In Figures 6 and 7 we compare output at the recording stations to the results obtained using the 2D ADCIRC

14

Dawson, Aizinger

model of Luettich et al [12] and UTBEST, our DG solver for the 2D depth-integrated equations.. Generally, all time series agree rather well, with only one exception being the x component of velocity at station 4, where we saw some significant disagreement. For sake of conciseness, we only include output from recording stations 3 and 4 noting that the results at the other two stations displayed good agreement similar to station 3. For readability, we plotted data for each recording station data on two separate graphs, the first one comparing results from both 2D models to the 3D LDG model with piecewise constants on various meshes and 3D LDG with piecewise linears on the coarse mesh. The second graph compares several higher order 3D LDG solutions. Commenting on the results of these numerical experiments, we first want to note that the 3D model is extremely sensitive to the choice of vertical eddy viscosity closure, whereas the 2D models are just as sensitive to the bottom friction formulation. Nevertheless, the results of 2- and 3D models displayed an excellent agreement in phase and a good agreement in amplitude (the last one with an exception of the anomaly at recording station 4). The plots for higher order approximations show that the piecewise linear solution on the coarse mesh is already sufficiently accurate so that refining the mesh or using piecewise quadratic polynomials changes the solution very little. The results of the simulation that employs piecewise constants for elevation and piecewise linears for velocity are also very close to the fully linear solution. This fact is very important since we want to have an option of using piecewise constants for problems with “rough” surface elevation. The graphs for piecewise constant polynomial spaces merit a closer look. Though rather similar to the piecewise linear solution in most cases, they can be a poor approximation at some points in space and time, and even refining the mesh does not guarantee convergence to the linear solution at every point (compare the left plot in the second row in Figure 7). The theory of LDG methods for linear convectiondiffusion problems proves convergence of order ∆xk for general triangulations, where k ≥ 1 is the order of the polynomial space [4]. Thus, there is no guarantee of convergence for piecewise constant approximations. Based on the results of the numerical experiments for piecewise constant approximation, we can conclude that our model displays behavior in agreement with the linear theory. 9. Conclusions. In this paper, we have described an implementation of the LDG method for the 3D shallow water equations, and given some preliminary numerical results. The methodology shows promise for modeling both high Froude number and tidal flows. More extensive testing of this algorithm and analysis of its convergence properties is underway and will be reported in future papers. REFERENCES [1] V. Aizinger, C. Dawson, B. Cockburn and P. Castillo, The local discontinuous Galerkin method for contaminant transport, Advances in Water Resources, 24, pp. 73-87, 2001. [2] V. Aizinger and C. Dawson, A discontinuous Galerkin method for two-dimensional flow and transport in shallow water, Advances in Water Resources, 25, pp. 67-84, 2002. [3] S. Chippada, C. N. Dawson, M. Martinez and M. F. Wheeler, A Godunov-type finite volume method for the system of shallow water equations, Comput. Meth. Appl. Mech. Engrg., 151, pp. 105-129, 1998. [4] B. Cockburn and C.-W. Shu, The local discontinuous Galerkin finite element method for convection-diffusion systems, SIAM J. Numer. Anal., 35, pp. 2440-2463, 1998. [5] B. Cockburn and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite

A DG method for the 3D SWE

[6]

[7] [8]

[9] [10]

[11] [12]

[13] [14] [15] [16] [17]

15

element method for scalar conservation laws II: General framework, Math. Comp., 52, pp. 411-435, 1989. B. Cockburn, G. Karniadakis and C.-W. Shu, The development of discontinuous Galerkin methods, in Discontinuous Galerkin Methods: Theory, Computation and Applications, B. Cockburn, G. Karniadakis and C.-W. Shu, editors, Lecture Notes in Computational Science and Engineering, volume 11, Part I: Overview, pp.3-50, Springer, 2000. A. M. Davies, A three-dimensional model of the Northwest European continental shelf, with application to the M4 tide, J. Phys. Oceanogr., 16(5), pp. 797-813, 1986. Ip, J.T.C. and D.R. Lynch, QUODDY3 User’s Manual: Comprehensive Coastal Circulation Simulation using Finite Elements: Nonlinear Prognostic Time-Stepping Model, Thayer School of Engineering, Dartmouth College, Hanover, New Hampshire, Report Number NML95-1, 1995. A. T. Ippen Mechanics of supercritical flow, Trans. of ASCE, 116, pp. 268-295, 1951. Johnson, B.H., Kim, K.W., Heath, R.E., Hsieh, B.B. and Butler, H.L., Development and verification of a three-dimensional numerical hydrodynamic, salinity and temperature model of Chesapeake Bay, technical report HL-91-7, U.S. Army Engineer Waterways Experiment Station, Vicksburg, MS, 1991. R.J. Le Veque, Numerical Methods for Conservation Laws, Basel, Birkh¨ auser, 1992. R.A. Luettich, J.J. Westerink, N.W. Scheffner, ADCIRC: An Advanced Three-Dimensional Circulation Model for Shelves, Coasts and Estuaries, Report 1, U.S. Army Corps of Engineers, Washington, D.C. 20314-1000, December 1991. Lynch, D. R., and F. E. Werner, Three-dimensional hydrodynamics on finite elements. Part II: Non-linear time-stepping model, Intl. J. Numer. Meths. Fluids., 12: 507-533, 1991. P.L. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes, J. Comput. Phys., 43, pp. 357-372, 1981. C. B. Vreugdenhil, Numerical Methods for Shallow-Water Flow, Kluwer Academic Publishers, 1994. T. Weiyan, Shallow Water Hydrodynamics, Elsevier Oceanography Series, 55, Elsevier, Amsterdam, 1992. O.C. Zienkiewicz and P. Ortiz, A split-characteristic based finite element model for the shallow water equations, Int. J. Num. Meth. Fluids, 20, pp. 1061-1080, 1995.

16

Dawson, Aizinger Elevation at recording station 3

Elevation at recording station 3

2D ADCIRC, coarse mesh 2D UTBEST, p.w. linear, coarse mesh 3D p.w. const, 1 layer, coarse mesh 3D p.w. const, 5 layers, refined mesh 3D p.w. linear, 1 layer, coarse mesh

0.3

0.3

0.2

elevation, m

0.2

elevation, m

3D p.w. linear, 1 layer, coarse mesh 3D p.w. linear, 5 layers, refined mesh 3D p.w. quadratic, 1 layer, coarse mesh 3D p.w. constant elevaton and linear velocity, 1 layer, coarse mesh

0.1

0.1

0

0

−0.1

−0.1

−0.2 8.6

8.8

9

9.2

9.4 9.6 time, sec

9.8

10

10.2

−0.2 8.6

10.4

8.8

9

x−velocity at velocity recording station 3

0.03

x−velocity, m/s

x−velocity, m/s

10.2

10.4 5

x 10

0.02

0

0.01

0

−0.01

−0.01

−0.02

−0.02

−0.03

−0.03 8.8

9

9.2

9.4 9.6 time, sec

9.8

10

10.2

10.4

8.6

8.8

9

9.4 9.6 time, sec

9.8

10

10.2

10.4 5

x 10

y−velocity at velocity recording station 3

2D ADCIRC, coarse mesh 2D UTBEST, p.w. linear, coarse mesh 3D p.w. const, 1 layer, coarse mesh 3D p.w. const, 5 layers, refined mesh 3D p.w. linear, 1 layer, coarse mesh

3D p.w. linear, 1 layer, coarse mesh 3D p.w. linear, 5 layers, refined mesh 3D p.w. quadratic, 1 layer, coarse mesh 3D p.w. constant elevaton and linear velocity, 1 layer, coarse mesh 0.1

y−velocity, m/s

0.1

9.2

5

x 10

y−velocity at velocity recording station 3

y−velocity, m/s

10

0.03

0.01

0

−0.1 8.6

9.8

3D p.w. linear, 1 layer, coarse mesh 3D p.w. linear, 5 layers, refined mesh 3D p.w. quadratic, 1 layer, coarse mesh 3D p.w. constant elevaton and linear velocity, 1 layer, coarse mesh

0.04

0.02

8.6

9.4 9.6 time, sec

x−velocity at velocity recording station 3

2D ADCIRC, coarse mesh 2D UTBEST, p.w. linear, coarse mesh 3D p.w. const, 1 layer, coarse mesh 3D p.w. const, 5 layers, refined mesh 3D p.w. linear, 1 layer, coarse mesh

0.04

9.2

5

x 10

0

−0.1 8.8

9

9.2

9.4 9.6 time, sec

9.8

10

10.2

10.4 5

x 10

8.6

8.8

9

9.2

9.4 9.6 time, sec

9.8

10

10.2

10.4 5

x 10

Fig. 6. Comparison of elevation and velocity approximations at Bahamas station 3 for days 11-12.

17

A DG method for the 3D SWE Elevation at recording station 4

0.25

0.25

0.2

0.2

0.15

0.15

0.1

0.05

0.1

0.05

0

0

−0.05

−0.05

−0.1

−0.1

−0.15 8.6

3D p.w. linear, 1 layer, coarse mesh 3D p.w. linear, 5 layers, refined mesh 3D p.w. quadratic, 1 layer, coarse mesh 3D p.w. constant elevaton and linear velocity, 1 layer, coarse mesh

0.3

elevation, m

elevation, m

Elevation at recording station 4

2D ADCIRC, coarse mesh 2D UTBEST, p.w. linear, coarse mesh 3D p.w. const, 1 layer, coarse mesh 3D p.w. const, 5 layers, refined mesh 3D p.w. linear, 1 layer, coarse mesh

0.3

−0.15 8.8

9

9.2

9.4 9.6 time, sec

9.8

10

10.2

10.4

8.6

8.8

9

x−velocity at velocity recording station 4

0.02

0.01

0.01

0

−0.01

−0.03

−0.03

−0.04

−0.04

9.2

9.4 9.6 time, sec

9.8

10

10.2

−0.05 8.6

10.4

8.8

9

9.8

10

10.2

10.4 5

x 10

0.2

0.15

0.1

0.1 y−velocity, m/s

y−velocity, m/s

9.4 9.6 time, sec

3D p.w. linear, 1 layer, coarse mesh 3D p.w. linear, 5 layers, refined mesh 3D p.w. quadratic, 1 layer, coarse mesh 3D p.w. constant elevaton and linear velocity, 1 layer, coarse mesh

0.25

0.15

0.05

0

0.05

0

−0.05

−0.05

−0.1

−0.1

−0.15

−0.15

−0.2 8.6

5

y−velocity at velocity recording station 4

2D ADCIRC, coarse mesh 2D UTBEST, p.w. linear, coarse mesh 3D p.w. const, 1 layer, coarse mesh 3D p.w. const, 5 layers, refined mesh 3D p.w. linear, 1 layer, coarse mesh

0.2

9.2

5

x 10

y−velocity at velocity recording station 4 0.25

10.4 x 10

−0.01

−0.02

9

10.2

0

−0.02

8.8

10

0.03

0.02

−0.05 8.6

9.8

3D p.w. linear, 1 layer, coarse mesh 3D p.w. linear, 5 layers, refined mesh 3D p.w. quadratic, 1 layer, coarse mesh 3D p.w. constant elevaton and linear velocity, 1 layer, coarse mesh

0.04

x−velocity, m/s

x−velocity, m/s

0.03

9.4 9.6 time, sec

x−velocity at velocity recording station 4

2D ADCIRC, coarse mesh 2D UTBEST, p.w. linear, coarse mesh 3D p.w. const, 1 layer, coarse mesh 3D p.w. const, 5 layers, refined mesh 3D p.w. linear, 1 layer, coarse mesh

0.04

9.2

5

x 10

−0.2 8.8

9

9.2

9.4 9.6 time, sec

9.8

10

10.2

10.4 5

x 10

8.6

8.8

9

9.2

9.4 9.6 time, sec

9.8

10

10.2

10.4 5

x 10

Fig. 7. Comparison of elevation and velocity approximations at Bahamas station 4 for days 11-12.