多年冻土约占世界陆地面积的25%，而在我国多年冻土约占国土面积的22.4%，主要分布在青藏高原和东北地区^［1-2］。青藏高原多年冻土对气候及环境变化敏感，为减少全球气候变暖和工程扰动所导致的融沉病害^［3-5］，以块石路基为主的“主动冷却路基”^［6］因其具有经济，易于施工等优点，可通过调控寒季强迫对流和暖季弱自然对流效应的方式将冷能储存到下伏多年冻土中来实现热防护^［7］，而被广泛应用于多年冻土路基工程中。对路基温度场的长期监测结果表明，块石结构路基作为一种有效的设计方案，其保持了路基工程完整性并最大限度减少了融沉与冻胀的影响，同时块石基底路基、U型块石路基对下伏多年冻土有着较好的降温效果^［8-10］。在传统块石结构路基模型传热的计算过程中，考虑了外界风的强迫对流效应，基于连续性方程、非达西流动量方程，以及假设多孔介质中固体骨架和临近孔隙流体温度时时相等的局部热平衡能量方程（LTE），建立出多年冻土区开放块石路基对流换热数值模型^［11-12］。这类模型的建立为系统评估块石结构长期换热效应奠定了坚实的基础。现有研究结果表明^［13］，多孔介质的局部热平衡理论的有效性受到流固有效导热系数比（k）、毕渥数（Biot）等一系列参数影响，当多孔介质中流体与固体骨架的热学性质差异逐渐增大时，同一尺度内固体骨架与邻近孔隙流体之间的温差也随之扩大，局部热平衡也将发展为局部非热平衡。青藏铁路块石基底路基所采用的块石粒径约20~30 cm，有效导热系数比可达1/20^［14］，块石固体骨架与孔隙内部空气之间不可避免地存在一定的温差，局部非热平衡效应不可忽视。在一些多孔介质传热研究过程中均发现了非热平衡效应的存在，如开孔泡沫^［15］、烧结金属等微孔隙-小粒径多孔介质内^［16］，亦或是热干岩储层热回收^［17］、液氮对煤矿防火灭火问题^［18-19］、废热回收系统^［20］等方面，使用LTNE能量方程将使模拟结果精度高于LTE模型。因此，块石这种大孔隙、高渗透率结构内部的局部热非平衡效应如何？在计算块石结构对下部土体热影响时，采用LTNE能量方程相比于传统LTE能量方程的计算精度如何？仍有待进一步研究。

如图1所示，与假设多孔介质内部同一时间流固温度局部相等的热平衡理论有所区别，局部非热平衡理论认为多孔介质中的固体骨架与孔隙流体存在一定温差（T_s ≠T_f ）并因此发生能量传递（q_sf ），而该能量传递由LTNE能量方程中的重要参数流固换热系数（h_sf ）与比表面积a控制^［21］。在多孔介质中，固体骨架温度与流体温度通常并非时刻相等。LTE作为LTNE的一种理想状态^［22］，在流固温度相差不大时使用LTE理论进行数值模拟可以有效精简计算量。因此当我们使用LTNE理论进行多孔介质的传热计算时，准确求得ah_sf 的表达式至关重要。科研人员^［23-26］已对此展开丰富的研究并求得了一些经典的h_sf 表达式，对多孔介质局部非热平衡模型试验的研究作出重要的贡献。Xu等^［27］研究了微小多孔介质中微尺度效应对h_sf 的影响，发现h_sf 随粒径（d_p ）的减小和努塞尔数的增大而减小。Alazmi等^［28］研究发现LTNE模型中不同h_sf 之间的差异并不大，通常在进口区域和实验不确定性范围内更为明显。但对于小孔隙度和低雷诺数的情况，Achenbach^［23］的求得的换热系数与Wakao^［24］和Dixon等^［25］求得的换热系数有所不同。路朗^［29］在研究低导热系数多孔介质换热时，采用了Achenbach^［23］所求得换热系数公式，而在高导热系数多孔介质换热过程中采用了Wakao^［24］所求得换热系数公式；而Amiri等^［20］在研究高孔隙率-大粒径岩石组成的废热回收填充床中，采用了Dixon等^［25］求得的换热系数公式。尽管前人已对各换热系数进行了一定对比分析，但是仍然缺少针对类似块石路基一类的大粒径、高孔隙率的多孔介质内部不同h_sf 系统性的对比与分析，以及会对多孔介质内部传热产生何种影响。

