考虑太阳辐射的寒区衬砌渠道水-热-力耦合冻胀模型与应用

江浩源，王正中，2，刘铨鸿，王 羿，谭志翔，孟晓栋

（1.西北农林科技大学 旱区寒区水工程安全研究中心 旱区农业水土工程教育部重点实验室，陕西 杨凌 712100；2.西北生态环境资源研究院 冻土工程国家重点实验室，甘肃 兰州 730000）

1 研究背景

我国水资源时空分布不均，南涝北旱、南暖北寒现象普遍，北方寒旱区的灌区及引调水工程建设是保障旱区农业用水及安全饮水的生命线工程［1］。作为主要输水方式的渠道工程在温-水-土-结构等渗-冻互馈机制作用下［2］，常发生鼓胀、断裂、脱空、失稳滑塌等破坏［3-5］，严重影响了渠道的正常输水功能及效益发挥。

太阳辐射是地表热量的主要来源，是影响渠道温度场的重要因素。因渠道走向、结构型式、所处位置等原因，使渠道阴阳坡受到的太阳辐射不同，引起了渠道横向温度场的非对称分布，产生了不稳定的冻融状态［6-8］。经调查，陕西冯家山水库总干退水渠［9］、新疆阜康灌区［10-11］及山东打渔张灌区［12］等部分东西走向渠段均存在显著的阴阳坡效应，阴阳坡冻深及冻胀变形差异明显。太阳辐射是冻融循环的前提基础，是引发渠道水分迁移、冰水相变等水-热-力耦合冻胀破坏的内在驱动力。太阳辐射导致冻结特征的不对称、不同步变化是渠道冻胀破坏的重要原因［7-8］，探明其变化规律对了解渠道冻胀破坏的机理具有重要意义。

季节性冻土区衬砌渠道的阴阳坡效应研究较少，王正中等［7-8］、安元等［6］分别将太阳辐射简化为等效气温边界条件，或采用HOTTEL 辐射模型建立了太阳辐射下渠道的温度场计算模型，而模型中并未考虑渠坡对太阳光的遮蔽作用或基土冻结产生的相变潜热、水分迁移等复杂的冻融循环作用。虽有学者建立了渠道水-热-力耦合模型，但并未考虑渠道在太阳辐射作用下的阴阳坡效应［2-3，13］，或仅是参考监测数据赋予衬砌板温度边界来近似考虑太阳辐射［12］，不能反映太阳辐射的时空变化；同样，多年冻土区的青藏公路、铁路等工程阴阳坡效应虽研究较多［14-15］，但亦存在上述问题。学者主要从现场监测数据入手，结合冻土的地温及路基变形，拟合太阳辐射量与地表温度的计算公式［16］，认为土体中的不均匀热状况导致了差异性的冻融过程和不均匀沉降，使路基产生冻害［17-19］，但与季冻区上建筑物破坏机理有所差异。现场监测虽可准确评价冻害是否发生，但因冻土区环境恶劣、监测成本高等原因，并不能全尺度的监测出渠道系统内温度、水分和变形的分布规律及阴阳坡太阳辐射的空间效应和昼夜温差的时间效应。从提升设计理论水平出发，需要建立季冻区考虑太阳辐射和昼夜温差的渠道水-热-力耦合冻胀模型，从理论上分析衬砌渠道的冻害发生原因，在已知外界环境下便可预测渠道基土的水、热、力耦合变化，对渠道冻胀破坏进行科学分析预测、防控和设计。

目前已有考虑太阳辐射的结构温度场计算模型，且得到验证。如Jin 等［20］基于光线追踪算法，计算了任一时刻拱坝的温度场及应力场分布；Liu 等［21］和Chen 等［22］基于太阳辐射下构件的非均匀温度场试验结果，建立了考虑太阳辐射的温度场计算模型，并对大型穹顶结构［23］和射电望远镜［24］的非均匀温度场进行了分析。上述模型仅需考虑结构吸收的太阳辐射能，而对于太阳辐射作用下的渠道冻胀分析而言，还需考虑温度场、水分场、变形场的耦合作用，模型更为复杂。事实上在太阳辐射作用下的渠基土存在“夜冻昼消”“冬冻春融”现象，冻融循环频繁，具有短时急变、分布不均匀、作用效应复杂等特点，因此建立考虑太阳辐射的渠道水-热-力耦合模型仍存在不少问题有待解决。

