平时做页岩气、致密气的生产动态模拟拿到现场给的压裂方案和测井数据最常见的需求就是“把这个井的降压开采过程算一遍看产量能到多少”。这类问题不少人会图省事用一块均匀渗透率的长方形储层代替压裂改造区跑出来的产量曲线倒是挺漂亮可拿给搞现场的人看对方第一句往往就是“你裂缝呢”。真正把压裂井写进模型之后你会发现裂缝、基质、井筒三者的压力传播速度完全不在一个量级这不是拍扁成均匀介质能糊弄过去的。这篇算例就是我基于COMSOL搭建的一套压裂井降压开采数值模拟模型主要用于熟悉多物理场耦合、低维裂缝域处理以及瞬态生产动态的完整建模流程适合地质力学、石油工程、水利和环境方向的研究生也适合刚接触COMSOL、想找一个能落地的工程案例来练手的工程师。整套算例不是我凭空捏的思路是照着页岩气水平井分段压裂的典型场景来的水平井筒贯穿储层横向压裂裂缝间隔分布井底流压从初始地层压力逐步降到生产压力气体从致密基质流进裂缝再由裂缝汇入井筒。COMSOL里主要用到地下水流模块中的达西定律接口配合必要的流固耦合表达式把裂缝处理成一维低维域基质用二维域表示这样既保留裂缝的高导流特征又不至于被精细几何拖垮网格。下面从物理背景、建模决策、实操步骤、结果验证到调试经验一层层拆开讲。1. 降压开采为什么值得单独建模而不是把储层拍平了事1.1 降压开采在数值模型里到底“算”的是什么降压开采这四个字听起来很简单就是把井底压力人为压低让地层压力高于井底形成压差驱动流体流入井筒。但搬到COMSOL里它实际是一个典型的非定常渗流过程数学上就是解达西方程∇·( - k/μ ∇p ) 0域内无源汇时你若只跑稳态解得到的只是某一时刻的压力场没法回答“生产十年后产量掉到多少”这种生产动态问题。所以算例用的是瞬态求解初始条件直接给整个储层一个均匀的原始地层压力井筒边界在t0时刻突然下降到生产压力。压力的下降会在储层内部形成不断向外扩散的“压降漏斗”这个漏斗的扩展速度就取决于储层渗透率、孔隙度、流体压缩系数以及裂缝的导流能力。这种模拟的关键点在于“压力波传播速度”。致密储层基质渗透率常常在0.0001 mD到0.01 mD这个量级压力波传播极慢几年时间可能也就影响井周围几十米而裂缝网络的渗透率能够达到几十甚至几千毫达西压力波可以沿着裂缝迅速传到数百米以外。如果不用裂缝域建模均匀介质模型会让压力波传播速率变成某种“平均”既低估了裂缝近井区的高速导流又高估了基质的供气能力产量曲线自然对不上现场实测。1.2 压裂井几何拓扑裂缝是“低维域”井是“线”页岩气水平井分段压裂后真正起导流作用的是数十到数百条激活的天然裂缝、次生裂缝网络以及少数几条主裂缝。把这些真实宽度只有毫米级、长度却有百米的裂缝全部用实体几何建模是不现实的——纵横比太悬殊网格数量会爆炸而且裂缝厚度远小于储层尺度时数值求解精度反而变差。COMSOL处理这类问题非常成熟的做法是降维裂缝变成一维边界井筒变成一维线在三维模型里是线在二维模型里是点或线基质保持二维。裂缝所在的边界上应用裂缝流动接口Fracture Flow它自动把裂缝开度b、渗透率k_f和实际三维储层厚度H折算成等效导流能力K_f·b/H。我在算例里就是这么处理的储层是1000 m × 800 m的二维平面井筒水平穿过5条横向裂缝以100 m间距垂直于井筒排列。每条裂缝就是一条长度200 m的线段开度取2 mm。这样网格主体是二维三角形只在裂缝附近加密计算量基本可控又不丢失裂缝的导流贡献。计算的物理背景明白之后你会发现真正决定算例成败的其实是在正式建模前的几个判断题这几个方向选错后面全部白搭。2. 建模前需要想清楚的四个判断选错方向后患无穷2.1 只跑渗流还是必须耦合岩石变形做降压开采模拟最常被问到的问题是“要不要把固体力学加上”。我的建议是分两种情况如果只需要看产量、压力分布、累计采出量这一类生产动态指标单独用达西定律接口就够了效率高、收敛稳、参数容易解释。但如果要研究出砂风险、裂缝闭合对产量的反馈影响或者注入开发时可能发生的水力压裂延伸问题就必须做流固耦合把有效应力对渗透率的反馈写进模型。所谓有效应力用的是基本的Terzaghi/Biot表达式σ_eff σ_total - α p在降压开采过程中孔隙压力p下降有效应力增大裂缝会被压得更紧渗透率随之下降。这个反馈机制对气井产量影响极大尤其在生产后期裂缝闭合明显。算例里我选择了一种折中方案不在COMSOL里完整解固体力学场而是采用工程上常用的“渗透率应力敏感经验关系”给裂缝渗透率添加一个随压力变化的系数。这样做的好处是仍然保留裂缝闭合的物理后果但不用引入额外的固体力学未知量求解规模小很多。如果你确实需要土体位移、井周应变这些结果再叠加固体力学接口方程上就是把达西方程和动量平衡方程联立求解。2.2 裂缝为什么要写成低维域而不是靠加密网格我曾经见过有人把裂缝区域画成一个很窄的矩形宽度取0.002 m长度200 m试图直接在这个薄条里划分网格。结果就是网格单元长宽比达到十万比一质量差到COMSOL直接警告即使强行求解裂缝两侧的压力梯度也无法被准确解析。低维域的核心逻辑是流动在裂缝宽度方向上的梯度不重要主要流动方向是沿裂缝延伸方向的。因此把宽度信息折算进等效渗透率用一维域的线来描述它既降低自由度又天然避开了长宽比问题。COMSOL里实现裂缝流动时会要求你填一个“厚度”参数这个厚度就是裂缝开度b。软件内部会自动把高度方向的流动效应折算进“单位面积导流系数”。如果你做的是三维模型裂缝是一个二维面厚度b写在截面的属性里即可。2.3 降压条件怎么给定压生产其实是最省事的选择降压开采在现场有两种实现方式一种是把井口压力控制在某个值让井底流压随产量动态调整另一种是固定井底流压看产量如何变化。对数值模型来说固定井底压力constant pressure drawdown是最简单也最实用的边界条件直接在裂缝汇入井筒的位置设定压力值即可。但也有一个细节新手容易犯t0时刻整个储层是均匀的初始压力20 MPa井筒边界突然设成5 MPa压力阶跃就产生了。达西方程本身对压力是连续的但这个突跳会在起步阶段产生数值振荡尤其是在裂缝导流能力很高的位置。解决的办法通常有三种一是把井底压力变化用斜坡函数ramp function在几分钟内平滑降到目标值二是把初始求解的步长压得极小三是求解器里把BDF方法的代数阶限制为1或2。这三个办法我后面调试章节会再展开讲。我在算例里就采用了斜坡过渡的方式过渡时间取0.1天既贴近真实开井过程也回避了刚开始就强行射流那种紧约束带来的收敛困难。2.4 单位体系混乱是最隐蔽的坑COMSOL允许你在输入参数时直接用带单位的数值比如渗透率写成1e-16[m^2]压力写成20[MPa]比较省事。但麻烦在于裂缝渗流接口里的“开度”和使用者自定义表达式之间的单位换算。打个比方如果你写裂缝渗透率时直接填10[darcy]COMSOL会抱怨。更常见的问题是把达西Darcy和毫达西mD弄混1darcy≈1e-12 m²1mD≈1e-15 m²。页岩基质0.0001 mD就是1e-19 m²差之毫厘谬以千里。我建议全程使用国际单位制但把产气量显示单位设置成m³/天或MMscf/天这样既保证内部计算不出错出图也符合工程习惯。版本方面我用COMSOL 6.1一直到6.4跑过这套算例底层物理接口设置基本一致只是6.3之后求解器默认选项有所变化。建议直接用6.4自动生成的网格控制在内存占用上更友好。判断做完接下来就是真正上手搭建模型了。下面这部分我会直接给出参数表以及每一步在COMSOL界面的操作路径照着做基本能跑通第一版。3. 算例参数与COMSOL实操步骤3.1 基础参数表一份可以直接抄作业的清单下表是一份典型的页岩气压裂井降压开采参数基于典型现场数据的量级你完全可以根据自己手头的地质资料替换数值。参数数值说明储层长度 L1000 m矩形区域x方向储层宽度 W800 m矩形区域y方向储层厚度 H30 m用于体积流量折算基质渗透率 k_m1e-19 m²约0.0001 mD基质孔隙度 φ_m0.08/裂缝渗透率 k_f1e-11 m²约10 darcy量级裂缝开度 b0.002 m垂直于流动方向裂缝间距100 m/裂缝半长100 m裂缝全长200 m左右各延伸100 m井底流压 p_wf5 MPa降压目标值初始地层压力 p_i20 MPa恒定初始条件流体黏度 μ2e-5 Pa·s简化取常数也可用干气黏度经验公式井筒位置y400 mx∈[150,850] m水平井沿x方向参数量级上稍微解释一下基质渗透率取1e-19 m²确实低到离谱但这是页岩基质的常态裂缝渗透率取1e-11 m²10 darcy看起来高其实压裂砂充填的主裂缝普遍在这个范围。这里把裂缝开度2 mm、渗透率1e-11 m²折算一下裂缝导流能力k_f·b2e-14 m²·m换算成常用的darcy·cm单位大约是2 darcy·cm欠压裂设计中目标导流能力普遍在1-5 darcy·cm之间说得通。3.2 几何建模和物理场设置路径打开COMSOL模型向导里选择二维空间维度物理场添加“达西定律”darcy研究选择“瞬态”。几何方面我用矩形表示储层井筒就是一条水平线段从x150 m延伸到x850 m用“线段”工具直接画即可5条裂缝垂直线段依次放在x200、300、400、500、600 m上与井筒垂直。这里有个几何层面容易出错的细节裂缝线段必须和井筒线段“相交”才能保证裂缝内的流体能汇入井筒。如果两线只是靠近但没有交到一起COMSOL需要额外设置“内部边界对”或“装配”中的连续性条件否则裂缝里的流体流到井筒附近就断头了。如果能在画几何绿地时把每条裂缝端点直接捕捉到井筒中心点上后续就省去手动耦合的设置非常省心。物理场设置里整个二维域选择达西定律裂缝线段需要在它上面再添加一个“裂缝”特征流体流动裂缝指定开度b、孔隙度和渗透率k_f。井筒本身一般不单独设置达西定律因为井筒是恒压边界直接在5个裂缝与井筒的交点处设定压力即可。边界上储层四周默认是零通量边界表示封闭区块弹性产状——如果没有补充能量这对应衰竭式开发。3.3 网格划分先粗后细裂缝和近井区单独处理达西问题本身对网格不像CFD那么苛刻但裂缝和井筒附近的压力梯度比较大网格太粗会把裂缝的导流作用平滑掉。我具体的做法是全局最大单元尺寸定50 m最大单元增长率1.2为矩形储层设定一个“尺寸”节点限定内部单元不粗于50 m再单独为5条裂缝线段设置“边尺寸”约束最大单元尺寸5 m井筒线段设置最大单元尺寸10 m。这样生成后裂缝附近自然形成局部加密带远离裂缝的区域网格稀疏总自由度数大概在2万到4万之间普通笔记本10分钟以内能算完10年生产史。3.4 求解器关键设置瞬态研究的默认求解器是BDF。默认的“自由时间步长”策略会把步长拉得很大而降压开发初期压力变化剧烈后期变化缓慢自适应时间步其实够用。你只需要在“时间步进”设置里把最大步长限制为3个月约0.25年防止COMSOL在后期直接把步长跳到一年以上导致产量曲线变形。求解器相对容差建议设成0.01页岩气瞬态生产这种问题不需要更严格0.001会明显拖慢速度但对于从事研究出图需求可以选0.001求稳。在算例里求解时间我设置了0 d、10 d、30 d、100 d、1 yr、5 yr、10 yr 这7个输出时间点。实际COMSOL内部会自动在中间插入很多计算步这7个点只是输出查看用的“存储时间”。模型跑通之后真正的收获在看结果这一步。同样的云图有人只看到漂亮的色带有人能看到产量递减的物理机制。这部分我重点讲三条应该盯住的曲线和一个验证思路。4. 结果里最该盯住的三个现象与验证思路4.1 产量递减曲线不是一条简单的指数衰减定压生产条件下井筒附近裂缝网络的产气量会先经历一个短暂的高峰然后迅速下降进入长尾递减。如果画出双对数坐标产量和时间的关系在早期应该接近1/√t的斜线这是裂缝性储层典型的线性流特征后期当压力波碰到相邻裂缝的影响半径或者储层边界后曲线会进一步变陡进入边界控制流阶段。这套算例跑出来我会建议你重点看两个量一是“井筒处裂缝累计质量流率”COMSOL里的边界积分结果二是“储层压力平均值随时间的下降”。后者可以帮助判断储层范围是否足够大如果十年后整个1000 m×800 m区域的平均压力都降到了10 MPa以下说明储层已经被“掏空”了一大部分模型外边界效应已经介入如果平均压力还接近20 MPa说明模型对这个产量目标来说过大关心近井响应即可。4.2 压降漏斗在裂缝和基质之间的巨大差异运行结束后在后处理里选择“压力”切片图逐时间步看压力场形状。你会非常直观地看到裂缝周围的压力很快降低形成细长条状的泄压带而基质内部还保持着接近原始压力的状态。这个对比是压裂井模拟最有说服力的可视化结果它直观说明了为什么体积压裂改造能提高产量——裂缝把泄压面积从井筒周围的几十米扩大到数百米的线性区域。同时可以计算每条裂缝的流量占比。一般中间裂缝和两端裂缝的贡献有差异这种差异来自裂缝间的“干扰”中间裂缝左右两侧都有相邻裂缝在竞相泄压压力波遇到相邻裂缝后泄压区域收缩贡献下降两端裂缝一侧没有邻居贡献相对大。对这种非均匀流量分配的理解是后面做压裂参数优化的基础。4.3 用Theis解和物质平衡做双重验证数值模拟最容易让人担心的就是“算得对不对”。COMSOL算完之后至少要过两道验证第一道是大致验证早期压力传播。把问题近似成无限大均匀储层的径向流利用Theis解析解可以估算距离井r处、生产时间t时的压降幅度。我一般取井筒附近10 m的位置对比数值解和解析解在1天、10天、100天三个时刻的压力值误差在10%以内就算合理。当然由于裂缝存在这个对比不能指望精确匹配异常偏离才说明网格或边界条件有问题。第二道是物质平衡检查。累计产气量应该等于储层初始可采气体量减去当前剩余气体量。可以做一次“积分”在全局计算里定义累计产量表达式同时在模型里积分当前压力对应的气体存储量看两者是否满足质量守恒。COMSOL求解器是稳定的如果你发现质量不平衡超过5%大概率是裂缝开度设置或者边界条件写错了。结果端没问题之后算例真正见功力的是跑“歪”之后的调试。这里我把实测中踩过的几个坑按顺序列出来照方抓药基本都能救回来。5. 调试教训启动振荡、边界反射与裂缝导流衰减的三个坑5.1 开井瞬间的压力阶跃为什么会把求解器逼疯第一种典型症状是t0附近出现明显的振荡锯齿或者求解器报“不收敛”直接卡死。根因就是井底压力从20 MPa瞬间跌到5 MPa裂缝里的达西速度瞬间冲到一个不合理的值。我调试这套算例时一开始也以为是自己裂缝渗透率给太高了试过下调到1e-12 m²振荡依然存在。后来想到给井底压力加了一个“阶跃函数变体”p(t) p_i - (p_i - p_wf) * min(1, t/ t_ramp)t_ramp取0.1天到1天。改成压力斜坡过渡后呕吐式的振荡立刻消失。只要你不是在研究“瞬间开井”的冲击问题都建议用斜坡过渡这样既不会人为改变稳态油气生产特征又显著提高求解稳定性和速度。5.2 裂缝导流能力衰减几何上改裂缝还是物性上改渗透率降开采后期有效应力增加裂缝会被压实导流能力下降。这个效应不建模算出来的十年产量往往偏高。我建议的处理方法是为裂缝渗透率定义一个依赖于局部压力p的表达式k_f(p) k_f0 * exp( -β * (p_i - p) / p_i )β是无量纲应力敏感系数典型值取0.3-0.5。写成这样以后COMSOL会在每个迭代步自动计算当前压力下的局部渗透率等效实现了“降压-有效应力增加-裂缝闭合-渗透率下降”的物理链条。需要特别提醒的是如果你用移动网格Deformed Geometry去“真实”地改裂缝几何开度同时又在达西定律里给了一个动态渗透率两者就会互相重复计算结果会不符。所以二选一轻量化研究用渗透率表达式中网格动网格变形这种交作业式的展示一定要关掉物性表达式防止双重影响。5.3 边界反射导致的假递减“假递减”是我见过最坑人的一个。模型区域取小了压力波在几年内就传到了外边界然后被零通量边界“弹”回来产量曲线会在某个时间点上突然加速变陡看起来像气藏快速枯竭。实际上现场储层边界往往没有这么近。排查方法很简单后处理里画一条远离井筒剖面的压力线如果某个时间点起最外缘的压力开始发生明显变化就要警惕边界效应进入了关注区。补救办法是扩大模型范围或者改用“无限大储层”的近似边界条件边界提取设置在足够远的位置。我对这套算例的存储区域1000 m×800 m其实是做过的“加宽”如果你把范围缩小一半跑同样参数下产量差异会非常显眼。坑填完该说这个算例能怎么扩展了。说实话一个算例的价值不在于它能出那一张云图而在于改建和扩展的灵活度。6. 在算例基础上做扩展批量扫描、三维化与向压裂过程延伸6.1 用MATLAB/Livelink批量跑参数扫描COMSOL的批处理能力很多人还没用起来。比如想扫描裂缝间距对累计产气量的影响完全可以在COMSOL Desktop里定义参数化扫描但更灵活的是用Livelink for MATLAB在MATLAB里写循环边改压强、裂缝半长参数边调用COMSOL求解。脚本思路很简单model mphopen(fracture_well.mph); for spacing [80, 100, 120, 150] model.param.set(spacing, spacing); model.study(std1).run(); res mphglobal(model, intop1(flux)); results(spacing) res; end这样做的好处是数据库、后处理曲线可以全部集成到MATLAB里适合做“参数敏感性图”“产量预测概率曲线”这类科研图表。如果你更习惯Java或PythonCOMSOL 6.x同样支持Java API和Python客户端只是MATLAB生态相对完整文献也多。6.2 往三维裂缝网络和压裂过程延伸这套算例是二维、5条主裂缝的简化版本。如果你手头有力学性质数据丰富可以把几何扩展到三维储层变成一个六面体裂缝变成二维面井筒变成三维线。三维模型自由度会上一个台阶网格需要几十万到上百万单元求解十年前史对工作站也有压力。想再往前走一步把水力压裂过程中裂缝的延伸与扩展也纳入模拟就没有必要继续用这个算例里的简化模型了——那种问题更适合用离散裂隙网络DFN方法或者岩土工程常用的颗粒流程序来模拟COMSOL也可以扩展XFEM但那是另一个量级的复杂度了。我的实际体会是这套算例最值得借鉴的部分不是某一项高级设置而是“如何把工程问题翻译成多物理场模型的思维链路”。先物理图景清楚再决定降维方式和边界条件最后才是操作COMSOL的细节。任何一项的此轻彼重都会在计算结果上留下痕迹。所以你拿到这份算例建议先别急着改参数先严格按照默认跑一遍然后逐项调整裂缝间距和渗透率看结果怎么变再动手“造”你自己的井。