企业官网建设流程全解析

1. 先搞清楚实验台里发生的物理链条：这不是简单的“把CO2灌进去”

如果你做过煤层气或页岩气的室内岩芯驱替实验，一定熟悉这个场景：把一块致密岩芯夹持在哈氏合金模具里，加围压、抽真空、饱和甲烷，然后在恒温箱里以恒定流量注入CO2，盯着出口组分检测器看CH4含量什么时候掉下来。整个过程看起来很简单——一种气体进去，另一种气体出来。可真要把它放到COMSOL里算明白，很多人第一步就走歪了，因为他们把问题简化成了“两种气体在孔隙里的混合”。

实际上，实验室岩芯里发生的是一个由三条物理链叠加的事件。

第一条链是渗流。气体在页岩或煤岩的纳米-微米级孔隙里流动，压力梯度是驱动力，渗透率是门卫。因为流速极慢、雷诺数远小于1，惯性项可以完全忽略，这就是达西流动的适用区间。第二条链是传质。注入的CO2和原有的CH4之间有分子扩散，又有对流裹挟，两者叠加成了前缘弥散。第三条链才是这个课题的核心——吸附/解吸。煤岩和页岩有机质对CH4和CO2都存在较强的吸附能力，而CO2的吸附亲和力普遍高于CH4，所以CO2分子会主动抢占吸附位，把原本吸附态的CH4“挤”进自由孔隙，这才是驱替与增采的微观来源。

三条链的时间尺度差别极大：压力传到稳定可能只需要几十秒，对流突破可能要几小时，而吸附平衡的局部响应介乎两者之间。这种多尺度耦合，恰恰是COMSOL里做多物理场仿真的价值所在——它能同时求解压力场、浓度场、吸附量场，还能把渗透率随吸附量变化这种“反过来影响第一条链”的反馈也闭环起来。

我不建议上来就建一个花里胡哨的三维模型。对于实验室尺度的岩芯驱替，COMSOL里最务实的选择是从一维或二维轴对称开始。为什么？因为实验夹具的设计通常就是轴向流动占绝对主导，径向梯度只在注入口附近有局部影响。一维模型能让你先把物理逻辑跑通，把吸附参数、边界条件这些变量解耦出来，再去扩展维度。我见过太多人第一步就建三维网格，节点几十万，最后计算发散或算一星期没结果，问题根本不在几何，而在物理场耦合关系没理顺。

所以在动手建模之前，先问自己三个问题：我要复现的是驱替过程中的哪个阶段？我的实验能提供哪些边界条件（流量还是压力）？我关注的是出口突破曲线、累计采气量，还是岩芯内部的饱和度分布？这三个答案直接决定模型复杂度和物理场选择。

2. 物理场选型与方程组耦合：达西定律、多组分输运和吸附源项的写法

2.1 为什么主控方程是达西定律，而不是N-S方程

很多从流体力学转过来的人总想用Brinkman方程或者直接上N-S方程，这在实验室岩芯尺度是完全没必要的。岩芯渗透率通常在0.1 mD到10 mD量级，孔隙喉道直径在几十纳米到几微米，流动速度按达西定律估算下来在10⁻⁵~10⁻³ m/s量级。折算成雷诺数Re = ρu√k/μ，大概在10⁻⁶到10⁻³之间，惯性效应比粘性效应低了好几个数量级。这种条件下用N-S方程不仅浪费算力，还会引入数值假扩散，收敛难度直线上升。

COMSOL里就用“达西定律（dl）”物理场，方程是经典的：

∂(ερ)/∂t + ∇·(ρu) = Qm

其中流速项 u = -(k/μ)∇p。这里有两个关键点容易被忽视。

第一，气体密度不能写成常数。CH4和CO2在4 MPa、313 K条件下，密度已经偏离低压稀薄状态很远。达西定律模块里要勾选可压缩流动选项，把密度定义成压力相关的理想气体形式 ρ = pM/(RT)。如果你图省事用不可压缩流动，压力解出来看着没错，但浓度输运方程里的对流项会严重失真，突破时间能差出百分之几十。

