考慮太陽輻射的寒區襯砌渠道水-熱-力耦合凍脹模型與應用

江浩源，王正中，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］，但與季凍區上建筑物破壞機理有所差異?，F場監測雖可準確評價凍害是否發生，但因凍土區環境惡劣、監測成本高等原因，并不能全尺度的監測出渠道系統內溫度、水分和變形的分布規律及陰陽坡太陽輻射的空間效應和晝夜溫差的時間效應。從提升設計理論水平出發，需要建立季凍區考慮太陽輻射和晝夜溫差的渠道水-熱-力耦合凍脹模型，從理論上分析襯砌渠道的凍害發生原因，在已知外界環境下便可預測渠道基土的水、熱、力耦合變化，對渠道凍脹破壞進行科學分析預測、防控和設計。

目前已有考慮太陽輻射的結構溫度場計算模型，且得到驗證。如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 考慮太陽輻射渠道凍脹模型的熱邊界選取及合理性評價

太陽輻射是渠道凍脹分析中能量的來源，影響著基土內部溫度場、水分場及變形場的耦合作用，是進行寒區渠道設計的重要因素，其中熱邊界合理選取至關重要?，F階段，凍脹數值模擬中常用以下兩種熱邊界：（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）前期設計時，宜采用多年氣溫擬合值和太陽輻射參數確定熱邊界；已建工程凍脹評價時，宜采用日氣溫值和太陽輻射參數確定熱邊界，并輔以監測值對其修正。