Turbulent Flow in Open Channels

Figure 147: Channel geometry

The turbulent flow in open channels can be approximated by one-dimensional network calculations. For the theoretical background the reader is referred to [17] and expecially [12] (in Dutch). [22] contains information on the solution of transient problems and (transient and steady state) analytical examples. The governing equation is the Bresse equation, which is a special form of the Bernoulli equation. For its derivation we start from Equation (613), which we write down for a flow line near the bottom of the channel in the form:

$$\frac{1}{g} \left [ \frac{\partial \Vert \boldsymbol{v} \Vert }{\... ...o g} t_{kl,l} \frac{dx_k}{ds} - \frac{1}{\rho g} \frac{dp}{ds} - \frac{dz}{ds}.$$ (621)

Assuming:

1. steady-state flow
2. each cross section is hydrostatic
3. the velocity is constant across each cross section
4. the velocity vector is perpendicular to each cross section,

one arrives at:

$$\frac{d}{ds} \left( \frac{Q^2}{2 g A^2} \right) = -S_f - \frac{dh}{ds} \sqrt{1-S_0^2} + S_0,$$ (622)

where (Figure 147) h is the water depth (measured perpendicular to the channel floor), s is the length along the bottom, $S_0= \sin(\phi)$, where $\phi$ is the angle the channel floor makes with a horizontal line, $S_f$ is a friction term (head loss per unit of length; results from the viscous stresses), $g$ is the earth acceleration, $Q$ is the volumetric flow (mass flow divided by the fluid density) and $A$ is the area of the cross section. This also amounts to:

$$\frac{Q}{gA^2} \frac{dQ}{ds} - \frac{Q^2}{gA^3} \left [ \left . \... ...}{\partial h} \frac{dh}{ds} \right ] + \frac{dh}{ds} \sqrt{1-S_0^2} = S_0 -S_f.$$ (623)

Assuming no change in flow (dQ/ds=0) and a trapezoidal cross section (for which $\partial A / \partial h = B$, where $B$ is the width at the free surface) one finally obtains (Bresse equation):

$$\frac{dh}{ds}=\frac{S_0 - S_f + \frac{1}{g} \frac{Q^2}{A^3} \frac {\partial A}{\partial s}}{\sqrt{1 - S_0^2} - \frac{Q^2 B}{g A^3}},$$ (624)

For $S_f$ several formulas have been proposed. In CalculiX the White-Colebrook and the Manning formula are implemented. The White-Colebrook formula reads

$$S_f=\frac{f}{8g} \frac{Q^2 P}{A^3},$$ (625)

where $f$ is the friction coefficient determined by Equation (163), and $P$ is the wetted circumference of the cross section. The Manning formula reads

$$S_f=\frac{n^2 Q^2 P^{4/3}}{A^{10/3}}$$ (626)

where $n$ is the Manning coefficient, which has to be determined experimentally.

In CalculiX, the channel cross section has to be trapezoidal (Figure 147). For this geometry the following relations apply:

$$A=h(b+h \tan\theta),$$ (627)

$$P=b+\frac{2h}{\cos\theta}$$ (628)

and

$$B=b+2h \tan\theta.$$ (629)

All geometry parameters are assumed not to change within an element(allowing a changing geometry within an element leads to complications, e.g. a non-constant width $b$ may lead to a fall (i.e. a transition from subcritical flow to supercritical flow) within one and the same element. In CalculiX. a changing width can be treated in a discontinuous way by using the Contraction element). Consequently:

$$\frac{\partial A}{\partial s} = 0.$$ (630)

and one obtains the Bresse equation in the form (for White-Colebrook):

$$\frac{dh}{ds}=\frac{S_0 - \frac{f}{8g} \frac{Q^2 P}{A^3}}{\sqrt{1 - S_0^2} - \frac{Q^2 B}{g A^3}}.$$ (631)

Recall that in the above formula $B$, $P$ and $A$ are a function of the depth $h$. The numerator has for positive $S_0$ exactly one root, which is called the normal depth. For this depth there is no change in $h$ along the channel. For zero or negative $S_0$ there is no root. The denominator has always exactly one root, called the critical depth. For this depth the slope is infinite. It is very weakly dependent on $S_0$. Notice that both the normal depth (if defined) and the critical depth are monotonically increasing functions of the volumetric fluid flow.

Let us for the time being assume that $S_0$ is positive. For $h$ close to zero both the denominator and numerator are negative, so the slope of $dh/ds$ is positive. For high enough values of $h$ both are positive, which also leads to a positive slope for $dh/ds$. Only for values of $h$ in between the normal and critical depth the slope $dh/ds$ is negative. For low values of $S_0$ the normal depth exceeds the critical depth and the corresopnding channel slope (slope of the bottom) is called weak. The corresponding water curves are denoted by A1, A2 and A3 depending on whether the curve is above the normal depth, in between normal depth and critical depth or below the critical depth, respectively. For high values of $S_0$ the critical depth exceeds the normal depth and the corresponding channel slope is called strong. The corresponding water curves are denoted by B1, B2 and B3 depending on whether the curve is above the critical depth, in between the critical depth and the normal depth or below the normal depth, respectively. Water curves below the critical depth are governed by upstream boundary conditions and are called frontwater curves. Water curves above the critical depth are governed by downstream boundary conditions and are called backwater curves.

