冷冲击作用下干热岩井井筒裂纹扩展数值模拟研究

摘要

干热岩（HDR）是一种清洁可再生能源，主要通过增强型地热系统（EGS）进行开发。在EGS工程中，无论是注入井还是生产井的建设都需要采用地热钻井技术，而高温高压钻井过程中地层破裂、井壁坍塌等现象是干热岩钻井施工面临的重要问题。在温差作用下，岩石的矿物颗粒由于其热膨胀、冷收缩特性的差异性，矿物颗粒之间产生温度应力，导致岩体发生热破裂。本文借助RFPA数值模拟软件，对干热岩开采过程中井筒花岗岩冷冲击作用下的裂纹扩展进行研究。结果表明：井筒模型在冷冲击过程中，岩石表面拉应力随着冷冲击时间的增加，先升高至峰值后缓慢下降。裂纹扩展可以大致分为前、中、后期。冷冲击前期，井筒周围出现环形拉应力区，开始出现均匀的微小裂纹。冷冲击中期，随着时间的增加，拉应力区逐渐向井筒外围扩展，裂纹随着拉应力区向外扩展。冷冲击后期，拉应力大小逐渐降低直至小于模型抗拉强度，裂纹扩展速度减缓直至停止扩展。围压、井径、温度对井筒围岩冷冲击时的破坏损伤效果影响显著，其中，温度对冷冲击裂纹扩展起到促进作用，围压对冷冲击裂纹扩展起到抑制作用，井径对冷冲击裂纹扩展起到促进作用。

关键词

干热岩; 井筒; 裂缝扩展; 数值模拟; 冷冲击; 增强型地热系统

0 引言

干热岩（HDR）一般是指在地球表面埋藏深度大于3000 m，内部不含或含少量流体，温度在180 ℃以上且渗透率极低的岩体^［1-2］。干热岩地热资源具有储量大、开采过程安全且受环境影响较小等优点。HDR储层岩石主要为变质岩和结晶岩，比较常见的有花岗岩、花岗闪长岩，具有低孔隙度、低渗透性、高强度的特点，因此HDR中热能的提取需要采用增强型地热系统（Enhanced Geothermal System， EGS），即通过水力压裂等工程手段在地层深处低渗透性干热岩体中形成人工地热储层，采出大量热能的人工地热系统^［3］。其工作原理是在干热岩中钻注水井，利用压裂液通过注水井对目标地层进行大排量压裂，采用加压致裂的方法建立裂隙体系即人工热储层。在储层中高温岩体与低温流体相互接触时产生较大热应力会使岩石产生裂隙，同时在高压力的作用下，裂隙会不断增多、延伸，相互贯通后形成裂隙网络。冷水流经热储层时渗入岩层裂隙吸收热量，与周围岩体发生热交换，会产生高温高压水和水汽，这些水和水汽通过生产井采出，用于地热发电和综合利用，冷却后的水可通过热交换系统进行循环利用^［4-5］。

在干热岩开采过程中时刻伴随着热交换过程，因热交换产生的热应力是导致岩石破裂的重要因素，在热应力的作用下岩体发生热破裂现象，会在一定程度上改变岩石内部孔隙结构，对岩体自身物理性质和岩体内的流体渗流产生影响^［6］。

在干热岩资源开采过程中，无论是开采井还是生产井的建设都需要采用地热钻井技术^［7］，而高温高压钻井过程中地层破碎、井壁坍塌等现象是干热岩钻井施工面临的重要问题。由于施工过程中泥浆的使用，使得井壁两侧岩石会持续接触低温流体^［8］，在钻进地层较深时，井壁岩石温度较高，此时与低温流体接触，会发生热交换作用导致井壁岩石温度发生变化，进而产生热应力，造成井壁岩石不同程度的损伤破坏，更易发生缩径和井筒坍塌事故。与此同时，在储层改造过程中，井筒周围岩石由于冷冲击作用产生的微裂纹可以提高压裂效率^［9-10］，提高造缝能力。因此研究井筒花岗岩冷冲击作用下的裂纹扩展对干热岩开采具有重要意义。