综上，本文基于冻土水-热-力耦合理论，考虑冰水相变及水分迁移作用，结合太阳辐射气象参数模型，建立了可反映阴阳坡太阳辐射空间效应和昼夜温差时间效应的渠道水-热-力三场耦合模型，并结合原型监测验证了数值模型的准确性。基于模型，通过系统分析渠道吸收的太阳辐射能、温度场的时空变化规律，量化了渠道断面的水、热、变形差异，分析了寒旱区渠道冻胀破坏的内在原因，并给出渠道冻胀数值计算中热边界的选取建议。该模型从理论上分析渠道冻胀的发展过程，不依赖建成后的现场监测数据，通用性较强，可为寒旱区渠道的设计、维护与研究提供参考。

2 太阳辐射下的渠道冻胀分析模型

季冻区衬砌渠道在空气热对流、太阳辐射、衬砌板间热辐射和环境热辐射等作用下产生复杂的水-热-力耦合冻胀现象。因此在水-热-力耦合模型和太阳辐射模型基础上，提出了渠道阴影计算方法，并基于辐射度算法对上述多种热辐射建立了渠道边界平衡方程，以作为温度场计算的热通量边界条件，从而建立了太阳辐射下的渠道冻胀分析模型。因渠道属线性工程，可按平面问题处理。

2.1 渠基土水-热-力耦合冻胀模型

2.1.1 水-热耦合控制方程 因土颗粒间隙较小而以考虑热传导为主，采用冰水相变修正后的傅里叶热传导方程如下［25-26］：

式中：T 为温度，℃；Lf为冰水相变潜热，kJ/kg；ρi为冰密度，kg/m3；θi为体积冰含量；Cp和λ分别为土体等效定压热容J/（kg·K）和等效导热系数W/（m·K），其值由下述半经验公式估算［2，9，27］：

式中：下标s、w、i、a 分别代表土颗粒、水、冰及气相。

变饱和多孔介质内的水分运动可用水头型Richards 方程描述，添加冰相的方程如下［9，13，27］：

式中：C 为比水容量；h 为基质势，m；k 为土体渗透系数，m/s；i 为重力项。

采用van Genuchten 模型来描述未冻水含量与基质势、渗透系数的关系，方程如下［9，27］：

式中：α、m 为试验拟合参数；θs、θr分别为饱和、残余含水量；Se 为等效饱和度；k（Se）、ks分别为非饱和土、饱和土的渗透系数，m/s。

引入冰阻抗系数I 来近似估算因冰的存在而引起的冻土渗透系数k 的降低，方程如下［25-26］：

因颗粒表面能作用，冻土中的未冻水含量始终与温度保持动态平衡，方程如下［28］：

式中：Wu为未冻水的质量含水量；a、b 为实验参数；T 为含水量为Wu时对应的冻结温度（Tf），℃。

2.1.2 应力-应变控制方程 基于增量弹塑性理论，冻土的应力-应变方程如下［18］：

式中：λp为塑性乘子；Q 为塑性势函数，满足相关联流动法则，等于屈服面函数F，取M-C 准则。

考虑冻土垂直和平行温度梯度的横观各向同性冻胀特征［9］，并依据温度梯度方向实时修正主冻胀方向［2］，坐标转换至整体坐标系，冻胀应变向量方程如下：

式中：ξ为冻胀方向的分配权重，取0.9；m、n 为温度梯度向量的方向余弦；n0为冻土初始孔隙率。

联立式（1）（4）（9）（10）构成渠基冻土的温度场、水分场和变形场的耦合控制方程。由式（1）（4）（9）可计算寒区渠道在外界环境作用下的温度场和冰、水含量分布，结合式（12）得到冻胀应变参数，采用式（10）便可求解得到渠道的应力变形分布，从而实现渠基冻土的水-热-力耦合冻胀分析。

