A characteristic particle method for the Saint Venant equations

Computers & Fluids 76 (2013) 58–72

Contents lists available at SciVerse ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

A characteristic particle method for the Saint Venant equations Yao-Hsin Hwang ⇑ Department of Marine Engineering, National Kaohsiung Marine University, No. 482, Jhong Jhou 3rd Road, Cijin District, Kaohsiung 805, Taiwan

a r t i c l e

i n f o

Article history: Received 10 October 2012 Received in revised form 25 January 2013 Accepted 30 January 2013 Available online 16 February 2013 Keywords: Characteristic particle method Saint Venant equations Dual-state particle Rankine–Hugoniot relation Open-channel flow Dam-break flow Tidal bore flow

a b s t r a c t A novel characteristic particle method is developed for the Saint Venant equations in the present study. Contrary to the conventional fixed-grid or moving particle methods, the present formulation is built by reallocating the computational particles along the characteristic curves. Both the right- and left-running characteristics are faithfully traced and imitated with their associated computational particles. The annoying numerical nuisance in realizing advection term for the fixed-grid arrangement can be effectively gotten rid of. Besides, special particles with dual states to the enforcement of the Rankine–Hugoniot relation are deliberately imposed to emulate the shock structure. Efficacy of this formulation is verified by solving some benchmark problems with significant transient effects. Computational results are meticulously compared with available analytical solutions. It is concluded that the proposed characteristic particle method will be a useful tool to replicate transient phenomena revealed by the Saint Venant equations. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction Although based on a simplified cross-sectional averaged derivation, the Saint Venant equations [1] are proven to be a useful model in coastal, estuarine and river engineering practices. Dynamics of water transportation phenomena on the ground such as tidal circulation in estuarine channel, flood routing and dam-break can be reasonably simulated with these equations [2–5]. Since the equation system consists of nonlinear partial differential equations for which analytical solution is essentially unavailable, search for an apposite numerical methodology becomes a paramount issue in realizing the physical phenomena divulged by these equations. In practical numerical respects, the major obstacle to simulate the Saint Venant equations lies in the dependence of wave celerity upon local flow condition and the nonlinear advection terms. A practical flow may cover from subcritical to supercritical regions such that discontinuities of flow depth and velocity may occur in the resulting solution. The commonly utilized method of characteristics (MOC) [5–7] for gradually-varied flows fails because the characteristic curves may intersect to form a shock wave (hydraulic bore or surge). In fact, multiple solutions must be taken into account in this circumstance for an accurate simulation. In general, there are two plausible approaches to secure numerical solutions for the Saint Venant equations: The fixed-grid and moving particle methods. The fixed-grid method adopts the Eulerian point of view by probing the solution at predefined locations. The working equations are discretized by taking the advantage of ⇑ Tel.: +886 0922731702; fax: +886 75716013. E-mail address: [email protected] 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2013.01.027

the given geometrical relations connected by the stationary nodes. Several effective numerical skills such as finite-difference (FD) [8–10], finite volume (FV) [11,12] and finite element (FE) [13,14] methods are invented in this regard. Numerical treatment of advection terms has been well accommodated and the shock wave can be successfully captured as long as the adopted grid spacing is sufficiently refined [3,5]. However, the fixed-grid methods still suffer from inevitable numerical annoyances in practical simulations. The stable first-order-accurate upwind scheme may introduce overwhelming numerical diffusivity and subsequently deteriorates the solution accuracy. On the other hand, adoption of a high-orderaccurate scheme may indeed eliminate the superfluously numerical diffusivity; however, unphysical oscillatory solution surfaces unless a supplementary limiter function is imposed to constraint the solution gradients [15–17]. On the contrary, the moving particle method is formulated based on the Lagrangian standpoint where the computational particles are maneuvered to move along the flow path. Realization of flow advection term is taken care by the particle moving strategy and the numerical drawbacks in a fixed-grid method can be effectively eradicated [18–21]. However, since the computational particles are designed to moves along with the flow field, they may become randomly distributed due to substantial flow velocity variation. As a consequence, it is unrealistic to expect a regular particle distribution to yield accurate representations of differential operators [22,23]. Meanwhile, solution characteristic features are basically transmitted along the flow characteristic curves instead of the flow path to formalize the compatibility relations. The simulation accuracy will be deteriorated if the computational particles are dealt with as material points. Therefore, to make a more

59

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

intimate portrayal, the conventional particle method must be modified to put up this crucial trait. Furthermore, discontinuity solution occurs when characteristic curves converge to form a shock. The differential governing equations become inappropriate to describe this peculiar multi-state situation. They must be transformed into their corresponding integral form to accommodate solution discontinuity. That is, the Rankine–Hugoniot relation must be incorporated in the numerical methodology for a reasonable simulation [3,5,17,24,25]. The objective of present study is to develop an accurate and effective numerical tactic for the Saint Venant equations. It is also formulated on the Lagrangian point of view. Nevertheless, in contrast with the conventional particle methods, the computational particles are repositioned to follow the characteristic curves rather than the flow path. This is one of the noticeable distinctions in the present formulation. It is also the reason why the present approach is coined as the characteristic particle method: Particles moving along with characteristic curves. In this manner, individual characteristics are unequivocally traced and replicated by their respective computational particles. Therefore, there are two families of particles taken up to tackle the unsteady features. The right-running characteristic will be handled by the right-running particle and the left-running characteristic by the left-running particle. Based on this distinctiveness, the present formulation is categorized as a dual-particle system, which is also a vital difference from a conventional particle method. Meanwhile, commencement of a shock (surge wave) can be identified as two particles in the same family collide, which is in accordance with the intersection of two characteristic curves. Once a shock is detected, a particular dual-state particle is imposed to mimic its intricacy. Its states are deduced to comply with the Rankine–Hugoniot relation across a shock. In this sense, this particle is designated as a shock particle. On the other hand, those following the characteristic curves to sort out the compatibility equations are regarded as regular ones. To the author’s limited knowledge, introduction of a dual-state particle has not been engineered in the existing literature. Validation of present formulation is conducted by solving some benchmark problems with significant transient effects. Those with negligible nonlinear advection term in the momentum equation suggested by Sobey [26,27] are computed to verify its feasibility for the gradually-varied flow. Numerical results are directly compared with available analytical solutions. Subsequently, cases with significant advection term are considered by solving dam-break and tidal bore problems. Both shock and depression waves are taken into account in this category. They cannot be well resolved by a conventional method of characteristics or linearized analysis. It will demonstrate the ability of present proposition in mimicking the evolvements of wave patterns. Comparisons with available exact solutions are also performed to depict its accuracy. Based on the derivation procedure and the obtained numerical results, it is concluded that the proposed characteristic particle method will be a useful tool to calculate the Saint Venant equations. 2. Characteristics of Saint Venant equations In this section, working equations and their essential features required for the present formulation are briefly outlined. It emphasizes the transmitting characteristics and Rankine–Hugoniot relation to depict an admissible shock. Numerical treatments for these equations will be detailed in Section 3. 2.1. Working equations Under certain assumptions, the governing continuity and momentum equations for a horizontal open-channel flow can be

described by the following dimensionless Saint Venant equations [1–5]:

@ g @q ¼0 þ @t @x  2 @q @ q @g þa þ ¼ j @t @x a @x

b

ð1aÞ ð1bÞ

