高温岩体地热开发系统热-水-力耦合模型的三维有限元分析

高温岩体地热开发系统热-水-力耦合模型的三维有限元分析

人类面临两个主题：能源和环境。新能源的研究和开发是能源发展的基本途径。干净的地热资源的开发和利用得到了国内外的广泛认可和广泛关注。20世纪70年代，美国洛斯阿拉莫斯国家实验室的科学家提出了从深度到内陆的热能量回蓄能力的假设。美国、德国等国家进行了大量的长期实验。一些高温岩石地热电站已经建成。中国在高温水域的研究还处于初步阶段。赵阳生在高温水体的热破裂和固流热耦合方面做了一些研究工作，康健等人用随机介质理论建立了岩石热裂隙数值模型和高温岩体热耦合模型。高温水体的热合作用是一个渐进的过程，它是从体平衡中借用多个激励的提议。自那以后，许多模型都出现了。例如，简单的热传导模型、水流动模型、水热耦合模型、热水力耦合模型、岩石介质的性质，包括均匀模型和裂缝介质模型。赵阳生和王瑞峰在国内建立了高温水体热合作用的块裂介质的固流热耦合数学模型。该模型根据离散介质的假设研究了热液流在干裂缝中的固流热耦合模型。然而，除了用提升方法形成的低序列小裂缝外，许多低序列小裂缝的基岩的渗透性较弱，但小裂缝分布密度大，供水性功能不足。为了研究高温水体的热态平衡系统，作者提出了双介质的理论，将人工热储存层的树干裂缝视为裂缝介质，并将许多低序列小裂缝的基岩视为具有适应空腔和双间隙的连续介质。根据双重介质（形状连续介质和裂缝介质）的假设，建立了双介质高温水体的热-水-力耦合模型，以研究高温水体的复杂热液流规律。1双介质热水力耦合数学模型1.1质单元模型1)将裂隙岩体(人工热储层)看作由主干裂隙介质和含众多的低序次的小裂隙的拟连续介质组成,裂缝介质服从节理单元模型,拟连续介质为均质同性的弹塑性体并符合应变叠加原理.2)渗流规律服从非线性达西定律.3)由于高温高压作用下,岩体为单相水所饱和.4)热质在固相、液相介质中传递方式以传导和强迫对流为主,各介质的比热容及热传导系数不随温度和压力而变化.5)拟连续介质的变形以空隙变形为主.1.2裂隙介质的渗流方程拟连续介质的渗流方程div(ρwqi)+ϕ∂ρw∂t+ρw∂e∂t=0,(1)ρw=ρw0[1-βΤ(Τ-Τ0)+βp(p-p0)],(2)qi=kij∂p∂xj,(3)kij=k0exp(-γΘe);(4)div(ρwqi)+ϕ∂ρw∂t+ρw∂e∂t=0,(1)ρw=ρw0[1−βT(T−T0)+βp(p−p0)],(2)qi=kij∂p∂xj,(3)kij=k0exp(−γΘe);(4)裂隙介质的渗流方程qi=kfi∂p∂si,(5)kfi=gb212ν(Τ),(6)ν(Τ)=0.1775/(1+0.037Τ+0.000221Τ2),(7)qi=kfi∂p∂si,(5)kfi=gb212ν(T),(6)ν(T)=0.1775/(1+0.037T+0.000221T2),(7)式中:kij,kfi为不同介质的渗透系数;e为体积变形;Θe为有效体积应力;ν为运动黏滞系数;b为裂隙宽度;ρw,ρw0分别为T与T0时刻水的密度;βp,βT分别为流体压缩系数和热膨胀系数;ϕ为空隙率;qi为比流量;p为水压.1.3连续介质和能量的迁移据热力学第一定律:外界传入系统的能量与内部热源产生的能量之和等于物质内能的增量与对外做功之和.建立分别考虑岩体(固相)与温度场、液相与温度场耦合的能量平衡方程,该方程了考虑液体热对流(强迫对流为主)和固相热传导而导致的能量迁移.拟连续介质ρrcpr∂Τr∂t-λrΤr,ii+ϕcpwρwkiip,iiΤwi+ϕΤwβΤβpk,iip,ii+12(1-ϕ)Τrβ∂∂t(μi,j+μj,i)δij=0.(8)ρrcpr∂Tr∂t−λrTr,ii+ϕcpwρwkiip,iiTwi+ϕTwβTβpk,iip,ii+12(1−ϕ)Trβ∂∂t(μi,j+μj,i)δij=0.(8)裂缝中主要介质为水,能量的迁移是通过裂隙水流的热传导、热对流(强迫对流为主)及岩体裂隙表面和裂隙水流之间的热传导的双过程实现的∂(ρwcpwΤw)∂t-λw∇2Τw+(ρwcpwΤwkfip,i)i-2λrb(Τrb-Τw)=0.(9)∂(ρwcpwTw)∂t−λw∇2Tw+(ρwcpwTwkfip,i)i−2λrb(Trb−Tw)=0.(9)1.4裂隙介质的力学方程拟连续介质的力学方程(λ+μ)uj,ji+μui,jj+Fi+(αpδij)i-(3λ+2μ)αsδijΤi=0,(10)(λ+μ)uj,ji+μui,jj+Fi+(αpδij)i−(3λ+2μ)αsδijTi=0,(10)裂隙介质的力学方程σn=knεn,σs=ksεs,(11)σ′n=σn-p-βΤδijΤi,σ′s=σs,(12)σn=knεn,σs=ksεs,(11)σ′n=σn−p−βTδijTi,σ′s=σs,(12)式中:λ,μ为拉梅常数;u为岩体位移;F为体力;αs为岩体热膨胀系数;Tr,Tw分别为岩体和水的温度;σ′n,σ′s为裂缝法向和切向有效应力;Kn,Ks为裂缝法向与切向刚度;εn,εs为裂缝法向与切向变形;λr,λw为岩体和水的热传导系数;ρr为岩体的密度;cpr,cpw分别为岩体和水的热容系数;Trb为岩体裂缝边缘的温度;α为等效孔隙压系数,0<α<1;s为切向座标.以上3组控制方程,再辅以初始和边界条件,构成了基于双重介质假设的裂隙岩体热-水-力耦合数学模型.2热水力耦合数学模型的数值解决方案2.1外载荷项t[Κ]{U}={F},(13)式中:[K]为总体刚度矩阵;{U}为位移矩阵;{F}分别为体积力、表面力、水压梯度、温度梯度引起的外载荷项.{F}={Fp}+{Fv}+{Fs}+{Ft}=n∑i=1{∫v[Ν]Τ{{Q}+{Fij}}dv+∫s[Ν]Τ{q}ds+∫ΩBΤDε0dΩ},(14)ε0=(3λ+2μ)αs(Τ-Τ0)×Τ,(15)式中:{Fp},{Fv},{Fs},{Ft}分别为体积力、表面力、水压梯度、温度梯度引起的外载荷项;[N]为插值形函数;[D]为弹性方程;Q为体力;q为面力;Fjj为水压梯度和温度梯度引起的附加荷载顶;αs为岩体的热膨胀系数.对(αpδij),i项的处理方法:渗透压力p作为均布体积力处理.按如下方法把渗透压力转化为等效节点载荷{Feij}=∫Ωe∂Νi∂xiαpδijdΩ=∫1-1∫1-1∫1-1∂Νi∂xiαpδij|J|dξdηdζ.(16)2.2单元排放系数法[Τ]{Ρ}+[S]∂Ρ∂t+{Η}=0.(17)双重介质热-水-力耦合模型中水-力耦合体现在介质的渗透系数是体积应力、裂缝宽度的函数,水-热耦合体现在流体的密度、黏度是温度的函数和水的热膨胀性上.拟连续介质传导矩阵[Τ]e=∫Ωe[B]Τ[Κ][B]|J|dΩ.(18)裂隙介质传导矩阵[Τ]e=gb312v∫Ωe[B]Τ[B]|J|dΩ.(19)单元储量矩阵[S]e=∫ΩeΦβΝiΝjdΩ.(20)单元列矢量[Η]e=∫ΩeΝi∂e∂tdΩ+∫ΓegΝidΓ,(21)∂e∂t=Θ(t+Δt)-Θ(t)Δt.(22)2.3双重介质热-水-力耦合模型的建立[R]{Τ}+[G]∂Τ∂t={Ρ},(23)式中:[R]是热容矩阵,[G]是热传导矩阵,{P}是温度载荷矩阵,{T}是结点温度列阵,温度对时间的离散采取隐式差分的形式.拟连续介质Geij=∫ΩeCprρrNiNjdΩ,(24)Reij=-∫Ωeλr(∂Νi∂x∂Νj∂x+∂Νi∂x∂Νj∂y+∂Νi∂z∂Νj∂z)dΩ+∫ΩeϕCpwρwΚiiΝj×(∂Νi∂x∂p∂x+∂Νi∂y∂p∂y+∂Νi∂z∂p∂z)dΩ,(25)Ρei=-∫Ωe12Νi(1-ϕ)β×∂(μi,j+μj,i)∂tδijΤr(t)dΩ.(26)裂隙介质Reij=-∫Ωeλw(∂Νi∂x∂Νj∂x+∂Νi∂x∂Νj∂y)dΩ+∫ΩeCpwρwΚfiΝj(∂Νi∂x∂p∂x+∂Νi∂y∂p∂y)dΩ,(27)Geij=∫ΩeCpwρwΝiΝjdΩ,(28)Ρei=∫ΩeΝi2λrb(Τrb-Τw(t))dΩ.(29)在双重介质热-水-力三维耦合计算中,先分别计算出裂隙介质和拟连续介质的有限元单元矩阵,然后将两类模型的有限元单元矩阵叠加,组装成同一总体静力离散方程、渗流离散方程、能量守恒离散方程,在耦合计算过程中,不同介质的固体变形、能量方程和渗流看作独立的子系统,具有耦合效应的参数值(温度、应力、水压、渗透率、密度、黏度、裂缝宽度等)在各子系统相互传递耦合求解.对时间域的离散采取用变时步隐式差分格式.开发了双重介质热-水-力耦合模型的三维有限元正分析程序DTHM.for.3模拟应用原则的高温岩体取样值3.1围岩应力方向x方向本文以云南腾冲热海热田为例,选一高温岩体地热区进行分析,模拟研究区域:地表4km以下,5.0km×5.0km×5.0km区域内的高温岩体地热开采过程,如图1所示.采用水平井、垂直裂缝的计算模型,如图1,2:沿垂直最小主应力方向(X方向)设置2条平行的主干裂缝,其坐标位置为x1=2.250km,x2=2.750km,即主干裂缝垂直间距500m,注入井和生产井垂直于裂缝而布置,岩体地热433K,应力边界条件如图1所示.计算参数见表1,在高温岩体地热开发系统中,热能的运移主要通过裂缝与两井的连接平面而发生,所以选取连接注入井与生产井的剖面为研究对象,且在此剖面上选择裂缝面点1#～11#,如图2.渗流场边界条件:在地表注入井给定压力27.3MPa,生产井压力9.8MPa,由于静水压力作用在图2所示的地下-6km左右的水平注入井将产生84.8MPa的注入水压,水平生产井会产出62.3MPa的产出水压.温度场边界条件:由于模型周围仍然是很大规模的热储,可以认为有热量补给,取热流边界条件,即:注入水给定温度30℃;地热梯度为50℃/km.3.2数值模拟和分析1采动覆岩裂缝宽度拉伸温度注入的低温水使得注入井附近围岩温度下降,围岩发生收缩变形,岩体应力场的重新分布,同时岩体渗流场扰动也将改变岩体应力场分布,裂隙单元形成应力降低区,系统运行3a后,注入井附近裂隙内应力约为104MPa.随开采时间的延长,岩体应力以注入井为中心的低应力区会逐渐扩大.单元的渗透系数随单元应力的降低而增大,因此低应力区的扩张将有利于于水的渗流.裂缝宽度变化总体上呈非线性增加的趋势,系统运行2.5a之内,各裂缝结点宽度急骤增长.注入井附近11#结点宽度达到162μm,其后,增长率逐渐降低.3a以后,11#,9#,7#,5#结点宽度几乎趋于稳定,而生产井附近1#,3#结点仍然随着岩体温度的降低而逐渐加宽,到9a左右与中段区域5#结点趋于一致.由此可见,随开采时间的延长,裂缝宽度会随着地热提取而增加,裂缝渗透阻力下降,渗透系数增加,在系统运行过程中适当可以降低注水压力,减少能耗.2主采系统稳定后的水压分布注入的低温水会在高水压梯度的作用下,沿人工热储层流向生产井的低水压区,在开始运行阶段,水大部分通过低应力区、强渗透性的主干裂缝而流动,主干裂缝内水压迅速上升,而随着渗流的发生,基岩水压逐渐滞后升高,在大约半年时间,几乎与主干裂缝水压趋于一致.1a左右,整个系统水压保持稳定,以固定的压力梯度循环流动,满足高温岩体地热开发系统的水循环,系统稳定后的水压分布如图5.图6给出了裂隙面各点水压随开采时间的变化规律,可见,在开采的前半年时间内,各点水压急骤增长,其后逐渐变得平缓,运行到1a后,系统保持稳定.结合图5,6可以看出,系统稳定渗流以后,在靠近注入井与生产井的区域压力梯度大,而在中间区域则比较小,固定井口从11#～8#,4#～1#压力降低约10MPa,平均水压梯度达0.046MPa/m;而中间区域4#～7#压力降低约2.5MPa,平均水压梯度约0.011MPa/m,这说明:由于固定井口压力的影响,使得井口附近区域水压下降快,水流速度快,有利于热能的输送.3不同开采时间下水温度变化规律在系统运行的开始,注入井周围岩体温度下降很快,形成一个相对低温区,并随开采时间的变化,系统温度逐渐降低,其低温区域逐渐扩大,在大约2a后,温度下降到364K并继续缓慢下降,到9a左右时间下降到330K,生产井区域573K的高温区,到系统运行9a后,只有405K.图7为9a时的系统温度场分布.图8为裂隙面各点的温度变化规律.9a内随开采时间温度整体呈下降趋势,在靠近注入井的区域内的点(6#～11#)其下降趋势呈时间的负指数关系,开始急骤下降后变得平缓.而在生产井附近区域1#～5#温度近乎线性下降.在系统运行大约一个季度,10#点温度下降很快,只有400K,并继续降低,而1#,2#,3#仍然保持570K的高温,在系统运行3a时间内,注入井附近区域温度下降的速率很快,远远超过生产井附近区域,而随开采时间的延长,注入井附近区域温度下降变的平缓,并逐渐与生产井趋于一致.可见随开采时间变化,其水温度升高的速率是逐渐减小的.4