2.2 渠道表面辐射换热模型

2.2.1 太阳辐射模型 太阳辐射具有明显的时空效应，其位置由赤纬δ（°）和时角ω（°）组成的赤道坐标系及由太阳高度角αs和方位角γs组成的地平坐标系决定，各参数计算公式如下［6，29］：

式中：n 为日序数，其中1月1日为1，12月31日为365，表征不同季节日期的参数。

式中：t 为从0 时开始的小时数，该值表征当日时间的参数，上午为负，下午为正。

式中：φ 为地理纬度；γs偏东为负，偏西为正。

对北方寒旱区而言，Hottel 晴空模型因具有较好的适用性和通用性而被优先采用［6，21，29］，可弥补国内站点稀疏、数据较少等不足。该模型中太阳射线到达大气层外切平面的太阳辐射强度G0由下式计算：

式中：Gsc为太阳常数，表示单位时间单位面积上所接收的太阳辐照度，取1367 W/m2。

任一斜面的太阳入射角i、太阳辐照度Gs（太阳直射辐照度Gtb和散射辐照度Gtd）公式如下［29-30］：

式中：β 为斜面倾角；γt为斜面方位角，指渠道表面法线在水平面上的投影与正南方向的夹角，面向东时为负，面向西时为正；斜面上的日照开始和结束时间以i=90°来计算［30］；α0、α1、k 为标准晴空大气常数，根据文献［6，29］计算。

2.2.2 渠道阴影计算方法 渠道边坡相互遮挡产生的阴影分布影响了渠道表面的太阳辐射分布。在太阳高度角较小时，阴坡的遮挡使渠道局部形成阴影，且此阴影的边界线与渠道的走向平行。基于此，可将渠道阴坡看作一根杆，将此杆在太阳下形成的阴影边界点做渠道走向的平行线，此为渠道阴坡产生的阴影边界。以北方东西走向渠道为例，为增加图幅立体感，将方向逆时针旋转，计算简图见图1。

图1 不同走向渠道阴影长度计算示意

图中OA 为细杆，高度为渠深h，l1和l4分别为渠坡顶在水平面上的投影。若此时太阳位于渠道正南方，细杆将产生OB 的阴影，渠道阴影区域为l2—l4，l1—l2为受光区；若将渠道顺时针旋转一定角度，此值为面向太阳的渠道表面方位角γt。为便于图幅描述，亦可表示为渠道走向不变，太阳逆时针旋转此角度。此时细杆OA 将产生OC 的阴影，渠道阴影区域为l3—l4，阴影边界离细杆距离为OD。因太阳高度角未发生变化，细杆产生的阴影长度一致，即OB=OC。综合考虑渠道走向和太阳位置变化，其产生的阴影长度OD 由下式计算：

进一步结合渠道断面进行阴影判定，如图2所示。A、E 均为太阳光线可穿过位置，为简化图幅，仅以A 点为例介绍。D1、D2、D3为不同太阳位置下产生的阴影点，结合阴影长度OD1、OD2、OD3和渠道断面几何关系，确定不同时刻下渠道的受光区R，计算公式如下：

图2 渠道阴影计算示意

联立式（13）—式（18）及（21）（22）得到不同地区、时间下的不同走向、断面大小的渠道表面太阳入射角和受光面积，二者共同决定渠道表面的日照时间，结合式（19）（20）得到各表面接受的太阳辐照度。2.2.3 衬砌板间辐射和环境热辐射 因衬砌板间温度及其与环境温度不同而产生热辐射，计算简图见图3，任一点P 得到的S ′面的辐照度（Gm：W/m2）及周围环境辐照度（Gamb：W/m2）见下式［31］：

式中：J ′、Jamb分别为S′面和环境（Samb）的辐射度，W/m2；n 和n′为外法线矢量；r 为面上两点的距离矢量；n 为折射率，不透明物体取1；σ为Stefan-Boltzmann 常数，取5.67×10-8W/（m2·K4）；Famb为环境角系数；Tamb为环境温度，K。2.2.4 基于辐射度方法的表面辐射换热方程 辐射度J 由自身辐射及对辐照度（Gm、Gs、Gamb）的反射辐射组成，而每一点Gm又是其他可见点J 的函数，据此得到如下辐射平衡方程［31］：

