CN109063383A - 基于微尺度重建模型的热-流-固多场耦合模拟方法 - Google Patents

Abstract

本申请公开了一种基于微尺度重建模型的热‑流‑固多场耦合模拟方法。该方法利用岩石微CT图像，基于重建的岩石孔隙尺度结构化网格模型，构建热‑流‑固三场耦合的数学模型，并通过数值模拟研究实现应力作用和温度变化对岩石孔隙结构演化和流动特性的预测。不同于传统宏观尺度的热‑流‑固耦合理论研究，本方法基于孔隙尺度的重建模型，从微观尺度深入揭示热‑流‑固多场耦合的相互作用机制。不仅有助于油气田开发过程中应力及温度参数对于开发效果影响的预测，同时，也能促进地下核废料处理、CO2封存以及地热开发等岩土工程领域内涉及多孔介质多场耦合问题的研究，更好的指导工程实践的开展。

Description

基于微尺度重建模型的热-流-固多场耦合模拟方法

技术领域

本发明涉及岩石物理技术领域，特别是涉及一种利用微CT图像，结合微尺度模型重建技术的岩石热-流-固多场耦合模拟方法。

背景技术

研究表明，自然界中的岩石作为一种强非均质性多孔介质，其中往往存在着气相、液相和固相相互作用的流-固耦合***。流-固耦合***对岩石的宏观物理力学特性及流体的输运特性起着重要的控制作用，同时这类流-固耦合***也往往受控于自然界中温度变化的影响。大量工程实践表明：诸如石油天然气开发、地热能利用、CO2地质封存和水电及隧道等与岩石相关的工程领域的建设中，均存在着应力场、渗流场和温度场等多场耦合的相互作用、相互影响的动态过程，从而使得工程建设过程的影响因素更加复杂化。因此，开展岩石热-流-固耦合机理的深入研究更有助于工程实践的实际需求。

总体来看，基于已有的研究成果，目前岩石热-流-固耦合理论主要是基于Terzaghi有效应力原理和Biot三维土体固结理论的基础，通过对原始基本假设的修正，使方程更接近于真实的岩石赋存环境。然而，现阶段对于岩石热-流-固耦合理论的研究往往基于宏观的连续介质思想，将岩石骨架和孔隙流体当做一个连续体元来表征多孔介质，局限于宏观尺度，从而忽略了岩石中大量存在的复杂无序的孔隙结构对于岩石力学性能、流体输运特性和热力学性能的影响。而这些恰恰是开展岩石热-流-固多场耦合机理研究的基础。由于真实环境的岩石中，固体骨架和孔隙流体各自占有一定空间，并且骨架在外部环境的作用下发生变形，从而导致孔隙形状的不断改变；同时在多相流体赋存的情况下，流体物性与分布形态也时时变化。因此，本申请公布了一种基于岩石CT扫描图像，利用孔隙尺度重建模型开展岩石热-流-固多场耦合理论研究的数值分析方法。该方法分别建立骨架和孔隙的孔隙尺度结构化网格模型，能更好的从微观的角度揭示多孔介质中热-流-固相互作用机制，为与岩石相关的工程实践提供科学的技术指导。

发明内容

本发明的目的在于提供一种基于微尺度重建模型的热-流-固多场耦合模拟方法。

为了达到上述目的，本申请采用下述技术方案来实现：

(一)以岩心微CT图像为基础数据，通过自主开发的模型重建算法构建孔隙尺度的结构化网格模型，分别构建固体骨架和孔隙流体的网格模型。

(二)在渗流力学和流-固耦合理论的基础上，通过引入与温度变化相关的参数对岩石变形与流体项进行修正，建立热-流-固耦合数学模型，作为后续数值模拟的控制方程。

(三)本方法采用双向耦合的计算方法，在重建的岩石骨架模型和孔隙结构模型上分别进行固体变形和流体场的求解，然后通过中间数据平台实现耦合变量的传递与求解控制。

(四)本方法中采用非耦合方程的形式进行求解，因此耦合效应的具体表现为：温度引起的热应变分别导致骨架模型和孔隙模型网格的变形与重新划分；骨架变形引起流体通道的改变，流体压力又导致骨架的变形；骨架的变形将影响温度场的分布情况，孔隙流体的流动以对流换热的形式改变骨架壁面处的温度，从而影响骨架的热传导。三者的相互影响通过在数值计算中不断的迭代收敛达到最终的平衡。

(五)由于本申请中采用的孔隙尺度模型，其物理尺寸较小，因此骨架模型上的温度梯度可以忽略不计，骨架与流体间的热传导也在极短的时间达到平衡，所以本方法主要研究不同温度对岩石变形与流体特性的影响。