第二，气体是混合物，粘度不是单一值。CO2的粘度在常温高压下比CH4高不少。严格做法是在模型里定义混合规则，比如COMSOL自带的质量分数平均粘度或者Hirschfelder近似。如果只做初步趋势研究，用常数粘度 2×10⁻⁵ Pa·s 我可以接受，但你要知道这是近似。

2.2 多组分输运：稀物质够用吗？别踩“浓度不可比”的坑

这是COMSOL建模里最容易埋雷的地方。甲烷和CO2在4 MPa、313 K下，各自浓度都在摩尔量级——按理想气体算，p/RT ≈ 4×10⁶/(8.314×313) ≈ 1500 mol/m³。注意，两组分浓度是同一个数量级，谁都不是谁的“痕量稀释物质”。

如果你用的是“多孔介质中的稀物质传递（tds）”，要清楚它的前提是溶质浓度远小于溶剂。CO2驱替CH4的后期，CO2浓度接近1，这时候稀物质假设直接崩了。解决方法有两种。

第一种方案：使用“浓物质传递（tdc）”物理场，COMSOL 6.x里这个模块对多组分气体混合物支持得比较完整，能定义二元和多元Fick扩散矩阵。代价是计算量和非线性程度都上去了。

第二种方案（我自己建模初期更常用）：用两个稀物质传递方程分别解浓度c_CH4和c_CO2，对流项共用同一个达西速度场，第三组分视为惰性平衡组分。只要控制模拟时间到主流突破前后，两组分的数值都在合理范围内，这个近似是可接受的。但要注意，扩散系数必须用有效扩散系数De，而不是分子扩散系数D。有效扩散系数要考虑孔隙度和迂曲度：

De = D·ε/τ

对于页岩或煤岩，迂曲度τ通常在3~5之间，多孔介质有效扩散系数会明显低于文献里查到的分子扩散系数。在4 MPa高压下，CH4-CO2二元分子扩散系数大约在1×10⁻⁶~3×10⁻⁶ m²/s量级，乘以孔隙度0.1再除以迂曲度4，De大概在10⁻⁸ m²/s量级。这个数值会直接影响Peclet数的判断和网格设计，我后面细说。

2.3 吸附源项：扩展Langmuir在COMSOL里的正确写法

吸附这个环节既是物理核心，也是数值陷阱。如果直接把吸附量当作平衡值代入浓度方程，会出现“瞬间吸附”假设，在驱替前缘会制造出一个剧烈的源项突变区，轻则收敛困难，重则解出负浓度。

更稳妥的做法是引入吸附动力学。用一个普通的常微分方程描述吸附量向平衡态的趋近：

dq_i/dt = k_i·(q_eq_i - q_i)

平衡吸附量用扩展Langmuir方程：

q_eq_CH4 = qm_CH4·b_CH4·p_CH4 / (1 + b_CH4·p_CH4 + b_CO2·p_CO2)

q_eq_CO2 = qm_CO2·b_CO2·p_CO2 / (1 + b_CH4·p_CH4 + b_CO2·p_CO2)

其中分压p_i = y_i·p = (c_i/c_tot)·p，c_tot = c_CH4 + c_CO2（如果按理想气体近似）。

COMSOL里实现这个有两种路径。一是直接在变量（Variables）节点写q_eq表达式，再用域常微分方程（Domain ODEs）的分布式ODE节点定义dq_i/dt=k_i·(q_eq_i-q_i)；二是在稀物质传递物理场的吸附子节点里自定义Langmuir型吸附，但那个子节点更适合单组分，多组分竞争吸附建议走手动ODE路线，可控性更强。

吸附动力学常数k_i怎么取？如果实验室的吸附等温线测试显示低压力下几分钟内就能达到平衡，那k_i取10⁻³~10⁻² s⁻¹量级合适。如果吸附速率很慢，比如受扩散控制，就要把k_i降到10⁻⁵~10⁻⁴ s⁻¹。这个参数对突破曲线形态影响极大，我建议做参数扫描对比，而不要拍脑袋定死。

2.4 三个物理场的耦合关系：在模型树里怎么组织

物理场之间的耦合并不复杂，但顺序要清楚。

