Hydraulic Design of a Barrage

import math
from scipy.optimize import brentq
from pprint import pprint

BHQ1 = 65.0    # m³/s  design flood 1 (standard case)
BHQ2 = 85.0    # m³/s  design flood 2 (control case)

# h_ue_zul: permissible overflow head above the weir crest at BHQ1 (= z_H - z_Kr, relative)
# Not to be confused with z_H (absolute headwater level in m a.s.l., see notation)
h_ue_zul = 1.10  # m
f        = 0.50  # m  freeboard (DIN 19700: f ≥ 0.5 m at BHQ1)

mu      = 0.64  # discharge coefficient (movable weir, sharp-crested)
c       = 1.0   # reduction coefficient (free overflow, initial value)
xi_pf   = 0.10  # pier contraction coefficient (rounded pier noses)

b_F  = 15.0              # m  chosen weir bay width (technically constrained)
b_pf = b_F * 0.225       # m  pier width (radial gate)

g = 9.81  # m/s²

print(f"Pier width b_pf = {b_pf:.2f} m")

c_pol = (2/3) * mu * math.sqrt(2 * g)  # combined Poleni coefficient without b and h_ue

def b_eff_from_Q(Q, h_ue, c_red=1.0):
    """Required effective weir width for a given discharge."""
    return Q / (c_red * c_pol * h_ue**1.5)

def h_ue_from_Q(Q, b_eff, c_red=1.0):
    """Overflow head for a given discharge and effective width."""
    return (Q / (c_red * c_pol * b_eff))**(2/3)

def b_eff_n_bays(n_active, h_ue):
    """Effective width of n_active adjacent weir bays (accounting for pier contraction)."""
    return n_active * (b_F - 2 * xi_pf * h_ue)

# The total b_eff for the (n-1) case must satisfy b_eff_req >= b_eff_from_Q(BHQ1, h_ue_zul)
b_eff_req = b_eff_from_Q(BHQ1, h_ue_zul, c)
b_F_eff   = b_F - 2 * xi_pf * h_ue_zul     # effective width per bay

n_minus1  = math.ceil(b_eff_req / b_F_eff)  # required number of (n-1) bays
n         = n_minus1 + 1                    # total number of weir bays

print(f"Combined Poleni coefficient c_pol = {c_pol:.4f}")
print(f"Required b_eff (n-1 case)      = {b_eff_req:.2f} m")
print(f"Effective bay width b_F_eff    = {b_F_eff:.2f} m")
print(f"Required (n-1) bays            = {n_minus1}")
print(f"Chosen number of weir bays n   = {n}")

n_piers    = n - 1               # intermediate piers
b_clear    = n * b_F             # total clear width (bays only)
b_piers    = n_piers * b_pf      # sum of pier widths
b_total    = b_clear + b_piers   # total width including piers

print(f"Number of intermediate piers : {n_piers}")
print(f"Total clear width            : {b_clear:.1f} m")
print(f"Sum of pier widths           : {b_piers:.2f} m")
print(f"Total structure width        : {b_total:.2f} m")

cases = [
    ("BHQ1, all n bays",      BHQ1, n),
    ("BHQ1, (n-1) bays",      BHQ1, n - 1),
    ("BHQ2, all n bays",      BHQ2, n),
]

print(f"{'Case':<28} {'n_act':>5} {'b_eff':>8} {'h_ü':>7} {'≤ h_ü,zul?':>10}")
print("-" * 62)
for label, Q, n_act in cases:
    h_ue_iter = h_ue_zul
    for _ in range(20):
        beff = b_eff_n_bays(n_act, h_ue_iter)
        h_ue_new = h_ue_from_Q(Q, beff, c)
        if abs(h_ue_new - h_ue_iter) < 1e-6:
            break
        h_ue_iter = h_ue_new
    ok = "OK" if h_ue_iter <= h_ue_zul else "!!"
    print(f"{label:<28} {n_act:>5} {beff:>8.2f} {h_ue_iter:>7.3f} {ok:>10}")

print(f"\n  h_ü,zul = {h_ue_zul:.2f} m  (max. permissible overflow head at BHQ1)")

# h_ue at BHQ1, (n-1) bays (governing case of the (n-1) rule)
h_ue_n1 = h_ue_zul
for _ in range(20):
    beff_n1 = b_eff_n_bays(n - 1, h_ue_n1)
    h_ue_n1 = h_ue_from_Q(BHQ1, beff_n1, c)

head_margin  = h_ue_zul - h_ue_n1   # > 0: permissible backwater level respected
crest_height = h_ue_zul + f