式中：ρd为漫反射系数；ε 为发射率；T 为渠道表面温度，K。

假设渠道表面为理想漫射灰体，发射率与吸收率α 相等，可计算得到渠道表面吸收的辐射量，并以此作为冻土水-热-力耦合中热模块的第二边界条件（热通量边界条件），方程如下：

图3 角系数计算示意

3 数值模型

以新疆阜康某灌区（87oE、44oN）E-W 走向“阴阳坡效应”显著的大型输水渠道为例，基于现场实测数据，联立上节考虑太阳辐射及阴影遮蔽的渠道水-热-力耦合方程，并采用COMSOL 求解因太阳辐射而导致的不均匀水、热、变形分布，分析渠道冻害发生原因，并给出热边界的选取建议。

3.1 基本概况该梯形渠道修建于季节性冻土区，属温带大陆性干旱气候，多年平均温度6.5 ℃，极端最低温度-25～-32.2 ℃，多年最大冻深为1.72 m，平均冻深为1.45 m。渠基土中粒径小于0.075 mm的土粒质量占比大于10%，为冻胀敏感性土。

渠道设计流量为21 m3/s，加大流量26 m3/s，渠深3.0 m，正常水深2.6 m，底宽4.0 m，坡比1∶1.5。采用C20 现浇混凝土衬砌，板厚12 cm。姜海波等［10-11］于2013年11月5日—2014年3月28日对该渠道的温度场（WS-4 型遥感土壤温度计，测量深度距基土表面5～120 cm）、冻深（DTM-2 型冻土器）、水分场（土钻取样后烘干，测量距基土表面5～120 cm）及衬砌板冻胀量（水准仪）等进行监测，监测点布设及渠道有限元网格如图4所示。冻深区域内网络高度最小为0.5 cm，最大为1 cm，最大时间步长120 s。其中A、B、C 为阳坡测点，分别位于1、2/3、1/3 设计水位处，D 为底板中心，E、F、G 为阴坡测点，位置同阳坡对称。

图4 有限元网格及数值模型示意（单位：m）

3.2 计算参数根据试验段土工试验［10-11］可知，浅层基土天然体积含水率为0.256～0.368，处于潮湿-饱和状态，液性指数为0.53～0.82，塑性指数为5.8～6.4，属低液限粉土。为简化计算，计算参数近似取试验均值，其中黏聚力为10.0 kPa，内摩擦角为20°，干密度为1.6 g/cm3；结合基土土质，其弹性参数与温度的关系采用式（27）（28）计算［2，13，18］；土体各相、漫反射系数和发射率、冻结曲线及土水特征曲线的参数见表1［2，9，24，28，32］，混凝土参数见表2［2，20］。

表1 土体材料计算参数

表2 混凝土材料计算参数

3.3 边界和初始条件渠道上表面温度边界包括热辐射（式（26））及热对流边界，其中对流换热方程如下：

式中：n 为渠道上表面法向向量；Tamb和T 分别为外界温度和渠道表面温度，℃；hc为对流换热系数，W/（m2·℃），由风速来近似确定［2］，公式如下：

图5 外界最高、最低温度变化曲线

根据现场实测数据，渠顶和渠底风速近似取5 m/s 和1 m/s，并采用二次抛物线函数过渡。

该渠道已修建多年，基土温度基本达到吞吐平衡，因此以年周期平衡温度场作为计算的初始温度场（11月1日）。根据气象站温度监测值，采用正弦函数拟合温度变化曲线［9，19］，公式如下：

式中：T0为年平均温度，取6.5 ℃；B 为年温度振幅，取22.5 ℃；tm为时间，月；α0为相角，取π。

选取试验段渠道岸边百叶箱每日温度实测极值［10］，整理如图5所示。日温度变化曲线结合每日极值温度采用正弦函数拟合［6］，以此作为计算分析的外界温度，公式如下：

渠基下表面恒温层温度近似取年平均温度6.5 ℃，距地表深度为11 m［33］。