where t is the temporal variable, x is the spatial variable, b is the width of channel section at free surface, g is the vertical distance of free surface or flow depth, q is the rate of discharge, a is the flow area and j is accounted for the friction term. The friction term can be estimated with the Darcy-Weisbach, Chezy or Manning models, which is generalized by the following form [26–28]:

j ¼ fq q

ð2Þ

Geometrically, both the sectional width and flow area are functions of flow depth:

b ¼ bðgÞ and a ¼ aðgÞ

ð3Þ

Furthermore, the flow area and sectional width are interrelated through the following relation:

bðgÞ ¼

daðgÞ dg

or aðgÞ ¼

Z

g

bðnÞdn

ð4Þ

0

Without loss of generality, the sectional width at referenced flow depth is assigned as:

bð1Þ ¼ 1

ð5Þ

If we introduce the cross-sectional averaged velocity:

v ¼ q=a

ð6Þ

the working equations (1a) and (1b) can be rearranged as:

rffiffiffi rffiffiffi rffiffiffi b @g b @g a @v þv þ ¼0 a @t a @x b @x @v @v @g þv þ ¼ fq v @t @x @x

ð7aÞ ð7bÞ

By adding and subtracting these equations will respectively lead to:

@ðw þ v Þ @ðw þ v Þ þ ðv þ cÞ ¼ fq v @t @x @ðw  v Þ @ðw  v Þ þ ðv  cÞ ¼ fq v @t @x

ð8aÞ ð8bÞ

with the wave celerity c and flow characteristic function w:

sffiffiffiffiffiffiffiffiffiffi aðgÞ cðgÞ ¼ bðgÞ

ð9aÞ

and

wðgÞ ¼

Z 0

g

dn ¼ cðnÞ

Z

g

0

sffiffiffiffiffiffiffiffiffi bðnÞ dn aðnÞ

ð9bÞ

Incorporating with the above relations, the Saint Venant equations (1a) and (1b) can be expressed by the following compatibility relations: (i) Along the characteristic curve of

dx ¼ kþ ; dt

þ

dw fq fq þ wþ ¼ w dt 2 2

ð10aÞ

(ii) Along the characteristic curve of

dx ¼ k ; dt



dw fq fq þ w ¼ wþ dt 2 2

ð10bÞ

60

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

with the characteristic transmitting velocities and variables:

k ¼ v  c

and w ¼ w  v

(i) For a right-running shock,

ð11Þ

In this way, the original partial differential equations (1a) and (1b) are transferred into ordinary ones along the characteristic curves. For discrimination, the first one with dx ¼ kþ is regarded as the dt right-running characteristic and the second one with dx ¼ k the dt left-running one (kþ > k ). It is noted that the original dependent variables can be secured from characteristic variables, v = (w+  w)/2 for instance. These compatibility equations are the bedrock for the method of characteristics [2–7].

As is implied by the above characteristic analysis, the characteristic curves may be divaricated to form depression waves or agglomerated to shocks. In the latter case, they will merge together such that it becomes impractical to determine a unique characteristic at the shock location. Therefore, relations (10a) and (10b) fail to describe the flow evolution and it must be substituted by the Rankine–Hugoniot relation to link the states across a shock [3,5,17,24,25]. For practice, the working equation system (1) should be rearranged as its equivalent conservative form:

@U @F þ ¼ SU @t @x

ð12aÞ

with the conserved variable, associated flux function and source vectors:

F¼

q q2 a

þ/

!

 and S U ¼

0



j

ð12bÞ

In the above equation, / is regarded as the driving function, which is arisen from the hydrostatic pressure term:

/ðgÞ ¼

Z

ð15aÞ

(ii) For a left-running shock,

wR  wL þ

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   1 1 ð/R  /L Þ ¼ wR  wL  aL aR

ð15bÞ

It is worth noting that an admissible right-running shock exists only þ   for wþ L > wR and a left-running one for wR > wL . 3. Numerical treatments

2.2. Rankine–Hugoniot relation

  a U¼ ; q

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   1 1 ð/L  /R Þ ¼ wþL  wþR wL  wR þ  aR aL

g

aðnÞdn

ð13Þ

0

This equation system (12) can be integrated over a region with solution jump to ensure mass and momentum conservation and yield the resulting Rankine–Hugoniot relation. For concreteness, by integrating this equation system over the space–time control volume depicted in Fig. 1 and taking the limit of infinitesimal control volume (dx, dt ? 0), its consequential integral form reads:

rðaL  aR Þ ¼ qL  qR rðqL  qR Þ ¼ q2L =aL þ /L  q2R =aR  /R

ð14aÞ ð14bÞ

where r is the shock speed, r ¼ and subscripts L and R respectively denote the left and right sides of a shock. Accordingly, the Rankine–Hugoniot relation can be derived by manipulating these two equations to eliminate the shock speed: dx dt



  2 1 1 qL qR  ð/L  /R Þ ¼  aR aL aL aR

For computational convenience, this relation is further categorized in the account of right- and left-running shocks,

Fig. 1. Space–time control volume for a moving shock.

In this section, we will describe the numerical tactics adopted in the present study to resolve the governing equations (1a) and (1b). Since a shock may emerge in the course of a simulation, both the regular and shock particles must be taken into consideration for a feasible numerical simulation. The regular particle is utilized to monitor solution advancement along characteristic curves; while the shock particle is used to emulate the states across a shock. It should be emphasized that these particles employed in the present study is not material ones adherent with specified mass as in conventional methods. They can be regarded as pure computational nodes moving along flow characteristics or shock waves. 3.1. Regular particles As was depicted in Eqs. (10a) and (10b), there are two types of characteristic curves with their corresponding compatibility equations to render the flow evolution. It is a natural preference by adopting two families of computational particles in tracking these characteristics. These particles are designed to move along their associated characteristic curves. The corresponding characteristic variables are computed according to the discretized expressions for the compatibility equations (10a) and (10b): (i) For the right-running characteristic, þ n xnþ1 pþ ¼ xpþ þ kpþ Dt

and wþ;nþ1 pþ

;n fq Dt=2 ¼ ðwþ;n þ w;n pþ  wpþ Þe pþ

ð16aÞ

(ii) For the left-running characteristic,  n xnþ1 p ¼ xp þ kp Dt

¼

ðw;n p



and w;nþ1 p

fq Dt=2 wþ;n p Þe

þ wþ;n p

ð16bÞ

In these equations, superscript n and n + 1 respectively denote the computational results at previous and current steps; Dt is the time step and subscripts p+ and p are employed to signify the rightand left-running computational particles, respectively. In this manner, the right-running characteristic of a right-running particle   wþ and the left-running characteristic of a left-running particle pþ    wp are updated by the compatibility equations. It is worth noting that these difference equations will yield an unconditionally stable scheme regardless of adopted time step.To settle the solution state and successfully proceed the computation, both the right- and left-running characteristic variables at a computational particle must be aptly updated. In general, a right-running particle may not happen to be coincident with a left-running one. As men tioned above, wþ pþ and wp are renewed with the compatibility relations (16a) and (16b), respectively. Therefore, the left-running   must be characteristic variable at a right-running particle w pþ secured with an appropriate interpolation scheme from those at   left-running particles w to settle the solution state and vice p

61

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

