多孔介质流固热三场全耦合数学模型及数值模拟.docx
- 文档编号:16755726
- 上传时间:2023-07-17
- 格式:DOCX
- 页数:20
- 大小:63.24KB
多孔介质流固热三场全耦合数学模型及数值模拟.docx
《多孔介质流固热三场全耦合数学模型及数值模拟.docx》由会员分享,可在线阅读,更多相关《多孔介质流固热三场全耦合数学模型及数值模拟.docx(20页珍藏版)》请在冰点文库上搜索。
多孔介质流固热三场全耦合数学模型及数值模拟
第25卷增1岩石力学与工程学报Vol.25Supp.12006年2月ChineseJournalofRockMechanicsandEngineeringFeb.,2006
多孔介质流–固–热三场全耦合
数学模型及数值模拟
盛金昌
(河海大学水利水电工程学院,江苏南京210098
摘要:
给出多孔岩体介质的流–固–热三场全耦合数学模型。
假定流体为单相流,固体介质为非沸腾的饱和、热
弹性多孔介质,该模型由流体物质守恒方程、力学平衡方程和能量守恒方程这3个相互耦合的方程组成,其中包
含了众多耦合项,并定义一系列的本构关系及耦合变量。
以FEMLAB工具为基础,将该数学模型转化成为一个统
一的偏微分方程组,在人机交互的环境下,实现流–固–热三场全耦合数值求解,一次解出渗流场、位移场和温
度场,给出更接近真实物理过程的数值解答,避免松散耦合法求解多场耦合问题带来的误差。
利用一个已知解析
解和数值解的算例来证明了耦合模型及求解方法的正确性。
最后模拟通过井孔向岩体中注入冷水时流–固–热
全耦合过程,详细地分析全耦合作用对井壁围岩应力的影响,计算结果表明:
流–固–热三场耦合作用对井壁的
稳定分析有非常重要的影响。
关键词:
岩石力学;多孔介质;流–固–热三场全耦合模型;数值模拟
中图分类号:
TU457文献标识码:
A文章编号:
1000–6915(2006增1–3028–06
FULLYCOUPLEDTHERMO-HYDRO-MECHANICALMODELOFSATURATEDPOROUSMEDIAANDNUMERICALMODELLING
SHENGJinchang
(CollegeofWaterConservancyandHydropowerEngineering,HohaiUniversity,Nanjing,Jiangsu210098,China
Abstract:
AfullycoupledThermo-Hydro-Mechanicalmodelispresentedinasingle-phase,non-boilinglinearthermoelasticmedium,whichincorporatescross-coupledfluidflowequation,energyconservationequationandmechanicalequilibriumequationwithmanycross-couplingterms.Aseriesofconstitutiverelationsandcross-couplingrelationsbetweenmaterialpropertiesandindependentvariablesaredefinedinthemodel.ThecoupledmultiphysicsmodelissimultaneouslysimulatedbyusingFEMLAB,thefirstengineeringtoolthatperformspartialdifferentialequation-basedmultiphysicsmodelinginaninteractiveenvironment,whichthemathematicalmodelistranslatedintoasetofpartialdifferentialequations.Theporepressure,displacementsandtemperature,whichshouldtheoreticallyapproachthemostrealisticresults,canbesolvedsimultaneouslybyusingFEMLAB,inwhichtheerrorsinothercouplingalgorithmscanbeavoided.Anexamplewithknownanalyticalandnumericalresultsisusedtovalidatethemultiphysicsmodel.Inparticular,coldwaterinjectionintowellboreismodeledwithrealistictimestepsindicatingthatthecoupledprocesseshavesignificanteffectsonthestressesofboreholewallandwellborestability.
Keywords:
rockmechanics;porousmedium;fullycoupledThermo-Hydro-Mechanicalmodel;numericalsimulation
收稿日期:
2005–06–13;修回日期:
2005–11–15
基金项目:
国家自然科学基金重点项目(50539030
第25卷增1盛金昌.多孔介质流–固–热三场全耦合数学模型及数值模拟•3029•
1引言
在众多的工程领域中,多孔岩体介质中的流、固、热多场耦合问题具有非常重要的意义和作用。
例如:
在油气田开采、CO2的地下储存、地热资源的开采、核废料的地下储存等领域中,都存在着多场耦合问题,多物理场之间的复杂耦合作用对工程具有非常重要的影响。
在这些耦合作用中,既包括线性物理作用,也包括材料的非线性耦合作用。
目前,国外在多场耦合分析方面已进行了大量的工作,如进行了两场(变形和流动或三场的耦合分析[1~3],其理论一般都假设岩体为多孔弹性介质。
这个耦合系统包括3个相互耦合的方程所组成的方程组以及一组物理参数之间的耦合关系,如:
饱和度与渗透系数之间的关系、应力与渗透系数之间的关系等。
地质体的耦合作用过程吸引了广泛的研究兴趣,目前在流、固、热等多场耦合分析方面已有不少研究成果[4~8]。
多孔岩体介质中的流–固–热多场耦合问题非常复杂,数值求解具有很高的难度。
求解多场耦合问题时多采用交叉迭代求解。
耦合问题求解有3种基本算法[3]:
单向耦合算法、松散耦合算法和全耦合算法。
对于单向耦合算法,即2组独立的方程在同一时间步内分开求解,求解时只是将其中的一个物理过程的计算结果作为另一个物理过程的输入,这种传递只是单向的。
例如:
由流动方程解出的孔隙压力作为荷载传给力学计算来求解应力和位移。
J.Fredrich等[4]的工作就是采用这一算法。
对于松散耦合算法,两组方程独立求解(和单向耦合算法一样,但是有关信息在指定的时间步内在两个求解器之间双向传递。
这一算法的优点是相对容易实现,而且还能反映较复杂的非线性物理过程,它接近全耦合算法[1,3,5~8]。
对于全耦合算法,需要推导出统一的一组全耦合方程组(通常是一个大型的非线性全耦合的偏微分方程组,这里面融合了所有的相关物理过程。
求解多物理耦合问题应该首选全耦合算法,因为在理论上它能给出最真实的结果。
然而由于同时求解多物理耦合过程异常复杂,到目前为止还没有这方面令人信服的研究成果。
本文的主要工作如下:
在J.Noorishad和C.F.Tsang[1]研究的基础上,首先建立一组多孔岩体介质瞬态的流、固、热三场全耦合方程组(其中包括固体变形、能量传输、流体流动3个相互耦合的过程,多物理场之间的交叉耦合还包括材料性质与独立变量之间的耦合关系,然后利用FEMLAB软件(基于偏微分方程组的多物理过程模拟工具作为平台,成功地求解了这一全耦合偏微分方程组,避免了松散耦合法求解多场耦合问题带来的误差,实现了同时求解多物理场耦合过程。
本文的求解方法是一个全耦合算法,在理论上它能给出最真实的结果。
通过对一个具有已知解析解和数值解的算例的计算分析来证明本文耦合模型及求解方法的正确性:
一维砂柱的等温固结和非等温固结问题。
最后模拟分析了通过井孔向岩体中注入冷水时流、固、热全耦合过程,详细地分析了全耦合作用对井壁围岩应力的影响。
2流–固–热三场全耦合数学模型
在J.Noorishad和C.F.Tsang[1]工作的基础上,下面给出流–固–热三场全耦合数学模型。
假设岩体介质为饱和的多孔弹性介质,并假设REV存在。
2.1流动方程
流体流动时质量平衡方程的欧拉形式如下:
Q
V
t
=
∇
+
∂
∂
l
l
l
(
φρ
φρ
(1
式中:
φ为岩体孔隙率,
l
ρ为流体的密度,t为时
间,
l
V为流体速度矢量,Q为流体的源汇项。
根据流体流动的动量方程可得Darcy定律:
(
l
l
l
g
P
k
Vρ
µ
−
∇
−
=(2
式中:
l
µ为流体的动力黏滞系数,k为孔隙介质的渗透率,P为孔隙压力,g为重力加速度矢量。
将式(2代入式(1,并加上固体骨架的变形项,经过一系列推导可得
−
∂
∂
+
∂
∂
+
∂
∂
t
T
t
P
tt
p
v
lφβ
φβ
ε
ρ
ρ
Q
g
P
k
=
⎥
⎦
⎤
⎢
⎣
⎡
−
∇
∇
(
l
l
lρ
µ
ρ
ρ
(3
式中:
ρ为流体的参考密度,
v
ε为岩体的体积应变,
T为温度,
t
β为流体的热体积膨胀系数,
p
β为流体的压缩系数。
式(3即为多孔热弹性介质中流动流体的控制方程的最后形式。
式(3等号左边的前3项分别表示
·3030·岩石力学与工程学报2006年
由于应变、流体压力和温度所引起的流体体积的变化,等号左边的最后1项代表由压力梯度和重力作用而引起的流体流量。
2.2能量守恒方程
固体骨架和流体共同存在于同一个体积空间,但它们具有不同的热动力学特性:
如比热容和热传导系数等。
因此固体骨架和流体的能量守恒方程需要分别定义。
固体骨架的能量守恒方程定义如下:
sssp1((1((1(qTKtT
cφφρφ−+∇∇−=∂∂−(4
式中:
sp(cρ为岩体骨架的热容,sK为岩体骨架的热传导张量,sq为岩体的热源强度。
对于流体,相应的能量守恒方程可定义如下:
llllplp((((qTKTVctT
cφφρρφ+∇∇=∇⋅+∂∂(5
式中:
lp(cρ,lK和lq分别为流体的热容、热传导张量和热源强度。
对于单相流,假设固体和流体之间总是处于热平衡状态,这样将式(4,(5迭加,并考虑到变形能,即可得到以下统一的能量守恒方程[9]
:
=∇⋅+∂∂−+∂∂TVct
TtTcpp((1((llv0tρεγφρ
tt(qTK+∇⋅⋅∇(6式中:
βλµγ32(+=,µ和λ为拉梅常数,β为各向同性固体的线性热膨胀系数;0T为无应力状态下的绝对温度;tq为充满了流体的多孔介质的热源汇项,定义为slt1(qqqφφ−+=。
tp(cρ和tK分别为充满了流体的多孔介质的比热容和热传导系数,且
splptp(1(((cccρφρφρ−+=
slt1(KKKφφ−+=2.3力学平衡方程
假设岩体为理想热弹性体,考虑流体的孔隙压力和热应力的本构关系如下:
PTijijklklijijijαδγδεδλδµεσ−−+=2(7
式中:
ijσ为应力分量,ijε为应变分量,α为Biot系数,ijδ为Kronecker数。
将式(7和应变位移关系代入静力平衡方程式,可以得到用位移表示的包含耦合项的修正的Navier平衡方程:
iiijijjjiFTPuu=−−+,,,,γαλµ2/((8
式中:
iF为体积力分量。
式(3,(6和(8一起完整地定义了岩体三场全耦合作用的数学模型。
施加一定的边界条件和初始条件,即可求解上述耦合控制方程组。
另外,渗透系数与应力之间的关系有很多表达式,本文采用以下关系:
⎪
⎭
⎪
⎬⎫
+=−=(21
exp(21vv0σσαpkk(9
3基于FEMLAB全耦合分析的实施
FEMLAB是一个在交互式环境下求解基于偏微分方程组的多物理耦合过程的有限元分析工具[10],它的功能非常强大,可以应用于广泛的科学和工程领域。
在FEMLAB中,不再需要编制复杂的偏微分方程组的求解器,它含有一些内嵌的物理模型,如:
化学工程模型、电磁学模型、结构力学模型等。
但是其中功能最强大、而且也是最灵活的还是它的偏微分方程组模式。
描述偏微分方程组有3个数学应用模式:
系数形式(coefficientform、通式
(generalform和弱形式(weakform。
其中弱形式最为灵活,但对大多数问题,通式模式已足够强大。
本文选择结构力学模型加通式模式微分方程组
(描述耦合热传导方程和流体流动方程来求解岩体三场全耦合问题,当然在结构力学模型和通式微分方程组模式之间还需定义很多交叉耦合项。
在这些问题正确定义出来之后,FEMLAB在求解时,将先把结构力学模型和通式微分方程组模式结合在一起转换成一个统一的通式形式的微分方程组,然后统一求解这个总的通式形式的微分方程组,同时解出位移场、渗流场和温度场,从而实现了三场的全耦合求解,避免了松散耦合法求解多场耦合问题带来的误差,给出了更接近真实物理过程的数值解答。
求解结果可以用多种方式来表达,如等势线、曲线、图像及动画等。
4流–固–热全耦合模型的验证
下面通过算例证明本文的三场全耦合数学模型以及基于FEMLAB的耦合数值模拟器的正确性和有效性。
M.A.Biot[2]利用一个有限元耦合分析程序
第25卷增1盛金昌.多孔介质流–固–热三场全耦合数学模型及数值模拟•3031•
(ROCMAS分析了一个砂柱的热弹性固结问题,这是一个典型的流固热耦合问题,有关研究[1
,2]
分别
给出了等温弹性固结解析解和非等温弹性固结的数值解。
所用的物理参数见表1。
用基于FEMLAB的耦合模型求解了这两个固结过程。
图1给出了本文解答与数值解(非等温弹性固结过程和Biot解析解
[2]
(等热固结过程的对比情况,由图可知,本文的解
答和相关研究[1
,2]
的结果非常接近,这说明本文的
数值模型非常精确可靠。
表1验证算例中所用的岩体材料参数
Table1Rockmassparametersusedinvalidationexample
弹性模量E/Pa泊松比νBiot系数α孔隙率φ渗透率k/m26.0×1030.4
1.00.2
4.0×10
-6
热传导率Kt/(kJ·m-
1·℃-
1
比热容(ρct/(kJ·m
-3
·℃-
1
线性热膨胀系数
β/℃
-18.36×10
-1
167.0
3.0×10
-
7
图中黑点代表本文模拟结果,实线代表已知解
图1本文模型的模拟结果与已知解析解和数值解的对比Fig.1Comparisonsofsimulationresultswiththeknown
solutions
5应用实例分析
在石油和天然气开采中,钻井过程中井壁稳定问题非常重要,实际上这一过程是一个流–固–热多场耦合问题。
注入液体的压力和温度将改变井壁围岩的应力分布,从而最终影响井壁的稳定[11]。
下面利用本文的全耦合模型求解该问题。
5.1模型的描述
图2给出了该算例的计算域和边界条件,这里计算域为10m×10m,井孔半径为0.1m。
水平和竖
图2计算域及边界条件示意图
Fig.2Schematicsimulationmodelanditsboundaryconditions
直方向的边界压力分别为30和20MPa。
假设注入液体压力比岩体中的大10MPa,注入液体是冷水,比岩体温度低50℃,液体温度低有利于井壁的稳定。
假设井壁孔隙压力和温度与注入液体的压力和温度完全相同。
计算域中的初始压力和温度差为0。
表2所示为应用实例中所用的岩体材料参数。
5.2计算工况
本文通过以下3种耦合和不耦合计算工况来说明耦合分析对井壁稳定分析的重要性,分别为:
工况1:
流–固–热全耦合计算分析;工况2:
把岩体看成不透水介质,即不考虑流动过程,但计算应力时把水压力加在井壁上。
计入温
度场引起的热应力;
工况3:
不考虑温度场和流场之间的耦合作用,即不考虑对流项对温度的影响和温度变化引起水压力的变化。
但考虑水压力场和热应力对变形的影响。
在FEMLAB求解时,本文选用四阶Lagrange型三角形单元。
下面本文将主要分析三场耦合作用对孔周围岩应力分布的影响。
5.3切向应力的对比分析
图3(a给出了在时刻为86400s时A-A剖面(见图2上不同工况下的切向应力对比图。
由图3(a可知,在井壁上,用全耦合算法得出的切向压应力(压为负值最大,
工况2的切向压应力值最小。
在点P,工况1的切向应力值为-25.14MPa,工况2的为-
17.78MPa,而工况3的为-23.79MPa,即不考虑耦合时点P的切向应力值比考虑耦合情况下的值小很多。
在岩体内部(离井壁距离大于0.2m,3种工况计算出来的结果逐渐趋向一致,它们之间的差别大
0.0
1
-0.1-0.2-0.3-0.4-0.5-0.6
2
3
4
5
6
固结时间/s
顶部沉降量/mm
·3032·岩石力学与工程学报2006年
表2应用实例中所用的岩体材料参数
Table2Rockmassparametersusedinapplicationexample
密度ρs/(kg·m-3
弹性模量E/MPa
热传导率Kt/
(kJ·m-1s-1·℃-
1热膨胀系数β/(℃-
1
比热容cps
/(kJ·kg-1
·
℃-
1
孔隙率φ渗透率k/m
2
Biot系数α泊松比ν密度ρl/(kg·m-3
压缩系数βp/
(GPa-
1动力黏滞系数(20℃η/
(N·s·m-
2热膨胀
系数βt/(℃-
1
比热容cpl/
(kJ·kg-
1·℃-
1
2.6×1035000
3.08×10
-3
1.0×10
-5
0.84
0.2
5.0×10
-13
1.0
0.25
1000.0
0.513
10
-3
3.17×10
-4
4.2
A-A剖面上自孔壁的距离/m(a孔周切向应力分布
A-A剖面上自孔壁的距离/m(b孔周径向应力分布
图3各种计算工况下孔周应力沿A-A剖面的分布Fig.3VariationofstressesalongA-Asectionfordifferent
cases
致在2.0MPa以内。
由此可知,考虑流动过程与热传导过程之间的耦合和不考虑它们之间的耦合对切向应力的影响。
图4(a给出了井壁上点A切向应力随时间的变化过程,图4中1,2,3的标注含义同图3。
由图4(a可知:
工况2和工况3条件下,该点的切向应力值随时间变化很小,而工况1条件下,点A切向应力则随着时间的增加而有一定幅度的减小。
这表明,考虑耦合作用后,井壁切向应力随着流动和热传导过程有一定的调整。
t/(104s
(a点A切向应力随时间的变化过程
t/(104
s
(b点B切向应力随时间的变化过程
图4各种计算工况下孔周应力随时间的变化过程Fig.4Variationofstresseswithtimefordifferentcases
5.4径向应力的对比分析
图3(b给出了在时刻为86400s时A-A剖面上不同工况下的径向应力对比图。
由图3可知,3种工况得出的径向应力的大小和分布大致相同,只是工况2(不耦合情况值要小一些,大约小1~2MPa。
到岩体内部后(大于0.2m,工况2比工况3的稍小一点。
在点B(岩体内部,工况1的径向应力值为
-24.66MPa,工况
2的为-23.38MPa,而工况3的
为-24.98MPa,可见3种工况下得出的岩体内部径向应力值差别不大。
图4(b给出了岩体内部点B径向应力随时间的变化过程,由图可知:
该点的径向应力在3种工况之间的差别随着流动和热传导过程的发展,它们之间的差别越来越大。
工况1条件下点B的径向应力
1—流–固–热全耦合过程
2—只有固–热耦合,岩体为不渗透介质,没有
流动过程,水压力作为面力加在井壁上3—只有流–热不耦合
切向应力/MPa
径向应力/MPa
切向应力/MPa
径向应力/MPa
第25卷增1盛金昌.多孔介质流–固–热三场全耦合数学模型及数值模拟3033在前5000s随时间的增加而减小,而此后则随时间增加而增加,这主要是因为热传导过程比较快,在前5000s时间内热应力对该点的径向应力的影响起主要作用,而此后孔隙压力对该点的径向应力的影响占主导作用。
工况2条件下点B的径向应力单调增大,而工况条件下则单调减小。
不过它们之间的差别不是太大,大约在2MPa之内。
以上分析表明:
采用耦合模型和采用非耦合模型,得出的孔周应力的大小和随时间的变化有比较大的差别,只有采用耦合模型,得出的孔周应力才是真实可靠的,这对井壁的稳定分析有重要价值。
明:
只有采用耦合模型,得出的孔周应力才是真实可靠的,这对井壁的稳定分析具有重要的参考价值。
参考文献(References:
[1]NoorishadJ,TsangCF.Coupledthermohydroelasticityphenomenainvariablysaturatedfracturedporousrocks-formulationandnumericalsolution[A].In:
StephanssonO,JingL,TsangCFed.CoupledThermo-Hydro-MechanicalProcessesofFracturedMedia[C].Amsterdam:
ElsevierSciencePublishers,1996.93–134.[2]BiotMA.Generaltheoryofthree-dimensionalconsolidation[J].JournalofAppliedPhysics,1941,12:
144–164.[3]MinkoffSE,StoneCM,BryantS,etal.Coupledfluidflowandgeomechanicaldeformationmodelling[J].JournalofPetroleumScienceandEngineering,2003,38:
37–5
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 多孔 介质 流固热三场全 耦合 数学模型 数值 模拟