化石能源属于不可再生能源，储量有限，且使用过程会导致众多环境污染问题，因而地热能、风能、太阳能等清洁能源日益被重视. 地热资源贮存在地表以下，储备十分丰富^[1]. 长期以来地热能源开发多以水热型浅层地热能为主，对于埋藏较深的干热岩的地热资源的开发利用较少^[2]. 深部岩石多为无天然裂隙的致密岩石^[3-5]，在开发过程中须通过人工刺激在岩层中开辟渗流通道，因此地热开发中存在裂隙岩体中的流固换热问题.

学者们多采用室内试验和理论分析对裂隙岩体中的流固换热问题进行研究. 在渗流换热室内试验方面，近年来国内外学者多制备粗糙裂隙岩样，探索应力状态、岩体温度、裂隙开度等变量与裂隙岩样渗流换热特性的关系^[6-8]. 赵坚^[9]在进行花岗岩裂隙的水热特性试验时发现裂隙几何与表面粗糙度在岩石裂隙热传导中起着重要作用. Bai等^[10]开展不同围压、温度、体积流量等条件下的换热试验，结果表明温度、流速、围压以及裂隙水力开度都会对换热效果造成影响. 李正伟等^[11]制备光滑平直裂隙与粗糙裂隙岩样开展花岗岩试样换热特性试验，研究结果表明水-岩之间的换热强度与体积流量、岩石温度正相关，且粗糙裂隙对流强度较光滑平直裂隙的有所提高. 黄奕斌等^[12]利用裂隙粗糙度（joint roughness coefficient，JRC）模板和3D打印技术，设置具有特定粗糙度的裂隙试样进行换热实验，发现裂隙面粗糙度会影响水-岩换热特性，在同等条件下顺水流方向粗糙度的增大会强化换热.

在渗流换热理论分析方面，学者们通过数值模拟进行渗流换热实验，探讨注入温度、体积流量、热物理参数等因素对于岩石换热特性的影响^[13-14]，研究成果与换热试验所得结论互为佐证. 换热系数作为表征岩石裂隙流体渗流换热特性的重要参数也是学者们研究的重点. 牛顿提出牛顿冷却定律，Chapman基于热边界层理论发展推导了换热系数的简单计算公式^[15]. 后人多利用局部热非平衡法、二维导热方程对于边界条件进行不同程度的简化，采用裂隙内水温线性分布假设或者平行板内温度分布解析模型，推导出裂隙换热系数不同形式的解析解，这也导致了不同的公式存在一定的适用条件，当超出一定条件时理论值与真实情况会存在较大偏差^[16-20]. Bai等^[16]以进出口水温均值近似代替进出口裂隙面温度均值，提出简单的换热系数计算公式. 朱家玲等^[17]采用局部热非平衡法和二维导热方程，推导出圆柱形试样裂隙换热系数的解，并对影响换热系数的因素进行敏感性分析，发现在同等条件下流速的影响大于裂隙开度的. Zhang等^[18]利用二维热传导方程提出计算圆柱形试样裂隙换热系数的公式. Jiang^[19]总结前人提出的几种计算换热系数公式，全面分析不同公式的适用性，对于公式中出现的数据异常，根据误差传播理论进行了解释. Zhao^[20]基于光滑平行板内温度分布解析模型，分别利用局部热平衡和热非平衡对边界条件进行简化，推导换热系数的2种解析解，可进行三维长方体模型的换热系数计算. 就渗流换热机理方面，研究表明，低温流体与其所接触的高温岩体裂隙壁面之间存在着渗流换热耦合作用. 一方面，流体的渗流形态会影响其热对流和热传导；另一方面，流固换热会导致流体热物理性质发生改变，如运动黏度、比热容，进而影响渗流形态^[21-22].

上述渗流换热试验多基于岩石三轴伺服系统改装而成，可进行不同温度以及围压下的渗流换热试验，但其换热介质在进入裂隙前后仍被加热，难以精确测定进出口水流温度和岩石表面温度，其试验数据存在一定误差，使得理论分析得出的结论也出现一定的偏差^[23-24].