head_criterion = 'OK' if head_margin >= 0 else 'Not fulfilled'

print(f"h_ü at BHQ1, (n-1) bays              : {h_ue_n1:.3f} m")
print(f"Margin up to h_ü,zul (head reserve)  : {head_margin:.3f} m")
print(f"Crest height (= h_ü,zul + f)         : {crest_height:.2f} m above weir sill")
print(f"Retention level check h_ü <= h_ü,zul : {head_criterion}")
print(f"Freeboard requirement f >= 0.5 m     : {'OK' if f >= 0.5 else 'Not fulfilled'}")

w_weir = 2.5   # m       weir height above the tailwater bed
k_st   = 30.0  # m^1/3/s Strickler coefficient of the tailwater reach (near-natural bed)
I_E    = 0.001 # --      energy slope of the tailwater reach (= bed slope for normal flow)
b_u    = 40.0  # m       width of the tailwater reach (rectangular cross-section)

# unit discharge for the stilling basin: conservative (n-1) bay case (DIN 19700)
n_sb     = n - 1     # active bays in the governing load case
h_ue_sb  = h_ue_zul  # initial value for the fixed-point iteration
for _ in range(20):
    beff_sb = b_eff_n_bays(n_sb, h_ue_sb)
    h_ue_sb = h_ue_from_Q(BHQ1, beff_sb, c)

q = BHQ1 / beff_sb   # m²/s - unit discharge in the (n-1) bay case

# h_ü at BHQ1 with all n bays (for the summary)
h_ue_n = h_ue_zul
for _ in range(20):
    beff_n = b_eff_n_bays(n, h_ue_n)
    h_ue_n = h_ue_from_Q(BHQ1, beff_n, c)

print(f"Governing load case: {n_sb} active bays, b_eff = {beff_sb:.2f} m, h_ü = {h_ue_sb:.3f} m")
print(f"Unit discharge q = BHQ1/b_eff     = {q:.3f} m²/s")

def Q_manning(h, b, k, I):
    A = b * h
    R = (b * h) / (b + 2 * h)
    return k * A * R**(2/3) * math.sqrt(I)

# brentq searches for h in [0.01, 20] m such that Q_manning(h) = BHQ1
h_u = brentq(lambda h: Q_manning(h, b_u, k_st, I_E) - BHQ1, 0.01, 20.0)

print(f"Tailwater depth h_u = {h_u:.3f} m")
print(f"Check: Q_Manning = {Q_manning(h_u, b_u, k_st, I_E):.2f} m³/s  (target: {BHQ1} m³/s)")

def stilling_basin_check(e):
    """
    Performs steps A-G of the stilling basin workflow for a given depression e.
    Returns (Fr1, epsilon, h1, h2, ok_Fr, ok_eps).
    """
    # Step B: h1 from Bernoulli (small, supercritical root);
    # energy head consistent with the load case of the unit discharge q
    H_ges = h_ue_sb + w_weir + e
    # search h1 < h_crit (critical depth) = (q**2/g)^(1/3)
    h_crit = (q**2 / g)**(1/3)
    h1 = brentq(lambda h: h + q**2 / (2 * g * h**2) - H_ges, 1e-4, h_crit)
    v1 = q / h1

    # Step C: Froude number
    Fr1 = v1 / math.sqrt(g * h1)

    # Step D: Fr1 check
    ok_Fr = 4.5 <= Fr1 < 9.0

    # Step F: sequent (conjugate) depth h2 (Belanger)
    h2 = (h1 / 2) * (math.sqrt(1 + 8 * Fr1**2) - 1)

    # Step G: submergence ratio
    eps = (h_u + e) / h2
    ok_eps = 1.05 <= eps <= 1.15

    return Fr1, eps, h1, h2, ok_Fr, ok_eps


# k-factor for the stilling basin length (cf. Peterka 1984 / USBR)
def k_factor(Fr1):
    if   Fr1 < 2.4:  return 4.8
    elif Fr1 < 4.0:  return 4.8 + (Fr1 - 2.4) / (4.0 - 2.4) * (5.8 - 4.8)
    elif Fr1 < 5.0:  return 5.8 + (Fr1 - 4.0) * (6.0 - 5.8)
    elif Fr1 < 6.0:  return 6.0 + (Fr1 - 5.0) * (6.13 - 6.0)
    elif Fr1 <= 11.0: return 6.13
    else:            return 6.0

e_test   = 1.0   # m initial value (step A)
e_min    = 0.0   # m lower limit of the search range
e_max    = 3.0   # m upper limit of the search range
iterations = []

