PREDICTION OF AFFLUX AT BRIDGE CONSTRICTIONS USING A

DEPTH AVERAGED NUMERICAL MODEL

Mohammad Reza Hadian(1), Amir Reza Zarrati(2)
(1)

Assistant Professor, Department of Civil Engineering, Yazd University, Safaiyeh, Yazd, Iran
phone: +98 0351 8211671; fax: +98 0351 8210699; e-mail: [email protected]
(2)
Associate Professor, Dept. of Civil Eng., Amirkabir University of Technology, 424, Hafez Ave., Tehran, Iran
phone: +98 021 64543002; fax: +98 021 66414213; e-mail: [email protected]

ABSTRACT
Calculation of afflux and backwater at bridge constrictions is one of the most important
factors which should be considered in hydraulic design of bridges over rivers. The 2D depth
averaged form of shallow water equations is employed in the present study to calculate
backwater upstream of a constriction. The equations are solved on unstructured blocks with
structured grids used in each block. The model uses collocated grid arrangement and
SIMPLEC like algorithm for velocity-water surface coupling. As a benchmark, a bridge
constriction was simulated in a laboratory flume and was tested for common range of Froude
numbers in engineering projects. The water surface profiles were measured at center line of
the channel and some cross sections upstream of the constriction and were compared with the
results of the numerical model. The comparison showed that both zero equation and k-
turbulence models which are implemented in the numerical model give accurate results in all
Froude numbers tested. It is also shown that the numerical model is more accurate than the
empirical method used even though the tests were conducted in a simple straight channel with
vertical abutments.
Keywords: afflux, bridge constriction, shallow water, numerical model, depth averaged.
1 INTRODUCTION
Bridges across waterways are grade separation structures like those built to carry traffic
over or under other highways. The important difference is that there is an interaction between
the river and bridge structure which should be considered in the bridge design. Hydraulic
design of bridges can be carried out by an engineer, which is trained in hydraulics, hydrology
and river behavior.
An important consideration in the hydraulic analysis of a proposed bridge crossing is the
amount of afflux and backwater which may occur by constricting the flow with the highway
crossing. Backwater can cause flooding upstream of the constriction and hamper the usage of
the highway.
Different methods can be employed in afflux prediction. The simplest method is
empirical equations derived by simplifying the classical hydraulic equations and applying
some coefficients to make them applicable for these sorts of calculations. There are many
tables and graphs to determine such coefficient which are found by experiment in laboratory
channels. Among different methods those presented by Chow (1959) and Bradley (1973) can
be mentioned here. However, because of many factors that can affect the afflux at a bridge
constriction, these equations and graphs are limited to cases which they are determined for.
With the rapid progress of computer science and technology, equations of motion can be
solved for any bridge crossing to determine afflux and backwater. In the present study the 2D
shallow water equations are used for this purpose.
In the present paper, details of the numerical model are briefly reviewed. Then the
experimental setup and measurements are explained. The results of the model are then
compared with experimental data and empirical method of Bradley (1973).
1

2 NUMERICAL MODEL
2.1 GOVERNING EQUATIONS
Conservative form of shallow water equations in general curvilinear coordinate system is
used in the present study to give the model the flexibility of dealing with complex geometries.
Neglecting wind shear stress, Coriolis acceleration, using Boussinesq approximation for
Reynolds stresses and applying partial transformation for mapping from Cartesian to
curvilinear coordinate with Cartesian velocities as primary variables, the equations will be as
(Hadian and Zarrati, 2005):

+ (Uh ) + (Vh ) = 0
t J

g 22 (h) g11 (h)

+ S
J (h) + (hU ) + ( hV ) =
+
t

J
J

(1)
(2)

Figure 1. Schematic control volume and notation.

In which, U and V are depth averaged contravariant velocity components in general

