计算无量纲压力（对应0-1无量纲时间）出现负值的求解方法咨询

我尝试计算无量纲压力并绘制其与0到1范围内无量纲时间的关系曲线，使用公式P_dim.append((2*pi*PP['k']*res['h']/(-Q[840]*pvt['muo']*30))*(3000-Pn_imp[840]))计算时，得到了接近-0.9的错误负值，希望了解正确的计算方式。

相关Python代码如下：

P = np.zeros((Nx,Ny))
m=len(t)          # time range for pressure calculation 
print(Pint)
print(np.matmul(Act,Pint) )
A_term_imp=Act + T + J    # first term of pressure equation using implicit method
B_term_imp = np.matmul(Act,Pint) + Q  
Pn_imp = np.linalg.solve(A_term_imp, B_term_imp)  
print(np.linalg.solve(A_term_imp, B_term_imp))
p_840=[]
Pwl=[]
time=[]
P_dim=[]
for j in range(1,m+1):
    A_term_imp=Act + T + J    # first term of pressure equation using implicit method
    B_term_imp = np.matmul(Act,Pint) + Q         # second term of pressure equation using implicit method   
    Pn_imp = np.linalg.solve(A_term_imp, B_term_imp)
    # Pwl.append(Pn_imp)
    # time.append(0.01*j)
    if j==1:
        print(Pn_imp)
    # average = Average(Pn_imp)
    # print(average)
    print('Pressure matrix at time step '+str(j*dt)+' days, implicit method is ',Pn_imp, ' psi', '\n')
    Pint = Pn_imp
    p_840.append(Pn_imp[840])
    
    P_dim.append((2*pi*PP['k']*res['h']/(-Q[840]*pvt['muo']*30))*(3000-Pn_imp[840])     #dimenssionless pressure 
 
            
    for j in range(1,m):
        P=Pn_imp
        P = np.reshape(P, (dy.shape[0], dy.shape[0]))

问题排查与修正 #

1. 符号逻辑错误 #

无量纲压力为负的核心原因是公式符号匹配错误：

• 如果Q[840]是生产井的产量（流体流出，定义为正值），公式中的-Q[840]会引入负号；而生产时井底压力Pn_imp[840]低于初始压力3000psi，3000-Pn_imp[840]为正，正负相乘得到负值。
• 标准无量纲压力公式中，生产井产量取正值，无需额外添加负号。

2. 公式参数与单位问题 #

油藏工程中无量纲压力的标准定义为：
P_D = (141.2 * k * h) / (μ * B * q) * (P_i - P_wf)
需要注意：

• 体积系数B：如果Q[840]是地面产量，必须乘以原油体积系数B（地下产量无需），你的原公式遗漏了该参数。
• 单位转换系数：当参数单位为k(md)、h(ft)、μ(cp)、q(bbl/d)、P(psi)时，转换系数为141.2，而非原公式中的30。

3. 代码嵌套循环冗余 #

外层循环内嵌套了另一个for j in range(1,m)循环，会导致每次时间步都重复执行多次矩阵重塑，浪费资源且可能引发变量混乱，需移除。

修正后的代码示例 #

import numpy as np

P = np.zeros((Nx,Ny))
m = len(t)          # 压力计算的时间步数
print(Pint)
print(np.matmul(Act, Pint))
A_term_imp = Act + T + J    # 隐式方法压力方程第一项
B_term_imp = np.matmul(Act, Pint) + Q  
Pn_imp = np.linalg.solve(A_term_imp, B_term_imp)  
print(np.linalg.solve(A_term_imp, B_term_imp))

p_840 = []
time = []
P_dim = []
# 获取原油体积系数，若为地下产量则取1.0
B = pvt.get('Bo', 1.0)  

for j in range(1, m+1):
    A_term_imp = Act + T + J    
    B_term_imp = np.matmul(Act, Pint) + Q       
    Pn_imp = np.linalg.solve(A_term_imp, B_term_imp)
    
    if j == 1:
        print(Pn_imp)
    
    print(f'隐式方法下，时间步{j*dt}天的压力矩阵为：{Pn_imp} psi\n')
    Pint = Pn_imp
    p_840.append(Pn_imp[840])
    
    # 修正后的无量纲压力计算
    numerator = 141.2 * PP['k'] * res['h']
    denominator = Q[840] * pvt['muo'] * B
    # 若Q定义为注入正、生产负，需在denominator前加负号
    p_dim_val = numerator / denominator * (3000 - Pn_imp[840])
    P_dim.append(p_dim_val)
    
    # 移除冗余嵌套循环，直接单次重塑矩阵
    P = np.reshape(Pn_imp, (dy.shape[0], dy.shape[0]))

内容的提问来源于stack exchange，提问作者Hana Bachi