地质力学学报  2008, Vol. 14 Issue (2): 141-148
引用本文
王连捷, 孙东生, 周春景, 王薇, 赵卫华, 王红才. ANSYS软件在求解地应力与流体耦合作用中的应用[J]. 地质力学学报, 2008, 14(2): 141-148.
WANG Lian-jie, SUN Dong-sheng, ZHOU Chun-jing, Wang Wei, ZHAO Wei-hua, WANG Hong-cai. APPLICATION OF THE ANSYS SOFTWARE TO SOLVING COUPLING OF IN-SITU STRESS AND FLUID[J]. Journal of Geomechanics, 2008, 14(2): 141-148.
ANSYS软件在求解地应力与流体耦合作用中的应用
王连捷 , 孙东生 , 周春景 , 王薇 , 赵卫华 , 王红才     
中国地质科学院地质力学研究所, 北京 100081
摘要:地应力与流体的耦合作用是近年来研究的热点, 是工程和地质研究中经常遇到的重要课题.因此, 寻找简便的求解流固耦合的方法具有重要的意义.本文根据流固耦合方程和热传导方程的对应关系, 找出了在一定边界条件下, 可用ANsYS中的结构力学模块和热传导模块求解流固耦合方程的简便方法.并将该方法所得结果与理论解和实际结果进行了对比, 结果显示他们具有很好的一致性, 说明简便方法是可行的.
关键词地应力    ANSYS软件    流固耦合    Biot方程    
文章编号:1006-6616(2008)02-0141-08
APPLICATION OF THE ANSYS SOFTWARE TO SOLVING COUPLING OF IN-SITU STRESS AND FLUID
WANG Lian-jie , SUN Dong-sheng , ZHOU Chun-jing , Wang Wei , ZHAO Wei-hua , WANG Hong-cai     
Institute of Geomechanics, Chinese Academy of Geological Sciences, Beijing 100081, China
Abstract: The coupling of the in-situ stress and fluid is the research focus in recent years and also an important subject in research on rock and soil engineering and geology.So it has prominent significance to find out a simple method to solve the stress-fluid coupling equation.In this paper, a simple method of solving the stress-fluid coupling equation at certain boundary conditions is proposed based on the corresponding relation of the seepage equation of stress-fluid coupling and heat conduction equation and using the structural mechanical module and heat conduction module in the ANSYS software.The results obtained from the method are in good agreement with the theoretical solution and actual outcomes.
Key words: in-situ stress    ANSYS software    solid-fluid coupling    Biot equation    
0 前言

地应力与流体的耦合作用(简称流固耦合)是近年来研究的热点, 是工程和地质研究中经常遇到的重要课题, 如岩土的固结、地面沉降、油气运移、油井出砂、油井套损、采油引起的井压变化、水力压裂岩石的破坏、煤矿瓦斯突出、煤层底板突水、边坡灾害、地震活动等现象都与流固耦合有关。因此, 寻找简便的求解流固耦合的方法具有重要的意义。

ANSYS软件是一种大型通用有限元分析软件, 有强大的前后处理功能, 具有广泛的应用性和普及性, 因而基于ANSYS软件基础上的求解流固耦合的简便方法具有更广泛的应用性。因为渗流场与温度场具有相似的微分方程, 所以可以用ANSYS软件中的温度模块求解渗流问题, 但由于流固耦合问题, 既要考虑渗流场, 还要考虑应力场, 所以流固耦合问题更加复杂。应用ANSYS软件简化方法只能在一定条件下求解流固耦合问题。本文将讨论在何条件下可用ANSYS软件简化方法对流固耦合问题求解, 以及求解方法。

这里所说的利用ANSYS软件简化方法求解是指:①利用ANSYS软件的结构模块解出应力, 然后求该应力作用下产生的加载时刻孔隙压力(初始孔隙压力); ②利用ANSYS软件中的热传导方程求解渗流方程, 解出孔隙压力和流速随时间的变化。

1 地应力与流体耦合方程[1~2]
(1)
(2)
(3)
(4)

θ=0时(在加载的瞬时, 流体尚未扩散, 此时θ=0), 由(4)式可得:

(5)

R =H时,

(6)
(7)

σ0称为体应力, 对有效应力来讲, 它称为总应力。

用位移表示的平衡方程:

(8)

耦合渗流方程[3~4] :

(9)

以上方程也称为Biot方程或流固耦合方程。

上述各式中p为孔隙压力; θ为单位体积上的孔隙介质中流体体积变化; 对三维问题, σii = (σx +σy +σz) = (σ1 +σ2 +σ3); 对二维问题σii = (σx +σz) = (σ1 +σ3), 对一维问题σii =σz =σ3; H, R, Q, α为Biot系数[2], δij为单位矩阵, ; kij为渗透系数。σij, εij为应力和应变, εv为体应变, εv112233。ε11, ε22, ε33为正应变; G为剪切模量; λ为拉梅常数, , ν为泊松比; k体积模量k=E为弹性模量; ρ为体积密度。

