Python水力学和水文学应用

Python(解非线性方程和线性方程)求水力学法向深度-浪涌高度速度及互连反应器中的浓度和流体分布

非线性方程

在水力学领域遇到的非线性方程的一个例子是通过长梯形通道寻找流动的法向深度 $y_n$ 。这样的流动深度出现在均匀流动区域，远离任何不均匀原因的影响，例如堰的上游。

法向深度 $y_n$ 可以通过求解以下方程获得：

$Q=\frac{1}{n} A R^{2 / 3} S_{0}^{1 / 2}$

浪涌高度和速度

非线性方程的另一种情况，也来自水利工程领域，是计算由于在其一端的流量突然增加或减少而在一维通道中发生的冲击波的高度。

图示了两种这样的情况，一种是由于下游端的流量减少，另一种是由于上游端的流量增加。两幅图都描绘了所谓的正涌浪，即行波后面的水深大于其前面的水深。第一种情况可能发生在水电运河中，其中控制进入涡轮机的水的下游闸门突然运行以减少排放。另一个可能在水力发电机组的出口通道中，由于涡轮转轮的启动，排放量突然增加。另一种类型的浪涌可能发生在带有闸门的通道中，该通道阻止了储备的水，而储备的水突然被撤回。然后产生一个负浪涌，其中波浪向上传播的水流深度小于其前面的水深。

狭窄和凸起通道部分的水流深度

线性方程

在求解联立方程时通常会出现一组包含 n 个未知变量的 n 个方程。下面介绍了这种“方程组”的一些代表性示例，这些示例取自水利工程和水文学领域，需要同时求解方程。

反应堆稳态分析系统

管网流量的稳态分布

单位水文的推导

Python求解上述方程

import numpy as np
import matplotlib.pylab as plt
b0 = 20.0
s = 2.0
s0 = 0.001
mn = 0.020
Q = 40
yinitial = 1.0
errorallow = 0.0001
def area(y):
 area=y*(b0+s*y)
 return area
def wetperi(depth):
 wetperi = b0+2*depth*np.sqrt(1+s*s)
 return wetperi
def hyrad(depth):
 hyrad = (b0+depth*s)*depth/(b0+2*depth*np.sqrt(1+s*s))
 return hyrad
def topwidth(depth):
 topwidth=b0+2*s*depth
 return topwidth
iter = 0
max_error = 1.0
yn = yinitial
while (max_error > errorallow):
 iter = iter+1
 Fyn = area(yn)**1.6667/wetperi(yn)**0.6667-mn*Q/
np.sqrt(s0)

源代码

Python(休恩法和龙格-库塔法解常微分方程)求水箱排水和溢洪道水流及水路生化需氧量和水位波动

休恩法

在数学和计算科学中，休恩方法可能是指改进的或修正的欧拉方法（即显式梯形法则），或类似的两阶段龙格-库塔方法。它以 Karl Heun 命名，是一种求解具有给定初始值的常微分方程 (ODE) 的数值过程。这两种变体都可以看作是欧拉方法对两阶段二阶龙格-库塔方法的扩展。

计算初值问题数值解的过程：

$y^{\prime}(t)=f(t, y(t)), \quad y\left(t_{0}\right)=y_{0}$

通过休恩方法，是先计算中间值 $\tilde{y}{i+1}$ ，然后在下一个积分点处的最终近似值 $y{i+1}$ 。

$\begin{aligned}&\tilde{y}{i+1}=y{i}+h f\left(t_{i}, y_{i}\right) \\&y_{i+1}=y_{i}+\frac{h}{2}\left[f\left(t_{i}, y_{i}\right)+f\left(t_{i+1}, \tilde{y}_{i+1}\right)\right]\end{aligned}$

其中 $h$ 是步长， $t_{i+1}=t_{i}+h$ 。

龙格-库塔法

在数值分析中，龙格-库塔方法是一系列隐式和显式迭代方法，其中包括欧拉方法，用于时间离散化联立非线性方程的近似解。