Channel flow can be supercritical or subcritical. For supercritical flow the velocity exceeds the propagation speed $c$ of a wave, which satisfies $c=\sqrt{g h}$. Defining the Froude number by $Fr=U/c$, where U is the velocity of the fluid, supercritical flow corresponds to $Fr>1$. Supercritical flow is controlled by upstream boundary conditions. If the flow is subcritical ($Fr<1$) it is controlled by downstream boundary conditions. In a subcritical flow disturbances propagate upstream and downstream, in a supercritical flow they propagation downstream only. The critical depth corresponds to $U=c$. Indeed, taking a rectangular cross section the denominator of the Bresse equation is zero if

$$U^2=gh \sqrt{1-S_0^2}.$$ (632)

For frontwater curves $h$ is less than the critical depth, consequently the velocity must exceed $c$ (conservation of mass) and is supercritical. For backwater curves $h$ exceeds the critical depth and the velocity is less than $c$, the flow is subcritical.

A transition from supercritical to subcritical flow is called a hydraulic jump, a transition from subcritical to supercritical flow is a fall. At a jump the following equation is satisfied [17]:

$$\dot{m}^2 + \rho^2 g \sqrt{1 - S_0^2} A_1^2 {y_G}_1 = \dot{m}^2 + \rho^2 g \sqrt{1 - S_0^2} A_2^2 {y_G}_2,$$ (633)

Figure 148: Conservation of momentum at a jump

where $A_1, A_2$ are the cross sections before and after the jump, ${y_G}_1$ and ${y_G}_2$ is the distance orthogonal to the channel floor between the fluid surface and the center of gravity of section $A_1$ and $A_2$, respectively, $\rho$ is the fluid density and $\dot{m}$ is the mass flow. This relationship can be obtained by applying the conservation of momentum principle to a mass of infinitesimal width at the jump. The conservation of momentum dictates that the time rate of change of the momentum must equal all external forces. In Figure 148 a mass of width $ds$ is shown at time $t$ crossing a jump. At time $t+dt$ this mass is moved to the right (width $ds'$). The change in momentum in s-direction amounts to

$$\lim_{dt \rightarrow 0} [(\rho A_2 U_2 dt) U_2 - (\rho A_1 U_1 dt) U_1]/dt = \rho Q (U_2 - U_1).$$ (634)

The forces are the hydrostatic forces on the right and left side of the mass:

$$\rho g ({y_G}_1 \sqrt{1-S_0^2}) A_1 - \rho g ({y_G}_2 \sqrt{1-S_0^2}) A_2.$$ (635)

All other forces such as gravity and wall friction disappear for $ds \rightarrow 0$. Equating both terms yields the jump equation. Notice that this relationship cannot be obtained by using the Bresse equation, since $h$ is discontinuous at the jump. The discrete form of the Bernoulli equation (615) cannot be used either, since it is obtained by integrating the differential form and $dp/ds = \rho g dh/ds$ is discontinuous. However, one can write down Equation (615) pro forma and deduce the head loss in a jump by formally substituting the jump equation. One obtains:

$$\Delta F =\frac{(h_2-h_1)^3}{4 h_1 h_2}.$$ (636)

Since the head loss must be positive, this also proves that a fall cannot occur in a prismatic channel (i.e. a channel with constant cross section). Therefore, a fall can only occur at discontinuities in the channel geometry, e.g. at a discontinuous increase of the channel floor slope $S_0$.

This approach opens up an alternative to using the conservation of momentum principle at discontinuities: if one knows the head loss (e.g. by performing experiments) one can apply the discrete form of the Bernoulli equation in the form:

$$z_1 + h_1 \cos \varphi + \frac{\dot{m}^2}{2 \rho^2 A_1^2 g} = z_2 + h_2 \cos \varphi + \frac{\dot{m}^2}{2 \rho^2 A_2^2 g} + \Delta F.$$ (637)

Figure 149: Allowable volumetric flow for a given specific energy

Defining the specific energy $E$ by:

$$E := h \cos \varphi + \frac{Q^2}{2 A^2 g},$$ (638)

one can write the above equation as $z_1 + E_1 = z_2 + E_2 + \Delta F$, from which it is clear that the total head $z+E$ can never increase in the direction of the flow, however, the specific energy can. From the definition of the specific energy one can derive the dependence of $Q$ on $h$, as shown in Figure 149:

$$Q = A \sqrt{2 g (E - h \cos \varphi)}.$$ (639)

To determine the maximum allowable volumetric flow for a given value of $E$ one has to set the first derivative of the above equation to zero, resulting in (recall that $\partial A / \partial h = B$):

$$2 B (E- h \cos \varphi) = A \cos \varphi.$$ (640)

Substituting this expression in Equation (639) yields:

$$Q_{max} = A \sqrt{g \frac{A \cos \varphi}{B} },$$ (641)

or

$$\frac{Q_{max}^2}{g A^3 B} = \cos \varphi.$$ (642)

This corresponds to the denominator of the Bresse equation, i.e. the depth for which the volumetric flow is maximum is the critical depth. This is illustrated in Figure 149. Curves corresponding to lower values of $E$ also go through the origin and are completely contained in the curve shown. These curves cannot intersect, since this would mean that the intersection point corresponds to different energy values.

For a volumetric flow lower than the maximal one ($Q_1$ in the figure), two depths are feasible: a subcritical one and a supercritical one. The transition from a supercritical one to a subcritical corresponds to a jump. At the location of the jump $z_1=z_2$, however, $\Delta F \ne 0$, so $E_2 < E_1$ and the subcritical height will be slightly lower than in the figure. For geometric discontinuities for which the head loss is known (e.g. for a contraction or an enlargement) the above reasoning can be used to obtain the fluid depth downstream of the discontinuity based on the specific energy upstream (or vice versa).

Available boundary conditions for channels are the sluice gate and the weir (upstream conditions) and the infinite reservoir (downstream condition). They are described in Section 6.6. Discontinuous changes within a channel can be described using the contraction, enlargement and step elements.

The elements used in CalculiX for one-dimensional channel networks are regular network elements, in which the unknowns are the fluid depth (in z-direction, i.e. not orthogonal to the channel floor) and the temperature at the end nodes and the mass flow in the middle nodes. The equations at our disposal are the Bresse equation in the middle nodes (conservation of momentum), and the mass and energy conservation (Equations 609 and 618, respectively) at the end nodes.

For channel elements the energy equation is used in its original form:

$$\sum_{j \in \text{in}} c_p (T_j) T_j \dot{m}_{ij} - c_p (T_i) T_i... ...\text{out}} \dot{m}_{ij} + \overline{h}(T_i, T) (T - T_i) + m_i h_i^{\theta}=0,$$ (643)

