水电站阀门关闭时的水团振荡

Question

我设法解决了这个任务，下面是解释。谢谢大家的努力。你知道吗

所以我需要解两个相关的微分方程，它代表了流体流量的变化（1），以及调压室（2）内上水位和水位之间的高差随时间的变化。你知道吗

这个特别的例子是一个水力发电厂（HPP）。如果他们关闭阀门，开关就在供水管的末端，启动压力管道，供水管内的水团就会开始振荡。因为在供水管道的末端有一个缓冲罐，它有一个功能，当振荡持续的时候，它可以把水团带到自己的体内，所以不会有任何严重的后果，比如水锤。你知道吗

然后将水电站分为三部分。第一部分为上水，第二部分为供水管道（通过大坝与上水相连），第三部分为调压井（与供水管道的上侧相连）。当振荡持续时，水团从上层水向调压井移动，反之亦然。你知道吗

所以这里有，一些关于这里计算的解释。。。你知道吗

所有的变量都在代码execpt中定义了当量直径、当量摩擦力和供应管道的长度，这些都是手工计算的。如果你想的话，我也可以在代码中为它们写方程。。。你知道吗

这个方程表示：运动方程（1）和连续性方程（2）。你知道吗

  dQt/dt = (g*At/Lt) * (yt - c*Qt*np.abs(Qt)/At**2) #(1)

  dy/dt = -Qt/As #(2)

流动的初始条件是：

  Qt = 18 m^3/s

高度为：

  y = 30 m

时间变化：从0 to 1000 s

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Variables for calculation:
g = 9.81 # gravity acceleration [m/s**2]
Dt = 2.776 # equivalent diameter of supply pipeline, from upper water to surge tank [m]
Ds = 3 # diameter of surge tank [m]
At = np.pi*Dt**2/4. # equivalent cross section of supply pipeline [m**2]
As = np.pi*Ds**2/4. # cross section of surge tank [m**2]

f_l = 0.011155 # equivalent friction coefficient of supply pipeline [/]
Lt = 9276 # lenght of supply pipeline, from upper water to surge tank  [m]

c = (1+f_l*Lt/Dt)/(2*g) # loss coefficient [s**2/m]

# Function si defined as:
def f(Qy, t):
    Qt = Qy[0]
    yt = Qy[1]

    dQdt = (g*At/Lt) * (yt - c*Qt*np.abs(Qt)/At**2)
    dydt = -Qt/As 

    return [dQdt, dydt]

# Where are:
# Qt - volumetric flow of fluid [m**3/s]
# yt - height difference between upper water and water level inside surge tank [m]

# Time interval:
t = np.linspace(0, 1000, 2000)

# Initial condition of volumetric flow, Qt and height difference between upper water and water level inside surge tank, yt:
Qy0 = [18, 30]

# Calling the function for time interval, t:
Q = odeint(f, Qy0, t)

# Graph drawing:
fig = plt.figure()
fig.set_facecolor("w")

plt.title('Oscillations')
plt.axhline(c='gray', lw=0.5)

plt.plot(t, Q[:,0], label='Change of flow in time')
plt.plot(t, Q[:,1], label='Change of height difference in time')

plt.xlabel('Time interval, t')
plt.ylabel('Flow, Qt \nheight difference, y')

legend = plt.legend(loc='lower right', shadow = False, fontsize = 'small')

plt.show()