龙格-库塔家族中最广为人知的成员通常被称为“RK4”、“经典龙格-库塔方法”或简称为“龙格-库塔方法”。

水箱排水

如图所示，我们考虑一个最初是空水箱的倒截锥形状，很像一个水桶。需要找出水箱的水深 $h$ 的变化，如果水箱充满水源 $Q_{in}$ ，并且底部出口也以可变流量排出一部分积聚水流输出。虽然可以控制 $Q_{in}$ ，但从水箱流出的 $Q_{out}$ 完全取决于桶中的水深 $h$ 。

我们可以假设罐是均匀变细的，横截面高度为 $h$ ， $A_h$ ，从等于底部的面积 $A_{bottom}$ 到罐顶部的 $A_{top}$ 变化。因此， $A_h$ 可以用以下形式的线性函数表示为 $A_{top}$ 和 $A_{bottom}$ ：

$A_{h}=A_{\text {bottom }}+\frac{A_{\text {top }}-A_{\text {bottom }}}{D} h$

溢洪道水流

溢洪道需要针对可能发生的最严重洪水进行设计，并根据该地区水文气象研究获得的集水区水文和降雨输入估算“设计洪水水位线”。水库问题需要已知溢洪道流出量对于给定的流入过程线。通过使用设计洪水过程线作为水库的入流，通过大坝溢洪道的出流，连同水库水位的上升，可通过以下连续性方程计算：

折叠目录

回水海拔剖面

溶解氧和生化需氧量

水位波动

回灌水位剖面

水质问题

Python解

水箱排水问题

import numpy as np
import matplotlib.pylab as plt
g = 9.81
Atop = 1.0
Abottom = 0.5
D = 0.75
Aoutlet = 0.005
Cd = 0.7

dt = 1 # Time step (in seconds
ntimes = int(time_simulation/dt)
print(ntimes)
h=np.zeros(ntimes) # Sequence of water depths in bucket
time = np.zeros(ntimes) # Time markers for plotting
h[0] = 0.0
for n in range (0,ntimes-1):
	 time[n+1] = time[n] + dt
	 if(time[n] <= time_control1):
		 Qin = Qin1
	 elif(time[n] > time_control1 and time[n] <=
		time_control2):
	 Qin = Qin2
	 else:
		 Qin = 0.0
	 Qout = Cd*Aoutlet*np.sqrt(2*g*h[n])
	 Area = Abottom + h[n]*(Atop-Abottom)/D
	 dhdt0 = (Qin - Qout)/Area
	 h1 = h[n]+dhdt0*dt
	 Qout1 = Cd*Aoutlet*np.sqrt(2*g*h1)
	 Area1 = Abottom + h1*(Atop-Abottom)/D
	 dhdt1 = (Qin - Qout1)/Area1
	 dhdt = 0.5*(dhdt0+dhdt1)
	 h[n+1] = dhdt*dt+h[n]
	 if(h[n+1] < 0.01): h[n+1] = 0.0

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(time, h)
ax.set_xlabel('Time (s)')
ax.set_ylabel('Water Depth (m)')
plt.show()

溢洪道水流

回水海拔剖面

溶解氧和生化需氧量

水位波动

回灌水位剖面

水质问题

源代码

Python(普赖斯曼方法解偏微分方程和运动波方程-圣维南流动方程)求棱形流道表面流体

如图所示，对于棱柱形通道的两个部分之间的流量范围，控制方程是从以下考虑推导出来的：

连续性方程，通过考虑截面内的净流入量、流出量和存储量得出，以及
动量方程，考虑流动方向的净力和动量。

尽管对于非棱柱形通道可能有一个更一般的方程可用，但这里给出了棱柱形通道方程，因为以下部分中处理的所有示例都符合这一假设。

我们可以用下面的一对方程写出连续性和动量方程，称为圣维南流动方程：

$\begin{gathered}B \frac{\partial y}{\partial t}+\frac{\partial Q}{\partial x}=q \\\frac{\partial Q}{\partial t}+\frac{\partial(Q u)}{\partial x}+g A \frac{\partial y}{\partial x}=V_x q+g A\left(S_0-S_f\right)\end{gathered}$

