多孔介质流固热三场全耦合数学模型及数值模拟.docx

1、多孔介质流固热三场全耦合数学模型及数值模拟第25卷增1岩石力学与工程学报V ol.25 Supp.1 2006年2月Chinese Journal of Rock Mechanics and Engineering Feb . ，2006多孔介质流固热三场全耦合数学模型及数值模拟盛金昌(河海大学水利水电工程学院，江苏南京 210098摘要：给出多孔岩体介质的流固热三场全耦合数学模型。假定流体为单相流，固体介质为非沸腾的饱和、热弹性多孔介质，该模型由流体物质守恒方程、力学平衡方程和能量守恒方程这3个相互耦合的方程组成，其中包含了众多耦合项，并定义一系列的本构关系及耦合变量。以FEMLAB 工具为

2、基础，将该数学模型转化成为一个统一的偏微分方程组，在人机交互的环境下，实现流固热三场全耦合数值求解，一次解出渗流场、位移场和温度场，给出更接近真实物理过程的数值解答，避免松散耦合法求解多场耦合问题带来的误差。利用一个已知解析解和数值解的算例来证明了耦合模型及求解方法的正确性。最后模拟通过井孔向岩体中注入冷水时流固热全耦合过程，详细地分析全耦合作用对井壁围岩应力的影响，计算结果表明：流固热三场耦合作用对井壁的稳定分析有非常重要的影响。关键词：岩石力学；多孔介质；流固热三场全耦合模型；数值模拟中图分类号：TU 457 文献标识码：A 文章编号：10006915(2006增1302806FULLY

3、COUPLED THERMO-HYDRO-MECHANICAL MODEL OF SATURATED POROUS MEDIA AND NUMERICAL MODELLINGSHENG Jinchang(College of Water Conservancy and Hydropower Engineering，Hohai University，Nanjing ，Jiangsu 210098，China Abstract ：A fully coupled Thermo-Hydro-Mechanical model is presented in a single-phase，non-boil

4、ing linear thermoelastic medium，which incorporates cross-coupled fluid flow equation，energy conservation equation and mechanical equilibrium equation with many cross-coupling terms. A series of constitutive relations and cross-coupling relations between material properties and independent variables

5、are defined in the model. The coupled multiphysics model is simultaneously simulated by using FEMLAB，the first engineering tool that performs partial differential equation-based multiphysics modeling in an interactive environment，which the mathematical model is translated into a set of partial diffe

6、rential equations. The pore pressure，displacements and temperature ，which should theoretically approach the most realistic results，can be solved simultaneously by using FEMLAB ，in which the errors in other coupling algorithms can be avoided. An example with known analytical and numerical results is

7、used to validate the multiphysics model. In particular，cold water injection into wellbore is modeled with realistic time steps indicating that the coupled processes have significant effects on the stresses of borehole wall and wellbore stability.Key words：rock mechanics；porous medium；fully coupled T

8、hermo-Hydro-Mechanical model；numerical simulation收稿日期：20050613；修回日期：20051115基金项目：国家自然科学基金重点项目(50539030第25卷增1 盛金昌. 多孔介质流固热三场全耦合数学模型及数值模拟 3029 1 引言在众多的工程领域中，多孔岩体介质中的流、固、热多场耦合问题具有非常重要的意义和作用。例如：在油气田开采、CO 2的地下储存、地热资源的开采、核废料的地下储存等领域中，都存在着多场耦合问题，多物理场之间的复杂耦合作用对工程具有非常重要的影响。在这些耦合作用中，既包括线性物理作用，也包括材料的非线性耦合作用。目前

9、，国外在多场耦合分析方面已进行了大量的工作，如进行了两场(变形和流动 或三场的耦合分析13，其理论一般都假设岩体为多孔弹性介质。这个耦合系统包括3个相互耦合的方程所组成的方程组以及一组物理参数之间的耦合关系，如：饱和度与渗透系数之间的关系、应力与渗透系数之间的关系等。地质体的耦合作用过程吸引了广泛的研究兴趣，目前在流、固、热等多场耦合分析方面已有不少研究成果48。多孔岩体介质中的流固热多场耦合问题非常复杂，数值求解具有很高的难度。求解多场耦合问题时多采用交叉迭代求解。耦合问题求解有3种基本算法3：单向耦合算法、松散耦合算法和全耦合算法。对于单向耦合算法，即2组独立的方程在同一时间步内分开求解，

10、求解时只是将其中的一个物理过程的计算结果作为另一个物理过程的输入，这种传递只是单向的。例如：由流动方程解出的孔隙压力作为荷载传给力学计算来求解应力和位移。J. Fredrich等4的工作就是采用这一算法。对于松散耦合算法，两组方程独立求解(和单向耦合算法一样 ，但是有关信息在指定的时间步内在两个求解器之间双向传递。这一算法的优点是相对容易实现，而且还能反映较复杂的非线性物理过程，它接近全耦合算法1，3，58。对于全耦合算法，需要推导出统一的一组全耦合方程组(通常是一个大型的非线性全耦合的偏微分方程组 ，这里面融合了所有的相关物理过程。求解多物理耦合问题应该首选全耦合算法，因为在理论上它能给出最

11、真实的结果。然而由于同时求解多物理耦合过程异常复杂，到目前为止还没有这方面令人信服的研究成果。本文的主要工作如下：在J. Noorishad和C. F. Tsang 1研究的基础上，首先建立一组多孔岩体介质瞬态的流、固、热三场全耦合方程组(其中包括固体 变形、能量传输、流体流动3个相互耦合的过程，多物理场之间的交叉耦合还包括材料性质与独立变量之间的耦合关系 ，然后利用FEMLAB 软件(基于偏微分方程组的多物理过程模拟工具 作为平台，成功地求解了这一全耦合偏微分方程组，避免了松散耦合法求解多场耦合问题带来的误差，实现了同时求解多物理场耦合过程。本文的求解方法是一个全耦合算法，在理论上它能给出最

12、真实的结果。通过对一个具有已知解析解和数值解的算例的计算分析来证明本文耦合模型及求解方法的正确性：一维砂柱的等温固结和非等温固结问题。最后模拟分析了通过井孔向岩体中注入冷水时流、固、热全耦合过程，详细地分析了全耦合作用对井壁围岩应力的影响。2 流固热三场全耦合数学模型在J. Noorishad和C. F. Tsang1工作的基础上，下面给出流固热三场全耦合数学模型。假设岩体介质为饱和的多孔弹性介质，并假设REV 存在。2.1流动方程流体流动时质量平衡方程的欧拉形式如下：QVt=+lll(1式中：为岩体孔隙率，l为流体的密度，t 为时间，lV 为流体速度矢量，Q 为流体的源汇项。根据流体流动的动

13、量方程可得Darcy 定律：(lllgPkV =(2式中：l为流体的动力黏滞系数，k 为孔隙介质的渗透率，P 为孔隙压力，g 为重力加速度矢量。将式(2代入式(1，并加上固体骨架的变形项，经过一系列推导可得+tTtPt tpvl QgPk=(lll (3式中：为流体的参考密度，v为岩体的体积应变，T 为温度，t为流体的热体积膨胀系数，p为流体的压缩系数。式(3即为多孔热弹性介质中流动流体的控制方程的最后形式。式(3等号左边的前3项分别表示3030 岩石力学与工程学报 2006年由于应变、流体压力和温度所引起的流体体积的变化，等号左边的最后1项代表由压力梯度和重力作用而引起的流体流量。 2.2

14、能量守恒方程固体骨架和流体共同存在于同一个体积空间，但它们具有不同的热动力学特性：如比热容和热传导系数等。因此固体骨架和流体的能量守恒方程需要分别定义。固体骨架的能量守恒方程定义如下：s s s p 1( ( 1( (1(q T K t Tc += (4式中：s p (c 为岩体骨架的热容，s K 为岩体骨架的热传导张量，s q 为岩体的热源强度。对于流体，相应的能量守恒方程可定义如下：l l l l p l p ( ( ( ( q T K T V c t Tc +=+(5式中：l p (c ，l K 和l q 分别为流体的热容、热传导张量和热源强度。对于单相流，假设固体和流体之间总是处于热平

15、衡状态，这样将式(4，(5迭加，并考虑到变形能，即可得到以下统一的能量守恒方程9：=+T V c tT t T c p p ( ( 1( (l l v 0t t t (q T K + (6 式中： 32(+=，和为拉梅常数，为各向同性固体的线性热膨胀系数；0T 为无应力状态下的绝对温度；t q 为充满了流体的多孔介质的热源汇项，定义为s l t 1( q q q +=。t p (c 和t K 分别为充满了流体的多孔介质的比热容和热传导系数，且s p l p t p (1( ( (c c c +=s l t 1( K K K += 2.3 力学平衡方程假设岩体为理想热弹性体，考虑流体的孔隙压力和

16、热应力的本构关系如下：P T ij ij kl kl ij ij ij +=2 (7式中：ij 为应力分量，ij 为应变分量，为Biot 系数，ij 为Kronecker 数。将式(7和应变位移关系代入静力平衡方程式，可以得到用位移表示的包含耦合项的修正的Navier 平衡方程：i i i ji j jj i F T P u u =+， 2/( (8式中：i F 为体积力分量。式(3，(6和(8一起完整地定义了岩体三场全耦合作用的数学模型。施加一定的边界条件和初始条件，即可求解上述耦合控制方程组。另外，渗透系数与应力之间的关系有很多表达式，本文采用以下关系：+= (21exp(21v v 0p

17、 k k (93 基于FEMLAB 全耦合分析的实施FEMLAB 是一个在交互式环境下求解基于偏微分方程组的多物理耦合过程的有限元分析工 具10，它的功能非常强大，可以应用于广泛的科学和工程领域。在FEMLAB 中，不再需要编制复杂的偏微分方程组的求解器，它含有一些内嵌的物理模型，如：化学工程模型、电磁学模型、结构力学模型等。但是其中功能最强大、而且也是最灵活的还是它的偏微分方程组模式。描述偏微分方程组有3个数学应用模式：系数形式(coefficient form、通式(general form和弱形式(weak form。其中弱形式最为灵活，但对大多数问题，通式模式已足够强大。本文选择结构力

18、学模型加通式模式微分方程组(描述耦合热传导方程和流体流动方程 来求解岩体三场全耦合问题，当然在结构力学模型和通式微分方程组模式之间还需定义很多交叉耦合项。在这些问题正确定义出来之后，FEMLAB 在求解时，将先把结构力学模型和通式微分方程组模式结合在一起转换成一个统一的通式形式的微分方程组，然后统一求解这个总的通式形式的微分方程组，同时解出位移场、渗流场和温度场，从而实现了三场的全耦合求解，避免了松散耦合法求解多场耦合问题带来的误差，给出了更接近真实物理过程的数值解答。求解结果可以用多种方式来表达，如等势线、曲线、图像及动画等。4 流固热全耦合模型的验证下面通过算例证明本文的三场全耦合数学模型

19、以及基于FEMLAB 的耦合数值模拟器的正确性和有效性。M. A. Biot2利用一个有限元耦合分析程序第25卷 增1 盛金昌. 多孔介质流固热三场全耦合数学模型及数值模拟 3031 (ROCMAS分析了一个砂柱的热弹性固结问题，这是一个典型的流固热耦合问题，有关研究1，2分别 给出了等温弹性固结解析解和非等温弹性固结的数值解。所用的物理参数见表1。用基于FEMLAB 的耦合模型求解了这两个固结过程。图1给出了本文解答与数值解(非等温弹性固结过程 和Biot 解析解2(等热固结过程 的对比情况，由图可知，本文的解答和相关研究1，2的结果非常接近，这说明本文的数值模型非常精确可靠。表1 验证算例

20、中所用的岩体材料参数Table 1 Rock mass parameters used in validation example弹性模量 E /Pa 泊松比 Biot 系数 孔隙率 渗透率k /m2 6.0103 0.41.0 0.24.0106热传导率 K t / (kJm 1 1比热容(c t /(kJm31线性热膨胀系数 / 1 8.36101167.03.0107 图中黑点代表本文模拟结果，实线代表已知解图1 本文模型的模拟结果与已知解析解和数值解的对比 Fig.1 Comparisons of simulation results with the knownsolutions5

21、应用实例分析在石油和天然气开采中，钻井过程中井壁稳定问题非常重要，实际上这一过程是一个流固热多场耦合问题。注入液体的压力和温度将改变井壁围岩的应力分布，从而最终影响井壁的稳定11。下面利用本文的全耦合模型求解该问题。 5.1 模型的描述图2给出了该算例的计算域和边界条件，这里计算域为10 m10 m，井孔半径为0.1 m。水平和竖图2 计算域及边界条件示意图Fig.2 Schematic simulation model and its boundary conditions直方向的边界压力分别为30和20 MPa。假设注入液体压力比岩体中的大10 MPa，注入液体是冷水，比岩体温度低50 ，

22、液体温度低有利于井壁的稳定。假设井壁孔隙压力和温度与注入液体的压力和温度完全相同。计算域中的初始压力和温度差为0。表2所示为应用实例中所用的岩体材料参数。 5.2 计算工况本文通过以下3种耦合和不耦合计算工况来说明耦合分析对井壁稳定分析的重要性，分别为：工况1：流固热全耦合计算分析； 工况2：把岩体看成不透水介质，即不考虑流动过程，但计算应力时把水压力加在井壁上。计入温度场引起的热应力；工况3：不考虑温度场和流场之间的耦合作用，即不考虑对流项对温度的影响和温度变化引起水压力的变化。但考虑水压力场和热应力对变形的影响。在FEMLAB 求解时，本文选用四阶Lagrange 型三角形单元。下面本文将

23、主要分析三场耦合作用对孔周围岩应力分布的影响。 5.3 切向应力的对比分析图3(a给出了在时刻为86 400 s时A -A 剖面(见图2 上不同工况下的切向应力对比图。由图3(a可知，在井壁上，用全耦合算法得出的切向压应力(压为负 值最大，工况2的切向压应力值最小。在点P ，工况1的切向应力值为25.14 MPa，工况2的为17.78 MPa，而工况3的为23.79 MPa，即不考虑耦合时点P 的切向应力值比考虑耦合情况下的值小很多。在岩体内部(离井壁距离大于0.2 m，3种工况计算出来的结果逐渐趋向一致，它们之间的差别大0.010.10.20.30.40.50.623456固结时间/s顶部沉

24、降量/m m3032 岩石力学与工程学报 2006年表2 应用实例中所用的岩体材料参数Table 2 Rock mass parameters used in application example密度s / (kgm 3弹性 模量 E /MPa热传导率K t /(kJm 1 s 1 1 热膨胀 系数/ (1比热容c ps/(kJkg 11孔隙率渗透率k /m2Biot 系数泊松比密度l /(kgm 3压缩系数p /(GPa1 动力黏滞系 数(20 /(Ns m 2 热膨胀系数t /(1比热容c pl /(kJkg 112.6103 5 0003.081031.01050.840.25.010

25、131.00.251 000.00.5131033.171044.2 A -A 剖面上自孔壁的距离/m (a 孔周切向应力分布 A -A 剖面上自孔壁的距离/m (b 孔周径向应力分布图3 各种计算工况下孔周应力沿A -A 剖面的分布 Fig.3 Variation of stresses along A -A section for differentcases致在2.0 MPa以内。由此可知，考虑流动过程与热传导过程之间的耦合和不考虑它们之间的耦合对切向应力的影响。图4(a给出了井壁上点A 切向应力随时间的变 化过程，图4中1，2，3的标注含义同图3。由图4(a可知：工况2和工况3条件下，

26、该点的切向应力值随时间变化很小，而工况1条件下，点A 切向应力则随着时间的增加而有一定幅度的减小。这表明，考虑耦合作用后，井壁切向应力随着流动和热传导过程有一定的调整。t / (104 s(a 点A 切向应力随时间的变化过程t / (104s(b 点B 切向应力随时间的变化过程图4 各种计算工况下孔周应力随时间的变化过程 Fig.4 Variation of stresses with time for different cases5.4 径向应力的对比分析图3(b给出了在时刻为86 400 s时A -A 剖面上不同工况下的径向应力对比图。由图3可知，3种工况得出的径向应力的大小和分布大致相

27、同，只是工况2(不耦合情况 值要小一些，大约小12 MPa。到岩体内部后(大于0.2 m，工况2比工况3的稍小一点。在点B (岩体内部 ，工况1的径向应力值为24.66 MPa，工况2的为23.38 MPa，而工况3的为24.98 MPa，可见3种工况下得出的岩体内部径向应力值差别不大。图4(b给出了岩体内部点B 径向应力随时间的变化过程，由图可知：该点的径向应力在3种工况之间的差别随着流动和热传导过程的发展，它们之间的差别越来越大。工况1条件下点B 的径向应力1流固热全耦合过程2只有固热耦合，岩体为不渗透介质，没有流动过程，水压力作为面力加在井壁上 3只有流热不耦合切向应力/M P a径向应

28、力/M P a切向应力/M P a径向应力/M P a第 25 卷 增1 盛金昌. 多孔介质流固热三场全耦合数学模型及数值模拟 3033 在前 5 000 s 随时间的增加而减小， 而此后则随时间 增加而增加，这主要是因为热传导过程比较快，在 前 5 000 s 时间内热应力对该点的径向应力的影响 起主要作用，而此后孔隙压力对该点的径向应力的 影响占主导作用。工况 2 条件下点 B 的径向应力单 调增大，而工况条件下则单调减小。不过它们之间 的差别不是太大，大约在 2 MPa 之内。 以上分析表明：采用耦合模型和采用非耦合模 型，得出的孔周应力的大小和随时间的变化有比较 大的差别，只有采用耦合

29、模型，得出的孔周应力才 是真实可靠的，这对井壁的稳定分析有重要价值。 明：只有采用耦合模型，得出的孔周应力才是真实 可靠的，这对井壁的稳定分析具有重要的参考价值。 参考文献(References： 1 Noorishad J，Tsang C F. Coupled thermohydroelasticity phenomena in variably saturated fractured porous rocks-formulation and numerical solutionA. In：Stephansson O，Jing L，Tsang CF ed. Coupled Thermo-Hy

30、dro-Mechanical Processes of Fractured MediaC. Amsterdam：Elsevier Science Publishers，1996. 93134. 2 Biot M A. General theory of three-dimensional consolidationJ. Journal of Applied Physics，1941，12：144164. 3 Minkoff S E，Stone C M，Bryant S，et al. Coupled fluid flow and geomechanical deformation modellingJ. Journal of Petroleum Science and Engineering，2003，38：375