江浩源，王正中，2，劉銓鴻，王 羿，譚志翔，孟曉棟

（1.西北農(nóng)林科技大學(xué) 旱區(qū)寒區(qū)水工程安全研究中心 旱區(qū)農(nóng)業(yè)水土工程教育部重點(diǎn)實(shí)驗(yàn)室，陜西 楊凌 712100；2.西北生態(tài)環(huán)境資源研究院 凍土工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，甘肅 蘭州 730000）

1 研究背景

我國(guó)水資源時(shí)空分布不均，南澇北旱、南暖北寒現(xiàn)象普遍，北方寒旱區(qū)的灌區(qū)及引調(diào)水工程建設(shè)是保障旱區(qū)農(nóng)業(yè)用水及安全飲水的生命線工程［1］。作為主要輸水方式的渠道工程在溫-水-土-結(jié)構(gòu)等滲-凍互饋機(jī)制作用下［2］，常發(fā)生鼓脹、斷裂、脫空、失穩(wěn)滑塌等破壞［3-5］，嚴(yán)重影響了渠道的正常輸水功能及效益發(fā)揮。

太陽(yáng)輻射是地表熱量的主要來(lái)源，是影響渠道溫度場(chǎng)的重要因素。因渠道走向、結(jié)構(gòu)型式、所處位置等原因，使渠道陰陽(yáng)坡受到的太陽(yáng)輻射不同，引起了渠道橫向溫度場(chǎng)的非對(duì)稱(chēng)分布，產(chǎn)生了不穩(wěn)定的凍融狀態(tài)［6-8］。經(jīng)調(diào)查，陜西馮家山水庫(kù)總干退水渠［9］、新疆阜康灌區(qū)［10-11］及山東打漁張灌區(qū)［12］等部分東西走向渠段均存在顯著的陰陽(yáng)坡效應(yīng)，陰陽(yáng)坡凍深及凍脹變形差異明顯。太陽(yáng)輻射是凍融循環(huán)的前提基礎(chǔ)，是引發(fā)渠道水分遷移、冰水相變等水-熱-力耦合凍脹破壞的內(nèi)在驅(qū)動(dòng)力。太陽(yáng)輻射導(dǎo)致凍結(jié)特征的不對(duì)稱(chēng)、不同步變化是渠道凍脹破壞的重要原因［7-8］，探明其變化規(guī)律對(duì)了解渠道凍脹破壞的機(jī)理具有重要意義。

季節(jié)性?xún)鐾羺^(qū)襯砌渠道的陰陽(yáng)坡效應(yīng)研究較少，王正中等［7-8］、安元等［6］分別將太陽(yáng)輻射簡(jiǎn)化為等效氣溫邊界條件，或采用HOTTEL 輻射模型建立了太陽(yáng)輻射下渠道的溫度場(chǎng)計(jì)算模型，而模型中并未考慮渠坡對(duì)太陽(yáng)光的遮蔽作用或基土凍結(jié)產(chǎn)生的相變潛熱、水分遷移等復(fù)雜的凍融循環(huán)作用。雖有學(xué)者建立了渠道水-熱-力耦合模型，但并未考慮渠道在太陽(yáng)輻射作用下的陰陽(yáng)坡效應(yīng)［2-3，13］，或僅是參考監(jiān)測(cè)數(shù)據(jù)賦予襯砌板溫度邊界來(lái)近似考慮太陽(yáng)輻射［12］，不能反映太陽(yáng)輻射的時(shí)空變化；同樣，多年凍土區(qū)的青藏公路、鐵路等工程陰陽(yáng)坡效應(yīng)雖研究較多［14-15］，但亦存在上述問(wèn)題。學(xué)者主要從現(xiàn)場(chǎng)監(jiān)測(cè)數(shù)據(jù)入手，結(jié)合凍土的地溫及路基變形，擬合太陽(yáng)輻射量與地表溫度的計(jì)算公式［16］，認(rèn)為土體中的不均勻熱狀況導(dǎo)致了差異性的凍融過(guò)程和不均勻沉降，使路基產(chǎn)生凍害［17-19］，但與季凍區(qū)上建筑物破壞機(jī)理有所差異。現(xiàn)場(chǎng)監(jiān)測(cè)雖可準(zhǔn)確評(píng)價(jià)凍害是否發(fā)生，但因凍土區(qū)環(huán)境惡劣、監(jiān)測(cè)成本高等原因，并不能全尺度的監(jiān)測(cè)出渠道系統(tǒng)內(nèi)溫度、水分和變形的分布規(guī)律及陰陽(yáng)坡太陽(yáng)輻射的空間效應(yīng)和晝夜溫差的時(shí)間效應(yīng)。從提升設(shè)計(jì)理論水平出發(fā)，需要建立季凍區(qū)考慮太陽(yáng)輻射和晝夜溫差的渠道水-熱-力耦合凍脹模型，從理論上分析襯砌渠道的凍害發(fā)生原因，在已知外界環(huán)境下便可預(yù)測(cè)渠道基土的水、熱、力耦合變化，對(duì)渠道凍脹破壞進(jìn)行科學(xué)分析預(yù)測(cè)、防控和設(shè)計(jì)。

目前已有考慮太陽(yáng)輻射的結(jié)構(gòu)溫度場(chǎng)計(jì)算模型，且得到驗(yàn)證。如Jin 等［20］基于光線追蹤算法，計(jì)算了任一時(shí)刻拱壩的溫度場(chǎng)及應(yīng)力場(chǎng)分布；Liu 等［21］和Chen 等［22］基于太陽(yáng)輻射下構(gòu)件的非均勻溫度場(chǎng)試驗(yàn)結(jié)果，建立了考慮太陽(yáng)輻射的溫度場(chǎng)計(jì)算模型，并對(duì)大型穹頂結(jié)構(gòu)［23］和射電望遠(yuǎn)鏡［24］的非均勻溫度場(chǎng)進(jìn)行了分析。上述模型僅需考慮結(jié)構(gòu)吸收的太陽(yáng)輻射能，而對(duì)于太陽(yáng)輻射作用下的渠道凍脹分析而言，還需考慮溫度場(chǎng)、水分場(chǎng)、變形場(chǎng)的耦合作用，模型更為復(fù)雜。事實(shí)上在太陽(yáng)輻射作用下的渠基土存在“夜凍晝消”“冬凍春融”現(xiàn)象，凍融循環(huán)頻繁，具有短時(shí)急變、分布不均勻、作用效應(yīng)復(fù)雜等特點(diǎn)，因此建立考慮太陽(yáng)輻射的渠道水-熱-力耦合模型仍存在不少問(wèn)題有待解決。