本研究研发可以解决上述问题的岩石裂隙水-岩渗流换热试验系统，制备裂隙花岗岩试样，开展粗糙裂隙花岗岩在不同温度和体积流量下的水-岩换热试验，对出口流体温度与裂隙试样温度及体积流量的关系、渗流规律进行分析，并开展花岗岩水-岩换热系数的影响因素研究及其敏感性分析. 本研究在裂隙试样温度为70~100 ℃时进行探索性试验，以期为深层地热资源开发提供一定参考.

裂隙水-岩换热试验系统如图1所示，主要分为循环注水、温控加热、数据采集3个子系统.

循环注水子系统：Isco双泵系统，可以控制裂隙内通过流体的渗流体积流量. 协同工作可控体积流量为0~204 mL/min.

温控加热子系统：包括温控装置、上下表面2块加热板，通过温控装置控制施加特定温度对岩样进行加热，最高加热温度可达200 ℃；

数据采集子系统：主要由温度探头以及压差计组成. 在试样前后表面安装温度探头，可实时传回岩石表面温度，以便实现对于热边界条件的精准控制；在进出水口位置装有温度探头以及压差传感器，可以记录进出口流体温度以及压差.

该系统换热介质在进入岩样前后不会被加热，同时岩样表面布置多支温度传感器，皆提升了数据的准确性，同时采用长方体试样更有利于控制裂隙的位置、走向.

裂隙水-岩换热试验流程图如图2所示.

1）采用100 mm×100 mm×200 mm的长方体花岗岩试样，即裂隙面宽度为100 mm，长度为200 mm，利用岩石三轴伺服多场耦合试验系统，采用巴西劈裂法制备裂隙试样.

2）针对预先制备好的花岗岩裂隙试样，采用OKIO 3M三维形貌扫描系统对岩样裂隙面进行扫描，而后计算裂隙面坡度的一阶导数均方根Z₂以及岩样的裂隙粗糙度JRC. 表达式如下：

(1) $ {Z}_{2}={\left[\frac{1}{(n-1){(\Delta y)}^{2}}\sum _{i=1}^{n-1}{({z}_{i+1}-{z}_{i})}^{2}\right]}^{1/2}, $

(2) $ \mathrm{J}\mathrm{R}\mathrm{C}=32.2+32.47\mathrm{l}\mathrm{g}\;{Z}_{2} . $

式中：n为每一个剖面线上全部的采样点数量，z_i为第i个采样点的粗糙高度，∆y为采样点之间的间隔.

裂隙试样H1~H4裂隙粗糙度分别为14.871、12.852、10.675、9.048.

1）处理好的岩样无错动贴合将导致无法过流，因此将其沿长边方向错动0.2~1.0 mm后固定以获取不同水力开度的裂隙岩样，采用耐高温的硅橡胶对花岗岩试样裂隙进行密封，并将其放置在试验盒中. 为了保证封水性以及隔热性，在进出口处布置密封垫和隔热板.

2）在裂隙试样进出水口、侧面布置共6支温度探头，实时记录岩样表面以及进出口温度，并通过温度控制箱与计算机相连接. 在进出口位置接入压差计，并通过压差计接入计算机，随时监测流体进出口压差数据.

3）将裂隙试样进水口与Isco泵相连，以体积流量为100 mL/min向裂隙内注入常温蒸馏水，如不存在明显渗漏现象，则说明裂隙密封良好.

4）通过温控装置调整加热板温度略高于70 ℃使岩样外壁加热温度升高，等到传感器温度维持在约70 ℃（30 min内变化不超过1 ℃），则认为裂隙试样温度已经达到恒温加热70 ℃的效果.

5）在恒温后启动Isco泵设置体积流量数值为10 mL/min，当流动和换热状态稳定之后，记录此时的进出口水温以及压差. 调整体积流量在10~80 mL/min逐步提高（8个体积流量等级），进行8次不同体积流量下的换热试验.

6）在各体积流量等级测量完成后，关闭Isco泵，调整岩样外壁加热温度至80、90 、100 ℃，重复上述步骤3）~5）. 待一块岩样完成试验后，更换岩样重复所有，待所有岩样完成试验后整理仪器.