本文基于多孔介质局部非热平衡理论，以北麓河非正线试验段开放块石基底路基信息为依托^［30］，通过数值模拟并结合一定量监测数据的方法，建立出块石基底路基局部热平衡和局部非热平衡数值模型，探求存在于块石层的随季节更替的局部非热平衡效应、分析LTE、LTNE模型的模拟结果中路基下部土体的温度变化趋势，以及LTNE模型内不同h_sf 在块石基底中的变化趋势及其对下部土体温度变化的影响。旨在建立出影响因素更全面的块石结构路基传热模型，从而更合理有效的预测青藏高原冻土路基的温度演化趋势，对于今后块石路基热质传递理论的研究起到积极的推动作用。

本文根据青藏铁路北麓河非正线试验段块石基底路基为原型进行建模^［30］。模型尺寸信息、土层分布及监测点布置方式如图2所示。路基顶面以及底面的宽度分别为7.6 m和18.7 m，路基高度为3.7 m，坡度为1∶1.5，块石层厚度为1.2 m。下伏土体从上至下分别为砂土、粉质黏土、强风化泥岩，其厚度分别为2.0 m、5.5 m以及22.5 m。通过仿真软件对瞬态计算结果记录、提取的功能实现“监测”的目的，在进行计算之前，首先通过在仿真软件中输入坐标的方式在块石夹层的左侧、中心和右侧布置监测点，其次在监测功能接口中设置不同的需要监测的变量后可获得瞬态计算结果中每个时间节点的多孔介质内流体温度、固体温度以及换热系数的变化趋势。与此同时，通过使用与上述相同的方法在块石层内三个监测点的正下方，地表以下2 m、5 m以及10 m处也相应设置了监测点，用于监测下伏多年冻土的温度变化趋势。

表1给出了数值模拟时使用到的控制方程，其中能量方程分为两种，即LTE与LTNE。LTE方程中，由于多孔介质固体骨架和流体在同一时间的温度局部相等，进而得到一个等效的“单相”能量守恒方程。而LTNE模型则相反，它对传热的计算基于两相之间的能量交换，从而存在两个能量方程：流体能量方程与固体能量方程。另外，本文探讨LTNE模型中几种换热系数h_sf 对块石结构换热的影响，其中所涉及到的h_sf 表达式均在此给出。

公式中涉及的重要参数在表2中给出。

块石、路基填土、下部土体以及空气的热力学参数在表3中给出。

考虑气候变暖影响，假设青藏高原在未来50年内平均气温将升高2.6 ℃^［31］，为精简模拟过程的计算量，根据附面层理论以及北麓河地区环境实测温度，本文使用的温度边界条件为基于实测温度拟合而成的温度正弦曲线［图3（a）］。可以看出整体上实测温度与拟合温度有着较好的契合性。但由于实际温度受环境因素影响存在一定程度波动，因此拟合得出的正弦形式的温度曲线与以天为单位测得的实测温度相比仍存在一定程度的误差。如图3（b）所示，可以发现在每年温度达到最高的7月以及最低的1月时间段内，拟合温度的误差偏高，其平均值在3.5~4.0 ℃之间，而在其他月份误差较低，平均值为2.69 ℃。北麓河地区的环境的风速变化规律可通过正弦函数表示，同时依据“综合幂次律”理论，环境中的风速随高度增加产生变化，模型右侧边界的风速可近似由指数形式描述。此外，在模型最底部，设置一个热流边界，其代表更深的土体带给模型底部的能量^［11］。具体边界条件在表4~5中列出，其中，A为年平均温度，B为温度变化幅值，在各个边界处的取值如表5所示。