綜上，本文基于凍土水-熱-力耦合理論，考慮冰水相變及水分遷移作用，結(jié)合太陽(yáng)輻射氣象參數(shù)模型，建立了可反映陰陽(yáng)坡太陽(yáng)輻射空間效應(yīng)和晝夜溫差時(shí)間效應(yīng)的渠道水-熱-力三場(chǎng)耦合模型，并結(jié)合原型監(jiān)測(cè)驗(yàn)證了數(shù)值模型的準(zhǔn)確性?；谀Ｐ?，通過(guò)系統(tǒng)分析渠道吸收的太陽(yáng)輻射能、溫度場(chǎng)的時(shí)空變化規(guī)律，量化了渠道斷面的水、熱、變形差異，分析了寒旱區(qū)渠道凍脹破壞的內(nèi)在原因，并給出渠道凍脹數(shù)值計(jì)算中熱邊界的選取建議。該模型從理論上分析渠道凍脹的發(fā)展過(guò)程，不依賴(lài)建成后的現(xiàn)場(chǎng)監(jiān)測(cè)數(shù)據(jù)，通用性較強(qiáng)，可為寒旱區(qū)渠道的設(shè)計(jì)、維護(hù)與研究提供參考。

2 太陽(yáng)輻射下的渠道凍脹分析模型

季凍區(qū)襯砌渠道在空氣熱對(duì)流、太陽(yáng)輻射、襯砌板間熱輻射和環(huán)境熱輻射等作用下產(chǎn)生復(fù)雜的水-熱-力耦合凍脹現(xiàn)象。因此在水-熱-力耦合模型和太陽(yáng)輻射模型基礎(chǔ)上，提出了渠道陰影計(jì)算方法，并基于輻射度算法對(duì)上述多種熱輻射建立了渠道邊界平衡方程，以作為溫度場(chǎng)計(jì)算的熱通量邊界條件，從而建立了太陽(yáng)輻射下的渠道凍脹分析模型。因渠道屬線性工程，可按平面問(wèn)題處理。

2.1 渠基土水-熱-力耦合凍脹模型

2.1.1 水-熱耦合控制方程 因土顆粒間隙較小而以考慮熱傳導(dǎo)為主，采用冰水相變修正后的傅里葉熱傳導(dǎo)方程如下［25-26］：

式中：T 為溫度，℃；Lf為冰水相變潛熱，kJ/kg；ρi為冰密度，kg/m3；θi為體積冰含量；Cp和λ分別為土體等效定壓熱容J/（kg·K）和等效導(dǎo)熱系數(shù)W/（m·K），其值由下述半經(jīng)驗(yàn)公式估算［2，9，27］：

式中：下標(biāo)s、w、i、a 分別代表土顆粒、水、冰及氣相。

變飽和多孔介質(zhì)內(nèi)的水分運(yùn)動(dòng)可用水頭型Richards 方程描述，添加冰相的方程如下［9，13，27］：

式中：C 為比水容量；h 為基質(zhì)勢(shì)，m；k 為土體滲透系數(shù)，m/s；i 為重力項(xiàng)。

采用van Genuchten 模型來(lái)描述未凍水含量與基質(zhì)勢(shì)、滲透系數(shù)的關(guān)系，方程如下［9，27］：

式中：α、m 為試驗(yàn)擬合參數(shù)；θs、θr分別為飽和、殘余含水量；Se 為等效飽和度；k（Se）、ks分別為非飽和土、飽和土的滲透系數(shù)，m/s。

引入冰阻抗系數(shù)I 來(lái)近似估算因冰的存在而引起的凍土滲透系數(shù)k 的降低，方程如下［25-26］：

因顆粒表面能作用，凍土中的未凍水含量始終與溫度保持動(dòng)態(tài)平衡，方程如下［28］：

式中：Wu為未凍水的質(zhì)量含水量；a、b 為實(shí)驗(yàn)參數(shù)；T 為含水量為Wu時(shí)對(duì)應(yīng)的凍結(jié)溫度（Tf），℃。

2.1.2 應(yīng)力-應(yīng)變控制方程 基于增量彈塑性理論，凍土的應(yīng)力-應(yīng)變方程如下［18］：

式中：λp為塑性乘子；Q 為塑性勢(shì)函數(shù)，滿足相關(guān)聯(lián)流動(dòng)法則，等于屈服面函數(shù)F，取M-C 準(zhǔn)則。

考慮凍土垂直和平行溫度梯度的橫觀各向同性?xún)雒浱卣鳎?］，并依據(jù)溫度梯度方向?qū)崟r(shí)修正主凍脹方向［2］，坐標(biāo)轉(zhuǎn)換至整體坐標(biāo)系，凍脹應(yīng)變向量方程如下：

式中：ξ為凍脹方向的分配權(quán)重，取0.9；m、n 為溫度梯度向量的方向余弦；n0為凍土初始孔隙率。

聯(lián)立式（1）（4）（9）（10）構(gòu)成渠基凍土的溫度場(chǎng)、水分場(chǎng)和變形場(chǎng)的耦合控制方程。由式（1）（4）（9）可計(jì)算寒區(qū)渠道在外界環(huán)境作用下的溫度場(chǎng)和冰、水含量分布，結(jié)合式（12）得到凍脹應(yīng)變參數(shù)，采用式（10）便可求解得到渠道的應(yīng)力變形分布，從而實(shí)現(xiàn)渠基凍土的水-熱-力耦合凍脹分析。