versa. Although much complicated schemes can be used, a linear approximation is employed in the present investigation for its simplicity. As demonstrated in Fig. 2a, wpþ;nþ1 at the objective rightþ running particle (p+) is determined from compatibility relation and w;nþ1 at two neighboring left-running par(16a); while w;nþ1 L R ticles (L and R) from Eq. (16b). The linear approximation between these two neighboring particles will lead to the interpolated characteristic variable

w;nþ1 pþ

¼

nþ1 xnþ1 R   xp þ

xnþ1 R 

w;nþ1 L xnþ1 L

þ

wp;nþ1 þ

expressed as:

nþ1 xnþ1 pþ  xL

xnþ1 R 

w;nþ1 R xnþ1 L

ð17aÞ

As a consequence, the resultant solution state at the right-running particle will be secured without vagueness. Similar treatments applied to the left-running particle (p) depicted in Fig. 2b will result in its right-running characteristic variable expressed as:

wþ;nþ1 ¼ p

xnþ1  xnþ1 p Rþ xnþ1 Rþ



xnþ1 Lþ

wþ;nþ1 þ Lþ

nþ1 xnþ1 p  xLþ

xnþ1  xnþ1 Rþ Lþ

wþ;nþ1 Rþ

ð17bÞ

To the author’s limited knowledge, induction of two-family particles to simulate the Saint Venant equations has not been practiced in the existing literature.Although numerical inaccuracy may arise from the interpolation procedure to secure interpolated variables þ (w pþ and wp ), it will not notoriously contaminate the subsequent calculations since these variables are secured from the accurate þ ones defined by the compatibility relations (w p and wpþ ). Indeed, inclusive solution state is required for specifying the transmitting velocity and realizing source terms. Nevertheless, the right-running particle carries the information of right-running characteristic variable, which is not affected by the interpolated left-running characteristic variable and the tempted inaccuracy therein. It is also the case for left-running characteristic variable on a left-running particle. That is, both characteristic variables can be accurately simulated but at different locations. The interpolation procedure is employed to reconstruct the solution distribution with the attained accurate characteristic variables. In fact, this interpolation procedure can be completely omitted for a linear problem without additional source terms. This is also one of the salient features distinguished from a conventional MOC method. In a conventional MOC method, both the characteristic variables at a computational node are obtained with an interpolation scheme if the Courant number deviates from unity. Therefore, its solution accuracy will be consequentially contaminated by the inherent numerical error. This is also the reason why the present proposition will yield more accurate results than a conventional MOC method.

tion of the shock structure. In our present formulation, the characteristic curve experienced by a computational particle can be palpably delineated with the particle-moving strategy. It will provide a simple but realistic route to make out the occurrence of a shock. As illustrated in Fig. 3a, a shock comes into existence if two successive computational particles of same family collide:

xnL þ kL Dt ¼ xnR þ kR Dt

ð18aÞ

In fact, it is also related to the time-step constraint to prevent occurrence of particle overtaking:

Dt 6

xnR  xnL

if

k;n  k;n L R

k;n > k;n L R

ð18bÞ

This condition will be taken as a criterion to restrict the computational time step. Once a shock has been detected, it will move in accordance with the shock speed:

xnþ1 ¼ xns þ rn Dt s

ð19Þ

where subscript s signifies a shock particle. The shock speed is derived from the Rankine–Hugoniot relation (14):

r¼

qL  qR aL  aR

r¼

or

q2L =aL þ /L  q2R =aR  /R qL  qR

ð20Þ

Therefore, the resulting trail of a shock particle may not happen to meet with a regular particle. The compatibility Eq. (16) is not applicable to determine its transmitted characteristic variable. Alternatively, special treatment must be adopted to determine the states across a shock. To illustrate a clearer picture on shock evolution, Fig. 3b shows the ploy suggested in the present study to determine the transmitted variable behind a right-running shock (sLn+1). As shown in this figure, the right-running shock moves from sLn to sLn+1 and its left neighboring regular particle from Ln to Ln+1. That is,

xsLnþ1 ¼ xsLn þ rDt

and xLnþ1 ¼ xLn þ kþLn Dt

A fictitious state sL is brought forth by assuming sLn moves across the shock as a regular particle:

xsL ¼ xsLn þ kþsLn Dt

  and wþsL ¼ wþsLn  wsLn efq Dt=2 þ wsLn

 It is noted that kþ sLn > r or xsL > xsLnþ1 for an admissible shock. The   is acquired from those transmitted variable after the shock wþ sLnþ1

at neighboring and fictitious particles (wþ and wþ sL ) with a linear Lnþ1 interpolation scheme:

3.2. Shock particles A shock is formed by the agglomeration of characteristic curves and, as a result, there does not exist a unique characteristic to pin down the transmitted characteristic variable. Meanwhile, multiple states must be taken into consideration for an accurate representa-

(a) Formation of a shock particle

(a) Right-running characteristic particle

(b) Left-running characteristic particle

Fig. 2. Interactions of characteristic particle.

(b) Left state of a shock particle

(c) Right state of a shock particle

Fig. 3. Characteristics of a shock particle.

62

wþsLnþ1 ¼

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

ðxsL  xsLnþ1 ÞwþLnþ1 þ ðxsLnþ1  xLnþ1 ÞwþsL xsL  xLnþ1

ð21aÞ

Similarly shown in Fig. 3c, the transmitted variable before the shock wþ will be obtained as: sRnþ1

wþsRnþ1 ¼

ðxsRnþ1  xsR ÞwþRnþ1 þ ðxRnþ1  xsRnþ1 ÞwþsR xRnþ1  xsR

ð21bÞ Fig. 4. Treatment of a boundary particle.

with the virtual right-side shock particle with subscript sR:

xsR ¼ xsRn þ kþsRn Dt

  and wþsR ¼ wþsRn  wsRn efq Dt=2 þ wsRn

Since the state ahead of a shock is not affected by the incoming shock, the interpolated variable is realized with the same comportment for its transmitted counterpart given in Eq. (21b):

wsRnþ1 ¼

ðxsRnþ1  xsR ÞwRnþ1 þ ðxRnþ1  xsRnþ1 ÞwsR xRnþ1  xsR

ð21cÞ

However, that behind the shock must be secured with the Rankine– Hugoniot relation (15a) to yield a correct shock speed and ensure mass and momentum conservation:

wsLnþ1  wsRnþ1

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi   1 1 ð/sLnþ1  /sRnþ1 Þ ¼ wþsLnþ1  wþsRnþ1 þ  asRnþ1 asLnþ1 ð22Þ

Evolution of a left-running shock can be handled with the same manner. Once the states at the shock particles have been determined, they can be employed as internal moving boundary conditions to carry on the simulation. In this sense, locations of shocks are explicitly made out without any intermediate nodes. Employment of dual-state shock particles is also a novel proposition in the present work.

3.3.1. Time step Unlike a fixed-grid method to deal with flow advection term, the particle method employing a Lagrangian viewpoint is free from numerical instability constraints on the adopted time step, unless particles in the same family collide to form a shock. Therefore, it is the physical constraint to determine the computational time step by preventing particle overtaking:

ð23aÞ

where u is the particle conveying velocity which may be transmitting velocity for a regular particle (k ) or shock speed of a shock particle (r). Subscript p denotes the minimum function is counted over all consecutive particles and L and R denote two consecutive neighboring particles on the left- and right-side, respectively (xR > xL). In practice, an additional constraint to confine particle moving range is also utilized:

d Dt2 ¼ Crminð Þ p up

ð23bÞ