• 达西定律输出压力场p，由p的空间导数生成达西速度u_darcy。
• 输运方程引用u_darcy作为对流速度，同时引用p来计算分压。
• 吸附ODE引用当地分压计算平衡吸附量和吸附速率。
• 吸附速率的累积反过来通过源项汇入输运方程：R_i = -ρ_b·(1-ε)·dq_i/dt。
• 如果还要考虑吸附引起的基质膨胀，那么吸附量还要反向影响达西定律里的渗透率k。

在COMSOL里，耦合变量直接写表达式就行，不需要手动指定“耦合链接”。但建议在“变量”节点里把单位对齐。浓度用mol/m³，压力用Pa，吸附量用mol/kg，骨架密度用kg/m³。单位错一个，源项就大错特错。我吃过一次亏：把吸附量单位写成了mol/m³，直接当源项用，结果计算出的产出量翻了好几倍，半天才查出是单位问题。

3. 几何、网格与边界条件：90%的收敛问题出在这里

3.1 一维模型还是二维轴对称？实验室几何怎么简化

实验室岩芯标准尺寸一般是直径25~50 mm、长度50~100 mm，更长的能达到300 mm。我建议起步阶段直接用一维模型：只沿长度方向建一条线，长度取岩芯实际长度。一维不仅能秒算，还能让你把吸附、传质这些物理逻辑彻底想清楚。等一维结果跟实验数据对上了，再扩展成二维轴对称模型研究注入口附近的径向效应。

二维轴对称模型有个好处，可以观察注入端附近的“指进”和径向浓度分布，但代价是非线性计算量大得多。如果只是看出口突破时间和累计采气量，一维精度已经足够。不要被“高维建模显得专业”这种想法绑架。

3.2 定压注入还是定流量注入？边界条件必须对应实验操作

这是边界条件设计里最实际的问题。

如果实验室用的是恒速泵（Syringe Pump）控制CO2流量，那么注入端边界条件就不要设固定压力。做法是设一个通量边界，流入的CO2摩尔通量等于泵流速换算而来的值。公式换算要留意：泵的体积流量Q_v（如1 mL/min）换算到岩芯端面面积A上的达西速度u = Q_v/A，再乘上摩尔浓度就是摩尔通量。换算完之后务必确认量纲。

如果实验是定压驱替——比如气瓶经减压阀稳定供给CO2、回压阀维持出口压力——那就设注入端压力p_inj，采出端压力p_out。COMSOL里注入端设压力约束，采出端也设压力约束，谁便宜谁就简单。

我见过不少人直接把入口设成固定浓度c_CO2 = 1。这在注入端是有物理依据的，但要注意：浓度边界会叠加在对通量格式上，如果同时设了流量边界，两种条件会互相矛盾，COMSOL会提示过约束。这里建议把问题拆分：入口要么给流量，要么给压力，浓度条件用“流入浓度等于注入气体浓度”这种方式放在通量边界里。

初始条件方面，岩芯先饱和CH4，那么初始压力p_init = 4 MPa，c_CH4由初始压力按理想气体算出来，c_CO2给一个很小的本底值比如0.001 mol/m³，避免数值除零。甲烷初始浓度不能给0，因为要和吸附方程里的分压关联。

3.3 Peclet数评估与网格时间步长匹配

这是影响求解稳定性的硬核问题。Peclet数定义为对流与扩散之比：

Pe = u_pore·L / De

其中u_pore = u_darcy/ε是孔隙实际流速。用我刚才给的参数估算：u_darcy = 1×10⁻⁵ m/s量级（定流量注入时），ε = 0.1，u_pore = 1×10⁻⁴ m/s，L = 0.3 m，De = 2×10⁻⁸ m²/s，算出来Pe = 1500。这是个非常大的数，意味着问题是对流主导的，前缘几乎是活塞式推进。

对流主导问题最怕两件事：网格太粗导致数值色散（前缘出现非物理振荡），时间步太长导致前缘跨过多个网格单元。COMSOL里稀物质传递模块有迎风稳定化选项，默认是开着的，但仅靠它是救不了粗网格的。