理想流体

二维深度平均流的控制方程

流体方程数值解(Python)

横向流入棱柱形通道中运动波动方程

三角通道中运动波近似路由洪水波

运动波逼近开书流域水文图

圣维南流动方程模拟通道中的非定常流动

理想流体流动方程求解

浅盆地二维深度平均流模拟

Python(有限差分法和杜普伊特假设)数值解非承压畜水层和堰底和板桩下稳定渗流

达西定律描述了计算通过多孔介质的流量的基本方程，在三个垂直坐标方向 x、y 和 z 上，可以用以下方式编写

$q_x=-K \frac{\partial h}{\partial x} ; q_y=-K \frac{\partial h}{\partial y} ; \text { and } q_z=-K \frac{\partial h}{\partial z}$

上式中，使用的符号分别为： $q_x$ 、 $q_y$ 、 $q_z$ ，三个方向单位面积的体积流量和 $K$ ，介质的水力传导率，为简单起见，假定为各向同性。等式中的空间导数是给定点 $h$ 上方土壤内测压头或饱和水柱深度相对于三个坐标轴的偏导数。方程中的表达式，将三个坐标方向上的速度与各自的压头梯度相关联，可以代入相应的稳态连续方程，得到以下公式：

$\frac{\partial q_x}{\partial x}+\frac{\partial q_y}{\partial y}+\frac{\partial q_z}{\partial z}=0$

非承压畜水层控制方程

承压畜水层控制方程

垂直面稳态渗流

数值解(Python)

解非承压畜水层中非定常一维和二维水流

堰底和板桩稳定渗流

PreviousPython微积分概率统计和图形论及几何学实现 NextPython主成分分析(简化分子线性输入)化学分子结构

Last updated 1 year ago

Was this helpful?

import numpy as np import matplotlib.pylab as plt b0 = 20.0 s = 2.0 s0 = 0.001 mn = 0.020 Q = 40 yinitial = 1.0 errorallow = 0.0001 def area(y): area=y*(b0+s*y) return area def wetperi(depth): wetperi = b0+2*depth*np.sqrt(1+s*s) return wetperi def hyrad(depth): hyrad = (b0+depth*s)*depth/(b0+2*depth*np.sqrt(1+s*s)) return hyrad def topwidth(depth): topwidth=b0+2*s*depth return topwidth iter = 0 max_error = 1.0 yn = yinitial while (max_error > errorallow): iter = iter+1 Fyn = area(yn)**1.6667/wetperi(yn)**0.6667-mn*Q/ np.sqrt(s0)

import numpy as np import matplotlib.pylab as plt g = 9.81 Atop = 1.0 Abottom = 0.5 D = 0.75 Aoutlet = 0.005 Cd = 0.7 dt = 1 # Time step (in seconds ntimes = int(time_simulation/dt) print(ntimes) h=np.zeros(ntimes) # Sequence of water depths in bucket time = np.zeros(ntimes) # Time markers for plotting h[0] = 0.0 for n in range (0,ntimes-1): time[n+1] = time[n] + dt if(time[n] <= time_control1): Qin = Qin1 elif(time[n] > time_control1 and time[n] <= time_control2): Qin = Qin2 else: Qin = 0.0 Qout = Cd*Aoutlet*np.sqrt(2*g*h[n]) Area = Abottom + h[n]*(Atop-Abottom)/D dhdt0 = (Qin - Qout)/Area h1 = h[n]+dhdt0*dt Qout1 = Cd*Aoutlet*np.sqrt(2*g*h1) Area1 = Abottom + h1*(Atop-Abottom)/D dhdt1 = (Qin - Qout1)/Area1 dhdt = 0.5*(dhdt0+dhdt1) h[n+1] = dhdt*dt+h[n] if(h[n+1] < 0.01): h[n+1] = 0.0 fig = plt.figure() ax = fig.add_subplot(1, 1, 1) ax.plot(time, h) ax.set_xlabel('Time (s)') ax.set_ylabel('Water Depth (m)') plt.show()