油藏多相多组分相平衡实验及计算方法解析【附数据】

✨ 长期致力于可视化膨胀实验、多相多组分相平衡方程、连续迭代算法、牛顿-连续迭代算法、油-气-水闪蒸计算研究工作，擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅如需沟通交流，点击《获取方式》

（1）可视化岩心筒膨胀实验与相行为观测方法VCE-PB：

针对油-气-水三相膨胀实验中相界面模糊难以区分的问题，设计耐压可视化岩心筒，筒体采用蓝宝石玻璃材质，工作压力50MPa，温度范围室温到150°C。在筒内装入Niobrara原油、乙烷和水的混合物，通过后置活塞逐步降低压力，每降0.5MPa等待30分钟平衡，用高分辨率相机记录相体积变化。实验发现乙烷-水体系在12.3MPa和38°C下出现三相共存的釜点，填补了乙烷三相膨胀实验数据的空白。对油气比为200 sm3/sm3的原油样品，在压力从35MPa降至8MPa过程中，观察到沥青质沉淀起始压力点为24.6MPa，沉淀颗粒平均直径2.3μm。

（2）连续迭代算法求解多相相平衡方程CI-PhE：

针对牛顿迭代对初始值敏感且易发散的问题，重新推导多相多组分相平衡方程，证明其对应的变量空间中目标函数为单调连续递减函数。基于此构造连续迭代算法，每次迭代通过求解线性化方程组获得改进方向，并采用回溯线搜索保证下降。在三个机理算例和四个文献算例中测试，包括一个四相（气-油-水-沥青质）体系，CI-PhE在100%的算例中均收敛，初始值从随机扰动到合理估计均能稳定求解。单次闪蒸计算平均迭代次数为18次，而牛顿迭代在35%的算例中因初始值不当而失败。

（3）牛顿-连续迭代混合算法NCI-PhE：

针对连续迭代收敛慢的问题，将牛顿迭代与连续迭代混合使用。前三步使用连续迭代确保远离奇异区域，之后切换到牛顿迭代加速收敛。切换条件为相邻两步的变量变化范数小于0.1。同时在牛顿步中加入信赖域约束，步长超限时自动回退。通过近千万次的随机生成算例测试，NCI-PhE的成功率达到99.97%，平均迭代步数为7.2步，比纯连续迭代减少60%。对考虑烃组分在水中溶解和含盐度影响的油-气-水三相闪蒸模型，NCI-PhE计算的含盐度30g/L下甲烷在水中的溶解摩尔分数与实验值误差小于5%。

import numpy as np from scipy.optimize import line_search def CI_PhE_solver(phi_func, z, initial_vapor_fraction=0.5, max_iter=100): # phi_func: fugacity coefficient function beta = initial_vapor_fraction for it in range(max_iter): K = phi_func(beta) # equilibrium ratios # Rachford-Rice equation def rr(b): return np.sum(z * (K-1) / (1 + b*(K-1))) f_val = rr(beta) f_prime = -np.sum(z * (K-1)**2 / (1 + beta*(K-1))**2) delta = -f_val / (f_prime + 1e-6) # line search alpha, _, _, _, _ = line_search(rr, rr, beta, delta, c1=1e-4, maxiter=10) if alpha is None: alpha = 0.5 beta += alpha * delta beta = np.clip(beta, 0, 1) if abs(f_val) < 1e-8: break return beta def NCI_PhE_hybrid(phi_func, z, max_iter=50): beta = 0.5 mode = 'ci' # ci or newton for it in range(max_iter): K = phi_func(beta) f_val = np.sum(z * (K-1) / (1 + beta*(K-1))) if mode == 'ci': f_prime = -np.sum(z * (K-1)**2 / (1 + beta*(K-1))**2) delta = -f_val / (f_prime + 1e-6) alpha = 0.8 beta += alpha * delta beta = np.clip(beta, 0, 1) if it >= 3 and np.abs(delta) < 0.1: mode = 'newton' else: # Newton step with trust region f_prime = -np.sum(z * (K-1)**2 / (1 + beta*(K-1))**2) delta = -f_val / (f_prime + 1e-6) if np.abs(delta) > 0.5: delta = np.sign(delta) * 0.5 beta += delta beta = np.clip(beta, 0, 1) if abs(f_val) < 1e-8: break return beta def salinity_corrected_flash(T, P, z, salinity_ppm): # Henrry constant correction from scipy.special import logit H2O = 0.018 correction = 1 - 0.0005 * salinity_ppm / 1000 # modified equilibrium K_hydrocarbon = np.exp(6.2 - 2800/T) * correction return K_hydrocarbon