引言

在现代矿业工程中，通过加热的方式来辅助破碎岩石已经有了超过一个世纪的历史^[1-3]。坚硬地层中，使用热裂解方法进行钻井，其机械钻速能够达到传统的旋转钻井方式的5$\sim$10倍^[4-6]。另外，采用热力射流的方法进行钻井作业，钻井工具并不会与岩石表面发生接触，所以在钻进时不会存在诸如传统的机械钻井过程中的钻具磨损和失效的问题^{[4, 7-8]}。热力射流的冲击作用之后，大量的热能由高温流体传递给岩石^[9-10]。由于岩石本身是由各种不同类型的矿物组成，而各种矿物组分的热导率存在很大差别，即相当于对岩石表面进行剧烈的局部加热行为，在此过程中，岩石表面会形成很大的温度梯度，其内部的局部受热部分会发生非均匀膨胀，由此产生压缩的温度应力，不仅会引起岩石内的原生裂缝不断扩展，而且导致岩石内部矿物颗粒之间出现新的裂纹，当其长度达到一个临界点，岩石表面受热部分发生屈曲从整块岩石上脱落^[11-12]。岩石在热冲击条件下的响应过程如图 1所示^[13]。

1 热裂解钻井模型建立及求解

目前，高温射流破岩主要基于连续管技术，其结构如图 2所示^[14]。因为热裂解破碎的时间很短，所以将由喷嘴喷出的超临界高温流体作为已知热源，将高温介质与岩石表面接触部分视为强制对流换热，忽略因为辐射等原因所造成的热量损失，将其简化为如图 3所示的模型。

高温射流喷嘴的半径及井筒半径分别设为5.0和25.4 mm^[15]，井筒以外区域为地层，在计算过程中可以作为无界区域进行处理。将模型中井筒部分与高温射流相接触的上表面设置为Robin边界，即第三类边界条件。其余边界都设置为绝热边界条件，这是因为岩石基质的导热率很低，在计算时间区间内，不会有热量扩散到岩石的外表面上^[16]。初始温度设置为岩石所处的地层温度。假定岩石物性不随温度变化，得到热传导控制方程

$ \rho {{C}_{\rm P}}(T)\dfrac{\partial T}{\partial t}=\dfrac{\partial }{\partial x}\left[ k(T)\dfrac{\partial T}{\partial x} \right]+\\{\kern 40pt}\dfrac{\partial }{\partial y}\left[ k(T)\dfrac{\partial T}{\partial y} \right]+\dfrac{\partial }{\partial z}\left[ k(T)\dfrac{\partial T}{\partial z} \right] $

(1)

式中：$\rho$—密度，kg/m$^{{\rm 3}}$；

$C_{\rm P}$—定压热容，J/K；

$T$—温度，℃；

$t$—时间，s；

$k$—导热系数，W/(m$\cdot$℃)；

$x$，$y$，$z$—$x$，$y$，$z$三个方向的坐标，m。

将井壁处理为绝热边界，设置为第二类边界条件，即Neumann边界条件，表示为

$ -k\left( T \right)\dfrac{\partial T}{\partial {{x}_{\rm i}}}{{\mathit{\boldsymbol{n}}}_{\rm i}}=0 $

(2)

式中：${{x}_{\rm i}}$—$x$向上的单位长度；

$\mathit{\boldsymbol{n}}_{\rm i}$—单位外法向向量。

在热力射流喷嘴的直径范围内，采用第三类边界条件，即Robin边界条件，设定高温流体与岩石表面的换热系数以及岩石所处的地层温度(可以使用初始温度表示)，表示如下

$ -k\left( T \right)\dfrac{\partial T}{\partial {{x}_{\rm i}}}{{\mathit{\boldsymbol{n}}}_{\rm i}}=h\left( {{T}_{\rm in}}-{{T}_{0}} \right) $

(3)

式中：$h$—高温流体与岩石表面的对流传热系数，$\times$10 kW/(m$^{{\rm 2}}\cdot$℃)^[9]；

