From 6366a2e1571c97c81972f829517b34091e5ee230 Mon Sep 17 00:00:00 2001 From: Poulard Christine <christine.poulard@irstea.fr> Date: Thu, 8 Jul 2021 16:29:29 +0000 Subject: [PATCH] Upload New File : better layout ; english translations --- subcritical_flow_channel_into_reservoir.py | 249 +++++++++++++++++++++ 1 file changed, 249 insertions(+) create mode 100644 subcritical_flow_channel_into_reservoir.py diff --git a/subcritical_flow_channel_into_reservoir.py b/subcritical_flow_channel_into_reservoir.py new file mode 100644 index 0000000..8373658 --- /dev/null +++ b/subcritical_flow_channel_into_reservoir.py @@ -0,0 +1,249 @@ +###################################################### +# Céline Berni, Christine Poulard, INRAE +# v1.02 juillet 2021 +# Démo pour le cours HySL2 Sciences de l'Eau ENTPE 1èreA +# +# /ENG/ DEMO : subcritical flow in a rectangular channel flowing into a reservoir +# +###################################################### + +import numpy as np # Numpy for array computing + +import matplotlib.pyplot as plt +from matplotlib.widgets import Slider # widgets du module matplotlib ! +from matplotlib.patches import Rectangle +from matplotlib.collections import LineCollection +import matplotlib.gridspec as gridspec + +## Calcul d'une ligne d'eau en fluvial pour un canal rectangulaire de grande largeur +# calcul d'aval en amont +# /ENG/ subcritical flow for a wide rectangular channel : computation from down- to up-stream + +## Données d'entrée /ENG/ Input data +# Q = débit /ENG/ Q = discharge +Q = 50 # m3/s/ENG/ +# Coefficient de Strickler // Strickler = 1 / Manning +K_ini = 60 # m**(1/3)/s +# pente // longitudinal slope +I = 0.0012 +# largeur /ENG/ width +b = 20 # m +# gravité /ENG/ gravity constant +g = 9.81 # m/s² + +## point de départ /ENG/ downstream point = starting point for computations (abscissis, water depth) +x0 = 0 # m +h0 = 1.5 # m +## pas de calcul /ENG/ computation step = delta water depth +delta_h_ini = 0.025 # m + +# FONCTIONS (pour être éventuellement réutilisées dans un autre code) /ENG/ FUNCTIONS + +## calcul de la hauteur normale /ENG/ normal depth +def hauteur_normale(Q, K, b, I): + hn = (Q / (K * b * np.sqrt(I))) ** (3 / 5) + return hn + + +## calcul de la hauteur critique /ENG/ critical depth +def hauteur_critique(Q, b): + hc = (Q ** 2 / (g * b ** 2)) ** (1 / 3) + return hc + + +def calcul_x_de_proche_en_proche(Q, K, b, I, h0, delta_h): + """ + step-by-step estimation of h(x), by searching for x successively for water depth h from h0 by steps of delta_h + """ + hn = hauteur_normale(Q, K, b, I) + hc = hauteur_critique(Q, b) + print(f"hauteur normale : {hn:.2f}") + print(f"hauteur critique : {hc:.2f}") + liste_h = None + liste_x = None + if hn < hc: + print('le régime uniforme doit être fluvial : calcul impossible // impossible : subcritical regime only.') + elif hn > h0: + print('la hauteur normale doit être inférieure à la hauteur imposée par la réservoir : calcul impossible.') + print('// impossible : downstream the channel, normal depth can not be higher than reservoir depth ') + else: + ## boucle de calcul pour calculer la hauteur d'eau par itération + # initialisation + hi = h0 + xi = x0 + + liste_h = [h0] + liste_x = [x0] + + while hi > hn: + dhi = min(delta_h, hi - hn) # pour ajuster le pas de calcul en dernière itération + hip1 = hi - dhi # hauteur d'eau au pas suivant i+1 + hm = hi - dhi / 2 # hauteur d'eau moyenne sur le segment i i+1 considéré + dhdxi = I * (1 - (hn / hm) ** (10 / 3)) / ( + 1 - (hc / hm) ** 3) # pente de la ligne d'eau (hyp grande largeur) + xip1 = xi - dhi / dhdxi # abscisse du point de hauteur d'eau hip1 + + # itération iplus1 + xi = xip1 + hi = hip1 + + print(f"itération : xi={xi:.3f} et hi={hi:.3f}") + + # stockage des valeurs de h et x dans un tableau + liste_h.append(hi) + liste_x.append(xi) + + return liste_h, liste_x + + +def update_from_sliders(val): + global plot_fb, ligne_normale + + # on prend comme valeur de K la valeur donnée par le slider, mais en convertissant en entier + K_slider = int(slider_K.val) + dh_slider = slider_dh.val + liste_h, liste_x = calcul_x_de_proche_en_proche(Q, K_slider, b, I, h0, dh_slider) + + try: + plot_fb.remove() + ligne_normale.remove() + except: + pass + hn = hauteur_normale(Q, K_slider, b, I) + xy_normale = (min_de_x, - I * min_de_x + hn), (max_de_x, - I * max_de_x + hn) + lcol_normale = LineCollection([xy_normale], color=["red"], lw=2, ls="--", label=f"Zf+ hn (pour K={K_slider})") + ligne_normale = ax_profile.add_collection(lcol_normale) + + if liste_h is not None and liste_x is not None: + # on a pu calculer + x = np.asarray(liste_x) + h = np.asarray(liste_h) + # altitude du fond (on a besoin de tous les points pour le tracé) + zf = - I * x + # altitude de la surface libre + zsl = zf + h + # charge + H = zf + h + Q ** 2 / (2 * g * b ** 2 * h ** 2) + + plot_zsl.set_data(x, zsl) + plot_H.set_data(x, H) + plot_fb = ax_profile.fill_between(x, y1=zf, y2=zsl, color="lightblue", alpha=0.75) + + fig.suptitle(f'profil en long // longitudinal profile, K={val:.2f}') + else: + # on n'a pas pu calculer + plot_zsl.set_data(x0, h0) + plot_H.set_data(x0, h0) + fig.suptitle(f'profil en long, K={val:.2f} ; calcul impossible') + + if hn < h_critique: + # print('le régime uniforme doit être fluvial : calcul impossible.') + plot_fb = ax_profile.fill_between([-1000, 0], y1=0, y2=h0, color="pink", alpha=0.75) + + elif hn > h0: + # print('la hauteur normale doit être inférieure à la hauteur imposée par la réservoir : calcul impossible.') + + plot_fb = ax_profile.fill_between([-1400, 0], y1=[1400 * I, 0], y2=[1400 * I + h0, h0], color="pink", alpha=0.75) + + # infos sur les sliders + + text_K.set_text(f'K={K_slider}') + text_dh.set_text(f'delta_h={dh_slider:.3f}') + + fig.canvas.draw_idle() + + +# CORPS DU PROGRAMME + +# initialisations + +liste_h, liste_x = calcul_x_de_proche_en_proche(Q, K_ini, b, I, h0, delta_h_ini) + +## pour la visualisation : on transforme les listes en vecteur numpy, pour ensuite travailler directement en VECTEURS +x = np.asarray(liste_x) +h = np.asarray(liste_h) +# vecteur des altitudes du fond +zf = - I * x # définition avec un point par point de x alors que 2 points suffisent pour définir la droite +min_de_x = min(x) +max_de_x = max(x) +xy_zfond = (min_de_x, - I * min_de_x) , (max_de_x, - I *max_de_x) +lcol_zfond = LineCollection([xy_zfond], color=["sienna"], lw=2, label='fond du canal \n Z_channel bottom') +# vecteur des altitudes de la surface libre +zsl = zf + h + +# vecteur des charges +H = zf + h + Q ** 2 / (2 * g * b ** 2 * h ** 2) + +h_critique = hauteur_critique(Q, b) +# Takes list of lines, where each line is a sequence of coordinates +xy_critique =(min_de_x, - I * min_de_x + h_critique), (max_de_x, - I * max_de_x+ h_critique) +lcol_critique = LineCollection([xy_critique], color=["orange"], lw=2, ls=":", label="Z_fond + hc (indépendant de K) \n Z_bottom + hc") + +h_normale = hauteur_normale(Q, K_ini, b, I) +xy_normale = (min_de_x, - I * min_de_x + h_normale), (max_de_x, - I * max_de_x+ h_normale) +lcol_normale = LineCollection([xy_normale], color=["red"], lw=2, ls="--", label = f"Zf+ hn (pour K={K_ini}) \n Z_bottom + hn(K={K_ini})") + + +## tracé du profil en long // PLOT + +# Grille // Layout (subplots) +gs = gridspec.GridSpec(5, 2, width_ratios = [1, 2], height_ratios = [1,10, 2, 1, 1], + left = 0, right = 1, bottom = 0.1, top = 1, hspace = 0, wspace = 0) + +fig = plt.figure() +fig.canvas.manager.set_window_title("Backwater, demo") + +# // profile plot ; items = free surface, head, bottom line, normal and critical depths +# on ajoute les courbes si possible dans l'ordre haut vers bas (cas de l'exercice) +ax_profile = plt.subplot(gs[1, 1], facecolor = 'white') # main plot : free surface profile +ax_profile.plot(x0, h0, 'v', label='cote aval imposée par le réservoir \n water level imposed by reservoir') +plot_H, = ax_profile.plot(x, H, '--', marker='.', label='ligne de charge \n head', color='green') + +plot_zsl, = ax_profile.plot(x, zsl, marker='.', label='surface libre \n free surface', color='blue') +ligne_normale = ax_profile.add_collection(lcol_normale) +# lignes fixes // fixed lines (bottom ad critical depth) +ax_profile.add_collection(lcol_critique) +ax_profile.add_collection(lcol_zfond) +plot_fb = ax_profile.fill_between(x, y1=zf, y2=zsl, color="lightblue", alpha=0.75) +rect = Rectangle((x0, -I * x0), 100, h0, color='darkblue') +ax_profile.add_patch(rect) +ax_profile.annotate("réservoir", (x0 + 10, h0 * 1.1), color="darkblue") +ax_profile.set_xlabel('x (m)') +ax_profile.set_ylabel('z (m)') +legend_profile = ax_profile.legend() +legend_profile.set_visible(False) +handles, labels = ax_profile.get_legend_handles_labels() + +ax_legend = plt.subplot(gs[1, 0], facecolor='lightcyan') # separate space devoted to legend +ax_legend.legend(handles, labels, ncol=1, prop={"size":8}, loc='lower left', bbox_to_anchor=(0.1, 0.1, 0.8, 0.8)) + +# K +ax_slider_K = plt.subplot(gs[3, 1], facecolor='lemonchiffon') # space for K slider (Strickler = 1 / Manning) +slider_K = Slider( + ax_slider_K, "Rugosité de Strickler \n roughness coefficient (1/Manning)", valmin=0.01, valmax=100, valfmt='%0.0f', valinit=K_ini, color="sienna") +slider_K.on_changed(update_from_sliders) + +# dh +ax_slider_dh = plt.subplot(gs[4, 1], facecolor = 'lemonchiffon') # space for delta_h slider +slider_dh = Slider( + ax_slider_dh, "Pas d'espace // discretisation step", valmin=0.001, valmax=0.5, valfmt='%0.3f', valinit=delta_h_ini, color="silver") +slider_dh.on_changed(update_from_sliders) + +text_K = ax_slider_K.text(0.5, 0.5, f'K={K_ini}', transform=ax_slider_K.transAxes, + fontsize=10, fontweight='bold', va='center', ha='center') +text_dh = ax_slider_dh.text(0.5, 0.5, f'delta_h={delta_h_ini:.3f}', transform=ax_slider_dh.transAxes, + fontsize=10, fontweight='bold', va='center', ha='center') + +for ax_s in [ax_legend, ax_slider_K, ax_slider_dh]: + ax_s.xaxis.set_visible(False) + ax_s.yaxis.set_visible(False) + +for pos in ['right', 'top', 'bottom', 'left']: + ax_legend.spines[pos].set_visible(False) + + +fig.suptitle(f'profil en long // longitudinal profile') + + + +plt.show() -- GitLab