渠道因防渗措施不当，渗漏较为严重，停水后，渠底下4.6 m 深处土体基本处于饱和状态，水分场基本稳定。此后按设计水位模拟通水一个月，并在冻结试验开始前一周排干水分，以此作为计算分析的初始水分场。

因复合土工膜较薄（<3 mm），采用无厚度土工膜渗阻计算，土工膜前区域为混凝土衬砌，膜后区域为渠基土，质量守恒方程建立如下：

式中：下标b、u、d 分别代表土工膜、膜前和膜后区域；n 为边界外法线矢量；k 为渗透系数，混凝土取1.2×10-9m/s，土工膜取1.0×10-11m/s；h 为压力水头，m；▽h为渗流域在边界处的压力梯度；R 为渗流阻抗，s；d 为厚度，m。

在渠道周边施加法向位移约束，见图4，满足初始地应力平衡。

4 结果分析

4.1 温度场分析

4.1.1 渠道阴影分布、太阳辐射及温度变化 因渠道为二维断面，采用衬砌板吸收的太阳辐射量沿渠周分布来定量描述不同时刻渠道的阴影分布，并分析由此而导致的温度差异。因东西走向渠道上、下午接收的太阳辐射量基本一致，且阜康地区比北京约晚2 h，因此取10∶00—14∶00 时段分析太阳辐射量；因衬砌板温度存在累积作用，取10∶00—18∶00 时段分析。以12月10日为例，结果如图6—8所示。其中，坐标原点从阳坡渠顶算起，下同。

由图6 可知，从早上太阳初升，随时间增加，太阳直射区域由阳坡逐渐扩大至渠底，阴影区域逐渐减少；随太阳高度角的增加，太阳辐射量逐渐增大，至14∶00 时最大达374 W/m2，照射近1/2 渠周；而阴坡始终处于阴影区域，仅吸收太阳散射值，太阳辐射远小于阳坡，辐射量差值达324.6 W/m2。

图6 不同时刻衬砌吸收的太阳辐射沿渠周分布曲线

由图7 可知，因太阳辐射区域及量值的增加，衬砌表面温度随之升高，阳坡温度远高于阴坡，最大差值为7.5 ℃，“阴阳坡效应”明显，且由太阳辐射引起的温度累积效应明显；同时衬砌板阴、阳坡温度均随外界气温的降低而降低，且二者差值逐渐减小。

图7 不同时刻衬砌表面温度沿渠周分布曲线

由图8 可进一步得知，衬砌表面温度与外界气温和太阳辐射存在较好的相关性，但滞后明显，二者共同造成了阴阳坡的横向温度差异。太阳辐射时段内，阳坡温度主要受太阳辐射影响，而阴坡则主要受外界气温热对流作用；非太阳辐射时段内，阴阳坡温度主要受外界气温影响。因此以往分析中不考虑阴阳坡或将太阳辐射包含在外界气温中或取差异性的渠道表面温度值，均不能科学合理的反映出阴阳坡温度场的时空分布差异。

图8 阴阳坡2/3 设计水位点温度、辐射分布曲线

4.1.2 渠道阴阳坡衬砌板表面温度变化 取监测期内阴阳坡衬砌板最低和日均温度值来分析二者的横向温度差异，如图9所示。

图9 阴阳坡衬砌板表面温度变化曲线

由现场实测值可知，衬砌板温度随外界气温的降低而降低，其变化规律与外界气温基本一致，但存在滞后性。阴、阳坡板最低温度分别为-24.6 和-21.5 ℃，温差最大为4.3 ℃。

由数值结果可知，最低温度值略高于现场实测值，可能是由于太阳辐射模型中未考虑空气透明度、阴天、降雪等因素引起的太阳辐射值偏高所致。同时温度变化规律同现场实测基本一致，且阴阳坡板日均温度差异明显，最大为3.5 ℃，进一步表明本模型可较好的分析由太阳辐射的时空效应和阴坡遮蔽作用而引起的衬砌板横向温度差异。

4.1.3 渠道阴阳坡冻深及基土温度场变化 渠道表面温度场的不对称导致基土温度场的差异性分布。以阴、阳坡2/3 设计水位点为例，将0 ℃等温线作为冻深判别依据，其发展如图10所示。

