• Author / Uploaded
• चन्द्र प्रकाश

Citation preview

   

 

Two‐dimensional modelling of free surface shallow  water flow         

Hydraulic Reference Manual  Iber v1.0  19.07.2014   

   

Two-dimensional modelling of free surface shallow water flow

REFERENCE MANUAL

1. PRESENTATION ...................................................................................... 5  2. HYDRODYNAMIC MODULE ...................................................................... 7  2.1. Introduction ........................................................................................7  2.2. Hydrodynamic Equations .......................................................................7  2.3. Bed Friction .........................................................................................8  2.4. Superficial Friction due to Wind ..............................................................9  2.5. Effective Stresses.................................................................................9  2.6. Hydrodynamic boundary Conditions ...................................................... 11  2.6.1. Closed Boundaries ........................................................................ 11  2.6.2. Open Boundaries .......................................................................... 13  2.7. Internal boundary conditions ............................................................... 14  2.7.1. Gates .......................................................................................... 15  2.7.2. Weirs .......................................................................................... 15  2.7.3. Gate-Weir Combination.................................................................. 16  2.7.4. Local Losses ................................................................................. 16  2.8. Infiltration ........................................................................................ 17  2.8.1. Green-Ampt ................................................................................. 17  2.8.2. Horton ........................................................................................ 18  2.8.3. Linear ......................................................................................... 19  2.9. Initial Abstraction............................................................................... 19  3. TURBULENCE MODULE .......................................................................... 21  3.1. Introduction ...................................................................................... 21  3.2. Scales of Turbulence in Shallow Waters ................................................. 22  3.3. Constant Turbulent Viscosity................................................................ 23 

Iber. Turbulent Flow in Shallow Waters

3.4. Parabolic Profile of turbulent Viscosity ................................................... 23  3.5. Mixing Length Model ........................................................................... 24  3.6. The Rastogi and Rodi K-E model (Rastogi & Rodi, 1978) .......................... 24  3.7. Dimensional Analysis of the turbulent terms in the shallow water equations 25  4. NON-STATIONARY SEDIMENT TRANSPORT MODULE ............................ 28  4.1. Sediment Conservation Equations......................................................... 28  4.2. Bed Transport.................................................................................... 29  4.2.1. Stress Partition ............................................................................. 29  4.2.2. Bedload Discharge ........................................................................ 29  4.2.3. Bed slope correction ...................................................................... 31  4.2.4. Avalanche model .......................................................................... 32  4.2.5. Considerations on a non-erodible layer ............................................ 32  4.3. 2D Turbulent Suspended Load Module................................................... 32  4.3.1. Turbulent Suspended Load Equation ................................................ 32  4.3.2. Entrainment/Deposition (E-D) term Calculation ................................. 33  4.3.3. Sedimentation Velocity .................................................................. 35  5. NUMERICAL SCHEMES .......................................................................... 37  5.1. Calculation Mesh ................................................................................ 37  5.2. Discretization in finite volumes for the 2D-SWE Equations ........................ 38  5.2.1. Discretization of the convective flux terms ........................................ 39  5.2.2. Discretization of the bed slope source term ...................................... 42  5.3. Discretization of the Transport Equations in the K-E Turbulence Model, and in the Suspended load Transport Model ........................................................... 43  5.3.1. Depth-averaged Transport Equation ................................................ 43  5.4. Discretization of the Exner Conservation of Sediment Equation ................. 45  5.4.1. Non-erodible bed layer .................................................................. 46  5.5. Managing Wet-Dry Fronts .................................................................... 47  Page 2 of 56

Iber. Turbulent Flow in Shallow Waters

6. REFERENCES ......................................................................................... 50  7. NOMENCLATURE ................................................................................... 52       

 

Page 3 of 56

Iber. Turbulent Flow in Shallow Waters

         

  IBER  Turbulent Flow in Shallow Waters       

PRESENTATION 

Page 4 of 56

Iber. Turbulent Flow in Shallow Waters

1. PRESENTATION  Iber  is  a  numerical  model  for  simulating  turbulent  free  surface  unsteady  flow  and  environmental  processes  in  river  hydraulics.  The  ranges  of  application  of  Iber  cover  river  hydrodynamics,  dam‐break  simulation, flood zones evaluation, sediment transport calculation and wave flow in estuaries.  At this moment Iber has 3 main computational modules: a hydrodynamic module, a turbulence module  and a sediment transport module. All of them work in a finite volume non‐structured mesh made up of  triangle or quadrilateral elements. In the hydrodynamic module, which constitutes the base of Iber, the  depth averaged two‐dimensional shallow water equations are solved (2D Saint‐Venant Equations). The  turbulent  module  allows  including  turbulent  stressess  in  the  hydrodynamic  calculation,  in  this  way  allowing its use for different turbulent models for shallow waters at different degrees of complexity. The  most recent version of Iber also includes a parabolic model, a mixing length model and a k‐ε model. The  sediment transport module solves the bedload and the turbulent suspended load transport equations,  based on the sediment mass balance evolution at the bed.  This  manual  describes  in  detail  the  equations  and  models  used  in  the  different  Iber  computational  modules, including the used numerical schemes.    In order to cite the Iber model in technical documents and publications the following reference should  be used:    Bladé, E., Cea, L., Corestein, G., Escolano, E., Puertas, J., Vázquez‐Cendón, M.E., Dolz, J., Coll, A. (2014).  "Iber:  herramienta  de  simulación  numérica  del  flujo  en  ríos".  Revista  Internacional  de  Métodos  Numéricos para Cálculo y Diseño en Ingeniería, Vol.30(1) pp.1‐10    

 

Page 5 of 56

Iber. Turbulent Flow in Shallow Waters

         

  IBER  Turbulent Flow in Shallow Waters       

HYDRODYNAMIC MODULE 

Page 6 of 56

Iber. Turbulent Flow in Shallow Waters

2. HYDRODYNAMIC MODULE 

2.1. Introduction  The hydrodynamic module solves the depth averaged Shallow Water Equations (2D‐SWE), also known as  the two‐dimensional St. Venant Equations. These equations consider a hydrostatic pressure distribution  and  a  relatively  uniform  distribution  of  the  depth  velocity.  The  hydrostatic  pressure  hypothesis  reasonably complies in river flows and estuary wave currents. Likewise, the hypothesis of uniform depth  velocity distribution usually complies in rivers and estuaries, even though there may be zones in which  this requirement is not fulfilled due to local three‐dimensional flows or saline wedges. In these cases it is  necessary to study the extension of these zones and its possible repercussion in the model results. At  this moment, the numerical models based on the 2D‐SWE are widely used in coastal and river dynamics  studies, flood zones evaluation and sediment and contamination transport calculations. 

2.2. Hydrodynamic Equations  The  hydrodynamic  module  solves  the  conservation  of  mass  and  momentum  equations  in  the  two  horizontal directions. 

h hU x hU y    MS t x y e Zs τ s,x τ b, x g h 2 ρ hU x hU 2x hU x U y  hτ exx  hτ xy    gh     2 Ω sinλ U y    MX t x y x x y   ρ 2 x

hU y t



hU x U y x



hU 2y y

 gh

 hτ exy  hτ eyy Zs τ s,y τ b, y g h 2 ρ     2 Ω sinλ U x    MY y x y   ρ 2 y  

Where h is depth, Ux, Uy are the depth averaged horizontal velocities, g is the gravity acceleration, Z is  the free layer elevation, τs is the free surface friction due to wind induced friction, τb is the bed friction, ρ  is the water density, Ω is the earth’s rotation angular velocity, λ is the latitude of the studied point, τexx,  τexy, τeyy are the effective horizontal tangential stressess, and Ms, Mx, My are respectively the terms of:  mass source/drain and momentum, which are used to model precipitation, infiltration and drainage.  The following terms are included in the hydrodynamic equations:   Hydrostatic pressure   Bed Slope  Page 7 of 56

Iber. Turbulent Flow in Shallow Waters

 Viscous and Turbulent tangential stresses   Bed Friction   Superficial friction due to wind   Rainfall (Precipitation)   Infiltration  Likewise, wet/dry fronts, both stationary and non‐stationary, that may appear in the area of study are  modelled.  These  fronts  are  fundamental  in  the  fluvial  and  coastal  flood  zones  modelling  (rivers  and  estuaries). In this way the possibility to evaluate the extension of flood zones in rivers is introduced, as is  the movement of wave fronts in estuaries and coastal zones. 

2.3. Bed Friction  The bed wields a friction force over the fluid that is equivalent to the friction with a wall, even though, in  general the bed’s roughness has a greater value in the case of rivers and estuaries.  Bed  friction  has  a  double  effect  on  the  flow  equations.  On  one  side  it  produces  a  friction  force  that  opposes  the  mean  velocity,  and  on  the  other  side,  it  produces  turbulence.  Both  effects  can  be  characterized  by  the  friction  velocity  uf,  which  is  a  form  of  expressing  the  bed’s  tangential  stress  in  velocity units: 

uf 

τb   ρ