where d is the initial particle spacing and Cr is the Courant number. This restriction is not imposed by numerical stability criterion but proposed to ensure an accurate simulation. Therefore, the adopted computational time step in the present study gives:

Dt ¼ minðDt1 ; Dt2 Þ

wþ;nþ1 B

    nþ1  wþ; wþ;nþ1 xnþ1  xnþ1 B L2 L1 þ xL1  xB L2   ¼ nþ1  xL 1  xL 2

ð24aÞ

On the other hand, the left-running characteristic variable at this boundary particle is derived from the prescribed boundary condition. For example, if the flow depth (gB) or flow velocity (vB) has been assigned as the physical boundary condition, the left-running characteristic becomes:

w;nþ1 ¼ 2wnþ1  wþ;nþ1 B B B

or w;nþ1 ¼ wþ;nþ1  2v nþ1 B B B

ð24bÞ

respectively. This treatment is analogous to that adopted in the existing literature [3,17] for non-reflective boundary conditions. Conditions for the particle at left boundary can be handled with the same manner. 4. Channel geometric functions

3.3. Time step and boundary treatment

  xR  xL with uL > uR Dt1 ¼ min p uL  uR

can be used to derive the numerical boundary condition for the boundary particle. As indicated in Fig. 4 for the right boundary, computational particle L1 will be removed since it moves across the boundary (xLL > xB ). Nevertheless, the right-running characteristic variable at boundary particle (wBþ;nþ1 ) is determined by assuming a linear distribution between nodes, L1 and L2:

ð23cÞ

3.3.2. Boundary treatment A particle is simply removed if it moves across the computational domain. However, its transmitted characteristic variable

As depicted in the previous derivations, several geometrical functions must be settled for a feasible calculation. These include the channel width b, flow area a, wave celerity c, driving function / and characteristic function w. In the present study, we consider a generalized trapezoidal channel:

bðgÞ ¼ b0 þ b1 g

ð25Þ

In this relation, b0 is the width at channel bottom and b1 corresponds to the slope of channel side-wall. When b0 = 0 and b0 = 1, the channel will be constructed with triangular and rectangular cross sections, respectively. Fig. 5a illustrates the corresponding channel geometry. With this generalized trapezoidal channel, all the related geometrical functions can be analytically derived. Therefore, the flow area is acquired from Eq. (4):

aðgÞ ¼

Z

g

bðnÞdn ¼ b0 g þ

0

b1 g2 2

ð26Þ

The wave celerity given in Eq. (9a) reads:

sffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   aðgÞ b0 g þ b1 g2 =2 pffiffiffi b1 g ¼ ¼ gcH cðgÞ ¼ bðgÞ b0 þ b1 g b0

ð27aÞ

with the celerity function:

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sþ2 cH ðsÞ ¼ 2s þ 2

ð27bÞ

Distribution of this function is displayed in Fig. 5b for reference. The driving function shown in Eq. (13) becomes:

/ðgÞ ¼

Z 0

g

aðnÞdn ¼

b0 g2 b1 g3 þ 2 6

ð28Þ

63

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

2.2 2.0 1.8

ψH

1.6 1.4 1.2 1.0

cH

0.8 0.6 10-2

10-1

100

101

102

103

104

b1η/b0 (b) Celerity and characteristic functions

(a) Cross section

Fig. 5. Geometry and characteristic functions of a trapezoidal channel.

The characteristic function defined in Eq. (9b) gives:

5. Numerical evaluations

sffiffiffiffiffiffiffiffiffi   g pffiffiffiffiffiffi bðnÞ b1 g dn ¼ 2gwH wðgÞ ¼ aðnÞ b0 0 Z

ð29aÞ

where wH(s) can be related to elliptic integrals [29,30]:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi nþ1

"rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðs þ 1Þðs þ 2Þ dn ¼ 2 2s s 0 n2 þ 2n     1 1þs 2þs 2 1þs 2þs Cw  RF ; ; 1  2 RD ; ; 1 þ pffiffiffi u 3s s s s s s

1 wH ðsÞ ¼ pffiffiffi

Z

u

ð29bÞ with Cw(1.19814) the integration constant to satisfy wH(0) = 0. In the above equation, RF(x, y, z) and RD(x, y, z) are the Carlson elliptic integrals:

1 2

Z

1

dt pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðt þ xÞðt þ yÞðt þ zÞ 0 Z 3 1 dt qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi RD ðx; y; zÞ ¼ 2 0 ðt þ xÞðt þ yÞðt þ zÞ3

RF ðx; y; zÞ ¼

ð30aÞ ð30bÞ

It is noted that this characteristic function p isffiffiffi a monotone increasing function and its value in the range of 2 6 wH ðsÞ 6 2. In fact, one can easily find that wH(s) complies with the following asymptotic relations:

limwH ðsÞ ¼ s!0

pffiffiffi 2 and

lim w ðsÞ s!1 H

¼2

ð31Þ

Fig. 5b also depicts the distribution of this characteristic function. Table 1 lists the definitions of all associated geometric functions considered in the present study.

Three categories of test problems with various unsteady flow characteristics are numerically solved to validate the present proposition. The first group consists of six problems with inconsequential flow advection term in the momentum equation. The working equations (1a) and (1b) will be reduced to a linear form, which can then be analytically solved with series solutions. In this sense, they can be served as benchmarks for computational flood and tide codes [26,27]. The second and third groups emphasize the effects of nonlinear advection terms which will bring on strong depression or shock waves. Computations are also verified with available exact solutions [31]. To perform an efficient simulation with reasonable accuracy, the computational particles are initially assigned with uniform distribution and managed to ensure sufficient number to resolve flow features. A particle is inserted when two successive ones are separated with distance larger than a prescribed value and is removed if it moves across the computational boundaries or internal shocks. For demonstration, Fig. 6 illustrates the evolution of computational particles in two computational time steps. In this figure, Rkj and Lkj respectively designate the jth right- and left-running regular particles implanted at the n + k time step and RB and LB are their associated boundary particles, which are remained motionless. The solid and dashed lines are employed to indicate the right- and left-running characteristic curves, respectively. As shown in this figure, particles R02 and R03 collide to yield a shock particle s at the n + 1 time step. Heavy solid line indicating the track of the shock is also displayed in this figure. Meanwhile, both R11 and L11 are inserted in the computational domain to increase solution accuracy.

Table 1 Geometric characteristic functions. Function

Definition

Rectangular

Triangular

Trapezoidal

b(g) a(g)

Given Rg 0 bðnÞdn q ffiffiffiffiffiffiffi

1

g g2/2 pffiffiffiffiffiffiffiffiffi g=2

b0 + b1g b0g + b1g2/2 pffiffiffi b1 g gcH b0

g3/6 pffiffiffiffiffiffi 2 2g

b0g2/2 + b1g3/6 pffiffiffiffiffiffi b g 2gwH b10

c(g)a /(g) w(g)b a

b

aðgÞ bðgÞ

Rg

0 aðnÞdn ffi R g qffiffiffiffiffiffi bðnÞ 0 aðnÞdn

g

pffiffiffi

g

2

g /2 pffiffiffi 2 g

qffiffiffiffiffiffiffiffiffiffiffi ðsþ2Þ 2ðsþ1Þ.

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C ðsþ1Þðsþ2Þ  1s RF ð1þs s ; 2þs s ; 1Þ  32s2 RD ð1þs s ; 2þs s ; 1Þ þ pwffiffis. wH ð s Þ ¼ 2 2s c H ðsÞ ¼