图10 阴、阳坡冻深发展过程线

因昼夜温差及太阳辐射作用，冻结前期阴、阳坡呈现出“夜冻昼消”的现象。由图可知，冻深随负积温的累积而逐渐向下发展，“夜冻昼消”现象逐渐消失，其中阴坡于11月25日消失，阳坡于12月10日消失，二者前后差异15 d。结合图5，至1月10日最低温时刻，冻深并未发展至最大，后期虽存在大幅度升温现象，但并未达到冻土的融化温度，冻深依然向下发展；至3月16日，冻深最大，为135 cm；随着温度的继续升高，阳坡于3月17日开始融化，阴坡于3月18日开始融化，二者前后差异1 d。外加太阳辐射作用，使得基土的融化速率远大于冻结速率，至监测期的最后一天3月28日，阴、阳坡基土表面分别融化40 和49 cm，冻深分别为97 和66 cm，差异明显。

为进一步分析太阳辐射作用下的基土温度场分布规律，取最大冻深时刻，如图11所示。

图11 3月16日渠道温度场分布曲线

由图可知，阴坡温度明显低于阳坡，温度场呈不对称分布；基土温度随外界气温（-8 ℃）变化存在滞后性，冻深以上温度变化剧烈，以下则较为缓慢。结合冻深监测结果可知，数值结果略小，最大差异为22 cm，主要原因在于太阳辐射模型未考虑阴天、降水对太阳光的阻挡所致。但反映出的冻深分布规律与监测结果分布规律基本一致，仅在数值上存在微小差异，表明本文模型可较好的模拟太阳辐射作用下的基土温度场分布规律。

4.2 水分场分析太阳辐射和环境温度共同作用引起渠道的不对称温度场，势必引起不对称的水分场分布。为此取11月5日（初始时刻）和1月20日（监测时刻）阳坡、渠底及阴坡的测点A、B、D、F、G 的总含水率沿基土深度分布结果，如图12所示。

图12 关键点渠基土体积含水率分布曲线

由图可知，相比于初始时刻，冻结区域内因外界负积温作用而使未冻区内的水分逐渐向上迁移，水分集聚现象明显。由数值结果可知，因太阳辐射和昼夜温度的叠加作用，冻土内存在多个水分集聚位置，而非只发生在冻结锋面处；阳坡、渠底及阴坡水分的集聚个数、位置及深度均不同，呈现出阴坡含水量最多而阳坡最少的分布规律。对比监测结果可知，因现场选取局部点来烘干，仅能反映出水分运移规律，并不能严格反映出含水量最高点或冻结锋面位置，但在冻深范围内呈现阴坡、渠底总含水量高而阳坡较少的分布规律，与数值结果一致。然数值解的表层含水量偏高，80～100 cm 深度处偏低，可能与未考虑阴天、蒸发等因素及现场局部地质引起的初始含水偏高所致；但其他位置点的含水量吻合度较高，满足工程要求，即此模型可模拟太阳辐射下的水分迁移、冻结相变等特征。

为进一步分析不对称的冰、水含量分布，取最大冻深时刻（3月16日）和最大融深时刻（3月28日）的水分场结果，如图13—14所示。

由图13 可知，冻结区域内的未冻水含量随温度降低而逐渐减少，最低为0.06，接近于基土的残余含水量；而相应位置的结冰量逐渐增加，冰含量呈条带状分布，无明显的空间分布规律［2-3，12］，但与基土冻深发展相呼应（图10）；因冻结初期基土表层常处于冻融循环状态，在多次温度抽吸力下表层冰含量较大；冰、水总含量最大为0.6，已远大于基土的孔隙率，将发生严重的冻胀变形。