in which $c_p$ for the incoming flow is determined at the correct upstream temperature $T_j$. $T$ is the temperature of the wall. The temperatures in the network are solved for as soon as the mass flow and fluid depth have been determined in the complete network. The above equation is applied on an element by element basis starting at the upstream In/Out elements and going in downstream direction. It can be reformulated as:

$$T_i = \frac{\sum_{j \in \text{in}} c_p (T_j) T_j \dot{m}_{ij}+ \o... ...theta}}{c_p (T_i) \sum_{j \in \text{out}} \dot{m}_{ij} + \overline{h}(T_i, T)}.$$ (644)

This is a slightly nonlinear equation in $T_i$ which is solved node by node in an iterative way. Convection with the wall can be defined using a *FILM card (forced convection), for the heat source *CFLUX is to be used on degree of freedom 0 (or, equivalently, 11).

Radiation to the environment can be included by modifying Equation (643) into:

$$\sum_{j \in \text{in}} c_p (T_j) T_j \dot{m}_{ij} - c_p (T_i) T_i \sum_{j \in \text{out}} \dot{m}_{ij} + \overline{h}(T_i, T) (T - T_i)$$

$$+ \epsilon(T_i) \sigma A (T^4-T_i^4) + m_i h_i^{\theta}=0,$$ (645)

where $A$ is the representative radiating surface for $T_i$, $\epsilon$ is the emissivity, $\sigma$ is the Stefan-Boltzmann constant and all temperatures have to be on an absolute scale. The radiating surface can be modeled using shell elements or solid elements covering the fluid surface. The above equation is a quartic equation, which can be solved using the conventional Newton-Raphson technique. Due to the independent processing of the energy equation after having solved the mass and momentum equation it is assumed that a change in temperature does not significantly influence the flow conditions.

Output variables are the mass flow (key MF on the *NODE PRINT or *NODE FILE card), the fluid depth (key PN -- network pressure -- on the *NODE PRINT card and DEPT on the *NODE FILE card) and the total temperature (key NT on the *NODE PRINT card and TT on the *NODE FILE card). These are the primary variables in the network. Internally, in network nodes, components one to three of the structural displacement field are used for the mass flow, the fluid depth and the critical depth, respectively. So their output can also be obtained by requesting U on the *NODE PRINT card. This is the only way to get the critical depth in the .dat file. In the .frd file the critical depth can be obtained by selecting HCRI on the *NODE FILE card. Notice that for liquids the total temperature virtually coincides with the static temperature (cf. previous section; recall that the wave speed in a channel with water depth 1 m is $\sqrt{10}$ m/s). If a jump occurs in the network, this is reported on the screen listing the element in which the jump takes place and its relative location within the element.