Where τb is the bed friction force module and ρ is the water density.  In  depth  averaged  models  it’s  not  possible  to  calculate  friction  velocity  by  means  of  standard  wall  functions,  as  it  is  done  in  wall  type  boundaries,  because  these  equations  are  solved  in  a  vertical  direction.  This  means  that  it’s  necessary  to  relate  friction  velocity  uf  with  depth  averaged  velocity  by  means of a friction coefficient. The bed stress can be expressed as: 

τ b  ρu f2  ρCf U

2

 

Where Cf is the bed friction coefficient. There are different expressions to approximate this coefficient.  Most of them assume uniform channel flow in a logarithmic profile of depth velocities. 

Page 8 of 56

Iber. Turbulent Flow in Shallow Waters

Unlike 1D models, in 2D models the hydraulic radius is not defined as the wet section area over the wet  perimeter, since in 2D models it makes no sense to define a transversal section. Taking a column of fluid  of length Δx and depth h, the hydraulic radius would be calculated as: 

Rh 

A h Δx  h  Pm Δx

Meaning that in 2D models hydraulic radius and depth are equivalent.  The  bed  friction  is  evaluated  using  the  Manning  formula,  which  uses  the  Manning  n  coefficient  as  parameter. The Manning formula uses the following roughness coefficient: 

n2 C f  g 1/3   h

2.4. Superficial Friction due to Wind  The friction force due to wind over the free surface can be calculated considering the wind velocity at 10  meters of height and a dragging coefficient, using the Van Dorn equation (Van Dorn, 1953). 

τs  ρ Cvd V10 2   Where ρ is the water density, V10 is the wind velocity at 10 meters of height and Cvd is the surface drag  coefficient. By default a drag coefficient of Cvd =2.5 10‐6 is considered.  

2.5. Effective Stresses  The horizontal effective stresses that appear in the hydrodynamic equations include the effects of the  viscous stresses, of the turbulent stresses and the dispersion terms due to the non‐homogeneity of the  depth velocity profile. 

τ ije  τ ijv  u'i u' j  D ij   v

Where  τij  are the viscous stresses,  u' i u' j  are the turbulent stresses (also called the Reynolds stresses),  and Dij are the lateral dispersion terms: 

D ij 

1 Zs U i  u i  U j  u j dz   h Z b

Page 9 of 56

Iber. Turbulent Flow in Shallow Waters

The dispersion terms are not considered in the 2D‐SWE equations (remember the hypothesis of uniform  depth velocity profile), because it is impossible to calculate them in a general way with a depth averaged  model. Their importance will be greater as less uniform is the depth velocity profile. A typical situation  in  which  these  terms  can  gain  importance  is  in  channels  with  elbows  or  small  curvature  radii,  like  in  confluence channels (Figure 1) 

Q1

Q

Q

3

 

2

 

Figure 1. Secondary flows (left) and depth velocity profile (right). Main causes of the dispersion terms.  The viscous stresses are calculated from the fluid’s kinematic viscosity,  ν  as: 

 U U j    τ ijv  ν i   x   x i   j In general, the order of magnitude of the viscous stresses is much lower than the rest of the terms that  appear in the hydrodynamic equations, except near the walls and in laminar flow.  The turbulent stresses are much larger in magnitude than the viscous stresses, especially in recirculation  zones, where turbulence production is elevated. In the case of two‐dimensional shallow equations the  turbulent stresses constitute three new unknowns to calculate, added to the depth and velocities Ux, Uy  they  mean  a  total  of  6  unknowns.  This  is  what  is  known  as  turbulence  closure  problem,  because  it  is  necessary to solve 3 equations with 6 unknowns. Due to this, it is necessary to use a turbulence model  that  allows  calculating  these  turbulent  stresses.  Most  of  the  models  calculate  the  turbulent  diffusive  terms from the following expression: 



 u'i u' j x j



  U i   νt   x j  x j 

Where  ν t  is the turbulent viscosity calculated in the turbulence module. The problem lies in the non‐ existence  of  a  universal  turbulence  model,  that  may  allow  calculating  in  a  precise  form  the  turbulent  stresses.  Through  time  some  models  have  been  developed  with  different  degrees  of  complexity.  The  Boussinesq formulation is used by all turbulence modules in Iber.  Page 10 of 56

Iber. Turbulent Flow in Shallow Waters

2.6. Hydrodynamic boundary Conditions  In a two‐dimensional problem it is necessary to distinguish two types of boundaries: open and closed.  The  closed  boundaries,  also  called  wall‐type  boundaries,  are  water‐resistant,  meaning  that  they  don’t  allow water through them.  

2.6.1. Closed Boundaries  The presence of a wall‐type boundary generates a lateral friction force in the fluid in a similar fashion  than the friction exerted by the bed’s roughness. At a closed boundary the following conditions can be  used:   Slip condition (null tangential stress)   Wall friction Condition (wall functions)  The slip condition is equivalent to undermining the friction stress generated by the wall‐type boundaries  over  the  fluid.  In  hydraulic  engineering,  and  especially  in  river  engineering,  the  surface  contact  with  lateral boundaries is much lower than the surface contact with the bed, due to the separation between  horizontal  and  vertical  scales,  so  the  friction  forces  in  the  wall  boundaries  can  be  undermined.  In  this  case the slip condition would hold.  In  problems  where  the  horizontal  and  vertical  dimension  are  similar  (narrow  section  channels)  this  roughness  force  can  have  some  importance  in  the  flow  development,  even  though  in  general  the  influence is small. If you want to consider the lateral friction, a wall‐type friction can be introduced; this  consists  in  imposing  a  tangential  force  in  the  opposite  direction  to  the  flow.  In  this  case,  Iber  distinguishes  between  smooth  turbulent  regime  and  rough  turbulent  regime  in  function  of  the  wall’s  roughness and the flow velocity close to the wall.  The wall friction velocity (u*) is defined in function of the wall friction ( τ w ) as: 

u* 

τw   ρ

The tangential velocity at the wall can be expressed as a function of the friction velocity, the roughness’s  depth and the wall distance as: 

u 

u* Ln  E  y      κ

y 

y u*   ν

Page 11 of 56

Iber. Turbulent Flow in Shallow Waters

Where y is the perpendicular distance to the wall, and E is a parameter which value depends on the flow  characteristics.  To  calculate  E,  Iber  considers  smooth  turbulent  flow  conditions,  rough  turbulent  flow  conditions and a transition between smooth and rough (Table 1). 

KS 

 Regimen type 

Smooth Turbulent  

K S  5  

 

Rough Turbulent 

KSu*   ν

u 

u* Ln  E  y     κ

E  9.0  

5 < KS  70  

E=

30   KS

KS  70  

E=

1   0.11 + 0.033  K S

Smooth‐Rough Transition  

Table 1. Wall Friction.  The turbulence is defined as smooth when this relationship complies: 

KS 

KSu*  5  ν

Where  Ks  is  the  wall’s  roughness  height,  a  measure  of  the  wall’s  roughness,  with  length  units.  In  this  condition the tangential velocity at the wall can be expressed as a function of friction velocity and the  kinematic viscosity, as: 

u 

u*  yu  Ln  9.0 *      κ 

Rough turbulent regime is applied when this relationship complies: 

K S 

K Su*  70   ν

In this condition, the tangential velocity at the wall can be expressed as a function of the friction velocity  and the bed’s roughness’s height, as 

u 

 u* y  Ln  30   κ  KS 

In the transition between regimes, the tangential velocity at the wall can be expressed as a function of  friction’s velocity, kinematic viscosity and roughness’s height, as: 

Page 12 of 56

Iber. Turbulent Flow in Shallow Waters

    u y   u  * Ln  κ  0.11 ν + 0.033  K  S   u*   2.6.2. Open Boundaries  In  open  bondaries,  different  types  of  boundary  conditions  can  be  introduced.  So  that  the  2D‐SWE  equations are correctly presented from a mathematical point of view, the number of conditions to use  in the open boundaries will depend if the boundary is an inlet or an outlet, it will also depend on the  type of regime at the boundary (fast/slow). In an inlet boundary it is necessary to introduce 3 boundary  conditions if the regime is supercritical (one for each of the 3 St. Venant equations), while if the regime  is  subcritical,  it  is  enough  to  present  2  conditions.  In  an  outlet  boundary,  it  is  enough  to  introduce  a  single condition if the regime is subcritical; meanwhile, it is not necessary to establish any condition if  the  regime  is  supercritical.  If  the  user  introduces  fewer  conditions  than  those  necessary,  from  a  mathematical point of view the equations will be undetermined and no correct solution will be found.  The conditions which can be entered by the user are height, velocity components, or a combination of  both. In Iber different options are considered to introduce boundary conditions, as explained in Table 2.  It is common in river hydraulics for the flow to be in subcritical regime in the contours of the modeled  section.  In  this  case,  what  is  often  done  is  to  introduce  free  surface  depth  or  water  level  in  the  downstream boundary. In the upstream boundary the most common is to enter a total inflow (in cubic  meters  per  second  cms)  in  the  direction  of  the  flow,  which  in  general,  lacking  more  precise  data,  is  assumed to be perpendicular to the inlet boundary. It is also possible to introduce upstream velocities  (in m/s) or specific flow (m2/s), though this occurs less often. In the case of total flow in the inlet, a unit  flow distribution (m2/s) is entered, according to the following expression: 