3月16日表层基土出现局部融化现象，含水量增加而冰含量减少。随着外界温度的逐渐升高，冻结层内冰透镜体开始从两端发生双向融化，局部区域形成过量孔隙水，随着时间推移，孔隙水逐渐扩散，形成随深度增加含水量逐渐增大的分布规律（图14（a）），因冻结区域内冰含量的阻水作用，使含水量较初始时刻明显增大。对比图13（b）和图14（b）可知，冻深范围外冰透镜体全部融化，且冻深线附近因温度升高至其融化温度，使冰含量值明显减少，而冻深范围内的其余区域因热流传输密度少和冻深线附近的冰融化吸热，从而使冰含量未发生较大变化。

图13 3月16日渠道水分场分布

图14 3月28日渠道水分场分布

由图13、图14 亦可发现，因太阳辐射和阴坡遮蔽而产生的热边界差异产生了明显的不对称水分场分布，阴坡总含水量明显大于阳坡。首先，阴阳坡的冻深差异特征形成了阴坡冰透镜体总体厚度较阳坡大的现象；其次，阴阳坡的热边界差异使阴坡含冰量大于阳坡含冰量；最后，因阳面接受太阳辐射多，温升快，从而使阳坡融化速度明显大于阴坡，冰含量快速减少。

4.3 变形场分析 由上节可知，阴、阳坡局部位置总含水量大于孔隙率，将会产生差异性冻胀变形，选取多个时刻下衬砌板的法向冻胀变形沿渠周分布曲线，如图15所示。

图15 衬砌板法向冻胀变形沿渠周分布曲线

由图可知，梯形渠道衬砌因基土冻胀作用，呈现出渠底向上隆起，坡脚受挤压明显，两坡板向渠内凸起，衬砌整体上抬的变形趋势。其中，渠底衬砌板在偏阳坡处的法向冻胀量最大，阴、阳坡板约在1/4～1/3 坡长位置处最大。在因太阳辐射而产生的差异冻胀力作用下，衬砌冻胀变形量由阴坡、渠底至阳坡依次递减，阴、阳坡板冻胀变形差异明显，整体与基土内冻深发展、水分分布规律一致，即渠道水热分布的不对称、不同步特征导致了阴、阳坡板不对称的变形分布。

由上图亦可知，随着冻结时间的发展，衬砌板法向冻胀变形均逐渐增大，至3月16日时，冻胀量最大，阳坡、渠底和阴坡最大分别为8.47、10.97 和12.48 cm，与现场监测值基本吻合，且衬砌板的整体变形趋势与监测值分布规律基本一致，同时，在后续融化过程中，衬砌板冻胀变形亦随基土内冰水总含量的减少而逐渐变小，表明本文的数值模型可较好的模拟衬砌结构在基土冻胀作用下的冻融变形特征。

5 考虑太阳辐射渠道冻胀模型的热边界选取及合理性评价

太阳辐射是渠道冻胀分析中能量的来源，影响着基土内部温度场、水分场及变形场的耦合作用，是进行寒区渠道设计的重要因素，其中热边界合理选取至关重要。现阶段，冻胀数值模拟中常用以下两种热边界：（1）以日或月为单位拟合环境温度，只采用热对流边界［2-3，9，13］，但无法考虑太阳辐射产生的阴阳坡效应；（2）结合现场实测数据，采用正弦函数拟合阴、阳坡温度变化曲线，并将其直接赋予结构边界［17-18］，但却无法充分反映出太阳辐射的空间效应和昼夜温差的时间效应。目前公认更为准确的是第二种以现场地温监测为依据的拟合温度边界，以此为参照，对比分析本文采用的太阳辐射、面-面辐射、环境辐射的热辐射边界和昼夜温差的热对流边界组合，对热边界选取及其合理性进行评价。拟合温度边界取上章仿真结果的阴坡、渠底和阳坡的日均表面温度数据（图9）来正弦拟合热边界，可弥补现场监测点少且数据不全的缺点。下面将从温度、水分、变形来分别阐述。

取阴、阳坡2/3 设计水位点的冻深发展过程结果，如图16所示。由图可知，拟合的温度边界使衬砌板温度缓慢变化，冻结速率发展较为均匀；而本文模型因考虑了太阳辐射空间效应和昼夜温度时间效应的耦合叠加作用，使得冻深发展在局部位置出现了短暂的平台状，冻结速率忽快忽慢，但与外界的环境变化相协调，更符合现场实际规律。