目前，利用数值模拟技术研究岩石热破裂问题已成为一种常用技术手段，并且优势明显，严成增^［11］在FDEM方法的基础上建立FDEM‑TH方法，对圆筒试样在两种不同温度边界条件下的热破裂情况进行分析。唐世斌等^［6］针对油气开发领域中低温液体的岩石破裂问题，运用热弹性理论获得平板模型，计算在模型遭受冷冲击时，其内部温度及应力分布理论解，结果表明，冷冲击边沿处的温度降低速度先快后慢，最终达到环境温度，拉应力也先急剧增加后逐渐降低。黄鑫等^［12］基于静态破裂剂破岩机理，利用数值模拟方法研究含孔岩样在热应力作用下的破坏过程，结果表明，在热应力和膨胀力的共同作用下，模型萌生两条垂直模型边缘的孔径方向的裂纹。包春燕^［13］建立层状模型，研究层状岩体的温度分布和裂纹扩张过程。高红梅等^［14］以岩石变形为切入点，研究热应力对花岗岩裂隙渗透系数和开度的关系，分析温度对水力开度的作用机理，数值模拟结果表明，温度越高热应力对花岗岩的渗透率影响越大。郑君等^［15］利用多场耦合数值模拟软件COMSOL，模拟二氧化碳爆破干热岩储层激发工艺的钻进过程，分析得出了冷却过程中炮孔附近储层的区域性温度分布特征及影响因素。肖鹏等^［16］利用COMSOL软件研究地热储层单裂隙岩体中渗流传热机理。研究发现，流体注入速度的提升会降低整个系统的寿命和出口法向总热通量；而流体注入温度的提升会增加系统的寿命和出口法向总热通量和总热量。Ali等^［17］利用离散元颗粒流数值模拟软件PFC软件，探究微波诱导的矿物颗粒热损伤的机理。模拟结果表明，岩石内部颗粒大小对数值模拟结果有很大的影响，使用颗粒流软件须在明确内部颗粒大小下进行参数校核和调整。Zhao等^［18］利用颗粒流数值模拟方法探究不同冷却方式对高温花岗岩巴西劈裂结果的影响。模拟结果表明随着模型所用颗粒尺寸的增大，常温花岗岩抗拉强度下降，而在花岗岩温度400 ℃以下，压缩强度无明显变化；在同一颗粒尺寸下，经过热处理的花岗岩的抗拉强度降低，并且温度越高，强度越低；同时随着颗粒尺寸的增加，产生的裂纹长度增加，但数量下降。沈宝堂等^［19］使用FRACOD软件，用迭代法进行力学数值计算和瞬态热流模拟，结果表明流体温度与围岩温度不同时，流体在裂隙中流动会引起热交换，在岩体内部产生温度压力，促进裂纹扩展。

对于岩石热破裂的问题，前人已经做了很多的研究，但目前的研究以热冲击作用为主，而冷冲击研究偏少，常以温度作用后的破坏结果为主要结论，缺少对微观损伤和裂隙扩展过程的展示，未能充分阐述岩石变形破坏行为及过程。此外，目前的研究主要考虑温度作用对岩石物理力学性质的影响，欠缺热力学参数对岩石冷冲击过程的温度、应力变化和破坏过程的深入探究。基于此，本文在考虑EGS实际工况的条件下，以传热学理论为基础，借助RFPA数值模拟软件，构建井筒冷冲击模型，探究冷冲击作用下井筒岩石的应力变化情况。研究冷冲击作用下井筒裂纹扩展过程，分析井筒围压、井径、围岩温度对井筒裂纹扩展的影响，揭示干热岩开采中井筒冷冲击裂纹扩展机理，为干热岩资源的开发利用提供理论参考。

1 本构方程

弹性材料的温度发生变化时，其各个位置都会随着温度变化而发生膨胀或收缩，由于弹性体受外界约束、材料单元间相互约束，膨胀或收缩又不能自由发生，因此产生温度应力，所以按照热弹性力学求出单元各点温度应力，首先要计算热传递和热传导公式，依据热力学知，单位条件下热力耦合模型热传导公式为^［20］：

