Modelling suspended sediment in environmental turbulent fluids
MModelling suspended sediment in environmentalturbulent fluids
Meng Cao ∗ A. J. Roberts † August 21, 2018
Abstract
Modelling sediment transport in environmental turbulent fluids is achallenge. This article develops a sound model of the lateral transport ofsuspended sediment in environmental fluid flows such as floods and tsunamis.The model is systematically derived from a 3D turbulence model based onthe Smagorinski large eddy closure. Embedding the physical dynamics intoa family of problems and analysing linear dynamics of the system, centremanifold theory indicates the existence of slow manifold parametrised bymacroscale variables. Computer algebra then constructs the slow manifoldin terms of fluid depth, depth-averaged lateral velocities, and suspendedsediment concentration. The model includes the effects of sediment erosion,advection, dispersion, and also the interactions between the sediment andturbulent fluid flow. Vertical distributions of the velocity and concentrationin steady flow agree with the established experimental data. Numericalsimulations of the suspended sediment under large waves show that thedeveloped model predicts physically reasonable phenomena.
Contents ∗ School of Mathematical Sciences, University of Adelaide, South Australia 5005. mailto:[email protected] or mailto:[email protected] † School of Mathematical Sciences, University of Adelaide, South Australia 5005. mailto:[email protected] a r X i v : . [ m a t h . D S ] J u l Cross-sectional structures of the flow and sediment 14
Environmental turbulent fluids, such as rivers, floods and tsunamis, always carryamounts of sediment. For example, Figure 1 shows the Yellow River in Chinawhich is famous for carrying large amounts of sediment. Modelling the sediment inthese environmental fluids is important for studying and predicting changes of themorphology and topography. We aim to develop a model to appropriately modelthe suspended sediment in turbulent fluid flows via systematic resolution of thephysical processes.Our modelling, which is based on dynamical systems theory instead of con-ventional depth-averaging, resolves out-of-equilibrium interactions between thevarying turbulence and suspended sediment. Most previous work studied suspendedsediment in uniform flows (Hunt 1954, van Rijn 1984, Celik & Rodi 1988, Fredsoe &Deigaard 1992, e.g.) or by depth-averaging flow and sediment equations (Wu et al.2000, Pittaluga & Seminara 2003, e.g.). We explore the implications of changingthe theoretical base from depth-averaging to a slow manifold of the turbulentSmagorinski large eddy closure.Cao & Roberts (2012) initially developed a 2D lateral model for environmentalfluid flows, derived from a 3D turbulence model based on Smagorinski large eddyclosure. This model includes the effects and interactions of inertia, advection, beddrag, gravitational forcing and turbulent dissipation with minimal assumptions.The innovation here is that the turbulent modelling and dynamics includes andinteracts with the suspended sediment transport (Cao 2014). A slow manifold isfound for the out-of-equilibrium dynamics of the coupled turbulent sediment system.We choose to parametrise the slow manifold in terms of emergent depth-averagedquantities. The evolution of these depth-averaged quantities on the slow manifoldgoverns the dynamics of the suspended sediment in the turbulent fluid flows.Consider a turbulent flow of depth h ( x , y , t ) flowing along a bed z = b ( x , y ) with a mean slope tan θ and carrying sediment. Sections 2.1–2.2 detail equationsof the Reynolds-averaged Navier–Stokes pde s, the advection-diffusion equationand boundary conditions on the free surface and the mean bed. Section 2.3 usesthe dynamical systems theory of centre manifolds (Roberts 1988, Potzsche &2igure 1: Scene of Yellow River turning in Shilou, China’s Shanxi ( http://news.xinhuanet.com ). This river carries vast amounts of suspended sediment.Rasmussen 2006, e.g.) to analyse the governing equations and derive the followingnon-dimensional suspended sediment model of the horizontal evolution of thedepth h ( x , y , t ) , the lateral depth-averaged velocities ¯ u ( x , y , t ) and ¯ v ( x , y , t ) anddepth-averaged concentration ¯ c ( x , y , t ) : ∂h∂t ≈ − ∂∂x ( h ¯ u ) − ∂∂y ( h ¯ v ) , (1a) ∂ ¯ u∂t ≈ − u ¯ qh + (cid:20) tan θ − ∂∂x ( h + b ) (cid:21) − u ∂ ¯ u∂x − v ∂ ¯ u∂y − ( s − ) h ∂ ¯ c∂x , (1b) ∂ ¯ v∂t ≈ − v ¯ qh − ∂∂y ( h + b ) − v ∂ ¯ v∂y − u ∂ ¯ v∂x − ( s − ) h ∂ ¯ c∂y , (1c) ∂ ¯ c∂t ≈ − w f h ( c − c ae ) − (cid:18) − w f ¯ q (cid:19) (cid:18) ¯ u ∂ ¯ c∂x + ¯ v ∂ ¯ c∂y (cid:19) , (1d)where ¯ q = √ ¯ u + ¯ v the mean local flow speed, constant w f is the falling velocityof the sediment, and constant c ae is an equilibrium reference concentration onthe mean bed z = b . The effective momentum equations (1b)–(1c), and the morerefined version (28b)–(28c), include the effects of gravitational forcing, bed drag, self-3dvection, turbulent dissipation, and sediment induced flow. The sediment pde (1d)includes the sediment erosion and deposition, and the advection. Although thismodel is expressed in terms of depth-averaged lateral velocities and concentration,they are derived not by depth-averaging, but instead by systematically accountingfor interaction between vertical profiles and horizontal gradients of the velocity, theconcentration, the stress, bed drag, lateral space variations and bed topography.That the systematically reduced model (1) has many terms so close to establishedmodels is a testament to the robustness of conservation of mass and momentumprinciples that underly traditional derivations. Additionally, our dynamical systemsapproach resolves finer microscale and out-of-equilibrium interactions that areactive in more physically complicated systems.Linear analysis in section 2.4 finds the spectrum supports the existence of aslow manifold. The computer algebra of Appendix A then constructs the slowmanifold of the system. Then section 2.5 derives the reduced model (1) of theturbulent flow and suspended sediment on the slow manifold.Section 3 discusses the predicted vertical distribution of the velocity and concen-tration fields in steady flow. Agreement with established experimental data (Schultz& Flack 2007, 2013, Celik & Rodi 1988, e.g.) indicates that our model is reason-able. Section 4 numerically explores the suspended sediment under large waves bythe comprehensive suspended sediment model (28d) coupled with the turbulencemodel (28a)–(28c). The numerical results indicate that our model reasonablydescribes the dynamics of the suspended sediment in turbulent flows. This section describes the derivation of the reduced models for turbulent flow andsuspended sediment. First, section 2.1 details the 3D continuity, Navier–Stokes andadvection-diffusion equations of the turbulent fluid flows and suspended sediment,whereas section 2.2 records the boundary conditions of the flow and the sediment onthe free surface and the mean bed. Second, section 2.3 embeds these equations in afamily of equations with modified tangential stresses on the free surface to establishthe existence of an appropriate slow manifold. The linear analysis of section 2.4supports the emergence of the slow manifold from the dynamics in the system. Thecomputer algebra of Appendix A handles the details of the construction of the slowmanifold model that is summarised in section 2.5.
Consider the turbulent flow flowing along a bed carrying sediment. This work onlyconsiders the suspended sediment and neglects the bed load transport on the meanbed. Figure 2 depicts the diagram of the suspended sediment in the turbulent4 yz fluid fluid h ( x , y , t ) q ( x , y , z , t ) g c ( x , y , z , t ) z = b ( x , y ) z = h + b Figure 2: This diagram depicts the suspended sediment in turbulent flowwith ( x , y , z ) coordinate system. The fluid of depth h ( x , y , t ) flows down thesloped bed at the turbulent mean velocity q ( x , y , z , t ) . The turbulent mean concen-tration of the suspended sediment is c ( x , y , z , t ) . Denote the mean bed z = b ( x , y ) ,the free surface z = h + b , and the gravity is g .flow. Define a coordinate system with x , y for the lateral directions, and z for thedirection normal to the mean slope.The fluid of depth h ( x , y , t ) flows down the sloped mean bed z = b ( x , y ) atthe turbulent mean velocity q ( x , y , z , t ) ; the velocity vector q = ( u , v , w ) in the ( x , y , z ) directions, respectively. The term ‘mean bed’ refers to the smooth averagebed over an ensemble of turbulence and bed roughness realisations: the physicalbed in any one realisation is envisaged to be rough, just like the physical fluidvelocity field would have rapid spatial variations in any one realisation. Denotethe turbulent mean pressure field by p ( x , y , z , t ) . The suspended sediment hasa turbulent mean concentration c ( x , y , z , t ) (volume fraction). The mean bed z = b ( x , y ) has an overall slope tan θ in the x -direction.We assume the particles of the suspended sediment have small sizes, and allthe particles of the suspended sediment have the same falling velocity, namely thatof a sphere of diameter d .The nondimensional governing partial differential equations for the incompress-ible, three dimensional, turbulent mean fluid fields are the Reynolds-averagedcontinuity and momentum equations, ∇ · q = ∂u∂x + ∂v∂y + ∂w∂z = ∂ q ∂t + q · ∇ q = − ∇ p + ∇ · τ + g , (3)5nd for the suspended sediment is the advection-diffusion equation, ∂c∂t + ∇ · ( q c ) = − ∇ · ( w f c n g ) + ∇ · ( ν ∇ c ) , (4)where w f is the falling velocity of the sediment particles, ν is the eddy viscosity,and the vector g = ( tan θ , 0, − ) is the direction of the nondimensional gravity.The suspended sediment influences the fluid turbulence. We assume the mixingdensity ρ mix of the fluid and suspended sediment satisfies1 ρ mix = ρ + c ( ρ m − ρ ) = ρ + c ( s − )= ρ (cid:2) − c ( s − ) + c ( s − ) (cid:3) + O ([ c ( s − )] ) , (5)where ρ m is the sediment density, ρ the density of the fluid and s = ρ m /ρ is therelative density.The pde s (2)–(5) are nondimensional and derived upon using the characteristicdepth H of the turbulent fluid as the length scale, the long wave speed √ g z H asthe velocity scale, and the fluid density ρ as the reference density. Thus, the mixingdensity ρ mix ≈ − c ( s − ) .The variable τ is the turbulent mean deviatoric stress tensor, which is ap-proximated using the eddy viscosity ν . We use the Smagorinski eddy closure toapproximate the turbulence stresses. Cao & Roberts (2012) and Georgiev et al.(2009) detailed this eddy closure and expressed the mean deviatoric stress tensor τ ij = ν (cid:18) ∂u i ∂x j + ∂u j ∂x i (cid:19) = c t h ˙ ε (cid:18) ∂u i ∂x j + ∂u j ∂x i (cid:19) , (6)where the magnitude of the second invariant of the strain-rate tensor satisfies | ˙ ε | = (cid:80) i , j ( ∂u i /∂x j + ∂u j /∂x i ) .The falling velocity w f is related to the mean particle effective size d , the relativedensity s and the gravity g . We set the falling velocity (Fredsoe & Deigaard 1992,e.g.) w f = (cid:115) ( s − ) gd c D . (7)The drag coefficient c D ≈ >
500 (Fredsoe & Deigaard 1992, e.g.). Typically, a quartz sedimenthas a relative density s = § We formulate boundary conditions on the mean bed z = b ( x , y ) and on thefree surface z = η ( x , y , t ) = h ( x , y , t ) + b ( x , y ) in terms of the turbulent mean6elocity field q ( x , y , z , t ) , turbulent mean concentration c ( x , y , z , t ) , the fluiddepth h ( x , y , t ) , and the turbulent mean pressure p ( x , y , z , t ) .First formulate the boundary conditions for the turbulent flows. • On the mean bed, no fluid penetrating the ground requires w = ub x + vb y on z = b . (8) • Positing a slip law on the mean bed to account for a negligibly thin turbulentboundary layer gives1 (cid:112) + b x ( u + wb x ) = c u h (cid:112) + b x + b y ( u + wb x ) n on z = b , (9)1 (cid:112) + b y ( v + wb y ) = c u h (cid:112) + b x + b y ( v + wb y ) n on z = b , (10)where the derivative ∂ n = − b x ∂ x − b y ∂ y + ∂ z and the constant c u ≈ • On the free surface (that is, on its turbulent mean position), the kinematiccondition that no fluid crosses the free surface is η t + uη x + vη y = w on z = η = h + b . (11) • Zero turbulent mean stress normal to the free surface indicates that on z = η − p + τ − η x τ − η y τ + η x τ + η x η y τ + η y τ + η x + η y = • No turbulent mean, tangential stresses at the free surface indicates that on z = η ( − η x ) τ + η x ( τ − τ ) − η y ( τ + η x τ ) = ( − η y ) τ + η y ( τ − τ ) − η x ( τ + η y τ ) = • On the free surface, the sediment flux normal to the surface is zero, whichrequires ( η x tan θ + ) w f c + ν (− η x c x − η y c y + c z ) = z = η . (15) This and the following boundary conditions are expressed in terms of ensemble mean quantities.Consequently, terms in the mean of the products of fluctuations might appear and a closure forthem invoked (Cao 2014, § On the mean bed, the upward net flux across the mean bed comes from theentrainment due to the fluid turbulence, that is ν ( b x c x + b y c y − c z ) = ( b x tan θ + ) w f c ae on z = b , (16)where, following the work of van Rijn (1984), the equilibrium referenceconcentration is approximated in terms of shear velocity u f and mean particlesize d as c ae ≈ u f d . (17)The nondimensional pde s (2)–(5), together with boundary conditions (8)–(16)govern the dynamics of the turbulent flow and suspended sediment. In order to provide theoretical support for the model redction, we embed the surfaceconditions (13)–(14) in a family of conditions that modify the tangential stressesto have an artificial forcing proportional to the square of the local, free surface,velocity: ( − η x ) τ + η x ( τ − τ ) − η y ( τ + η x τ )= ( − γ ) √ c t ( + c u )( + c u ) u √ u + v on z = η = h + b , (18) ( − η y ) τ + η y ( τ − τ ) − η x ( τ + η y τ )= ( − γ ) √ c t ( + c u )( + c u ) v √ u + v on z = η = h + b . (19)When evaluated at parameter γ = γ =
0, and when the mean slope andthe lateral derivatives are negligible (tan θ = ∂ x = ∂ y = ν∂u/∂z = ν ( u/η ) and ν∂v/∂z = ν ( v/η ) . In thiscase, two neutral modes of the dynamics are the lateral shear flow ( u , v ) ∝ z/h + c u . Analogously we embed the sediment boundary conditions (15) and (16) in thefamily ν (− η x c x − η y c y + γ c c z ) + ( − γ c ) ν ch The Euler parameter of a toy problem suggests introducing a factor ( − γ/ ) into the left-handside of the tangential stress boundary conditions (13)–(14) in order to improve convergence inthe parameter γ when evaluated at the physically relevant γ = (cid:104) + ( − γ c ) w f (cid:105) ( + η x tan θ ) w f c = z = h + b , (20) ν [− b x c x − b y c y + ( − γ c ) c z ] + ( − γ c ) ν ch + (cid:104) + ( − γ c ) w f (cid:105) ( + b x tan θ ) w f c ae = z = b . (21)Upon setting the embedded parameter γ c = γ c = w f / O ( w f ) . Withoutsuch term, the model only conserves sediment to errors O ( w f ) . The reason forimplementing such high order errors is that the information about different fallingvelocities comes into the model in O ( w f ) terms.When the artificial parameter γ c =
0, the lateral gravity and lateral derivativesare negligible (tan θ = ∂ x = ∂ y = w f =
0, boundarycondition (20) requires the concentration c = ∂c∂z + ch = z = b .Together with the sediment pde (5), these imply a neutral mode of the sedimentdynamics is c ∝ − z/h when artificial parameter γ c = γ = tan θ = ∂ x = ∂ y = γ c = w f =
0, a four parameter subspace ofequilibria exists corresponding to some uniform lateral shear, turbulent mean, flow, ( u , v ) ∝ z/h + c u , some turbulent mean concentration c ∝ − z/h , on a fluid ofany constant fluid depth h . For large enough lateral length scales, these equilibriaoccur independently at each location x and y (Roberts 1988, 2008, e.g.) and hencethe subspace of equilibria are in effect parameterised by ¯ u ( x , y ) , ¯ v ( x , y ) , ¯ c ( x , y ) and h ( x , y ) .Provided we can treat lateral derivatives ∂ x and ∂ y as perturbing influences,that is provided solutions vary slowly enough in x and y , centre manifold the-orems (Roberts 1988, Chicone 2006, e.g.) assure us of three vitally importantproperties: • this subspace of equilibria are perturbed to a slow manifold, whereon theevolution is slow, that exists for a finite range of gradients ∂ x and ∂ y , andparameters γ , γ c , tan θ and w f , and which may be parameterised by thedepth-averaged lateral velocities ¯ u ( x , y , t ) and ¯ v ( x , y , t ) , the depth-averagedconcentration ¯ c ( x , y , t ) , and the local fluid depth h ( x , y , t ) ; • the slow manifold attracts solutions from all nearby initial conditions providedthe spectrum of the linearised dynamics is suitable;9 and that a formal power series in the parameters γ , γ c , tan θ , w f andgradients ∂ x and ∂ y approximate the slow manifold to the same order of erroras the order of the residuals of the governing differential equations.That is, the theorems support the existence, emergence, and construction of slowmanifold models such as (1). This support occurs in a finite domain in parameterspace, and we assume the finite domain is big enough to include interesting casesof the physically relevant γ = γ c = θ and falling velocity w f . The linear dynamics of the system support the application of centre manifoldtheory. For the flat bed of b = constant, and with tan θ = γ = γ c = w f =
0, thebase problem (2)–(21) has the equilibrium of a shear flow which is, in terms of thestretched vertical coordinate Z = ( z − b ) /h , h = constant , u = u Z + c u + c u , v = v Z + c u + c u , w = p = h ( − Z ) , c = c ( − Z ) , ν = c t h ˙ ε = c t h ¯ q √ + c u , (23)where ¯ u and ¯ v are the depth-averaged lateral velocities, ¯ c is the depth-averagedconcentration, and ¯ q = √ ¯ u + ¯ v is the mean speed. Environmental turbulentflows have eddy viscosity variations. In this linear analysis, we assume the eddyviscosity ν is effectively constant.Then we consider the dynamics of the pde s (2)–(5) linearised in the smallperturbation fields ( h (cid:48) , u (cid:48) , v (cid:48) , w (cid:48) , c (cid:48) , p (cid:48) ) about each of these equilibria. Becauseenvironmental turbulent fluids have very large lateral scales compared with thedepth the lateral variations are very slow. As the lateral variations vary slowlythey do not affect the dominant linear process. Thus we also treat the lateralderivatives ∂ x and ∂ y as ‘small operators’ in this linearisation. Thus, the pde s (2)–(5) and the boundary conditions (8)–(21), in the linearisation that effectively ∂ x = ∂ y = tan θ = γ = γ c = w f =
0, result in the linear problem ∂w (cid:48) ∂z = ∂u (cid:48) ∂t + w (cid:48) ∂u∂z = ν ∂ u (cid:48) ∂z (24b) ∂v (cid:48) ∂t + w (cid:48) ∂v∂z = ν ∂ v (cid:48) ∂z , (24c) ∂w (cid:48) ∂t = − ∂p (cid:48) ∂z + ν ∂ w (cid:48) ∂z , (24d) This general type of argument has recently been made rigorous in one lateral dimension byRoberts (2013). c (cid:48) ∂t + w (cid:48) ∂c∂z = ν ∂ c (cid:48) ∂z , (24e) w (cid:48) = Z = u (cid:48) = c u h (cid:48) ∂u∂z + c u h ∂u (cid:48) ∂z on Z = v (cid:48) = c u h ∂v (cid:48) ∂z + c u h (cid:48) ∂v∂z on Z = ∂h (cid:48) ∂t = w (cid:48) on Z = − p (cid:48) + ν ∂w (cid:48) ∂z = Z = ν ∂u (cid:48) ∂z = √ c t ¯ q ( + c u )( + c u ) u (cid:48) on Z = ν ∂ν (cid:48) ∂z = √ c t ¯ q ( + c u )( + c u ) v (cid:48) on Z = c (cid:48) = Z = ∂c (cid:48) ∂z + c (cid:48) h = Z = w (cid:48) =
0. Equa-tion (24i) implies the free surface perturbation h (cid:48) = constant, which correspondsto the freedom already in (22) so without loss of generality we here set h (cid:48) = p (cid:48) = pde s (24b), (24c) and (24e) indicate fields u (cid:48) , v (cid:48) and c (cid:48) have solutions in the form ( u (cid:48) , v (cid:48) , c (cid:48) ) ∝ [ A sin ( kz ) + B cos ( kz )] exp ( λt ) , (25)where k is a nondimensional vertical wavenumber. Substitute these solutionforms (25) into pde s (24b), (24c) and (24e), and obtain λ = − νk . Substitutingthe solution forms (25) into boundary conditions (24k)–(24n) leads to two separateconditions for the velocity fields and concentration fields, respectivelytan k = k + c u ( + c u ) k and tan k = hk . (26)The characteristic equations (26) have the non-zero wavenumbers k > π , whichimplies the leading non-zero eigenvalue − νk < − νπ . In addition to these negativeeigenvalues, there are zero eigenvalues corresponding to the freedom to vary the fluiddepth h , depth-averaged velocities ¯ u and ¯ v and depth-averaged concentration ¯ c .Thus, there is a spectral gap between the eigenvalues λ = λ < − νπ . Centremanifold theory (Roberts 1988, Potzsche & Rasmussen 2006, e.g.) then supportsthe existence and emergence of a slow manifold of large lateral scale in the threedimensional turbulent fluid and sediment system.11 .5 Reduced model of the flow and sediment dynamics This section focusses on interpreting the slow manifold of the leading order modelof the turbulent flow and suspended sediment.Instead of depth-averaging equations, the centre manifold approximation theo-rem underlies an iterative construction of a slow manifold that resolves the turbulentand sediment interactions within the fluid layer. The computer algebra programlisted in Appendix A constructs the slow manifold of the turbulent flow and sedimentsystem. The program derives evolution equations for the water depth h ( x , y , t ) , thedepth-averaged lateral velocities ¯ u ( x , y , t ) and ¯ v ( x , y , t ) , and the depth-averagedconcentration ¯ c ( x , y , t ) .The order of error in the construction is phased in terms of the small parameters.Here the small parameters are the lateral derivatives ∂ x and ∂ y , the small meanslope tan θ , the falling velocity w f , and the artificial parameters γ and γ c . Generallywe report results to errors O ( ∂ p/ x + ∂ p/ y + tan p/ θ + w pf + γ p ) for some prescribedexponent p (where, for example, a term with factor w mf γ n is O (cid:0) w pf + γ q (cid:1) if m/p + n/q (cid:62) γ c is introduced tointegrate the sediment dynamics into the theoretical support for a slow manifold.In order to ensure the sediment dynamics are accurate we construct the slowmanifold to higher orders in the parameter γ c .Cao (2014) [ § γ converge quicklyin the turbulent fluid flow systems. Computations with the sediment equationsindicate that the dependence upon γ c also converges reasonably quickly. Truncatingto errors O ( ∂ / x + ∂ / y + tan / θ + γ + w f , γ c ) , the computer algebra in Appendix Aderives the evolution of the depth-averaged concentration ¯ c ( x , y , t ) : ∂ ¯ c∂t = · · · − (cid:18) ¯ u ∂ ¯ c∂x + ¯ v ∂ ¯ c∂y (cid:19) (cid:0) + γ c + γ c + γ c − γ c (cid:1) (27a) + ¯ ch (cid:18) ¯ u ∂h∂x + ¯ v ∂h∂y (cid:19) ( − γ c − γ c − γ c + γ c ) (27b) + O ( ∂ / x + ∂ / y + tan / θ + γ + w f , γ c ) .These power series converge quickly enough in γ c to reasonably evaluate at γ c = γ = γ c = h ( x , y , t ) , depth-averaged lateral velocities ¯ u ( x , y , t ) and ¯ v ( x , y , t ) and depth-averaged concentration ¯ c ( x , y , t ) . The equations hereare complicated due to the methodology systematically resolving all the intricate12icroscale physical interactions. ∂h∂t ≈ − (cid:18) ∂h ¯ u∂x + ∂h ¯ v∂y (cid:19) , (28a) ∂ ¯ u∂t ≈ − u ¯ qh + (cid:20) tan θ − ∂ ( h + b ) ∂x (cid:21) − u ∂ ¯ u∂x − v ∂ ¯ u∂y − (cid:18) ¯ u h ∂h∂x − ¯ u ¯ vh ∂h∂y (cid:19) + qh (cid:20) ∂∂x (cid:18) h ∂ ¯ u∂x (cid:19) + ∂∂y (cid:18) h ∂ ¯ u∂y (cid:19)(cid:21) + u − ¯ v h ¯ q (cid:20) ∂∂x (cid:18) h ∂ ¯ u∂x (cid:19) − ∂∂y (cid:18) h ∂ ¯ u∂y (cid:19)(cid:21) + ( s − ) ¯ u ¯ c ¯ qh − ( s − ) h ∂ ¯ c∂x , (28b) ∂ ¯ v∂t ≈ − v ¯ qh − ∂ ( h + b ) ∂y − v ∂ ¯ v∂y − u ∂ ¯ v∂x − (cid:18) ¯ u ¯ vh ∂h∂x − ¯ v h ∂h∂y (cid:19) + qh (cid:20) ∂∂x (cid:18) h ∂ ¯ v∂x (cid:19) + ∂∂y (cid:18) h ∂ ¯ v∂y (cid:19)(cid:21) + u − ¯ v h ¯ q (cid:20) ∂∂x (cid:18) h ∂ ¯ v∂x (cid:19) − ∂∂y (cid:18) h ∂ ¯ v∂y (cid:19)(cid:21) + ( s − ) ¯ v ¯ c ¯ qh − ( s − ) h ∂ ¯ c∂y , (28c) ∂ ¯ c∂t ≈ − w f h (cid:18) + w f ¯ q (cid:19) ¯ c + w f h (cid:18) − w f ¯ q (cid:19) c ae − (cid:18) − w f ¯ q (cid:19) (cid:18) ¯ u ∂ ¯ c∂x + ¯ v ∂ ¯ c∂y (cid:19) + qh (cid:20) ∂∂x (cid:18) h ∂ ¯ c∂x (cid:19) + ∂∂y (cid:18) h ∂ ¯ c∂y (cid:19)(cid:21) + u − ¯ v h ¯ q (cid:20) ∂∂x (cid:18) h ∂ ¯ c∂x (cid:19) − ∂∂y (cid:18) h ∂ ¯ c∂y (cid:19)(cid:21) (28d)Equation (28a) is a direct consequence of conservation of fluid. The momentumequations (28b)–(28c) include the effects of drag ¯ u ¯ q/h , advection, such as ¯ u∂ ¯ u/∂x ,turbulent dissipation, gravitational forcing tan θ − ∂ ( h + b ) /∂x , and pressure gra-dients established by the suspended sediment h∂ ¯ c/∂x . The sediment concentrationequation (28d) contains the equilibration of vertical sediment distribution due toterms such as w f ¯ c/h , advection such as ¯ u∂ ¯ c/∂x , and turbulent dispersion effects.Although equations (28a)–(28d) are expressed in terms of depth-averaged lat-13ral velocities and depth-averaged concentration, they are derived not by depth-averaging, but instead by systematically accounting for interaction between verticalprofiles of the velocity and concentration, the stress, bed drag, lateral space varia-tions and bed topography. The form of coefficients in equations (28a)–(28d) aresupported by dynamical systems theory: the detail in the equations reflects thata slow manifold is in principle composed of exact solutions of the Smagorinskidynamics (3) and convection-diffusion equation (5), and hence accounts for allinteractions up to a given order of analysis no matter how small the numericalcoefficient in the interactions.The sediment model (28d) includes all the dominated terms in established mod-elling (Wu 2004, Duan 2004, Duan & Nanda 2006, e.g.). For example, Duan (2004)derived the following depth-averaged advection-diffusion equation of suspendedsediment: ∂ ¯ c∂t = − w f h ( ¯ c − c ae ) − ¯ u ∂ ¯ c∂x − ¯ v ∂ ¯ c∂y + h ¯ q (cid:18) ∂ ¯ c∂x + ∂ ¯ c∂y (cid:19) , (29)which consists of the effects of vertical distribution w f /h ( ¯ c − c ae ) , advection¯ u∂ ¯ c/∂x and ¯ v∂ ¯ c/∂y , and dispersion h ¯ q∂ ¯ c/∂x and h ¯ q∂ ¯ c/∂y . But the de-rived model (28d) also includes more subtle effects, which could be importantfor suspended sediment in complex flow regimes. The model (28d) further givesmodifications in presence of the ratio w f / ¯ q due to different distributions of sedi-ment in the vertical for the different levels of turbulent mixing characterised bythe mean flow speed ¯ q . The coefficients in equation (28d) are a little differentto the established model (29). Take the advection term ¯ u∂ ¯ c/∂x for example,the established model (29) has the coefficient 1, whereas our coefficient of suchterms is 1.01 − w f / ¯ q . Physically, a higher falling velocity w f means sedimentconcentrates more near bed where the mean advection velocity is lower, and hencenet transport will be slower. Our model (28d) has smaller dispersion coefficentcompared with the model (29). The centre manifold approach does not impose a specific cross-sectional velocitydistribution as done by other methods, instead our approach empowers the sedimentand Smagorinski turbulent equations to determine the cross-sectional structures inout-of-equilibrium dynamics. Recall that the locally stretched vertical coordinate Z = ( z − b ) /h . This section concentrates on the vertical distribution of the lateralvelocity u ( Z ) and concentration c ( Z ) in steady flow, and compare our predictionof the lateral velocity u ( Z ) with analogous published experimental data (Schultz& Flack 2007, 2013, Celik & Rodi 1988, e.g.).14 .1 Distribution of the suspended sediment We concentrate on the vertical distribution of the suspended sediment c ( Z ) inthe slow manifold of the system. Computer algebra in Appendix A derives thephysical field of sediment concentration c in terms of the depth h ( x , y , t ) , depth-averaged velocities ¯ u ( x , y , t ) and ¯ v ( x , y , t ) , depth-averaged concentration ¯ c ( x , y , t ) and stretched local normal coordinate Z = ( z − b ) /h on the slow manifold, evalu-ating at γ = γ c = c ( Z ) = ¯ c (cid:0) + Z − Z − Z (cid:1) (30a) + ¯ c w f ¯ q (cid:0) − Z − Z (cid:1) (30b) + ¯ c w f ¯ q (cid:0) − − Z + Z − Z (cid:1) (30c) + c ae w f ¯ q (cid:0) − Z + Z + Z (cid:1) (30d) + c ae w f ¯ q (cid:0) − + Z − Z (cid:1) (30e) + h ¯ q (cid:18) ¯ u ∂ ¯ c∂x + ¯ v ∂ ¯ c∂y (cid:19) (cid:0) + Z − Z + Z (cid:1) (30f) + h ¯ c ¯ q ∂ ¯ u∂x (cid:0) − + Z − Z − Z (cid:1) (30g) + h ¯ c ¯ q ∂ ¯ v∂y (cid:0) − + Z + Z − Z (cid:1) (30h) + ¯ c ¯ q (cid:18) ¯ u ∂h∂x + ¯ v ∂h∂y (cid:19) (cid:0) − + Z − Z − Z (cid:1) (30i) + ¯ c ¯ q (cid:18) ¯ u ∂b∂x + ¯ v ∂b∂y (cid:19) (cid:0) − Z .549 Z + Z (cid:1) (30j) + ¯ c ¯ q (cid:18) ¯ u ∂b∂x + ¯ v ∂b∂y (cid:19) (cid:0) − Z .549 Z + Z (cid:1) (30k) − tan θ h ¯ u ¯ c ¯ q (cid:0) − Z + Z + Z (cid:1) (30l) + O (cid:0) ∂ / x + ∂ / y + tan / θ + w f + γ , γ c (cid:1) .The vertical sediment distribution (30) describes the low-order shape of the slowmanifold in state space. Physically this equation describes the details of thesuspended sediment concentration in out-of-equilibrium flow. The terms in equa-tion (30) have physical interpretations. For example, the line (30a) approximatesthe mean concentration, together with the lines (30b)–(30e), which forms thebasic distribution of the concentration c ( Z ) in the vertical in the presence of thedepth-averaged concentration ¯ c , the equilibrium bed concentration c ae , the fallingvelocity w f , and the mean fluid speed ¯ q . The lines (30f)–(30i) describe the effect15y advection on the vertical distribution of the concentration. Lines (30j)–(30k)describe modifications due to the change of the bed topography. The line (30l)describes modifications due to lateral flow down a slope in its affect on the verticaldistribution of the suspended sediment.Now we explore the distribution of the suspended sediment c ( Z ) in steady flows.Consider the turbulent flow of constant depth H = θ ; that is, the mean bed b =
0. Weconsider the concentration fraction is small, ¯ c < u ¯ c ¯ q/h and ¯ v ¯ c ¯ q/h are negligible. Then the evolution equations (28a)–(28c) predicts the equilibrium U = / θ and V =
0, so the mean speed¯ q = U = / θ .The nondimensional equilibrium reference concentration c ae on the mean bed z = b in steady flow is determined by the particle size and the slope. We ap-proximate the nondimensional shear velocity u f = ¯ q/C (cid:48) , where C (cid:48) =
18 log ( /d ) is the Chezy coefficient (van Rijn 1984, e.g.). Thus, equation (17) gives thenondimensional equilibrium reference concentration c ae = θ d ( − log d ) , (31)which only depends on the nondimensional mean slope tan θ and the nondimensionalmean particle size d .For the steady flow, the concentration c ( Z ) in equation (30) then reduces to c ( Z ) = ¯ c (cid:0) + Z − Z − Z (cid:1) + ¯ c w f ¯ q (cid:0) − Z − Z (cid:1) + ¯ c w f ¯ q (cid:0) − − Z + Z − Z (cid:1) + c ae w f ¯ q (cid:0) − Z + Z + Z (cid:1) + c ae w f ¯ q (cid:0) − + Z − Z (cid:1) . (32)Equation (32) shows that the concentration c ( Z ) depends on the vertical coordi-nate Z , the falling velocity w f , the depth-averaged concentration ¯ c , the equilibriumreference concentration c ae and the mean flow speed ¯ q . The falling velocity w f andthe equilibrium reference concentration c ae vary with the mean particle size d andthe mean slope tan θ according to equation (7) and (31), respectively. Thus, theconcentration profile c ( Z ) depends on the vertical coordinate Z , the mean particlesize d and the mean slope tan θ .Figure 3 depicts the profiles of the nondimensional suspended sediment con-centration c ( Z ) in the vertical for four different nondimensional mean particlesize d . The bed has a mean slope tan θ = (cid:239) c ( Z ) d = × − d = × − d = × − Figure 3: Profiles (line curves) of the suspended sediment concentration c ( Z ) inthe vertical for three different nondimensional mean particle size d according toequation (32). The mean slope tan θ = c ae = c ae ≈ d > × − . Forsmall nondimensional mean particle size d , the nondimensional concentration c ( Z ) is approximately linear in the vertical coordinate Z . In this steady flow, theapproximation of the eddy diffusivity is (cid:15) s ( Z ) ≈ ¯ q ( − Z − Z )+ tan θ q ( − Z + Z ) . (33)For an indicative comparison, we integrate equation (5) from the bottom to thefree surface in the steady flow, and obtain an approximation c ( Z ) ≈ c ae (cid:20) + Z ( − Z ) (cid:21) − w f / ¯ q . (34)The dash curves in Figure 3 depict the approximation (34). When the nondimen-sional mean particle size d is small, the computed vertical distribution (32) is17 ¯ c/c ae Figure 4: Vertical distribution of the suspended sediment: (blue curve) fromequation (32); (red circle) the numerical prediction by Celik & Rodi (1988); and(green stars) the corresponding experimental data used by Celik & Rodi (1988).The nondimensional mean particle size d = × − and then the fallingvelocity w f = d increases, there is a difference between the computed verticaldistribution (32) and the approximation (34) at the upper flow.Figure 4 plots the vertical distribution of the suspended sediment from equa-tion (32) (blue curve), the numerical prediction (red circles) by Celik & Rodi (1988),and the corresponding experimental data (green stars) used by Celik & Rodi (1988).Celik & Rodi (1988) calculated the suspended sediment transport in unidirectionalchannel flow, where they used the nondimensional variables of fluid depth H = U ≈ d = × − and falling velocity w f = .2 Vertical distribution of the velocity in steady flow This section reports on the vertical distribution of the lateral velocity u ( Z ) inthe flow and sediment system, v ( Z ) is similar. Computer algebra in Appendix Aderives the following physical flow field of lateral velocity u in term of stretchedlocal normal coordinate Z = ( z − b ) /h on the slow manifold, evaluated at thephysical γ = γ c = u ( Z ) = ¯ u ( + Z − Z − Z − Z − Z )+ tan θ h ¯ q ( + Z − Z + Z − Z + Z + Z + Z )+ h ¯ u ¯ q ∂ ¯ u∂x ( + Z − Z + Z + Z + Z + Z + Z + Z )+ h ¯ v ¯ q ∂ ¯ u∂y ( + Z − Z + Z + Z + Z + Z + Z + Z )+ h ¯ q (cid:18) ∂h∂x + ∂b∂x (cid:19) (− − Z + Z − Z + Z − Z − Z − Z ) , + ( s − ) ¯ c ¯ u (cid:0) + Z − Z + Z + Z (cid:1) + O (cid:0) ∂ / x + ∂ / y + tan / θ + w f + γ , γ c (cid:1) . (35)Physically, equation (35) describes the vertical details of the lateral velocity u ( Z ) in terms of the vertical position Z , fluid depth h , depth-averaged velocities ¯ u and ¯ v ,depth-averaged concentration ¯ c , and the slope tan θ of the mean bed b .Schultz & Flack (2013) experimentally studied the smooth-wall turbulent chan-nel flow with the Reynolds number Re up to 300 000, and showed that the meanflow is approximately independent of the Reynolds number. In earlier work, Schultz& Flack (2007) showed when the relative roughness, the ratio of the roughnessheight and the boundary-layer thickness, is small, the mean velocity shape for therough and smooth walls are similar in the outer layer.However, compared with large roughness in the environmental flows underconsideration, the roughnesses in the experiments of Schultz & Flack (2007) aresmall—the ratio between the roughness height and the fluid depth was ≈ × − .To compare with the experimental data (Schultz & Flack 2007, 2013, e.g.), wederive the equilibria profiles and evaluate the lateral velocity at these equilibria.Consider the turbulent flow with suspended sediment flowing on a flat meanbed of constant slope tan θ ; that is, the mean bed b =
0. Consider the suspendedsediment in steady flow of depth H =
1. We consider the concentration fraction19 u ( Z ) /u ( ) Figure 5: Vertical distribution of the lateral velocity from: (blue curve) theapproximation (36); and (red circle curve) the experimental measurements (Schultz& Flack 2013, Fig. 2). The ratio u ( Z ) /u ( ) is independent of the slope tan θ ,where u ( ) is the lateral velocity at the level Z = c ae , ¯ c < u ¯ c ¯ q/h and ¯ v ¯ c ¯ q/h arenegligible. Then the evolution equations (28a)–(28c) predicts the equilibrium U = / θ and V =
0, so the mean speed ¯ q = U = / θ . For thissteady flow, the lateral velocity (35) reduces to u ( Z ) (cid:112) tan θ/c t ≈ + Z − Z − Z − Z − Z − Z − Z . (36)Figure 5 compares the profile of the lateral velocity from the approximation (36) (blueline curve) with the experimental measurements (red circle curve) by Schultz &Flack (2013). From the approximation (36), the velocity ratio u ( Z ) /u ( ) is indepen-dent of the slope tan θ . Our prediction of the lateral velocity in the vertical agreesreasonably the experimental lateral velocity by Schultz & Flack (2013), except nearthe mean bed. In environmental flows we expect that the large roughness of stones,roots and debris to typically break up any log boundary layer. Thus, we do notresolve the turbulent log layer, and are interested in the dynamics determined bythe relatively large scale of the fluid depth.Shear stress arises in the turbulent fluid. In this steady flow, we predict the20 x z ( Z ) / t a n θ Z Figure 6: Shear stress profile of: (blue curve) our approximation τ z ; and (redtriangle curve) the experimental measurements of flow over rough bed (Schultz &Flack 2007, Fig. 9). The error bars show the ±
9% uncertainty in their experiments.shear stress τ xz has the near linear profile of τ xz ( Z ) tan θ ≈ − Z + Z − Z + Z + Z + Z . (37)Figure 6 shows that our approximation of shear stress τ xz (blue curve) approachesthe experimental measurements (red triangle curve) by Schultz & Flack (2007).The shear stress is approximately straight due to the need to dissipate the nearconstant forcing of gravity. A difference occurs near the bed, because we do notresolve the log boundary layer. This section explores the depth-averaged concentration in waves on an inclinedrippled bed by the suspended sediment model (28d), coupled with the modifiedmomentum equations (28a)–(28c). Numerical results are qualitatively comparedwith the experimental measurements of the suspended sediment (Zedler & Street2006, Kos’Yan et al. 2007, e.g.). 21 e l o c i t y ¯ u (cid:239) c o n ce n tr a t i o n ¯ c t Figure 7: Time series of the depth-averaged lateral velocity ¯ u (blue curves) and thedepth-averaged concentration ¯ c (red curves) at the trough x =
50 (line curves) andat the crest x =
60 (dash curves) in Figure 8. The mean particle size d = × − ,so the falling velocity w f = c ae = h ( x , y , t ) with suspended sediment of depth-averagedconcentration ¯ c ( x , y , t ) flowing down a rippled bed. The fluid has the depth-averaged lateral velocities ¯ u ( x , y , t ) and ¯ v ( x , y , t ) along the bed. Let the bed havenondimensional length L x =
100 and width L y =
10. The mean slope of the bedis tan θ = v = H =
1, then the mean equilibrium depth-averaged lateral velocity U ≈ c ≈ ( π/L x x ) to the equilibrium.Recall the Froude number U/ √ gH = > x and y directions for both the flow and bed. Figure 7 plots the time22 & ¯ u (cid:239) x Figure 8: Plots of the depth h (black) and depth-averaged velocity ¯ u (blue) of thefluid flowing over the rippled bed (green) at time t =
180 in Figure 7. The dashline represents the zero mean bed. This figure shows supercritical flow arises.series of the depth-averaged lateral velocity ¯ u (blue curves) and the depth-averagedconcentration ¯ c (red curves) at a trough x =
50 (line curves) and at a crest x =
60 (dash curves). The periodic depth-averaged velocity ¯ u indicates that largeroll waves are generated on the free surface (Balmforth & Mandre 2004, e.g.).The depth-averaged velocity ¯ u is bigger at the trough x =
50 (blue line curve)than at the crest x =
60 (blue dash curve), which indicates enhanced turbulentmixing arises at the trough (Zedler & Street 2001, e.g.). Then the turbulent mixingproduces slightly bigger depth-averaged concentration ¯ c at the trough x =
50 (redline curve) than at the crest x =
60 (red dash curve). Zedler & Street (2006), intheir calculation of suspended sediment over rippled beds, found that far from thebed the concentration at the crest and trough are approximately π out of phase.They commented that this phase lag is due to the vortex produced by the ripplenear the trough. Kos’Yan et al. (2007) experimentally and mathematically plottedthe time series of the horizontal velocity and concentration under waves over sandybottom, which showed that the concentration slightly lags the horizontal velocityto reach maximum. However, in our simulation, no significant lag happens whichappears due to our ripple not producing a vortex near the trough.Figure 8 plots the water depth h (black) and depth-averaged velocity ¯ u (blue)of the fluid flowing over the rippled bed (green) in the x direction at time t = e p t h h (cid:239) c o n ce n tr a t i o n ¯ c x Figure 9: Plots of the depth h (black) and the depth-averaged concentration ¯ c (red)in the x direction at time t =
180 in Figure 7. The depth-averaged concentration ¯ c is ahead to reach maximum over a ripple.exhibits the supercritical flow as the fluid flowing over each ripple on the bed. Thedepth h rises at the crest and the depth-averaged velocity ¯ u declines at the crest,which corresponds to the depth-averaged velocity ¯ u reaches minimum at the crestin Figure 7. Figure 9 plots the depth h and the depth-averaged concentration ¯ c in x -direction at time t =
180 in Figure 7. The depth-averaged concentration ¯ c isapproximate π/ h . That is because the strong turbulentmixing at the troughs makes the concentration peak quickly.There is only one significant peak in one period in Figure 8–9. Zedler &Street (2006), who reported numerical results of large eddy simulation of the flowand suspended sediment over sinusoidal ripples, found three peaks on the timeseries of concentration at the crest and trough. Their ripples have a height towavelength ratio of 0.1, which is five times steeper than our ripples. Zedler & Street(2006) commented that these peaks are mainly due to the vortex, shear stress andadvection near the ripple. In our simulation, the only significant peak in a periodis due to the turbulent mixing.Figure 10 compares at the time t =
180 the depth h (black) and depth-averaged concentration ¯ c (red) of the flow over ripples with different heights. Thedash curves represent the depth h and depth-averaged concentration ¯ c for theripple height 0.4, while the line curves are for the ripple height 0.6. The fluid24 e p t h h (cid:239) c o n ce n tr a t i o n ¯ c x Figure 10: Plots of the depth h (black curves) and depth-averaged concen-tration ¯ c (red curves) of the fluid flowing over the rippled bed with rippleheight 0.4 (dash curves) and 0.6 (line curves) in the x direction at time t = h is usually bigger at the crest of the steeper ripple, but the depth-averagedconcentration ¯ c becomes smaller for the steeper ripples. However, the laboratoryexperiment by Osborne & Vincent (1996), whose ripples have an approximate heightto wavelength ratio of 0.2, verifies that steep asymmetric ripples under shoalingwaves produce greater concentrations higher in the water column than low steepnessripples. Such difference is possibly because the ripples with small height in oursimulation do not produce strong vortices that enhance the pick up of sedimentinto suspension. The phenomenon of smaller depth-averaged concentration forsteeper ripples is because the increased fluid depth for steeper ripples producessmall depth-averaged concentration according to the erosion and deposition w f ¯ c/h in the governing equation (28d). This work derives a suspended sediment model (28a)–(28d) to simulate the in-teractions between suspended sediment and turbulent flows. The concentrationequation (28d) consists of the effects of sediment erosion, advection, and dispersion.Section 2.3 embedded the physical boundary conditions on the free surface andon the mean bed in a family of problems to access a slow manifold in the system.25he parameter γ = γ c = γ = γ c = tan θ = w f = ∂ x = ∂ y =
0, a four parameter familyof equilibria exists to support the existence of a slow manifold in the system, aslow manifold describing the large lateral structures. Computer algebra detailedin Appendix A leads to the evolution equations in the field of depth h ( x , y , t ) ,depth-averaged lateral velocities ¯ u ( x , y , t ) and ¯ v ( x , y , t ) , and depth-averaged con-centration ¯ c ( x , y , t ) . It is reassuring that the dominant terms in our model (28d)agree with established modelling (Wu 2004, Duan 2004, Duan & Nanda 2006,e.g.). Then our model includes more subtle effects, that could be important forsuspended sediment in complex flow regimes. The trends of the suspended sedimentconcentration corresponds to the published experimental measurements (Cellino &Graf 1999, Yoon & Kang 2005, e.g.).Section 4 implemented numerical simulations of the suspended sediment underlarge waves by the suspended sediment model (28a)–(28d). The time series ofthe depth-averaged suspended sediment concentration rises fast and falls slowly.Supercritical flow arises when the fluid flowing over the rippled bed. The plots ofthe depth-averaged concentration ¯ c in space show that high concentration arises atthe troughs and low concentration arises at the crests. References
Balmforth, N. J. & Mandre, S. (2004), ‘Dynamics of roll waves’,
J. Fluid Mech. , 1–33. doi:10.1017/S0022112004009930 .Cao, M. (2014), Modelling environmental turbulent fluids and multiscale modellingcouples patches of wave-like system, PhD thesis, School of Mathematical Sciences,University of Adelaide.Cao, M. & Roberts, A. J. (2012), Modelling 3D turbulent floods based upon theSmagorinski large eddy closure. Proceedings of the 18th Australasian FluidMechanics Conference 3rd-7th December 2012 Edited by P. A. Brandner andB. W. Pearce Published by the Australasian Fluid Mechanics Society. .Celik, I. & Rodi, W. (1988), ‘Modeling suspended sediment transport in nonequi-librium situations’,
Journal of Hydraulic Engineering , 1157–1191. doi:10.1061/(ASCE)0733-9429(1988)114:10(1157) .Cellino, M. & Graf, W. H. (1999), ‘Sediment-laden flow in open channels un-der noncapacity and capacity conditions’,
Journal of Hydraulic Engineering (5), 455–462. doi:10.1061/(ASCE)0733-9429(1999)125:5(455) .Chanson, H. (2004),
Hydraulics of Open Channel Flow , 2 edn, Butterworth–Heinemann. 26hicone, C. (2006),
Ordinary Differential Equations with Applications , Vol. 34of
Texts in Applied Mathematics , Springer New York. doi:10.1007/0-387-35794-7 .Duan, J. G. (2004), ‘Simulation of flow and mass dispersion in meandering chan-nels’,
Journal of Hydraulic Engineering , 964–976. doi:10.1061/(ASCE)0733-9429(2004)130:10(964) .Duan, J. G. & Nanda, S. (2006), ‘Two-dimensional depth-averaged model simulationof suspended sediment concentration distribution in a groyne field’,
Journal ofHydrology (3–4), 426 – 437. doi:10.1016/j.jhydrol.2005.11.055 .Fredsoe, J. & Deigaard, R. (1992),
Mechanics of coastal sediment transport. , Vol. 3of
Advanced Series on Ocean Engineering , World Scientific Publication.Georgiev, D. J., Roberts, A. J. & Strunin, D. V. (2009), Modelling turbulent flowfrom dam break using slow manifolds, in G. N. Mercer & A. J. Roberts, eds, ‘Pro-ceedings of the 14th Biennial Computational Techniques and Applications Confer-ence, CTAC-2008’, Vol. 50 of
ANZIAM J. , pp. C1033–C1051. http://anziamj.austms.org.au/ojs/index.php/ANZIAMJ/article/view/1466 [September 3,2009].Hunt, J. N. (1954), ‘The turbulent transport of suspended sediment in openchannels’,
Proc. R. Soc. Lond. A (1158), 322–335. doi:10.1098/rspa.1954.0161 .Kos’Yan, R., Divinskiy, B., Krylenko, M. & Vincent, C. (2007), ‘Modelling of thevertical distribution of suspended sediment concentration under waves with agroup structure’,
Ocean 2007–Europe (art. no. 4302351). https://ueaeprints.uea.ac.uk/28240/ .Osborne, P. D. & Vincent, C. E. (1996), ‘Vertical and horizontal structure issuspended sand concentrations and wave-induced fluxes over bedforms’,
MarineGeology (3–4), 195–208. doi:10.1016/0025-3227(95)00002-X .Pittaluga, M. & Seminara, G. (2003), ‘Depth-integrated modeling of suspendedsediment transport’,
Water Resources Research (1137), 11. doi:10.1029/2002WR001306 .Potzsche, C. & Rasmussen, M. (2006), ‘Taylor approximation of integral manifolds’, Journal of Dynamics and Differential Equations (2), 427–460. doi:10.1007/s10884-006-9011-8 .Roberts, A. (2008), ‘The inertial dynamics of thin film flow of non-newtonian fluids’, Physics Letters A (10), 1607–1611. doi:10.1016/j.physleta.2007.10.014 .27oberts, A. J. (1988), ‘The application of centre-manifold theory to the evolutionof systems which vary slowly in space’,
J. Austral. Math. Soc. Ser. B , 480–500. doi:10.1017/S0334270000005968 .Roberts, A. J. (2013), Macroscale, slowly varying, models emerge from the mi-croscale dynamics in long thin domains, Technical report, [ http://arxiv.org/abs/1310.1541 ].Roberts, A. J., Georgiev, D. J. & Strunin, D. V. (2008), Model turbulent floodswith the Smagorinsky large eddy closure, Technical report. http://arxiv.org/abs/0805.3192 .Schultz, M. P. & Flack, K. A. (2007), ‘The rough-wall turbulent boundary layer fromthe hydraulically smooth to the fully rough regime’, Journal of Fluid Mechanics , 381–405. doi:10.1017/S0022112007005502 .Schultz, M. P. & Flack, K. A. (2013), ‘Reynolds-number scaling of turbulent channelflow’,
Physics of Fluids . http://dx.doi.org/10.1063/1.4791606 .van Rijn, L. C. (1984), ‘Sediment transport, part II: Suspended load transport’, Journal of Hydraulic Engineering (11), 1613–1641. doi:10.1061/(ASCE)0733-9429(1984)110:11(1613) .Wu, W. (2004), ‘Depth-averaged two-dimensional numerical modeling of unsteadyflow and nonuniform sediment transport in open channels’,
Journal of HydraulicEngineering (10), 1013–1024. doi:10.1061/(ASCE)0733-9429(2004)130:10(1013) .Wu, W., Rodi, W. & Wenka, T. (2000), ‘3D numerical modeling of flow andsediment transport in open channels’,
Journal of Hydraulic Engineering . doi:10.1061/(ASCE)0733-9429(2000)126:1(4) .Yoon, J.-Y. & Kang, S.-K. (2005), ‘A numerical model of sediment-laden turbulentflow in an open channel’, Canadian Journal of Civil Engineering (1), 233–240. doi:10.1139/l04-089 .Zedler, E. A. & Street, R. L. (2001), ‘Large-eddy simulation of sediment transport:currents over ripples’, J. Hydraulic Engineering (6), 444–452. doi:10.1061/(ASCE)0733-9429(2001)127:6(444) .Zedler, E. A. & Street, R. L. (2006), ‘Sediment transport over ripples in oscillatoryflow’,
J. Hydraul. Eng. (2), 180–193. doi:10.1061/(ASCE)0733-9429(2006)132:2(180) . 28
Ancillary computer algebra program
This appendix lists the computer algebra to construct the slow manifold model ofthe suspended sediment in turbulent flow.Denote the fluid depth h ( x , y , t ) by h , depth-averaged lateral velocities ¯ u ( x , y , t ) and ¯ v ( x , y , t ) by uu and vv , depth-averaged suspended sediment concentration¯ c ( x , y , t ) by cc , and their time derivatives by h t = gh , ¯ u t = gu , ¯ v t = gv and¯ c t = gc . Denote the mean bed b ( x , y ) by b . The coefficients of lateral and normalgravitational forcing are represented by grx , gry and grz:=1 . Use qq representthe mean flow speed ¯ q = √ ¯ u + ¯ v and rqq for the reciprocal of this mean speed.Use the operator h(m,n) to denote the various lateral derivatives of the fluiddepth, ∂ mx ∂ ny h . Similarly use the operators uu(m,n) , vv(m,n) , cc(m,n) to denotelateral derivatives of the depth-averaged lateral velocities ¯ u ( x , y , t ) and ¯ v ( x , y , t ) ,and the depth-averaged concentration ¯ c ( x , y , t ) . These operators depend upontime and lateral space. Then the lateral derivative ∂ x h ( m , n ) = ∂ m + ∂ ny h , and thetime derivative ∂ t h ( m , n ) = ∂ mx ∂ ny gh , for example. Define readable abbreviationsfor h ( x , y , t ) and its first spatial derivatives. We use d to count the number oflateral derivatives so we can easily truncate the asymptotic expansion.Define the operators for the mean flow speed ¯ q and its reciprocal. The lastsimplification rule for rqq breaks the symbolic symmetric between the two lateraldirections. However the benefit of canonical representation outweighs the cost ofloss of symbolic symmetry.The key to the correctness of this program is that the residuals of the governingequations are computed correctly, and that the algorithm only terminates whenthese residuals are zero to the specified error.
09 resc:=df(u,x)+df(v,y)+df(w,z);110 resa:=sub(zz=0,w-u*df(b,x)-v*df(b,y));111 write length_resc:={length(resc),length(resa)};112 ok:=if {resc,resa}={0,0} then ok else 0;113 w:=w+(dw:=-h*wsolv(resc,zz))-resa;114 ezz:=ezz+df(dw,zz)/h;115 tzz:=tzz+2*r2*ct/(1+2*cu)*qq*df(dw,zz);116 % update thickness evolution117 gh:=sub(zz=1, w-u*hx-v*hy-u*df(b,x)-v*df(b,y));118 %solve vertical momentum and normal stress119 resw:=df(w,t)+u*df(w,x)+v*df(w,y)+w*df(w,z)+exc*df(p,z)120 +grz+exc*(-df(txz,x)-df(txy,y)-df(tzz,z));121 restn:=sub(zz=1,-p*(1+(hx+df(b,x))^2+(hy+df(b,y))^2)122 +tzz-2*(hx+df(b,x))*txz-2*(hy+df(b,y))*tyz123 +(hx+df(b,y))^2*txx+2*(hx+df(b,x))124 *(hy+df(b,y))*txy+(hy+df(b,y))^2*tyy);125 write length_resw:={length(resw),length(restn)};126 ok:=if {resw,restn}={0,0} then ok else 0;127 % update the pressure128 p:=p+h*psolv(resw,zz)+restn;129 %Smagorinski large eddy closure130 exx:=df(u,x);131 eyy:=df(v,y);132 ezz:=df(w,z);133 exz:=(df(u,z)+df(w,x))/2;134 exy:=(df(u,y)+df(v,x))/2;135 eyz:=(df(v,z)+df(w,y))/2;136 rese:=exx^2+ezz^2+eyy^2+2*exz^2+2*exy^2+2*eyz^2-ros^2;137 write length_rese:=mylength(rese);138 ok:= if rese=0 then ok else 0;139 ros:=ros+rese*h*(cu+1/2)/r2*rqq;140 txx:=2*ct*h^2*ros*exx;141 tyy:=2*ct*h^2*ros*eyy;142 tzz:=2*ct*h^2*ros*ezz;143 txz:=2*ct*h^2*ros*exz;144 txy:=2*ct*h^2*ros*exy;145 tyz:=2*ct*h^2*ros*eyz;146 %solve lateral momentum147 resu:=df(u,t)+u*df(u,x)+v*df(u,y)+w*df(u,z)+exc*df(p,x)148 -grx+exc*(-df(txx,x)-df(txy,y)-df(txz,z));149 resv:=df(v,t)+u*df(v,x)+v*df(v,y)+w*df(v,z)+exc*df(p,y)150 -gry+exc*(-df(tyy,y)-df(txy,x)-df(tyz,z));151 resbu:=sub(zz=0,(-u-w*df(b,x))*(1-df(b,x)^2/2)152 +cu*h*(1-df(b,x)^2/2-df(b,y)^2/2)*(153 -df(u+w*df(b,x),x)*df(b,x)
54 -df(u+w*df(b,x),y)*df(b,y)155 +df(u+w*df(b,x),z)));156 resbv:=sub(zz=0,(-v-w*df(b,y))*(1-df(b,y)^2/2)157 +cu*h*(1-df(b,x)^2/2-df(b,y)^2/2)*(158 -df(v+w*df(b,y),x)*df(b,x)159 -df(v+w*df(b,y),y)*df(b,y)160 +df(v+w*df(b,y),z)));161 write length_resuv:={length(resu),length(resv)162 ,length(resbu),length(resbv)};163 ok:=if {resu,resv,resbu,resbv}={0,0,0,0} then ok else 0;164 % update lateral mean velocities165 u:=u+(du:=resbu*(1-2*zz)/(1+2*cu)166 +h*(1+2*cu)*r2/(4*ct)*rqq^3*usolv(167 +(qq^2+vv(0,0)^2)*resu168 -uu(0,0)*vv(0,0)*resv ,zz));169 v:=v+(dv:=resbv*(1-2*zz)/(1+2*cu)170 +h*(1+2*cu)*r2/(4*ct)*rqq^3*usolv(171 +(qq^2+uu(0,0)^2)*resv172 -uu(0,0)*vv(0,0)*resu ,zz));173 ros:=ros+(uu(0,0)*df(du,zz)+vv(0,0)*df(dv,zz))/(r2*h)*rqq;174 % compute tangential stresses175 exz:=(df(u,z)+df(w,x))/2;176 eyz:=(df(v,z)+df(w,y))/2;177 txz:=2*ct*h^2*ros*exz;178 tyz:=2*ct*h^2*ros*eyz;179 resttu:=(-sub(zz=1,180 (1-0*gamm)*((1-(hx+df(b,x))^2)*txz181 +(hx+df(b,x))*(tzz-txx)182 -(hy+df(b,y))*(txy+(hx+df(b,x))*tyz))183 -(1-gamm)*r2*ct/(cu+1)/(2*cu+1)*u*qq) );184 resttv:=(-sub(zz=1,185 (1-0*gamm)*((1-(hy+df(b,y))^2)*tyz186 +(hy+df(b,y))*(tzz-tyy)187 -(hx+df(b,x))*(txy+(hy+df(b,y))*txz))188 -(1-gamm)*r2*ct/(cu+1)/(2*cu+1)*v*qq) );189 write length_restt:={length(resttu),length(resttv)};190 ok:=if {resttu,resttv}={0,0} then ok else 0;191 % update lateral evolutions192 gu:=gu-3*(1+2*cu)*(1+cu)193 /2/h/(3+11*cu+12*cu^2)/(1+3*cu+3*cu^2)/(3+4*cu)194 *(((1+5*cu+8*cu^2)*uu(0,0)^2*rqq^2195 -(9+45*cu+80*cu^2+48*cu^3))*resttu196 +(1+5*cu+8*cu^2)*uu(0,0)*vv(0,0)*rqq^2*resttv);197 gv:=gv-3*(1+2*cu)*(1+cu)198 /2/h/(3+11*cu+12*cu^2)/(1+3*cu+3*cu^2)/(3+4*cu)
99 *(((1+5*cu+8*cu^2)*vv(0,0)^2*rqq^2200 -(9+45*cu+80*cu^2+48*cu^3))*resttv201 +(1+5*cu+8*cu^2)*uu(0,0)*vv(0,0)*rqq^2*resttu);202 % solve the suspended sediment concentration.203 es:=ct*h^2*ros;204 sed:=df(c,t)+df(c*u,x)+df(c*v,y)+df(w*c,z)205 -df(wf*c,z)*the+grx*df(wf*c,x)206 -df(es*df(c,x),x)-df(es*df(c,y),y)-df(es*df(c,z),z);207 % boundary conditions on the free surface and bed.208 sedf:=sub(zz=1,209 (1+(1-gamc)*wf/6)*wf*c*(1+grx*df(h+b,x))210 +es*(-df(h+b,x)*df(c,x)211 -df(h+b,y)*df(c,y)212 +df(c,z)*gamc) +2*(1-gamc)*es*c/h);213 sedb:=sub(zz=0,214 (1+(1-gamc)*wf/6)*wf*(cae)*(1+grx*df(b,x))215 +es*(-df(b,x)*df(c,x)-df(b,y)*df(c,y)+(2-gamc)*df(c,z))216 +2*(1-gamc)*es*c/h217 );218 write length_sed:={length(sed),length(sedf),length(sedb)};219 ok:=if {sed,sedf,sedb}={0,0,0} then ok else 0;220 % update the concentration.221 gc:=gc+(dc:=-mean(3/2*(1-zz)*sed,zz)+3/4*(sedb-sedf/ct/r2*(1+2*cu)*rqq)/h);222 c:=c+csolv(sed+2*(1-zz)*dc,zz)*h/ct/r2*(1+2*cu)*rqq223 -sedf*(zz-1/2)/ct/r2*(1+2*cu)*rqq;224 showtime;225 if ok then write iter:=100000+iter;226 end;227 %write results228 r2:=sqrt(2)$ gz:=1$229 on rounded; print_precision 4;230 write dhdt:=gh; write dudt:=gu;231 write dvdt:=gv; write dcdt:=gc;232 end;99 *(((1+5*cu+8*cu^2)*vv(0,0)^2*rqq^2200 -(9+45*cu+80*cu^2+48*cu^3))*resttv201 +(1+5*cu+8*cu^2)*uu(0,0)*vv(0,0)*rqq^2*resttu);202 % solve the suspended sediment concentration.203 es:=ct*h^2*ros;204 sed:=df(c,t)+df(c*u,x)+df(c*v,y)+df(w*c,z)205 -df(wf*c,z)*the+grx*df(wf*c,x)206 -df(es*df(c,x),x)-df(es*df(c,y),y)-df(es*df(c,z),z);207 % boundary conditions on the free surface and bed.208 sedf:=sub(zz=1,209 (1+(1-gamc)*wf/6)*wf*c*(1+grx*df(h+b,x))210 +es*(-df(h+b,x)*df(c,x)211 -df(h+b,y)*df(c,y)212 +df(c,z)*gamc) +2*(1-gamc)*es*c/h);213 sedb:=sub(zz=0,214 (1+(1-gamc)*wf/6)*wf*(cae)*(1+grx*df(b,x))215 +es*(-df(b,x)*df(c,x)-df(b,y)*df(c,y)+(2-gamc)*df(c,z))216 +2*(1-gamc)*es*c/h217 );218 write length_sed:={length(sed),length(sedf),length(sedb)};219 ok:=if {sed,sedf,sedb}={0,0,0} then ok else 0;220 % update the concentration.221 gc:=gc+(dc:=-mean(3/2*(1-zz)*sed,zz)+3/4*(sedb-sedf/ct/r2*(1+2*cu)*rqq)/h);222 c:=c+csolv(sed+2*(1-zz)*dc,zz)*h/ct/r2*(1+2*cu)*rqq223 -sedf*(zz-1/2)/ct/r2*(1+2*cu)*rqq;224 showtime;225 if ok then write iter:=100000+iter;226 end;227 %write results228 r2:=sqrt(2)$ gz:=1$229 on rounded; print_precision 4;230 write dhdt:=gh; write dudt:=gu;231 write dvdt:=gv; write dcdt:=gc;232 end;