curvilinear and directions respectively (Figure 1), is transport parameter such as U and
V, h= water depth, =water density, = depth averaged viscosity which is equal to e for
momentum equations and equal to e for any other transported parameters. e = depth
averaged effective viscosity ( e = + t ), = kinematic viscosity of fluid, t = depth averaged
turbulent viscosity, g=gravitational acceleration, =water surface elevation ( = h + Z b ),
Z b =bed elevation, J= jacobian of transformation ( J = x y x y ) where x , x , y and y
are metrics, bx and by =bed shear stresses in x and y directions which can be calculated
using Mannings equation as:

bx g n 2 u u 2 + v 2
=
1

h3
by g n 2 v u 2 + v 2
=
1

h3

(3)
(4)

where u and v are depth averaged Cartesian velocities. The depth averaged contravariant

velocity components, which are tangential to curvilinear coordinate and respectively, are
defined as:

U = uy vx , V = vx uy

(5)

Other parameters are defined as below:

g 22 = x2 + y 2 , g12 = x x + y y , g11 = x2 + y2

(6)

(Uh )
(Uh )
(hu )
J bx + e y

y
y
g12
SU = gh y

J

(Vh )
(Vh )
(hu )

+ e y
y
g12
J

e (Uh )


(Uh )
(hv )
J by

SV = gh x
x
+ x
g12
-x

J

(Vh )
(Vh )
(hv )

+ e -x
+ x
g12
J

(7)

(8)

The depth averaged turbulent viscosity can be considered as a constant value,

calculated by a zero-equation models or more advanced models like depth-averaged k-. Here
the zero equation model of Jia and Wang (1999) and depth-averaged k- model are used. The
former can be written as:
t =

u* h
6

(9)

In which u* = bed shear velocity and is the von Karman constant (=0.4).
In depth averaged k- model t is calculated as (Rastogi and Rodi, 1978):
t = C

k2

(10)

The equations of k and can be expressed in general form of transport equations (2)
when is substituted by k or and with using the following source terms:
g12 (hk ) t g12 (hk )

+
S k = J (G + Pk h ) t

k
k

2 g (h ) t g12 (h )
+

S = J C1 G + P C 2 h t 12
k
k J J

t
2T12 + 2T22 + T32
h
1 (uh )
(uh )

y
T1 = y

J
G=

(11)
(12)
(13)
(14)

T2 =
T3 =

1
(vh )
(vh )
x

+ x

(15)

1
(vh )
(vh )
(uh )
(uh )
x

y
+ y
+ x

(16)

u*3
Pk =
cf
P = C

(17)

u*4
h

(18)

Other coefficients of k- equations are tabulated in Table 1(Launder and Spalding, 1974).
Table 1. constant coefficients of k- equations.

C1

C2

0.09

1.43

1.92

1.3

2.2 SOLUTION ALGORITHM

The equations are solved on structured grid using control volume method in which the
equations are integrated on control volumes (cells) in time and space. (Versteeg and
Malalasekera, 1995). Applying Gauss diversion method to convert surface integrals to
boundary integrals and the first-order implicit Euler method for time integration, the
momentum equations are discretized in the form of:
AP P = AE E + AW W + AN N + AS S + Apo op + S

(19)

In which the Ai s are coefficients and the indices are defined considering Figure 1. These
coefficients depend on the applied method for discretizing convection-diffusion terms. Here
the Power-law scheme is used in numerical model. The east and west factors for power-law
schemes can be calculated as:
0.1 F 5
e

AE = De max 0,1
(20)
+ max( Fe ,0)

D
e

5
0.1 F
w

AW = Dw max 0 ,1
(21)
+ max(Fw ,0)

D
w

Fe = (hU )e , Fw = (hU )w
(22)
g 22 h
De =
J
and finally:

g 22 h
, Dw =
e
J

(23)

F
644
47
444
8
AP = A + AE + AW + AN + AS S (Fe Fw + Fn Fs )
o
p

A po = J h po

1

t

(24)
(25)

where the superscript o denotes the values at previous time step.

For better convergence of equations the under-relaxation factors are used by replacing
equation (19) with:
1
AP
P = AE E + AW W + AN N + AS S + S + Apo op +
AP **p
(26)