2.2 渠道表面輻射換熱模型

2.2.1 太陽(yáng)輻射模型 太陽(yáng)輻射具有明顯的時(shí)空效應(yīng)，其位置由赤緯δ（°）和時(shí)角ω（°）組成的赤道坐標(biāo)系及由太陽(yáng)高度角αs和方位角γs組成的地平坐標(biāo)系決定，各參數(shù)計(jì)算公式如下［6，29］：

式中：n 為日序數(shù)，其中1月1日為1，12月31日為365，表征不同季節(jié)日期的參數(shù)。

式中：t 為從0 時(shí)開(kāi)始的小時(shí)數(shù)，該值表征當(dāng)日時(shí)間的參數(shù)，上午為負(fù)，下午為正。

式中：φ 為地理緯度；γs偏東為負(fù)，偏西為正。

對(duì)北方寒旱區(qū)而言，Hottel 晴空模型因具有較好的適用性和通用性而被優(yōu)先采用［6，21，29］，可彌補(bǔ)國(guó)內(nèi)站點(diǎn)稀疏、數(shù)據(jù)較少等不足。該模型中太陽(yáng)射線到達(dá)大氣層外切平面的太陽(yáng)輻射強(qiáng)度G0由下式計(jì)算：

式中：Gsc為太陽(yáng)常數(shù)，表示單位時(shí)間單位面積上所接收的太陽(yáng)輻照度，取1367 W/m2。

任一斜面的太陽(yáng)入射角i、太陽(yáng)輻照度Gs（太陽(yáng)直射輻照度Gtb和散射輻照度Gtd）公式如下［29-30］：

式中：β 為斜面傾角；γt為斜面方位角，指渠道表面法線在水平面上的投影與正南方向的夾角，面向東時(shí)為負(fù)，面向西時(shí)為正；斜面上的日照開(kāi)始和結(jié)束時(shí)間以i=90°來(lái)計(jì)算［30］；α0、α1、k 為標(biāo)準(zhǔn)晴空大氣常數(shù)，根據(jù)文獻(xiàn)［6，29］計(jì)算。

2.2.2 渠道陰影計(jì)算方法 渠道邊坡相互遮擋產(chǎn)生的陰影分布影響了渠道表面的太陽(yáng)輻射分布。在太陽(yáng)高度角較小時(shí)，陰坡的遮擋使渠道局部形成陰影，且此陰影的邊界線與渠道的走向平行?；诖?，可將渠道陰坡看作一根桿，將此桿在太陽(yáng)下形成的陰影邊界點(diǎn)做渠道走向的平行線，此為渠道陰坡產(chǎn)生的陰影邊界。以北方東西走向渠道為例，為增加圖幅立體感，將方向逆時(shí)針旋轉(zhuǎn)，計(jì)算簡(jiǎn)圖見(jiàn)圖1。

圖1 不同走向渠道陰影長(zhǎng)度計(jì)算示意

圖中OA 為細(xì)桿，高度為渠深h，l1和l4分別為渠坡頂在水平面上的投影。若此時(shí)太陽(yáng)位于渠道正南方，細(xì)桿將產(chǎn)生OB 的陰影，渠道陰影區(qū)域?yàn)閘2—l4，l1—l2為受光區(qū)；若將渠道順時(shí)針旋轉(zhuǎn)一定角度，此值為面向太陽(yáng)的渠道表面方位角γt。為便于圖幅描述，亦可表示為渠道走向不變，太陽(yáng)逆時(shí)針旋轉(zhuǎn)此角度。此時(shí)細(xì)桿OA 將產(chǎn)生OC 的陰影，渠道陰影區(qū)域?yàn)閘3—l4，陰影邊界離細(xì)桿距離為OD。因太陽(yáng)高度角未發(fā)生變化，細(xì)桿產(chǎn)生的陰影長(zhǎng)度一致，即OB=OC。綜合考慮渠道走向和太陽(yáng)位置變化，其產(chǎn)生的陰影長(zhǎng)度OD 由下式計(jì)算：

進(jìn)一步結(jié)合渠道斷面進(jìn)行陰影判定，如圖2所示。A、E 均為太陽(yáng)光線可穿過(guò)位置，為簡(jiǎn)化圖幅，僅以A 點(diǎn)為例介紹。D1、D2、D3為不同太陽(yáng)位置下產(chǎn)生的陰影點(diǎn)，結(jié)合陰影長(zhǎng)度OD1、OD2、OD3和渠道斷面幾何關(guān)系，確定不同時(shí)刻下渠道的受光區(qū)R，計(jì)算公式如下：

圖2 渠道陰影計(jì)算示意

聯(lián)立式（13）—式（18）及（21）（22）得到不同地區(qū)、時(shí)間下的不同走向、斷面大小的渠道表面太陽(yáng)入射角和受光面積，二者共同決定渠道表面的日照時(shí)間，結(jié)合式（19）（20）得到各表面接受的太陽(yáng)輻照度。2.2.3 襯砌板間輻射和環(huán)境熱輻射 因襯砌板間溫度及其與環(huán)境溫度不同而產(chǎn)生熱輻射，計(jì)算簡(jiǎn)圖見(jiàn)圖3，任一點(diǎn)P 得到的S ′面的輻照度（Gm：W/m2）及周?chē)h(huán)境輻照度（Gamb：W/m2）見(jiàn)下式［31］：

式中：J ′、Jamb分別為S′面和環(huán)境（Samb）的輻射度，W/m2；n 和n′為外法線矢量；r 為面上兩點(diǎn)的距離矢量；n 為折射率，不透明物體取1；σ為Stefan-Boltzmann 常數(shù)，取5.67×10-8W/（m2·K4）；Famb為環(huán)境角系數(shù)；Tamb為環(huán)境溫度，K。2.2.4 基于輻射度方法的表面輻射換熱方程 輻射度J 由自身輻射及對(duì)輻照度（Gm、Gs、Gamb）的反射輻射組成，而每一點(diǎn)Gm又是其他可見(jiàn)點(diǎn)J 的函數(shù)，據(jù)此得到如下輻射平衡方程［31］：