为求解块石路基LTNE模型，在建立几何模型、剖分网格的过程结束后，通过仿真软件将模型中的块石域设置为多孔介质，在多孔介质窗口中对块石层的参数进行设置后开启局部非热平衡接口，此时仿真软件会在多孔介质域生成一组新的网格，在局部非热平衡接口中输入a、h_sf 等计算流固传热所需要的参数，仿真软件将依此对二者之间的换热进行计算。如图4所示，此时在块石域将同时存在两副网格，其中一副通过手动绘制生成，另一副则由软件自动复制生成，其分别代表块石域与空气域。两组网格剖分方式相同，但左右两侧边界形式有所区别：其中流体域网格为开放边界，与外部空气域相连，其在仿真软件中显示为绿色；而固体域网格为流固边界，在仿真软件中显示为灰色。在模拟结果中，可以根据需要去选择显示不同计算域的计算结果。图5给出了2013年1月和7月LTNE模型中块石路基的固体温度与流体温度，可以看出，在同一时间节点内，块石路基内部块石与空气的温度并不相等：在1月15日［图5（a）、（b）］），低温空气从右至左流经块石路基，温度从-14.8 ℃升高至-13.9 ℃；而块石的最低温为-13.5 ℃，最高温为-12.0 ℃，分别处在路基的内部以及左边缘处。在7月15日［图5（c）、（d）］，高温空气在流经路基后，温度从10.7 ℃降低至9.9 ℃；块石最低温和最高温分别8.9 ℃、13.0 ℃，同样处在块石路基的内部和左边缘处。

由于两个计算域同时存在于块石路基中，在同一监测点上也可以同时读取到两个计算域的模拟结果。图6给出了LTNE模型中块石路基中心位置监测点上多孔介质块石固体骨架和空气随时间变化的温度曲线，以及二者的温差。可以看出，空气的温度变化幅度大于块石，且块石温度变化相比于空气存在一定的迟滞效应。一年内，二者的温度在两个时间点达到相同水平，分别为12月底以及4月底。在4—12月期间，空气温度高于块石，最大温差为1.54 ℃；而在其他时间段，块石温度高于空气温度，最大温差达到0.80 ℃。二者在夏季的温差大于冬季。文献［32-33］中指出，受风速和风向的影响，块石路基以冬半年通风为主的强迫对流与夏半年传导为主的自然对流对多年冻土进行保护，且夏季块石路基顶底板的温差大于冬季。这与模拟得出的结论相同。

通过以上对模拟结果的分析可知，局部非热平衡效应的确存在于块石基底中，其不仅体现在相同时间内块石与空气温度的不相等，也体现二者温度随时间的变化趋势上的差别。两相传热的迟滞效应导致块石的温度变化慢于空气，进而导致了多孔介质中固体与流体的温度不平衡^［20］。

为探求LTNE模型是否可以提高模拟结果的准确性，本节将对开放块石基底路基的模拟过程中同时采用两种多孔介质传热模型进行计算，即LTE与LTNE。在LTNE模型中，使用的换热系数（h_sf ）由Dixon等^［25］在1979年通过理论推导结合试验求得，已有学者在研究加拿大偏远地区岩石废热回收系统时使用该h_sf 公式并取得理想结果^［20］。除此之外，模型的参数保持一致。模拟的结果将与实测值进行对比分析，以证明模拟结果的有效性。

研究表明，块石路基下冻土上限往往在4月达到最高，在10月达到最低^［34］。因此，图7给出了2013年1月、4月、7月、10月共4个月份两个模型的模拟结果中，块石路基以及下部土体温度等值线图，其中蓝色等值线代表LTE模型的温度，红色等值线代表LTNE模型温度。由图可知，10月15日，冻土上限达到一年之内最低，此时块石路基模型的土体温度从上至下降低，且LTE模型的土体温度高于LTNE模型：在此期间，同为-1.5 ℃、-2.0 ℃等值线，LTE模型的等值线位置低于LTNE模型［图7（d）］；而在4月15日，冻土上限达到最高，土体温度从上至下升高，LTE模型中-5 ℃、-2 ℃等值线所覆盖区域大于LTNE模型，在此时LTE模型土体温度整体低于LTNE模型［图7（b）］。由此可知，土层与块石路基之间的能量交换高于LTNE模型，这导致LTE模型中土体温度受环境温度影响比LTNE模型大。

