0 引言
PF的相变过程表明，温度是自支撑相变压裂成功实施的关键[5]。在相变过程中，一方面释放热量促使温度升高，另一方面，温度较低的PFFS对地层、裂缝进行降温，改变地层、裂缝温度场。因此，有必要对相变热作用下的裂缝温度场进行预测，为PFFS配方和施工参数的优化提供依据。关于压裂过程中的裂缝温度场，学者们展开了大量研究，李木坤等[6]在岩体热传导方程基础上，建立了热-流-固耦合模型，研究认为CO2在裂缝中流动发生温降而产生的温度应力有助于主裂缝周围复杂缝网的形成。Seth等[7]研发了计算水力压裂和关井期间裂缝温度的数值模型，研究了裂缝长度与岩石物理性质对温度的影响。Guo等[8]考虑温度、压力和CO2体积功的影响，建立了适用于碳酸盐岩储层酸压的温度分布计算模型，计算表明酸-岩相变热可导致裂缝内的温度升高约15 ℃。
PF相变涉及的化学反应是其内部各组分之间的有机化学反应，现有关于酸-岩相变热的研究成果无法直接应用于自支撑相变压裂技术，有必要建立PF相变反应动力学模型，明确相变热与温度、时间的关系，阐明相变热对裂缝温度场的影响规律。为此，本文首先利用差示扫描量热仪（DSC）的测试结果建立反应动力学模型，在此基础上，根据能量守恒定律建立压裂过程中地层和裂缝温度分布的数学模型，提出求解方法，并验证模型的可靠性，最后通过数值实验探讨了相变热、PF体积分数和比热容对裂缝温度的影响规律。
1 反应动力学模型和温度场模型
利用DSC实验测得的热通量曲线，建立相变反应动力学模型，以此建立放热速率、温度和时间之间的数学关系。在此基础上，构建考虑化学相变热的地层和裂缝温度场数学模型。
1.1 PF的反应动力学模型
采用DSC Q20测试PF在相变过程中的热通量。实验在高纯度氮气环境下进行，流量为50 mL/min。实验步骤为：①称取PF样品约10 mg放入铝坩埚中，快速密封后放入DSC样品室；②在DSC样品室放置一个密封的空铝坩埚作为参考；③以特定的升温速率（5，10，20 K/min）将测试样品从313 K（40 ℃）加热至473 K（200 ℃），并记录不同时间的热通量。
热通量为正值，其随温度的变化曲线接近于正态分布曲线（见图1 ），每条热通量曲线只有一个比较对称的放热峰，表明PF在相变过程中放出热量，且化学反应过程是一个单一的主化学反应，不存在其他副化学反应。从热通量曲线中可提取出各条曲线显示的开始相变温度、相变峰值温度以及相变热（见表1 ）。
$ \alpha(t)=1-\left[1-1.8 \times 10^{8} \exp \left(-\frac{8530}{T}\right) t\right]^{12.09}$
根据（1）式可以计算出不同温度和时间下PF的转化率，进而计算任意时间放出的热量。
1.2 裂缝温度场模型
如图2 所示，取裂缝中一个长、宽、高分别为dx、dy、dz的微元体进行分析，其热量传递包括4个部分：①裂缝内的流体不断向前运动，发生热对流；②缝内存在温度梯度，相邻微元体之间发生热传导；③缝内压裂液与岩体之间形成热交换；④PF在相变过程中释放一定的热量。因次，对于裂缝内的流体，需要考虑热对流、热传导、热交换、相变放热4个方面的影响。
包含热对流、热传导、热交换、相变热等因素的缝内温度场模型为：
$ \begin{array}{l} \lambda_{\mathrm{w}}\left(\frac{\partial^{2} T_{\mathrm{w}}}{\partial x^{2}}+\frac{\partial^{2} T_{\mathrm{w}}}{\partial y^{2}}+\frac{\partial^{2} T_{\mathrm{w}}}{\partial z^{2}}\right)- \\ c_{\mathrm{w}} \rho_{\mathrm{w}} \frac{K_{\mathrm{f}}}{\mu}\left(\frac{\partial p}{\partial x} \frac{\partial T_{\mathrm{w}}}{\partial x}+\frac{\partial p}{\partial y} \frac{\partial T_{\mathrm{w}}}{\partial y}+\frac{\partial p}{\partial z} \frac{\partial T_{\mathrm{w}}}{\partial z}\right)+ \\ \frac{\lambda_{\mathrm{r}}}{w / 2}\left(T_{\mathrm{r}}-T_{\mathrm{w}}\right)+m_{\mathrm{pf}} \Delta H_{\mathrm{R}} \frac{\mathrm{~d} \alpha}{\mathrm{~d} t} \frac{1}{\mathrm{~d} x \mathrm{~d} y \mathrm{~d} z}=c_{\mathrm{w}} \rho_{\mathrm{w}} \frac{\partial T_{\mathrm{w}}}{\partial t} \end{array}$
1.3 岩体温度场模型
对于岩体温度场，建立如图3 所示的地层岩体温度场物理模型，取岩体中的一个长、宽、高分别为dx、dy、dz的微元体进行分析，该微元体的热量受相邻微元体之间的热传导以及基质内流体的对流传热影响。
建立地层岩体温度场模型需满足的前提条件包括：①假设岩体和压裂液的热力学参数均质且恒定；②忽略摩擦和动能对换热的影响；③基质岩体和孔隙流体满足局部热平衡条件；④地层的热力学参数取岩
体骨架和孔隙流体的加权平均值。满足上述假设条件时，地层岩体的温度场模型为：
$ \begin{aligned} \lambda_{x} \frac{\partial}{\partial x}\left(\frac{\partial T_{\mathrm{r}}}{\partial x}\right)+ & \lambda_{y} \frac{\partial}{\partial y}\left(\frac{\partial T_{\mathrm{r}}}{\partial y}\right)+\lambda_{z} \frac{\partial}{\partial z}\left(\frac{\partial T_{\mathrm{r}}}{\partial z}\right)- \\ & c_{\mathrm{mav}} \rho_{\mathrm{mw}} \phi\left(v_{x} \frac{\partial T_{\mathrm{r}}}{\partial x}+v_{y} \frac{\partial T_{\mathrm{r}}}{\partial y}+v_{z} \frac{\partial T_{\mathrm{r}}}{\partial z}\right) \\ & m_{\mathrm{pm}} \Delta H_{\mathrm{R}} \frac{\mathrm{~d} \alpha}{\mathrm{~d} t} \frac{1}{\mathrm{~d} x \mathrm{~d} y \mathrm{~d} z}=c_{\mathrm{r}} \rho_{\mathrm{r}} \frac{\partial T_{\mathrm{r}}}{\partial t} \end{aligned}$
2 模型求解与验证
2.1 求解方法
初始时刻，整个物理模型中的温度等于原始储层温度，即：
$ T_{\mathrm{r}}(x, y, z, t)=T_{\mathrm{res}}(t=0)$
边界条件为：
$ \left\{\begin{array}{ll} q(x, t)=q_{\text {in }} & \forall t, x=0 \\ T_{\mathrm{r}}(x, t)=T_{\mathrm{wf}} & \forall t, x=0 \\ T_{\mathrm{f}}(x, y, z, t)=T_{\mathrm{res}} & \forall x, \forall z, \forall t, y \rightarrow y_{\max } \\ T_{\mathrm{w}}=T_{\mathrm{b}} & 0<x<L_{\mathrm{f}}, y=w / 2,-H_{\mathrm{f}} / 2<z<H_{\mathrm{f}} / 2 \end{array}\right.$
（5）式中，第1个表达式表示任意时刻缝口注入速度恒定为qin；第2个表达式表示任意时刻缝口温度都等于井底温度Twf，井底温度Twf利用井筒温度场模型进行计算[5]；第3个表达式表示物理模型外边界温度等于原始储层温度Tres；第4个表达式表示裂缝壁面内侧的流体温度Tw与裂缝壁面外侧的基质岩体温度Tb相等。
采用有限元方法对裂缝、地层岩体温度场进行求解。将求解域离散为单元体，对基质岩体采用四节点四面体单元，对裂缝采用三角形单元。将裂缝单元体的温度场Tw,e表示为节点温度的插值关系[10]，即：
$ T_{\mathrm{w}, \mathrm{e}}(x, y, z)=\boldsymbol{N}(x, y, z) \boldsymbol{q}_{\mathrm{T} \mathrm{P}}$
形函数矩阵的计算式为[11]：
$ N_{i}=\frac{1}{6 V}\left(a_{i}+b_{i} x+c_{i} y+d_{i} z\right)$
采用伽辽金法[12]对（2）式中的裂缝温度场进行离散，则有：
$ \boldsymbol{R} \boldsymbol{T}_{\mathrm{w}}+\boldsymbol{S} \frac{\partial \boldsymbol{T}_{\mathrm{w}}}{\partial t}=\boldsymbol{F}$
矩阵R、S、F的元素由相应的单元矩阵元素组装而成，单元矩阵元素为：
$ \begin{aligned} R_{\mathrm{e}, i, j}= & \int_{\Omega_{e}} \lambda_{\mathrm{w}}\left(\frac{\partial N_{i}}{\partial x} \frac{\partial N_{j}}{\partial x}+\frac{\partial N_{i}}{\partial y} \frac{\partial N_{j}}{\partial y}+\frac{\partial N_{i}}{\partial z} \frac{\partial N_{j}}{\partial z}\right) \mathrm{d} \Omega+ \\ & \int_{\Omega_{\mathrm{e}}} c_{\mathrm{w}} \rho_{\mathrm{w}} K_{\mathrm{f}} N_{j}\left(\frac{\partial N_{i}}{\partial x} \frac{\partial p}{\partial x}+\frac{\partial N_{i}}{\partial y} \frac{\partial p}{\partial y}+\frac{\partial N_{i}}{\partial z} \frac{\partial p}{\partial z}\right) \mathrm{d} \Omega \end{aligned}$
$ S_{\mathrm{e}, i, j}=\int_{\Omega} c_{\mathrm{w}} \rho_{\mathrm{w}} N_{i} N_{j} \mathrm{~d} \Omega $
$ F_{\mathrm{e}, i}=\int_{\Omega} N_{i}\left[\frac{\lambda_{\mathrm{s}}}{\delta}\left(T_{\mathrm{r}}-T_{\mathrm{w}}\right)+m_{\mathrm{pf}} \Delta H_{\mathrm{R}} \frac{\mathrm{~d} \alpha}{\mathrm{~d} t} \frac{1}{\mathrm{~d} x \mathrm{~d} y \mathrm{~d} z}\right] \mathrm{d} \Omega$
（8）式中对时间求偏导的目的是从tn时刻对应的Tw,n推出tn+1时刻对应的未知量Tw,n+1，Tw,n在tn时刻为已知，并且在较短的时间步长Δt内，F也为已知，故对（8）式进行处理可得：
$ \begin{array}{l} {\left[\frac{2 \boldsymbol{S}\left(\boldsymbol{T}_{\mathrm{w}, n+1}, t_{n+1}\right)}{\Delta t}+\boldsymbol{R}\left(\boldsymbol{T}_{\mathrm{w}, n+1}, t_{n+1}\right)\right] \boldsymbol{T}_{\mathrm{w}, n+1}=} \\ \boldsymbol{S}\left(\boldsymbol{T}_{\mathrm{w}, n+1}, t_{n+1}\right)\left(\frac{2}{\Delta t} \boldsymbol{T}_{\mathrm{w}, n}+\frac{\boldsymbol{T}_{\mathrm{w}, n+1}-\boldsymbol{T}_{\mathrm{w}, n}}{\Delta t}\right)+\boldsymbol{F}_{n+1} \end{array}$
采用同样的方法对（3）式进行离散得：
$ \boldsymbol{G} \boldsymbol{T}_{\mathrm{r}}+\boldsymbol{H} \frac{\partial \boldsymbol{T}_{\mathrm{r}}}{\partial t}=\boldsymbol{E}$
将（13）式进行处理得：
$ \begin{array}{l} {\left[\frac{2 \boldsymbol{H}\left(\boldsymbol{T}_{\mathrm{r}, n+1}, t_{n+1}\right)}{\Delta t}+\boldsymbol{G}\left(\boldsymbol{T}_{\mathrm{r}, n+1}, t_{n+1}\right)\right] \boldsymbol{T}_{\mathrm{r}, n+1}=} \\ \quad \boldsymbol{H}\left(\boldsymbol{T}_{\mathrm{r}, n+1}, t_{n+1}\right)\left(\frac{2}{\Delta t} \boldsymbol{T}_{\mathrm{r}, n}+\frac{\boldsymbol{T}_{\mathrm{r}, n+1}-\boldsymbol{T}_{\varepsilon, n}}{\Delta t}\right)+\boldsymbol{E}_{n+1} \end{array}$
2.2 模型验证
由于本文模型对基质岩体采用的是四节点四面体单元，对裂缝采用的是三角形单元，单元节点难以完全准确落在裂缝中轴线上；而一维温度场模型采用差分求解，各个单元的边相互垂直，单元节点可以全部分布在中轴线上。因此，在提取本文模型裂缝中轴线上的温度值时，裂缝高度方向上的坐标可能会略微偏移中轴线，导致本模型温度曲线不如一维模型的温度曲线平滑。减小单元尺寸，可以让温度曲线更加平滑，考虑到三维模型计算量较大，本文在增加曲线平滑度和减少计算工作量之间寻求平衡，在保证计算精度的前提下，选择合适的单元尺寸。对比两种模型裂缝中轴线上的温度分布发现，两种模型计算结果存在一定的差异（见图6 ），但差异不大，说明本文模型具有较高的可靠性。
3 模型应用与案例分析
采用表2 所示的参数计算不同条件下裂缝温度场，分析不同因素对裂缝温度场的影响规律。
3.1 相变热对裂缝温度的影响
相变热指单位质量的PF放出的热量（见表1 ），在升温速率为5，10，20 K/min条件下测试得到的相变热分别为32.07，31.40，32.10 J/g，平均值为31.86 J/g。为便于分析相变热对裂缝温度场的影响规律，计算相变热分别为0，20，40，60 J/g条件下的温度分布（见图7 ）。结果表明：相变热越大，缝内温度相对越高，且不同相变热条件下缝口附近温度差值较大，越靠近裂缝远端，温度差值越小。注液时间为5，20，40，60 min时，不同相变热条件下，中轴线上距离缝口约0，7，10，10 m处出现最大温度差值。出现这种现象的原因在于：①注液时间为5 min时，缝内温度整体较高，PF在缝口大量相变，相变热造成缝口温度升高（见图7a ）。②注液时间为20 min时，缝口温度较低，在缝长方向上温度逐渐升高，缝口温度低于PF相变的临界温度值，因此，缝口处的PF相变率较低，缝口温度上升幅度有限，而在缝内7 m位置处，裂缝温度上升至足以使PF大量相变的临界值，在此处PF大量相变，相变热引起裂缝温度大幅上升，超过7 m范围后，残留的PF较少，相变体积也较少，因此，相变热引起的温度变化幅度降低（见图7b ）。持续注入较低温度的压裂液不会造成缝内温度无限制降低，而是逐渐达到稳定[17]。③注液时间为40，60 min时，不同相变热条件下缝内温度差距不大；缝口温度较低，PF相变量较少，相变热引起的温度升高幅度较小，在缝内10 m位置处，裂缝温度才能上升至足以使PF大量相变，超过10 m范围后，残留的PF较少，相变体积也较少，因此，相变热引起的温度变化幅度降低（见图7c 、图7d ）。
为验证上述推论，导出了相变热为40 J/g时不同时间PF相变转化率（见图8 ），即每个网格单元位置处PF转化率的绝对值，该值直接由（1）式计算得到。同时导出了相变热为40 J/g时不同时间PF转化率的变化率（见图9 ），该值由图8 中长度方向相邻两个单元的总转化率相减得到。图8 表明，注液时间为5，20，40，60 min时，缝内中轴线上距离缝口约15，40，48，52 m处PF的转化率接近1.0。裂缝中轴线上温度越低，PF转化率为1.0的等值线离缝口越远。图9 所示的PF转化率的变化率反映了在某一位置发生相变的PF总体积，也反映了在该位置处释放的热量和温度升高幅度。从转化率的变化率来看，注液时间为5，20，40，60 min时，裂缝中轴线上转化率的变化率最大值出现在距离缝口约0，7，10，10 m处，因此，在这些位置的温度升高幅度最大，这验证了图7 所示的温度计算结果。
3.2 PF体积分数对裂缝温度的影响
PF与NPF共同组成PFFS，在现场压裂施工过程中，PF和NPF分别通过两条管线注入到井筒中，针对不同的储层，PF与NPF的体积分数不同[18]，为研究PF体积分数对裂缝温度场的影响，基于表2 所示的数据，将PF体积分数分别设置为10%，20%，30%与40%，计算分析PF体积分数对裂缝中轴线温度的影响（见图12 ）。由图可知，PF体积分数越高，PF相变时释放的热量越大，因此缝内温度升高幅度也就越高。注液时间为5 min时，PF体积分数对温度差值的影响主要出现在缝口，PF体积分数为10%，20%，30%和40%时，缝口温度分别为88.1，89.9，91.7，93.5 ℃。注液时间为20，40，60 min时，温度差值主要出现在靠近缝口的位置，在裂缝远端，缝内温度等于储层原始温度。
3.3 比热容对裂缝温度的影响
PFFS的比热容决定了相同热量变化情况下的温度变化幅度。基于表2 数据，将PFFS比热容分别设置为1 700，2 200，2 700，3 200 J/（kg·K），其中基准值采用表2 中所示的2 200 J/（kg·K），计算分析比热容对裂缝温度场的影响（见图15 ）。由图可见，温度较低的PFFS注入地层后，一方面，PFFS从地层吸收热量，温度升高，另一方面，PFFS相变放热，导致温度上升。因此，比热容越小，缝内温度越高。注液时间为5 min时，比热容为1 700，2 200，2 700，3 200 J/（kg·K）条件下，缝口温度分别为93.4，91.8，90.8，90.1 ℃，随着注液时间的增加，缝口温度逐渐降低。
定义比热容2 200 J/（kg·K）为基准值，将比热容为1 700，2 700，3 200 J/（kg·K）条件下的裂缝中轴线上的温度与基准值条件下的温度进行对比（见图16 ）。因注入压裂液的温度比储层原始温度低，压裂液在储层岩石传递热量和相变放热的影响下温度升高，比热容越大，温度升高幅度越小，因此，当比热容为1 700 J/（kg·K）时，裂缝中轴线上压裂液的温度高于基准条件下的温度，温度差值为正数；而比热容为2 700，3 200 J/（kg·K）时，裂缝中轴线上压裂液的温度低于基准条件下的温度，温度差值为负数。
4 结论
相变热对裂缝温度的影响不可忽略，相变过程主要受温度控制，准确预测裂缝温度是相变自支撑压裂成功的关键。与未考虑相变热条件下的温度相比，当相变热为20，40，60 J/g时，注液结束时刻，缝内温度的上升最大值分别为2.1，4.2，6.2 ℃。相变热和PF体积分数与裂缝温度变化呈正相关，比热容与温度变化呈负相关。
缝内不同位置、不同时间温度受注入冷流体的冷却效应和PF相变过程中的放热效应交替主导。压裂液注入初期，缝内温度较高、相变速度较快，相变放热对储层岩石温度影响较大；注入后期，缝内温度降低，相变放热速度减慢、压裂液对储层岩石的冷却效应增强。在注入40，60 min时，最大相变速率和最大温差出现在距缝口约10 m处。
符号注释：
a，b，c，d——与节点几何位置相关的系数，无因次；cmw——基质内流体的比热容，J/（kg·K）；cr——岩体的比热容，J/（kg·K）；cw——缝内压裂液的比热容，J/（kg·K）；E——岩体温度场的温度载荷向量；F——裂缝温度场的温度载荷向量；Fe,i——单元裂缝温度场温度载荷向量的元素；G——岩体温度场的热传导矩阵；Hf——裂缝高度，m；H——岩体温度场热容矩阵；ΔHR——PF固化反应放热量，J/kg；i，j——序号；Kf——裂缝渗透率，m2；Lf——裂缝长度，m；mpf——六面体裂缝微元体中包含的PF的质量，kg；mpm——六面体岩体微元体中包含的PF的质量，kg；n——时间步；Ni——矩阵N中的元素；N(x,y,z)——插值矩阵；p——裂缝内流体压力，Pa；q——边界注入速度，m3/min；qin——恒定注入速度，m3/min；qT,e——节点温度矩阵，K；R——裂缝温度场热传导矩阵；Re,i,j——单元裂缝温度场热传导矩阵的元素；S——裂缝温度场热容矩阵；Se,i,j——单元裂缝温度场热容矩阵的元素；t——时间，s；T——温度，K；Tb——边界温度，K；Tr——岩体温度，K；Tres——储层原始温度，K；Tr——岩体温度场矩阵；Tw——裂缝壁面内侧温度，K；Tw——裂缝温度场矩阵；Tw,e——裂缝单元体的温度，K；Twf——井底温度，K；vx，vy，vz——沿不同坐标方向的渗流速度，m/s；V——单元体体积，m3；w——裂缝宽度，m；x，y，z——直角坐标系，m；ymax——缝宽方向无穷远处位置，m；α——转化率，无因次；$ \alpha(t)$——t时刻的转化率，无因次；δ——裂缝半宽，m；λr——压裂液与岩体之间的换热系数，W/（m2·K）；λx，λy，λz——沿不同坐标方向的岩石导热系数，W/（m·K）；λw——导热系数，W/（m·K）；μ——缝内流体黏度，Pa·s；ρmw——基质内液体的密度，kg/m3；ρr——岩体的密度，kg/m3；ρw——压裂液密度，kg/m3；ϕ——基质岩体孔隙度，无因次；Ω——积分域；Ωe——单元积分域。