式中：ρd為漫反射系數(shù)；ε 為發(fā)射率；T 為渠道表面溫度，K。

假設(shè)渠道表面為理想漫射灰體，發(fā)射率與吸收率α 相等，可計(jì)算得到渠道表面吸收的輻射量，并以此作為凍土水-熱-力耦合中熱模塊的第二邊界條件（熱通量邊界條件），方程如下：

圖3 角系數(shù)計(jì)算示意

3 數(shù)值模型

以新疆阜康某灌區(qū)（87oE、44oN）E-W 走向“陰陽(yáng)坡效應(yīng)”顯著的大型輸水渠道為例，基于現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)，聯(lián)立上節(jié)考慮太陽(yáng)輻射及陰影遮蔽的渠道水-熱-力耦合方程，并采用COMSOL 求解因太陽(yáng)輻射而導(dǎo)致的不均勻水、熱、變形分布，分析渠道凍害發(fā)生原因，并給出熱邊界的選取建議。

3.1 基本概況該梯形渠道修建于季節(jié)性?xún)鐾羺^(qū)，屬溫帶大陸性干旱氣候，多年平均溫度6.5 ℃，極端最低溫度-25～-32.2 ℃，多年最大凍深為1.72 m，平均凍深為1.45 m。渠基土中粒徑小于0.075 mm的土粒質(zhì)量占比大于10%，為凍脹敏感性土。

渠道設(shè)計(jì)流量為21 m3/s，加大流量26 m3/s，渠深3.0 m，正常水深2.6 m，底寬4.0 m，坡比1∶1.5。采用C20 現(xiàn)澆混凝土襯砌，板厚12 cm。姜海波等［10-11］于2013年11月5日—2014年3月28日對(duì)該渠道的溫度場(chǎng)（WS-4 型遙感土壤溫度計(jì)，測(cè)量深度距基土表面5～120 cm）、凍深（DTM-2 型凍土器）、水分場(chǎng)（土鉆取樣后烘干，測(cè)量距基土表面5～120 cm）及襯砌板凍脹量（水準(zhǔn)儀）等進(jìn)行監(jiān)測(cè)，監(jiān)測(cè)點(diǎn)布設(shè)及渠道有限元網(wǎng)格如圖4所示。凍深區(qū)域內(nèi)網(wǎng)絡(luò)高度最小為0.5 cm，最大為1 cm，最大時(shí)間步長(zhǎng)120 s。其中A、B、C 為陽(yáng)坡測(cè)點(diǎn)，分別位于1、2/3、1/3 設(shè)計(jì)水位處，D 為底板中心，E、F、G 為陰坡測(cè)點(diǎn)，位置同陽(yáng)坡對(duì)稱(chēng)。

圖4 有限元網(wǎng)格及數(shù)值模型示意（單位：m）

3.2 計(jì)算參數(shù)根據(jù)試驗(yàn)段土工試驗(yàn)［10-11］可知，淺層基土天然體積含水率為0.256～0.368，處于潮濕-飽和狀態(tài)，液性指數(shù)為0.53～0.82，塑性指數(shù)為5.8～6.4，屬低液限粉土。為簡(jiǎn)化計(jì)算，計(jì)算參數(shù)近似取試驗(yàn)均值，其中黏聚力為10.0 kPa，內(nèi)摩擦角為20°，干密度為1.6 g/cm3；結(jié)合基土土質(zhì)，其彈性參數(shù)與溫度的關(guān)系采用式（27）（28）計(jì)算［2，13，18］；土體各相、漫反射系數(shù)和發(fā)射率、凍結(jié)曲線及土水特征曲線的參數(shù)見(jiàn)表1［2，9，24，28，32］，混凝土參數(shù)見(jiàn)表2［2，20］。

表1 土體材料計(jì)算參數(shù)

表2 混凝土材料計(jì)算參數(shù)

3.3 邊界和初始條件渠道上表面溫度邊界包括熱輻射（式（26））及熱對(duì)流邊界，其中對(duì)流換熱方程如下：

式中：n 為渠道上表面法向向量；Tamb和T 分別為外界溫度和渠道表面溫度，℃；hc為對(duì)流換熱系數(shù)，W/（m2·℃），由風(fēng)速來(lái)近似確定［2］，公式如下：

圖5 外界最高、最低溫度變化曲線

根據(jù)現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)，渠頂和渠底風(fēng)速近似取5 m/s 和1 m/s，并采用二次拋物線函數(shù)過(guò)渡。

該渠道已修建多年，基土溫度基本達(dá)到吞吐平衡，因此以年周期平衡溫度場(chǎng)作為計(jì)算的初始溫度場(chǎng)（11月1日）。根據(jù)氣象站溫度監(jiān)測(cè)值，采用正弦函數(shù)擬合溫度變化曲線［9，19］，公式如下：

式中：T0為年平均溫度，取6.5 ℃；B 為年溫度振幅，取22.5 ℃；tm為時(shí)間，月；α0為相角，取π。

選取試驗(yàn)段渠道岸邊百葉箱每日溫度實(shí)測(cè)極值［10］，整理如圖5所示。日溫度變化曲線結(jié)合每日極值溫度采用正弦函數(shù)擬合［6］，以此作為計(jì)算分析的外界溫度，公式如下：

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

渠道因防滲措施不當(dāng)，滲漏較為嚴(yán)重，停水后，渠底下4.6 m 深處土體基本處于飽和狀態(tài)，水分場(chǎng)基本穩(wěn)定。此后按設(shè)計(jì)水位模擬通水一個(gè)月，并在凍結(jié)試驗(yàn)開(kāi)始前一周排干水分，以此作為計(jì)算分析的初始水分場(chǎng)。

因復(fù)合土工膜較薄（<3 mm），采用無(wú)厚度土工膜滲阻計(jì)算，土工膜前區(qū)域?yàn)榛炷烈r砌，膜后區(qū)域?yàn)榍?，質(zhì)量守恒方程建立如下：