$T_{\rm in}$—高温流体的入射温度，℃；

$T_0$—初始温度，℃。

采用Crank-Nicolson差分格式^[17]对式(1)进行离散，可得

$ {T^{{{(n + 1)}^1}}} \!-\! {T^n} \!=\! \dfrac{{\alpha ({T^{{{\left( {n + 1} \right)}^1}}})\Delta t}}{2}\dfrac{{{\partial ^2}{T^{{{(n + 1)}^1}}}}}{{\partial {x^2}}}\! +\! \dfrac{{\alpha ({T^n})\Delta t}}{2} \cdot \\[2pt]{\kern 40pt} \dfrac{{{\partial ^2}{T^n}}}{{\partial {x^2}}} + \dfrac{{\Delta t}}{2}\dfrac{1}{{\rho {C_{\rm{P}}}({T^{{{(n + 1)}^1}}})}}\dfrac{{\partial k({T^{n + 1}})}}{{\partial x}}\dfrac{{\partial {T^{{{(n + 1)}^1}}}}}{{\partial x}} + \\[2pt]{\kern 40pt} \dfrac{{\Delta t}}{2}\left[ {\dfrac{1}{{\rho {C_P}(T)}}\dfrac{{\partial k(T)}}{{\partial x}}} \right]\dfrac{{\partial T}}{{\partial x}} + \alpha ({T^n})\Delta t\dfrac{{{\partial ^2}{T^n}}}{{\partial {y^2}}} + \\[2pt]{\kern 40pt} \dfrac{{\Delta t}}{{\rho {C_P}({T^n})}}\dfrac{{\partial k({T^n})}}{{\partial y}}\dfrac{{\partial {T^n}}}{{\partial y}} + \alpha ({T^n})\Delta t\dfrac{{{\partial ^2}{T^n}}}{{\partial {z^2}}} + \\[2pt]{\kern 40pt} \dfrac{{\Delta t}}{{\rho {C_{\rm{P}}}({T^n})}}\dfrac{{\partial k({T^n})}}{{\partial y}}\dfrac{{\partial {T^n}}}{{\partial y}} $

(4)

$ {{T}^{{{(n+1)}^{2}}}}\!-\!{{T}^{{{(n+1)}^{1}}}}\!=\!\dfrac{\alpha ({{T}^{{{\left( n+1 \right)}^{2}}}})\Delta t}{2}\dfrac{{{\partial }^{2}}{{T}^{{{(n+1)}^{2}}}}}{\partial {{y}^{2}}}- \\[2pt]{\kern 40pt} \dfrac{\alpha ({{T}^{n}})\Delta t}{2}\dfrac{{{\partial }^{2}}{{T}^{n}}}{\partial {{y}^{2}}}+ \dfrac{\Delta t}{2}\dfrac{1}{\rho {{C}_{\rm P}}({{T}^{{{\left( n+1 \right)}^{2}}}})} \cdot \\[2pt]{\kern 40pt} \dfrac{\partial k({{T}^{{{\left( n+1 \right)}^{2}}}})}{\partial y}\dfrac{\partial {{T}^{{{\left( n+1 \right)}^{2}}}}}{\partial y}- \dfrac{\Delta t}{2}\dfrac{1}{\rho {{C}_{\rm P}}({{T}^{n}})}\cdot \\[2pt]{\kern 40pt} \dfrac{\partial k({{T}^{n}})}{\partial y}\dfrac{\partial {{T}^{n}}}{\partial y} $

(5)