经验法则：前缘宽度至少要覆盖3~5个网格单元。驱替前缘的宽度跟纵向弥散系数相关，弥散系数综合了分子扩散和机械弥散。机械弥散系数通常近似为α·u_pore，其中α是纵向弥散度，实验室岩芯通常取0.001~0.01 m。代入之后D_disp = 0.005×1×10⁻⁴ = 5×10⁻⁷ m²/s，远大于分子扩散的De。这么算下来，前缘宽度量级在厘米级。300 mm长的岩芯，网格总数300~500个就够，再细分意义不大。

时间步长方面，建议限制CFL条件：每个时间步内前缘移动不超过一个网格单元。如果网格尺寸1 mm，u_pore = 1×10⁻⁴ m/s，那CFL限制的步长约10 s。COMSOL的自适应时间步长通常能处理，但如果你设了太BT的参数扫描，偶尔会被不连续点卡住。这时手动设初始步长0.1 s，最大步长60 s，往往能平稳跑完。

4. 求解器与参数扫描：非线性不收敛的排查链路

4.1 报错之前先检查这五个地方

“求解器未收敛”或“找不到一致的初始值”这类红色报错，90%的原因不在求解器本身，而在模型设置的某个角落。我按排查优先级列一下，你照着做能省大量时间。

第一，初始条件自洽性。压力初始值、浓度初始值、吸附量初始值三者必须同时代入方程后不产生巨大的初始残差。最常见的坑是把压力初始值设为某个值，但浓度初始值按另一个压力算的，导致分压和吸附平衡量在t=0就矛盾。解决方法是先跑一个稳态求解（只开达西定律，关闭输运和吸附），把压力场算出来之后，再以稳态解作为瞬态的初始条件。这一步几乎免费，但能消除绝大多数初始震荡。

第二，时间尺度不匹配。压力扩散的特征时间可能只有几十毫秒，而驱替过程是几小时。COMSOL在全耦合求解时会同时推进两个时间尺度，如果压力方程的存储项设置不当，会在初始几步产生极小的步长。解决思路有两个——要么在达西定律里忽略压力存储项做准稳态假设；要么用分离式求解器，先把压力场在每一时间步内快速收敛，再去更新浓度场。

第三，吸附平衡突变。扩展Langmuir表达式里，如果分母接近零（也就是所有分压都趋近零），会造成除零。在COMSOL里可以用平滑函数给分母加一个微小常数，比如1e-6 MPa，这种处理对结果影响可以忽略，但能让雅可比矩阵稳定下来。

第四，入口条件的阶跃突变。实验上你不可能瞬间把注入端从封闭切换到6 MPa，COMSOL里如果用阶跃函数，初始时刻就是一个强不连续，这会让全耦合求解器在第一轮就被打爆。正确做法是用平滑阶跃（cosh或双曲正切形式的过渡段），让注入压力或流量在1~10秒内缓慢升高。

第五，材料参数单位出错。COMSOL会在求解前做单位检查，但不会检查物理合理性。渗透率如果按SI写成了毫达西数值没换算（1 mD = 9.87×10⁻¹⁶ m²)，算出来的达西速度会对，但时间尺度会天差地别。这个错误我见过太多次。

4.2 参数扫描怎么设才不浪费时间

COMSOL的“参数化扫描（Parametric Sweep）”功能很强大，可以配合“辅助扫描”做二维参数组合。做这个课题时，我建议第一轮扫描固定其他参数，只扫注入流量：0.5、1、2、5 mL/min。你会发现突破时间随流量基本呈线性缩短，但注入流量太高时出口CH4曲线的“拖尾”现象会加重，这是因为局部CO2浓度高但吸附置换还没充分完成。

第二轮再扫Langmuir参数。重点是b_CO2对b_CH4的比率，这个是选择性的核心。实验室文献里CO2/CH4的吸附选择性通常在2~6之间。选择性越高，突破越晚、采出效率越高。把选择性这个参数扫一遍，你能直观看到“置换能力”对曲线的决定作用。

第三轮看渗透率的影响。我建议在参数扫描里专门设一档k0的基准值分别乘以0.1、1、10，然后画出累计采出量一系列曲线。你很快会发现：渗透率影响的是驱替过程的时间尺度，不影响最终的累计采出量（在吸附参数不变的前提下）。这能帮你区分哪些参数决定“时间”，哪些参数决定“总量”。

