Hydraulische Durchbildung einer Staustufe

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

BHQ1 = 65.0    # m³/s  Bemessungshochwasser 1 (Regelfall)
BHQ2 = 85.0    # m³/s  Bemessungshochwasser 2 (Kontrollfall)

# h_ue_zul: zulässige Überfallhöhe über Wehrkrone bei BHQ1 (= z_H - z_Kr, relativ)
# Nicht zu verwechseln mit z_H (absolutem Oberwasserstand in m ü. NHN, vgl. Notation)
h_ue_zul = 1.10  # m
f        = 0.50  # m  Freibord (DIN 19700: f ≥ 0,5 m bei BHQ1)

mu      = 0.64  # Abflussbeiwert (bewegliches Wehr, scharfkantig)
c       = 1.0   # Abminderungsbeiwert (vollkommener Überfall, Anfangswert)
xi_pf   = 0.10  # Einschnürungsbeiwert (abgerundete Pfeilerköpfe)

b_F  = 15.0              # m  gewählte Wehrfeldbreite (technisch bedingt)
b_pf = b_F * 0.225       # m  Pfeilerbreite (Drucksegmentverschluss)

g = 9.81  # m/s²

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

c_pol = (2/3) * mu * math.sqrt(2 * g)  # kombinierter Poleni-Koeffizient ohne b und h_ue

def b_eff_von_Q(Q, h_ue, c_abd=1.0):
    """Erforderliche effektive Wehrbreite für gegebenen Durchfluss."""
    return Q / (c_abd * c_pol * h_ue**1.5)

def h_ue_von_Q(Q, b_eff, c_abd=1.0):
    """Überfallhöhe bei gegebenem Durchfluss und effektiver Breite."""
    return (Q / (c_abd * c_pol * b_eff))**(2/3)

def b_eff_n_felder(n_aktiv, h_ue):
    """Effektive Breite für n_aktiv nebeneinanderliegende Felder (Einschnürung berücksichtigt)."""
    return n_aktiv * (b_F - 2 * xi_pf * h_ue)

# Gesamte b_eff für (n-1)-Fall muss b_eff_erf >= b_eff_von_Q(BHQ1, h_ue_zul) erfüllen
b_eff_erf = b_eff_von_Q(BHQ1, h_ue_zul, c)
b_F_eff   = b_F - 2 * xi_pf * h_ue_zul     # wirksame Feldbreite je Feld

n_minus1  = math.ceil(b_eff_erf / b_F_eff)  # benötigte Anzahl (n-1)-Felder
n         = n_minus1 + 1                     # Gesamtanzahl Wehrfelder

print(f"Kombinierter Poleni-Koeffizient c_pol = {c_pol:.4f}")
print(f"Erforderliche b_eff (n-1-Fall) = {b_eff_erf:.2f} m")
print(f"Wirksame Feldbreite b_F_eff    = {b_F_eff:.2f} m")
print(f"Benötigte (n-1)-Felder         = {n_minus1}")
print(f"Gewählte Wehrfeldanzahl n      = {n}")

n_pfeiler  = n - 1               # Zwischenpfeiler
b_licht    = n * b_F             # lichte Gesamtbreite (nur Felder)
b_pfeiler  = n_pfeiler * b_pf    # Summe Pfeilerbreiten
b_gesamt   = b_licht + b_pfeiler # Gesamtbreite inkl. Pfeiler

print(f"Anzahl Zwischenpfeiler : {n_pfeiler}")
print(f"Lichte Gesamtbreite    : {b_licht:.1f} m")
print(f"Summe Pfeilerbreiten   : {b_pfeiler:.2f} m")
print(f"Bauwerksbreite gesamt  : {b_gesamt:.2f} m")

faelle = [
    ("BHQ1, alle n Felder",   BHQ1, n),
    ("BHQ1, (n-1) Felder",    BHQ1, n - 1),
    ("BHQ2, alle n Felder",   BHQ2, n),
]

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

print(f"\n  h_ü,zul = {h_ue_zul:.2f} m  (max. zulässige Überfallhöhe bei BHQ1)")