qn 

h 5/3 5/3  h dy

Q 

Where qn is the normal specific discharge (m2/s) at each point of the inlet, and Q is the total inlet flow at  mentioned boundary. The integral in the denominator is extended along the length of the boundary.  Besides  height,  other  possibilities  can  be  considered  at  the  outlet,  such  as  weirs  or  rating  curves.  The  weir condition establishes a relationship between outflow and height at each point of the contour: 

q = Cd (ZS  ZW )1.5   Page 13 of 56

Iber. Turbulent Flow in Shallow Waters

Where Cd is  the weir’s discharge  coefficient, Zs is the water level upstream  of the weir and Zw is the  water level on top of the weir. The user must introduce Cd and Zw.  The  rating  curve  contour  condition  establishes  a  general  relationship  between  the  outflow  and  the  water  level  en  each  contour  point.  This  relationship  is  introduced  by  the  user  using  a  table  defining  values of specific discharge and water level.  The set of conditions implemented in Iber in open contours are shown in Table 2.    Boundary 

 

Regime 

Boundary conditions 

Subcritical / Critical 

Total  discharge  in  normal  direction  to  the  boundary 

Supercritical 

Total  discharge  in  normal  direction  to  the  boundary and average velocity  

Subcritical / Critical 

Specific  discharge  in  normal  direction  to  the boundary 

Discharge 

Inlet   

a) Specific discharge in normal direction to  the boundary and water depth 

Unit discharge  Supercritical 

b) Specific discharge in normal direction to  the  boundary  and  water  surface  elevation  a) Water depth  b) Water surface elevation 

Subcritical  c) Weir (elevation and discharge coefficient) 

Outlet 

d) Rating curve  Supercritical / Critical 

Not necessary to introduce any condition 

Table 2. Contour Conditions implemented in open contours. 

2.7. Internal boundary conditions  The internal conditions are used to model hydraulic structures like gates, weirs or bridges which change  the conditions of the system. 

Page 14 of 56

Iber. Turbulent Flow in Shallow Waters

The  internal  boundary  conditions  implemented  in  Iber  can  be  used  to  model  the  following  flow  conditions:   Flow below a gate   Flow over a free falling weir   Weir‐Gate combination   Local loss 

2.7.1. Gates  The equation used is the gate drainage, which works on free or submerged flow. The data to insert are  the discharge coefficient, the gate’s bed level, the gate opening and its width. By default, the discharge  coefficient will have a value of Cd=0.6. 

ZU ZD

h ZB Figure 2.  Scheme and equations of the gate internal contour condition. 

(ZD  ZB ) / (ZU  ZB )  

Discharge Equations 

Free Flow 

0.00 – 0.67 

Q = Cd B h 2g (Z U  ZB )  

Transition 

0.67 – 0.80 

Q = Cd B h 6g (ZU  ZD )  

Submerged Flow 

0.80 – 1.00 

Q = Cd B h 2g (ZU  ZD )  

 

 

2.7.2. Weirs  The rectangular weir discharge equation is used, which can work for free flow and submerged flow. The  data to input are the weir’s upper level, the discharge coefficient and the weir’s length. By default the  discharge coefficient has a value of Cd=1.7. 

Page 15 of 56

Iber. Turbulent Flow in Shallow Waters

ZU Zw

ZD

ZB Figure 3.  Scheme and equations of the weir internal contour condition. 

(ZD  ZW ) / (ZU  ZW )  

Discharge Equation 

Free Flow 

 0.67 

Q = 2.598 Cd B (ZD  ZW ) (ZU  ZW )0.5

 

 

2.7.3. Gate‐Weir Combination  In this case the condition includes the two previous ones, so both parameters must be used, the gates  and weirs. The total discharged flow is taken as the sum below the gate and the flow above the weir. 

ZU Zw

ZD

Q  Qcompuerta  Qvertedero   h ZB Figure 4. Scheme and equations of the gate‐weir internal contour condition. 

2.7.4. Local Losses  In  this  case,  there  is  a  local  energy  loss  in  the  flow  transference  between  two  finite  volumes  of  value  H=

v2/2g.  The  St.  Venant  Equations  are  the  mathematical  expressions  of  the  conservation  of  mass 

and momentum, a way to consider this loss of energy is over the friction slope. For this, an additional  term is added to the friction slope, Sf, through the edge of the finite volume, this added term is equal to  H/V, where V is the element’s volume. In this way, the energy loss through the two elements will be 

Page 16 of 56

Iber. Turbulent Flow in Shallow Waters

equal  to  H+  Sf∙L,  where  L  is  the  distance  between  element’s  centre  on  both  sides  of  the  boundary  where the local loss is implemented. 

 

v2 H     2g H S 'f  S f  V

 

Figure 5. Scheme and equations of the local loss internal boundary condition. 

2.8. Infiltration  To simulate rainfall processes, it can be necessary to consider water infiltration in non‐saturated terrain  for  calculating  runoff.  Modeling  infiltration  of  surface  waters  is  especially  important  in  rainfall‐runoff  models.  Infiltration is considered in the model using a negative term in the conservation of mass equation (loss  of water): 

h hU x hU y   =-i  t x y Where i is the real infiltration rate, measured as the minimum between the potential infiltration rate f  (soils infiltration capacity at every instant, which depends on soil characteristics and conditions), and the  amount of readily available surface water for infiltration. 

i = min ( f ,

h )  Δt

To  calculate  the  potential  infiltration  3  infiltration  common  models  are  used:  The  Green‐Ampt  model,  the Horton model and linear model. 

2.8.1. Green‐Ampt  The infiltration rate, expressed in m/s, is calculated in each cell using the Green‐Ampt formula (Chow,  Maidment,  &  Mays,  1988),  in  it;  it  is  assumed  that  there  exists  a  saturated  front  that  separates  the  saturated  soil  region,  immediately  below  the  terrain,  and  the  non‐saturated  region,  where  there  is  suction. 

Page 17 of 56

Iber. Turbulent Flow in Shallow Waters

As the water infiltrates, the saturated front descends and the depth of the saturated region L becomes  bigger. The infiltration rate f is calculated as: 

  h  Ψ  Δθ  f  k s 1   L0  Δθ  F  

t

F   f dt 0

L  L0 +

F Δθ

Δθ    θi  

Where ks the saturated permeability of the soil, h is the depth, ψ is the suction in the non‐saturated soil  region, Δθ the humidity content change of the soil as the saturated front increases, θi the soil’s initial  humidity content,   the soil’s total porosity, and L is the depth of the saturated region. The infiltration  rate is equal to the potential infiltration rate as long as there is enough surface water to infiltrate.  The parameters to introduce by the user for this module are:  

The soil’s saturated permeability (ks) 



The non‐saturated region’s suction (ψ) 



Soil’s effective porosity (drainable) (θe) 



Soil’s initial effective saturation (Se), defined as: 

Se =

θi - θ r θe

 

Where  θ r  is the soil’s retention capacity (non‐drainable humidity) and  θi the soil’s initial humidity. The  soil’s porosity    is equal to the drainable porosity plus the soil’s retention capacity ( 

= θe  θr

). 

Starting  with  the  effective  porosity  and  the  soil’s  initial  effective  saturation,  the  change  in  humidity  content can be calculated as the saturation front advances with the formula: 

Δθ =  - θi =  - θr - θe  Se  θe  1  Se    All the parameters of the Green‐Ampt equation can be introduced as spatial variables (different for each  element of the mesh). 

2.8.2. Horton  In the Horton model, the potential infiltration rate is calculated as: 

f  f c   f 0 - f c   exp   k  t   

Page 18 of 56

Iber. Turbulent Flow in Shallow Waters

Where  t  is  time,  starting  with  the  beginning  of  the  precipitation.  The  user  must  introduce  as  model  parameters  an  initial  infiltration  rate  (݂௢ ),  the  infiltration  rate  at  infinite  time  (݂௖ )  and  the  constant  k,  which defines the temporal variation of the potential infiltration rate. All the parameters of the Horton  infiltration equation are spatially variable (calculated for each element of the mesh). 

2.8.3. Linear  The  linear  model  considers  an  initial  abstraction  P0  (volume  per  unit  of  area),  and  afterwards  some  continuous constant losses (volume per unit of area and unit of time). The value of the initial abstraction  and the continuous losses can vary per element. 

  Figure 6. Temporal evolution of the infiltration rate according to the linear model 

2.9. Initial Abstraction  If the Green‐Ampt or Linear infiltration models are used to calculate the infiltration losses, there is also  the  possibility  to  consider  an  initial  abstraction.  This  initial  abstraction  can  represent  other  processes  like the canopy retention (because of vegetation) or land depressions, or the infiltration capacity of dry  soils with an elevated porosity.  The initial abstraction is defined as a volume per unit of area, and as that has length unit. This value is  subtracted from the water that reaches the soil, being from rainfall or from surface runoff. That way it  can be used either in areas with or without rainfall. 

 

 