(六)基于重建的孔隙尺度骨架模型和孔隙模型，结合构建的热-流-固耦合数学模型，分别采用ANASYS和CFX软件作为固体变形和流体流动的求解平台，MPCCI作为耦合数据的传递平台开展多场耦合的数值模拟计算。

与现有热-流-固多场耦合理论研究相比，本发明的有益效果在于：

本方法通过构建骨架和孔隙网格模型，分别形成固体域和流体域两套独立的网格***，通过定义的耦合边界实现耦合物理量的传递，通过数值模拟手段揭示热-流-固多场耦合的微观作用机制。

利用CFX集成于ANASYS中的优点，本方法中提出的双向耦合计算不需要第三方耦合软件的辅助，只需借助不具备求解功能的MPCCI软件作为耦合数据的传递平台即可实现耦合计算的求解，耦合求解实现过程简单，数值模拟软件易操作。

附图说明

为了更清楚地说明本发明方法的技术方案，以下结合附图和具体实施方式来进一步说明本发明方法。

图1为本申请方法的理论基础概念图。

图2为本申请方法实施流程示意图。

图3为本方法中采用的岩心原始微CT图像。

图4为基于微CT图像重建的孔隙尺度网格模型，(a)为固体骨架模型，(b)为流体孔隙模型。

图5为模型边界条件设置及耦合边界定义(a)骨架模型边界设置；(b)孔隙模型边界设置。

图6为骨架模型与孔隙模型的几何模型。

图7为模型孔隙度随围压的变化规律。

图8为模型渗透率随围压的变化规律。

图9为模型孔隙度随温度的变化规律。

图10为模型渗透率随温度的变化规律。

图11为模型渗透率随温度和有效压力的变化规律，(a)x方向；(b)y方向；(c)z方向。

图12为模型相对渗透率随围压的变化规律。

图13为模型相对渗透率随温度的变化规律。

具体实施方式

为了使本发明实现的基础理论、技术手段和达成目的易于说明，下面将结合本申请实施例中的附图，以及具体的计算案例对本申请实施例中的技术方案进行清楚、完整地描述。需要知悉的是，所描述的实施例仅是本申请一部分实施例，而不是全部实施例。基于本申请的实施例，本领域的其他技术人员在没有其他创新性劳动的前提下获得的所有其他实施例，都应当属于本申请保护范围。

本发明方法利用岩石微CT图像，基于微尺度重建模型开展岩石热-流-固多场耦合的数值分析方法实施例结合附图详述如下：

具体如图1所述，本申请所述的岩石热-流-固多场耦合理论中，温度场、应力场和变形场之间相互作用、相互影响的关系为：流体对温度的影响主要体现在热传导上，温度对流场的改变主要通过影响对流项；温度的改变会引起固体的热应变，固体的变形又会导致温度场的重分布；固体的变形会直接引起渗流通道的改变，流体压力又将影响有效应力从而改变固体变形量。因此，这三个物理场之间的相互作用和影响控制着岩石的宏观物理特性，从而决定着工程建设安全开发。

本方法的具体实施流程如图2所示：主要包含孔隙尺度模型重建、多场耦合数学模型构建、数值模拟计算三个主要的步骤。

S1：如图3所示为本实施例中采用的原始微CT图像数据，图中所示岩石样本为人造砂岩(MS1)，其余的岩石样本还有：碳酸盐岩(C1)、贝雷砂岩(B1)和疏松砂岩(S6)，本发明实施例的岩芯微观CT实验主要采用西南石油大学油气藏地质及开发工程国家重点实验室的蔡司Xradia MicroXCT-400微观成像***。

利用获取的岩心微CT图像数据，采用自主开发的模型重建技术实现孔隙尺度结构化网格模型的构建，为了减少重建模型的网格数量，提高模型重建的生成和后续数值计算的收敛速率，构建如图4所示的立方体模型(一般为400³～600³像素尺寸)。

S2：本申请中提出的方法有别于传统多场耦合理论研究，不是将岩石骨架和孔隙流体视作一个连续的体元来表征岩石多孔介质，而是分别构架骨架模型(图4a)和孔隙模型(图4b)实现变形场和流场的单独计算，再通过耦合参数的传递实现多场耦合的分析，计算更加精确，也更符合实际。

本文所述方法采用双向耦合的方式分别采用ANASYS和CFX软件求解固体变形场和渗流场，利用MPCCI软件作为数据交换平台实现中间耦合物理量的传递。整个***每一个步骤的计算称作一次大迭代，每一次大迭代中固体场和流体场分别进行迭代运算，单场运算收敛再进行耦合物理量传递，直到整个***的计算收敛为止。

S3：本方法中流体场的数值模拟计算采用CFX软件，因此采用N-S方程组作为基础控制方程来描述孔隙尺度的渗流。以油水两相流为例，不混相流体连续性方程如公式(1)所示：

其中α_i为计算单元中i_th项流体的体积分数，ρ_i和是i_th项流体的密度和速度，S_ai为源汇项，可为用户自定义的质量源项。