图16 不同热边界下的冻深发展过程线

阴、阳坡含水量分布规律基本类似，取最大冻深时刻下阴坡2/3 设计水位点的总含水量沿基土深度的分布曲线，如图17所示。由图可知，因拟合温度边界下冻结速率变化较为均匀，除浅层基土外，含水率由上至下逐渐增大，至冻结峰面处最大，与单向冻结结果基本一致。而本文模型因考虑太阳辐射和昼夜温差作用，冻结初期的“夜冻昼消”往复冻融循环作用使得基土表层含水增加；水分分布与冻土的冻结速率相对应，冻结速率慢时水分迁移多，反之亦然，更能反映出外界环境作用引起的水热耦合作用。

图17 不同热边界下的水分场分布

总含水率和冻深差异导致了不同的冻胀变形，取阴、阳坡2/3 设计水位点，其衬砌板法向冻胀变形随时间的分布规律如图18所示。由图可知，衬砌板的起始冻胀时刻、冻胀量发展速度和量级，以及融化速度等变量均与热边界有关，拟合温度边界方法得到的冻胀量不能反映出太阳辐射和昼夜温差作用所引起的冻结速率差异而导致的水热分布规律，导致结构冻融破坏计算不准确。

图18 不同工况下衬砌板法向冻胀量随时间变化曲线

结合温度、水分和变形结果可知，本文数值模型的结果与目前公认的以地温监测数据为基础的拟合温度边界所得规律基本一致，结合上一节结果，可进一步说明本文模型的准确性；同时本文模型综合考虑了太阳辐射的空间效应和昼夜温度的时间效应，可反映出外界环境作用下寒区渠道水、热、力耦合的冻胀破坏本质，特别是能减少对现场监测值的依赖性，且具有理论通用性，相比于目前使用较多的拟合温度边界，结果更加合理、准确，且可辅助设计单位进行前期渠道设计，满足精细化设计、复核的要求。

6 结论

考虑太阳辐射和昼夜温差的时空变化，提出了渠道的阴影计算方法，基于冻土的横观各向同性冻胀特征，建立了“太阳辐射-环境温度-冻土-工程”相互作用的水-热-力耦合冻胀模型，并探讨了不同热边界对渠道水-热-力耦合计算的影响，结合工程监测及数值模拟结果，得到如下结论：（1）考虑太阳辐射的渠道冻胀模型计算的温度、水分和变形结果与现场实测值较为吻合，分布规律一致，表明本文模型可分析太阳辐射空间效应和昼夜温差时间效应下的寒区渠道冻结过程，理论通用性强，且可减少对现场监测的依赖，可为寒区渠道的冻胀模拟、预测及灾害防控设计提供理论依据。（2）探明了渠道表面太阳辐射的时空分布规律，阴、阳坡太阳辐射时长、辐射量存在显著差异，其表面温度由外界气温和太阳辐射量共同决定，计算的阴、阳坡日均温差最大为3.5 ℃，导致了阴坡早于阳坡15 d 冻结，晚于阳坡1 d 融化，冻深发展差异明显，最大相差31 cm，渠道阴阳坡温度场分布不对称且不同步。（3）冻结期渠道阴坡结冰区域及量值明显大于阳坡，“夜冻昼消”的冻融循环作用使基土浅层含水量明显增多，而融化期阳坡冰融化速度大于阴坡，水分场呈不对称、不同步、层状分布等特点。（4）渠坡板和底板的最大冻胀变形分别发生在1/4～1/3 坡长位置和偏阳坡位置，阴、阳坡变形随时间的发展过程呈不均匀、不同步、不对称性，计算的最大变形分别为12.48 和8.47 cm。（5）渠道阴阳坡效应产生的本质在于太阳辐射和昼夜温差引起的渠道能量收支差异、温度梯度和冻结速率变化所致，即热边界引起的渠道基土的水-热-力动态耦合作用。（6）前期设计时，宜采用多年气温拟合值和太阳辐射参数确定热边界；已建工程冻胀评价时，宜采用日气温值和太阳辐射参数确定热边界，并辅以监测值对其修正。