如图8所示，考虑到实测值在2005年12月—2013年6月较为完整，模拟结果也截取相同时间段的土体2 m、5 m以及10 m深度的中央监测点处土体温度随时间变化趋势。在路基建成的前几年内（2004—2007年）为快速降温的阶段，两种模型的模拟结果显示在不同深度处均存在一定程度的降温趋势^［34］。其中，地表以下2 m和5 m处的温度随着季节气温变化而产生明显的周期性波动，而在地表以下10 m处，温度波动幅度受气温变化影响则明显变小。在2005—2006年间，LTE模型在2 m、5 m、10 m三处的年平均地温分别下降了0.73 ℃、0.42 ℃、0.04 ℃；LTNE模型则分别下降了0.85 ℃、0.56 ℃、0.05 ℃；而实测分别下降0.34 ℃、0.28 ℃、0.08 ℃。在此之后，路基的快速降温阶段结束，土体的降温幅度变小，趋于稳定。2007—2011年，LTE模型在三个深度的降温趋势为0.01、0.03、0.04 ℃·a^-1；LTNE模型为0.04、0.07、0.04 ℃·a^-1。除此之外可以发现，LTE模型在各个深度的土体温度的波动幅度都比LTNE模型大。以拥有全年完整实测温度的最后一个年份（2012年）为例，在地表以下2 m处，LTE、LTNE模型的最低温度为 -6.18 ℃、-4.82 ℃，最高温达到-0.87 ℃、-1.26 ℃；在5 m处，二者最低温度达到了-3.46 ℃、-2.88 ℃，最高温为-1.51 ℃、-1.67 ℃；而在10 m处，二者最低温度达到了-2.13 ℃、-1.91 ℃，最高为 -1.99 ℃、-1.82 ℃。在2012全年，LTE模型在三个土层深度的地温波动幅值分别比LTNE模型大1.75 ℃、0.74 ℃以及0.05 ℃。

为更直观地看出两种模型的模拟结果与实测值的贴近程度，图9给出了不同深度处模拟结果与实测数据的月平均温差［图9（a）］以及年误差百分比［图9（b）］。结合图8可以发现，由于环境温度的不稳定，近地表处（2 m深度）的实测地温也存在一定的不稳定性。除此之外，两个模型在该深度温度变化相位存在一定差异，这最终导致了模拟与实测的温差呈现一定的不规律性。在2 m深度，LTE与LTNE两种模型年误差平均值分别为19%以及25%。在5 m处，气温影响以及模拟结果中温度变化相位差逐渐减小，模拟与实测的误差随着月份变化呈现较强的规律性，在每年的1—5月，模拟结果显示该深度地温通常低于实测数据，而在其他月份则相反。由于LTE模型的温度变化幅值更大，因此在相同时间段LTE模型的温差往往大于LTNE模型。此外，由于LTNE模型在快速降温阶段的降温幅度小于LTE模型，导致其2005—2006年的误差百分数大于LTE模型。在2007年之后，LTNE模型误差便逐渐缩小，这说明此时LTNE模型的模拟结果更为准确。二者在该深度年平均误差分别为13%和12%。10 m深度的土体受环境气温影响变得很小，因此模拟结果与实测地温的差值随月份变化的波动也极其微弱。可以发现，两种模型在该深度的误差均随时间逐渐增大，该误差主要因为10 m深度的土体实测地温呈现微弱的上升趋势，而两个模型的模拟结果中的升温趋势并不明显或者呈现微弱的降温趋势。其中，LTE模型在2006—2013年中的月平均地温逐渐低于LTNE模型0.05~0.16 ℃。在该深度，二者年误差分别为7%、15%。

