Effectuer une analyse isocurrent des canaux en eau libre avec Python et matplotlib
introduction
Un isocourant est un écoulement dans lequel la forme de la section transversale, le gradient du lit de la rivière, la profondeur de l'eau et le débit ne changent pas dans le temps et le lieu, et dans les relations fluviales, il est souvent analysé à l'aide de la formule de vitesse d'écoulement moyenne de Manning. La formule suivante est utilisée pour l'analyse.
- $ v $: vitesse d'écoulement moyenne
- $ n $: Coefficient de rugosité
- $ R $: rayon hydraulique
- $ I $: Pente du lit de la rivière
- $ A $: Coupe transversale de l'eau courante
- $ S $: Junbe
- $ Q $: débit
Étant donné que la section transversale de l'eau courante et du côté humide sont des fonctions de la profondeur de l'eau, la méthode consistant à supposer d'abord la profondeur de l'eau $ h $ (ou l'altitude du niveau d'eau) et à calculer le débit descendant $ Q $ à cette profondeur d'eau est un processus de calcul simple. y a-t-il.
Script Shell pour l'exécution
Dans le script d'exécution, les trois programmes suivants sont exécutés.
- py_eng_uniflow.py: Corps du programme d'analyse de flux non uniforme
- py_eng_uni_figHQ.py: Dessin de la courbe niveau-débit d'eau
- py_eng_uni_figsec.py: Dessin en coupe d'analyse
Le texte du script d'exécution est le suivant.
a_py_uni.txt
python3 py_eng_uniflow.py inp_sec1.txt > out_PSHQ1.txt
python3 py_eng_uni_figHQ.py out_PSHQ1.txt fig_PSHQ1.png
python3 py_eng_uni_figsec.py inp_sec1.txt fig_sec1_PS.png n=0.04,i=0.001 79.4
Des données d'entrée
1ère ligne: coefficient de rugosité de la 1ère colonne $ n $, pente du lit de la 2ème colonne $ I $ 2e ligne et suivantes: 1ère colonne coordonnée x du lit de la rivière (de la rive gauche à la rive droite), coordonnée y (élévation)
inp_sec1.txt
0.04 0.001
-138 100
-114 90
-70 80
-45 70
-32 62
0 61.5
32 62
60 70
80 80
98 84
120 84
des données de sortie
La séquence des données de sortie est la suivante. Pour l'analyse isocinétique uniquement, il suffit de connaître la relation entre le débit et l'altitude ou la profondeur de l'eau, mais dans le but de créer des conditions d'entrée pour l'analyse d'écoulement inégal, la section transversale de l'eau qui coule, le côté humide et la largeur de la surface de l'eau sont également sorties. À propos, la largeur de la surface de l'eau est l'information nécessaire pour calculer la profondeur limite de l'eau.
- 1ère rangée: débit
- Deuxième rangée: élévation du niveau d'eau
- 3e rangée: profondeur de l'eau
- 4ème rangée: section transversale de l'eau courante
- 5e rangée: Junbe
- 6ème rangée: largeur de la surface de l'eau
out_PSHQ1.txt
Q WL h A S wsw
0.000 61.500 0.000 0.000 0.000 0.000
0.069 61.600 0.100 0.640 12.802 12.800
0.436 61.700 0.200 2.560 25.603 25.600
1.285 61.800 0.300 5.760 38.405 38.400
2.768 61.900 0.400 10.240 51.206 51.200
5.019 62.000 0.500 16.000 64.008 64.000
8.760 62.100 0.600 22.426 64.563 64.513
.......... ..........
Programme d'analyse isocurrent
Dans le programme, l'aire transversale $ A $ et la marge $ S $ sont calculées. Cela se fait en suivant le segment de ligne représentant le terrain d'entrée depuis la rive gauche (côté gauche) en référence à la figure ci-dessous, et en calculant la zone et la longueur de la ligne de terrain de la partie entourée par la ligne d'élévation (el) et la ligne de terrain spécifiées. J'essaye de le faire.
py_eng_uniflow.py
import numpy as np
from math import *
import sys
import matplotlib.pyplot as plt
def INTERSEC(x1,y1,x2,y2,el):
xx=1.0e10
if 0.0<abs(y2-y1):
a=(y2-y1)/(x2-x1)
b=(x2*y1-x1*y2)/(x2-x1)
xx=(el-b)/a
return xx
# Main routine
param=sys.argv
fnameR=param[1]
# Input data
x=[]
y=[]
fin=open(fnameR,'r')
text=fin.readline()
text=text.strip()
text=text.split()
rn=float(text[0]) # Manning roughness coefficient
gi=float(text[1]) # Gradient of river bed
while 1:
text=fin.readline()
if not text: break
text=text.strip()
text=text.split()
x=x+[float(text[0])]
y=y+[float(text[1])]
fin.close()
nn=len(x)
print('{0:>10s} {1:>10s} {2:>10s} {3:>10s} {4:>10s} {5:>10s}'.format('Q','WL','h','A','S','wsw'))
ymin=min(y)
ymax=max(y)
dy=0.1
yy=np.arange(ymin,ymax+dy,dy)
for el in yy:
S=0.0
A=0.0
xxl=0.0
xxr=0.0
for i in range(0,nn-1):
x1=x[i]
y1=y[i]
x2=x[i+1]
y2=y[i+1]
xx=INTERSEC(x1,y1,x2,y2,el)
dS=0.0
dA=0.0
if y1<el and y2<el:
dS=sqrt((x2-x1)**2+(y2-y1)**2)
dA=0.5*(2.0*el-y1-y2)*(x2-x1)
if x1<=xx and xx<=x2:
if y2<=el and el<=y1:
dS=sqrt((x2-xx)**2+(y2-el)**2)
dA=0.5*(x2-xx)*(el-y2)
xxl=xx
if y1<=el and el<=y2:
dS=sqrt((xx-x1)**2+(el-y1)**2)
dA=0.5*(xx-x1)*(el-y1)
xxr=xx
S=S+dS
A=A+dA
if 0.0<S:
R=A/S
v=1.0/rn*R**(2/3)*sqrt(gi)
Q=A*v
else:
R=0.0
v=0.0
Q=0.0
h=el-ymin
wsw=xxr-xxl
print('{0:10.3f} {1:10.3f} {2:10.3f} {3:10.3f} {4:10.3f} {5:10.3f}'.format(Q,el,h,A,S,wsw))
Programme de dessin de courbe niveau-débit d'eau
Ceci suit la procédure standard de création de graphe bidimensionnel.
py_eng_uni_figHQ.py
import numpy as np
import matplotlib.pyplot as plt
import sys
def DINP(fnameR):
data=np.loadtxt(fnameR,skiprows=1,usecols=(0,1))
QQ=data[:,0]
EL=data[:,1]
return QQ,EL
# Parameter setting
param=sys.argv
fnameR=param[1] # input filename
fnameF=param[2] # output image file name
qmin=0.0;qmax=10000.0;dq=1000
emin=60.0;emax=90.0;de=2.0
# Input data
QQ,EL=DINP(fnameR)
# Plot
fig=plt.figure()
plt.plot(QQ,EL,color='black',lw=2.0,label='PS (n=0.04, i=0.001)')
plt.xlabel('Discharge (m$^3$/s)')
plt.ylabel('Elevation (EL.m)')
plt.xlim(qmin,qmax)
plt.ylim(emin,emax)
plt.xticks(np.arange(qmin,qmax+dq,dq))
plt.yticks(np.arange(emin,emax+de,de))
plt.grid()
plt.legend(shadow=True,loc='upper left',handlelength=3)
plt.savefig(fnameF,dpi=200,, bbox_inches="tight", pad_inches=0.2)
plt.show()
plt.close()
Image de sortie
Programme de dessin de coupe calculée
Bien que ce soit une question d'apparence, la surface sous-marine est remplie de bleu clair en utilisant fill_between.
py_eng_uni_figsec.py
import numpy as np
#from math import *
import sys
import matplotlib.pyplot as plt
# Main routine
param=sys.argv
fnameR=param[1] # input filename
fnameF=param[2] # output image file name
strt=param[3] # graph title (strings)
fwl=float(param[4]) # target flood water level
# Input data
x=[]
y=[]
fin=open(fnameR,'r')
text=fin.readline()
text=text.strip()
text=text.split()
rn=float(text[0])
gi=float(text[1])
while 1:
text=fin.readline()
if not text: break
text=text.strip()
text=text.split()
x=x+[float(text[0])]
y=y+[float(text[1])]
fin.close()
xx=np.array(x)
yy=np.array(y)
# plot
xmin=-150
xmax=150.0
ymin=50.0
ymax=100
fig = plt.figure()
ax=fig.add_subplot(111)
ax.plot(xx,yy,color='black',lw=1.0,label='ground surface')
ax.fill_between(xx,yy,fwl,where=yy<=fwl,facecolor='cyan',alpha=0.3,interpolate=True)
ax.text(0,fwl,'10,000 years flood (EL.'+str(fwl)+')',fontsize=10,color='black',ha='center',va='top')
ax.set_title(strt)
ax.set_xlabel('Distance (m)')
ax.set_ylabel('Elevation (EL.m)')
ax.set_xlim(xmin,xmax)
ax.set_ylim(ymin,ymax)
aspect=1.0*(ymax-ymin)/(xmax-xmin)*(ax.get_xlim()[1] - ax.get_xlim()[0]) / (ax.get_ylim()[1] - ax.get_ylim()[0])
ax.set_aspect(aspect)
plt.savefig(fnameF,dpi=200,, bbox_inches="tight", pad_inches=0.2)
plt.show()
plt.close()
Image de sortie
c'est tout
Recommended Posts