Page 19 of 56

Iber. Turbulent Flow in Shallow Waters

          IBER  Turbulent Flow in Shallow Waters       

TURBULENCE MODULE 

Page 20 of 56

Iber. Turbulent Flow in Shallow Waters

3. TURBULENCE MODULE 

3.1. Introduction  A great number of studies in hydraulic engineering implies the analysis of open channel flows, many of  which can be considered shallow flows, understanding this as low depth with small vertical‐horizontal  relationship. Practically all open channel flows are turbulent. While looking at a river, small swirls appear  and  disappear  with  an  almost  chaotic  movement,  showing  the  complexity  of  the  turbulent  motion.  These swirls are the main cause of mixing processes, so they play an important role in the diffusion of  soluble substances and suspended solids, among others.  Despite  the  fact  that  all  flows  in  river  engineering  are  turbulent,  in  many  cases  the  turbulence  is  not  large enough to have a notorious influence in the mean velocity. Such is the case of rivers, estuaries and  in  general  in  coastal  zones  with  a  smooth  enough  geometry  so  that  no  recirculation  zones  appear.  Nevertheless, even in this type of situations it’s important to consider turbulence while modelling these  water systems, since it plays an important role in the processes of transport and mixing of contaminants  and sediments. The diffusion of heat, of a solute, or suspended sediment is produced basically because  of turbulence, except in laminar flow, which generally does not occur in hydraulic engineering, even less  in rivers and estuaries. The turbulent diffusion coefficient is several times higher in magnitude than the  molecular  diffusion  coefficient.  Because  of  this,  it  is  necessary  to  evaluate  previously  the  turbulent  kinetic energy to calculate the diffusive flow.  One of the main characteristics of Iber is the inclusion of diverse RANS type turbulence models, which  are solved in the turbulence module. The following turbulence models are included for shallow waters,  enumerated by order of complexity:   Constant Turbulent Viscosity   Parabolic Model   Mixing Length Model   k‐ε (ka‐epsilon) model from Rastogi ‐ Rodi (Rastogi y Rodi, 1978)  Using models with different levels of complexity allows selecting the most adequate for each case study,  considering  the  complexity  of  the  flow  and  the  model.  In  general  the  mixing  length  model  provides  satisfactory results in rivers and estuaries, this way discouraging the use of other turbulence models in  these  cases.  In  hydraulic  structures  such  as  open  flow  channels  with  marked  elbows  and  recirculation 

Page 21 of 56

Iber. Turbulent Flow in Shallow Waters

zones, it may be necessary to use at least one mixing length model and one k‐ε model to test the results.  The  choice  of  the  turbulence  model  that  best  fits  each  case  is  based  on  user  experience,  considering  that  the  more  complex  the  model  is,  the  longer  it  takes  to  calculate  and  more  complex  the  equation  solution is.  The  objective  of  the  turbulence  models  is  to  calculate  the  Reynolds  stresses.  Models  based  on  the  Boussinesq hypothesis (the ones enumerated before, used in Iber), the Reynolds stresses are evaluated  from the following expression: 

 U i U j  2   k δ ij     uiu j  ν t   x   j x i  3 The turbulence model calculates the turbulent viscosity to be used in the previous expression. 

3.2. Scales of Turbulence in Shallow Waters  One of the main characteristics of shallow flows is the split‐up between horizontal and vertical scales,  this  happens  because  the  vertical  extension  of  the  fluid  (limited  by  depth)  is  much  smaller  than  the  horizontal extension. This scale split‐up is applicable in the spatial dimension, velocities, and of course  the  turbulence.  Speaking  of  turbulence,  its  main  effect  is  the  fragmentation  between  tri‐dimensional  (swirls) and bi‐dimensional turbulent structures. The 3D turbulence spatial scale is limited by its depth,  and  thus  its  structures  are  much  smaller  than  those  associated  to  the  2D  turbulence,  which  are  only  limited  by  the  horizontal  scale.  The  3D  turbulence  is  originated  mainly  by  the  bed  friction,  while  the  origins of the 2D turbulence are the horizontal velocity gradients   It’s important for the turbulence model to include not only the effects of the 3D turbulence, formed by  the bed friction, but also the 2D turbulence, formed by the horizontal velocity gradients. In the shallow  water  models,  the  bi‐dimensional  characteristics  of  the  flow  are  implicitly  considered  in  the  transport  equations  when  considering  a  homogenous  depth  velocity  profile,  while  the  tri‐dimensional  causes  of  turbulence are usually included using another term which depends on the bed tangential stress. In the  same  way,  even  when  considering  3D‐SWE,  the  turbulence  model  should  consider  turbulence  anisotropy in the horizontal and vertical directions.  Iber  uses  the  following  turbulence  models.  All  of  them  are  shallow  water  depth  averaged  turbulence  models. 

Page 22 of 56

Iber. Turbulent Flow in Shallow Waters

3.3. Constant Turbulent Viscosity  The order of magnitude of the turbulent viscosity can approximated depending on the considered flux.  Several  published  articles  propose  these  turbulent  viscosity  approximated  values  in  function  of  the  considered flux. This focus is very simple, and cannot be considered as an adequate, nor realistic in any  case, turbulence model, since it does not consider spatial changes of turbulent viscosity. It’s important  to  highlight  the  importance  of  spatially  variable  turbulence  viscosity  to  the  mean  velocities.  Furthermore, the existing tables only provide approximated values. Because of these reasons, the use of  this method is discouraged, since it can lead to excessively large errors, generally coming from the use  of large values of turbulent viscosity and from not considering its spatial variability.  Another  important  downfall  of  the  method  is  produced  while  modeling  non‐stationary  flow  regimes,  since in these cases; turbulence varies through space and time. 

3.4. Parabolic Profile of turbulent Viscosity  This model assumes a parabolic distribution of turbulent viscosity in depth, this way considering a mean  viscosity in depth which can be calculated from the following expression: 

ν t  0.068 u f h   Where h is the total depth and uf is the bed friction velocity, calculated from the bed tangential stress,  from the expression: 

uf 

τf   ρ

If the Manning formula is used to calculate the bed friction, the following formula for turbulent viscosity  can be used: 

ν t  0.068 g n U h 5/6   This means that the turbulent viscosity depends on the local depth, the depth averaged velocity and the  manning  coefficient.  For  better  results,  a  coefficient  is  introduced  to  adjust  the  turbulent  viscosity,  which can be handled by the user. 

ν t  C m 0.068 g n U h 5/6  

Page 23 of 56

Iber. Turbulent Flow in Shallow Waters

3.5. Mixing Length Model  In  the  shallow  waters  mixing  length  model,  the  turbulent  viscosity  is  calculated  based  on  the  local  characteristics of the flow using the following expression:  2

ν t  min 0.267 κ h, κ d wall 

2

 u  2S ij S ij   2.34 f    κh 

Where  κ=0.41  is  the  von  Karman  constant.  It’s  a  relatively  simple  algebraic  model  that  produces  acceptable  results  in  flows  where  the  turbulence  is  locally  generated  mainly  from  bed  friction.  It  considers  the  turbulence  caused  by  horizontal  velocity  gradients,  but  it  does  not  consider  convective  transport or turbulence dissipation. Using the model in flows with strong recirculation zones, the results  have large errors. 

3.6. The Rastogi and Rodi K‐E model (Rastogi & Rodi, 1978)  In  this  model  k  represents  the  turbulent  kinetic  energy  and  ε  the  rate  of  dissipation  of  the  turbulent  energy.  The  model  considers  turbulence  caused  by  bed  friction,  velocity  gradients,  dissipation  and  convective transport. It solves the following shallow water equations for k‐ε model: 

k U x k U y k     x y x j t  ε U x ε U y ε    y x j t x

ν t  cμ

k2 ε

    ν  ν t  k  σ k  x j 

 ν  ν  t  σε 

c k  c f1/2

 ε   x j

3    2ν t S ij S ij  c k u f  ε    h 

4 2    c ε1 ε 2ν t S ij S ij  c ε u f  c ε2 ε    k k h2 

1/2 c ε  3.6c 3/2 k c ε2 c μ

cf 

τb 1   ρ U2

With the following constants: 

c μ  0.09

c ε1  1.44

c ε2  1.92

σ k  1.0

σ ε  1.31  

Where  k  is  the  kinetic  turbulent  energy,  ε  is  the  turbulence  dissipation  rate,  Sij  is  the  deformation  sensor. The terms which include the bed friction velocity uf  are responsible for modeling the forming of  turbulence due to bed friction.  The  k‐ε  model  is  relatively  sophisticated.  The  results  for  shallow  water  turbulent  flows  are  relatively  good, being one of the most used models, whenever turbulence levels are important. Nevertheless, its  Page 24 of 56

Iber. Turbulent Flow in Shallow Waters

level of complexity does not guarantee correct results for all types of flows. Like any other turbulence  model, the results from the k‐ε model must be critically analyzed, for this, user experience in turbulent  flow is essential. 