岩样在不同温度下的出口水流温度与体积流量之间的关系如图3所示. 图中， θ₀为裂隙试样温度，q_V为体积流量，Q_out为裂隙出口水流温度. 以岩样H1为例，体积流量从10 mL/min增大至80 mL/min，当裂隙试样温度为70 ℃时，出口水流温度下降了29.4%，与裂隙试样温度为80 、90、100 ℃时相比，分别下降了31.3%、32.8%、33.8%. 在同一体积流量下，随着裂隙试样温度升高，出口流体温度也会越大. 不同温度下的出口水流温度与体积流量关系线性拟合如图3所示，可知拟合后R²皆为0.99，线性程度较高. 在相同的裂隙试样温度下，体积流量增大，出口水流温度随之减小. 这是因为体积流量越大，渗流水与岩石接触的时间越少，即热传递时间越少，出口水流温度越低. 在同一体积流量下随着裂隙试样温度的升高，换热量越高，相同时间内流体温度升高越明显.

其余岩样的出口水流温度与体积流量间关系基本与H1的表现一致，皆为拟合关系较好的线性关系. 不过不同岩样之间表现出一定差异性，随着试样裂隙粗糙度的降低，出口水流温度降低. 这是因为裂隙粗糙度降低，换热接触面积减少，换热量降低.

Forchheimer方程可以较好地描述岩石裂隙中流体压力梯度和体积流量的非线性关系^[25-27]，其表达式如下：

(3) $ - \nabla {{p}} = A{q_V}+B{q_V}^2. $

式中： $ \nabla p $为压力梯度，A、B分别表示流体黏性影响和惯性影响系数. 当流体的体积流量足够小时，此时由惯性影响造成的压力损失较小，即Bq_V²项较小，相比黏性影响可以忽略，此时Forchheimer方程即和立方定律等价.

不同温度下裂隙梯度随体积流量的变化如图4所示. 图中，拟合曲线是通过Forchheimer方程对数据点拟合得到的. 可以看出，随着体积流量增大，压力梯度增大，流体非线性程度增大. 在现有研究中多引入非线性因子α=Bq_V²/（Aq_V+Bq_V²）定量表征裂隙内流体非线性程度，其定义为流体流动过程中由惯性效应造成的压降损失与总的压降损失的比值，当α>0.1时即为非线性流^[28]. 随着体积流量增大，非线性因子逐渐增大. 本研究中当体积流量大于20 mL/min时，非线性因子α>0.1，流态转为非线性.

换热系数反映在流固接触单位温度差下，单位换热面积上在单位时间内所能传递的热量大小，是岩体裂隙换热评价的重要分析参数. Zhao^[20]采用解析方法，根据局部热非平衡理论结合换热方程和相应边界条件，推导出岩体裂隙壁面温度的分布情况，然后基于岩体裂隙壁面温度进一步推出换热系数的表达式：

(4) $ h=-\dfrac{{\rm{ln}}\;\dfrac{{\theta }_{\mathrm{o}\mathrm{u}\mathrm{t}}-{\theta }_{0}}{{\theta }_{\mathrm{i}\mathrm{n}}-{\theta }_{0}}{C}_{\mathrm{w}}{\rho }_{\mathrm{w}}u{e}_{\mathrm{h}}{\lambda }_{\mathrm{r}}/2}{{\rm{ln}}\;\dfrac{{\theta }_{\mathrm{o}\mathrm{u}\mathrm{t}}-{\theta }_{0}}{{\theta }_{\mathrm{i}\mathrm{n}}-{\theta }_{0}}{C}_{\mathrm{w}}{\rho }_{\mathrm{w}}u{e}_{\mathrm{h}}R/2+{\lambda }_{\mathrm{r}}L}. $