式中：下標(biāo)b、u、d 分別代表土工膜、膜前和膜后區(qū)域；n 為邊界外法線矢量；k 為滲透系數(shù)，混凝土取1.2×10-9m/s，土工膜取1.0×10-11m/s；h 為壓力水頭，m；▽h為滲流域在邊界處的壓力梯度；R 為滲流阻抗，s；d 為厚度，m。

在渠道周邊施加法向位移約束，見(jiàn)圖4，滿足初始地應(yīng)力平衡。

4 結(jié)果分析

4.1 溫度場(chǎng)分析

4.1.1 渠道陰影分布、太陽(yáng)輻射及溫度變化 因渠道為二維斷面，采用襯砌板吸收的太陽(yáng)輻射量沿渠周分布來(lái)定量描述不同時(shí)刻渠道的陰影分布，并分析由此而導(dǎo)致的溫度差異。因東西走向渠道上、下午接收的太陽(yáng)輻射量基本一致，且阜康地區(qū)比北京約晚2 h，因此取10∶00—14∶00 時(shí)段分析太陽(yáng)輻射量；因襯砌板溫度存在累積作用，取10∶00—18∶00 時(shí)段分析。以12月10日為例，結(jié)果如圖6—8所示。其中，坐標(biāo)原點(diǎn)從陽(yáng)坡渠頂算起，下同。

由圖6 可知，從早上太陽(yáng)初升，隨時(shí)間增加，太陽(yáng)直射區(qū)域由陽(yáng)坡逐漸擴(kuò)大至渠底，陰影區(qū)域逐漸減少；隨太陽(yáng)高度角的增加，太陽(yáng)輻射量逐漸增大，至14∶00 時(shí)最大達(dá)374 W/m2，照射近1/2 渠周；而陰坡始終處于陰影區(qū)域，僅吸收太陽(yáng)散射值，太陽(yáng)輻射遠(yuǎn)小于陽(yáng)坡，輻射量差值達(dá)324.6 W/m2。

圖6 不同時(shí)刻襯砌吸收的太陽(yáng)輻射沿渠周分布曲線

由圖7 可知，因太陽(yáng)輻射區(qū)域及量值的增加，襯砌表面溫度隨之升高，陽(yáng)坡溫度遠(yuǎn)高于陰坡，最大差值為7.5 ℃，“陰陽(yáng)坡效應(yīng)”明顯，且由太陽(yáng)輻射引起的溫度累積效應(yīng)明顯；同時(shí)襯砌板陰、陽(yáng)坡溫度均隨外界氣溫的降低而降低，且二者差值逐漸減小。

圖7 不同時(shí)刻襯砌表面溫度沿渠周分布曲線

由圖8 可進(jìn)一步得知，襯砌表面溫度與外界氣溫和太陽(yáng)輻射存在較好的相關(guān)性，但滯后明顯，二者共同造成了陰陽(yáng)坡的橫向溫度差異。太陽(yáng)輻射時(shí)段內(nèi)，陽(yáng)坡溫度主要受太陽(yáng)輻射影響，而陰坡則主要受外界氣溫?zé)釋?duì)流作用；非太陽(yáng)輻射時(shí)段內(nèi)，陰陽(yáng)坡溫度主要受外界氣溫影響。因此以往分析中不考慮陰陽(yáng)坡或?qū)⑻?yáng)輻射包含在外界氣溫中或取差異性的渠道表面溫度值，均不能科學(xué)合理的反映出陰陽(yáng)坡溫度場(chǎng)的時(shí)空分布差異。

圖8 陰陽(yáng)坡2/3 設(shè)計(jì)水位點(diǎn)溫度、輻射分布曲線

4.1.2 渠道陰陽(yáng)坡襯砌板表面溫度變化 取監(jiān)測(cè)期內(nèi)陰陽(yáng)坡襯砌板最低和日均溫度值來(lái)分析二者的橫向溫度差異，如圖9所示。

圖9 陰陽(yáng)坡襯砌板表面溫度變化曲線

由現(xiàn)場(chǎng)實(shí)測(cè)值可知，襯砌板溫度隨外界氣溫的降低而降低，其變化規(guī)律與外界氣溫基本一致，但存在滯后性。陰、陽(yáng)坡板最低溫度分別為-24.6 和-21.5 ℃，溫差最大為4.3 ℃。

由數(shù)值結(jié)果可知，最低溫度值略高于現(xiàn)場(chǎng)實(shí)測(cè)值，可能是由于太陽(yáng)輻射模型中未考慮空氣透明度、陰天、降雪等因素引起的太陽(yáng)輻射值偏高所致。同時(shí)溫度變化規(guī)律同現(xiàn)場(chǎng)實(shí)測(cè)基本一致，且陰陽(yáng)坡板日均溫度差異明顯，最大為3.5 ℃，進(jìn)一步表明本模型可較好的分析由太陽(yáng)輻射的時(shí)空效應(yīng)和陰坡遮蔽作用而引起的襯砌板橫向溫度差異。

4.1.3 渠道陰陽(yáng)坡凍深及基土溫度場(chǎng)變化 渠道表面溫度場(chǎng)的不對(duì)稱(chēng)導(dǎo)致基土溫度場(chǎng)的差異性分布。以陰、陽(yáng)坡2/3 設(shè)計(jì)水位點(diǎn)為例，將0 ℃等溫線作為凍深判別依據(jù)，其發(fā)展如圖10所示。

圖10 陰、陽(yáng)坡凍深發(fā)展過(guò)程線

