平稳模型到动态模型的转换

Question

我尝试使用以下方程式在Python中模拟水库的流量：

def dRho(rho_y, t):
    P_r = 10e5
    rho_r = 900
    L = 750
    H = 10
    W = 150
    A = H * W
    V = A * L
    fi = 0.17

    k = 1.2e-13
    c = 12.8e-9
    mu = 2e-3

    N = 1
    dV = V/N
    dx = L/N

    P_in = P_r
    rho_in = rho_r

    P_w = 1e5
    rho_w = rho_r* np.exp(c*(P_w-P_r))

    P = P_r + (1/c) * np.log(rho_y/rho_r)
    dP_in = P - P_in
    dP_uit = P_w - P

    Q_in = (-A*k/mu)*(dP_in/dx)
    Q_uit = (-A*k/mu)*(dP_uit/dx)
    return (Q_in*rho_in - Q_uit*rho_y)/dV

t0 = np.linspace(0,1e9, 1e9/100)
rho0 = 900
solve = odeint(dRho, rho0, t0)
plt.plot(t0,solve, '-')
plt.show()

到目前为止，我能够使用以下公式求解方程的平稳形式： 我曾经尝试过使用numpy数组和切片将静态脚本转换为动态版本，但是我没有得到任何接近于动态求解方程的地方。我主要在使用/实现带有下标的变量（即I、I+1和I-1）时遇到问题。你知道吗

应用动力学方程的尝试：

def dRho(rho_y, t):
    # reservoir 
    P_r = 10e5
    rho_r = 900
    L = 750
    H = 10
    W = 150
    A = H * W
    V = A * L
    fi = 0.17
    # fluid 
    k = 1.2e-13
    c = 12.8e-9
    mu = 2e-3
    # components
    N = 10
    dV = V/N
    dx = L/N

    # inlet
    P_in = P_r
    rho_in = rho_r
    # outlet
    P_w = 1e5    
    rho_w = rho_r* np.exp(c*(P_w-P_r))


    rhoi = rho_y[1:-1] # i
    rhop = rho_y[2:]   # i + 1
    rhom = rho_y[:-2]  # i - 1

    # compressibility
    Pi = P_r + (1/c) * np.log(rhoi/rho_r)  # i
    Pm = P_r + (1/c) * np.log(rhom/rho_r) #i - 1

    # Darcy
    Qi = (-A*k/mu)*((Pi)/dx) # i
    Qp = (-A*k/mu)*((Pm-Pi)/dx) # i + 1

    return (Qp*rhop - Qi*rhoi)/dV

t0 = np.linspace(0,1e9, 1e9/100)
rho0 = np.ones(10)*900
solve = odeint(dRho, rho0, t0)
plt.plot(t0,solve, '-', label='dRho')
plt.legend(loc='upper right')
plt.show()

编辑1： 我添加了一些脚本中缺少的初始值。我的新代码是：

def dRho(rho_y, t):
    # reservoir 
    P_r = 10e5
    rho_r = 900
    L = 750
    H = 10
    W = 150
    A = H * W
    V = A * L
    fi = 0.17
    # fluid 
    k = 1.2e-13
    c = 12.8e-9
    mu = 2e-3
    # components
    N = 8
    dV = V/N
    dx = L/N

    # inlet
    P_in = P_r
    rho_in = rho_r
    # outlet
    P_w = 1e5    
    rho_w = rho_r* np.exp(c*(P_w-P_r))


    rhoi = rho_y[1:-1] # i
    rhop = rho_y[2:]   # i + 1
    rhom = rho_y[:-2]  # i - 1

    # compressibility
    Pi = P_r + (1/c) * np.log(rhoi/rho_r) # i
    Pm = P_r + (1/c) * np.log(rhom/rho_r) #i - 1

    # Darcy
    Q0 = 0
    Qi = (-A*k/mu)*((Pi)/dx) # i
    Qp = (-A*k/mu)*((Pm-Pi)/dx) # i + 1

    out = np.empty(N+1)
    out[0] = 0
    out[1:N] = (Qp*rhop - Qi*rhoi)/dV
    out[N] = 0
    return out

t0 = np.linspace(0,1e9, 1e9/200)
N = 8
rho0 = np.ones(N+1)*900
solve = odeint(dRho, rho0, t0)
plt.plot(t0,solve, '-', label='dRho')
plt.legend(loc='upper right')
plt.show()

不幸的是，这个脚本的输出不是我想要的，老实说，我没有看到我的错误。你知道吗

电流输出：

虽然所需的输出应该更像这样：

我通过对I=8的p\u I、Q\u I和dRho/dt进行硬编码来实现所需的输出。你知道吗

编辑2：

我通过创建一系列的3个for循环来获得所需的输出，每个for循环对应于前面提到的等式。然而，我的脚本需要相当长的时间来执行，因为大O国=n^3

这是我的密码：

def dRho(rho_y, t, N):
    # reservoir 
    P_r = 10e5
    rho_r = 900
    L = 750
    H = 10
    W = 150
    A = H * W
    V = A * L
    fi = 0.17
    # fluid 
    k = 1.2e-13
    c = 12.8e-9
    mu = 2e-3
    # components
    N = N
    dV = V/N
    dx = L/N

    # inlet
    P_in = P_r
    rho_in = rho_r
    # outlet
    P_w = 1e5    
    rho_w = rho_r* np.exp(c*(P_w-P_r))

    P = np.empty(N+1)*10e5
    Q = np.ones(N+1)
    out = np.empty(N+1)

    P[0] = P_w
    for i in range(1,N):
        P[i] = P_r + (1/c) * np.log(rho_y[i]/rho_r)

    P[N] = P_r + (1/c) * np.log(rho_y[N]/rho_r)

    Q[0] = 0
    for i in range(1,N):
        Q[i] = (-A*k/mu)*((P[i-1] - P[i])/dx)
    Q[N] = (-A*k/mu)*((P[N]-P_r)/dx)

    out[0] = 0
    for i in range(1,N):
        out[i] = ((Q[i+1]*rho_y[i+1] - Q[i]*rho_y[i])/dV) 
    out[N] = 0

    return out

t0 = np.linspace(0,1e9, int(1e9/200))
N = 100
rho0 = np.ones(N+1)*900
solve = odeint(dRho, rho0, t0, args=(N,))
print(solve[:,0])


rho_i = np.ones(len(t0))*900
plt.plot(t0,solve[:,1:len(rho0)-1], '-', label='dRho')
plt.plot(t0,rho_i, '-', label='dRho_initial')
# plt.legend(loc='upper right')
plt.show()

有没有一种方法可以在合理的时间内执行时仍然自动执行。你知道吗