3.7. Dimensional Analysis of the turbulent terms in the shallow water equations  If  the  shallow  water  equations  are  turned  dimensionless,  the  following  non‐dimensional  numbers  are  obtained: 

F

U gh

T

1 H Cf L

Rl 

UL ν

Rt 

UL   νt

Which  respectively  refer  to  the  relationship  between  water  mass  inertia  (convective  forces)  and  pressure  forces  (F),  bed  friction  force  (T),  laminar  tangential  stresses  (Rl)  and  turbulent  tangential  stresses  (Rt).  The  relative  importance  of  these  associated  processes  to  each  dimensionless  number  is  inversely proportional to the value of each number. For example, the larger a dimensionless number is,  the less important it will be for the process it represents. That way, for a large laminar Reynolds value,  the flow is turbulent and thus the laminar forces loose importance in the flow development. Likewise,  the  importance  of  turbulent  stresses  in  the  mean  velocity  will  depend  on  the  value  of  the  turbulent  Reynolds number, which depends on the turbulent viscosity. A parabolic model can be used to estimate  the order of magnitude of the turbulent viscosity, with the following expression: 

νt 

1 1 κ u f h  κ g n h 5/6 U  0.21 n h 5/6 U   6 6

The Manning formula has been used to estimate the bed friction velocity uf. This estimation will be more  precise in cases where the turbulence is generated mainly by bed friction, as is the case in rivers, and  will  provide  some  errors  in  cases  where  the  turbulence  is  originated  from  horizontal  shear  stress,  for  example  in  cases  of  recirculation  zones.  In  any  case,  using  this  approximation,  the  turbulent  Reynolds  number can be expressed as: 

Rt 

UL 4.8 L    νt n h 5/6

This expression can be used as a first evaluation of the importance of turbulent stress in the velocity and  depth.  For  example,  if  we  model  a  river  section  with  a  depth  of  10m,  a  width  of  400m,  a  manning  coefficient of 0.025 and a mean velocity of 0.5m/s, the approximate turbulent viscosity is 0.02m2/s and 

Page 25 of 56

Iber. Turbulent Flow in Shallow Waters

the  turbulent  Reynolds  number  is  close  to  Rt  ~  11000,  which  is  a  pretty  large  value,  so  the  turbulent  stresses are expected to have a small effect on the mean flow.  According  to  the  expression  for  Tb,  the  importance  of  bed  friction  grows  in  shallow  water  flows,  and  looses importance as the relationship between depth and the horizontal dimension grows.   

 

Page 26 of 56

Iber. Turbulent Flow in Shallow Waters

           

  IBER  Turbulent Flow in Shallow Waters       

SEDIMENT TRANSPORT MODULE 

Page 27 of 56

Iber. Turbulent Flow in Shallow Waters

4. NON‐STATIONARY SEDIMENT TRANSPORT MODULE  The sediment transport module solves the non‐cohesive sediment non‐stationary transport equations.  The  equations  include  the  bedload  transport  equations  and  the  suspended  sediment  transport  equations,  coupling  the  bedload  and  the  suspended  load  through  a  sedimentation‐rise  term.  The  sediment transport module uses the velocity, depth and turbulence fields from the hydrodynamic and  turbulence modules. The bedload is calculated using an empirical formula, chosen from the Meyer‐Peter  Muller  or  the  Van  Rijn  Methods.  The  suspended  load  transport  is  modeled  from  a  depth  averaged  turbulent transport equation.  Hydrodynamic + Turbulence

Bed-Load

Suspended Load

Sediment Conservation

Bed variations

 

Figure 7. Non‐Stationary sediment transport scheme. 

4.1. Sediment Conservation Equations  The bed level variation is calculated with the Exner sediment conservation equation: 

1  p 

Z b q sb,x q sb,y   D -E  t x y

Where p is the porosity of the sediment that forms the bed layer, Zb is the bed level, qsb,x and qsb,y are the  two solid flow components. The difference D‐E is the balance between the bedload and suspended load. 

Page 28 of 56

Iber. Turbulent Flow in Shallow Waters

4.2. Bed Transport  4.2.1. Stress Partition  The total bedload stress in the river bed is generated by the grain’s roughness (which is proportional to  the sediment diameter) and the shape of the bed layer (ripples, dunes, antidunes). Only the grain stress  contributes  the  bedload  movement.  Thus,  before  calculating  the  bed  sediment  flow  it  is  necessary  to  estimate this grain stress. To do this, Iber uses the Einstein formula based on the total stress:  1.5

n  τ = τ  s   n  * bs

* b

ns 

K1/6 s (m)

K s  2  3 D s   

25

Where n is the total Manning coefficient, ns is the equivalent Manning coefficient due to the grain, Ds is  the sediment’s diameter, Ks the grain’s roughness’s height (measured from the grain’s diameter),  τ b is  * * the  bed  total  stress,  τ bs is  the  bed  stress  due  to  the  grain,  τ b , τ bs are  dimensionless  total  and  grain 

stress, calculated with: 

τ*b =

τb  ρs  ρ  g Ds

τ*bs =

τ bs  ρs  ρ  g Ds

 

Where  ρs is the sediment’s density, and  ρ is the water density. Iber uses a value of  Ks  2.5 Ds  

4.2.2. Bedload Discharge  The Bedload discharge is calculated based on empirical formulas. In the latest version of the model, two  well regarded and well known formulas are used:   Meyer‐Peter Müller   Van Rijn   Meyer‐Peter Müller (Meyer‐Peter & Muller, 1948)  The  original  Meyer‐Peter  &  Müller  equation,  deduced  from  gravel  beds  of  diameters  up  to  30mm,  calculates the Bedload discharge with the following expression: 

q*sb = 8   τ*bs - τ*c  = 8   τ*bs - τ*c  3/2

3/2

 

Where the dimensionless discharge is calculated as:  Page 29 of 56

Iber. Turbulent Flow in Shallow Waters

q*sb 

q sb

 

 ρs  3  ρ  1 g Ds  

In  the  case  of  flat  beds,  a  dimensionless  critical  bed  stress  is  considered  of   c



 0.047 .  On  the 

contrary, it is necessary to make a bed slope correction. Detailed in section 4.2.3  Other corrections have been made to this formula, (Wong, M, 2003), (Wong & Parker, 2006), resulting  in: 

q*sb = 3.97   τ*bs - τ*c 

3/2

 

If  the  bed  is  flat,   c  0.0495 ,  is  considered.  On  the  contrary,  it  is  necessary  to  make  a  bed  slope  

adjustment, detailed in section 4.2.3. This corrected version is included in Iber.   Van Rijn (van Rijn, L.C., 1984)  In the van Rijn Formula the bedload discharge is calculated from the following expressions: 

T 2.1 T  0.3  q  0.053  0.3 D* * sb

T1.5 T  0.3  q  0.100  0.3 D*

 

* sb

Where  T  is  a  dimensionless  parameter  that  measures  the  excess  of  bed  friction  above  a  critical  value  which defines the movement threshold. 

τ*bs - τ*c T τ*c

 

The dimensionless diameter is defined as:  1/3

gR D*  Ds   2   ν 

con R 

γs -γ   γ

Page 30 of 56

Iber. Turbulent Flow in Shallow Waters

4.2.3. Bed slope correction  Whenever the bed is not flat, the equations must be corrected to include the effects of gravity, either to  increase  the  bedload  discharge  in  a  positive  slope,  or  to  reduce  it  in  an  adverse  slope.  The  correcting  factor  is  applied  over  the  critical  stress  of  movement  threshold,  as  detailed  by  Apsley  and  Stansby  (Apsley & Stansby, 2008), where they present and generalize the works of many other authors, citing at  least (Dey, 2003) and (Wu, 2004).  To consider the bed slope both in the movement threshold and in the bedload discharge, the sediment’s  weight,  due  to  the  bed’s  slope,  is  combined  in  a  vector  form  with  the  bed  stress  to  get  an  effective  stress. Si b is a unit vector in the direction of the maximum slope; the dimensionless effective stress is  defined as: 

τ*bs,eff = τ*bs + Do  sinβ  b   Where     is  the  maximum  slope,  with  the  horizon,  and  D0  is  a  parameter  depending  on  the  particle  shape. So that the movement begins in the absence of flow, whenever    is equal to the internal friction  angle of the material (  ), the parameter D0 is defined as: 

Do 

τ*c,0 tan 

 

*

Where  τ c,0  is the dimensionless critical stress for the flat bed. On the other side, the effective critical  stress is reduced proportionally to the gravity component normal to the bed slope. 

τ*eff,crit = τ*c,0  cosβ   Where 

τ*c,0   is  the  same  as  before.  Starting  at  this  point,  the  bedload  discharge  equations  presented 

above  are  used,  substituting  the  (critical  and  bed)  stresses  by  effective  stresses,  and  getting  the  sediment  discharge,  which  is  a  function  of  the  fluid’s  stress  and  the  bed  slope  in  each  of  the  x,  y  directions.  The formulas used are entirely in vector form for the bedload discharge, and thus able to consider any  flow direction with respect to the maximum slope 