In the above equation the ** superscript denotes the values of last iteration.
In which is the under-relaxation factor for equation of .
The water surface is calculated by a SIMPLEC like algorithm following Zhou (1995)
and Weerakoon et al. (2003), which was developed for non-orthogonal curvilinear coordinate
system by Hadian and Zarrati (2005). Therefore the corrections of contravariant velocities are
calculated as:

g h
, B=
(27)
U = B
g 22
Ap

Anb
v
g h

, C=
(28)
V = C
g11
Ap

Anb
v
The corrections of water surface elevation are calculated by solving the following
convection-diffusion type equations:
*
*

*
m
+
J
U +
V = Bh*
(29)
+ Ch
t

* * *
* *
(30)
+
U h +
V h
t

Though curvilinear coordinate system is a robust tool to handle complex geometry

domains, it results in a very skewed mesh in some regions when the flow domain is very
complex as it is in the case of bridge constriction. The multi-block technique is one of the
most powerful methods which can solve this problem. Here the method of halo cells is used
(Lien et al,, 1996; Apsley and Hu, 2003; Hadian and Zarrati, 2005) and the equations are
solved in unstructured blocks with one layer thick halo cells.
For the sake of simplicity in programming and especially for implementing the multiblock method, collocated grid has been used, so all the variables are stored at cell centers. To
avoid the checkerboard water surface fluctuation, it is necessary to use momentum
interpolation method for calculation of contravariant velocities at cell faces in solution of
equation (29) and coefficients of transport equation (2). In the present study the method of
Lien and Leschziner (1994) was developed for shallow water equations (Hadian and Zarrati,
2006).
In the test cases of present study the flow is sub-critical in both inlet and outlet. So, the
water surface elevation is used at outlet as boundary condition while the discharge is
introduced at inlet. The wall function method is used at solid walls to bridge the boundary
layer and the necessary source terms are added to the equations (Versteeg and Malalasekera,
1995).
At the beginning of calculation, the water surface level at the outlet section was used for
all the flow domain as initial condition and velocity in the flow domain was calculated using
known discharge and water surface level.
The procedure for solution of equations at every time step can be summarized as follows:
1. Solve the momentum equations for all blocks assuming fixed velocities at halo cells
m=J

2.
3.
4.
5.
6.
7.

(frozen halo cells) in each block and then update the velocities in the entire flow
domain including in halo cells.
Solve the depth correction equation (29) for all blocks using frozen halo cells and then
update in halo cells.
Correct the water surface elevation, water depth and all velocities at all blocks.
Calculate eddy viscosity by Equation (9) or by solving k- equations at every block
and update the halo cells.
Considering the inlet discharge correct the velocities at inlet section.
Repeat the above steps till convergence is achieved.
Go to the next time step.

Inlet tank

Tail gate

Bridge Axis

Tailwater control section

3 EXPERIMENTAL SETUP
The experiments were conducted in hydraulic laboratory of Amirkabir University of
Technology. The experimental flume was 6 m long, 0.513 m wide and 0.75 m deep and fed
with 2 centrifugal pumps. The bottom of flume was made of painted steel and walls were
Perspex material. The bridge constriction was simulated by 2 abutments with length of 10 cm
and thickness of 4.2 cm, made up of painted wood. By installing these abutments the bridge
opening ratio of 39% was formed in the flume. Water pumped into an inlet tank and flowed
into flume after passing some screens that produced almost uniform flow across the channel
width. There was a tailgate at the end of the flume which was used to adjust the downstream
tail water depth (Figure 2).

C6
C5
C4
C3 C2
C1

0.513

2.00

3.20
6.00

Figure 2. Layout of the experimental flume and transverse sections position.

Experiments were performed for 6 cases with different flow discharge and tail waters as
it is summarized in Table 2. With discharges and depths selected, the range of Froude number
was between 0.15 to 0.3 similar to mild slope rivers. Figure 3 shows flow condition at the
constriction in one of these cases.
Table 2. Main parameters in the experiments.