在整个流体计算域内采用单一动量方程控制，如公式(2)所示。域内速度场共用，表现形式为控制体内各相物性参数(密度和黏度)取体积平均。

其中ρ和μ分别为密度和粘度，p代表压力梯度，代表体积力，孔隙尺度模型由于尺寸太小，可以忽略不计，项描述Laplace压力效应，其值仅在流体相界面处非0。与域内各项体积分数相联系，密度和黏度的体积平均可以分别表示为：

其中水的黏度随温度的变化可以表示为：

油的黏度随温度的变化关系式与油的组分相关，但一般参照如下关系：

lgμ₀＝A+BT (6)

其中系数A和B可以通过油的黏温曲线拟合得到。

计算域内各项共用能量方程，如公式(7)所示：

其中为i项流体有效导热系数，S_h为热辐射及其他体积热源的贡献，E和T分别为各相能量和温度的加权平均，分别如下式所示：

S4：本方法中的岩石变形场采用ANASYS软件进行数值模拟，岩石采用等向理想塑性强化模型，弹性阶段的骨架变形场数学模型主要由平衡方程。几何方程和物理方程构成，分别如公式(10)，(11)和(12)所示：

其中σ_ij为应力张量，f_i为体积力。

其中ε_x，ε_y和ε_z为法向应变，γ_xy，γ_yz和γ_zx为剪应变，u和v和ω为位移分量

其中E和μ分别为弹性模量和泊松比。

变形场中，骨架的热传导遵循傅里叶热传导定律：

其中q_n″为热流密度，k为导热系数。

岩石骨架和孔隙流体之间的热对流现象则遵循牛顿准则：

q_n″＝h(T_s-T_f) (14)

其中h为对流换热系数，T_s为骨架表面温度，T_f为周围流体温度。

S4：进行数值模拟计算之前，需要进行模型边界条件的设定与耦合边界的定义。固体变形通过ANASYS进行求解，流体流动通过CFX进行求解，并且以Workbench平台为基础进行模型边界条件添加。如图5(a)所示，骨架模型四个侧面分别施加围压，图5(b)所示为流体场上下两面分别施加进出口压力，同时定义固体与流体的交界面为流-固耦合界面。图6所示为在Workbench中建立的骨架与孔隙的几何模型文件。

由于本实施例中模型尺寸很小，因此骨架模型中的温度梯度可以忽略不计；同时由于固体与流体间的热交换也在极短时间内达到热平衡，所以实施例中主要考虑了不同的温度下热应变导致的孔隙结构变化及对渗流的影响，因此对骨架模型设置恒定温度边界条件。耦合计算从渗流场开始求解，并将计算得到的孔隙压力传递给固体场求解器；固体场计算得到的边界节点变形量传递给流体场的壁面，如此直至***迭代收敛，完成计算。

S5：基于重建的模型及构建的热-流-固多场耦合数学模型，实施例1通过数值模拟的方法研究了应力和温度对岩石孔隙结构演化和绝对渗透率影响。

其中模拟采用的微观岩石力学参数由岩石微米压痕实验测得，岩石的导热系数和膨胀系数采用经验值，如表1所示为本实施例中采用的不同岩样的物性参数。

表1.岩石物性参数表

如图7所示，在孔隙压力不变的情况下，模型孔隙度随围压增大而减小，当围压增大到30MPa时下降趋势变缓，因为岩样内部的塑性区域增多。图8所示为渗透率随围压的变化规律，随着围压增大渗透率变小，主要是渗流通道压缩变形导致，此外随着塑性区的增大，渗透率下降幅度变缓，因为本文将岩样假设为各向同性强化弹塑性材料，因此塑性区的应变不会继续增大。

在假设围压和孔隙压力不变的情况下，同时不考虑热应变引起的裂纹萌生等因素(设定温度范围20℃～100℃)。图9和图10所示分别为孔隙度和渗透率随温度的变化关系。可以看出，温度升高，随着骨架膨胀导致孔隙度降低，但是此温度变化范围内孔隙度的下降幅度有限，仅为2％左右。通过不同岩样模型渗透率与温度的变化曲线可以看出，温度对膨胀系数较大和原始渗透率较低的模型影响更大。图11(a)、(b)和(c)可以看出岩心渗透率在空间上的非均质性，三个方向渗透率随应力和温度的变化情况说明：孔隙尺寸小、连通性差和基准渗透率越低的岩样对应力和温度的敏感性越高。

S6：实施例2研究了应力和温度对于水驱油效果的影响。选用了孔喉尺寸较大，连通性好，且空间渗透率非均质性较弱的人造砂岩样品MS1。

应力对相对渗透率影响的研究在恒定温度和孔隙压力的假设下进行。数值模拟中采用的油水物性参数如表2所示。