Page 31 of 56

Iber. Turbulent Flow in Shallow Waters

4.2.4. Avalanche model  Apsley  and  Stansby  (Apsley  &  Stansby,  2008)  also  proposed  including  an  avalanche  sliding  model  to  avoid those slopes higher than the bed material’s friction angle. For that, if an angle    between two  finite  volumes  is  higher  than  the  friction  angle  ,  then  a  unit  bedload  discharge  appears  from  the  upstream to the downstream element equal to: 

0.5  L2  (tan β  tan  ) q aval =(1  p)    cos β  t Where L is the maximum horizontal dimension of the adjacent finite volumes. 

4.2.5. Considerations on a non‐erodible layer  While calculating bedload dragging and changes in the bed layer section, the possibility to include a rock  layer, or non‐erodible surface, below which the bed erosion cannot go further. 

4.3. 2D Turbulent Suspended Load Module  4.3.1. Turbulent Suspended Load Equation  The suspended load transport equation uses the velocity, depth and turbulence fields calculated in the  hydrodynamic  and  turbulence  modules.  The  suspended  load  is  modelled  using  a  depth‐averaged  equation. The equation used in the Iber is: 

ν hC hU x C hU y C          t t x y x j   Sc,t

 C  Dsx Dsy   E  D     h  y  x j  x

Where  C  is  the  depth‐averaged  concentration  of  suspended  solids,  Ux,  Uy  are  the  horizontal  depth‐ averaged velocity components,  ν t  is the turbulent viscosity, Г is the molecular diffusion coefficient for  suspended  solids,  and  Sc,t  is  the  Schmidt  number,  which  relates  the  moment  turbulent  diffusion  coefficient with the suspended turbulent diffusion coefficient.  The  terms  Dsx,  Dsy  model  the  suspended  sediment  dispersion  due  to  the  non‐homogeneous  vertical  velocity  profile  and  sediment  concentration.  Normally  its  effect  is  undermined  in  2D  shallow  water  models, despite its importance may be relevant when the concentrations and velocities vary along the  depth of an element, as is the case in open channels with slam curvature radii.  Page 32 of 56

Iber. Turbulent Flow in Shallow Waters

The terms E and D respectively model those bedload grains which become suspended (entrainment) and  deposition  of  suspended  sediments  to  the  bed  layer.  The  difference  represents  a  balance,  and  thus  a  coupling of the bedload and suspended load. 

4.3.2. Entrainment/Deposition (E‐D) term Calculation   Three  formulas  are  implemented  for  calculating  the  entrainment/deposition  term  (E‐D):Van  Rijn    (van  Rijn, 1987), Smith (Smith, J.D., 1977) and Ariathurai and Arulanandan (Ariathurai & Arulanandan, 1978).  The first two are valid on sand bed layers, while the Ariathurai is valid for cohesive bed layers. The three  formulas are especially recommended in the latest ASCE Manual of Sediment Transport, among them,  the most widely used is the Van Rijn Formula.   Van Rijn  In the Van Rijn formula  (van Rijn, 1987), the term E‐D is evaluated from the following expression 

E  D  Ws  c*a  ca   α Ws  C*  C    Where  α  is a coefficient that relates the mean suspended particle concentration and the river bedload  concentration,  which  value  comes  from  the  Rouse  profile  for  the  depth  sediment  concentration  distribution, Ws is the sediment grain velocity, C is the depth‐averaged suspended load concentration, 

C* is  the  depth‐averaged  suspended  load  concentration  at  equilibrium  condition  (suspended  load  transport  capacity), 

c a and c*a are  respectively  the  instantaneous  concentration  and  the  equilibrium   

concentration at a height z=a above the river’s bed layer, where a is the height of the layer in which the  transport  is  produced  (theoretical  separation  limit  between  the  bedload  and  suspended  load).  This  height must be approximated from the grain diameter. The coefficient  α  is calculated from the vertical  concentration distribution (Rouse profile), from the following integral: 

h-a

α=



h

a

 h-z a     z h-a 

ws

k u *

a = 3  D50   dz

Where κ=0.41 is the von Karman constant.  The equilibrium concentration near the bed layer proposed by Van Rijn (van Rijn, 1987) is: 

Page 33 of 56

Iber. Turbulent Flow in Shallow Waters

c*a  0.015

a = ks

D50  T1.5   a  D*0.3

ks  3  Ds  

1/3

gR D*  D   2     ν 

 

 Smith  This formula is similar to the Van Rijn, its only difference lies in the formula to calculate the equilibrium  concentration, for which Smith (Smith, J.D., 1977) proposes the following formula: 

1.56 103  T c    1 + 2.4 103  T * a

a = 26.3   τ*s - τ*c   Ds + k s

k s  3  Ds  

 Ariathurai y Arulanandan  For  cohesive  soils  the  formula  proposed  by  Ariathurai  and  Arulanandan  (Ariathurai  &  Arulanandan,  1978) is used, in it the erosion depends on the difference between the tangential stress and the critical  tangential stress to start the erosion  ce , it’s also dependant of a value M representative of the erosion  rate (equivalent to an erosion rate when  τ b =

2  τ ce ): 

τ  E = M   b  1    τ ce  In  cohesive  soils  a  modification  is  introduced  for  calculating  D,  which  considers  the  deposition  critical  tangential stress  cd  . In this case: 

D  P  α Ws C   With:  

 τ  P  1  b        if     τ b < τ cd   and  P  0  otherwise   τ cd 

Page 34 of 56

Iber. Turbulent Flow in Shallow Waters

4.3.3. Sedimentation Velocity  The sedimentation velocity is calculated based on its diameter (van Rijn, 1987):  2 R  g  D50 18  υ 10  υ Ws = 1+0.01  D*3 -1 D50

Ws =



Ws =1.1 R  g  D50



 D50 < 10-4 m  10-4 m < D50 0

 

Page 41 of 56

Iber. Turbulent Flow in Shallow Waters

5.2.2. Discretization of the bed slope source term   Iber uses a centred discretization of all the source terms except of the bed slope term. The main cause  of using an upwind discretization of the bed slope instead of a centred discretization is that it calculates  exactly the hydrostatic solution in an irregular bathymetry, this way avoiding spurious oscillations in the  free surface waters and in the velocities. These oscillations are in general small, but they can become  larger  in  magnitude  in  irregular  bathymetry  problems,  such  is  the  common  case  in  river  and  coastal  hydraulics.  Whenever  the  upwind  numerical  schemes  are  used  to  discretize  convective  flow,  the  upwind  discretization of the bed slope provides has better properties and provides more precise results than the  classic centred schemes.  The used discretization of the bed slope source term in a finite volume Ci can be expressed as: 

Si =  S dA = Ci

S jKi

ij

 

Where Sij is an upwind discretization of the bed slope source term in each edge of the considered finite  volume and is calculated as: 

 0  | nij | hL  hR  -1 -1  S ij = -g zb , R  zb , L   I  X | D | D X   n x ,ij     2 2  n   y ,ij  The same as in the convective flux, to implement this discretization in Iber, it is decomposed using the  vectors em as:  3

Sij =   mem   m 1

1 2

1C   cij zij | n ij |  2C  0

 

1 2

 3C  cij zij | n ij |

Page 42 of 56

Iber. Turbulent Flow in Shallow Waters

5.3. Discretization of the Transport Equations in the K‐E Turbulence Model, and  in the Suspended load Transport Model  The  conservation  equations  used  in  the  k‐ε  turbulence  model  and  in  the  suspended  load  transport  model can be expressed as: 

 Fj   S  t x j Where   = C h   is  the  conserved  variable,  C  is  the  non‐conserved  variable,  corresponding  to  the  considered  model  (sediment  concentration,  kinetic  turbulent  energy,  turbulent  dissipation  rate),  S  are  the source terms which model generation or destruction processes of the conserved variable, and F is  the convective and diffusive flux, which may be expressed as: 

Fj  C h u j   e h

C   x j

Where Γe is the effective diffusivity coefficient, which includes the molecular and turbulent diffusion. In  general,  turbulent  diffusion  is  many  times  larger  than  the  molecular  diffusion,  meaning  that  this  term  can be omitted.  In  the  k‐ε  turbulence  model  it  is  necessary  to  solve  2  transport  equations;  one  due  to  the  turbulent  kinetic energy (k) and the other one for the turbulence dissipation rate (ε). 

5.3.1. Depth‐averaged Transport Equation   To spatially discretize the transport equation using the finite volume method the differential equation is  integrated  in  each  cell  of  the  computational  mesh.  This  process  is  especially  advantageous  for  solving  the  conservation  equations,  since  they  are  being  solved  in  an  integral  form,  which  allows  the  formulation  of  natural  conservative  methods.  Integrating  the  convection‐diffusion  equation  in  a  two‐ dimensional cell we obtain: 

 C h i

n 1

  C h i

Δt

n

A i      C h  u  dA       e h  C  dA   Ai

Ai

Ai

E - D

dA  