$$\frac{\partial }{\partial x}\left(\lambda \frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left(\lambda \frac{\partial T}{\partial y}\right)+\frac{\partial }{\partial z}\left(\lambda \frac{\partial T}{\partial z}\right)+Q=\rho c\frac{\partial T}{\partial t}$$

（1）

$$Q=qA$$

（2）

式中：T——物体瞬态温度，K；$\lambda$——材料导热系数，W/（m·K）；$Q$——物体表面换热量，J；$q$——物体单位面积换热量，J；$A$——传热面比表面积，m/g；$\rho$——材料密度，kg/m³；$c$——比热，J/（kg·K）；$t$——时间，s。

边界条件指物体与外界接触界面的特征，本文采用第三类边界条件，即热对流换热边界，此时物体周围介质温度和物体表面的换热规律为：

$$h\left(T_{\mathrm{f}}-T_{\mathrm{w}}\right)=-\lambda {\left(\partial T/\partial n\right)}_{\mathrm{w}}=q$$

（3）

式中：$T_{\mathrm{w}}$——物体表面温度，K；（$\partial T$/$\partial n$）_w——物体温度梯度；$T_{\mathrm{f}}$——周围介质温度，K；$h$——介质与材料对流换热系数，W/（m²·K）。

由上式可得温度场：

$$T=T(x,y,z,t)$$

（4）

热弹性理论平衡方程^［21］，与弹性理论平衡方程类似，但应力分量包含温度应力，具体公式如下：

$$\left\{\begin{array}{c}\frac{\partial {\sigma }_x}{\partial x}+\frac{\partial {\sigma }_{yx}}{\partial y}+\frac{\partial {\sigma }_{zx}}{\partial z}+X=\mathrm{0}\\ \frac{\partial {\tau }_{xy}}{\partial x}+\frac{\partial {\sigma }_y}{\partial y}+\frac{\partial {\sigma }_{zy}}{\partial z}+Y=\mathrm{0}\\ \frac{\partial {\tau }_{xz}}{\partial x}+\frac{\partial {\sigma }_{yz}}{\partial y}+\frac{\partial {\sigma }_z}{\partial z}+Z=\mathrm{0}\end{array}\right.$$

（5）

式中：${\sigma }_x$、${\sigma }_y$、${\sigma }_z$——x、y、z方向上的正应力；${\tau }_{xy}$、${\sigma }_{yx}$、${\tau }_{xz}$——剪应力；X、Y、Z——分别为x、y、z方向的体积应力。

进一步推导可得热弹性理论本构方程为：

$$\left\{\begin{array}{l}{\sigma }_x=\mathrm{\lambda }\text{'}e+\mathrm{2}G{\epsilon }_x-\beta ∆T\\ {\sigma }_y=\mathrm{\lambda }\text{'}e+\mathrm{2}G{\epsilon }_y-\beta ∆T\\ {\sigma }_z=\mathrm{\lambda }\text{'}e+\mathrm{2}G{\epsilon }_z-\beta ∆T\\ {\tau }_{xy}=G{\gamma }_{xy}\\ {\tau }_{xz}=G{\gamma }_{xz}\\ {\tau }_{yz}=G{\gamma }_{yz}\end{array}\right.$$

（6）

其中：$$\lambda \text{'}=$$$$E\mu /[(\mathrm{1}+\mu )(\mathrm{1}-\mathrm{2}\mu )]$$

（7）

$$G=E/[\mathrm{2}(\mathrm{1}+\mu )]$$

（8）

$$\beta =E\alpha /(\mathrm{1}-\mathrm{2}\mu )=(\mathrm{3}\lambda +\mathrm{2}G)\alpha$$

（9）

式中：$∆T$——温差；$\lambda \text{'}$——拉梅系数；$G$——剪切模量；$\beta$——热应力系数；$\mu$——泊松比；$E$——弹性模量。

RFPA数值计算中，当单元应力状态满足破坏准则时，单元发生损伤并累积，此时单元的弹性模量损伤累积为^［22］：

$$E=(\mathrm{1}-D)E_{\mathrm{0}}$$

（10）

式中：E——损伤后的弹性模量；E₀——初始弹性模量；D——损伤变量。

（1）拉伸状态下的损伤，以拉应力为判断依据，此时岩石单元的损伤变量为：

$$D=\left\{\begin{array}{ll}\mathrm{0},& (\epsilon <{\epsilon }_{\mathrm{0}})\\ \mathrm{1}-{\sigma }_r/(\epsilon E_{\mathrm{0}}),& ({\epsilon }_{\mathrm{0}}<\epsilon <{\epsilon }_{\mathrm{u}})\\ \mathrm{1},& (\epsilon >{\epsilon }_{\mathrm{u}})\end{array}\right.$$

（11）

式中：${\sigma }_{\mathrm{r}}$——残余强度；${\epsilon }_{\mathrm{0}}$——最大弹性拉应变；${\epsilon }_{\mathrm{u}}$——极限拉应变。

因而拉伸状态下的本构关系为：

$${\sigma }_{ij}=\left\{\begin{array}{ll}\lambda \text{'}{\epsilon }_k{\delta }_{ij}+\mathrm{2}G{\epsilon }_{ij}-\beta ∆T{\delta }_{ij},& (\epsilon <{\epsilon }_{\mathrm{0}})\\ {\sigma }_{\mathrm{r}}/(\epsilon E_{\mathrm{0}})(\lambda \text{'}{\epsilon }_k{\delta }_{ij}+\mathrm{2}G{\epsilon }_{ij}-\beta ∆T{\delta }_{ij}),& ({\epsilon }_{\mathrm{0}}<\epsilon <{\epsilon }_{\mathrm{u}})\\ \mathrm{0},& (\epsilon >{\epsilon }_{\mathrm{u}})\end{array}\right.$$

（12）

（2）压缩状态下的损伤，以摩尔-库伦准则判断：

$${\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}(\mathrm{1}+\mathrm{s}\mathrm{i}\mathrm{n}\phi )/(\mathrm{1}-\mathrm{s}\mathrm{i}\mathrm{n}\phi )\ge {\sigma }_{\mathrm{c}\mathrm{o}}$$

（13）

式中：${\sigma }_{\mathrm{1}}$——最大主应力；${\sigma }_{\mathrm{c}\mathrm{o}}$——岩石抗压强度；$\phi$——岩石内摩擦角。

此时压缩状态下的本构关系为：

$${\sigma }_{ij}=\left\{\begin{array}{ll}\lambda \text{'}{\epsilon }_k{\delta }_{ij}+\mathrm{2}G{\epsilon }_{ij}-\beta ∆T{\delta }_{ij},& (\epsilon <{\epsilon }_{\mathrm{c}})\\ \frac{{\sigma }_{\mathrm{r}}}{\epsilon E_{\mathrm{0}}}(\lambda \text{'}{\epsilon }_k{\delta }_{ij}+\mathrm{2}G{\epsilon }_{ij}-\beta ∆T{\delta }_{ij}),& ({\epsilon }_{\mathrm{0}}<\epsilon <{\epsilon }_{\mathrm{u}})\end{array}\right.$$

（14）

式中：${\sigma }_{\mathrm{r}}$——残余强度；${\epsilon }_{\mathrm{c}}$——最大弹性压应变。

岩石材料的弹性模量、密度、泊松比等参数具有较大离散性，这也是岩石类材料的主要特征之一。这种离散型可以利用分布函数来表达，RFPA使用Weibull分布来考虑材料非均质性，将宏观的线性破坏问题转变为细观线性破坏问题，RFPA中的非均质性表达式为^［23］：

$$\phi \left(\alpha \right)=\frac{m}{{\alpha }_{\mathrm{0}}}{(\frac{\alpha }{{\alpha }_{\mathrm{0}}})}^{m-\mathrm{1}}{e}^{-{(\frac{\alpha }{{\alpha }_{\mathrm{0}}})}^m}$$

（15）

式中：$\phi \left(\alpha \right)$——非均质函数；$\alpha$——细观单元的材料参数；${\alpha }_{\mathrm{0}}$——材料参数平均值；$m$——RFPA数值模型中定义的材料均质度系数，代表材料性质的均匀性，m越高表示材料越均质，离散程度越低。

2 数值模型

为了更好的分析井筒高温花岗岩受冷冲击作用后裂纹的扩展情况，因此利用RFPA软件建立模型进行模拟，如图1所示。该模型为边长300 mm的正方形，内部正中央设有半径50 mm的圆形井筒；模型的初始温度为300 ℃，冷冲击流体的温度为20 ℃，网格划分为500×500；此外该模型考虑围压的影响，参考前人的研究^［24］，将模型四周围压设置一致，同时受限于客观实验因素，因此围压设为3 MPa；各边界设置为绝热边界，且冷冲击边界采用第三类热交换边界；模型的其他参数条件依据花岗岩力学实验的结果^［25］及相应文献研究进行赋值，对应参数如表1所示。

图1  井筒模型

Fig.1  Wellbore model

表1  模型参数

Table 1  Model parameters

参数	取值
弹性模量/MPa	18130
抗压强度/MPa	118
抗拉强度/MPa	12
密度/（kg·m-3）	2600
泊松比	0.25
均质度	50
拉压比	1/10
内摩擦角/(°)	30
导热系数/[W·(m·K)-1]	1$\times$10-5
热膨胀系数/ºC-1	1$\times$106
热容量/[(J·(m3·K)-1]	1.5
对流换热系数/[W·(m2·K)-1]	1000
低温流体介质/ºC	20
${\sigma }_x$、${\sigma }_y$/MPa	3、3

由弹性力学知，地应力作用时的井筒周围弹性体应力为^［26］：

$${\sigma }_{\mathrm{r}\mathrm{s}}=\frac{\left(K+\mathrm{1}\right)}{\mathrm{2}}\sigma \left(\mathrm{1}-\frac{a^{\mathrm{2}}}{r^{\mathrm{2}}}\right)+\frac{\left(K-\mathrm{1}\right)}{\mathrm{2}}\sigma \left(\mathrm{1}-\mathrm{4}\frac{a^{\mathrm{2}}}{r^{\mathrm{2}}}+\mathrm{3}\frac{a^{\mathrm{4}}}{r^{\mathrm{4}}}\right)\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\theta )$$

（16）

$${\sigma }_{\mathrm{\theta }\mathrm{s}}=\frac{\left(K+\mathrm{1}\right)}{\mathrm{2}}\sigma \left(\mathrm{1}+\frac{a^{\mathrm{2}}}{r^{\mathrm{2}}}\right)-\frac{\left(K-\mathrm{1}\right)}{\mathrm{2}}\sigma \left(\mathrm{1}+\mathrm{3}\frac{a^{\mathrm{4}}}{r^{\mathrm{4}}}\right)\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\theta )$$

（17）

式中：${\sigma }_{\mathrm{r}\mathrm{s}}$、${\sigma }_{\mathrm{\theta }\mathrm{s}}$——分别为径向和周向应力，MPa；$K$——侧压系数；$\sigma$——竖向地应力，MPa，$a$——半径，m；$r$、$\theta$——坐标系坐标。

对于半无限大固体，第三类边界条件的瞬态温度场分布为^［27］：

$$T\left(r,t\right)=T_{\mathrm{i}}+(T_{\mathrm{0}}-T_{\mathrm{i}})\left[erfc\left(\frac{r}{\mathrm{2}\sqrt[]{\alpha t}}\right)-\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{hr}{\lambda }+\frac{h^{\mathrm{2}}\alpha t}{{\lambda }^{\mathrm{2}}}\right)erfc(\frac{r}{\mathrm{2}\sqrt[]{\alpha t}})\right]$$

（18）

$$erfc(r)$$=1-$$erf(r)$$

（19）

$$erf(r)$$=$$\frac{\mathrm{2}}{\sqrt[]{\mathrm{\pi }}}{\int }_a^r{e}^{-y^{\mathrm{2}}}\mathrm{d}y$$

（20）

式中：r——计算处的半径，m；a——圆孔模型内径，m；$\lambda$——材料热传导系数，W/（m·K）；$erfc\left(*\right)$——高斯互补误差函数；$h$——介质与材料对流换热系数，W/（m²·K）；$T_{\mathrm{i}}$——初始温度，K；$T_{\mathrm{0}}$——低温介质温度，K。

圆孔周围任意一点的应力计算公式为^［28］：

$${\sigma }_{\mathrm{r}\mathrm{t}}=\frac{{\alpha }_{\mathrm{t}}E}{r^{\mathrm{2}}}\left[(\frac{r^{\mathrm{2}}-a^{\mathrm{2}})}{b^{\mathrm{2}}-a^{\mathrm{2}}}){\int }_a^{b}T\left(r,t\right)rdr-{\int }_a^{r}T\left(r,t\right)r\mathrm{d}r\right]$$

（21）

$${\sigma }_{\mathrm{\theta }\mathrm{t}}=\frac{{\alpha }_{\mathrm{t}}E}{r^{\mathrm{2}}}\left[(\frac{r^{\mathrm{2}}+a^{\mathrm{2}})}{b^{\mathrm{2}}-a^{\mathrm{2}}}){\int }_a^{b}T\left(r,t\right)rdr+{\int }_a^{r}T\left(r,t\right)r\mathrm{d}r-T\left(r,t\right)r^{\mathrm{2}}\right]$$

（22）

式中：${\sigma }_{\mathrm{r}\mathrm{t}}$、${\sigma }_{\mathrm{\theta }\mathrm{t}}$——分别径向和轴向热应力，MPa；$E$——弹性模量，MPa；${\alpha }_{\mathrm{t}}$——热膨胀系数；r——计算处的半径，m；a、b——分别为圆孔模型内、外半径，m；$T\left(r,t\right)$——温度场函数。

此时两种情况叠加，可以得到冷冲击作用下井筒周围岩石应力状态分布：

$${\sigma }_{\mathrm{r}}={\sigma }_{\mathrm{r}\mathrm{s}}+{\sigma }_{\mathrm{r}\mathrm{t}}$$

（23）

$${\sigma }_{\mathrm{\theta }}={\sigma }_{\mathrm{\theta }\mathrm{s}}+{\sigma }_{\mathrm{\theta }\mathrm{t}}$$

（24）

利用以上公式结合本构方程可以得到温度-应力耦合作用下各时刻温度场和应力场分布。本模型设置围压固定，侧压系数为1，考虑径向温度和应力分布，考虑径向温度场应力场，$\theta$取0，此时温度和应力分布如图2。

图2  不同时刻径向温度、应力分布解析解与数值解对比

Fig.2  Comparison of analytical and numerical

solutions of radial temperature and stress

；

distribution at different moments

由图2可以看出随着冷冲击时间的增加，井筒周围径向温度逐渐降低，应力场先增加后降低。温度场分布和应力场分布的数值解与解析解存在一定偏差，但总体上趋势相近，误差在合理范围内，因此利用RFPA研究井筒冷冲击过程的数值计算方法可靠。

3 结果与讨论

3.1 冷冲击作用下井筒应力分析

为了探究井筒岩石仅在冷冲击作用下的应力变化趋势，因此不考虑岩石的裂纹扩展情况，将模型的抗压强度设置为1000 MPa，其余条件不变。

图3是冷冲击过程中不同时刻模型最小主应力曲线，首先基于RFPA软件的特性，所有力学分析部分，规定压应力为正值，拉应力为负值，从图3中的最小主应力曲线可以明显观察到岩石分为两个部分，靠近模型内部部分受到拉应力影响，而靠近边界部分则受到压应力影响。其次，根据曲线的变化规律，可以分为前期和后期。前期：如图3（a）所示，随着降温时间的增长，模型冷冲击边界的拉应力不断增加，对应图中最小主应力大小朝负值方向增加，此时拉应力从10 s时的22.04 MPa增大到40 s时的31.01 MPa，并且拉应力也从冷冲击边界向模型内部扩张，此时的拉应力已经大于花岗岩的抗拉强度12 MPa，这也可以解释井筒花岗岩冷冲击边界起裂现象。随着r不断增大，最小主应力从负值逐渐变为正值，即拉应力不断减小，当r大于一定值时，拉应力几乎减小为0，甚至随着r增大，而产生了压应力，此时模型大部分区域没有发生温度变化。后期如图3（b）所示，随着冷冲击时间的增长，模型降温区域继续扩大，但靠近边界区域拉应力始终保持较高水平（远大于抗拉强度）。100 s时拉应力35.63 MPa，在t=300、500、1000 s，拉应力依次降低为30.42、28.51、26.14 MPa。造成冷冲击过程中前期、后期的应力变化先增后减的原因是温度梯度的变化。前期冷冲击边界温度降低快，且降温范围小，因此温度梯度不断增大，此时拉应力不断增大；后期降温范围扩大，而冷冲击边界温度与低温介质相近，温度降低幅度小，降温区域不断扩大，此时温度梯度小，所以拉应力不断降低。整体上看，拉应力呈现先增后减的趋势，且模型拉应力分布为靠近冷冲击边界处拉应力最大，离冷冲击边界越远，拉应力越小。

图3  冷冲击过程中不同时刻模型最小主应力

Fig.3  Minimum principal stress curve at each

moment under cold shock

依据对冷冲击过程中应力场的分析可以总结出：冷冲击作用时，岩石表面拉应力随着冷冲击时间的增加，先升高至峰值后缓慢下降。300 ℃花岗岩在接受20 ℃低温流体冷冲击过程中产生的拉应力大小普遍大于一般花岗岩的抗拉强度，这可以解释冷冲击作用对花岗岩具有致裂效果，并且上述数值计算的模型应力分布是假设材料均质度高的基础上进行的，数据离散程度低，而实际上花岗岩内部颗粒分布不均，并且含有较多原生裂纹和孔隙，在冷冲击时会造成裂隙尖端产生应力集中效应，更加有利于裂纹的萌生和扩展。

3.2 冷冲击作用下井筒裂纹扩展过程

图4为模型声发射（Acoustic Emission， AE）事件数量随时间变化图，由图可见，岩石在经受冷冲击作用后裂纹扩展过程可以大致分为3个时期；前期（0~20 s）、中期（20~350 s）、后期（350~2000 s）。

图4  AE数量随时间变化曲线

Fig.4  Curve of AE number over time

在冷冲击前期，由冷冲击引起的热应力是在为裂纹的产生积蓄应变能，当储存的应变能达到能够使岩石单元间的接触断开时，就会产生微裂纹^［29-30］，此时声发射事件数量没有明显的增加。在冷冲击中期随着热交换的持续进行，冷冲击进程继续推进，降温区域继续向模型内部扩展，同时裂纹的产生伴随应力释放和消耗，但随着冷冲击时间的不断增大，应变能不断增大，能够使更多的岩石单元间接触断开，伴随着AE数量的急剧增大，裂纹开始大规模产生。冷冲击后期，冷冲击区域向模型内部扩展，此时模型整体温度与冷冲击流体温度趋于一致，因此热交换效率低，温度梯度降低，降温区域逐渐覆盖整个模型，此时裂纹缓慢延伸直至停止，因此声发射事件数量发展缓慢。

图5、图6分别为各时刻井筒裂纹扩展过程的温度场和应力场。可以看出300 ℃井筒模型在遭受20 ℃流体冷冲击作用后裂纹扩展可以分为3个阶段。前期、中期以及后期。冷冲击前期，井筒周围温度迅速降低，形成温度梯度，产生很大的拉应力，此时井筒周围产生环形拉应力区，在井筒周围开始出现单元破坏，随着时间的增加，降温范围逐渐向井筒外围扩大，单元破坏情况向外延伸，且大体均匀分布在井筒周围。冷冲击中期，随着井筒边界冷冲击时间的增加，降温区域逐渐向井筒外围扩展，环形拉应力区也持续向外扩展，裂纹随着拉应力区沿径向朝外扩展。冷冲击后期，降温范围持续扩展直至整个模型温度发生变化，拉应力区继续向外部扩展，但此时的拉应力大小逐渐降低直至小于模型抗拉强度，裂纹基本停止扩展，形成井筒冷冲击裂纹扩展的最终形态。由此可以看出，在冷冲击作用下，裂纹扩展主要是由拉应力决定的，拉应力是由于降温区域处产生的温度梯度引起的，而井筒冷冲击产生的裂纹大体上均匀分布在井筒周围，随着环形拉应力区向外移动而扩展。

图5  各时刻井筒裂纹扩展过程的温度场

Fig.5  Temperature field of the wellbore fracture extension process at different moments

图6  各时刻井筒裂纹扩展过程的应力场

Fig.6  Stress field of the wellbore fracture extension process at different moments

3.3 冷冲击作用时裂纹扩展的影响因素

针对于干热岩储层，不同地层深度对应的温度和围压不同，在相同的温度和围压下，不同干热岩工程中的井径不同，本节考虑围压、井径、温度等因素，对井筒冷冲击过程中的裂纹扩展进行研究，分析各因素对冷冲击作用时裂纹扩展的影响。

3.3.1 围压

为讨论围压对井筒冷冲击效果的影响，在所述模型的基础上将围压依次设置为1、3、5、10、20 MPa，其它参数条件保持不变。

图7为不同围压下井筒冷冲击裂纹扩展的最终形态，图8为不同围压下裂纹扩展半径和数量。可以看出，随着围压由1 MPa依次增大到3、5、10、20 MPa，井筒模型在冷冲击过程中产生的裂纹扩展半径由80.6 mm依次降低为58.2、31.3、15.2、5.7 mm，裂纹数量也由44条依次减小到29、23、17、6条，裂纹数量和扩展半径几乎呈线性减小，此时的裂纹扩展区域面积呈指数减小。与此同时，冷冲击产生的裂隙开度随着围压的增大而减小。另外，随着围压的增大，冷冲击产生的裂纹出现分支更少。因此可以得出结论，井筒在经受冷冲击作用时，井筒水平围压越大，冷冲击作用产生的裂纹数量、裂纹扩展半径、裂隙开度越小，冷冲击致裂效果随着水平围压的增大显著减小。

图7  不同围压裂纹扩展最终形态

Fig.7  Final morphology of fracture extension at

different confining pressure

图8  不同围压裂纹扩展半径和数量

Fig.8  Diffusion radius and number of fracture at different confining pressure

3.3.2 井径

为讨论井径对井筒冷冲击效果的影响，依旧在所述模型的基础上，围压设置为3 MPa，井筒半径依次设置为30、50、70 mm，其它参数不变。

图9为不同井径裂纹扩展最终形态，图10为不同井径裂纹有效扩展半径和数量，由图可以看出，随着井筒井径的增大，诞生于井筒周围的裂纹数量明显增多，裂纹扩展半径没有显著变化，但是总体来看，随着井径由30 mm依次增大到50、70 mm时，冷冲击下的井筒模型裂纹扩展面积明显增大，由174.2 cm²依次增加到289.2、332.1 cm²。由以上分析可以看出，井筒在经受冷冲击作用时，井筒井径越大，井筒周围产生的裂纹数量越多，裂纹扩展半径无显著变化，但裂纹扩展面积增大。

图9  不同井径裂纹扩展最终形态

Fig.9  Final morphology of fracture extension

for different well diameters

图10  不同井径裂纹扩展半径和数量

Fig.10  Radius and number of fracture extention in different well diameters

3.3.3 温度

为探究温度对井筒模型冷冲击效果的影响，在所述模型的基础上，改变模型温度，依次设置为200、300、400、500 ℃，其它参数保持不变。

图11为不同模型温度裂纹扩展最终形态，图12为不同模型温度裂纹扩展半径和数量。由图可知，随着模型温度由200 ℃依次增大为300、400、500 ℃时，裂纹扩展半径明显增大，冷冲击产生的裂纹数量显著增多，其中裂纹扩展半径由200 ℃时的20 mm增大到500 ℃时的81.2 mm，裂纹数量也由22条增加到39条。同时由图10可以看出，伴随着模型温度的升高，冷冲击产生的裂隙开度增大。由此可以得出，温度对干热岩井筒围岩冷冲击效果的影响显著，随着围岩温度的升高，冷冲击产生的裂纹数量、扩展半径、开度增大。

图11  不同模型温度裂纹扩展最终形态

Fig.11  Final morphology of fracture extension

at different model temperatures

图12  不同模型温度裂纹扩展半径和数量

Fig.12  Fracture extention radius and number at different model temperatures

4 结论

本文通过数值模拟方法，使用RFPA数值模拟软件构建井筒冷冲击模型，对冷冲击作用下井筒裂纹扩展进行分析，同时考虑围压、井径以及温度对冷冲击效果的影响，得出的结论如下：

（1）冷冲击作用时，井筒岩石表面出现拉应力，同时随着冷冲击时间的增加，拉应力先升至峰值后缓慢下降。在实际工况下，产生的拉应力大小普遍大于一般花岗岩的抗拉强度，会使井筒岩石产生裂纹。

（2）井筒模型在冷冲击过程中，裂纹扩展可以大致分为前、中、后期，冷冲击前期，井筒周围温度迅速降低，井壁产生环形拉应力区，此时井筒围岩出现单元破坏，随着冷冲击时间的增加，单元破坏情况向外延伸，诞生均匀分布在井筒周围的微小裂纹；冷冲击中期，随着时间的增加井筒围岩温度继续降低，环形拉应力区继续向井筒外围扩展，裂纹随着拉应力区沿径向朝外扩展；冷冲击后期，降温范围持续扩展直至整个模型温度发生变化，拉应力大小逐渐降低直至小于模型抗拉强度，裂纹扩展速度减缓直至停止扩展，形成裂纹扩展的最终形态。

（3）围压、井径、温度对井筒围岩冷冲击时的破坏损伤效果影响显著，其中，围压对井筒围岩冷冲击时的破坏损伤效果影响显著，围压对冷冲击裂纹扩展起到抑制作用，围压越大，裂纹扩展半径、数量、开度越低。井径对冷冲击裂纹扩展起到促进作用，随着井径的增大，裂纹数量、裂纹扩展面积增大，但对裂纹开度和裂纹扩展半径没有显著影响。温度对冷冲击裂纹扩展起到促进作用，随着井筒围岩温度增大，冷冲击产生的裂纹数量、扩展半径以及开度明显增大。