Test No
T1 T2 T3 T4 T5 T6
Discharge (lit/s)
10 10 15 15 20 20
Downstream tail depth (cm) 10 12 10 12 12 15

Figure 3. Flow at constriction with Discharge=15 lit/s and tail water =12 cm ( (view looking upstream).

4 RESULTS AND COMPARISONS

To reduce the CPU time, only half of the channel was modeled and the symmetric
condition was applied at the channel axis. As it is presented in Figure 4, five blocks was used
to generate the mesh in the flow domain.

Figure 4. Mesh applied in numerical model.

The water surface profile along the centerline of channel is presented in Figure 5 for
different test cases. As the figure shows, the model could predict accurately the water surface
location at the upstream side of the bridge. Both the zero equation and k- model could
predict the water surface elevations with almost the same accuracy. Downstream of the
constriction, the pressure distribution deviates from hydrostatic and the results of the shallow
water model deteriorate. Deviation from experimental results downstream of the constriction
is less in lower values of afflux.
Figure 6 shows water surface profile at transverse sections C1 to C6 (Figure 2) for the
test with discharge of 15 lit/s and tail water level of 10 cm. Predictions of the numerical
model is very good except at section C1 which is the downstream edge of the abutment where
pressures deviate from hydrostatic distribution (Figure 3). It should be mentioned here that
empirical relations are unable to give these details on transverse distribution of afflux
upstream of the constriction
The empirical method of Bradley (1973) was used here for calculation of afflux in the 6
test cases and its results are compared with the numerical model. The results are presented in
Figure 7. In this figure the error in calculating afflux normalized by the tail water depth is
plotted verses Froude number. As the figure shows, the accuracy of the present model is much
better than the empirical equation especially at higher Froude numbers. Accuracy of both zero
equation and k- model is almost the same, although zero equation model is much simpler.
While the error of empirical equation increases with Froude number, the numerical model

error is similar in the range of Froude number tested.

Figure 5. Comparison of water surface profile along channel axis.

Figure 6. water surface profile at transverse sections.

Afflux Error (%)

16
14
12
10
8
6
4
2
0

zero equation model

k- model
empirical (Bradley, 1973)

0.1

0.2
Fr

0.3

0.4

Figure 7.Error in afflux calculation vs. Froude number.

5 CONCLUSIONS
A numerical model based on shallow water equations is developed in the present study to
simulate afflux upstream of constrictions. The model solves the transformed equation to
general curvilinear coordinate system by control volume method on unstructured Multi-Block
grids and structured grid in each block. Halo cell technique is used for data transformation
and collocated grids are applied to store the variables. The water surface is calculated using a
SIMPLEC like algorithm and both the zero equation and depth averaged k- model are
implemented for calculating turbulent eddy viscosity. Some experiments were carried out for
vertical abutments in a straight channel and water surface profiles were measured along center
line of the channel and some cross sections upstream of the constriction.
The model could predict the afflux accurately in all the test cases. The accuracy of the
model with k- and zero equation turbulence models was similar. Comparison of the present
model and empirical method of Bradley (1973) showed that the present numerical model is
more accurate even in the simple geometry used in the experiments which these equations
were derived for. Moreover, the numerical model gives more details about water surface
profiles in transverse direction upstream of the constriction. The accuracy of empirical
equations decreased as the Froude number increased however the error of the numerical
model was almost similar for different Froude numbers.
ACKNOWLEDGMENTS
The authors would like to appreciate Dr. David Apsleys comments during development
of the numerical model. Also, for experimental setup and measurements, kindly assist of Mr.
Karimi from the hydraulic laboratory of Amirkabir University of Technology is
acknowledged.
REFERENCES
Apsley, D. and Hu,W. (2003), CFD Simulation of Two- and three-Dimensional Free Surface
Flow , International Journal of Numerical Methods in Fluids, 42, pp. 465-491.
Bradley, J.N. (1973), Hydraulics of bridge waterways (Hydraulic design series), U.S. Federal
Highway Administration.
Choi, S.K. (1999), Note on the Use of Momentum Interpolation Method for Unsteady Flows,
Numerical Heat Transfer, Part A, 36, pp. 545-550.
Chow, V.T. (1959), Open-Channel Hydraulics, McGraw-Hill Book Company.
Hadian, M.R. and Zarrati, A.R. (2005), Simulation of Shallow Flows in a Circular Basin
Using Collocated Multi-Block Method, in Proceedings of XXXI IAHR Congress, Seoul,
9

Korea.
Hadian, M.R. and Zarrati, A.R. (2006), Study on the Effect of Momentum Interpolation
Method for Solving Shallow Water Equations, Amirkabir Journal of Science and
Technology, Vol.16, No.63- C(CMM), pp. 31-41 (in Persian).
Jia, Y. and Wang, S.S.Y. (1999), Numerical Model for Channel Flow and Morphological
Channel Studies, Journal of Hydraulic Engineering, ASCE, 125(9), pp. 924-933.
Lien, F.S. and Leschziner, M.A. (1994), A General Non-Orthogonal Collocated Finite
Volume Algorithm for Turbulent Flow at all Speeds Incorporating Second-Moment
Turbulence-Transport Closure, Part 1: Computational Implementation, Computer
Methods in Applied Mechanics and Engineering, 114, pp. 123-148.
Lien, F.S., Chen, W.L. and Leschziner,M.A. (1996). A Multiblock Implementation of NonOrthogonal Collocated Finite Volume Algorithm for Complex Turbulent Flows,
International Journal for Numerical Methods in Fluids, 23, pp. 567-588.
Majumdar,S., (1988), Role of Underrelaxation in Momentum Interpolation for Calculation of
Flow with Nonstaggered Grids, Numerical Heat Transfer, 13, pp. 125-132.
Miller, T.F. and Schmidt, F.W. (1988), Use of Pressure-Weighted Interpolation Method for
the Solution of the Incopressible Navier-Stokes Equations on a Nonstaggered Grid System,
Numerical Heat Transfer, 14, pp. 213-233.
Rastogi, A,K. and Rodi, W. (1978), Prediction of Heat and Mass Transfer in Open Channels,
Journal of Hydraulics Division, ASCE, 104(3), pp. 397-420.
Rhie,C.M. and Chow,W.L. (1983), Numerical Study of the Turbulent Flow Past an Airfoil
Trailing Edge Separation, AIAA Journal, 21(11), pp. 1525-1532.
Yu, B., Kawaguchi, Y., Tao, W.Q. and Ozoe, H. (2002a), Checkerboard Pressure Predictions
Due to the Underrelaxation Factor and Time Step Size for a NonStaggered Grid with
Momentum interpolation Method, Numerical Heat Transfer, Part B, 41, pp. 85-94.
Yu, B., Tao, W.Q., Wei, J.J., Kawaguchi, Y., Toshio, T. and Ozoe, H. (2002b), Discussion on
Momentum Interpolation Method for Collocated Grids of Incompressible Flow, Numerical
Heat Transfer, Part B, 42, pp. 141-166.
Versteeg, H.K. and Malalasekera, W. (1995), An Introduction to Computational Fluid
Dynamics, The Finite Volume Method, Longman Book Publisher.
Wang,Y. and Komori, S. (1999), Comparison of Using Cartesian and Covariant Velocity
Components on Non-orthogonal Collocated Grids, International Journal for Numerical
methods in Fluids, 31, pp. 1265-1280.
Weerakoon, S.B., Tamai N. and Kawahara Y. (2003), Flow Computation at a River
Confluences, Journal of Water and Maritime Engineering, Inst. of Civil Engineers, 1, pp.
73-83.
Zhou, Jian Guo (1995), Velocity-Depth Coupling in Shallow-Water Flows, Journal of
Hydraulic Engineering, ASCE, 121(10), pp. 717-724.

10