因晝夜溫差及太陽(yáng)輻射作用，凍結(jié)前期陰、陽(yáng)坡呈現(xiàn)出“夜凍晝消”的現(xiàn)象。由圖可知，凍深隨負(fù)積溫的累積而逐漸向下發(fā)展，“夜凍晝消”現(xiàn)象逐漸消失，其中陰坡于11月25日消失，陽(yáng)坡于12月10日消失，二者前后差異15 d。結(jié)合圖5，至1月10日最低溫時(shí)刻，凍深并未發(fā)展至最大，后期雖存在大幅度升溫現(xiàn)象，但并未達(dá)到凍土的融化溫度，凍深依然向下發(fā)展；至3月16日，凍深最大，為135 cm；隨著溫度的繼續(xù)升高，陽(yáng)坡于3月17日開(kāi)始融化，陰坡于3月18日開(kāi)始融化，二者前后差異1 d。外加太陽(yáng)輻射作用，使得基土的融化速率遠(yuǎn)大于凍結(jié)速率，至監(jiān)測(cè)期的最后一天3月28日，陰、陽(yáng)坡基土表面分別融化40 和49 cm，凍深分別為97 和66 cm，差異明顯。

為進(jìn)一步分析太陽(yáng)輻射作用下的基土溫度場(chǎng)分布規(guī)律，取最大凍深時(shí)刻，如圖11所示。

圖11 3月16日渠道溫度場(chǎng)分布曲線

由圖可知，陰坡溫度明顯低于陽(yáng)坡，溫度場(chǎng)呈不對(duì)稱(chēng)分布；基土溫度隨外界氣溫（-8 ℃）變化存在滯后性，凍深以上溫度變化劇烈，以下則較為緩慢。結(jié)合凍深監(jiān)測(cè)結(jié)果可知，數(shù)值結(jié)果略小，最大差異為22 cm，主要原因在于太陽(yáng)輻射模型未考慮陰天、降水對(duì)太陽(yáng)光的阻擋所致。但反映出的凍深分布規(guī)律與監(jiān)測(cè)結(jié)果分布規(guī)律基本一致，僅在數(shù)值上存在微小差異，表明本文模型可較好的模擬太陽(yáng)輻射作用下的基土溫度場(chǎng)分布規(guī)律。

4.2 水分場(chǎng)分析太陽(yáng)輻射和環(huán)境溫度共同作用引起渠道的不對(duì)稱(chēng)溫度場(chǎng)，勢(shì)必引起不對(duì)稱(chēng)的水分場(chǎng)分布。為此取11月5日（初始時(shí)刻）和1月20日（監(jiān)測(cè)時(shí)刻）陽(yáng)坡、渠底及陰坡的測(cè)點(diǎn)A、B、D、F、G 的總含水率沿基土深度分布結(jié)果，如圖12所示。

圖12 關(guān)鍵點(diǎn)渠基土體積含水率分布曲線

由圖可知，相比于初始時(shí)刻，凍結(jié)區(qū)域內(nèi)因外界負(fù)積溫作用而使未凍區(qū)內(nèi)的水分逐漸向上遷移，水分集聚現(xiàn)象明顯。由數(shù)值結(jié)果可知，因太陽(yáng)輻射和晝夜溫度的疊加作用，凍土內(nèi)存在多個(gè)水分集聚位置，而非只發(fā)生在凍結(jié)鋒面處；陽(yáng)坡、渠底及陰坡水分的集聚個(gè)數(shù)、位置及深度均不同，呈現(xiàn)出陰坡含水量最多而陽(yáng)坡最少的分布規(guī)律。對(duì)比監(jiān)測(cè)結(jié)果可知，因現(xiàn)場(chǎng)選取局部點(diǎn)來(lái)烘干，僅能反映出水分運(yùn)移規(guī)律，并不能?chē)?yán)格反映出含水量最高點(diǎn)或凍結(jié)鋒面位置，但在凍深范圍內(nèi)呈現(xiàn)陰坡、渠底總含水量高而陽(yáng)坡較少的分布規(guī)律，與數(shù)值結(jié)果一致。然數(shù)值解的表層含水量偏高，80～100 cm 深度處偏低，可能與未考慮陰天、蒸發(fā)等因素及現(xiàn)場(chǎng)局部地質(zhì)引起的初始含水偏高所致；但其他位置點(diǎn)的含水量吻合度較高，滿足工程要求，即此模型可模擬太陽(yáng)輻射下的水分遷移、凍結(jié)相變等特征。

為進(jìn)一步分析不對(duì)稱(chēng)的冰、水含量分布，取最大凍深時(shí)刻（3月16日）和最大融深時(shí)刻（3月28日）的水分場(chǎng)結(jié)果，如圖13—14所示。

由圖13 可知，凍結(jié)區(qū)域內(nèi)的未凍水含量隨溫度降低而逐漸減少，最低為0.06，接近于基土的殘余含水量；而相應(yīng)位置的結(jié)冰量逐漸增加，冰含量呈條帶狀分布，無(wú)明顯的空間分布規(guī)律［2-3，12］，但與基土凍深發(fā)展相呼應(yīng)（圖10）；因凍結(jié)初期基土表層常處于凍融循環(huán)狀態(tài)，在多次溫度抽吸力下表層冰含量較大；冰、水總含量最大為0.6，已遠(yuǎn)大于基土的孔隙率，將發(fā)生嚴(yán)重的凍脹變形。

3月16日表層基土出現(xiàn)局部融化現(xiàn)象，含水量增加而冰含量減少。隨著外界溫度的逐漸升高，凍結(jié)層內(nèi)冰透鏡體開(kāi)始從兩端發(fā)生雙向融化，局部區(qū)域形成過(guò)量孔隙水，隨著時(shí)間推移，孔隙水逐漸擴(kuò)散，形成隨深度增加含水量逐漸增大的分布規(guī)律（圖14（a）），因凍結(jié)區(qū)域內(nèi)冰含量的阻水作用，使含水量較初始時(shí)刻明顯增大。對(duì)比圖13（b）和圖14（b）可知，凍深范圍外冰透鏡體全部融化，且凍深線附近因溫度升高至其融化溫度，使冰含量值明顯減少，而凍深范圍內(nèi)的其余區(qū)域因熱流傳輸密度少和凍深線附近的冰融化吸熱，從而使冰含量未發(fā)生較大變化。

圖13 3月16日渠道水分場(chǎng)分布

圖14 3月28日渠道水分場(chǎng)分布