Where  Φ=C.h  is  the  average  value  at  the  cell  and  S=E‐D  is  the  source  term  due  to  generation/destruction  processes  of  the  considered  variable.  Applying  the  divergence  theorem  to  the  second and third term we obtain: 

Page 43 of 56

Iber. Turbulent Flow in Shallow Waters

 C h i

n 1

  C h i

n

Δt

A i   C h  u  n dL    e h   C   n dL   Li

Li

Ai

E - D 

dA  

Where  the  integrals  around  the  cell’s  area  are  extended  to  the  cell’s  edges.  This  equation  is  the  conservation  equation  expressed  in  an  integral  form.  The  integrals  that  appear  in  the  integral  form  conservation equation are solved discretely, this way getting: 

 C h i

n 1

  C h i

Δt

n

Ai 

 C h  u  n  L j K i

ij

ij



   h  C   n  L   E - D  A   j K i

ij

ij

i

i

Where the summations are extended to all the edges that form the boundary of the cell. The sub‐index  ij identifies the common face of the cells i and j. Each term of the summations represent a flux of the  considered  variable  that  enters  the  cell  through  the  corresponding  face,  that  way  the  sum  of  all  the  fluxes through all the faces that form the contour of the cess are equal to the balance of what goes out  minus what goes in, for example; the net flux outside the cell. In function of the used interpolation to  calculate the flux through the cell’s contours, especially the convective flux, different numerical schemes  can be found.  The  diffusive  flux  between  cells  is  calculated  using  a  centred  discretization  of  second  order,  without  numerical stability problems. The convective flux discretization is more problematic from the numerical  stability point of view. In the next paragraphs some implemented numerical schemes are presented in  the code to discretize these terms.   First Order Upwind scheme   In the upwind schemes the direction of the information flow is considered to make the discretization of  the convective flux. The diffusive flux is discretized using a centred scheme of second order. In the case  of the convective flux, the information travels in the same direction as the water velocity, for example  upstream  to  downstream.  The  simplest  upwind  method  is  the  first  order  upwind  scheme,  calculated  using: 

 C h  u  n ij   u  n ij  C h ij  un,ij   C h ij     The  velocity  is  discretized  in  a  centred  manner,  while  the  value  of  C.h  is  discretized  in  an  upwind  manner, considering the upstream node value: 

Page 44 of 56

Iber. Turbulent Flow in Shallow Waters

u n ,ij  α u n ,i  (1  α) u n ,j

 C h ij   C h i  C h ij   C h  j

si u n ,ij  0   si u n ,ij  0

Where  α  is a linear interpolation coefficient. For an equidistant discretization the value of  α  is 0.5.  The fact of upwinding the convective term in the first order scheme is equivalent from the mathematical  point of view as adding a diffusive term (usually called numerical diffusion) with a diffusivity (numerical)  coefficient  Γn  proportional  to  the  mesh  size.  The  numerical  diffusion  is  a  consequence  of  the  used  technique for numerical stabilization. The desired situation in a numerical scheme is for it to be stable  with the introduction of as little as possible numerical diffusion. This is why it’s so convenient to use fine  meshes to reduce the error introduced by numerical diffusion   Second Order Upwind scheme  The  previously  mentioned  scheme  is  of  first  order  because  the  numerical  diffusion  introduced  in  the  discretization of the convective flux. Despite this it is widely used in Computational Fluid Dynamic (CFD)  codes as the default scheme, due to its numerical stability. When more precision is required in a regular  size mesh, it is necessary to appeal to higher order schemes. To do this a second order scheme of type  MUSCL (Monotonic Upstream Scheme for Conservative Laws) has been implemented Iber, which makes  a linear reconstruction of each mesh element of non‐conserved variable.  Once this linear reconstruction of the non‐conserved variable is made, the convective flux is calculated  at each edge of the computational mesh with: 

 C h  u  n ij   u  n ij  C h ij  un,ij   C h ij u n ,ij  α u n ,i  (1  α) u n ,j

 C h ij   C h Ij  C h ij   C h iJ

si u n ,ij  0

 

si u n ,ij  0

Where   C h  I j  and   C h i J  the values of the variable at edge Lij taken from the linear reconstruction  of cells Ci y Cj respectively. 

5.4. Discretization of the Exner Conservation of Sediment Equation  The Exner equation of conservation of sediment can be as:  Page 45 of 56

Iber. Turbulent Flow in Shallow Waters

(1  p)

z b q sb,x q sb,y    DE  t x y

The  Exner  equation  is  discretized  in  a  similar  way  as  the  suspended  load  sediment  conservation  equation. The integral of the Exner equation in each two‐dimensional cell can be written as: 

(1  p)

n 1 n  zb,i zb,i

Δt

q sb,y  q A i    sb,x + Ai y  x

  dA   D - E i A i   

Applying the divergence theorem to the second term of the equation, we get: 

(1  p)

n 1 n zb,i  zb,i

Δt

Ai 

 q n  L   D - E  A   jK i

* sb

ij

ij

i

i

The value of the bedload in each edge of the mesh   q *sb  is calculated in an upwind form as:  ij

q  q 

* sb ij

 q sb,i

si  q*sb n   0

* sb ij

 q sb,j

si  q n   0

ij

* sb

 

ij

The erosion/deposition source term is discretized in a centred form in each mesh cell. 

5.4.1. Non‐erodible bed layer  There is the possibility to include a rock layer, meaning a non‐erodible bed layer. The implementation in  the code of this non‐erodible bed layer has been done using a 2 stage calculation (predictor‐corrector)  of the numerical bedload fluxes through the element’s contours.  Predictor: searches a new bed layer elevation due mainly to the outflow in an element. 

z ip = z in +

Δt (1  p)  A i

  max 0, q jK i

* sb



n  Lij   ij

Corrector:  if  at  the  predictor  the  element  has  eroded  below  the  rock  layer,  the  bedload  outflow  is  corrected en each of the edges, and afterwards it is updated to the new bed layer elevation.  p

Si (zi

< z roca ) y (  q*sbn  > 0) :  ij

Page 46 of 56

Iber. Turbulent Flow in Shallow Waters

q  = q  *,c sb ij

q 

* sb ij

*,c sb ji

(1  p)

n 1 n z b,i  zb,i

Δt

Ai 

z in -z roca z in -z ip

= -  q*sb 

 

ij  

 q n  L   D - E  A   jK i

*,c sb

ij

ij

i

i

5.5. Managing Wet‐Dry Fronts  Modeling  flooding  zones,  as  is  the  movement  of  wave  fronts  in  estuaries  and  coastal  zones  is  fundamental  in  environmental  hydraulics.  In  Iber  wet‐dry  fronts,  both  stationary  and  non‐stationary,  that may appear in the domain are modelled with a fixed finite volume mesh, allowing each element to  have water or not depending on the flow conditions. Among the volumes which do not have water and  those which do have, there is a dry‐wet front which is necessary to treat appropriately from a numerical  point  of  view  to  avoid  instabilities  and  non‐physical  oscillations  in  the  results.  For  this  wet‐dry  front,  either in a flood front or a wave front, there is a wet‐dry tolerance εwd, so if the water depth is lower  than εwd it is considered that the cell is dry and is not included in the calculations. The wet‐dry tolerance  can be changed by the user to a value close to zero, even though in irregular bathymetry problems, as is  very  usual  in  river  an  coastal  hydraulics,  it  is  advisable  to  use  values  in  the  order  of  1cm  or  1mm  to  increase the stability of the results without spoiling their precision. In any case, the water height is never  forced  to  be  equal  to  zero,  with  the  purpose  of  avoiding  mass  losses  in  the  calculation  domain.  The  numerical scheme used to solve the wet‐dry front is stable and non‐diffusive.  The management the wet‐dry fronts in Iber is stable, conservative and non‐diffusive, meaning that the  fronts  are  adequately  solved  without  instabilities  of  the  numerical  type,  even  when  these  occur  in  strong bed slopes. Every finite element has an associated bed level. Schematically shown in Figure 9. 

 

Page 47 of 56

Iber. Turbulent Flow in Shallow Waters

Figure 9. Schematic representation of the wet‐dry bed level management. 

Between two volumes with different bed level one of the following situations may occur, as presented  in Figure 10: 

   

   

 

Figure 10. Several situations of the water levels between two adjacent cells. 

In  the  first  figure  both  volumes  have  water,  so  no  front  is  produced  and  thus  no  special  treatment  is  necessary. In the other two cases there is a wet‐dry front. The difference is that in the second case, the  free water level in the wet cell is higher than the bed level of the dry cell, while in the third case it is  lower. Only in the third case it is necessary to use a special treatment, which consists of redefining the  bed slope to impose a reflexive condition in the front. In this case the bed slope is redefined as: 

h i - h j Δz b,ij =   z b,j - z b,i

si

h j  z b,j - z b,i

si

h j > z b,j - z b,i

 

The reflexive condition is imposed as: 