表2.油水物性参数表

如图12所示，在孔隙压力一定的情况下，相对渗透率曲线随着围压的增大逐渐向右下方移动，等渗点右移，两相共渗面积减小，且油水两相渗透率均呈下降趋势，残余油饱和度升高。这主要是由于随着模型有效应力增大，孔隙结构受压变形，部分狭窄喉道收缩甚至闭合，模型渗透率下降，油水流动能力下降，导致更多的残余油无法被驱替。

温度对相对渗透率影响的研究在恒定围压和孔隙压力的假设下进行。数值模拟中采用的油水物性参数如表3所示。

表3.油水物性参数表

模拟流体介质采用地层水和原油，地层水的黏温曲线采用水的黏温曲线；模型近壁面出的黏性系数采用微流边界层理论进行修正，即模型中水的黏性系数有下式给出：

其中μ_w为地层水在273K时的动力黏度，当忽略温度对为流边界层的影响时，可以采用不随温度变化的作用系数Φ。

本实施例中选用的岩样埋深较浅，地层温度低于100℃，因此以100℃以内的普通稠油样品的黏温曲线为基础开展模拟研究，忽略温度对微流边界层的影响。流体温度由骨架温度确定，并计算得到相应的油水黏性系数。依据原油黏温曲线结合宾汉流体基本公式，综合考虑微流边界层效应，可以得到微流通道中的流变关系时：

其中τ₀和μ_p分别代表极限切应力和结构黏度，可以通过原油流变曲线进行拟合。

如图13所示，随着温度的升高，曲线右端点及等渗点均发生右移，油相流动能力提高，相对渗透率增大。同时原油黏度随温度升高而降低，油水流动性能均增强，束缚水饱和度升高，残余油饱和度降低。油水流度比减小，弱化了水驱过程的指进效应，提高了最终采收率。虽然温度的升高引起孔喉收缩，但是黏度的降低带来的水驱效果的提升明显高于热应变带来的负面影响。

以上所述仅为本发明的较佳实施例，用于描述本发明的基本原理、特征和主要优点，并不用以限制本发明，凡是依据本发明的技术实质对以上实施方式所做的任何修改，等同变化与修饰，均属于本发明技术方案的范围内。

Claims (9)

1.基于微尺度重建模型的热-流-固多场耦合模拟方法，其特征在于：

(1)以岩心微CT图像为基础重建孔隙尺度结构化网格模型；

(2)基于两相流体力学和热-流-固耦合理论，建立了热-流-固两相渗流数学模型；

(3)以ANASYS和CFX软件作为耦合数值计算平台，分别实现变形域和流体域的数值求解；以MPCCI软件作为耦合数据交换平台实现数据交换和求解器控制。

2.按照权利要求1所述基于微尺度重建模型的热-流-固多场耦合模拟方法，其特征在于：所述步骤(1)中，通过数字图像处理技术对原始CT图像进行处理，分别提取岩石骨架和孔隙结构，并利用自主开发的模型重建算法构建岩石的孔隙尺度结构化网格模型。

3.按照权利要求1所述基于微尺度重建模型的热-流-固多场耦合模拟方法，其特征在于：所述步骤(2)中，基于传统的渗流力学和流-固耦合理论，通过引入温度项对流体场和变形场进行修正，构建热-流-固多场耦合数学模型。

4.按照权利要求1所述基于微尺度重建模型的热-流-固多场耦合模拟方法，其特征在于：所述步骤(3)中，以Workbench平台为基础，分别对岩石骨架网格模型和孔隙网格模型进行边界条件添加和耦合边界定义。

5.按照权利要求3所述基于微尺度重建模型的热-流-固多场耦合模拟方法，其特征在于：将构建的数学模型利用UDF(User Define Function)在数值模拟软件中实现编译。

6.按照权利要求1所述基于微尺度重建模型的热-流-固多场耦合模拟方法，其特征在于：所述步骤(3)中，将岩石骨架模型在ANASYS软件中进行数值模拟，同时引入热-固耦合模块的计算，实现固体变形域及温度影响的求解与分析。

7.按照权利要求1所述基于微尺度重建模型的热-流-固多场耦合模拟方法，其特征在于：所述步骤(3)中，将孔隙模型在CFX软件中进行数值模拟，同时引入热-流耦合模块的计算，实现流体域和温度影响的求解与分析。

8.按照权利要求1所述基于微尺度重建模型的热-流-固多场耦合模拟方法，其特征在于：所述步骤(3)中，利用MPCCI软件作为数据交换和求解器控制平台，实现固体域和流体域计算参数的传递。

9.按照权利要求8所述基于微尺度重建模型的热-流-固多场耦合模拟方法，其特征在于：通过对数值计算结果的后处理，分析应力条件及温度参数对岩石孔隙结构演化及输运特性的影响。