1. CO₂驱油模型概述与核心原理CO₂驱油技术作为提高原油采收率EOR的重要手段近年来在油气开采和碳封存领域获得广泛应用。这项技术的本质是通过注入超临界CO₂来降低原油粘度、改善流度比同时实现部分碳封存。在数值模拟层面我们需要同时考虑三个关键物理过程多孔介质流动遵循达西定律渗透率场呈现高度非均质性物质传输包含对流-扩散机制和相间质量交换相变行为超临界CO₂与原油的混相/非混相过程实际工程中这三个过程通过以下控制方程耦合∇·(ρu) 0 // 质量守恒 u -(k/μ)∇p // 达西定律 ∂(φc)/∂t ∇·(uc) ∇·(D∇c) q // 组分输运其中φ为孔隙度k为渗透率张量D为扩散系数张量。在COMSOL中这三个物理场通过共享变量如压力p、浓度c自动耦合。注意初学者常犯的错误是直接使用默认的物理场接口。正确的做法是在多物理场节点显式添加耦合项特别是相变导致的源项q需要手动定义。2. 几何建模与材料定义2.1 储层几何构建采用二维简化模型时建议使用20m×5m的矩形区域代表储层剖面。岩石颗粒的生成有几种策略随机分布法适合快速验证num_grains 50; // 颗粒数量 x_pos rand(num_grains,1)*20; // X坐标 y_pos rand(num_grains,1)*5; // Y坐标 radii 0.3 rand(num_grains,1)*0.7; // 半径0.3-1m规则排列法适合理论分析[x_pos,y_pos] meshgrid(1:2:19, 1:2:4); radii ones(size(x_pos))*0.8;实际岩心CT扫描数据导入需配合LiveLink for MATLAB实操技巧颗粒体积分数建议控制在30%-40%过高会导致孔隙连通性差收敛困难。可通过调整radii系数控制。2.2 材料参数设置渗透率场的定义需要区分岩石骨架和孔隙空间if (sqrt((x-x_pos)^2(y-y_pos)^2) radii) { k 1e-15; // 岩石骨架渗透率(m²) } else { k 2e-13; // 孔隙空间渗透率 }更专业的做法是采用平滑过渡函数k 1e-15 (2e-13-1e-15)/(1 exp(-10*(sqrt((x-x_pos)^2(y-y_pos)^2)-radii)))关键材料参数参考值参数岩石骨架孔隙空间渗透率(m²)1e-152e-13孔隙度0.020.25压缩系数(1/Pa)1e-104e-103. 多物理场耦合设置3.1 达西流与输运耦合在多物理场节点添加多孔介质流与浓物质传递耦合接口时需要特别注意密度耦合勾选质量守恒-密度耦合选项非等温效应如果考虑温度场需添加热力学属性相变源项通过弱贡献形式添加典型的注入井边界条件设置// 入口边界 P_in 10e6; // 10MPa注入压力 c_CO2 1.0; // 纯CO₂注入 // 出口边界 P_out 5e6; // 生产井压力3.2 状态方程处理CO₂物性参数计算推荐两种方法简化模型适合快速计算// 压缩因子近似 Z 1 0.3*(P/1e6) - 0.1*(T/50); rho_CO2 (P*44.01)/(Z*8.314*(T273.15));精确模型需自定义函数// 调用Span-Wagner方程 rho_CO2 CO2_density_SW(P,T);避坑指南当压力超过7.38MPa、温度超过31.1℃时必须使用超临界状态方程否则会出现物性突变导致求解失败。4. 求解器配置技巧4.1 稳态求解策略建议采用分步求解策略先求纯流动场冻结浓度场再求输运场固定流场最后全耦合求解求解器参数设置方法PARDISO直接求解器 非线性方法牛顿迭代 容差1e-4可放宽至1e-3初算4.2 瞬态求解优化对于时间相关的问题使用向后差分公式(BDF)初始时间步长设为总时间的1%启用自动时间步进典型设置时间范围0到1e6秒 步长控制严格 最大步长1e5秒5. 后处理与结果分析5.1 波及效率计算定义波及系数有两种方法网格遍历法适合教学演示double swept_area 0; for (int i0; imesh.numElements; i) { if (c_CO2[i] 0.1) { swept_area mesh.elementVolume(i); } } swept_eff swept_area/total_area;积分算子法工程推荐swept_eff intop1((c_CO20.1)*1)/intop1(1);5.2 典型结果可视化必看的后处理内容压力场云图观察驱动效果CO₂浓度场动画分析前缘推进生产井流量监测评估采收率高级分析技巧在流线图中添加粒子追踪计算局部流速矢量场导出特定位置的参数随时间变化6. 常见问题排查手册6.1 收敛问题处理现象可能原因解决方案初始步不收敛初始猜测不合理使用物理场初始化中间步发散物性突变减小时间步长周期性震荡网格质量差重构边界层网格6.2 物理场异常诊断当出现反常识结果时按以下顺序排查网格质量检查最大长宽比100最小单元质量0.3边界层足够密集边界条件验证注入/生产量守恒压力梯度方向合理无冲突边界设置材料属性复核单位制一致性参数数量级正确温度/压力适用范围7. 模型扩展方向热力耦合添加能量方程考虑温度场影响组分模型定义油相多组分相平衡地质力学耦合孔隙变形效应井筒模型使用Pipe Flow模块耦合对于想深入研究的同学建议尝试使用参数化扫描研究注入压力影响添加随机裂隙网络增强真实性调用外部MATLAB函数实现复杂本构关系这个模型虽然简化但已经包含了CO₂驱油的核心物理过程。我在实际项目中发现把注入压力控制在最小混相压力(MMP)附近时采收率会出现跃升这个现象值得用参数化扫描深入研究。另外建议输出驱替前缘的实时位置数据这对优化井网布置很有帮助。 CO₂驱油技术：原理、建模与数值模拟实践
1. CO₂驱油模型概述与核心原理CO₂驱油技术作为提高原油采收率EOR的重要手段近年来在油气开采和碳封存领域获得广泛应用。这项技术的本质是通过注入超临界CO₂来降低原油粘度、改善流度比同时实现部分碳封存。在数值模拟层面我们需要同时考虑三个关键物理过程多孔介质流动遵循达西定律渗透率场呈现高度非均质性物质传输包含对流-扩散机制和相间质量交换相变行为超临界CO₂与原油的混相/非混相过程实际工程中这三个过程通过以下控制方程耦合∇·(ρu) 0 // 质量守恒 u -(k/μ)∇p // 达西定律 ∂(φc)/∂t ∇·(uc) ∇·(D∇c) q // 组分输运其中φ为孔隙度k为渗透率张量D为扩散系数张量。在COMSOL中这三个物理场通过共享变量如压力p、浓度c自动耦合。注意初学者常犯的错误是直接使用默认的物理场接口。正确的做法是在多物理场节点显式添加耦合项特别是相变导致的源项q需要手动定义。2. 几何建模与材料定义2.1 储层几何构建采用二维简化模型时建议使用20m×5m的矩形区域代表储层剖面。岩石颗粒的生成有几种策略随机分布法适合快速验证num_grains 50; // 颗粒数量 x_pos rand(num_grains,1)*20; // X坐标 y_pos rand(num_grains,1)*5; // Y坐标 radii 0.3 rand(num_grains,1)*0.7; // 半径0.3-1m规则排列法适合理论分析[x_pos,y_pos] meshgrid(1:2:19, 1:2:4); radii ones(size(x_pos))*0.8;实际岩心CT扫描数据导入需配合LiveLink for MATLAB实操技巧颗粒体积分数建议控制在30%-40%过高会导致孔隙连通性差收敛困难。可通过调整radii系数控制。2.2 材料参数设置渗透率场的定义需要区分岩石骨架和孔隙空间if (sqrt((x-x_pos)^2(y-y_pos)^2) radii) { k 1e-15; // 岩石骨架渗透率(m²) } else { k 2e-13; // 孔隙空间渗透率 }更专业的做法是采用平滑过渡函数k 1e-15 (2e-13-1e-15)/(1 exp(-10*(sqrt((x-x_pos)^2(y-y_pos)^2)-radii)))关键材料参数参考值参数岩石骨架孔隙空间渗透率(m²)1e-152e-13孔隙度0.020.25压缩系数(1/Pa)1e-104e-103. 多物理场耦合设置3.1 达西流与输运耦合在多物理场节点添加多孔介质流与浓物质传递耦合接口时需要特别注意密度耦合勾选质量守恒-密度耦合选项非等温效应如果考虑温度场需添加热力学属性相变源项通过弱贡献形式添加典型的注入井边界条件设置// 入口边界 P_in 10e6; // 10MPa注入压力 c_CO2 1.0; // 纯CO₂注入 // 出口边界 P_out 5e6; // 生产井压力3.2 状态方程处理CO₂物性参数计算推荐两种方法简化模型适合快速计算// 压缩因子近似 Z 1 0.3*(P/1e6) - 0.1*(T/50); rho_CO2 (P*44.01)/(Z*8.314*(T273.15));精确模型需自定义函数// 调用Span-Wagner方程 rho_CO2 CO2_density_SW(P,T);避坑指南当压力超过7.38MPa、温度超过31.1℃时必须使用超临界状态方程否则会出现物性突变导致求解失败。4. 求解器配置技巧4.1 稳态求解策略建议采用分步求解策略先求纯流动场冻结浓度场再求输运场固定流场最后全耦合求解求解器参数设置方法PARDISO直接求解器 非线性方法牛顿迭代 容差1e-4可放宽至1e-3初算4.2 瞬态求解优化对于时间相关的问题使用向后差分公式(BDF)初始时间步长设为总时间的1%启用自动时间步进典型设置时间范围0到1e6秒 步长控制严格 最大步长1e5秒5. 后处理与结果分析5.1 波及效率计算定义波及系数有两种方法网格遍历法适合教学演示double swept_area 0; for (int i0; imesh.numElements; i) { if (c_CO2[i] 0.1) { swept_area mesh.elementVolume(i); } } swept_eff swept_area/total_area;积分算子法工程推荐swept_eff intop1((c_CO20.1)*1)/intop1(1);5.2 典型结果可视化必看的后处理内容压力场云图观察驱动效果CO₂浓度场动画分析前缘推进生产井流量监测评估采收率高级分析技巧在流线图中添加粒子追踪计算局部流速矢量场导出特定位置的参数随时间变化6. 常见问题排查手册6.1 收敛问题处理现象可能原因解决方案初始步不收敛初始猜测不合理使用物理场初始化中间步发散物性突变减小时间步长周期性震荡网格质量差重构边界层网格6.2 物理场异常诊断当出现反常识结果时按以下顺序排查网格质量检查最大长宽比100最小单元质量0.3边界层足够密集边界条件验证注入/生产量守恒压力梯度方向合理无冲突边界设置材料属性复核单位制一致性参数数量级正确温度/压力适用范围7. 模型扩展方向热力耦合添加能量方程考虑温度场影响组分模型定义油相多组分相平衡地质力学耦合孔隙变形效应井筒模型使用Pipe Flow模块耦合对于想深入研究的同学建议尝试使用参数化扫描研究注入压力影响添加随机裂隙网络增强真实性调用外部MATLAB函数实现复杂本构关系这个模型虽然简化但已经包含了CO₂驱油的核心物理过程。我在实际项目中发现把注入压力控制在最小混相压力(MMP)附近时采收率会出现跃升这个现象值得用参数化扫描深入研究。另外建议输出驱替前缘的实时位置数据这对优化井网布置很有帮助。