从以上分析可知，无论是LTE模型还是LTNE模型，在路基初建成的几年内，都有着很好的降温效果，虽然后续效果减弱，但仍存在一定的降温效应。这说明块石路基作为广为使用的冻土路基保护措施，其防治效果显而易见。通过模拟结果可以发现，LTE模型的土体温度波动幅度更大，且降温趋势比LTNE模型更明显。原因在于两种模型在数值模拟期间，对块石路基多孔介质内不可渗透固体边界（即块石基底与土体的交界面）计算方式不同：LTE模型中，块石与空气共用一套计算网格，构成的多孔介质将作为一个整体与路基下方土体进行能量交换的数值计算；而LTNE模型中，块石路基域将同时存在两组网格分别对应块石固体骨架以及空气流体，二者以不同的边界条件与下方土体换热。由图5可知，块石温度波动幅度小于空气，所以LTNE模型中土体与路基的能量交换弱于LTE模型，进而导致了二者在模拟结果上的差异。通过与实测温度的对比分析过程中可以发现，LTNE模型的模拟结果与实测值更为接近，其在地下2 m、10 m处的误差比LNE模型平均误差百分数分别小6%以及8%，而在地下5 m处相差1%。这与之前预测的一样，LTNE模型有着更高的准确性，因为它可以更准确地计算多孔介质内的传热。

流固换热系数h_sf 作为LTNE方程中的重要参数，其控制着多孔介质固体骨架与流体之间的能量交换。前人对此展开了丰富的研究。然而由于多孔介质内部结构的不均匀性、各向异性以及研究方法的多样性，h_sf 也存在多种形式。本节将对四种LTNE模型h_sf 的计算公式进行数值模拟，它们均在非微观尺度的多孔介质模型中有着较为广泛的应用，模型的编号如表6所示，各h_sf 的具体公式在表1中给出。

模拟结果显示块石路基中h_sf 的值以年为周期呈现规律性波动，根据其变化趋势，提取2014年2月、5月、8月、11月这4个月份的块石路基中h_sf 的分布云图（图10）。根据云图可以发现块石内部的各个模型h_sf 的变化规律相似，在块石路基的边缘位置，即其与路基填土和地表接触处以及与外界空气交汇处h_sf 的数值高于块石内部的h_sf 数值，这代表在这些区域块石与邻近孔隙空气的能量交换高于其他区域。在一年的大部分时间内，块石路基中的换热处于较为平稳的状态，直到每年的5月，在块石路基中央的底部，h_sf 开始增大，直到8月，在路基内部中央靠左的位置上，迎来一年变化的峰值。从图中还可得知，不同模型的h_sf 数值大小有所区别。在a、b、c、d四个模型的块石内部，h_sf 年平均值分别为0.22、0.35、0.09 W⋅m^-2⋅℃^-1以及0.68 W⋅m^-2⋅℃^-1。

根据表1所列出的h_sf 的公式可知，各个h_sf 虽在形式上有较大差别，但其中均涉及雷诺数（Re）与普朗特数（Pr）。受篇幅所限，不能将每个换热系数中的Re、Pr的变化趋势一一列出。本文以模型c为例，展示了2014年8月15日块石内部的Re以及Pr的分布云图（图11）。可以看出，二者虽在数值大小上有很大差别，其中Re在0~400之间；而Pr在0.85~0.87之间。但二者的分布方式与h_sf 分布方式保持高度一致性。由此可见，块石多孔介质内的换热系数变化与雷诺数（Re）、普朗特数（Pr）的变化相关。

为进一步研究换热系数的变化趋势，图12给出了块石路基内部从左至右三个监测点在2012年7月至2014年7月h_sf 的变化趋势。基于图10可进一步得出，监测点上的h_sf 值随时间产生规律性波动，2012—2014年内波动的范围没有变化，且块石内部三个监测点从左至右波动的峰值依次减小。然而，各个模型的h_sf 波动峰值仍有明显区别。其中，模型b的波动峰值最大，三个监测点从左至右分别为2.89、1.23、0.41 W⋅m^-2⋅℃^-1；而模型c的波动峰值最小，分别为1.16、1.49、0.15 W⋅m^-2⋅℃^-1；模型a和d分别为2.11、0.89、0.28 W⋅m^-2⋅℃^-1以及2.29、1.32、0.68 W⋅m^-2⋅℃^-1。