4.3 完整排查案例：连续三个小时“求解器未收敛”

我讲一个真实经历。当时设了定流量边界，在入口给了2 mL/min的CO2注入，初始条件设为压力4 MPa、甲烷分压4 MPa，打开瞬态就报“没有找到一致的初始值”。

一个排查动作就找到了根因：我设的流体条件里保留了达西定律的“存储项S_p”，默认值是1e-4 1/Pa，这对应于刚性骨架。但气体压缩性跟压力非线性强相关，在4 MPa下气体压缩系数远远超过骨架压缩系数，存储项就出现了几个数量级的动态变化，让初始雅可比矩阵条件数大到爆炸。

处理办法很简单：把达西定律改用“理想气体”可压缩选项，存储项由软件自动计算；同时把输运方程的求解器从“全耦合”改成“分离式”，让压力和浓度分步迭代。之后求解器很顺利地跑完了10小时的瞬态模拟。后来我把同样的做法迁移到另一台新装的COMSOL 6.4上，连Linux环境下的批处理也没出问题。

5. 进阶方向：渗透率动态变化与移动网格的坑

5.1 吸附膨胀/收缩怎么反馈到渗透率

搞CO2驱替CH4，尤其是煤岩，最不能忽略的反馈就是吸附膨胀。CO2吸附量升高会让基质膨胀，孔喉变窄，渗透率下降，反过来阻碍后续注入。这是实验室里能实际观测到的现象，也是模型的灵魂所在——不写这个反馈，模型就只是个“两种气体混合”的空壳。

最常用的渗透率变化模型是Shi-Durucan或Palmer-Mansoori型。Shi-Durucan的简化版本是这样：

k = k0·exp(3·c_p·Δp - 3·ε_s/φ0·Δq)

其中Δp是压力变化，Δq是总吸附量变化，ε_s是最大基质应变系数，φ0是初始孔隙度。在COMSOL里，把k写成达西定律材料属性里的表达式即可，比如k_var = k0exp(3cp*(p-pinit)-3es/phi0(q_CH4+q_CO2-q_CH4_init-q_CO2_init))。

这一步实现起来不难，难在参数标定。ε_s通常要从体积应变实验数据去拟合，没有实验数据时可暂取0.01~0.05。c_p（孔隙压缩系数）大致在10⁻⁶~10⁻⁵ Pa⁻¹量级，具体值与岩性关系密切。初次建模时建议把渗透率变化的幅度用参数扫描包住，先看渗透率衰减50%和衰减十倍的差异，再决定要不要做实验标定。

5.2 移动网格做岩芯变形：谨慎，不是所有版本都好用

COMSOL有“移动网格（Moving Mesh）”功能，很多人在做基质膨胀仿真时乐呵呵地打开这个物理场，指望岩芯几何能自动跟着吸附量变形。我的建议是：除非你研究的是宏观应变与应力耦合，否则别碰移动网格。

为什么？因为实验室岩芯被夹持在刚性哈氏合金模具里，径向变形几乎被完全约束，轴向变形也受围压和端部摩擦限制。几何变形量级在微米到数十微米级别，对宏观渗流路径的影响完全可以折算到渗透率表达式里。你打开移动网格之后，网格会被扭曲，极化到输运方程和达西定律里的雅可比矩阵，时间步长至少会缩小一个数量级，计算效率显著下降。

我前两年在COMSOL 6.2上试过一次带移动网格的全耦合煤岩模型，岩芯膨胀0.2 mm，网格扭曲集中在进口端2 cm区域，即便用了自动重划分，最终曲线跟固定网格版本相比差异不超过3%。从那以后我的默认方案就是“渗透率反馈”而不是“几何变形”。

如果你的研究目的确实涉及应力场，比如模拟水力压裂后的裂缝开度变化，再考虑引入固体力学模块。否则，请把我的经验记下来：移动网格是锦上添花，不是雪中送炭。

6. 把模型和实验数据对齐：参数标定与关键输出的定义

6.1 从纯组分吸附等温线到二元扩展Langmuir的参数整理