$ {{T}^{n+1}}-{{T}^{{{(n+1)}^{2}}}}=\dfrac{\alpha ({{T}^{n+1}})\Delta t}{2}\dfrac{{{\partial }^{2}}{{T}^{n+1}}}{\partial {{z}^{2}}}- \dfrac{\alpha ({{T}^{n}})\Delta t}{2} \cdot \\[2pt]{\kern 40pt} \dfrac{{{\partial }^{2}}{{T}^{n}}}{\partial {{z}^{2}}}+ \dfrac{\Delta t}{2}\dfrac{1}{\rho {{C}_{\rm P}}({{T}^{n+1}})}\dfrac{\partial k({{T}^{n+1}})}{\partial z}\dfrac{\partial {{T}^{n+1}}}{\partial z}- \\[2pt]{\kern 40pt} \dfrac{\Delta t}{2}\dfrac{1}{\rho {{C}_{\rm P}}({{T}^{n}})}\dfrac{\partial k({{T}^{n}})}{\partial z}\dfrac{\partial {{T}^{n}}}{\partial z} $

(6)

式中：$\alpha$—线膨胀系数，无因次；

$\Delta t$—时间步长，s。

对于Robin边界条件，如图 4所示。取$T(m, n)$为控制节点，在控制体积内，入射流体传入的热量加上由周边节点通过热传导传出的热量等于零，通过能量守恒定律，可将其写成式(7)形式的差分格式。

$ h\Delta y\left( {{T}_{\rm in}}-{{T}_{m, n}} \right)+k\Delta y\dfrac{{{T}_{m-1, n}}-{{T}_{m, n}}}{\Delta x}+\\{\kern 40pt}\dfrac{k\Delta x}{2}\dfrac{{{T}_{m, n+1}}-{{T}_{m, n}}}{\Delta y}+ \\{\kern 40pt} \dfrac{k\Delta x}{2}\dfrac{{{T}_{m, n-1}}-{{T}_{m, n}}}{\Delta y}=0 $

(7)

式中：$\Delta x$，$\Delta y$—$x$，$y$向的步长，m；

${T}_{m, n}$—网格($m$，$n$)处的温度，℃。

设$\Delta x=\Delta y$，则可以将式(7)改写成

$ {{T}_{m, n}}=\dfrac{\dfrac{h\Delta x}{k}{{T}_{\rm in}}+\left( \dfrac{2{{T}_{m-1, n}}+{{T}_{m, n+1}}+{{T}_{m, n-1}}}{2} \right)}{\left( 2+\dfrac{h\Delta x}{k} \right)} $

(8)

将岩石种类分别设置为砂岩、页岩和花岗岩这3种在钻井过程中常会碰到的岩石种类，其热物理性质和力学参数见表 1所示^[18]。

通过计算得到的温度场，当物体的温度由$T_1$升高到$T_2$时，则其热应变可以由下式进行计算

$ \varepsilon = \alpha \left( T_2-T_1 \right)=\alpha \Delta T $

(9)

式中：$\varepsilon$—应变，无因次；

$T_1$，$T_2$—热应变前后的温度，℃；

$\Delta T$—温度变化量，℃。

将应变用应力和温度梯度表示，可以得到