qn,ij = q x,ij n x,ij + q y,ij n y,ij  0   The  use  of  the  mentioned  conditions  provides  can  reproduce  the    hydrostatic  solution  for  any  bathymetry, without streatching the front and without generating any spurious oscillations in the free  surface.  This  type  of  wet‐dry  front  treatment  has  been  successfully  used  in  the  modeling  of  both  stationary  and  non‐stationary  processes,  being  particularly  useful  in  the  simulation  of  flood  zones  in  rivers and coastal zones, and in calculating the evolution of wave fronts.  A more detailed description of the implementation in Iber can be found in the references provided in  section 6.     

 

Page 48 of 56

Iber. Turbulent Flow in Shallow Waters

         

  IBER  Shallow waters turbulent flow       

REFERENCES 

Page 49 of 56

Iber. Turbulent Flow in Shallow Waters

6.  REFERENCES  Apsley, D., & Stansby, P. (2008). Bed‐Load Sediment Transport on large slopes: Model Formulation and  Implementation within a RANS solver. Journal of Hydraulic Engineering, ASCE , vol 134 (10).  Ariathurai, R., & Arulanandan, K. (1978). Erosion rate of cohesive soils. Journal of the Hydraulics Division.  ASCE , 104 (HY2), 279‐283.  Bermudez,  A.,  Dervieux,  A.,  Desideri,  J.,  &  Vazquez‐Cendon,  M.  (1998).  Upwind  schemes  for  the  two  dimensional  shallow  water  equations  with  variable  depth  using  unstructured  meshes.  Computational  Methods Applied for Mechanical Engineering , Vol 155.  Bladé,  E.,  &  Gómez‐Valentín,  M.  (2006).  Modelación  del  flujo  en  lámina  libre  sobre  cauces  naturales.  Análisis integrado en una y dos dimensiones. Monografía. Barcelona: CIMNE.  Bladé,  E.,  Gómez‐Valentín,  M.,  Sánchez‐Juny,  M.,  &  Dolz,  J.  (2008).  Preserving  steady  state  in  Finite  Volume Computations of River Flow. Journal of Hydraulic Engineering, ASCE , Vol 134 (9).  Cea, L. (2005). Unstructured finite volume model for unsteady turbulent shallow water flow with wet‐dry  flows: numerical solver ans experimental validation. Phd thesis. A Coruña: Universidad de A Coruña.  Cea,  L.,  French,  J.,  &  Vázquez‐Cendón,  M.  (2006).  Numerical  modelling  of  tidal  flows  in  complex  estuaries  including  turbulence:  An  unstructured  finite  volume  solver  and  experimental  validation.  International Journal for Numerical Methods in Engineering , vol 67 (13).  Cea, L., Puertas, J., & Vázquez‐Cendón, M. (2007). Depth averaged modelling of turbulent shallow water  flow with wet‐dry fronts. Archives of Computational Methods in Engineering, State of the art reviews ,  vol 14 (3).  Chiew,  Y.,  &  Parker,  G.  (1994).  Incipient  Sediment  Motion  on  Non‐Horizontal  Slopes.  Journal  of  Hydraulics Research , 649‐660.  Chow, V., Maidment, D., & Mays, L. (1988). Applied Hydrology. McGraw‐Hill.  CIMNE. (2009). GiD The personal pre and post‐processor. Barcelona: www.gidhome.com.  Corestein, G., Bladé, E., Gómez‐Valentín, M., Dolz, J., Oñate, E., & Piazzese, J. (2004). New GiD interface  for Ramflood‐DSS project hydraulic simulation code. Proceeding of the 2nd Conference on advances and  applications of GiD. Barcelona: CIMNE.  Dey.  (2003).  Threshold  of  sediment  motion  on  combined  transverse  and  longitudinal  sloping  beds.  Journal of Hydraulic Research , #4 405‐415.  Page 50 of 56

Iber. Turbulent Flow in Shallow Waters

García, M., & Parker, G. (1991). Entrainment of Bed Sediment into suspension. ASCE Journal of Hydraulic  Engineering , #4 414‐435.  Meyer‐Peter, E., & Muller, R. (1948). Formulas for bedload transport. Proceeding of the 2nd Conference  (pp. 39‐64). Stockholm: IAHR.  Rastogi,  A.,  &  Rodi,  W.  (1978).  Predictions  of  heat  and  mass  transfer  in  open  channels.  Journal  of  Hydraulics Division , 104 (3) 397‐420.  Smith,  J.,  &  McLean,  S.  (1977).  Spatially  averaged  flow  over  a  wavy  surface.  Journal  of  Geophysical  Research , 82(12) 1735‐1746.  Smith,  J.D.  (1977).  Modelling  of  sediment  transport  on  continental  shelves.  The  Sea:  ideas  and  observations on progress in the study of the seas. New York: John Wiley and Sons.  Van Dorn, W. (1953). Wind Stress on an artificial pond. Journal of Marine Research , Vol 12.  van  Rijn,  L.  (1987).  Mathematical  modelling  of  morpholigical  processes  in  the  case  of  suspended  sediment transport. Delft: Delft Hydraulics Communications #382.  van  Rijn,  L.C.  (1984).  Sediment  Transport,  part  I:  Bedload  transport.  Journal  of  Hydraulic  Engineering  ,  vol 110 (10).  Wong, M. (2003). Does the bedload equation of Meyer‐Peter and Muller fit its own data? 30th Congress  (pp. 73‐80). Thessaloniki: International Association of Hydraulic Research.  Wong, M., & Parker, G. (2006). Reanalysis and Correction of Bedload relation of Meyer‐Peter and Muller  using their own database. Journal of Hydraulic Engineering , 132(11) 1159‐1168.  Wu, W. (2004). Depth‐averaged two‐dimensional numerical modeling of unsteady flow and nonuniform  sediment transport in Open Channels. Journal of Hydraulic Engineering , 130(10) 1013‐1024. 

7.

Page 51 of 56

Iber. Turbulent Flow in Shallow Waters

NOMENCLATURE  a width of the layer in which bedload transport occurs 

c a ,  c*a  instant concentration and equilibrium concentration of suspended solids at a height z=a above  the river bed level   C depth‐averaged suspended load  

C* Depth‐averaged suspended load in equilibrium conditions   Cd discharge coefficient at a gate or a weir  Cf bed friction coefficient  Cm adjust coefficient for the turbulent viscosity in a parabolic model  dwall wall distance  Dij lateral dispersion terms  Ds sediment diameter  D50 mean sediment diameter 

D*  Non‐dimensional sediment diameter  Dsx,  Dsy  Suspended  load  dispersion  due  to  the  non‐homogeneity  of  the  velocity  profiles  and  of  the  vertical direction sediment concentration  D‐E balance between the bedload and suspended load  f potential infiltration rate  g gravity acceleration  h depth  i real infiltration rate  k turbulent kinetic energy  ks soil’s saturated permeability  Ks grain’s roughness height   kvd wind drag coefficient 

Page 52 of 56

Iber. Turbulent Flow in Shallow Waters

L0 Initial depth of the soil’s saturated region  L Depth of the soil’s saturated region   M2 parameter that measures the erosion rate because of re‐suspension   Ms mass source/drain term   Mx, My momentum source/drain terms  n Manning coefficient  ns equivalent Manning coefficient due to the grain  p porosity of the bedload sediments  qn unit normal flow at the inlet contour  qsb,x, qsb,y components of the bedload discharge  Q discharge/Flow  Rh hydraulic radius  R non‐dimensional submerged specific weight   Sc,t Schmidt number  Se Effective initial soil saturation  Sij deformation tensor   T non‐dimensional parameter which measures the excess of friction at the bed layer above the critical  value that defines the start of movement  uf friction velocity due to bed roughness  

u'i u' j  Turbulent stresses or Reynolds stresses   Ux, Uy depth‐averaged horizontal velocities  |U| depth‐averaged mean velocity  V10 Wind velocity at 10 meters height  Ws sedimentation velocity of solid particles  Zs Free layer level  Zb Bed layer level 

Page 53 of 56

Iber. Turbulent Flow in Shallow Waters

α coefficient  that relates the mean  concentration of suspended solids and the concentration near the  river bed  β Bed slope  Г Molecular diffusion of suspended solids coefficient  δij Kronecker delta  ε turbulence dissipation rate  Δθ changes in the soil’s humid content as the saturation advances  θi initial soil’s humidity content   θe soil’s effective porosity (drainable)   θr retention capacity (irreducible humidity or non‐drainable)  κ=0.41 von Karman constant  λ latitude of the considered point  μc bed material’s internal friction angle   ν fluid’s cinematic viscosity  νt turbulent viscosity  ρ water density 

ρs Sediment density   τ b Total stress due to bed friction    τ bs  Bed stress due to grain  τ c  Critical bed stress  τ*b , τ*bs  Total grain stresses (non‐dimensional)  τ*c  Critical bed stress (non‐adimensional)  τs free surface friction due to wind induced roughness  τexx, τexy, τeyy tangential effective horizontal stresses 

τijv  Viscous stresses   Page 54 of 56

Iber. Turbulent Flow in Shallow Waters

 Soil total porosity  ψ soil non‐saturated región suction   Ω earth’s rotation angular velocity  Ws sediment velocity     

Page 55 of 56