print(f"{'Iter':>4}  {'e [m]':>7}  {'Fr1':>6}  {'Fr-OK':>6}  {'h1 [m]':>8}  {'h2 [m]':>8}  {'ε':>6}  {'ε-OK':>6}")
print("-" * 65)

for i in range(1, 15):
    Fr1, eps, h1, h2, ok_Fr, ok_eps = stilling_basin_check(e_test)
    iterations.append((e_test, Fr1, eps, h1, h2, ok_Fr, ok_eps))

    fr_sym  = "OK" if ok_Fr  else "Not fulfilled"
    eps_sym = "OK" if ok_eps else "Not fulfilled"
    print(f"{i:>4}  {e_test:>7.3f}  {Fr1:>6.2f}  {fr_sym:>6}  {h1:>8.4f}  {h2:>8.4f}  {eps:>6.3f}  {eps_sym:>6}")

    if ok_Fr and ok_eps:
        print("\n >>> Criteria fulfilled: iteration finished.")
        break

    # adjustment strategy (bisection logic, cf. lecture slide):
    # larger e -> larger energy head -> smaller h1 -> larger Fr1
    if not ok_Fr:
        # Fr1 < 4.5 >>> e too small (increase e); Fr1 >= 9.0 >>> e too large (decrease e)
        if Fr1 < 4.5:
            e_min = e_test
        else:  # Fr1 >= 9.0
            e_max = e_test
    else:
        # ok_Fr = True, but ok_eps = False
        if eps > 1.15:   # too much submergence >>> decrease e
            e_max = e_test
        else:            # eps < 1.05 >>> too little submergence >>> increase e
            e_min = e_test

    e_test = 0.5 * (e_min + e_max)  # bisection

# store final values
e_opt = e_test
Fr1_opt, eps_opt, h1_opt, h2_opt, ok_Fr_opt, ok_eps_opt = stilling_basin_check(e_opt)
if not (ok_Fr_opt and ok_eps_opt):
    print("\n >>> No solution in the search range [e_min, e_max]: criteria cannot be fulfilled simultaneously.")
    print("     Remedy: baffle elements (blocks, chute blocks, end sill) or check boundary conditions (cf. lecture).")

k   = k_factor(Fr1_opt)
l_T = k * h2_opt
l_K = 3.5 * l_T

print(f"Optimized depression      e   = {e_opt:.3f} m")
print(f"Supercritical depth       h1  = {h1_opt:.4f} m")
print(f"Froude number             Fr1 = {Fr1_opt:.2f}")
print(f"Sequent depth             h2  = {h2_opt:.3f} m")
print(f"Tailwater depth           h_u = {h_u:.3f} m")
print(f"Submergence ratio         ε   = {eps_opt:.3f}  (target: 1.05-1.15)")
print()
print(f"k-factor (Fr1={Fr1_opt:.1f})    k   = {k:.2f}")
print(f"Stilling basin length     l_T = {l_T:.2f} m")
print(f"Scour protection length   l_K = {l_K:.2f} m")

summary = {
    "Weir_hydraulics": {
        "BHQ1 (m3/s)":                        BHQ1,
        "BHQ2 (m3/s)":                        BHQ2,
        "h_ue_zul (m)":                       h_ue_zul,
        "Freeboard (m)":                      f,
        "n (weir bays)":                      n,
        "bF (m)":                             b_F,
        "Clear width (m)":                    round(b_clear, 2),
        "Total width (m)":                    round(b_total, 2),
        "hü at BHQ1 with n bays (m)":         round(h_ue_n, 3),
        "hü at BHQ1 with (n-1) bays (m)":     round(h_ue_n1, 3),
        "Head margin (n-1) case (m)":         round(head_margin, 3),
        "Check h_ü <= h_ü,zul":               head_criterion,
        "Crest height (m)":                   round(crest_height, 2),
    },
    "Stilling_basin": {
        "Weir height w (m)":                  w_weir,
        "Discharge per meter weir width (m²/s)": round(q, 3),
        "hu (m)":      round(h_u, 3),
        "eopt (m)":    round(e_opt, 3),
        "h1 (m)":      round(h1_opt, 4),
        "h2 (m)":      round(h2_opt, 3),
        "Fr1":         round(Fr1_opt, 2),
        "epsilon":     round(eps_opt, 3),
        "k":           round(k, 2),
        "lT (m)":      round(l_T, 2),
        "lK (m)":      round(l_K, 2),
        "Energy dissipation in the stilling basin": 'OK' if (ok_Fr_opt and ok_eps_opt) else 'Not fulfilled',
    },
}
pprint(summary)