Fig. 6. Evolution of computational particles.

64

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

At the n + 2 step, particle R01 collides with shock particle and is removed from subsequent calculation. Also both the particles L05 and L06 are removed since they will otherwise move across the shock. Similarly, particles R21 ; R22 ; R23 and L21 are added in the regions with deficient particle distribution. It is noted that particle L21 is inserted in the rear region of the shock since the left-running particles ahead of the shock (i.e. L05 and L06 ) cannot move across the shock.

equation is omitted. Consequently, they can then be analyzed with series solutions based on the principle of superposition [26,27]. Series solution with 200 terms is employed as the analytical reference to compare with numerical results furnished by present formulation. All calculations are carried out in the dimensionless region of 0 6 x 6 1 with initially uniform particle distribution of d = 1/128. Related characteristic quantities for length (L0), area (A0), discharge rate (Q0), time (s0), friction factor (F0) and frequency (X0) adopted in these test problems are listed in Table 2 for reference. Titles affixed to these test problems are directly borrowed from those in Sobey’s work [27].

5.1. Problems without advection term in momentum equation The governing equations (1a) and (1b) will be simplified to become a linear form if the nonlinear advection term in momentum

5.1.1. AT1: Steady flow and backwater curves The initial condition (IC) and boundary condition (BC) for this case are given as:

Table 2 Characteristic quantities in test problems. L0 (m) 5000

A0 (m2)

Q0 (m3/s)

50

s0 (s)

250

0.20 0.15

a

F0 (s1) 103

1000

X0 (s1)) 103

IC : aðx; 0Þ ¼ a0 þ m0 x and qðx; 0Þ ¼ q0 þ Dq

ð32aÞ

BC : að0; tÞ ¼ 0 and qð1; tÞ ¼ q0

ð32bÞ

-0.20

t=0.2 t=0.4 t=0.6 t=0.8 t=1.0

-0.25

0.10

q

-0.30 t=0.2 t=0.4 t=0.6 t=0.8 t=1.0

0.05

-0.35

0.00 -0.05 0.0

0.2

0.4

0.6

0.8

-0.40 0.0

1.0

x

0.10

1.0

0.15

a

0.05

0.00

0.00

-0.05

-0.05

-0.10

-0.10 -0.15

0.5

1.0

1.5

x=0.2 x=0.4 x=0.6 x=0.8

0.10

0.05

2.0

2.5

3.0

3.5

4.0

0

(c) Flow area at various locations 0.00

0.00 -0.05

x=0.2 x=0.4 x=0.6 x=0.8

-0.15

40

50

60

-0.25

-0.30

-0.30

-0.35

-0.35 2.0

80

90

100

90

100

-0.20

-0.25

1.5

70

x=0.2 x=0.4 x=0.6 x=0.8

-0.15

q

1.0

30

-0.10

-0.20

0.5

20

(e) Long-term response of flow area

-0.05 -0.10

10

t

t

-0.40 0.0

0.8

0.20 x=0.2 x=0.4 x=0.6 x=0.8

0.15

q

0.6

(b) Discharge rate at various times

0.20

-0.15 0.0

0.4

x

(a) Flow area at various times

a

0.2

2.5

3.0

3.5

4.0

-0.40

0

10

20

t

30

40

50

60

70

80

t

(d) Discharge rate at various locations

(f) Long-term response of discharge rate

Fig. 7. AT1: Steady flow and backwater curves.

65

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72



ð34aÞ

  x  t  x0 qðx; tÞ ¼ 0:25a0 1 þ cos p w

  x þ t  x0  0:25a0 1 þ cos p w

ð34bÞ

There are two characteristic waves move in the upward and downward directions, respectively. The computational results as well as analytical solutions (depicted with solid lines) for flow area at various times and locations are displayed in Fig. 8a and b, respectively. The original wave is split into two characteristic waves which can be clearly identified in Fig. 8a. Simulation results with the proposed formulation imitate this phenomenon very accurately. Interestingly, the flow area evolutions at x = 0.2 and x = 0.4 respectively coincide with those at x = 0.8 and x = 0.6, which can be observed both in the computational results and analytical solutions.

5.1.2. AT2: Free propagation This problem is proposed to simulate free propagation of a flow depth wave in a frictionless channel (fq = 0). The initial condition (IC) and boundary condition (BC) for this case are given as:

IC : aðx; 0Þ  0 0:5a0 ½1 þ cosðp xx Þ x0  w 6 x 6 x0 þ w w and qðx; 0Þ ¼ 0 ¼ 0 otherwise ð33aÞ ð33bÞ

with a0 = 0.2, x0 = 0.5 and w = 0.15. This is a simple sinusoidal wave centered at x0 with wave length 2w. The analytical solution before the propagation waves reach boundaries (t 6 0.35) read:

5.1.3. AT3: Boundary reflection This problem explores the wave propagation in a frictional channel as well as its reflection from the boundaries. The associated initial and boundary conditions are set as:

0.10

0.10

0.08

0.08

0.06

0.06

a

a 0.04

0.2

0.4

0.6

x=0.2 x=0.4 x=0.6 x=0.8

0.04

t=0.1 t=0.2 t=0.3

0.02 0.00 0.0



x  t  x0 aðx; tÞ ¼ 0:25a0 1 þ cos p w

  x þ t  x0 þ 0:25a0 1 þ cos p w

with a0 = 0.08, m0 = 0.1, q0 = 0.2, Dq = 0.04 and the friction coefficient is assigned as fq = 0.1. The computational results for flow area and discharge rate at various computational times are respectively shown in Fig. 7a and b. Fig. 7c and d respectively illustrates the evolution of flow area and discharge rate at various locations. As compared with the analytical solutions (depicted with solid lines), it is shown that the proposed characteristic particle method is capable of providing an accurate simulation for this problem. Both the locations and strengths of propagation waves are reproduced in highfidelity. Fig. 7e and f show the long-term response of flow area and discharge rate, respectively. As expected, the flow will approach steady state (q = 0.2) with diminishing oscillation amplitude.

BC : að0; tÞ ¼ 0 and qð1; tÞ ¼ 0



0.02

0.8

0.00 0.00

1.0

0.05

0.10

x

0.15

0.20

0.25

0.30

t

(a) Flow area at various times

(b) Flow area at various locations

Fig. 8. AT2: Free propagation.

0.09 t=0.4 t=0.8 t=1.2

0.07 0.05

0.08

0.04

0.03

a

0.01

a

0.00

-0.01 -0.03

-0.04

-0.05 -0.07 0.0

0.2

0.4

0.6

0.8

1.0

-0.08 0.0

x (a) Flow area at various times Fig. 9. AT3: Boundary reflection.

x=0.2 x=0.4 x=0.6 x=0.8 0.3

0.6

0.9

t (b) Flow area at various locations

1.2

66

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

a

a 0.0

0.0 t=0.25T t=0.5T t=0.75T

-0.1 -0.2 -0.3 0.0

0.2

0.4

0.6

x=0.2 x=0.4 x=0.6 x=0.8

-0.1 -0.2

0.8

-0.3 0.0

1.0

0.2

0.4

0.6

x

0.8

1.0

t/T

(a) Flow area at various times

(b) Flow area at various locations

Fig. 10. AT4: Transition to steady state tidal circulation in partially-closed basin.

