简介面向石油工程研究人员与学生提供一套基于MATLAB的黑油模型油藏数值模拟实现核心代码涵盖黑油模型、两相油水模型、三相模型及水热模型可辅助理解离散化、状态方程、Darcy流动方程与时间推进等关键环节。包体共67个文件以.m脚本为主体另含txt说明、mat数据及Eclipse格式的smspec/unsmry输出文件压缩包约173KB目录按components、utils、examples等模块划分便于检索和对照学习。示例包含SPE1、SPE9、SPE10、Egg及Norne等经典模拟案例能够从组件构建到实例仿真完整展示黑油模型应用流程用于预测油藏压力分布、产量与采收率并配有README与多种工具脚本适合有一定MATLAB基础的学习者快速上手。已有697人学习下载对深入理解油藏数值模拟与黑油模型实现具有较好的参考价值。1. 从 MATLAB 代码包说起黑油模型为什么值得专门拆一套拿到 ad-blackoil 这个资源包第一反应是它把 MATLAB Reservoir Simulation ToolboxMRST里黑油模型相关的那条线单独抽了出来。MRST 本身是个庞然大物但真正做油藏数值模拟的人都知道日常最常用的核心就是黑油模型这一套从单相水驱到三相黑油再到带热水驱、聚合物这些扩展项。这个包把 model 层、components 层、utils 层和 examples 层分得清清楚楚相当于把 MRST 的模型开发框架以黑油模型为主线重新组织了一遍。对 5 年以上从业者来说这个包的价值不在于能跑——MRST 官方仓库本来是能跑的——而在于它把黑油模型的实现路径单独暴露出来了。比如 ThreePhaseBlackOilModel.m 和 TwoPhaseOilWaterModel.m 的继承关系、computeFlashBlackOil.m 里溶解气比的计算逻辑、getFluxAndPropsOil_BO.m 里油相达西流量的组装方式这些都是做二次开发时真正要改的地方。对新手而言SPE1、SPE9、SPE10 三个标准算例加上 Egg 和 Norne 两个实际模型足够把黑油模型从理论到工程完整过一遍。这篇文章就按模型体系 → 标准算例 → 时间步与性能 → 复杂模型 → 验证手段这条线来拆。2. 黑油模型体系拆解从 WaterModel 到 ThreePhaseBlackOilModel 的递进关系2.1 黑油模型的核心假设与组分逻辑黑油模型之所以叫黑油是因为它把地下流体简化成三个拟组分油、气、水。天然气在处理上有两个去处自由气Free Gas和溶解气Dissolved Gas溶解气溶于油相中落泡点压力以下时从油相中析出这就是黑油模型与组分模型最关键的区别——黑油不追踪组分摩尔分数只追踪溶解气油比Rs这个状态量。2.2 资源包里的五个模型类各自边界在哪里打开 models 目录能从类的命名直接看出 MRST 的模型设计思路。五个核心模型类分别是WaterModel、TwoPhaseOilWaterModel、ThreePhaseBlackOilModel、GenericBlackOilModel和WaterThermalModel。这里的继承关系不是随意叠的每多一个类就多一组物理过程模型类相态适用场景关键方程文件WaterModel单相水示踪剂测试、水驱前缘分析equationsWater.mTwoPhaseOilWaterModel油水不考虑溶解气的水驱油藏equationsOilWater.mThreePhaseBlackOilModel油气水溶解气驱、气顶油藏equationsBlackOil.mGenericBlackOilModel油气水组分扩展聚合物、低盐度水驱updateStateBlackOilGeneric.mWaterThermalModel水温度热水驱、地热equationsWaterThermal.m以两相油水模型为例它不计算气相所以getFluxAndPropsOil_BO.m里只需要组装油相和水相的相对渗透率而三相模型必须在每个网格单元里判断自由气是否存在——如果压力高于泡点自由气饱和度为零此时方程里气相流动项直接消掉但溶解气比仍然是一个未知量。2.3 equationsBlackOil.m 的方程组装逻辑从残差到雅可比真正做二次开发的人最关心的应该是equationsBlackOil.m。它把每个网格单元的质量守恒方程组织成残差形式然后交给 MRST 的 Newton 求解器迭代。核心代码逻辑如下function eqs equationsBlackOil(state, model, dt, drivingForces) % 计算每个网格单元的达西流量 [flux, fluxFace, fluxUpstream] getFluxAndPropsOil_BO(model, state); % 组装油、气、水三相的质量守恒残差 eqs assembleMassBalanceEquations(model, state, flux, ... state.flux, drivingForces); % 计算溶解气比 Rs 的约束方程 if model.disgas rsEq computeDissolvedGasEquation(model, state); eqs [eqs; rsEq]; end % 加上井源汇项 eqs addWellContributions(model, eqs, state, drivingForces); end这段代码里最容易被忽略的是model.disgas这个标志位。它决定黑油模型是否启用溶解气计算——关闭时三相模型退化成带自由气的油水模型方程数量直接少一组。实际项目里如果遇到不收敛先检查这个标志位是否匹配油藏实际流体性质。2.4 毛管压力与相对渗透率的处理方式getCapillaryPressureBO.m处理的是油水两相的毛管压力。MRST 的默认做法是读取state.s里的饱和度然后调用model.rock里的Pc函数句柄。这个函数句柄既可以是解析式如 Brooks-Corey也可以是查表插值。注意黑油模型默认把气相毛管压力设为零如果有气水毛管压力数据需要自己扩展这个文件。3. SPE1、SPE9 实战标准算例的正确跑法与结果解读3.1 SPE1 三层气顶油藏案例的参数矩阵SPE1 是黑油模型绕不开的基准算例。它模拟一个带气顶的三层油藏顶部注气、底部采油。资源包里的simulateSPE1.m把整个流程封装好了但作为工程师不能只点运行要能说清每步在干什么% 读取 SPE1 网格和属性10x10x3笛卡尔网格 deck readEclipseDeck(SPE1.DATA); % 转换为 MRST 网格结构激活有效单元 G computeGeometry(processGRDECL(deck.GRID)); % 初始化三相黑油模型 model ThreePhaseBlackOilModel(G, deck.ROCK, deck.FLUID, ... disgas, true, surfstd, true); % 设置初始状态气顶区域自由气饱和度非零 state0 initResSol(G, deck, rs, deck.INITIAL.RS);这里的surfstd参数容易被忽略它控制溶解气比 Rs 是否按地面标准条件折算。SPE1 的流体数据是标准的所以开不开启结果差异不大但对实际油藏数据油、气的地面体积系数定义方式会直接影响产量曲线的数值。3.2 SPE9 的坑高含水率下的时间步崩溃SPE9 比 SPE1 难收敛得多。它包含 25 口生产井含水率快速上升很多初学 MRST 的人跑 SPE9 会在中后期遇到时间步截断问题。资源包里的blackoilTutorialSPE9.m和simulateSPE9.m是两种不同的处理路线% simulateSPE9.m 采用固定时间步策略 timesteps repmat(30*day, 1, 120); % 30 天一步共 120 步 [wellSols, states] simulateSchedule(state0, model, ... schedule, OutputLevel, 1); % blackoilTutorialSPE9.m 采用自适应时间步 schedule.step.control ones(numel(timesteps), 1); schedule.step.val timesteps;实际跑下来固定步长在 2000 天附近会卡住因为气饱和度在局部快速变化Newton 迭代不收敛。自适应步长能自动把 30 天切成 1 天甚至更小代价是总步数接近翻倍。在真实项目中我的建议是先用blackoilTutorialSPE9.m跑通流程再用simulateSPE9.m做精度验证两个结果对比看离散误差。3.3 从 WellSols 提取生产指标跑完算例后wellSols结构里存了每口井每个时间步的产量数据。提取指标时注意单位换算% 读取第一口生产井的产油速率 qOs cellfun((ws) ws(1).qOs, wellSols, UniformOutput, false); qOs cell2mat(qOs); % 单位m^3/day地面条件 % 累加得到累计产油量 cumO cumsum(qOs .* timesteps / day); % 按井筒流压约束检查是否达到产能极限 bhp cellfun((ws) ws(1).bhp, wellSols, UniformOutput, false); bhp cell2mat(bhp);这里最常犯的错是把qOs当成地下体积流量。MRST 里带s后缀的字段都是地面标准条件下的流量如果直接拿去做物质平衡会差一个体积系数倍率。3.4 SPE10 模型与 upscaling 的必要性SPE10 模型有 120 万网格60x220x85直接跑全模型对内存和时间都不友好。资源包里的tenthCSP_Model_I.m是 SPE10 上半部分的渗透率场取自第 1 到第 35 层Tarbert 组。实测在 16GB 内存的机器上跑 5 年生产史大约需要 2-3 小时。如果嫌慢常见做法是空间粗化upscaling比如把 60x220 粗化到 30x110渗透率用 Harmonic Mean 做等效。4. 时间步控制与性能优化从 timestepControlDemo 到 Mex 加速4.1 自适应时间步策略的核心逻辑timestepControlDemo.m值得仔细研究。MRST 的默认时间步控制器是SimpleTimeStepSelector它的逻辑是Newton 迭代时如果遇到过大的残差变化或超过最大迭代次数就把时间步减半退回重试如果收敛得很快则线性增大时间步。% timestepControlDemo.m 中的关键配置 model.controlTimeStepSelector TimeStepSelector(maxGrowth, 1.5, ... maxSteps, 10, minSteps, 2);maxGrowth设为 1.5 表示时间步最多增大 50%maxSteps和minSteps控制每步内的 Newton 迭代次数上下限。实际调试时如果油藏在注水突破后含水率快速攀升建议把maxGrowth降到 1.2收敛稳定性明显改善。4.2 Mex 加速的适用边界blackoilTutorialMexAcceleration.m把核心物理计算编译成 Mex 文件。启用方法简单mrstModule add mex; mexSetup();但要注意Mex 加速只对getFluxAndPropsOil_BO.m这类稠密计算有效对模型组装和井计算反而可能更慢。实测下来SPE9 这种 2 万网格规模的算例加速比约 2-3 倍SPE10 全模型能到 5 倍以上。如果是刚入门调试不建议第一时间开 Mex因为出错了很难定位。4.3 内存优化与稀疏矩阵复用黑油模型的性能瓶颈集中在雅可比矩阵组装上。MRST 每次 Newton 迭代都会重建雅可比矩阵这占了总耗时的一半以上。常见的优化手段是复用稀疏矩阵的稀疏模式symbolic factorization reuse% 预先计算稀疏模式后续迭代只更新数值值 Jac model.getLinearJacobian(state, state0, dt, drivingForces); % 在循环外定义求解器选项 ls LinearSolver(type, ilu, reuse, true);LinearSolver的reuse参数决定是否复用预处理矩阵。对于长期生产模拟几百个时间步、几千个 Newton 迭代打开这个开关能省掉大量重复计算。实测 SPE1 全周期模拟从 30 分钟缩短到 12 分钟左右。5. 复杂场景落地Egg 模型、Norne 模型与多段井模型5.1 Egg 模型注采优化的理想测试平台Egg 模型是一个 101x20x7 的 synthetic 油藏101 口井10 注 91 采专门用来测试注采优化算法。资源包里的eggExample.m演示了最基本的定液量生产方案。这个模型的价值在于它对异质性敏感——渗透率场是层状分布的注水推进不均匀天然适合测试闭环优化策略。5.2 Norne 模型真实油田数据的复杂之处fieldModelNorneExample.m加载的是 Norne 油田的公开数据。这个数据集的难点在于它包含断层和复杂井轨迹processGRDECL的断层自动处理有时候会产生活性网格与死网格的邻接问题。跑这个模型需要注意两点一是确保mrstModule add deckformat加载了完整的 Eclipse 数据解析模块二是在输出结果时善用plotWellSols辅助函数否则多口井的产量曲线会挤在一起。5.3 Multisegment Well 模型的实现逻辑multisegmentWellExample.m是多段井模型MSW的示例。它把水平井或复杂结构井分成多段每段单独计算压降和流量比传统直井模型的井指数精度高得多。MRST 中启用多段井的代码结构% 读入 Eclipse 风格的井定义 wellSpecs getEclipseWells(deck, G, multisegment, true); % 装配多段井模型 model addMultisegmentWells(model, wellSpecs); % 求解 [ws, states] simulateSchedule(state0, model, schedule);多段井的核心是每段之间的流量平衡方程。MRST 在处理水平井时把井筒沿轨迹离散成若干节点每个节点有自己的 BHP 和流量段间压降用 Fanning 摩擦因子公式计算。实际应用时需要注意段长不要太短否则摩擦压降会在段间反复迭代拖慢收敛速度。5.4 完井射孔层位设定与边界条件处理井的射孔层位设定也是个容易出问题的地方。MRST 的getEclipseWells默认读取 Eclipse 数据文件中的COMPDAT关键字如果没有这个关键字就需要手动指定% 手动定义一口射开 1-3 层的生产井 W addWell(G, rock, 50, 50, 1:3, ... Name, P1, comp_i, 50, comp_j, 50, ... sign, -1, Type, bhp, Val, 200*barsa);这里sign-1表示生产井注采方向Typebhp表示定井底流压生产。注意黑油模型里射孔层位跨越多个网格时各层的流量分配是按各层 Kh 值自动加权的不需要手动设置分层产液比例。6. compare 目录的用武之地回归测试与批量基准对比资源包里的examples/compare目录是很多人会忽略但实际价值很高的部分。compare目录下存放了不同模型在同一算例上的运行结果对比比如 SPE1 在ThreePhaseBlackOilModel与GenericBlackOilModel下的差异。这为模型开发提供了一条标准化的回归测试路径。常规做法是写一个测试脚本把修改后的模型输出与基线结果做差检查最大误差是否在可接受范围内% 读取两个模型在相同时间点的油藏状态 state_ref states_ref{end}; % 基线模型末时刻状态 state_mod states_mod{end}; % 修改后模型末时刻状态 % 计算压力场的最大偏差 pDiff max(abs(state_ref.pressure - state_mod.pressure)) / barsa; assert(pDiff 0.01, 压力场偏差超过 0.01 bar); % 计算饱和度场的最大偏差 sDiff max(abs(state_ref.s(:,1) - state_mod.s(:,1))); assert(sDiff 1e-3, 油饱和度偏差超过 0.1%%);断言阈值不能拍脑袋定。SPE1 这类小规模算例压力偏差应该在 1e-2 bar 以内SPE9 这种强非线性问题放宽到 1e-1 barEgg 模型由于井数多、推进明显饱和度偏差放宽到 1e-2 是合理的。另一个实用技巧是用getReportTimings.m分析各模块耗时占比。它能输出每个时间步里模型组装、线性求解、井计算各自花的时间。线性求解占比超过 70% 时优先考虑换预处理器组装时间占比高时考虑开 Mex 或优化网格邻接结构。这个脚本是判断性能瓶颈方向最直接的工具。最后说一个验证黑油模型结果是否物理合理的标准做法物质平衡检查。跑完模拟后用calculateHydrocarbonsFromStatusBO.m计算油藏内油、气、水的地下体积配合累计产出和注入数据验证守恒性。黑油模型的物质平衡误差通常应该在 1e-6 量级如果差到 1e-3 以上且误差来源不是井数据插值那就要回到模型配置检查disgas标志位或 PVT 表插值是否合理了。本文还有配套的精品资源点击获取