# h_ue bei BHQ1, (n-1)-Felder (maßgebender Fall der (n-1)-Regel)
h_ue_n1 = h_ue_zul
for _ in range(20):
    beff_n1 = b_eff_n_felder(n - 1, h_ue_n1)
    h_ue_n1 = h_ue_von_Q(BHQ1, beff_n1, c)

stau_marge  = h_ue_zul - h_ue_n1   # > 0: zulässiger Aufstau eingehalten
kronenhoehe = h_ue_zul + f

stau_kriterium = 'OK' if stau_marge >= 0 else 'Nicht erfüllt'

print(f"h_ü bei BHQ1, (n-1)-Felder        : {h_ue_n1:.3f} m")
print(f"Marge bis h_ü,zul (Stau-Reserve)  : {stau_marge:.3f} m")
print(f"Kronenhöhe (= h_ü,zul + f)        : {kronenhoehe:.2f} m über Wehrschwelle")
print(f"Stauziel-Nachweis h_ü <= h_ü,zul  : {stau_kriterium}")
print(f"Freibordanforderung f >= 0,5 m    : {'OK' if f >= 0.5 else 'Nicht erfüllt'}")

w_wehr = 2.5   # m       Wehrhöhe über UW-Sohle
k_st   = 30.0  # m^1/3/s Strickler-Beiwert UW-Strecke (naturnahe Sohle)
I_E    = 0.001 # --      Energieliniengefälle UW-Strecke (= Sohlgefälle bei Normalabfluss)
b_u    = 40.0  # m       Breite UW-Strecke (Rechteckprofil)

# Einheitsdurchfluss für Tosbecken: konservativer (n-1)-Felder-Fall (DIN 19700)
n_tos    = n - 1     # aktive Felder im maßgebenden Lastfall
h_ue_tos = h_ue_zul  # Startwert für Fixpunkt-Iteration
for _ in range(20):
    beff_tos = b_eff_n_felder(n_tos, h_ue_tos)
    h_ue_tos = h_ue_von_Q(BHQ1, beff_tos, c)

q = BHQ1 / beff_tos   # m²/s – Einheitsdurchfluss im (n-1)-Felder-Fall

# h_ü bei BHQ1 mit allen n Feldern (für die Zusammenfassung)
h_ue_n = h_ue_zul
for _ in range(20):
    beff_n = b_eff_n_felder(n, h_ue_n)
    h_ue_n = h_ue_von_Q(BHQ1, beff_n, c)

