Abstract
We reformulate the mathematical model for the 2D sedimentation in an estuary as a coupled nonlinear differential system. Combining the mass-conservation character of the discontinuous Galerkin method and the jump-capturing property of Lesaint-Raviart upwind technique, we design an upwind discontinuous Galerkin finite element method, which obeys the local mass conservation and possesses good stability. Our theoretical analysis shows that there exists a unique solution to the numerical procedure and the discrete solution permits convergence rate. Numerical experiments are conducted to verify our theoretical findings. This may provide a theoretical principle for better understanding of the mechanism and morphological characters of sedimentation at estuaries.
1. Introduction
This article is initiated from our project “Numerical Simulation for the Sedimentation at Yellow River’s Estuary.” It is known that the estuary of Yellow River is a weak-tidal, arenose, frequently swinging, and accumulated firth. Its sedimentation is heavily influenced by the reduced amount of water but a relatively large amount of sand from the upstream and the oceanic dynamics such as the tide, current, and wave in the littoral region. Numerical modeling and simulation for the flow and sedimentation are very important means of quantitatively predicting the formation and development of the delta and the evolution of the fluvial bed, as well as the silt deposition transport, and then providing a theoretical principle for better understanding of the mechanism and morphological characters of the sedimentation. This will be of great significance for the further development of Yellow River delta.
Based upon the literature of the sedimentation studies [1–5] for estuaries, theoretically we reformulate the mathematical model for 2D sedimentation as a combination of 4 governing equations: the flow continuity equation, the flow motion equation, the sediment continuity equation, and the bed deformation equation, in which the first two equations describe the water dynamics of the littoral region and the last two are for the dynamics of sediment deposition.
From the point of view of numerical and applicable issue, there are two important ingredients which should be strengthened in designing an ideal numerical procedure for this model. The first one is that the procedure should obey the local mass/momentum conservation law to reflect the physical characteristics of the dynamics of the flow and sediment, and the second is that the procedure should recognize the flow direction so that the procedure can capture the jumps or interfaces and maintain the good stability since the velocity of the fluid is much larger than that of the sediment deposition. For this purpose, we recall the mass-conservation character of the discontinuous Galerkin method and the jump-capturing property of upwind technique and attempt to use these merits so as to design a strongly stable and locally conservative numerical procedure for the mathematical model.
The goals of this article are as follows: () in the flow continuity equation, we reformulate the mathematical model by taking the bed height varying with time so that the flow continuity equation is naturally coupled with the sediment continuity equation and the bed deformation equation. By doing so, the mathematical model can much better describe the practical sediment process compared to the simplified case of -constant [6–8] and () by introducing the Broken Sobolev spaces and applying the discontinuous Galerkin method to discrete those space-derivative terms to design a nonsymmetric interior penalty Galerkin (NIPG) method. Considering the fact that the flow velocity is much larger than that of the sediment deposition, we pay our special attention to those terms involving by inserting Lesaint-Raviart upwind technique [9]. Therefore, an upwind discontinuous Galerkin finite element method is proposed for the reformulated mathematical model. () With the help of the induction argument, we prove that the upwind discontinuous Galerkin finite element method possesses good stability, from which the existence and uniqueness of the discrete solution follow. () Combining the proof that the discrete depth of the fluid possesses a positive lower bound, we prove that the approximate solutions for the depth of fluid , the velocity , the sediment concentration , and the bed height are of convergence rates, respectively. () Numerical experiments are conducted to verify the theoretical findings, stability, local conservation property, and convergence rates.
The achievement of these five goals constitutes the main content of the article, presented in Sections 2–6, respectively. For simplicity of presentation, we will use to denote the usual Sobolev spaces that provided the norm and the seminorm for any 2D domain . When and , the subscript and will be omitted. We also let and stand for the inner products on and on the edges , respectively, and denotes their -norms. is a generic constant and is a small positive constant, which may take different values at different places but independent of the triangulation parameters and . We define , where if , the integral is replaced with the essential supreme.
2. Reformulation of the Mathematical Model
Let be a bounded polygonal domain in , and let be a time interval. Then, the sediment transportation model at an estuary can be described by the following 4 governing equations defined over the space-time domain ; here we ignore the Coriolis force, since it does not affect the correctness of the model.
The flow continuity equation is
the flow motion equation is
the sediment continuity equation is
and the bed deformation equation iswhere the unknowns and denote the depth of water, the water level, the sediment concentration, the vertical mean velocity, and the bed height, respectively; is the sediment carrying capacity, is the gravitational acceleration, , denotes the Chezy coefficient, denotes Manning’s roughness coefficient, is the viscosity coefficient, is the sediment diffusion coefficient, is the coefficient of saturation recovery, is the settling velocity of sediment, and is the sediment dry bulk density.
In the uncoupled sedimentation model [5] and the shallow water system (SWS) [10, 11], (1) is replaced by . It means that the flow continuity equation is not coupled with the bed deformation equation; in other words, the bed height was assumed to be independent of time in the flow continuity equation. Therefore, the shallow water equations (1) and (2) and the sediment deposition models (3) and (4) can be solved separately. We remark that this assumption will lead to a simple calculation, but it is unreasonable in the real process of sediment transport. In this paper, we assume changes with time in shallow water equations, so that the flow continuity equation is much naturally coupled with the sediment continuity equation and the bed deformation equation.
Without loss of generality, the boundary of the domain is decomposed into three parts: the inflow part , the rigid boundary , and the outflow part (see Figure 1), and the initial-boundary conditions are prescribed aswhere the initial values and the boundary conditions are given data, and n denotes the unit outward normal vector to .
3. The Upwind Discontinuous Galerkin Finite Element Procedure
In this section, we shall revisit the Broken Sobolev spaces and modify those terms involving by Lesaint-Raviart upwind technique and then apply the discontinuous Galerkin method to discrete those space-derivative terms in (1)–(4) to design an upwind nonsymmetric interior penalty Galerkin (NIPG) method.
3.1. Revisiting the Broken Sobolev Spaces
We begin this subsection by recalling some notations of the triangulation for the given domain . We let be a quasi-uniform subdivision of with element ; let , denote the diameter of , and . Further we need to describe the edges of an element and the functions defined on them. We denote, by , any edge of an element , and then belongs to either its interior or boundary of . We also denote, by , the set of interior edges. If , the unit normal vector of is taken to be outward to , or if is a common edge of two neighbor elements and , we understand that is oriented from to . For any function of , we denote by its traces on and and define its average and jump by and , respectively. When , we naturally define With the help of these notions, we present the definitions of the Broken Sobolev spaces and the functionals on the whole edges. The Broken Sobolev spaces for any real number are defined by equipped with the norm Here with are being the -order generalized derivatives of . The functionals with respect to the averages and jumps on the whole edges are defined by where is the length of and the selected nonnegative real numbers are called penalty parameters.
3.2. Variational Formulation for the Mathematical Model
Let be any function in We multiply (1) by suitable test functions on each element , apply the Green formula, and then sum the resulting equations over the domain to obtain the variational form for the flow continuity equationSimilarly, we obtain the variational forms for the rest of the three equations (2)–(4): Consequently, the variational formulation corresponding to models (1)–(4) is defined as to find the map , satisfying (10)-(11).
3.3. The Upwind Discontinuous Galerkin Finite Element Scheme
Assume that is a positive integer and denotes the space of polynomials of degree less than or equal to on . Then, the finite element spaces used later are chosen to be the subspace of defined by For time discretization we divide the time interval [0, T] by the nodes with time step and let be the node-function value of .
Considering the dominance of convection in a real sedimentation process, we hope that the numerical procedure should recognize the direction of the flow, so as to design a stable algorithm. For this purpose, we introduce the upwind values for and on the interior edge by and , respectively. First we define the upwind values of and on the interior edge by and , respectively as follows:
In discretization, we use backward Euler’s scheme to approximate the time derivatives and use the nonsymmetric interior penalty Galerkin method to the diffusion part, in which the average values of on interior edge are replaced by and the nonlinear convective term is handled by the Lesaint-Raviart upwind technique. By doing so, the desired stability can be guaranteed. For easy computation we use, for example, instead of to linearize the nonlinear terms.
Upon these ideas, we propose our fully upwind discontinuous Galerkin finite element procedure (UWDG) as follows: find the map satisfying The initial values also are clarified by their -orthogonal projections: Here , , and denote the interior and exterior traces of on . If , then .
4. Stability Analysis
In this section, we shall analyze the stability of the upwind DG scheme and prove the existence and uniqueness of the discrete solution by borrowing some lemmas from the literature on finite elements and making some reasonable hypotheses.
4.1. Some Lemmas and Hypotheses
To conduct convergence analysis in the next section we need to define the -orthogonal projections for the unknowns These projections are maps satisfying
Following the practical background of a real sedimentation process, we shall make some reasonable hypotheses on the exact solutions of the model.(H1)There exists a unique solution to the models (1)–(4). Then, one can deduce, at least, the existence of the variational formulations (10)-(11).(H2)There is a positive constant such that .(H3), , , and .(H4), , , and .(H5), , , and .(H6) and .
Next we introduce the following lemmas which are crucial to our convergence analysis.
Lemma 1 (see [12–14]). Assume that (H4)–(H6) hold and and are defined in (19). Then, the following error estimates hold for and :
Lemma 2 (trace theorem [12, 15]). Assume that and Then, there exists a constant such that In particular, if , then
Lemma 3 (discrete Gronwall inequality [15, 16]). Assume that the nonnegative sequences and satisfy Then, we have
Lemma 4 (inverse inequalities [12, 16]). Let Then, there is a constant such that
Lemma 5 (boundedness). Assume that is sufficiently small. Then, there is a constant such that, for ,
Proof. Combining the definition of , the hypotheses (H2) and (H4), and Lemma 1, we obtain Let ; then if , we have On the other hand, by the triangle inequality, we can deduce Let ; then if , we have
4.2. Stability Estimate
Preceding the proof of the stability for (14)–(18), we first present a bounded estimate for , and . Since its proof involves an induction argument with the error bounds of these unknowns we here only give its presentation in Lemma 6, whose detailed proof can be found at the last part of the next section.
For the concise exposition we decompose the errors in the following manner: let
Lemma 6 (boundedness of approximate solution). There exists a small positive constant such that when and , for ,
Theorem 7 (stability estimate). There exists a constant independent of and such that for all ,
Proof. Setting in (14), using Green’s formula for the second term on the left-hand side, we have Using , we see that Combining (34) and (35), we rewrite (14) as the following form: Next we estimate the values of the five terms . By Lemma 6 and Young’s inequality, we obtain Using the trace theorem, we have where is the bilinear form defined in Section 3.1.
Similarly, Combining and (37)–(39), (36) becomes Similarly, setting in (15), in (16), and in (17) and using the same techniques, we have Multiplying (40)-(41) by , summing both sides with respect to over , and using the discrete Gronwall inequality (Lemma 3), we complete the proof of Theorem 7.
From Theorem 7, we see that if the initial-boundary conditions (5)-(6) are zero, then are all equal to zero. This implies that the solution of the UWDG scheme (14)–(18) is unique and thus the existence.
Theorem 8 (existence and uniqueness). For , there exists a unique solution to the UWDG scheme (14)–(18).
5. Convergence Analysis
In this section, we shall conduct error estimates for the UWDG solution and prove the convergence rate based upon the prior-estimate techniques and induction argument. Combining these error estimates and applying the induction argument again, we give the detailed proof for the boundedness of the approximate solution stated in Lemma 6.
5.1. Error Equations
First, we derive the error equation of the flow continuous equation (1). Choosing the test function in (10) and (14), subtracting (10) from (14), and using the definition of in (31), we derive the error equation for the flow continuous equation as follows: Taking into account that there is no diffusion term in the flow continuous equation (1), we apply Green’s formula to the second term on the left-hand side Using , we obtain Similarly, we use Green’s formula for the second term on the right-hand side; thus Combining (43)–(45), we rewrite the error equation (42) into the following form: Similarly, we obtain the other three error equations
5.2. Error Estimates
We begin this subsection by deriving the error estimate for the water depth stated in (46). For this we denote the right-hand terms of (46) by and estimate them term by term.
For , we use Cauchy-Schwarz’s inequality and Young’s inequality to give By noting that and Lemma 6 (boundedness of approximate solution), we split as the sum of and estimate them separately; in which, Collecting the estimates for these 3 terms, we obtain the bound of : Similarly, for on the interior edge, , so we split as the sum of and estimate them separately: in which Collecting the estimates for these 2 terms, we obtain the bound of Similarly, we obtain the bounds for , and For , we use the identity and the triangle inequality to give For and , we apply the equality to have Combining the bounds (48)–(59) for to obtain the estimate for the right-hand terms then applying to the first term on the left-hand side of (46), we finally have the estimate for the water depth
Combining the hypotheses (H2)–(H6) and Lemmas 1–6 and using the same method and technique as these used in (48)–(59), we can obtain the error estimates of for the rest of the error equations (47). Here we only give the results:
Considering (1)–(4) is a coupled sedimentation model, we need to combine (60) and (61) to get the error estimates of . For sufficiently small , we have Multiplying both sides of (62) by , summing from to , and using Lemmas 1 and 3, we have Here, we use to replace to denote the generic constant in order to prove Lemma 6.
By using the triangle inequality and Lemma 1, we can get the main conclusion for the norm error estimate of the UWDG scheme. But the process of (63) needs Lemma 6 (boundedness of approximate solution); now we use an induction argument to prove Lemma 6 for .
Proof. We shall use the induction argument to Lemma 6. For , considering the definitions of and Lemma 1, it is easy to get that the lemma holds valid, that is to say, Assuming Lemma 6 holds for , then repeat the derivation process of (63) above; for , we obtain By the converse inequality (Lemma 4), we have Here is the generic constant in (63) and is the constant in Lemma 4.
Let ; then when satisfy and , we obtain On the other hand, combining Lemma 5 and the triangle inequality, we find Therefore, Lemma 6 holds and (63) is true for . Considering the above analysis result and combining the triangle inequality and Lemma 1, we can prove that the approximate solutions for the depth of fluid , the velocity , the sediment concentration , and the bed height are of convergence rates, respectively.
Theorem 9. Assume that the hypotheses (H1)–(H6) hold, , , and are sufficiently small, there is constant independent of and , and , such that the following suboptimal a priori error estimate holds:
6. Algorithm Flow and Computational Cost
In dealing with practical issues, there are three steps: collecting the river terrain data, selecting the hydrologic parameter, and program calculation. These steps are carefully displayed in the following text and Figure 2.
Step 1. Collect the river terrain data and determine the computational domain, the initial-boundary conditions, and the time range.
Step 2. Obtain the hydrologic parameters by properly selecting the sediment carrying capacity, Manning’s roughness coefficient, the viscosity coefficient, the sediment diffusion coefficient, the coefficient of saturation recovery, the settling velocity of sediment, and the sediment dry bulk density.
Step 3. Program calculation step includes(1)determining the time and space step to get triangular subdivision;(2)calculating , and the unknowns on the time level, by inserting , and the knowns on the time level into the fully upwind discontinuous Galerkin finite element procedure (UWDG);(3)judging whether the scheduled time can be reached, if it is reached, output the results and then the program ends.
Now we analyze the computational cost of the UWDG scheme. Because of the lack of continuity constraints between mesh elements for the test functions, we must define the basic functions for each triangle element. We usually use the monomial functions as the basic functions; in 2D, the basic functions are Let denote the total number of elements ; if and , then the total degrees of freedom is , and , respectively. can be understood as the scale of the computation, and in 2D; here denote the maximum diameter of elements.
From Step 3 above, in order to calculate the unknowns on the time level of , we need the knowns on the time level of . Thus for , the computational complexity has
There are three steps to calculating the time level: generate meshes, form the total coefficient matrix, and solve linear equations. Let , and denote the computational complexity of three steps above; if we use Delaunay triangulation algorithm to generate meshes and use iterative methods to solve linear equations, the computational complexity has
From Theorem 9, the choice of the time step must satisfy to ensure the convergence of the UWDG scheme. Let ; then , and we have Combining (70)–(72), the computational cost is , for example, , that is to say, , the computational cost is ; , , the computational cost is .
7. Numerical Experiments
In this section, we perform two numerical experiments, one is for unsteady state case and the other is for stationary-state case, to verity the theoretically proven convergence results.
Example 1 (unsteady-state case). Let the domain be the unit square and the time interval be . The analytic and smooth solution is prescribed to beThe boundary conditions, initial conditions, and the source terms are then completely determined. We also let , , and the coefficient of saturation recovery , settling velocity of sediment , and the sediment dry bulk density are all equal to one.
We choose the piecewise quadratic polynomials () and the piecewise linear polynomials () in the finite element space , respectively. The errors and convergence orders in norm are displayed in Tables 1 and 3. The norm and convergence orders between the approximate solution and the analytic solution of and are displayed in Table 2. The numbers in brackets are the theoretically predicted convergence orders.
By the data in Tables 1–3, we see that the UWDG scheme works well and the convergence order is greater than the results theoretically predicted by (69), even for which is not included in Theorem 9.
For simplicity of presentation, we use Figures 3 and 4 to describe the contours of the sediment concentration and the velocity field at , respectively, and omit those of the remaining unknowns.
(a)
(b)
(a)
(b)
Example 2 (stationary-state case). We further consider nonconstant bathymetry but stationary case. We take the domain to be the unit square and the time interval to be . The analytic and smooth solution is prescribed to beThe boundary conditions, initial conditions, and the source terms are then completely determined by the analytic and smooth solution. We let and and also let the coefficient of saturation recovery , settling velocity of sediment , and the sediment dry bulk density be all equal to one. Further we assume that Manning’s roughness coefficient is equal to zero; this forces the ted shear stress to be equal to zero and (2) turns to be a vortex without friction. The numerical results are displayed in Tables 4–6. The numbers in brackets are the theoretically predicted convergence orders.
Our numerical tests show that the convergence accuracy of the proposed UWDG scheme in this paper is at least as sharp as ; the results are theoretically predicted by (69), even for which is not included in Theorem 9.
8. Conclusions
In this paper, we design an upwind discontinuous Galerkin finite element method for the 2D sedimentation in an estuary as coupled nonlinear differential system, which obeys the local mass conservation and possesses good stability. There exists a unique solution to the numerical procedure and the discrete solution permits convergence rate. Numerical experiments are conducted to verify our theoretical findings.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper
Acknowledgments
This work is supported in part by National Natural Science Foundation of China under Grant nos. 10971254, 11301311, 11471196, 91130010, 11471194, 11571115, and 11171193, by National Science Foundation under Grant nos. EAR-0934747 and DMS-1216923, by the OSD/ARO MURI under Grant no. W911NF-15-1-0562, and by Shandong Province Natural Science Foundation under Grant no. ZR2014AZ001.