- Email: [email protected]
A three-dimensional model for suspended sediment transport based on the compact discontinuous Galerkin method
Author's Accepted Manuscript
A three-dimensional model for suspended sediment transport based on the compact discontinuous Galerkin method Zhao Zhang-Yi, Zhang Qing-He, Zhao Hong-Bo, Yang Hua
www.elsevier.com/locate/ijsrc
PII: DOI: Reference:
S1001-6279(15)00043-8 http://dx.doi.org/10.1016/j.ijsrc.2014.06.001 IJSRC31
To appear in:
International Journal of Sediment Research
Received date: 4 February 2013 Revised date: 17 June 2014 Accepted date: 29 June 2014 Cite this article as: Zhao Zhang-Yi, Zhang Qing-He, Zhao Hong-Bo, Yang Hua, A three-dimensional model for suspended sediment transport based on the compact discontinuous Galerkin method, International Journal of Sediment Research, http://dx.doi.org/10.1016/j.ijsrc.2014.06.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A three-dimensional model for suspended sediment transport based on the compact discontinuous Galerkin method Zhao Zhang-Yi1, 2, Zhang Qing-He1, Zhao Hong-Bo2, Yang Hua2 1
State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin
University, Tianjin 300072, China 2
Key Laboratory of Engineering Sediment, Ministry of Transport, Tianjin Research
Institute for Water Transport Engineering, M. O. T., Tianjin 300456, China Corresponding author: Zhang Qing-He Tel: +86-22-27404266; fax: +86-22-27404266. E-mail address: [email protected]
A three-dimensional model for suspended sediment transport based on the compact discontinuous Galerkin method Zhao Zhangyi a, b, Zhang Qinghe a,, Zhao Hong-Bo b, Yang Hua b a
State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University, Tianjin 300072, China
b
Key Laboratory of Engineering Sediment, Ministry of Transport, Tianjin Research Institute for Water Transport Engineering, M. O. T., Tianjin 300456, China Abstract A three-dimensional numerical model was developed in the present paper, with the goal of simulating suspended sediment transport in discontinuous situations. In the model, the advection-diffusion equation of suspended sediment was solved using the Compact Discontinuous Galerkin (CDG) method based on unstructured meshes. The evaluation of the model revealed good agreement between the proposed model and existing theories. The model has second-order accuracy. Moreover, the model was applied to an idealized case of suspended sediment transport caused by rip currents based on the field spatial scale. Key Words: advection-diffusion equation; model evaluation; discontinuity; rip currents; unstructured meshes; slope limiter
1 Introduction Numerical modeling has become a powerful tool in the contemporary study of sediment transport in rivers, estuaries and coastal areas. Based on various purposes, a number of related studies have been completed. For example, Wang and Adeff (1986) proposed a 3-D model for suspended sediment transport based on the finite element method to study river sedimentation. Lin and Falconer (1996) established a layer-integrated 3D model for the estuarine and coastal waters around the UK using a finite difference method. Wu et al. (2000) developed a 3-D model to investigate the transport properties of both the suspended sediment and the bed load based on the finite volume method. Khosronejad et al. (2007) developed a 3-D model for sediment transport in curved open channels based on the finite volume method. Recently, the suspended module of sediment transport has become a regular component of an environmental model, including both commercial models and community models, e.g., EFDC (Environmental Fluid Dynamics Code) (Hamrick, 1992), FVCOM (Finite Volume Coastal Ocean Model) (Chen et al., 2003), ROMS (Regional Ocean Modeling System) (Warner et al., 2008), and Delft-3D (Deltares, 2010). Most of the above-mentioned models are only suitable for continuous
Corresponding author. Email address: [email protected]
sediment transport. However, in discontinuous or large gradient flows, such as dam-breaking flow and tidal bore flow, the discontinuous transportation of sediment may exist. In this case, standard numerical methods, such as the finite element or finite volume method, may generate solutions with numerical diffusion and/or non-physical oscillations (Younes & Ackerer 2008). Therefore, there is practical-significance to develop a high-resolution suspended sediment transport model using a discontinuity-capturing method; a candidate for such a model is the Discontinuous Galerkin (DG) method (Reed & Hill, 1973). Currently, in addition to Gourgue’s work (2011) on depth-averaged 2-D model applied in the Scheldt estuary, there are only a few published papers regarding the suspended sediment transport model based on the DG method although there have been some research reports on the application of the DG method for advection-diffusion problems (Cockburn & Shu, 1998b; Gómez et al., 2007; Zhang et al., 2013). Therefore, we will develop a 3-D suspended sediment transport model based on the DG method with unstructured meshes. The DG method used in this study is based on the so-called compact discontinuous Galerkin (CDG) method, which was first presented by Peraire and Persson (2008). This method is an extension of the Runge-Kutta DG method (Cockburn & Shu, 1998a; 2001) and is closely related to the local discontinuous Galerkin (LDG) method (Cockburn & Shu, 1998b). Peraire and Persson (2008) proved that the CDG method inherits all of the attractive features of the LDG method. In addition, the CDG method is more compact and stable than the LDG method in situations involving multiple dimensions. The outline of the paper is as follows. The mathematical formulation of 3-D suspended sediment transport is briefly introduced in Section 2. The numerical scheme and the spatial discretization of the CDG method as well as its temporal discretization are described in Section 3. Some numerical test calculations are presented for evaluating the performance of the model in Section 4. The results of a numerical example for suspended sediment transport are presented to demonstrate the applicability of the proposed model in Section 5. Finally, some conclusions are presented in Section 6. 2 Mathematical formulation of suspended sediment transport 2.1 Governing equation The distribution of the suspended sediment concentration in water is governed by the following advection-diffusion equation: c t (U c)
K ( c) 0
(1) where c is the suspended sediment concentration, ( x , y , z) , U = (ux, uy, uz-ωs) , ux, uy and T
uz are the velocity components, ωs is the settling velocity of suspended sediment, and K ( x
0 0), (0 y
0), (0 0 z )
T
where εx, εy and εz are the mixing coefficients of the sediment in the x-, y- and z-directions, respectively. In applying Eq. (1) to rivers, estuaries and coastal areas where the ratio of the vertical length scale to the horizontal length scale is small, the vertical diffusivity terms are usually several orders of magnitude larger than the horizontal diffusivity terms (Lin & Falconer, 1996). Therefore, it is more important to exactly calculate the vertical mixing coefficients than the horizontal ones. In the present paper, the horizontal mixing coefficients are assumed to be constant, while the vertical mixing coefficient was prescribed using the method given by van Rijn (1984).
2.2 Boundary conditions Eq. (1) is solved with the following boundary conditions: at the free surface, the vertical sediment flux is zero, and hence, a no-flux boundary condition is used for the sediment concentration: z
c s c 0 z
(2) At the bottom boundary of the suspended sediment layer, the sediment flux is the difference between the rates of deposition and erosion: z
c s c Eb Db z
(3)
The deposition rate at the interface is Db sc
(4) The erosion rate is calculated as described by Ariathurai and Arulanandan (1978): Eb E0 (1 Pb ) b c 1
(5) where E0 is the erosion flux parameter, Pb is the bottom porosity, τb is the bottom shear stress, and τc is the critical shear stress for erosion. 3 Numerical discretization of the governing equation by the DG method 3.1 Definition The following definition is used in the proposed method c g D , on D
(6a) Kc n g N , on N
(6b) where Ω is a bounded domain with boundary
D
, ∂Ω N
D
N
and ∂Ω are the Dirichlet and
Neumann boundaries, respectively, gD and gN are specified values, and n is the outward unit normal to the boundary of Ω. 3.2 The spatial discretization For the sake of simplicity, Eq. (1) is rewritten as the two respective parts to be discretized. c t (Uc) f
(7a) (Kc) f
(7b) 3.2.1 Discretization of the advection terms The bounded domain Ω is partitioned into N non-overlapping elements. To apply the DG method to Eq. (7a), the equation is first multiplied by the trial or test function φ. The resulting equation is integrated over an element e with the face ∂e, and after integration by parts, the equation can be written as
c
( t ( U)c)de
e
e
ˆ cˆd(e) f de n U e
(8) where the hatted terms denote the numerical fluxes across the element faces. Here, the averaged velocity on each side of the face is an approximation to Uˆ , and cˆ is taken to be the upwind value (Fig. 1). If we use cupw to denote the value of c on the upwind side of the face, then by summing Eq. (8) over all elements and considering Eq. (6), we have the following global expression:
(
c ˆ c)d U D n Uˆ ( gD c)d() s n Uˆ (cupw c)ds f d t
(9) where s is the union of the interior faces of all elements, s- is the inflow part of s, and
D
and D
are the outflow and inflow Dirichlet boundaries, respectively.
(a) when the flow is outflow.
(b) when the flow is inflow. Fig. 1
Illustration of the upwind value
3.2.2 Discretization of the diffusion terms Following Cockburn and Shu (1998a), we introduce the new variable ξ and rewrite Eq. (7b) as a first order system of equations ξ f
(10a) ξ Kc
(10b) Introducing a scalar test function φ and a vector test function λ, multiplying Eqs. (10a) and (10b) by the respective test functions, integrating over an element e, and then integrating by parts, we obtain ξ d e n ξˆd(e) f de e
e
e
(11) λKcde n λKcˆd(e) λ ξde e
e
e
(12) where the numerical fluxes ξˆ and cˆ are approximations to ξ and c, respectively, on the face of the element. Following the notation of Arnold et al. (2002), we consider two neighboring elements e+ and e- with the common face k e
e , and assume that n denotes the unit normal vectors to the element
boundary e . Let λ and denote the test functions on k. The average and jump operators are defined as follows
{λ} (λ λ ) 2, {} ( ) 2
(13a) [λ] λ n λ n , [ ] n n
(13b) By adding Eqs. (11) and (12) over all elements, the following global equations are obtained
ξ d f d ξˆ [ ]ds ξˆ nd()
s
(14)
ξ λd Kc λd Kcˆ[λ ]ds Kcˆλ nd()
s
(15) Next, after performing integration by parts on the first term on the right side of Eq. (15), we obtain
ξ λd λ (Kc)d K([c] λ {cˆ c}[λ])ds K(cˆ c)λ nd()
s
(16) The numerical fluxes ξˆ and cˆ for the interior faces in Eqs. (15) and (16) are given by ξˆ {ξ e} C11[c] C12[ξ e ]
(17) cˆ {c} C12 [c]
(18) and the numerical fluxes for the boundary faces are ξˆ ξe C11 (c g D )n, cˆ g D , on D
(19a) ξˆ g N n, cˆ c, on N
(19b) e
e
where C11 is a positive constant, C12 is a vector, and ξ c ξ . Making use of Eqs. (18) and (19), Eq. (16) can be written as
ξ λd λ (Kc)d K([c] {λ} C12 [c][λ])ds
D
s
K( g D c)λ nd()
(20) To obtain an expression for ξ as a function of c, we follow Peraire and Persson (2008) and introduce the local lifting operators:
r e (ξ) λ ξ {λ}, e s e
(21a)
r e (c) λ c[λ ], e s e
(21b)
r (c) λ cλ n, e D
e D
e
(21c)
Setting λ in Eq. (20), Eq. (14) can be rewritten as
(Kc)d K([c] {} C12 [c][ ])ds
D
s
K( g D c) nd() f d ξ [ ]ds ξ nd()
s
(22) According to Eqs. (17), (19) and (21), the terms involving ξˆ in Eq. (22) can be written as
ξˆ [ ]ds (({c} C s
12
s
[c])[ ] C11[c] [ ])ds
(r e ([ ]) r e (C12 [ ])) (r e ([c]) r e (C12 [c]) rDe (c))d (r e ([ ]) r e (C12 [ ])) rDe ( g D )d es
es
(23)
ξˆ nd()
D
(c n C11 (c g D ))d()
eD
eD
e
rDe ( g D ) nd(e)
N
g N d()
r ( ) (r e ([c]) r e (C12 [c]) rDe (c))d()
e D
(24) Therefore, Eqs. (9), (22)~(24) comprise the weak formulation of the advection-diffusion equation. Because the second-order accurate numerical model will be developed as described by Cockburn and Shu (2001), the DG method with pth (p = 0, 1, 2...) order polynomials is always (p+1)th order accurate. The linear solution c on a tetrahedron element e in the weak formulation of the advection-diffusion equation can be expressed in terms of local basis functions i and degrees of freedom ci: Nloc
c cii i 1
(25) where Nloc is the number of vertices of the element. 3.3 Temporal discretization Making use of Eq. (25), the weak formulation of the advection-diffusion equation produces a semi-discrete equation that can be written in matrix form as follows M
dc A(U)c Dc r dt
(26) where M is the mass matrix, A(U) is the advection matrix, D is the diffusion matrix and r is the right-hand side vector containing boundary, source and absorption terms. When solving the 3D advection-diffusion equation for realistic problems, it is advantageous to use a large simulation time step to achieve fast CPU times, and it is suitable to use an implicit scheme without a strong restriction of the size of the time step (Hartmann & Houston, 2005). However, a fully implicit scheme leads to nonlinear algebraic equations at each time level, which is rather expensive to solve (Dolejší et al., 2011). Therefore, only the advection terms in the advection-diffusion equation are treated explicitly, while the diffusion terms are treated implicitly in the proposed model. To achieve this, Eq. (26) is discretized in two steps as follows: M c* cn t A(Un )cn rDn
(27)
M cn 1 c* t Dc n 1 rNn 1 rsn 1
(28) where c* is the preliminary value, rD, rN, and rs are the Dirichlet, Neumann boundary and source terms, respectively. 3.4 Slope limiter The high-order method often produces spurious oscillations and instability in the vicinity of discontinuities (Burbeau et al., 2001). To ensure the stability of the high-order scheme and eliminate the non-physical oscillation, a slope limiter, which can limit the solution gradient near discontinuities, is usually used. In recent years, a number of slope limiting techniques are available for the DG method (Cockburn & Shu, 1998b; Burbeau et al., 2001; Krivodonova et al, 2004; Kuzmin, 2010). However, the slope limiter may introduce a certain amount of numerical diffusion while suppressing oscillation and often requires parameter adjustment (Kuzmin, 2010). The modified minmod limiter and the vertex-based limiter, with only one parameter or no user-adjusted parameters, will be compared in Section 4. A brief introduction regarding the two slope limiters for the second-order DG method will first be presented here. 3.4.1 Modified minmod limiter The modified minmod limiter (Cockburn & Shu, 1998b) is defined as follows: m(a1 , a2 ,
a1 , if a1 M x 2 , an ) otherwise m(a1 , a2 , , an ),
(29) where M is a positive number, and m is the minmod limiter function (Seby, 1984). More detailed discussion about how sensitive the numerical model results are to the parameter M will be provided in Section 4. m(a1 , a2 ,
, an ) s min ai if sgn(a1 ) 1 i n
sgn(an ) s
(30) 1,
x0
where sgn( x) 0, x 0 . 1, x 0
To construct the modified slope limiter for tetrahedral elements, we proceed as follows. Considering the tetrahedrons in Fig. 2, where m1 is the barycenter of the face of e0, and bi denotes the barycenter of the tetrahedron ei (i = 0, 1, 2, 3, 4).
Fig. 2
The tetrahedron element e0 and its four neighbors
Because we have 4
m1 b0 i (bi b0 ) i 1
(31) where αi (i = 1, 2, 3, 4) depends on m1 and the geometry, we can obtain for any linear function c 4
c(m1 ) c(b0 ) i (c(bi ) c(b0 )) i 1
(32) and cei
cde
ei c(bi ), i 0,1,2,3,
ei
(33) we have 4
c (m1 , e0 ) c(m1 ) ce0 i (cei ce0 ) c (m1, e0 ) i 1
(34) Supposing c is a piecewise linear function, and let mi (i = 1, 2, 3, 4) be the four barycenters of the faces of the e0, we obtain 4
4
i 1
i 1
c( x, y, z ) c(mi ) i ( x, y, z ) ce0 c(mi , e0 ) i ( x, y, z), ( x, y, z) e0
(35) where βi (i = 1, 2, 3, 4) depends on the geometry. To use the slope limiter for solution, the following quantities are computed first: i m(c(mi , e0 ), c (mi , e0 )) (36) where m is the modified minmod limiter defined in Eq. (29), and υ is the limiter factor (Cockburn & Shu, 1998a), taking υ = 1.1 in our model. If
4
i 1
i
0 , then the limiter is only applied to c* in Eq. (27), and we have 4
c* ( x, y, z ) ce0 i i ( x, y, z ), ( x, y, z ) e0 i 1
(37) If
4
i 1
i
0 , then the following variables are computed: 4
4
i 1
i 1
pos= max(0, i ), neg= max(0, i )
(38) + = min 1,neg pos , min 1,pos neg
(39) Next, the limiter is applied to c* in Eq. (27), and we have
4
c* ( x, y, z ) ce0 ˆ i i ( x, y, z ), ( x, y, z ) e0 i 1
(40) where ˆ i max(0, i ) max(0, i ) .
3.4.2 Vertex-based limiter The vertex-based limiter was introduced by Kuzmin (2010) and is defined as
min 1, cimax ce ci ce , ci ce 0 ve min 1, ci ce 0 i min min 1, ci ce ci ce , ci ce 0
where c
max i
max(ce , c
max i
), c
min i
(41) min(ce , c ) , ci is the value on the node of the tetrahedron e, and ce min i
denotes the value on the barycenter of the tetrahedral element e. Note that cimax and cimin are initialized by small and large constants, respectively. The vertex-based limiter is only applied to the DG field c* in Eq. (27). For linear elements, c* in any element may be written as c* ( x, y, z) ce ve (ci ce ), ( x, y, z) e
(42) 4 Numerical tests 4.1 Accuracy test The accuracy of the second-order CDG scheme is tested by solving c (uc) (vc) 2 2 2 ( uc) 2 ( vc) 0 t x y x y
(43) in a three-dimensional grid system. The domain is defined on ( x, y, z) [0, 1]×[0, 1] ×[0, 1] with velocity u = v = 1.0 and periodic boundary conditions. The initial condition is c0 x, y sin 2 x y . The exact solution is given by c x, y, t exp 8 2 t sin 2 x y 2t . The numerical convergence rate is evaluated by estimating the errors at five different grid levels. For the coarsest and the finest meshes, the horizontal domain is discretized by 126 and 32,766 triangle elements, respectively. The vertical domain is discretized into three layers for all five meshes. We take ε = 0.0001 and the final time t = 0.1, with time step Δt = 0.005 at the coarsest mesh and halved at the next finer level. The results of the CDG scheme (see Table 1) show that the second order of accuracy is achieved in the proposed model. The model predictions are found to be in very close agreement with the exact solution. Table 1
Two dimensional accuracy test for CDG CDG
Horizontal elements
2
L error
Convergence rate
126
1.91e-1
---
510
4.00e-2
2.26
2046
9.33e-3
2.10
8190
2.19e-3
2.09
32,766
5.02e-4
2.13
4.2 Advection of a square-shaped scalar field In this section, the test case is the advection of a square-shaped scalar field which was studied by Tamamidis and Assanis (1993). This case is selected because it can test the ability to capture the discontinuity of the numerical schemes. In this case, the two-dimensional equation c (uc) (vc) 0 t x y
(44) is solved in a three-dimensional grid system. The domain is defined on ( x, y, z) [0, 6]×[0, 6] ×[0, 6] with velocity u = v = 1.0, while the initial size of the square scalar field is [0.5, 2]×[0.5, 2] on the horizontal plane. The 3D view of the initial field is plotted in Fig. 3. The initial value in the square field is 10.0, and it is 0 in the other areas. The field is advected a distance of 3.95 from its initial location. The horizontal domain is discretized by 7198 triangle elements, and the vertical domain is discretized into three layers. With the irregular mesh, we test the robustness of the scheme presented here against the finite volume method developed by Chen et al. (2003). The maximum and minimum values and the root mean square (RMS) error of the scalar field after the advection are presented in Table 2. The 3D views of the concentration fields are shown in Figs. 4(a)-4(g). The concentration profiles at x = 5 are shown in Fig. 5. From these results, the finite volume method is found to destroy the accuracy (Fig. 4(a)), and the result from the CDG without the slope limiter produces spurious oscillations around discontinuities (Fig. 4(b)). The results from the CDG with the modified minmod limiter are sensitive to the parameter M, i.e., too small a value of M implies higher local dissipation and order reduction (Fig. 4(c)), whereas too high a value of M reintroduces the oscillations (Fig. 4(e) and Fig. 4(f)). Although the solution from the CDG with the modified minmod limiter (M = 5) noticeably improves the results and more exactly captures the maximum and minimum field values than the former, it introduces heavier numerical diffusion at sharp fronts (Fig. 4(d)) and exhibits a larger RMS error. In contrast, the solution from the CDG with the vertex-based limiter achieves rather satisfactory results in maintaining the field gradient and exactly capturing the maximum and minimum field values (Fig. 4(g)) and exhibits the smallest RMS error. The conclusion can be drawn that the CDG scheme with the vertex-based limiter is the most accurate and performs well for solutions with sharp fronts in this test problem. The results indicated that the proposed model has the ability to capture a discontinuity and that the vertex-based limiter should be the first choice for the suspended sediment model. Table 2
Performance of the different schemes Maximum scalar value
Minimum scalar value
RMS error
Finite volume method
10.0
0.0
0.171
CDG without slope limiter
10.904
-0.576
0.131
CDG with modified minmod limiter (M =
10.077
-0.0649
0.144
5) CDG with vertex-based limiter
Fig. 3
10.0
3D view of the initial scalar field
0.0
0.095
(a) Finite volume method.
(b) CDG method without the slope limiter.
(c) CDG method with the modified minmod limiter (M = 1).
(d) CDG method with the modified minmod limiter (M = 5).
(e) CDG method with the modified minmod limiter (M = 10).
(f) CDG method with the modified minmod limiter (M = 50).
(g) CDG method with the vertex-based limiter. Fig. 4 3D view plots of the advected square scalar field
Fig. 5 Comparisons of the results given by different numerical methods at x = 5 (─ Finite volume method, ○ CDG method without the slope limiter, □ CDG method with the modified minmod limiter (M = 5), △ CDG method with the vertex-based limiter)
5 A numerical example for high-gradient suspended sediment transport In a coastal area, a strong, narrow and channelized seaward flow, known as rip current, may occur at certain conditions. The rip current, which exhibits concentrated sediment in the rip channel, may be responsible for significant offshore transport of sediment (Thorpe et al., 2013). In this section, the developed model will be used for the description of high-gradient suspended sediment transport in a field-scale rip current. The domain is defined on ( x, y) [0,30km] [0,1km] , the bed slope is set to 1:2000, and the minimum water depth is 2.0 m at x = 0. The horizontal domain is discretized with an unstructured mesh of 6475 nodes and 12,176 elements. The height of each vertical layer is 2.0 m. The flow is one-dimensional and remains constant (flow per unit width is 4.0 m2/s). The initial condition is zero concentration. The boundary conditions of sediment for the problem are of Dirichlet type at the inflow with 0, c(0, y, z ) 10kg/m3 , 0,
y [0,250m) y [250m,750m] y (750m,1000m]
(45) To form the high-gradient concentration front, the advection-dominated sediment transport (εx =εy = εz = 0.02 m2/s) is considered, and the deposition and erosion processes of the sediment are ignored. The concentration solutions and the concentration profiles at y = 500 m of this problem are shown at times t = 3600 s, 7200 s, 10,800 s and 14,400 s (Fig. 6 and Fig. 7). Note that the concentration front is well described and that the
proposed model maintains the concentration gradient and captures the maximum and minimum field values at different times.
Fig. 6 Simulation of the advection-dominated sediment transport at times t = 3600 s, 7200 s, 10,800 s and 14,400 s (top to bottom) Fig. 7 The concentration profiles at times t = 3600 s (─), 7200 s (○), 10,800 s (□) and 14,400 s (△)
6 Conclusions In the present paper, a three-dimensional model was developed to simulate the transport characteristics of the suspended sediment in situations with large horizontal gradients (discontinuous situations) based on the compact discontinuous Galerkin method for unstructured mesh systems. In the model, the slope limiter is introduced to eliminate numerical oscillations around discontinuities of the suspended sediment. The evaluation of the performance of the model against theory revealed good agreement. The model is further compared to the commonly used algorithms, such as the finite volume method (FVM). Compared to the true results, the RMS (root mean square) from the model introduced in this paper is approximately half of that obtained with the FVM. The model is thus able to achieve second-order accuracy. The model was also applied to an idealized field-scale case: a discontinuous distribution of suspended sediment caused by rip currents. The results of the application of the model indicated that the model is capable of capturing the high-gradient features of the suspended sediment. Acknowledgements The study was supported by the Natural Science Foundation of Tianjin, China (Grant No. 12JCZDJC30200), the Research Fund for the Doctoral Program of Higher Education of China (RFDP) (Grant No. 20120032130001) and the Science Fund for Creative Research Groups of the National Natural Science Foundation of China (Grant No. 51321065). The authors gratefully acknowledge the support of High Performance Computing Center of Tianjin University. We would like to thank Prof. Christopher C Pain of Imperial College in UK for providing software FLUIDUTY, which played a guide role in developing our model. References Ariathurai, R., & Arulanandan, K. (1978). Erosion rates of cohesive soils. Journal of the Hydraulics Division, ASCE, 104(2), 279-283. Arnold, D. N., Brezzi, F., Cockburn, B., & Marini, L. D. (2002). Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM Journal on Numerical Analysis, 39(5), 1749-1779. Burbeau, A., Sagaut, P., & Bruneau, C. H. (2001). A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods. Journal of Computational Physics, 169(1), 111-150. Chen, C. H., Liu, R. C., & Beardsley, R. C. (2003). An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: application to coastal ocean and estuaries. Journal of Atmospheric and Oceanic Technology, 20(1), 159-186. Cockburn, B., & Shu, C. W. (1998a). The Runge-Kutta discontinuous Galerkin method for conservation laws V: Multidimensional systems. Journal of Computational Physics, 141(2), 199-224. Cockburn, B., & Shu, C. W. (1998b). The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM Journal on Numerical Analysis, 35(6), 2440-2463. Cockburn, B., & Shu, C. W. (2001). Runge-Kutta discontinuous Galerkin methods for convection-dominated problems. Journal of Scientific Computing, 16(3), 173-261. Deltares. (2010). Delft3D-FLOW: Simulation of multi-dimensional hydrodynamic flows and transport phenomena, including sediments-User Manual(Version 3.04). Delft, The Netherlands. Dolejší, V., Holík, M., & Hozman, J. (2011). Efficient solution strategy for the semi-implicit discontinuous Galerkin discretization of the Navier-Stokes equations. Journal of Computational Physics, 230(11), 4176-4200.
Gómez, H., Colominas, I., Navarrina, F., & Casteleiro, M. (2007). A discontinuous Galerkin method for a hyperbolic model for convection-diffusion problems in CFD. International Journal for Numerical Methods in Engineering, 71(11), 1342-1364. Gourgue, O. (2011). Finite element modeling of sediment dynamics in the Scheldt (Doctora dissertation), Université Catholique de Louvain. Hamrick, J. M. (1992). A three dimensional environmental fluid dynamics computer code: Theoretical and computational aspects. Special report on marine science and ocean engineering No.317. Virginia Institute of Marine Science, The college of William and Mary, Virginia. Hartmann, R., & Houston, P. (2006). Symmetric interior penalty DG methods for the compressible Navier-Stokes equations I: Method formulation. International Journal of Numerical Analysis and Modeling, 3(1), 1-20. Khosronejad, A., Rennie, C. D., Neyshabouri, S. A. A. S., & Townsend, R. D. (2007). 3D numerical modeling of flow and sediment transport in laboratory channel bends. Journal of Hydraulic Engineering, ASCE, 133(10), 1123-1134. Krivodonova, L., Xin, J., Remacle, J. F., Chevaugeon, N., & Flaherty, J. E. (2004). Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws. Applied Numerical Mathematics, 48(34), 323-338. Kuzmin, D. (2010). A vertex-based hierarchical slope limiter for p-adaptive discontinuous Galerkin methods. Journal of Computational and Applied Mathematics, 233(12), 3077-3085. Lin, B. L., & Falconer, R. A. (1996). Numerical modelling of three-dimensional suspended sediment for estuarine and coastal waters. Journal of Hydraulic Research, 34(4), 435-456. Peraire, J., & Persson P. O. (2008). The compact discontinuous Galerkin (CDG) method for elliptic problems. SIAM Journal of Scientific Computing, 30(4), 1806-1824. Reed, W. H., & Hill T. (1973). Triangular mesh method for the neutron transport equation. Los Alamos Report La, UR-73-479. Retrieved from http://www.osti.gov/scitech/servlets/purl/4491151 Seby, P. K. (1984). High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM Journal on Numerical Analysis, 21(5), 995-1011. Tamamidis, P., & Assanis, D. N. (1993). Evaluation of various high-order-accuracy schemes with and without flux limiters. International Journal for Numerical Methods in Fluids, 16(10), 931-948. Thorpe, A., Miles, J., Masselink, G., Russell, P., Scott, T., & Austin, M. (2013). Suspended sediment transport in rip currents on a macrotidal beach. Journal of Coastal Research, 118(3), 1880-1885. Rijn, L. C. V. (1984). Sediment transport, part II: suspended load transport. Journal of Hydraulic Engineering, 110(11), 1613-1641. Wang, S. Y., & Adeff, S. E. (1986). Three-dimensional modeling of river sedimentation process. Proceedings of the Third International Symposium on River Sedimentation, University of Mississippi, Mississippi. 3, 1496-1505. Warner, J. C., Sherwood, C. R., Signell, R. P., Harris, C. K., & Arango, H. G. (2008). Development of a three-dimensional, regional, coupled wave, current, and sediment-transport model. Computers and Geosciences, 34(10), 1284-1306. Wu, W., Rodi, W., & Wenka, T. (2000). 3D numerical modeling of flow and sediment transport in open channels. Journal of Hydraulic Engineering, ASCE, 126(1), 4-15. Younes, A., & Ackerer, P. (2008). Solving the advection-dispersion equation with discontinuous Galerkin and multipoint flux approximation methods on unstructured meshes. International Journal for Numerical Methods in Fluids, 58(6), 687-708.
Zhang, Y., Zhang, X., & Shu, C. W. (2013). Maximum-principle-satisfying second order discontinuous Galerkin schemes for convection-diffusion equations on triangular meshes. Journal of Computational Physics, 234(2), 295-316.
Table 1. Two dimensional accuracy test for CDG CDG
Horizontal elements
L2 error
Convergence rate
126
1.91e-1
---
510
4.00e-2
2.26
2046
9.33e-3
2.10
8190
2.19e-3
2.09
32766
5.02e-4
2.13
Table 2. Performance of the different schemes Maximum
Minimum
RMS
scalar value
scalar value
error
Finite volume method
10.0
0.0
0.171
CDG without slope limiter
10.904
-0.576
0.131
10.077
-0.0649
0.144
10.0
0.0
0.095
CDG with modified minmod limiter (M=5) CDG with vertex-based limiter
A three-dimensional model for suspended sediment transport based on the compact discontinuous Galerkin method Zhao Zhang-Yi, Zhang Qing-He, Zhao Hong-Bo, Yang Hua
www.elsevier.com/locate/ijsrc
PII: DOI: Reference:
S1001-6279(15)00043-8 http://dx.doi.org/10.1016/j.ijsrc.2014.06.001 IJSRC31
To appear in:
International Journal of Sediment Research
Received date: 4 February 2013 Revised date: 17 June 2014 Accepted date: 29 June 2014 Cite this article as: Zhao Zhang-Yi, Zhang Qing-He, Zhao Hong-Bo, Yang Hua, A three-dimensional model for suspended sediment transport based on the compact discontinuous Galerkin method, International Journal of Sediment Research, http://dx.doi.org/10.1016/j.ijsrc.2014.06.001 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
A three-dimensional model for suspended sediment transport based on the compact discontinuous Galerkin method Zhao Zhang-Yi1, 2, Zhang Qing-He1, Zhao Hong-Bo2, Yang Hua2 1
State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin
University, Tianjin 300072, China 2
Key Laboratory of Engineering Sediment, Ministry of Transport, Tianjin Research
Institute for Water Transport Engineering, M. O. T., Tianjin 300456, China Corresponding author: Zhang Qing-He Tel: +86-22-27404266; fax: +86-22-27404266. E-mail address: [email protected]
A three-dimensional model for suspended sediment transport based on the compact discontinuous Galerkin method Zhao Zhangyi a, b, Zhang Qinghe a,, Zhao Hong-Bo b, Yang Hua b a
State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University, Tianjin 300072, China
b
Key Laboratory of Engineering Sediment, Ministry of Transport, Tianjin Research Institute for Water Transport Engineering, M. O. T., Tianjin 300456, China Abstract A three-dimensional numerical model was developed in the present paper, with the goal of simulating suspended sediment transport in discontinuous situations. In the model, the advection-diffusion equation of suspended sediment was solved using the Compact Discontinuous Galerkin (CDG) method based on unstructured meshes. The evaluation of the model revealed good agreement between the proposed model and existing theories. The model has second-order accuracy. Moreover, the model was applied to an idealized case of suspended sediment transport caused by rip currents based on the field spatial scale. Key Words: advection-diffusion equation; model evaluation; discontinuity; rip currents; unstructured meshes; slope limiter
1 Introduction Numerical modeling has become a powerful tool in the contemporary study of sediment transport in rivers, estuaries and coastal areas. Based on various purposes, a number of related studies have been completed. For example, Wang and Adeff (1986) proposed a 3-D model for suspended sediment transport based on the finite element method to study river sedimentation. Lin and Falconer (1996) established a layer-integrated 3D model for the estuarine and coastal waters around the UK using a finite difference method. Wu et al. (2000) developed a 3-D model to investigate the transport properties of both the suspended sediment and the bed load based on the finite volume method. Khosronejad et al. (2007) developed a 3-D model for sediment transport in curved open channels based on the finite volume method. Recently, the suspended module of sediment transport has become a regular component of an environmental model, including both commercial models and community models, e.g., EFDC (Environmental Fluid Dynamics Code) (Hamrick, 1992), FVCOM (Finite Volume Coastal Ocean Model) (Chen et al., 2003), ROMS (Regional Ocean Modeling System) (Warner et al., 2008), and Delft-3D (Deltares, 2010). Most of the above-mentioned models are only suitable for continuous
Corresponding author. Email address: [email protected]
sediment transport. However, in discontinuous or large gradient flows, such as dam-breaking flow and tidal bore flow, the discontinuous transportation of sediment may exist. In this case, standard numerical methods, such as the finite element or finite volume method, may generate solutions with numerical diffusion and/or non-physical oscillations (Younes & Ackerer 2008). Therefore, there is practical-significance to develop a high-resolution suspended sediment transport model using a discontinuity-capturing method; a candidate for such a model is the Discontinuous Galerkin (DG) method (Reed & Hill, 1973). Currently, in addition to Gourgue’s work (2011) on depth-averaged 2-D model applied in the Scheldt estuary, there are only a few published papers regarding the suspended sediment transport model based on the DG method although there have been some research reports on the application of the DG method for advection-diffusion problems (Cockburn & Shu, 1998b; Gómez et al., 2007; Zhang et al., 2013). Therefore, we will develop a 3-D suspended sediment transport model based on the DG method with unstructured meshes. The DG method used in this study is based on the so-called compact discontinuous Galerkin (CDG) method, which was first presented by Peraire and Persson (2008). This method is an extension of the Runge-Kutta DG method (Cockburn & Shu, 1998a; 2001) and is closely related to the local discontinuous Galerkin (LDG) method (Cockburn & Shu, 1998b). Peraire and Persson (2008) proved that the CDG method inherits all of the attractive features of the LDG method. In addition, the CDG method is more compact and stable than the LDG method in situations involving multiple dimensions. The outline of the paper is as follows. The mathematical formulation of 3-D suspended sediment transport is briefly introduced in Section 2. The numerical scheme and the spatial discretization of the CDG method as well as its temporal discretization are described in Section 3. Some numerical test calculations are presented for evaluating the performance of the model in Section 4. The results of a numerical example for suspended sediment transport are presented to demonstrate the applicability of the proposed model in Section 5. Finally, some conclusions are presented in Section 6. 2 Mathematical formulation of suspended sediment transport 2.1 Governing equation The distribution of the suspended sediment concentration in water is governed by the following advection-diffusion equation: c t (U c)
K ( c) 0
(1) where c is the suspended sediment concentration, ( x , y , z) , U = (ux, uy, uz-ωs) , ux, uy and T
uz are the velocity components, ωs is the settling velocity of suspended sediment, and K ( x
0 0), (0 y
0), (0 0 z )
T
where εx, εy and εz are the mixing coefficients of the sediment in the x-, y- and z-directions, respectively. In applying Eq. (1) to rivers, estuaries and coastal areas where the ratio of the vertical length scale to the horizontal length scale is small, the vertical diffusivity terms are usually several orders of magnitude larger than the horizontal diffusivity terms (Lin & Falconer, 1996). Therefore, it is more important to exactly calculate the vertical mixing coefficients than the horizontal ones. In the present paper, the horizontal mixing coefficients are assumed to be constant, while the vertical mixing coefficient was prescribed using the method given by van Rijn (1984).
2.2 Boundary conditions Eq. (1) is solved with the following boundary conditions: at the free surface, the vertical sediment flux is zero, and hence, a no-flux boundary condition is used for the sediment concentration: z
c s c 0 z
(2) At the bottom boundary of the suspended sediment layer, the sediment flux is the difference between the rates of deposition and erosion: z
c s c Eb Db z
(3)
The deposition rate at the interface is Db sc
(4) The erosion rate is calculated as described by Ariathurai and Arulanandan (1978): Eb E0 (1 Pb ) b c 1
(5) where E0 is the erosion flux parameter, Pb is the bottom porosity, τb is the bottom shear stress, and τc is the critical shear stress for erosion. 3 Numerical discretization of the governing equation by the DG method 3.1 Definition The following definition is used in the proposed method c g D , on D
(6a) Kc n g N , on N
(6b) where Ω is a bounded domain with boundary
D
, ∂Ω N
D
N
and ∂Ω are the Dirichlet and
Neumann boundaries, respectively, gD and gN are specified values, and n is the outward unit normal to the boundary of Ω. 3.2 The spatial discretization For the sake of simplicity, Eq. (1) is rewritten as the two respective parts to be discretized. c t (Uc) f
(7a) (Kc) f
(7b) 3.2.1 Discretization of the advection terms The bounded domain Ω is partitioned into N non-overlapping elements. To apply the DG method to Eq. (7a), the equation is first multiplied by the trial or test function φ. The resulting equation is integrated over an element e with the face ∂e, and after integration by parts, the equation can be written as
c
( t ( U)c)de
e
e
ˆ cˆd(e) f de n U e
(8) where the hatted terms denote the numerical fluxes across the element faces. Here, the averaged velocity on each side of the face is an approximation to Uˆ , and cˆ is taken to be the upwind value (Fig. 1). If we use cupw to denote the value of c on the upwind side of the face, then by summing Eq. (8) over all elements and considering Eq. (6), we have the following global expression:
(
c ˆ c)d U D n Uˆ ( gD c)d() s n Uˆ (cupw c)ds f d t
(9) where s is the union of the interior faces of all elements, s- is the inflow part of s, and
D
and D
are the outflow and inflow Dirichlet boundaries, respectively.
(a) when the flow is outflow.
(b) when the flow is inflow. Fig. 1
Illustration of the upwind value
3.2.2 Discretization of the diffusion terms Following Cockburn and Shu (1998a), we introduce the new variable ξ and rewrite Eq. (7b) as a first order system of equations ξ f
(10a) ξ Kc
(10b) Introducing a scalar test function φ and a vector test function λ, multiplying Eqs. (10a) and (10b) by the respective test functions, integrating over an element e, and then integrating by parts, we obtain ξ d e n ξˆd(e) f de e
e
e
(11) λKcde n λKcˆd(e) λ ξde e
e
e
(12) where the numerical fluxes ξˆ and cˆ are approximations to ξ and c, respectively, on the face of the element. Following the notation of Arnold et al. (2002), we consider two neighboring elements e+ and e- with the common face k e
e , and assume that n denotes the unit normal vectors to the element
boundary e . Let λ and denote the test functions on k. The average and jump operators are defined as follows
{λ} (λ λ ) 2, {} ( ) 2
(13a) [λ] λ n λ n , [ ] n n
(13b) By adding Eqs. (11) and (12) over all elements, the following global equations are obtained
ξ d f d ξˆ [ ]ds ξˆ nd()
s
(14)
ξ λd Kc λd Kcˆ[λ ]ds Kcˆλ nd()
s
(15) Next, after performing integration by parts on the first term on the right side of Eq. (15), we obtain
ξ λd λ (Kc)d K([c] λ {cˆ c}[λ])ds K(cˆ c)λ nd()
s
(16) The numerical fluxes ξˆ and cˆ for the interior faces in Eqs. (15) and (16) are given by ξˆ {ξ e} C11[c] C12[ξ e ]
(17) cˆ {c} C12 [c]
(18) and the numerical fluxes for the boundary faces are ξˆ ξe C11 (c g D )n, cˆ g D , on D
(19a) ξˆ g N n, cˆ c, on N
(19b) e
e
where C11 is a positive constant, C12 is a vector, and ξ c ξ . Making use of Eqs. (18) and (19), Eq. (16) can be written as
ξ λd λ (Kc)d K([c] {λ} C12 [c][λ])ds
D
s
K( g D c)λ nd()
(20) To obtain an expression for ξ as a function of c, we follow Peraire and Persson (2008) and introduce the local lifting operators:
r e (ξ) λ ξ {λ}, e s e
(21a)
r e (c) λ c[λ ], e s e
(21b)
r (c) λ cλ n, e D
e D
e
(21c)
Setting λ in Eq. (20), Eq. (14) can be rewritten as
(Kc)d K([c] {} C12 [c][ ])ds
D
s
K( g D c) nd() f d ξ [ ]ds ξ nd()
s
(22) According to Eqs. (17), (19) and (21), the terms involving ξˆ in Eq. (22) can be written as
ξˆ [ ]ds (({c} C s
12
s
[c])[ ] C11[c] [ ])ds
(r e ([ ]) r e (C12 [ ])) (r e ([c]) r e (C12 [c]) rDe (c))d (r e ([ ]) r e (C12 [ ])) rDe ( g D )d es
es
(23)
ξˆ nd()
D
(c n C11 (c g D ))d()
eD
eD
e
rDe ( g D ) nd(e)
N
g N d()
r ( ) (r e ([c]) r e (C12 [c]) rDe (c))d()
e D
(24) Therefore, Eqs. (9), (22)~(24) comprise the weak formulation of the advection-diffusion equation. Because the second-order accurate numerical model will be developed as described by Cockburn and Shu (2001), the DG method with pth (p = 0, 1, 2...) order polynomials is always (p+1)th order accurate. The linear solution c on a tetrahedron element e in the weak formulation of the advection-diffusion equation can be expressed in terms of local basis functions i and degrees of freedom ci: Nloc
c cii i 1
(25) where Nloc is the number of vertices of the element. 3.3 Temporal discretization Making use of Eq. (25), the weak formulation of the advection-diffusion equation produces a semi-discrete equation that can be written in matrix form as follows M
dc A(U)c Dc r dt
(26) where M is the mass matrix, A(U) is the advection matrix, D is the diffusion matrix and r is the right-hand side vector containing boundary, source and absorption terms. When solving the 3D advection-diffusion equation for realistic problems, it is advantageous to use a large simulation time step to achieve fast CPU times, and it is suitable to use an implicit scheme without a strong restriction of the size of the time step (Hartmann & Houston, 2005). However, a fully implicit scheme leads to nonlinear algebraic equations at each time level, which is rather expensive to solve (Dolejší et al., 2011). Therefore, only the advection terms in the advection-diffusion equation are treated explicitly, while the diffusion terms are treated implicitly in the proposed model. To achieve this, Eq. (26) is discretized in two steps as follows: M c* cn t A(Un )cn rDn
(27)
M cn 1 c* t Dc n 1 rNn 1 rsn 1
(28) where c* is the preliminary value, rD, rN, and rs are the Dirichlet, Neumann boundary and source terms, respectively. 3.4 Slope limiter The high-order method often produces spurious oscillations and instability in the vicinity of discontinuities (Burbeau et al., 2001). To ensure the stability of the high-order scheme and eliminate the non-physical oscillation, a slope limiter, which can limit the solution gradient near discontinuities, is usually used. In recent years, a number of slope limiting techniques are available for the DG method (Cockburn & Shu, 1998b; Burbeau et al., 2001; Krivodonova et al, 2004; Kuzmin, 2010). However, the slope limiter may introduce a certain amount of numerical diffusion while suppressing oscillation and often requires parameter adjustment (Kuzmin, 2010). The modified minmod limiter and the vertex-based limiter, with only one parameter or no user-adjusted parameters, will be compared in Section 4. A brief introduction regarding the two slope limiters for the second-order DG method will first be presented here. 3.4.1 Modified minmod limiter The modified minmod limiter (Cockburn & Shu, 1998b) is defined as follows: m(a1 , a2 ,
a1 , if a1 M x 2 , an ) otherwise m(a1 , a2 , , an ),
(29) where M is a positive number, and m is the minmod limiter function (Seby, 1984). More detailed discussion about how sensitive the numerical model results are to the parameter M will be provided in Section 4. m(a1 , a2 ,
, an ) s min ai if sgn(a1 ) 1 i n
sgn(an ) s
(30) 1,
x0
where sgn( x) 0, x 0 . 1, x 0
To construct the modified slope limiter for tetrahedral elements, we proceed as follows. Considering the tetrahedrons in Fig. 2, where m1 is the barycenter of the face of e0, and bi denotes the barycenter of the tetrahedron ei (i = 0, 1, 2, 3, 4).
Fig. 2
The tetrahedron element e0 and its four neighbors
Because we have 4
m1 b0 i (bi b0 ) i 1
(31) where αi (i = 1, 2, 3, 4) depends on m1 and the geometry, we can obtain for any linear function c 4
c(m1 ) c(b0 ) i (c(bi ) c(b0 )) i 1
(32) and cei
cde
ei c(bi ), i 0,1,2,3,
ei
(33) we have 4
c (m1 , e0 ) c(m1 ) ce0 i (cei ce0 ) c (m1, e0 ) i 1
(34) Supposing c is a piecewise linear function, and let mi (i = 1, 2, 3, 4) be the four barycenters of the faces of the e0, we obtain 4
4
i 1
i 1
c( x, y, z ) c(mi ) i ( x, y, z ) ce0 c(mi , e0 ) i ( x, y, z), ( x, y, z) e0
(35) where βi (i = 1, 2, 3, 4) depends on the geometry. To use the slope limiter for solution, the following quantities are computed first: i m(c(mi , e0 ), c (mi , e0 )) (36) where m is the modified minmod limiter defined in Eq. (29), and υ is the limiter factor (Cockburn & Shu, 1998a), taking υ = 1.1 in our model. If
4
i 1
i
0 , then the limiter is only applied to c* in Eq. (27), and we have 4
c* ( x, y, z ) ce0 i i ( x, y, z ), ( x, y, z ) e0 i 1
(37) If
4
i 1
i
0 , then the following variables are computed: 4
4
i 1
i 1
pos= max(0, i ), neg= max(0, i )
(38) + = min 1,neg pos , min 1,pos neg
(39) Next, the limiter is applied to c* in Eq. (27), and we have
4
c* ( x, y, z ) ce0 ˆ i i ( x, y, z ), ( x, y, z ) e0 i 1
(40) where ˆ i max(0, i ) max(0, i ) .
3.4.2 Vertex-based limiter The vertex-based limiter was introduced by Kuzmin (2010) and is defined as
min 1, cimax ce ci ce , ci ce 0 ve min 1, ci ce 0 i min min 1, ci ce ci ce , ci ce 0
where c
max i
max(ce , c
max i
), c
min i
(41) min(ce , c ) , ci is the value on the node of the tetrahedron e, and ce min i
denotes the value on the barycenter of the tetrahedral element e. Note that cimax and cimin are initialized by small and large constants, respectively. The vertex-based limiter is only applied to the DG field c* in Eq. (27). For linear elements, c* in any element may be written as c* ( x, y, z) ce ve (ci ce ), ( x, y, z) e
(42) 4 Numerical tests 4.1 Accuracy test The accuracy of the second-order CDG scheme is tested by solving c (uc) (vc) 2 2 2 ( uc) 2 ( vc) 0 t x y x y
(43) in a three-dimensional grid system. The domain is defined on ( x, y, z) [0, 1]×[0, 1] ×[0, 1] with velocity u = v = 1.0 and periodic boundary conditions. The initial condition is c0 x, y sin 2 x y . The exact solution is given by c x, y, t exp 8 2 t sin 2 x y 2t . The numerical convergence rate is evaluated by estimating the errors at five different grid levels. For the coarsest and the finest meshes, the horizontal domain is discretized by 126 and 32,766 triangle elements, respectively. The vertical domain is discretized into three layers for all five meshes. We take ε = 0.0001 and the final time t = 0.1, with time step Δt = 0.005 at the coarsest mesh and halved at the next finer level. The results of the CDG scheme (see Table 1) show that the second order of accuracy is achieved in the proposed model. The model predictions are found to be in very close agreement with the exact solution. Table 1
Two dimensional accuracy test for CDG CDG
Horizontal elements
2
L error
Convergence rate
126
1.91e-1
---
510
4.00e-2
2.26
2046
9.33e-3
2.10
8190
2.19e-3
2.09
32,766
5.02e-4
2.13
4.2 Advection of a square-shaped scalar field In this section, the test case is the advection of a square-shaped scalar field which was studied by Tamamidis and Assanis (1993). This case is selected because it can test the ability to capture the discontinuity of the numerical schemes. In this case, the two-dimensional equation c (uc) (vc) 0 t x y
(44) is solved in a three-dimensional grid system. The domain is defined on ( x, y, z) [0, 6]×[0, 6] ×[0, 6] with velocity u = v = 1.0, while the initial size of the square scalar field is [0.5, 2]×[0.5, 2] on the horizontal plane. The 3D view of the initial field is plotted in Fig. 3. The initial value in the square field is 10.0, and it is 0 in the other areas. The field is advected a distance of 3.95 from its initial location. The horizontal domain is discretized by 7198 triangle elements, and the vertical domain is discretized into three layers. With the irregular mesh, we test the robustness of the scheme presented here against the finite volume method developed by Chen et al. (2003). The maximum and minimum values and the root mean square (RMS) error of the scalar field after the advection are presented in Table 2. The 3D views of the concentration fields are shown in Figs. 4(a)-4(g). The concentration profiles at x = 5 are shown in Fig. 5. From these results, the finite volume method is found to destroy the accuracy (Fig. 4(a)), and the result from the CDG without the slope limiter produces spurious oscillations around discontinuities (Fig. 4(b)). The results from the CDG with the modified minmod limiter are sensitive to the parameter M, i.e., too small a value of M implies higher local dissipation and order reduction (Fig. 4(c)), whereas too high a value of M reintroduces the oscillations (Fig. 4(e) and Fig. 4(f)). Although the solution from the CDG with the modified minmod limiter (M = 5) noticeably improves the results and more exactly captures the maximum and minimum field values than the former, it introduces heavier numerical diffusion at sharp fronts (Fig. 4(d)) and exhibits a larger RMS error. In contrast, the solution from the CDG with the vertex-based limiter achieves rather satisfactory results in maintaining the field gradient and exactly capturing the maximum and minimum field values (Fig. 4(g)) and exhibits the smallest RMS error. The conclusion can be drawn that the CDG scheme with the vertex-based limiter is the most accurate and performs well for solutions with sharp fronts in this test problem. The results indicated that the proposed model has the ability to capture a discontinuity and that the vertex-based limiter should be the first choice for the suspended sediment model. Table 2
Performance of the different schemes Maximum scalar value
Minimum scalar value
RMS error
Finite volume method
10.0
0.0
0.171
CDG without slope limiter
10.904
-0.576
0.131
CDG with modified minmod limiter (M =
10.077
-0.0649
0.144
5) CDG with vertex-based limiter
Fig. 3
10.0
3D view of the initial scalar field
0.0
0.095
(a) Finite volume method.
(b) CDG method without the slope limiter.
(c) CDG method with the modified minmod limiter (M = 1).
(d) CDG method with the modified minmod limiter (M = 5).
(e) CDG method with the modified minmod limiter (M = 10).
(f) CDG method with the modified minmod limiter (M = 50).
(g) CDG method with the vertex-based limiter. Fig. 4 3D view plots of the advected square scalar field
Fig. 5 Comparisons of the results given by different numerical methods at x = 5 (─ Finite volume method, ○ CDG method without the slope limiter, □ CDG method with the modified minmod limiter (M = 5), △ CDG method with the vertex-based limiter)
5 A numerical example for high-gradient suspended sediment transport In a coastal area, a strong, narrow and channelized seaward flow, known as rip current, may occur at certain conditions. The rip current, which exhibits concentrated sediment in the rip channel, may be responsible for significant offshore transport of sediment (Thorpe et al., 2013). In this section, the developed model will be used for the description of high-gradient suspended sediment transport in a field-scale rip current. The domain is defined on ( x, y) [0,30km] [0,1km] , the bed slope is set to 1:2000, and the minimum water depth is 2.0 m at x = 0. The horizontal domain is discretized with an unstructured mesh of 6475 nodes and 12,176 elements. The height of each vertical layer is 2.0 m. The flow is one-dimensional and remains constant (flow per unit width is 4.0 m2/s). The initial condition is zero concentration. The boundary conditions of sediment for the problem are of Dirichlet type at the inflow with 0, c(0, y, z ) 10kg/m3 , 0,
y [0,250m) y [250m,750m] y (750m,1000m]
(45) To form the high-gradient concentration front, the advection-dominated sediment transport (εx =εy = εz = 0.02 m2/s) is considered, and the deposition and erosion processes of the sediment are ignored. The concentration solutions and the concentration profiles at y = 500 m of this problem are shown at times t = 3600 s, 7200 s, 10,800 s and 14,400 s (Fig. 6 and Fig. 7). Note that the concentration front is well described and that the
proposed model maintains the concentration gradient and captures the maximum and minimum field values at different times.
Fig. 6 Simulation of the advection-dominated sediment transport at times t = 3600 s, 7200 s, 10,800 s and 14,400 s (top to bottom) Fig. 7 The concentration profiles at times t = 3600 s (─), 7200 s (○), 10,800 s (□) and 14,400 s (△)
6 Conclusions In the present paper, a three-dimensional model was developed to simulate the transport characteristics of the suspended sediment in situations with large horizontal gradients (discontinuous situations) based on the compact discontinuous Galerkin method for unstructured mesh systems. In the model, the slope limiter is introduced to eliminate numerical oscillations around discontinuities of the suspended sediment. The evaluation of the performance of the model against theory revealed good agreement. The model is further compared to the commonly used algorithms, such as the finite volume method (FVM). Compared to the true results, the RMS (root mean square) from the model introduced in this paper is approximately half of that obtained with the FVM. The model is thus able to achieve second-order accuracy. The model was also applied to an idealized field-scale case: a discontinuous distribution of suspended sediment caused by rip currents. The results of the application of the model indicated that the model is capable of capturing the high-gradient features of the suspended sediment. Acknowledgements The study was supported by the Natural Science Foundation of Tianjin, China (Grant No. 12JCZDJC30200), the Research Fund for the Doctoral Program of Higher Education of China (RFDP) (Grant No. 20120032130001) and the Science Fund for Creative Research Groups of the National Natural Science Foundation of China (Grant No. 51321065). The authors gratefully acknowledge the support of High Performance Computing Center of Tianjin University. We would like to thank Prof. Christopher C Pain of Imperial College in UK for providing software FLUIDUTY, which played a guide role in developing our model. References Ariathurai, R., & Arulanandan, K. (1978). Erosion rates of cohesive soils. Journal of the Hydraulics Division, ASCE, 104(2), 279-283. Arnold, D. N., Brezzi, F., Cockburn, B., & Marini, L. D. (2002). Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM Journal on Numerical Analysis, 39(5), 1749-1779. Burbeau, A., Sagaut, P., & Bruneau, C. H. (2001). A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods. Journal of Computational Physics, 169(1), 111-150. Chen, C. H., Liu, R. C., & Beardsley, R. C. (2003). An unstructured grid, finite-volume, three-dimensional, primitive equations ocean model: application to coastal ocean and estuaries. Journal of Atmospheric and Oceanic Technology, 20(1), 159-186. Cockburn, B., & Shu, C. W. (1998a). The Runge-Kutta discontinuous Galerkin method for conservation laws V: Multidimensional systems. Journal of Computational Physics, 141(2), 199-224. Cockburn, B., & Shu, C. W. (1998b). The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM Journal on Numerical Analysis, 35(6), 2440-2463. Cockburn, B., & Shu, C. W. (2001). Runge-Kutta discontinuous Galerkin methods for convection-dominated problems. Journal of Scientific Computing, 16(3), 173-261. Deltares. (2010). Delft3D-FLOW: Simulation of multi-dimensional hydrodynamic flows and transport phenomena, including sediments-User Manual(Version 3.04). Delft, The Netherlands. Dolejší, V., Holík, M., & Hozman, J. (2011). Efficient solution strategy for the semi-implicit discontinuous Galerkin discretization of the Navier-Stokes equations. Journal of Computational Physics, 230(11), 4176-4200.
Gómez, H., Colominas, I., Navarrina, F., & Casteleiro, M. (2007). A discontinuous Galerkin method for a hyperbolic model for convection-diffusion problems in CFD. International Journal for Numerical Methods in Engineering, 71(11), 1342-1364. Gourgue, O. (2011). Finite element modeling of sediment dynamics in the Scheldt (Doctora dissertation), Université Catholique de Louvain. Hamrick, J. M. (1992). A three dimensional environmental fluid dynamics computer code: Theoretical and computational aspects. Special report on marine science and ocean engineering No.317. Virginia Institute of Marine Science, The college of William and Mary, Virginia. Hartmann, R., & Houston, P. (2006). Symmetric interior penalty DG methods for the compressible Navier-Stokes equations I: Method formulation. International Journal of Numerical Analysis and Modeling, 3(1), 1-20. Khosronejad, A., Rennie, C. D., Neyshabouri, S. A. A. S., & Townsend, R. D. (2007). 3D numerical modeling of flow and sediment transport in laboratory channel bends. Journal of Hydraulic Engineering, ASCE, 133(10), 1123-1134. Krivodonova, L., Xin, J., Remacle, J. F., Chevaugeon, N., & Flaherty, J. E. (2004). Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws. Applied Numerical Mathematics, 48(34), 323-338. Kuzmin, D. (2010). A vertex-based hierarchical slope limiter for p-adaptive discontinuous Galerkin methods. Journal of Computational and Applied Mathematics, 233(12), 3077-3085. Lin, B. L., & Falconer, R. A. (1996). Numerical modelling of three-dimensional suspended sediment for estuarine and coastal waters. Journal of Hydraulic Research, 34(4), 435-456. Peraire, J., & Persson P. O. (2008). The compact discontinuous Galerkin (CDG) method for elliptic problems. SIAM Journal of Scientific Computing, 30(4), 1806-1824. Reed, W. H., & Hill T. (1973). Triangular mesh method for the neutron transport equation. Los Alamos Report La, UR-73-479. Retrieved from http://www.osti.gov/scitech/servlets/purl/4491151 Seby, P. K. (1984). High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM Journal on Numerical Analysis, 21(5), 995-1011. Tamamidis, P., & Assanis, D. N. (1993). Evaluation of various high-order-accuracy schemes with and without flux limiters. International Journal for Numerical Methods in Fluids, 16(10), 931-948. Thorpe, A., Miles, J., Masselink, G., Russell, P., Scott, T., & Austin, M. (2013). Suspended sediment transport in rip currents on a macrotidal beach. Journal of Coastal Research, 118(3), 1880-1885. Rijn, L. C. V. (1984). Sediment transport, part II: suspended load transport. Journal of Hydraulic Engineering, 110(11), 1613-1641. Wang, S. Y., & Adeff, S. E. (1986). Three-dimensional modeling of river sedimentation process. Proceedings of the Third International Symposium on River Sedimentation, University of Mississippi, Mississippi. 3, 1496-1505. Warner, J. C., Sherwood, C. R., Signell, R. P., Harris, C. K., & Arango, H. G. (2008). Development of a three-dimensional, regional, coupled wave, current, and sediment-transport model. Computers and Geosciences, 34(10), 1284-1306. Wu, W., Rodi, W., & Wenka, T. (2000). 3D numerical modeling of flow and sediment transport in open channels. Journal of Hydraulic Engineering, ASCE, 126(1), 4-15. Younes, A., & Ackerer, P. (2008). Solving the advection-dispersion equation with discontinuous Galerkin and multipoint flux approximation methods on unstructured meshes. International Journal for Numerical Methods in Fluids, 58(6), 687-708.
Zhang, Y., Zhang, X., & Shu, C. W. (2013). Maximum-principle-satisfying second order discontinuous Galerkin schemes for convection-diffusion equations on triangular meshes. Journal of Computational Physics, 234(2), 295-316.
Table 1. Two dimensional accuracy test for CDG CDG
Horizontal elements
L2 error
Convergence rate
126
1.91e-1
---
510
4.00e-2
2.26
2046
9.33e-3
2.10
8190
2.19e-3
2.09
32766
5.02e-4
2.13
Table 2. Performance of the different schemes Maximum
Minimum
RMS
scalar value
scalar value
error
Finite volume method
10.0
0.0
0.171
CDG without slope limiter
10.904
-0.576
0.131
10.077
-0.0649
0.144
10.0
0.0
0.095
CDG with modified minmod limiter (M=5) CDG with vertex-based limiter