0.4

0.4

0.2

0.2

0.0

0.0

a -0.2 -0.4 -0.6

a

x=0.2 x=0.4 x=0.6 x=0.8

-0.2 t=0.2T t=0.4T t=0.6T t=0.8T

0.0

-0.4

0.2

0.4

0.6

0.8

1.0

-0.6 0.0

0.2

0.4

0.6

0.8

1.0

x

t/T

(a) Flow area at various times

(b) Flow area at various locations

6 x=0.2 4 2

a 0 -2 -4 6 -6 x=0.4

t/T

4 2

a 0 -2 -4 -6 0

5

10

15

20

25

30

35

40

41

42

43

44

t/T

45

46

t/T

(c) Long-term response of flow area Fig. 11. AT5: Transition to steady state tidal circulation in open reach.

47

48

49

50

67

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

IC : aðx; 0Þ ¼



0 0:5a0 ½1 þ cosðp xx Þ x0  w 6 x 6 x0 þ w w and qðx; 0Þ ¼ 0 0 otherwise ð35aÞ

BC : að0; tÞ ¼ 0 and qð1; tÞ ¼ 0

ð35bÞ

with a0 = 0.2, x0 = 0.5 and w = 0.15. This is a simple sinusoidal wave centered at x0 with wave length 2w. The friction coefficient is assigned as fq = 1.0. Fig. 9a and b illustrates the computational results of flow area along with analytical solutions (shown with solid lines) at various times and locations, respectively. Although the inclusions of friction factor and boundary reflection may make the wave propagation procedure more complicated, very accurate results are still secured with the proposed formulation. Before the waves reach the boundaries (t 6 0.35), evolutions of flow area at x = 0.2 and x = 0.4 almost coincide with those at x = 0.8 and x = 0.6, respectively. However, it is not the case for the reflected waves. This phenomenon is closely reproduced by the present characteristic particle method shown in Fig. 9b. 5.1.4. AT4: Transition to steady state tidal circulation in partiallyclosed basin This is a typical start up problem of tidal circulation. The associated initial and boundary conditions read:

IC : aðx; 0Þ ¼ 0 and qðx; 0Þ ¼ 0

ð36aÞ

BC : að0; tÞ ¼ a0 sinðxtÞ and qð1; tÞ ¼ 0

ð36bÞ

respectively. Very complicated flow area distributions can be found in these figures. Nevertheless, the proposed characteristic particle method still yields quite accurate results in this problem. Fig. 11c demonstrates the long-term response of flow area at x = 0.2 and x = 0.4. It clearly shows the initial transients are damped by channel friction and the flow response is transitioning to a sinusoidal distribution induced by the boundary condition (37b). 5.1.6. AT6: Flood routing This problem investigates a typical flood routing phenomenon. The associated initial and boundary conditions read:

IC : aðx; 0Þ ¼ m0 ðx  1Þ and qðx; 0Þ ¼ q0 BC : qð0; tÞ ¼



ð38aÞ

q0 þ 0:5qp ½1  cosð2pt=t0 Þ 0 6 t 6 t 0 and að1; tÞ ¼ 0 q0 otherwise ð38bÞ

with m0 = 0.004, q0 = 0.04, qp = 0.2, t0 = 0.2 and fq = 0.1. It is noted that m0 = fqq0 will yield an exact uniform flow balance in the ini-

with the driving amplitude at the left end (x = 0) being set as a0 = 0.2 and its frequency as x = p. It is starting from an initially quiescent condition. The friction coefficient is assigned as fq = 0.1. The computational results as well as the analytical solutions (shown with solid lines) for flow area at various times and locations are displayed in Fig. 10a and b, respectively. With the imposed frequency, the corresponding dimensionless period becomes: T = 2. Again, quite accurate prediction on this problem is produced the proposed characteristic particle method.

(a) Initial flow depth

5.1.5. AT5: Transition to steady state tidal circulation in open reach This is a start up problem from an initially quiescent condition. It is a typical tidal problem driven by flow area variation at both ends but with a phase difference. The associated initial and boundary conditions read:

IC : aðx; 0Þ ¼ 0 and qðx; 0Þ ¼ 0 BC : að0; tÞ ¼ a0 sinðxtÞ and að1; tÞ ¼ a0 sinðxt þ /L Þ

ð37aÞ ð37bÞ

with a0 = 0.16, x = p, /L = 1 and fq = 0.1. Fig. 11a and b illustrates the computational results of flow area along with analytical solutions (shown with solid lines) at various times and locations,

0.20

0.15

(b) Wave pattern and corresponding flow depth Fig. 13. Dam-break problem.

0.20

x=0.2

0.15

t=0.25 0.5 0.75

a 0.10

0.4

a 0.10

1.0

0.6 0.8

0.05

0.00 0.0

0.05

0.2

0.4

0.6

0.8

1.0

0.00 0.0

0.2

x

0.4

0.6

0.8

t

(a) Flow area at various times

(b) Flow area at various locations

Fig. 12. AT6: Flood routing.

1.0

68

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

tial condition. The computational results as well as the analytical solutions (shown with solid lines) for flow area at various times and locations are displayed in Fig. 12a and b, respectively. A flood propagating toward downstream is obviously observed in these figures and the computational results delivered by the proposed characteristic particle method yield a very accurate imitation.

shock. A feasible numerical tool must be verified to yield faithful results in this circumstance in which both the gradually- and rapidly-varied flow features are accommodated [3,5]. In the present study, we consider an initially quiescent channel separated by a gate into two regions with different flow depths, which is schematically illustrated in Fig. 13a. Without loss of generality, the flow depth in the left region is assumed to be larger than the right region (gL > gR or gR/gL 6 1). As the gate is suddenly removed, there will form a shock and depression wave running into the right and left regions, respectively. During this wave propagation procedure, an intermediate state emerges between these two waves. Fig. 13b schematically depicts this scenario and the corresponding distribution of flow depth. For comparison, channels with generalized trapezoidal cross sections shown in Fig. 5 are considered. All simulations were conducted with initially uniform particle

5.2. Dam-break problems From the computational results for linearized problems shown above, it is concluded that the proposed characteristic particle method is a very accurate tool if the advection term in momentum equation is insignificant. Cases with substantial advection term will be investigated with the dam-break problem. In the dambreak problem, the characteristic lines will converge to form a

1.2

1.2 exact

1.0

exact

1.0

η

η

v

v

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.0

0.0

-0.2 0.0

0.2

0.4

0.6

0.8

-0.2 0.0

1.0

0.2

0.4

x

0.6

0.8

1.0

x

(a) η R / η L = 0.5

(b) η R / η L = 0.1

2.5 exact

η

2.0

v

1.5 1.0 0.5 0.0 0.0

0.2

0.4

0.6

0.8

1.0

x

(c) η R / η L = 0 Fig. 14. Computational results of a rectangular channel at t = 0.2.

2.0

2.0

exact

1.5

SPH (n=1000), Chang et al. SPH (n=5000), Chang et al. present (n=64)

exact FV (n=120), Ying et al. present (n=64)

1.5

v

v 1.0

1.0

0.5

0.5

0.0 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.0 0.2

0.4

0.6

x

x

(a) SPH

(b) FV Fig. 15. Comparisons with SPH and FV methods.

0.8

1.0

69

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