文献中查来的吸附参数往往是纯组分等温线拟合结果，比如纯CH4的Langmuir体积V_L和Langmuir压力p_L。转成COMSOL模型里扩展Langmuir的参数，要小心映射关系。纯组分Langmuir形式是：

q = qm·b·p/(1+b·p) 或者写作 q = V_L·p/(p_L+p)

两种表达式的参数换算关系是 qm = V_L，b = 1/p_L。这里的p_L是压力常数。比如一块煤样测得CH4的V_L = 1.2 mol/kg，p_L = 2 MPa，那么b = 0.5 1/MPa。CO2的典型V_L大一些，p_L小一些，这意味着更高的吸附容量和更强的亲和力——所以能置换CH4。

扩展到二元体系时需要同一套qm和b参数。不同文献的数值会互相矛盾，最稳的做法是：优先采用同一篇文献里同时给出CH4和CO2参数的实验数据，不要拼接两篇文章的数值。拼接很容易造成选择性颠倒——比如把两个不同温压条件下的数据进行混搭，计算结果和实验趋势对不上，这时候你会以为是建模问题，其实是数据本身有问题。

6.2 出口组分浓度、累计采气量怎么定义和输出

COMSOL里有一堆内置算子。出口端组分浓度的平均值用aveop或者intop算子。我习惯在出口端定义一个边界探针（Probe），在线监视c_CH4随时间的演化。突破曲线就是从c_CH4 = 1（初始纯甲烷）下降到0的过程，读出c_CH4降到0.1的时刻，就是通常意义上的“突破时间”。

累计采出量定义成出口边界上CH4通量对时间的积分。在COMSOL里可以用intop算子配合时间积分表达式，或者直接用“全局ODE”开辟一个变量跟踪累积量。要注意的是，如果你用定流量入口，岩芯内的气体总有压缩和吸附存储量，千万别直接用“注入体积乘以结束浓度”这种粗暴计算，那样的误差可达百分之十几。老老实实做边界积分，最稳妥。

模拟结束后，把实验测的出口浓度曲线和模拟曲线叠在一张图里，看两个特征：突破时间对不对，突破后的拖尾形态对不对。如果只修参数就能让两条线基本重合，说明物理场框架是对的；如果无论如何都套不上，回头检查吸附动力学常数和有效扩散系数。

6.3 哪些参数最敏感？我的判断供你参考

做了三轮参数扫描之后，我对这个模型的敏感性排序如下，从高到低分别是：CO2/CH4吸附选择性、吸附动力学常数、注入流量、有效扩散系数、渗透率绝对值。第一梯队是两个吸附参数，它们直接决定置换能力和拖尾形状；第二梯队是注入流量，它决定时间尺度；第三梯队是扩散系数，在低流量对流传质不占绝对优势的时候才明显。

这套排序意味着：如果你只有有限的实验数据做校准，优先校准吸附选择性，其次调吸附动力学常数。渗透率数值就算错了一个数量级，只要反馈机制还在，突破曲线形态只是整体平移，不会产生本质偏差。

话题回到COMSOL操作层面。这个模型用LiveLink for MATLAB或LiveLink for Python跑参数扫描会非常顺手。尤其是在COMSOL 6.4里用Python批量提交参数扫描，可以自动输出每个样本的突破时间和累计产量表格，配合matplotlib直接画敏感性图。我自己已经很少手动在GUI里一个个点扫描了，因为那太耗时，而且容易在等待结果时打乱思路。

写在最后的实操体会

COMSOL模拟CO2驱替甲烷这个课题，难点从来不在软件操作，而在“你清不清楚自己在模拟什么”。把达西流动、多组分传质、竞争吸附三条线先拆开，再逐条闭合，模型就不会是黑箱。我自己从初版模型到和实验数据对齐，大概花了三周，其中一半时间都耗在参数标定和收敛排查上，但思路一旦清晰，后面的各种参数扫描和模型扩展都是水到渠成的事。

如果这个模型下一步要往实际工程方向推，我建议加多孔介质传热场，把吸附放热和温度变化纳入进来。实验室恒温假设在矿场尺度和高压差条件下并不成立，温度会影响Langmuir参数和粘度，那是另一个维度的课题了。