print(f"Maßgebender Lastfall: {n_tos} aktive Felder, b_eff = {beff_tos:.2f} m, h_ü = {h_ue_tos:.3f} m")
print(f"Einheitsdurchfluss 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 sucht h in [0.01, 20] m mit 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"Unterwassertiefe h_u = {h_u:.3f} m")
print(f"Probe: Q_Manning = {Q_manning(h_u, b_u, k_st, I_E):.2f} m³/s  (Soll: {BHQ1} m³/s)")

def tosbecken_check(e):
    """
    Führt Schritte A-G des Tosbecken-Workflows für gegebene Eintiefung e durch.
    Gibt (Fr1, epsilon, h1, h2, ok_Fr, ok_eps) zurück.
    """
    # Schritt B: h1 aus Bernoulli (kleine, schießende Wurzel);
    # Energiehöhe konsistent mit dem Lastfall des Einheitsdurchflusses q
    H_ges = h_ue_tos + w_wehr + e
    # Suche h1 < h_grenz (kritische Tiefe) = (q**2/g)^(1/3)
    h_krit = (q**2 / g)**(1/3)
    h1 = brentq(lambda h: h + q**2 / (2 * g * h**2) - H_ges, 1e-4, h_krit)
    v1 = q / h1

    # Schritt C: Froude-Zahl
    Fr1 = v1 / math.sqrt(g * h1)

    # Schritt D: Fr1-Check
    ok_Fr = 4.5 <= Fr1 < 9.0

    # Schritt F: konjugierte Tiefe h2 (Belanger)
    h2 = (h1 / 2) * (math.sqrt(1 + 8 * Fr1**2) - 1)

    # Schritt G: Einstaugrad
    eps = (h_u + e) / h2
    ok_eps = 1.05 <= eps <= 1.15

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


# k-Faktor für Tosbeckenlänge (vgl. Peterka 1984 / USBR)
def k_faktor(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 Startwert (Schritt A)
e_min    = 0.0   # m Suchbereich untere Grenze
e_max    = 3.0   # m Suchbereich obere Grenze
iterationen = []

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 = tosbecken_check(e_test)
    iterationen.append((e_test, Fr1, eps, h1, h2, ok_Fr, ok_eps))

    fr_sym  = "OK" if ok_Fr  else "Nicht erfüllt"
    eps_sym = "OK" if ok_eps else "Nicht erfüllt"
    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 >>> Kriterien erfüllt: Iteration beendet.")
        break

    # Anpassungsstrategie (Bisektionslogik, vgl. Vorlesungsfolie):
    # größeres e -> größere Energiehöhe -> kleineres h1 -> größeres Fr1
    if not ok_Fr:
        # Fr1 < 4.5 >>> e zu klein (e erhöhen); Fr1 >= 9.0 >>> e zu groß (e verringern)
        if Fr1 < 4.5:
            e_min = e_test
        else:  # Fr1 >= 9.0
            e_max = e_test
    else:
        # ok_Fr = True, aber ok_eps = False
        if eps > 1.15:   # zu viel Einstau >>> e verringern
            e_max = e_test
        else:            # eps < 1.05 >>> zu wenig Einstau >>> e erhöhen
            e_min = e_test

    e_test = 0.5 * (e_min + e_max)  # Bisektion

# Endwerte speichern
e_opt = e_test
Fr1_opt, eps_opt, h1_opt, h2_opt, ok_Fr_opt, ok_eps_opt = tosbecken_check(e_opt)
if not (ok_Fr_opt and ok_eps_opt):
    print("\n >>> Keine Lösung im Suchbereich [e_min, e_max]: Kriterien nicht gleichzeitig erfüllbar.")
    print("     Abhilfe: Schikanen (Blöcke, Zähne, Endschwelle) oder Randbedingungen prüfen (vgl. Vorlesung).")

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

print(f"Optimierte Eintiefung   e   = {e_opt:.3f} m")
print(f"Schießende Tiefe        h1  = {h1_opt:.4f} m")
print(f"Froude-Zahl             Fr1 = {Fr1_opt:.2f}")
print(f"Konjugierte Tiefe       h2  = {h2_opt:.3f} m")
print(f"Unterwassertiefe        h_u = {h_u:.3f} m")
print(f"Einstaugrad             ε   = {eps_opt:.3f}  (Ziel: 1,05-1,15)")
print()
print(f"k-Faktor (Fr1={Fr1_opt:.1f})  k   = {k:.2f}")
print(f"Tosbeckenlänge          l_T = {l_T:.2f} m")
print(f"Kolkschutzlänge         l_K = {l_K:.2f} m")

summary = {
    "Hydraulik_Wehr": {
        "BHQ1 (m3/s)":                       BHQ1,
        "BHQ2 (m3/s)":                       BHQ2,
        "h_ue_zul (m)":                      h_ue_zul,
        "Freibord (m)":                       f,
        "n (Wehrfelder)":                     n,
        "bF (m)":                            b_F,
        "Lichte Breite (m)":                  round(b_licht, 2),
        "Gesamtbreite (m)":                   round(b_gesamt, 2),
        "hü bei BHQ1 mit n Feldern (m)":      round(h_ue_n, 3),
        "hü bei BHQ1 mit (n-1) Feldern (m)":  round(h_ue_n1, 3),
        "Stau-Marge (n-1)-Fall (m)":          round(stau_marge, 3),
        "Nachweis h_ü <= h_ü,zul":            stau_kriterium,
        "Kronenhöhe (m)":                     round(kronenhoehe, 2),
    },
    "Tosbecken": {
        "Wehrhöhe w (m)":                w_wehr,
        "Abfluss je Meter Wehrbreite (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),
        "Energieumwandlung im Tosbecken": 'OK' if (ok_Fr_opt and ok_eps_opt) else 'Nicht erfüllt',
    },
}
pprint(summary)