distribution with d = 1/64 in the range of 0 6 x 6 1 without taking the friction term into account (fq = 0). Two distinct quiescent states with different flow depths are separated by a gate at xm = 0.5. First, a channel with rectangular cross section (b0 = 1) is considered. Fig. 14a–c shows the computational results of flow depth and averaged velocity at t = 0.2 for gR/gL = 0.5, 0.1 and 0, respectively. Exact solutions [31] are also included in these figures for comparison. It should be noted that the case of dry bed (gR/gL = 0) is actually calculated with gR/gL = 1010 to prevent numerical nuisance induced by computer round-off error. As clearly demonstrated in

2.0

exact

η*

1.5

v* xs (t=0.2) 1.0

0.5

0.0 0.0

0.2

0.4

0.6

0.8

1.0

ηR /ηL Fig. 16. Intermediate states for a rectangular channel.

these figures, the proposed characteristic particle method can yield very accurate simulations irrespective of initial flow depth ratio. It is worth emphasizing that the simulation accuracy does not deteriorated for the case of gR/gL = 0, which is difficult to be accurately predicted with a fixed-grid or conventional particle methods [10,21]. Fig. 15a and b shows the comparisons with SPH [21] and FV [10] methods in this case, respectively. In the SPH method [21], the slice water particles are employed to solve the onedimensional shallow water equations and a variable smoothing length strategy is used to economize the computations. In the FV method [10], a first-order accurate upwind conservative scheme is developed for the flow advection and the flow driving term is treated as a source term. That is, it is the flow velocity rather than the flow characteristics to work out the conservative flux function. It is clearly illustrated in these figures that the present characteristic particle method can provide more accurate results than SPH and FV methods even with much less computational particles. The SPH method can yield very steep shock front but with erroneous shock location; while the FV method predict the correct shock location but with rather smeared front. These shortcomings are effectively eliminated by the present characteristic particle method. To highlight the numerical performance of present formulation, the resulting intermediate state of flow depth and averaged velocity (g and v) and the shock location at t = 0.2 in terms of initial flow depth ratio (gR/gL) are exhibited in Fig. 16. Again, these quantities are precisely reproduced by the proposed method. Since the case with dry bed (gR/gL = 0) incurs severe numerical difficulties, we further calculate the dam-break problems with triangular (b0 = 0) and trapezoidal (b0 = 0.5) cross sections. The

1.2

3.0

exact t=0. t=0.05 t=0.10 t=0.15

1.0 0.8

η

0.6

exact t=0. t=0.05 t=0.10 t=0.15

2.5 2.0

v

0.4

1.5 1.0

0.2 0.5

0.0 0.0

-0.2 0.2

0.4

0.6

0.8

0.2

1.0

0.4

0.6

0.8

1.0

x

x

(a) Flow depth at various times

(b) Averaged velocity at various times

Fig. 17. Computational results of a triangular channel with gR/gL = 0.

2.5

1.2 exact t=0. t=0.05 t=0.10 t=0.15

1.0 0.8

η

0.6

exact t=0. t=0.05 t=0.10 t=0.15

2.0 1.5

v

0.4

1.0

0.2

0.5

0.0 0.0

-0.2 0.2

0.4

0.6

0.8

x

(a) Flow depth at various times

1.0

0.2

0.4

0.6

0.8

1.0

x

(b) Averaged velocity at various times

Fig. 18. Computational results of a trapezoidal channel with gR/gL = 0 .

70

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

3.0 1.0 2.5

0.8

η 0.6

2.0

1.25 1.75

t=0

0.0 0.8 0.2

0.4

0.6

0.8

0.6

1.0

t=0 t=0.25 t=0.5 t=0.75 t=1.0 t=1.25 t=1.5 t=1.75

v

0.2 0.0

0.5

t=0

-0.2

1.5 1.75

predicted distributions of flow depth and averaged velocity as well as the exact solutions for the triangular channel at various times are demonstrated in Fig. 17a and b; while those for the trapezoidal channel are illustrated in Fig. 18a and b. As shown in these figures, the present formulation delivers simulation results in a very good agreement with exact solutions. Apart from the above simulations, we also investigate the intermediate states of gR/gL = 0 for channel with generalized trapezoidal cross section. The computational results are shown in Fig. 19 along with the shock locations at t = 0.15. The exact values of these parameters can be derived from the flow characteristics and Rankine–Hugoniot relation [31]:

  1  b0 and xs ¼ xm þ v  t b0

0.75 1.0

1.25

0.25

Fig. 19. Intermediate states of generalized trapezoidal channels with gR/gL = 0.

pffiffiffi

0.5

x

0.4

b0

g ¼ 0; v  ¼ 2wH

0.25

0.2

1.0

0.5 0.0

0.75

0.4

exact v* xs

1.5

1.0

1.5

-0.4 -0.6 0.0

0.2

0.4

0.6

0.8

1.0

x

(a) Computational results at various times

ð39Þ

Again, the proposed characteristic particle method produces computational results in a very agreement with exact solutions, which is evidently depicted in this figure.Finally, we demonstrate a simulation of dam-break problem with gR/gL = 0.2 but with closed ends at simulation boundaries (x = 0 and x = 1). Accordingly, the corresponding boundary conditions read:

qð0; tÞ ¼ qð1; tÞ ¼ 0

ð40Þ

It will be a close resemblance for the propagation and reflection procedure of a hydraulic bore. Fig. 20a shows the distribution of flow depth and averaged velocity at various times. In this figure, solid and dashed lines are employed to signify the shock waves moving in the right- and left-ward directions, respectively. Besides the propagations of depression and shock waves, reflections from the solid walls are also observed in this simulation. It clearly demonstrates the numerical treatment of the solid walls based on the flow characteristics given in Eq. (24) is a non-reflective boundary condition. No noticeable spurious wiggle in the reflected waves delivered by the computational results can be observed. These observations are further verified by the wave pattern of flow depth in Fig. 20b, where the locus of a propagating shock (hydraulic bore) is clearly identified.

(b) Wave pattern of flow depth Fig. 20. Dam-break problem with closed ends.

5.3. Tidal bore problem The initial condition is imposed with constant flow discharge rate: A tidal bore problem in a horizontal frictionless channel with divergent rectangular channel is simulated to verify the present proposition. The associated channel width and flow area are related to the flow depth,

bðg; xÞ ¼ a þ b tanh½cðx  0:7Þ and aðg; xÞ ¼ gfa þ b tanh½cðx  0:7Þg

ð41Þ

qðx; 0Þ ¼ qin

ð42aÞ

Besides, the initial flow depth is derived from the following relation to yield a steady state distribution:

q2in 1 q2 ðx; 0Þ 1 g ðx; 0Þ ¼ þ þ gL 2 a2 ðgðx; 0Þ; xÞ 2 a2 ðgL ; 0Þ

ð42bÞ

71

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

where v B0 is the initial velocity derived from Eq. (42b) at x = 1. Calculation is carried out with the following parameters: fq = 0, qin = 0.5, gL = 1.0, tB = 0.2 and vB = 0.1. The distribution of channel width depicted in Fig. 21a is determined by the following conditions:

The bore is triggered by the downstream boundary condition at x = 1,



vB vB

0

þ ðv B  v B0 Þt=t B

t 6 tB

ð43Þ

t > tB

5 4

α+β tanh[γ (x-0.7)]

3 2 1

b 0 -1 -2 -3 -4 -5 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

x