$ \left\{ \begin{array}{l} {{\varepsilon }_{x}}=\dfrac{1}{E}\left[ {{\sigma }_{x}}-\mu \left( {{\sigma }_{y}}+{{\sigma }_{z}} \right) \right]+\alpha {{T}_{\upsilon }} \\[3pt] {{\varepsilon }_{y}}=\dfrac{1}{E}\left[ {{\sigma }_{y}}-\mu \left( {{\sigma }_{x}}+{{\sigma }_{z}} \right) \right]+\alpha T_{\upsilon }^{'} \\[3pt] {{\varepsilon }_{z}}=\dfrac{1}{E}\left[ {{\sigma }_{z}}-\mu \left( {{\sigma }_{y}}+{{\sigma }_{z}} \right) \right]+\alpha T_{\upsilon }^{''} \\[3pt] {{\gamma }_{xy}}=\dfrac{2(1+\mu )}{E}{{\tau }_{xy}} \\[3pt] {{\gamma }_{xz}}=\dfrac{2(1+\mu )}{E}{{\tau }_{xz}} \\[3pt] {{\gamma }_{yz}}=\dfrac{2(1+\mu )}{E}{{\tau }_{yz}} \end{array} \right. $

(10)

式中：$\varepsilon _x$，$\varepsilon _y$，$\varepsilon _z$—$x$，$y$，$z$向的应变，无因次；

$E$—杨氏模量，MPa；

${\sigma }_{x}$，${\sigma }_{y}$，${\sigma }_{z}$—$x$，$y$，$z$向的的正应力，MPa；

$\mu$—泊松比，无因次；

$T_ {\upsilon }$，$T_ {\upsilon }^{\prime}$，$T_ {\upsilon }^{''}$—$x$，$y$，$z$向的温度差，℃；

${\gamma }_{xy}$，${\gamma }_{xz}$，${\gamma }_{yz}$—$xy$，$xz$，$yz$面上的剪应力，MPa；

${\tau }_{xy}$，${\tau}_{xz}$，${\tau }_{yz}$—$xy$，$xz$，$yz$面上的剪应变，无因次。

边界条件采用位移边界条件，设其处于无限远地层当中，则可以写成

$ \left\{ \begin{array}{l} {{u}_{x}}=0 \\ {{u}_{y}}=0 \\ {{u}_{z}}=0 \\ \end{array} \right. $

(11)

式中：$u_{x}$，$u_{y}$，$u_{z}$—$x$，$y$，$z$向的位移，m；

对式(10)按位移进行求解，得到应力用应变及温差表示的表达式为

$ \left\{ \begin{array}{l} {{\sigma }_{x}}=\dfrac{E}{1+\mu }\left( \dfrac{\mu }{1-2\mu }e+{{\varepsilon }_{x}} \right)-\dfrac{E\alpha {{T}_{\upsilon }}}{1-2\mu } \\ {{\sigma }_{y}}=\dfrac{E}{1+\mu }\left( \dfrac{\mu }{1-2\mu }e+{{\varepsilon }_{y}} \right)-\dfrac{E\alpha T_{\upsilon }^{'}}{1-2\mu } \\ {{\sigma }_{z}}=\dfrac{E}{1+\mu }\left( \dfrac{\mu }{1-2\mu }e+{{\varepsilon }_{z}} \right)-\dfrac{E\alpha T_{\upsilon }^{''}}{1-2\mu } \\ {{\tau }_{xy}}=\dfrac{E}{2\left( 1+\mu \right)}{{\gamma }_{xy}} \\ {{\tau }_{yz}}=\dfrac{E}{2\left( 1+\mu \right)}{{\gamma }_{yz}} \\ {{\tau }_{xz}}=\dfrac{E}{2\left( 1+\mu \right)}{{\gamma }_{xz}} \\ \end{array} \right. $

(12)

式中：$e$—体积应变，无因次。

2 温度场分析

求解以上公式，得到井眼中心温度随时间变化情况如图 5所示。

从图 5可见，当射流介质刚接触到岩石表面时，因为周围环境与岩石表面存在巨大的温度梯度，在0$\sim$0.1 s内，井眼中心的温度急剧增大，超过了300 ℃，曲线斜率几乎接近90°；在1 s内，岩石表面的最高温度就会达到高温流体温度的90%。在10 $\sim$20 s时间内，曲线斜率逐渐趋于平缓，温度只增加了30 ℃左右。这表明在岩石受热过程中，井眼中心温度随着射流时间的增加而不断增大，但温度增大的速率却不断减小。这是由于岩石基质的导热性能较差，随着入射热量的不断积累，降低了入射的高温流体与岩石表面之间的热交换效率，从而导致岩石表面温度增长的速率越来越慢。

取$xy$平面，分别计算热裂解钻井过程中径向和轴向的温度分布，见图 6和图 7。

随着岩石表面受热时间的不断增加，喷嘴直径区域内的温度迅速增大，同时温度开始由井眼中心向四周传递。但是由于岩石本身的导热性差，因此温度传播的速率很低。由图 6可见，在径向上，井眼中心温度最高，距离井眼中心的距离越远，温度逐渐降低；而随着受热时间的不断增加，沿径向上的温度梯度也逐渐增大，即温度波及范围有限，由此能够在岩石内部形成较大的温度梯度和温度应力，造成岩石裂解破碎。

由图 7可见，在轴向上，离井眼表面的距离越远，温度越低。随着岩石表面受热时间增加，温度传播的距离也越远，且温度在轴向上能够产生比径向上更大的温度梯度。

岩石类型是高温热裂解钻井效率的最大影响因素^[19]。不同种类的岩石井底温度在1 s内的变化曲线见图 8。

由图 8可见，岩石种类对温度场的影响较大。岩石的比热越大，其吸收热量的能力也越强，在相同受热时间内，表面的温度也越高；岩石的热传导系数越大，其传播热量的能力也越强，因此，温度的波及范围越大，而岩石的导热系数越小，温度传播速度越慢，曲线斜率越大，能形成更大的温度梯度。总的说来，图中所示3种岩石的热物理性质相对比较接近，但温度场仍具有很大差异。因此，岩石种类对热裂解破岩效果的影响不能忽略。

3 温度应力分析

在热裂解钻井过程中，岩石的表层受到热冲击的影响会发生屈曲，最终从岩石表面剥离。%在柱坐标条件下，对岩石的受力情况进行研究。图 9为不同时间热裂解钻井沿井眼半径径向应力分布情况。

热裂解钻井过程中，岩石基质在径向上受到压应力的作用，且井眼中心处的压应力最大。这主要是因为岩石受热后由于其导热性差而造成温度梯度大，高温区域体积迅速膨胀，受到周围岩石的挤压所造成的。温度越高的区域，体积膨胀程度越大，受到的挤压也越严重，因此压应力越大，离受热区域越远，压应力迅速下降。而随着受热时间的增加，井眼中心处在径向上受到的压应力迅速增加；在1 s内压应力即可以达到最大值142.5 MPa。此后，压应力的峰值几乎不再增加，仅仅是随着热量的传递，最大压应力的区域扩大。由于模型具有对称性，沿井眼半径的环向应力随时间变化与径向应力变化具有相同的特征(图 10)。

图 11为不同时间热裂解钻井井眼平面内剪切应力沿井眼半径分布曲线。由图 11可见，热裂解钻井过程中，在井眼平面内，岩石基质受到的剪切力沿井眼半径呈对称分布，但剪切力的绝对值很小，不到0.4 MPa。这是因为在岩石受热过程中，假定岩石基质的物性是均匀且各向同性的，则岩石在沿径向的各个方向上其体积应变都是相等的，所以在井眼平面上几乎不受到剪应力的作用。

图 12为不同时间热裂解钻井井眼半径轴向剪切应力分布曲线。

由图 12可以明显看出，岩石在轴向上受到很大的剪切作用。这是因为在岩石受热过程中，中心处高温区域的岩石基质受热膨胀后受到挤压会发生屈曲，产生如图 13所示的沿$z$轴，即垂直于井眼平面的应变。将应力状态与材料强度进行比较，将花岗岩的平均抗剪强度9 MPa作为其强度极限^[13]，可以看出在0.1 s内岩石就会发生破坏。

由图 13可以看出，沿径向上，几乎只会在受热区域表面及其附近发生屈曲，远离受热区域则几乎没有位移发生。数值模拟得到的结果与文献中岩石的受力情况和屈曲响应具有很好的一致性^[20]。

4 结论

(1) 热裂解钻井过程中，由于岩石基质的导热性差，形成非均匀温度场，会在径向和轴向上产生温度梯度，从而形成温度应力。

(2) 热裂解钻井过程中，岩石基质在径向上受到压应力的作用(不考虑围压)，最大值达到142.5 MPa，且井眼中心处的压应力最大。

(3) 岩石中心受热部分会在高温作用下发生屈曲，产生垂直于井眼平面方向的应变，并且该方向上受到接近20 MPa的剪应力(不考虑围压)，超过花岗岩的抗剪强度，造成岩石表面在高温作用下从岩石基质剥离脱落，这与文献中的实验结果相符。