图13为使用四种不同h_sf 的LTNE模型地下土体在2 m、5 m、10 m深度的温度变化的模拟结果。可以看出，在地下2 m和5 m处，四个模型的模拟结果差异不大，温度曲线基本重合，仅在温度达到每年最低值以及最高值时可以看出差别，且温度最低值处的差异大于最高值处的差异。以2006年为例，2 m处的温度在该年4月达到最低值，分别为 -3.18 ℃、-3.00 ℃、-3.65 ℃以及-2.91 ℃，四个模型在该处总温差为0.74 ℃，平均温差为0.19 ℃；同年5月，5 m处温度达到该年最低值，分别为 -1.97 ℃、-1.89 ℃、-2.00 ℃、-1.84 ℃，总温差0.16 ℃。平均温差0.04 ℃。而在10 m处地温受环境影响变小，温度波动范围也随之减小，四条温度曲线不再发生交叉，在2005—2013年间，a、b、c、d四种模型在该深度的温度平均值分别为-1.69 ℃、 -1.76 ℃、-1.80 ℃以及-1.65 ℃，总温差为0.15 ℃，平均温差为0.04 ℃。除此之外值得注意的是，在这三个不同土体深度中，均是模型c的模拟结果土体温度最低，而模型d温度最高。

结合2.3.1节不难发现，尽管不同模型块石路基内部的h_sf 平均值以及峰值存在一定差异，但是各h_sf 的整体变化趋势十分相近。由于h_sf 峰值存在时间较短（仅在每年8月附近）且范围很小（仅在块石层中心位置），而在其他绝大多数时间中块石基底中的换热系数均处在一个很低的状态，除此之外，计算过程在其他相同模型参数以及相同边界条件的限制下，不同模型的块石与邻近孔隙之间的总体换热差异而导致的土体温差也并不明显，因此不同h_sf 对下部土体温度变化的影响有限。总体而言，路基下部土体的降温幅度随块石夹层内部h_sf 平均值的增大而减小。其中，h_sf 平均值最大的模型d的土体温度最高，而h_sf 平均值最小的模型c的土体温度最低。这可能与h_sf 在块石内部的变化趋势有关。由于h_sf 代表着多孔介质中固体与气体之间的热交换，h_sf 值越大，能量交换也越活跃。7—9月正处夏季，环境温度高，在此时较高的h_sf 使块石路基内部的能量交换更为剧烈，块石从空气中交换的热量较多，从而使块石路基内部温度升高，进而影响土体的降温效果。

本文通过数值模拟的方法，建立了块石基底路基局部热平衡（LTE）模型与局部非热平衡（LTNE）模型，并对LTNE模型中的重要参数——流固换热系数（h_sf ）的几种经典形式进行了数值模拟并将结果对比分析，得出主要结论如下：

（1）受环境温度变化幅度大、块石与空气的有效导热系数差异较大等因素影响，局部非热平衡效应存在于开放块石基底路基中：在同一时间，块石温度与空气温度并不相等，块石温度变化相对于空气存在一定的迟滞效应，且空气的年温差大于块石。空气与块石在夏季和冬季的最大温差可达1.54 ℃和0.80 ℃。

（2）模拟结果显示，LTE模型下伏多年冻土的温度受环境温度变化的影响大于LTNE模型。在与实测温度的对比分析过程中，LTNE模型的模拟结果与实测温度呈现更好的匹配性；在地下2 m和10 m处，LTNE模型误差百分数分别比LTE模型的小6%和8%，在5 m处二者误差百分数相近；LTE模型略高估了块石基底路基对下部土体的降温效果。

（3）LTNE模型中，不同的h_sf 在块石层中的变化趋势相似，且与雷诺数（Re）、普朗特数（Pr）的变化密切相关。但是各个h_sf 的平均值以及峰值仍存在一定的差异。

（4）由于不同h_sf 在块石层中的峰值存在时间短且出现范围很小，因此虽然不同h_sf 平均值和波动峰值存在一些差异，但其对土体温度的影响有限。总体而言，土体温度随着h_sf 平均值的增大而略有升高。这可能与其峰值所处月份（8月）有关，夏季更活跃的换热将使块石与空气的能量交换更高，从而影响块石路基的降温效果。