(a) Distribution of channel width 1.8

1.8

1.6

1.6

η 1.4

1.4

η

1.2

1.2 1.0

1.0 t=0.2 t=0.4 t=0.6 t=0.8 t=1.0

0.8 0.6 0.4

0.6 0.4

v 0.2 0.0

0.0 -0.2 0.0

x=0.3 x=0.5 x=0.7 x=0.9

0.8

v 0.2

0.2

0.4

0.6

0.8

-0.2 0.0

1.0

0.1

0.2

0.3

0.4

x

0.9

1.45

1.45

1.5

0.8 1.5

0.7 0.6

1.4 1.35

0.5 0.4

1.3 1.025

0.3

1.25 1.05 1.2

0.2 1.075

0.1

1.15

1.1

0

0.6

0.7

0.8

0.9

1.0

(c) Flow depth and velocity at various locations

1

0

0.5

t

(b) Flow depth and velocity at various times

t

v ð1; tÞ ¼

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x

(d) Wave pattern of flow depth Fig. 21. Computational results of a tidal bore problem.

1

72

bðg; 0Þ ¼ 1:0;

Y.-H. Hwang / Computers & Fluids 76 (2013) 58–72

bðg; 1Þ ¼ 10:0 and c ¼ 10:0

ð44Þ

It indicates rapid cannel width enlargement in the range of 0.6 6 x 6 0.8. Owing to the effects of non-prismatic channel, the characteristic equations (8a) and (8b) should be modified as:

  @ðw þ v Þ @ðw þ v Þ v @a þ ðv þ cÞ ¼ fq v  @t @x bc @x g   @ðw  v Þ @ðw  v Þ v @a þ ðv  cÞ ¼ fq v  @t @x bc @x g

ð45aÞ ð45bÞ

Nevertheless, essential computational tactics exhibited in the present study is applicable to solve this equation system. Computational results for flow depth and velocity at various times and locations are depicted in Fig. 21b and c, respectively. It clearly shows the establishment and propagation of a tidal bore. The flow depth (tidal elevation) will gradually increase as it travels in the upstream direction. Fig. 21d demonstrates the wave pattern of flow depth in the x–t plane. It also manifests the locus of a propagating tidal bore with increasing elevation. 6. Conclusions An accurate numerical method for the Saint Venant equations is proposed and validated in the present work. It is built based on the compatibility equations derived from the governing equation system. Contrary to the conventional fixed-grid or moving particle methods, two families of particles are employed to take care of two transmitted characteristics. The computational particles are managed to follow their individual characteristic curves. Therefore, the present proposition is coined as a characteristic particle method. Besides the regular particles to tackle flow characteristics, auxiliary shock particles with dual states are imposed to imitate the admissible shock structure. They are introduced to the fulfillment of Rankine–Hugoniot relation. The shock location and its propagation speed can then be specified without numerical ambiguity. Numerical difficulties associated with the fixed-grid and moving particle methods can be effectively eliminated. Validations of the present formulation are performed by solving some benchmark problems with significant transient effects. The computational results are compared with available analytical solutions to depict its accuracy. From the derivation procedure and the obtained numerical results, it is conclude that the present characteristic particle method will be a useful tool to simulate one-dimensional open-channel flows modeled by the Saint Venant equations. References [1] De Saint-Venant B. Theorie du mouvement non permanent des eaux, avec application aux crues de rivieras at a l’introduction des marces dans leur lit. Compt Rend Acad Sci 1871;73:147–54. 237–40.

[2] Henderson FM. Open channel flow. New York: McGraw-Hill; 1966. [3] Cunge JA, Holly FM, Verwey A. Practical aspects of computational river hydraulics. London: Pitman Publishing Co.; 1980. [4] Chaudhry MH. Open-channel flow. 2nd ed. New York: Springer; 2008. [5] Dronkers JJ. Tidal computations in rivers and coastal waters. Amsterdam (The Netherlands): North-Holland Publ.; 1964. [6] Abbott MB. An introduction to the method of characteristics. London: Thame and Hudson; 1966. [7] Wylie EB, Streeter VL. Fluid transients in systems. Englewood Cliffs: PrenticeHall; 1993. [8] Preissmann A, Cunge J. Calcul du mascaret sur machines electroniques. La Houille Blanche 1961;5:588–96. [9] Fennema RJ, Chaudhry MH. Simulation of one-dimensional dam-break flows. J Hydraul Res 1987;25:41–51. [10] Wang JS, Ni HG, He YS. Finite-difference TVD scheme for computation of dambreak problems. J Hydraul Eng 2000;126:253–62. [11] Mingham CG, Causon DM. High-resolution finite-volume method for shallow water flows. J Hydraul Eng 1998;124:605–14. [12] Ying X, Khan AA, Wang SY. Upwind conservative scheme for the Saint Venant equations. J Hydraul Eng 2004;130:977–87. [13] Hicks FE, Steffler PM. Characteristic dissipative Galerkin scheme for open channel flow. J Hydraul Eng 1992;118:337–52. [14] Berger RC, Stockstill RL. Finite-element model for high-velocity channels. J Hydraul Eng 1995;121:710–6. [15] VanLeer B. Towards the ultimate conservative difference scheme. V. A second order sequel to Godunov’s method. J Comput Phys 1979;32:101–36. [16] Sweby PK. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM J Numer Anal 1984;21:995–1011. [17] Hirsch C. Numerical computation of internal and external flows. Computational methods for inviscid and viscous flows, vol. 2. New York: Wiley; 1990. [18] Monaghan JJ. Simulating free surface flows with SPH. J Comput Phys 1994;110:399–406. [19] Koshizuka S, Nobe A, Oka Y. Numerical analysis of breaking waves using the moving particle semi-implicit method. Int J Numer Methods Fluids 1998;26:751–69. [20] Ata K, Soulaimani A. A stabilized SPH method for inviscid shallow water flows. Int J Numer Methods Fluids 2005;47:139–59. [21] Chang T-J, Kao H-M, Chang K-H, Hsu M-H. Numerical simulation of shallowwater dam break flows in open channels using smoothed particle hydrodynamics. J Hydrol 2011;408:78–90. [22] Quinlan NJ, Basa M, Lastiwka M. Truncation error in mesh-free particle methods. Int J Numer Methods Eng 2006;66:2064–85. [23] Hwang Y-H. Smoothing difference scheme in a moving particle method. Numer Heat Transfer B 2011;60:203–34. [24] LeVeque RJ. Finite volume methods for hyperbolic problems. Cambridge: Cambridge Press; 2002. [25] Shapiro AH. The dynamics and thermodynamics of compressible fluid flow. New York: John Wiley & Sons; 1953. [26] Sobey RJ. Analytical solution of non-homogeneous wave equation. Coast Eng J 2002;44:1–24. [27] Sobey RJ. Analytical solutions for flood and tide codes. Coast Eng J 2002;44:25–51. [28] Sobey RJ. Evaluation of numerical models of flood and tide propagation in channels. J Hydraul Eng 2001;127:805–24. [29] Carlson BC. Computing elliptic integrals by duplication. Numer Math 1979;33:1–16. [30] Carlson BC, Notis EM. Algorithms for incomplete elliptic integrals. ACM Trans Math Softw 1981;7:398–403. [31] Ritter A. Die Fortpflanzung der Wasserwellen. Zeitschrift des Vereines Deutscher Ingenieure 1892;36:947–54.