由圖13、圖14 亦可發(fā)現(xiàn)，因太陽(yáng)輻射和陰坡遮蔽而產(chǎn)生的熱邊界差異產(chǎn)生了明顯的不對(duì)稱(chēng)水分場(chǎng)分布，陰坡總含水量明顯大于陽(yáng)坡。首先，陰陽(yáng)坡的凍深差異特征形成了陰坡冰透鏡體總體厚度較陽(yáng)坡大的現(xiàn)象；其次，陰陽(yáng)坡的熱邊界差異使陰坡含冰量大于陽(yáng)坡含冰量；最后，因陽(yáng)面接受太陽(yáng)輻射多，溫升快，從而使陽(yáng)坡融化速度明顯大于陰坡，冰含量快速減少。

4.3 變形場(chǎng)分析 由上節(jié)可知，陰、陽(yáng)坡局部位置總含水量大于孔隙率，將會(huì)產(chǎn)生差異性?xún)雒涀冃危x取多個(gè)時(shí)刻下襯砌板的法向凍脹變形沿渠周分布曲線，如圖15所示。

圖15 襯砌板法向凍脹變形沿渠周分布曲線

由圖可知，梯形渠道襯砌因基土凍脹作用，呈現(xiàn)出渠底向上隆起，坡腳受擠壓明顯，兩坡板向渠內(nèi)凸起，襯砌整體上抬的變形趨勢(shì)。其中，渠底襯砌板在偏陽(yáng)坡處的法向凍脹量最大，陰、陽(yáng)坡板約在1/4～1/3 坡長(zhǎng)位置處最大。在因太陽(yáng)輻射而產(chǎn)生的差異凍脹力作用下，襯砌凍脹變形量由陰坡、渠底至陽(yáng)坡依次遞減，陰、陽(yáng)坡板凍脹變形差異明顯，整體與基土內(nèi)凍深發(fā)展、水分分布規(guī)律一致，即渠道水熱分布的不對(duì)稱(chēng)、不同步特征導(dǎo)致了陰、陽(yáng)坡板不對(duì)稱(chēng)的變形分布。

由上圖亦可知，隨著凍結(jié)時(shí)間的發(fā)展，襯砌板法向凍脹變形均逐漸增大，至3月16日時(shí)，凍脹量最大，陽(yáng)坡、渠底和陰坡最大分別為8.47、10.97 和12.48 cm，與現(xiàn)場(chǎng)監(jiān)測(cè)值基本吻合，且襯砌板的整體變形趨勢(shì)與監(jiān)測(cè)值分布規(guī)律基本一致，同時(shí)，在后續(xù)融化過(guò)程中，襯砌板凍脹變形亦隨基土內(nèi)冰水總含量的減少而逐漸變小，表明本文的數(shù)值模型可較好的模擬襯砌結(jié)構(gòu)在基土凍脹作用下的凍融變形特征。

5 考慮太陽(yáng)輻射渠道凍脹模型的熱邊界選取及合理性評(píng)價(jià)

太陽(yáng)輻射是渠道凍脹分析中能量的來(lái)源，影響著基土內(nèi)部溫度場(chǎng)、水分場(chǎng)及變形場(chǎng)的耦合作用，是進(jìn)行寒區(qū)渠道設(shè)計(jì)的重要因素，其中熱邊界合理選取至關(guān)重要。現(xiàn)階段，凍脹數(shù)值模擬中常用以下兩種熱邊界：（1）以日或月為單位擬合環(huán)境溫度，只采用熱對(duì)流邊界［2-3，9，13］，但無(wú)法考慮太陽(yáng)輻射產(chǎn)生的陰陽(yáng)坡效應(yīng)；（2）結(jié)合現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)，采用正弦函數(shù)擬合陰、陽(yáng)坡溫度變化曲線，并將其直接賦予結(jié)構(gòu)邊界［17-18］，但卻無(wú)法充分反映出太陽(yáng)輻射的空間效應(yīng)和晝夜溫差的時(shí)間效應(yīng)。目前公認(rèn)更為準(zhǔn)確的是第二種以現(xiàn)場(chǎng)地溫監(jiān)測(cè)為依據(jù)的擬合溫度邊界，以此為參照，對(duì)比分析本文采用的太陽(yáng)輻射、面-面輻射、環(huán)境輻射的熱輻射邊界和晝夜溫差的熱對(duì)流邊界組合，對(duì)熱邊界選取及其合理性進(jìn)行評(píng)價(jià)。擬合溫度邊界取上章仿真結(jié)果的陰坡、渠底和陽(yáng)坡的日均表面溫度數(shù)據(jù)（圖9）來(lái)正弦擬合熱邊界，可彌補(bǔ)現(xiàn)場(chǎng)監(jiān)測(cè)點(diǎn)少且數(shù)據(jù)不全的缺點(diǎn)。下面將從溫度、水分、變形來(lái)分別闡述。

取陰、陽(yáng)坡2/3 設(shè)計(jì)水位點(diǎn)的凍深發(fā)展過(guò)程結(jié)果，如圖16所示。由圖可知，擬合的溫度邊界使襯砌板溫度緩慢變化，凍結(jié)速率發(fā)展較為均勻；而本文模型因考慮了太陽(yáng)輻射空間效應(yīng)和晝夜溫度時(shí)間效應(yīng)的耦合疊加作用，使得凍深發(fā)展在局部位置出現(xiàn)了短暫的平臺(tái)狀，凍結(jié)速率忽快忽慢，但與外界的環(huán)境變化相協(xié)調(diào)，更符合現(xiàn)場(chǎng)實(shí)際規(guī)律。

圖16 不同熱邊界下的凍深發(fā)展過(guò)程線