求解流固耦合问题需要联立求解方程(8)和方程(9)并进行迭代, 以便得出应力、孔隙压力、流速以及他们随时间的变化[5]。迭代求解是比较麻烦的, 因而我们寻求简便方法。

上述方程可退化为一维和二维问题的方程。

2 温度场方程与耦合渗流方程的对应关系

ANSYS中温度场方程如下[6] :

(10)

ρ为密度, C为比热, T为温度, K为热传导系数, Qt为热源密度。

比较式(9)和(10)可知, ρc对应, Qt对应, 孔隙压力ρ与温度T对应。用ANSYS求解流固耦合方程时, 必须满足:

(11)
(12)

式(11)容易满足, 但式(12)得到满足比较困难, 只有当或为常数时, 式(12)才满足, 下面我们观察, 在何种情况下, 式(12)得到满足。

3 一维问题

一维情况下流体的流动是简单的, 假定一维问题是z方向, 在此方向作用有固定荷载F, 则总应力为σ0 =σz =F/s, s为荷载作用的面积。由式(7)可知, 载荷F作用瞬时的孔隙压力p为:p =σ0 =σz =F/s。当孔隙压力耗散终了时, 得到的应力仍为σz =σ0 =F/s, 这与弹性解是一样的。可见, 总应力是不随孔隙压力变化而变化的, 且可由弹性解求出。因为总应力σ0是不变的, 即, 故可用ANSYS简便方法求解, 此时Biot方程为[1] :

(13)

式中c为固结系数。此方程也称为太沙基方程, 上述方程在给定初始条件和边界条件后, 可用ANSYS中的热传导模块求解。而初始孔隙压力p可由ANSYS软件的线弹性求解。由ANSYS软件的结构模块求出应力σz, 由下式求出初始孔隙压力p, 由式(13)解出孔隙压力p随时间的变化。

(14)
4 三维问题

我们通过一个模型来研究不同边界条件时, 三维情况下总应力及其随流体耗散的变化, 观察在何条件下为零或为常数, 即适合于用ANSYS简便方法求解的条件。模型如图 1所示, 模型各表面的名称如下:面1485称为x1面; 面2376称为x2面; 面1265称为y1面; 面5678称为y2面; 面1234称为z1面; 面5678面称为z2面; 下面给出了几种边界条件, 如表 1

图 1 三维计算模型 Figure 1 3D model for calculation

表 1 三维模型边界条件 Table 1 Boundary conditions of the 3D model

表 1中边界条件的说明如下:①xy方向无变形, z方向加力, 即εxy =0, σz ≠0;②xy方向自由变形, z方向加力, 即σx =σy =0, σz ≠0;③x方向无变形, y方向自由, z方向加力, 即εx =0, σy =0, σz ≠0;④x方向自由, y方向无变形, z方向加力, 即σx =0, εy =0, σz ≠0; ⑤xyz方向都加力, 即σx ≠0, σy ≠0, σz ≠0;⑥x方向自由变形, yz方向加力, 即σx =0, σy ≠0, σz ≠0;⑦ xz方向加力, y方向自由变形, 即σx ≠0, σy =0, σz ≠ 0;⑧x方向无变形, y方向自由变形, z方向加力, 即εx =0, σy =0, σz ≠0;⑨x方向自由变形, y方向无变形, z方向加力, 即σx =0, εy =0, σz ≠0;

让我们求出加载瞬时的总应力、孔隙压力和孔隙压力耗散终了时的应力, 对两种情况进行比较, 观察总应力是否随时间变化。

加载瞬时, 由于流体尚未耗散, 因而流体的体积变化为零, 即θ=0, 故由式(4)可得:

(15)

根据边界条件, 由式(15)及(3), 可求解加载瞬时的总应力和孔隙压力。

当孔隙压力耗散终了时, 即孔隙压力为0时, 式(3)变为:

(16)

可以看出, 上式恰好为弹性本构方程。利用式(16)可计算孔隙压力耗散终了时的应力, 而这时的应力恰好就是线弹性应力。表 2为不同边界条件的计算结果。

表 2 三维模型不同边界条件的计算结果 Table 2 Results of calculation of the 3D model for different boundary conditions

由上表可以看出:

1.在三维情况下, 只要在任一方向上不允许有变形, 即εij =0 (i =j), 如边界条件1、3、4、8、9, 则加载时刻的总应力与孔隙压力耗散终了时刻的总应力是不同的, 即总应力是变化的, , 所以不能用ANSYS软件的简便方法求解;

2.在边界条件2、5、6、7情况下, 即在(x, y, z)任一方向上允许有变形, 或施加了载荷, 此时的加载瞬时和流体耗散终了时的总应力是相同的, 不随耗散时间变化, 即 =0, 此时可用ANSYS软件求解。也就是用线弹性方程求解总应力, 按式(7)求出孔隙压力, 也是初始孔隙压力。

3.得到初始孔隙压力和边界条件后, 用ANSYS软件的热传导模块求解渗流方程, 得出流体压力、压力梯度、流速、流量等随时间的变化及空间分布。

5 二维问题

对二维问题我们同样通过不同边界条件的模型来观察。图 2为不同边界条件的二维模型。

