1. 煤层瓦斯抽采数值模拟为什么会走到三场耦合这一步做煤矿瓦斯治理的同行应该都有体会早期的抽采设计基本靠经验公式加现场实测。钻孔怎么布置、抽采半径取多少大多参照《煤矿安全规程》里的推荐值再根据本矿的瓦斯含量和透气性系数做些修正。这套办法能用但有个很尴尬的问题——现场条件稍微复杂一点比如煤层赋存深度变化大、地应力场不对称、或者抽采负压导致煤体变形明显经验值就开始漂了。我接触过不少矿井按经验半径布孔结果抽采达标时间比设计晚了三分之一后期要靠加密钻孔来补救。p这也是为什么近些年三维煤层瓦斯渗流与传热耦合的数值模拟越来越受重视。核心诉求很直接把瓦斯在煤层里的流动过程、煤体在应力作用下的变形过程、以及温度变化带来的热效应放到同一个三维模型里联合求解让抽采半径、应力、位移这些原本靠估的参数能从一个共同的物理模型里算出来而不是各算各的最后靠人为拼凑。p凡是在COMSOL里做过这类模型的同行应该都清楚真正的难点不在于软件操作而在于方程要不要自己写、怎么写、写完了怎么跟固体力学模块稳定地耦合在一起。COMSOL自带的Darcy模块能用但做瓦斯抽采研究时往往不够用——需要考虑Klinkenberg效应、吸附膨胀引起的渗透率动态变化、瓦斯解吸放热带来的温度场扰动。这些效应自带的模块没法直接表达必须自定义渗流方程。这篇文章就把我实际把自定义渗流方程和固体力学方程在一个三维模型里联合求解的完整思路写出来包括方程改写的细节、耦合操作链路、以及后处理时怎么从应力场和位移场里提取抽采半径的判据。p适合谁来参考一是做瓦斯抽采设计、需要论证钻孔布置方案的工程技术人员二是做煤层气开发、热流固耦合研究的科研人员三是正在用COMSOL做多物理场耦合、卡在自定义PDE和固体力学耦合这一步的研究生。文里涉及的操作以COMSOL Multiphysics为背景但方程改写和耦合思路换到其他有限元软件同样成立。p2. 从单场模型到三场耦合驱动逻辑是什么2.1 单场模型为什么算不准抽采半径先看最简单的单场模型把煤体当成刚性骨架瓦斯流动用达西定律描述渗透率取实验室测的常数。这种模型算出来的压力分布很平滑抽采半径呈规则椭圆扩展看起来挺像回事。但对照井下实测数据就会发现模型预测的抽采影响范围普遍偏大特别是抽采中后期实测瓦斯含量下降曲线往往比模拟滞后不少。p原因其实不复杂瓦斯抽采不是单纯的气体流动而是一个瓦斯解吸-渗流-煤体变形的连锁过程。瓦斯被抽走后煤层压力下降有效应力升高煤体被压缩裂隙闭合同时瓦斯解吸导致煤基质收缩裂隙张开再加上抽采造成的温度变化三方面共同影响渗透率。实测数据里那个滞后很大程度上就是应力敏感效应在起作用。你抽得越猛压力降得越快煤体压缩越明显渗透率下降得越快后续瓦斯就越难抽出来。刚性模型完全忽略了这个负反馈自然高估抽采效率。p2.2 应力场、渗流场、温度场之间的三条耦合路径要建立真正有用的三维模型至少要把三条耦合路径放进方程体系里p第一条应力对渗流的影响——煤体变形改变孔隙度和裂隙开度渗透率跟着变。这是三场耦合里最核心的一条工程上常用负指数形式的应力敏感模型来描述渗透率随有效应力增大呈指数下降。p第二条渗流对应力的影响——瓦斯压力本身就是孔隙压力直接参与有效应力计算。抽采过程中孔隙压力下降有效应力升高煤体发生压缩变形。这条路径在固体力学模块里通过有效应力原理天然耦合进去也是应力、位移计算结果的主要驱动因素。p第三条温度对两者的影响——瓦斯解吸是吸热过程抽采时煤体局部温度下降温度变化一方面改变瓦斯气体的黏度和吸附平衡另一方面产生热应力叠加到煤体变形里。反过来煤体变形和气体流动也会改变热量传递路径。这三条路径交织在一起就是典型的**热-流-固三场耦合THM耦合**问题。p2.3 自己写渗流方程而不是用自带模块差别在哪这时候就能看出为什么必须自定义渗流方程了。COMSOL自带的Darcy模块使用的是固定孔隙度和固定渗透率或者允许你用表达式驱动渗透率但吸附膨胀、Klinkenberg效应这类非线性项要反复改写源项操作很不顺手。更关键的是做科研写论文时审稿人一定会问你的控制方程是什么边界条件怎么给的参数怎么标定的用自带模块黑盒计算这些问题很难回答清楚。p自定义渗流方程的实际做法是把质量守恒方程自己写出来明确每一项的物理意义。通常采用气体压力作为求解变量方程写为p其中第一项是瓦斯累积量的变化率第二项是对流项右边是源汇项。把孔隙度、渗透率、气体密度都写成压力、应力和温度的函数后这个方程就同时包含了应力耦合和温度耦合。在COMSOL里用General Form PDE或Coefficient Form PDE模块把方程系数一项一项填进去控制权完全在自己手里。p3. 三维模型怎么搭几何、材料参数与定解条件3.1 几何模型的选择策略三维煤层瓦斯抽采模型几何上不需要做得很花哨。我常用的做法是建立长方体煤体区域尺寸根据钻孔间距来定走向长度取相邻两个抽采钻孔间距的一半利用对称性倾向长度取钻孔控制范围高度为煤层厚度。利用对称面做半模型或四分之一模型能显著降低网格量和计算时间。p钻孔的处理有两种思路。第一种是几何简化把抽采钻孔处理成圆柱形空腔钻孔壁面设置为定压边界负压值如-20 kPa表压换算成绝对压力约81 kPa。第二种是源汇等效不建钻孔几何直接在钻孔位置施加一个源项。第一种更直观后处理时可以看钻孔周围的应力集中和渗流形态第二种胜在网格友好、收敛性好。我在做抽采半径系统分析时更倾向于第一种因为要同时观察钻孔周围的应力、位移分布特征几何简化更有说服力。p3.2 渗透率动态演化模型的选择渗透率模型是整个自定义渗流方程里最敏感的一个环节。文献里常见的模型有Palmer-Mansoori模型、Shi-Durucan模型、以及各种形式的应力敏感指数模型。我的建议是不要盲目追求复杂优先选参数能通过现场数据标定的模型。p我自己用下来最顺手的组合是这样渗透率随有效应力的变化用指数形式吸附膨胀引起的渗透率变化用附加项表达Klinkenberg效应在低压力段单独修正。整体渗透率写成pk k0·exp(-α·Δσ_eff) Δk_sw其中Δσ_eff是有效应力变化量Δk_sw是吸附膨胀引起的渗透率增量。α叫应力敏感系数MPa⁻¹量级一般通过实验室三轴渗流实验拟合Δk_sw通过煤的吸附膨胀实验测定。把这些都写成压力和应力的显式函数后在COMSOL里以变量Variables形式全局定义方程引用时直接调用即可。p3.3 边界条件与初始条件怎么给才不飘边界条件的选择直接影响结果能不能站得住脚。我习惯这样设置p渗流场煤体四周外边界设为无流动边界钻孔壁面设为定压边界取抽采负压对应的绝对压力煤层顶底板界面按实际情况处理如果顶底板透气性差设为无流动边界如果顶底板是采动裂隙发育区则设为定压边界。p固体力学场底部固定四侧施加水平应力约束大小按地应力实测值折算顶部施加垂直应力模拟上覆岩层自重。钻孔壁面为自由面不额外加压。p温度场初始温度按地温梯度计算煤体外边界设为恒温钻孔壁面与抽采气体之间考虑对流换热。p初始条件方面渗流场初始压力取实测煤层瓦斯压力换算成绝对压力应力场初始状态建议先做一次地应力平衡计算把重力荷载和构造应力加载进去得到初始应力场和初始位移场应接近零再以此为初值开始抽采模拟。很多新手漏掉这一步结果位移场里出现离谱的整体刚体位移后面应力分析全乱套。p4. COMSOL里的实际操作自定义方程与固体力学耦合的完整链路4.1 选择哪个PDE接口合适COMSOL里自定义渗流方程可以用Coefficient Form PDE系数型偏微分方程或General Form PDE一般型偏微分方程。两者的差别在于系数型PDE适合能整理成标准形式da·∂u/∂t ∇·(-c∇u - αu γ) β·∇u au f的方程一般型PDE更自由允许任意形式的通量项和源项。p对于瓦斯渗流方程我推荐用Coefficient Form PDE省内存、收敛性好而且COMSOL对系数型的求解器优化做得更好。关键是把方程整理成标准形式。例如简化的瓦斯渗流方程可以整理成pda φ·μ/(k·p)对应累积项系数 c p/μ对应扩散系数 f Q_m源汇项p实际写的时候da、c、f都是表达式里面可以引用全局定义的变量渗透率k、孔隙度φ、黏度μ等都是空间坐标和时间t的函数。这样写的好处是方程形式一目了然检查和调试都方便。p4.2 固体力学模块的关键设置固体力学模块用COMSOL内置的Solid Mechanics即可。需要额外处理的是初始应力场的施加和孔隙压力的传递。p孔隙压力传递是最关键的耦合操作在固体力学模块中把瓦斯压力作为孔隙压力加入有效应力计算。具体操作是在固体力学接口的实体设置里勾选孔隙弹性或直接在体载荷里施加孔隙压力项载荷大小引用渗流场的求解变量p。这样渗流方程每算一步的压力分布都会实时更新到固体力学方程的载荷项里。p初始应力场我建议用两次计算的策略第一次计算只加载重力和边界力不激活渗流场得到初始应力状态然后把第一次计算的应力场结果作为第二次计算完整耦合计算的初始值。这样能避免初始位移过大导致的收敛问题也符合地下工程数值模拟的一般做法。p4.3 热场模块的自定义与耦合热场可以直接用COMSOL的Heat Transfer in Solids模块但要把瓦斯解吸吸热和对流换热两项处理清楚。解吸吸热是负的源项和单位时间解吸量成正比。单位时间解吸量怎么算需要从渗流方程里的累积项变化率推算——在自定义渗流方程中累积项通常是压力、温度的函数它对时间的偏导数就是单位时间的瓦斯解吸速率。把这个量提取出来乘上解吸热一般取40-60 kJ/mol作为热场方程的负源项施加。p对流换热项则通过渗流速度矢量传递把达西速度u_darcy计算出来在热场方程里添加对流项。COMSOL的Heat Transfer in Solids默认只有传导需要手动添加对流项或使用Heat Transfer in Porous Media接口。我个人的经验是如果瓦斯流速很小达西速度在10⁻⁶ m/s量级对流换热量通常可以忽略温度场主要受解吸吸热控制。这个判断对计算收敛性影响很大——加了强对流项之后非线性程度上升收敛难度明显增加。p4.4 求解器设置与时间步长控制三维THM耦合模型最大的实际困难是收敛性。我自己总结几个行之有效的做法p第一先做稳态地应力平衡再做瞬态抽采。不要一上来就全耦合瞬态求解那样很少能收敛。p第二时间步长用对数增长。抽采刚开始时压力梯度大、变化快时间步取1秒甚至0.1秒后期压力场趋于平稳时间步可以放到1天甚至10天。COMSOL里直接设置时间序列比如range(0,1,100) range(100,10,1000) range(1000,100,10000)这种梯度序列。p第三开启自适应网格细化特别是钻孔壁面附近。瓦斯压力梯度在钻孔周围最剧烈网格不够细时压力场会出现非物理振荡。我在钻孔壁面处设置边界层网格第一层厚度取钻孔半径的1/20层数5-6层效果明显。p第四使用分离式求解器。渗流场和固体力学场一次全耦合求解牛顿迭代矩阵规模大很容易不收敛。改用分离式求解器Segregated Solver先求渗流场再求应力场反复迭代到收敛稳定性好很多代价是每个时间步的迭代次数增加了。实际算下来分离式求解的总时间未必更长但很少因为不收敛而中途崩溃。p5. 抽采半径、应力重分布与位移场后处理怎么提取判据5.1 抽采半径的标准定义与三维可视化抽采半径的国标定义主要依据《煤矿瓦斯抽采达标暂行规定》抽采后煤层残余瓦斯压力降到0.74 MPa以下或残余瓦斯含量降到8 m³/t以下的范围。模拟时一般用压力阈值判据计算结束后提取钻孔周围压力低于0.74 MPa绝对压力的区域边界该边界到钻孔中心的距离就是抽采半径。p三维模型里做这件事有几种方法。最简单的是用COMSOL的Cut Plane功能在煤层中切出中面或不同高度截面画等值线图找到0.74 MPa等值线量距离。更精确的可以用Volume Selection提取整个低压区看它的三维形态。实际抽采钻孔是空间弯曲的煤矿井下定向钻孔低压区形态并不是标准圆柱对称这时候三维可视化就非常有必要——它能清楚地看出抽采半径在走向、倾向、垂向三个方向上的差异。p5.2 应力重分布的三个典型区域抽采过程中煤体应力发生重分布三维模型后处理时我一般把钻孔周围分成三个典型区域来看p卸压区紧邻钻孔壁面瓦斯压力降幅大有效应力升高但总应力卸除煤体发生向钻孔方向的位移。这个区域扩容裂隙发育是瓦斯流动的主要通道。p应力集中区在卸压区外围由于卸压区煤体向钻孔方向收缩外围煤体承担了额外的载荷出现应力集中。应力集中系数和钻孔直径、抽采负压、地应力大小都有关。三维模型算出来的应力集中区往往呈环形围绕钻孔但形态受层理面和原始裂隙影响而变得不规则。p原岩应力区距钻孔足够远的地方应力不受抽采扰动影响保持原始地应力状态。p判断这三个区域的边界最直观的指标是最大主应力分布和位移矢量场。应力集中区的最大主应力明显高于原岩应力卸压区的位移矢量指向钻孔中心应力集中区的位移矢量指向卸压区。p5.3 位移场怎么用从位移反推卸压圈位移场是很多人算了但不用的一项输出实际上它的信息量比应力场更丰富。有两个典型用途p用途一评估卸压圈范围。煤体向钻孔方向的径向位移达到一定阈值比如5 mm的区域基本对应裂隙充分发育的卸压圈。这个卸压圈大小对判断二次增透措施如水力压裂、CO₂致裂的布孔位置很有参考价值。p用途二判断抽采对巷道稳定性的影响。如果钻孔距离巷道太近抽采引起的煤体收缩位移可能造成巷道变形。三维模型算出巷道周边的位移分布后可以提前评估是否需要加强支护。p5.4 不同抽采负压和布孔方案下的对比分析模型标定好之后最有工程价值的应用是做方案对比改变抽采负压、钻孔直径、钻孔间距观察抽采半径和应力位移场的变化规律。我做过一个具体案例在相同地质条件下把抽采负压从20 kPa提高到40 kPa抽采半径在初期增长明显因为压力梯度大了瓦斯流速加快但90天后的抽采半径增幅却有限——原因是高负压导致煤体压缩加剧渗透率下降抵消了一部分负压增大的增益。这个结论在刚性模型里永远算不出来只有三场耦合模型才能捕捉到。p这种负压不是越高越好的分析结果在工程设计中非常有用。实际做方案时可以做成参数化扫描负压、孔径、孔间距三个变量各取几个水平组合计算最后输出一个抽采效果矩阵供设计人员直接查表选取。p6. 调试与验证三维耦合计算最常见的几个坑6.1 收敛性卡死的定位思路三维THM耦合模型算到一半突然不收敛是每个人都会遇到的事情。我的排查顺序是固定的p第一步看发散发生在哪些物理场。分离式求解器会把每个物理场的迭代残差单独显示先看是谁先发散的。如果是渗流场发散检查时间步长是否过大、压力边界是否合理如果是固体力学场发散检查材料参数是否出现了负刚度、网格是否严重畸变。p第二步看发散的部位。用结果图滚动显示最后几步的分布通常发散点就是问题点。我遇到过好几次发散点都在钻孔底部尖端处——那里几何尖锐应力奇异网格稍粗就过不去。解决办法是把钻孔底部做成圆角或者局部加密网格。p第三步看参数是否越界。自定义方程里如果有分母项比如渗透率作为分母检查压力或应力是否让分母趋近于零。COMSOL里可以在变量定义里加约束函数如flc2hs平滑截断防止参数越界。p6.2 渗透率模型的参数标定陷阱渗透率动态演化模型里应力敏感系数α是一个极其敏感的参数。α差一个数量级抽采半径预测结果可以差到一倍以上。但α又不是一个纯物理参数它受煤阶、裂隙发育程度、围压条件影响很大。我踩过的坑是直接从文献里拿一个α值来用结果模型算出来的抽采半径比实测小得多。p我的建议是用现场抽采数据反演标定α。具体做法先用一组钻孔的实测瓦斯压力下降数据作为目标把α作为待定参数手动或借助优化模块反复计算直到模型的压力下降曲线和实测曲线拟合上。这个过程比较耗时每一组参数都要跑一次完整的三维计算但值得做——标定完的模型才能用于方案预测。p6.3 温度场耦合要不要每次都加温度场耦合是很耗计算资源的如果研究问题不涉及明显的温度效应可以适当简化。我的判断依据是如果煤层瓦斯含量低、抽采时间短比如预抽时间小于90天解吸吸热引起的温度降幅一般只有几摄氏度对瓦斯流动的影响可以忽略这时可以关掉温度场只做流-固两场耦合计算速度能快50%以上。p但如果做的是深部煤层地温高、低透气性煤层抽采周期长、或者研究热激励增透这类方法温度效应就必须完整考虑。深部煤层原始温度高瓦斯解吸吸热比例也高温度变化可能达到10℃以上这时候热应力、气体黏度变化都不敢忽略。p6.4 网格无关性验证真的不能省三维模型计算量大大家都想用粗网格快速出结果但网格无关性验证不能省否则你算出来的抽采半径到底是物理规律还是网格效应自己心里都没底。我常用的办法取三套网格——粗网格钻孔壁面边界层6层、中等网格10层、细网格15层网格数大概差2-3倍算同一个工况对比钻孔中心压力下降曲线和抽采半径。如果细网格和中等网格的结果偏差小于5%基本可以认为中等网格够用。p实际经验是钻孔壁面附近的压力梯度和应力梯度最大那里的网格密度决定结果准确性远场网格影响不大。所以重点是边界层的厚度和层数而不是整体加密。p7. 从模型结果到工程设计还有一段路要走模型算完之后数据怎么落到实际工程里是另一个阶段的活。我通常把模拟结果输出成三个层面的东西p设计参数层面——抽采半径随抽采时间的变化曲线。这条曲线可以直接用于钻孔间距设计两排钻孔间距取2倍抽采半径留一定安全系数抽采时间取达标所需的预抽期。p方案对比层面——不同布孔方式下的达标时间对比。比如顺层平行布孔和交叉布孔在相同孔数下哪个达标更快用模拟可以定量回答。p风险预警层面——把应力集中区和位移较大区域标注在采掘工程平面图上提醒现场注意这些区域的突出危险性和支护薄弱点。p我个人实际操作中比较深的一点体会三维耦合模型不是算完就结束了真正的问题往往是算得出来但用不上——因为你算的结果和你拿到的现场数据对不上。所以越是做实用型研究越要把模型标定放在第一位。先花时间把历史抽采数据整理清楚再调模型比起先算一堆花样结果再回头找数据效率高得多。抽采半径、应力分布、位移场这些常规输出经过标定的模型算一遍就够没有标定的模型算十遍也是在原地转圈。