陰、陽(yáng)坡含水量分布規(guī)律基本類(lèi)似，取最大凍深時(shí)刻下陰坡2/3 設(shè)計(jì)水位點(diǎn)的總含水量沿基土深度的分布曲線，如圖17所示。由圖可知，因擬合溫度邊界下凍結(jié)速率變化較為均勻，除淺層基土外，含水率由上至下逐漸增大，至凍結(jié)峰面處最大，與單向凍結(jié)結(jié)果基本一致。而本文模型因考慮太陽(yáng)輻射和晝夜溫差作用，凍結(jié)初期的“夜凍晝消”往復(fù)凍融循環(huán)作用使得基土表層含水增加；水分分布與凍土的凍結(jié)速率相對(duì)應(yīng)，凍結(jié)速率慢時(shí)水分遷移多，反之亦然，更能反映出外界環(huán)境作用引起的水熱耦合作用。

圖17 不同熱邊界下的水分場(chǎng)分布

總含水率和凍深差異導(dǎo)致了不同的凍脹變形，取陰、陽(yáng)坡2/3 設(shè)計(jì)水位點(diǎn)，其襯砌板法向凍脹變形隨時(shí)間的分布規(guī)律如圖18所示。由圖可知，襯砌板的起始凍脹時(shí)刻、凍脹量發(fā)展速度和量級(jí)，以及融化速度等變量均與熱邊界有關(guān)，擬合溫度邊界方法得到的凍脹量不能反映出太陽(yáng)輻射和晝夜溫差作用所引起的凍結(jié)速率差異而導(dǎo)致的水熱分布規(guī)律，導(dǎo)致結(jié)構(gòu)凍融破壞計(jì)算不準(zhǔn)確。

圖18 不同工況下襯砌板法向凍脹量隨時(shí)間變化曲線

結(jié)合溫度、水分和變形結(jié)果可知，本文數(shù)值模型的結(jié)果與目前公認(rèn)的以地溫監(jiān)測(cè)數(shù)據(jù)為基礎(chǔ)的擬合溫度邊界所得規(guī)律基本一致，結(jié)合上一節(jié)結(jié)果，可進(jìn)一步說(shuō)明本文模型的準(zhǔn)確性；同時(shí)本文模型綜合考慮了太陽(yáng)輻射的空間效應(yīng)和晝夜溫度的時(shí)間效應(yīng)，可反映出外界環(huán)境作用下寒區(qū)渠道水、熱、力耦合的凍脹破壞本質(zhì)，特別是能減少對(duì)現(xiàn)場(chǎng)監(jiān)測(cè)值的依賴(lài)性，且具有理論通用性，相比于目前使用較多的擬合溫度邊界，結(jié)果更加合理、準(zhǔn)確，且可輔助設(shè)計(jì)單位進(jìn)行前期渠道設(shè)計(jì)，滿足精細(xì)化設(shè)計(jì)、復(fù)核的要求。

6 結(jié)論

考慮太陽(yáng)輻射和晝夜溫差的時(shí)空變化，提出了渠道的陰影計(jì)算方法，基于凍土的橫觀各向同性?xún)雒浱卣鳎⒘恕疤?yáng)輻射-環(huán)境溫度-凍土-工程”相互作用的水-熱-力耦合凍脹模型，并探討了不同熱邊界對(duì)渠道水-熱-力耦合計(jì)算的影響，結(jié)合工程監(jiān)測(cè)及數(shù)值模擬結(jié)果，得到如下結(jié)論：（1）考慮太陽(yáng)輻射的渠道凍脹模型計(jì)算的溫度、水分和變形結(jié)果與現(xiàn)場(chǎng)實(shí)測(cè)值較為吻合，分布規(guī)律一致，表明本文模型可分析太陽(yáng)輻射空間效應(yīng)和晝夜溫差時(shí)間效應(yīng)下的寒區(qū)渠道凍結(jié)過(guò)程，理論通用性強(qiáng)，且可減少對(duì)現(xiàn)場(chǎng)監(jiān)測(cè)的依賴(lài)，可為寒區(qū)渠道的凍脹模擬、預(yù)測(cè)及災(zāi)害防控設(shè)計(jì)提供理論依據(jù)。（2）探明了渠道表面太陽(yáng)輻射的時(shí)空分布規(guī)律，陰、陽(yáng)坡太陽(yáng)輻射時(shí)長(zhǎng)、輻射量存在顯著差異，其表面溫度由外界氣溫和太陽(yáng)輻射量共同決定，計(jì)算的陰、陽(yáng)坡日均溫差最大為3.5 ℃，導(dǎo)致了陰坡早于陽(yáng)坡15 d 凍結(jié)，晚于陽(yáng)坡1 d 融化，凍深發(fā)展差異明顯，最大相差31 cm，渠道陰陽(yáng)坡溫度場(chǎng)分布不對(duì)稱(chēng)且不同步。（3）凍結(jié)期渠道陰坡結(jié)冰區(qū)域及量值明顯大于陽(yáng)坡，“夜凍晝消”的凍融循環(huán)作用使基土淺層含水量明顯增多，而融化期陽(yáng)坡冰融化速度大于陰坡，水分場(chǎng)呈不對(duì)稱(chēng)、不同步、層狀分布等特點(diǎn)。（4）渠坡板和底板的最大凍脹變形分別發(fā)生在1/4～1/3 坡長(zhǎng)位置和偏陽(yáng)坡位置，陰、陽(yáng)坡變形隨時(shí)間的發(fā)展過(guò)程呈不均勻、不同步、不對(duì)稱(chēng)性，計(jì)算的最大變形分別為12.48 和8.47 cm。（5）渠道陰陽(yáng)坡效應(yīng)產(chǎn)生的本質(zhì)在于太陽(yáng)輻射和晝夜溫差引起的渠道能量收支差異、溫度梯度和凍結(jié)速率變化所致，即熱邊界引起的渠道基土的水-熱-力動(dòng)態(tài)耦合作用。（6）前期設(shè)計(jì)時(shí)，宜采用多年氣溫?cái)M合值和太陽(yáng)輻射參數(shù)確定熱邊界；已建工程凍脹評(píng)價(jià)時(shí)，宜采用日氣溫值和太陽(yáng)輻射參數(shù)確定熱邊界，并輔以監(jiān)測(cè)值對(duì)其修正。