a.左、右边界x方向约束, 底边界z方向约束; b.左边界x方向约束; 底边界z方向约束 图 2 二维模型边界条件 Figure 2 Boundary conditions of a 2D model

图 2a为边界条件1, 图 2b为边界条件2。假定边界力为σz =1, σx =2, 计算结果如表 3

表 3 二维模型不同边界条件的计算结果 Table 3 Results of the 2D model for different boundary conditions

表 3可以看出, 对二维问题, 边界条件1时, 总应力随耗散时间而变化, 边界条件边2时, 总应力是不变化的。因此, 对边界条件2的情况, 可用ANSYS软件简化方法求解。

6 实例

通过实例检验ANSYS软件简便方法求解流固耦合的合理性。

(1) 一维空间中的固结问题求解

假定模型厚度z =10米(图 3), 其初始条件和边界条件为t =0时刻, z =0处作用荷载σz=1MPa, 在模型底部z =10处, 位移u =0, 且流体是不可渗透的, 即。在顶部z = 0处为透水边界, 即p =0。

图 3 一维固结模型 Figure 3 1D consolidation model

用ANSYS软件进行了求解, 并与解析解的结果进行了对比如图 4图 5。解析解的表达式如下:

图 4 底部z =10米处孔隙压力随时间的变化 Figure 4 Variation of pore pressures at bottom (z =10m) with time

图 5 20天时孔隙压力和有效应力的分布 Figure 5 Pore pressures and effective stresses in 20 days
(17)

式中c为固结系数, c =kK, k为渗透率, K为压缩模量。

σz =1MPa, k =0. 02米/天, K =107MPa, 计算了底部z =10米处孔隙压力随时间的变化及固结20天后, 孔隙压力和有效应力的分布(图 34)。

图 4图 5可以看出, 随着时间的增加, 孔隙压力减小。从分布来看, 模型的上部靠近透水部位孔隙压力小, 下部大。有效应力的分布与之相反。ANSYS简便方法计算结果与解析解吻合很好。

(2) 油田油井抽油时的压力下降和停产后井底压力的恢复

油井生产时, 井底压力随开采过程下降以及停产后井压恢复的规律对油田生产有重要意义。这也是流固耦合的渗流问题。我们对胜利油田某采油厂的实例, 利用二维问题求解。该问题的边界条件可使总应力不变, 可利用ANSYS线弹性模块求解总应力, 由式(7)求出加载瞬时的孔隙压力, 并以此压力做为初始孔隙压力, 利用式(9)求出孔隙压力随时间的变化。

由ANSYS软件计算的地层压力, 即初始孔隙压力为23MPa, 与实测的初始孔隙压力一致。由ANSYS软件的热传导方程求解的油井生产时的井底压力下降曲线和停产后井底压力恢复曲线如图 6。图中给出了实测结果。可以看出, 计算曲线与实测曲线吻合较好。

图 6 井底压力下降曲线和停产后井底压力恢复曲线 Figure 6 Bottom-hole pressure decline curve during oil extraction and bottom-hole build-up curve after oil production stop

图 6指出: (1)油井生产后, 井压随时间下降, 开始时刻井压下降较快, 随着时间的增加, 下降速度变缓, 曲线趋于平缓。(2)油井停产后, 压力开始恢复, 开始时恢复较快, 后来曲线变缓, 逐步恢复到开采前水平。

7 结论

(1) 在一维情况下, 可以用ANSYS软件中的结构模拟和热传导模块求解流体耦合方程。

(2) 在二维和三维情况下, 如果其边界条件使得任一方向上的变形受到约束。总应力随时间而变化时, 不能用ANSYS软件的方法求解流固耦合问题。其他边界条件下, 总应力不随时间而变化, 可以用ANSYS软件求解。

(3) 用ANSYS软件求解, 就是当总应力不随时间而变化时, 可用线弹性方程(或弹塑性方程)求解总应力和孔隙压力, 即初始孔隙压力。用ANSYS软件的热传导模块求解渗流方程, 得出流体压力、压力梯度、流速、流量等随时间的变化及空间分布。

(4) 用ANSYS软件求解的流固耦合方程, 其解答具有足够的精度, 方法简单、快速。是求解流固耦合方程的一种有效方法。

参考文献
[1]
Jaeger JC, Cook NG. W. 岩石力学基础. 中国科学院工程力学研究所译. 科学出版社, 1981. 253~ 267.
[2]
杨天鸿, 唐春安, 徐涛, 等. 岩石破裂过程的渗流特性[M]. 科学出版社, 2004: 47-48.
[3]
SMITH. M. 有限元方法编程[M]. 王菘译, 北京: 电子工业出版社, 2003. 301~ 325.
[4]
华东水力学院土力学研究室. 土工原理与计算[M]. 北京: 水力电力出版社, 1985: 215-218.
[5]
王连捷, 王薇, 张利容. 地应力驱动油气运移基本方程及有限元模拟[J]. 地质力学学报, 1999(3): 27-31.
[6]
王国强. 实用工程数值模拟技术及基础在ANSYS上的实践[M]. 西北工业大学出版社, 1999: 30-33.