式中：R为裂隙面宽度的一半，R=50 mm；L为特征长度，通常取岩样流动方向的长度，取200 mm； $ {\lambda }_{{\rm{r}}} $为岩样的导热系数，为2.9 W/(m·K）；θ₀为裂隙试样温度；θ_in为裂隙进口水温；u为流体流速；e_h为裂隙水力开度； $ {\rho }_{\mathrm{w}} $为水的密度，为1000 kg/m³； $ {C}_{\mathrm{w}} $为水的比热容，为4200 J/(kg·℃).

裂隙经典立方定律表达式 ^[29]如下：

(5) $ {q_V} = - \frac{{w{e_{\text{h}}}^3}}{{12\mu }}\nabla p . $

式中：w为裂隙面宽度，为100 mm；μ为流体的动力黏滞系数.

立方定律将岩石裂隙假设为光滑平行板，且仅考虑流体流动状态为层流，实际上并不存在这样的理想化裂隙. 但是天然岩石裂隙粗糙度不一以及开度分布不均匀使得真实开度难以得到精确的测量，因此学者们提出利用裂隙水力开度进行分析. 其定义为与岩石天然粗糙裂隙拥有相同过水能力时理想光滑平行板裂隙的开度^[29-32]，可以通过式（5）进行反算，计算可知裂隙试样H1~H4裂隙水力开度分别为72.2、157.9、235.2、345.4 μm.

在探究裂隙试样温度和渗流体积流量与换热系数之间关系的过程中，以试样H1为例得到不同的裂隙试样温度下换热系数h随渗流体积流量的变化，如表1所示.

不同体积流量下换热系数随裂隙试样温度的变化如图5所示. 可以看出，当裂隙试样温度相同时，换热系数随着体积流量的增大而增大，70、80、90、100 ℃条件下体积流量从10 mL/min增大到80 mL/min，换热系数分别增加了551.5%、568.4%、469.3%和411.5%. 当渗流体积流量相同时，换热系数将随着温度升高而增加. 与70 ℃相比，80、90、100 ℃条件下体积流量为10 mL/min时增幅分别为25.2%、51.8%、82.2%. 根据热边界层理论，流体的运动黏滞系数减小以及流速的增大均会影响边界层厚度，使其减小，从而导致裂隙流体换热的阻力降低，换热强度提高，换热系数增大.

为了进一步考察温度对对流换热系数的影响，绘制表1中80、90、100 ℃下的对流换热系数与70 ℃下的对流换热系数的差值（即Δh(θ)=h(θ₀)−h(70 ℃)）随渗流体积流量的变化曲线，如图6所示. 可以看出，换热系数差值均为正值，且同体积流量下裂隙试样温度越高差值越大，进一步验证了对流换热系数与裂隙试样温度的正相关关系. 同温度时该差值与渗流体积流量表现为非严格线性，具有一定波动，但整体上换热系数差值与渗流体积流量近似为线性关系，表明同条件下体积流量越大，对换热系数产生的影响越大.

试样H2~H4换热系数规律与试样H1的基本一致，因此不再讨论.

归一化换热系数（h/h₀）和归一化渗流流速（u/u₀）间的关系可以表征换热系数与渗流流速之间的关系，研究表明其近似符合幂函数关系^[17]，表达式如下：

(6) $ h/{h}_{0}=a{(u/{u}_{0})}^{b}. $

式中： ${u}_{0} $和 $ {h}_{0} $分别为初始（体积流量为10 mL/min）流速和初始（体积流量为10 mL/min）对流换热系数，a、b为通过拟合得到的参数.

如图7所示为不同裂隙试样的归一化换热系数和归一化渗流流速之间的关系. 可以看出，拟合曲线较好地满足了如式（5）所示的幂函数关系，参数a的取值为0.90~1.14，参数b的取值为0.81~0.90，拟合曲线相关系数均大于0.90，拟合效果较好. 拟合曲线参数b的取值均小于1.0，说明换热系数与渗流流速之间并非简单的线性关系，而是随着流速的增大，换热系数一直增大，但增幅越来越小.

对换热问题中的重要无量纲特征数努塞尔数Nu、普朗特数Pr、雷诺数Re进行无量纲分析.

研究表明，一般情况下可以用幂函数形式来表示各无量纲数之间的关系^[33-34]：

(7) $ Nu/P{r^{{\text{1/3}}}} = CR{e^m}. $

研究Nu/Pr^1/3与Re之间的关系，拟合曲线如图8所示. 拟合曲线相关数皆大于0.92，表明拟合效果良好. 如图8所示，Re增大，Nu也会随之增大，这是因为渗流流速越大，热传导越强，这进一步验证了归一化换热系数（h/h₀）和归一化渗流流速（u/u₀）的分析结论.

以裂隙试样温度为100 ℃时得到的数据为例，分析该裂隙试样温度下不同体积流量时对流换热系数随水力开度的变化规律，如图9所示. 可以看出，在相同裂隙试样温度和渗流体积流量下，随着裂隙水力开度的增大，总体上对流换热系数基本表现为线性减小，但是变化的幅度不大，水力开度由72.2 μm增至345.4 μm，减小幅度小于20%. 这是因为在相同渗流体积流量下，裂隙水力开度的增大导致裂隙过水断面面积增大，从而使渗流流速降低. 根据热边界层理论，流速降低使热边界层厚度增大，对流换热热阻增大，换热系数减小.

以裂隙试样温度为100 ℃时得到的数据为例，分析该裂隙试样温度下不同体积流量时对流换热系数随JRC的变化规律，如图10所示. 可以看出，随着裂隙粗糙度JRC的增大，换热系数基本表现为线性增大关系. 相关性系数R²皆大于0.80，拟合效果较好. 但是随着体积流量的增大，尤其是当体积流量大于60 mL/min后，线性拟合系数会低于0.90，换热系数与裂隙粗糙度JRC线性相关性降低，这是因为当体积流量大于一定数值时，换热介质的流态非线性程度较高，使得换热系数与岩石裂隙JRC数值之间不再是单纯的线性关系.

对换热系数进行敏感性分析的基本原理如下：对于一个系统的特性Z，其主要取决于n个影响因素 $ \alpha =\left\{{\alpha }_{1},{\alpha }_{2},{\alpha }_{3},\cdots ,{\alpha }_{{n}}\right\} $，即 $Z=f\left\{{\alpha }_{1}, {\alpha }_{2},{\alpha }_{3},\cdots ,{\alpha }_{{n}}\right\}$. 系统在某基准状态 $ {\alpha }^{*} $下 ${Z}^{*}=f\left\{{{\alpha }_{1}}^{*},{{\alpha }_{2}}^{*}, {{\alpha }_{3}}^{*},\cdots , {{\alpha }_{{n}}}^{*}\right\}$. 敏感性分析即是基于此基准状态，分析当影响特性Z的某一因素发生变化时，其偏离基准状态 $ {Z}^{*} $的趋势以及程度^[17]. 定义敏感性函数如下：

(8) $ S\left({\alpha }_{{k}}\right)=\left|\frac{{\rm{d}}Z\left({\alpha }_{{k}}\right)}{{\rm{d}}{\alpha }_{{k}}}\right|\frac{{\alpha }_{{k}}}{Z}. $

通过对不同因素下的特性进行拟合得到敏感性函数，而后将不同因素的基准状态代入到式（8）中即可得到不同因素的敏感性因子. 敏感性因子越大，表明该因素对特性的影响越显著.

首先对渗流体积流量q_V和裂隙水力开度 $ {e}_{\mathrm{h}} $进行敏感性分析，即Z=f{q_V，e_h}，特性Z即为换热系数h.

以裂隙试样温度为70 ℃时为例. 裂隙体积流量q_V变化区间为10 ~80 mL/min，裂隙水力开度e_h变化区间为0.0722~0.3454 mm，选择基准状态为 $q_V^* $=50 mL/min， $e_{\rm{h}}^* $=0.1579 mm. 当q_V= $q_V^* $=50 mL/min时，拟合h=f (e_h)；当e_h= $e_{\rm{h}}^* $=0.1579 mm时，拟合h=f (q_V)利用origin软件进行多项式拟合，拟合结果如图11、12所示.

拟合方程分别如下：

$ \begin{gathered} h({q_V}) = - 132.295+29.971{q_V} - 0.401{q_V}^2+0.002\;18{q_V}^3, \\ h\left( {{e_{\rm{h}}}} \right) = 857.728 - 2.282{e_{\rm{h}}}+0.006\;2{e_{\rm{h}}}^2 - 6.65 \times {10^{ - 6}}{e_{\rm{h}}}^3. \\ \end{gathered} $

依据式（8）计算得到q_V、e_h分别对应的敏感性函数分别为

$ S\left( {{{{q}}_{{V}}}} \right) = \left| {\frac{{29.971{{{q}}_{{V}}} - 0.802{{{q}}_{{V}}}^2+0.006\;54{{{q}}_{{V}}}^3}}{{ - 132.295+29.971{{{q}}_{{V}}} - 0.401{{{q}}_{{V}}}^2+0.002\;18{{{q}}_{{V}}}^3}}} \right|, $

$ S({e_{\rm{h}}}) = \left| {\frac{{ - 2.282\;2{e_{\rm{h}}}+0.012\;46{e_{\rm{h}}}^2 - 1.995 \times {{10}^{ - 5}}{e_{\rm{h}}}^3}}{{857.728 - 2.282{e_{\rm{h}}}+0.006\;2{e_{\rm{h}}}^2 - 6.65 \times {{10}^{ - 6}}{e_{\rm{h}}}^3}}} \right|. $

代入基准状态得到θ₀=70 ℃条件下，S( $q_V^* $)=0.489，S( $e_{\rm{h}}^* $)=0.204. 按照同样的方法得到θ₀=80、90、100 ℃时换热系数关于体积流量和水力开度的敏感性因子，如表2所示. 可以看出，随着温度的升高，体积流量的敏感性越来越大，温度升高导致流体运动黏度降低，使得热边界层厚度增大影响换热系数. 水力开度的敏感性会先增大而后减小，这可能是因为当温度上升到一定程度以后，温度对于换热系数的影响会远大于水力开度造成的影响，从而导致换热系数会出现突然的增大. 而在其他条件一致的情况下，换热系数关于渗流体积流量的敏感性因子S( $q_V^* $)始终大于关于水力开度的敏感性因子S( $e_{\rm{h}}^* $)，说明换热系数受渗流体积流量的影响大于受水力开度的.

考虑温度的影响，继续对裂隙体积流量q_V和裂隙试样温度θ₀进行敏感性分析，选择基准状态为 $q_V^* $=50 mL/min， $\theta _0^* $=100 ℃. 得到H1~H4裂隙换热系数h关于体积流量和裂隙试样温度的敏感性因子，结果如表3所示. 随着裂隙水力开度的增大，换热系数对于渗流体积流量和温度的敏感性均会提高. 这是因为裂隙水力开度增大，裂隙过水断面面积增大，过水能力增强，流速相应降低. 据热边界层理论，流速降低使热边界层厚度增加，对流换热热阻增大，换热系数减小. 在同条件下，S( $\theta _0^* $)基本都大于S( $q_V^* $)，则认为换热系数受裂隙试样温度的影响较受渗流体积流量的更大.

综上，就q_V、θ₀、e_h三个因素而言，换热系数h受裂隙试样温度θ₀的影响最大，其次是裂隙体积流量q_V，影响最小的是水力开度e_h.

（1）当体积流量与裂隙水力开度条件保持不变时，裂隙试样温度升高，出口水流温度线性升高；换热系数随温度升高而线性增大.

（2）在裂隙试样温度与裂隙水力开度保持不变的条件下，体积流量增大，出口水流温度随之线性减小；换热系数随体积流量的增大而增大，表现为指数小于1的幂函数关系.

（3）在相同裂隙试样温度和渗流体积流量下，随着裂隙水力开度的增大，换热系数随之线性减小，水力开度由72.2 μm增至345.4 μm，裂隙试样温度为100 ℃时减幅小于20%. 岩石裂隙JRC增大，换热系数线性增大.

（4）利用敏感性函数对换热系数进行影响因素的敏感性分析，换热系数受裂隙试样温度的影响最大，其次是裂隙渗流体积流量，影响最小的是水力开度.

（5）实际地热开发中深部岩石温度较高，本研究仅针对较低温条件展示试验，后续将对于实验仪器进一步改进，进而开展较高温下的裂隙渗流换热试验.