吐魯番市典型灌區輸水渠道冬季冰蓋多場耦合下力學響應特性分析研究

李 鳴

(吐魯番市高昌區水管總站，新疆 吐魯番 838000)

1 概 述

為提升農業灌溉效率，我國許多地區都建有長距離輸水渠道，調控水資源至農業生產用水[1-3]。但在西北等冬季高寒地區，冬季輸水勢必會形成河冰及較厚的冰蓋，對渠道安全輸水帶來較大影響，因而開展輸水渠道冰蓋體研究具有重要作用[4-5]。已有許多農業水利工程師針對河冰等流冰體開展過渠道結冰預警等研究，還有一些學者通過現場安裝監測設備，研究冰蓋在升溫過程中內部溫度場以及應力變形等變化，探討各參數變化規律，為渠道穩定性設計提供重要參考[6-9]。針對多物理場條件，以有限元數值軟件作為計算手段，研究不同工況下冰蓋體力學以及溫度場特征，從數值分析角度解釋冰蓋體對渠道安全運營的威脅[10-12]，為渠道安全建設提供理論依據。

2 工程概況

吐魯番市典型代表灌區內有5條河流，最大年流量達3.31 m3/s，區域內春冬季溫差較大，年風速超過1.5 m/s，工程結構需耐冷熱凍脹效應影響。灌區內水資源豐枯水期層次性顯著，可供應水量之比為2∶1，其中地表水主要來源于5條河流，為流域內農業生產提供灌溉用水。區域內已建有長度超過240 km的輸水渠道，均采用格賓石籠針對性襯砌防滲結構，其中干渠有9條，分別支援不同鄉鎮農業生產用水；支渠超過35條，渠道斷面形式為梯形，渠底寬度為40 m，水深設計為4 m，上覆蓋層冰蓋體平均50 cm，年灌溉水資源超過24 000×104m3，滿足2.933 3×104hm2農業土地用水需求。局部地區水資源利用效率較低，灌溉效率每年以2%速率下降，這主要受地下水資源近幾年有節制調度使用，而灌溉用水中地下水資源又占比較大。另一方面，灌區輸水渠道襯砌結構在低溫冬季易形成冰蓋體，當溫度變化時，冰蓋體與渠道之間互相作用，降低水資源在渠道內輸送效率，而實質上輸水渠道與冰蓋體為多物理場耦合作用下(溫度場與固體場)，其力學特征關乎水資源輸送效率。

為準確評判輸水渠道在冬季低溫環境下運營效率，針對渠道開展多場耦合數值試驗分析。根據現場地質踏勘發現，灌區內地勢并不較高，處于盆地低高程區域，高差100～150 m，地質構造主要在東南部丘陵地帶可見背斜等地質現象。灌區內鉆孔取樣資料表明，覆蓋土層以人工雜填土、粉質壤土及砂礫土為主，其中人工雜填土在表面分布范圍較廣，主要為一些種植填土與碎石填土，密實性較差，析水性較強，分布厚度約為2.5 m；粉質壤土主要存在于東南部丘陵進入灌區內的支渠周圍，是部分支渠的渠地基，中等承載力，防滲性較好，含水量為18%；砂礫土現場取樣得知，最大粒徑超過6 mm，級配較差，與下臥砂巖層界面錯動接觸，松散型較大，覆蓋土層地質年代基本均屬于第四系。巖層以二疊統砂巖以及角礫石為主，強風化作用，完整性較差。地下水位最深處達16 m，基巖內孔隙水分布較少。另一方面調查發現，冬季輸水渠道常出現冰蓋體，部分渠道區段內襯砌結構在冬季易受損壞，冰蓋體在春季溫度升高時轉換成流體運動，冰蓋體顆粒尺寸不超過20 mm，厚度不等，淺層冰蓋體內部具有較多的氣泡，常見的冰蓋體微觀示意圖見圖1，冰蓋體產生的冰壓力分布是渠道襯砌結構安全考慮的重要方面。

圖1 冰蓋體微觀示意圖

3 多場耦合力學響應理論

冰蓋體是一種特殊狀態下固體結構，本文在分析過程中認定冰蓋體為各向同性天然材料，在冬季至春季溫度變化過程中處于流變變形，而描述流變變形的重要手段即是蠕變本構模型。本文以非線性Burger模型作為其蠕變本構模型，其具體表述可用下式：

(1)

式中:cice指常數參數；dref、dice為冰蓋體尺寸；Asinha為黏滯系數。

式(1)中，表述了冰蓋體所在流變過程中發生的彈性變形與塑性不可逆變形，直至失穩破壞狀態。常用的冰蓋體蠕變耦合模型還有以冰體流變速率表述的流變本構模型，其表達式為：

(2)

式中：Anor為流變實系數；n為常數參數，取3；Q為冰蓋體自身能量；T、R分別為物理場參數。

在輸水渠道中，冰蓋溫度變化過程中產生的冰壓力會對襯砌結構產生一定破壞影響，考慮此，應結合溫度場-應力場開展渠道冰蓋體力學特征分析，其中熱力耦合下控制方程為：

(3)

其中應力場變形服從以下方程式：

(4)

冰蓋體在耦合場中會出現長期蠕變變形與熱變形，其總變形可用下列3項式子表述：

(5)

式(5)中，各向同性橫向變形取值參數以矩陣形式表述為：

(6)

上述理論為輸水渠道中冰蓋體熱力耦合下應力分析基本理論，其中還包括有流變應力，要求解上述材料的本構模型方程，在有限元數值軟件中需要建立耦合場，本文引入COMSOL Multiphyisics有限元數值軟件，通過多物理場的微分方程進行模擬研究環境，并求解應力變形等特征參數。而在水利工程中，常見的多物理場微分方程一般呈非線性，其邊界條件及方程式可用下式表述：

(7)

n·(cu+αu-γ)=g-qu+hTμ

(8)

式中：u為物理變化中的因變量；ea、da、c、α、β、γ、f均為耦合場中特定物理參數。

本文計算時忽略渠道凍脹變形影響，且認為冰蓋體受熱效應影響僅在融點內，忽略流變變形對體積應變張量不產生影響。

4 力學特征分析

4.1 研究工況與建模

模擬工況上下冰蓋體溫差為10℃，表面溫度為-10℃，流變時間為6 d，所選取渠道特征斷面見圖2。分析計算溫度從低溫至融點過程中，即升溫過程中冰蓋體應力響應特征以及渠道結構參數改變對冰蓋體應力特征影響。

圖2 渠道特征斷面

4.2 多場耦合下冰蓋體應力特征

根據COMSOL數值有限元軟件獲得圖3所示冰蓋體在不同時間段時溫度場分布云圖。從圖3中可看出，冰蓋體在流變過程中表面層溫度逐漸增大，初始狀態時處于-10℃，在流變時間為3 d時，溫度增大1.5倍，達-4℃；在流變時間為第5 d時，表面層溫度已達到0℃，即溫度由最底層完成傳遞至冰蓋體頂層。從整體溫度場變化分布來看，在未完成溫度傳遞之前，即第5 d之前，冰蓋體整體溫度分布呈從底部至頂部，逐漸減小，但在冰蓋體完成自上而下的溫度傳遞后，頂部與底部溫度為最高，均已達到0℃。但在冰蓋體中間較厚區域仍然存在負溫，即冰蓋體在長期流變過程中溫度場分布并不是均勻狀態，即使已完成溫度熱傳遞，這也表明溫度場與應力場耦合效應下，冰蓋體應力特征變化過程中，內部溫度場的溫差作用持續發生。

圖3 不同時間段時溫度場分布云圖

圖4為計算獲得流變升溫過程中冰蓋體壓應力分布云圖。從圖4中可看出，初始狀態下拉應力為最大，達150 kPa，最大壓應力為241 kPa，位于距離冰蓋體頂部5 cm區域，而拉應力區域集中在該冰蓋體頂部，即初始狀態下冰蓋體頂部實質上是出于收縮變形才致拉應力出現；隨著升溫與流變時間推移，頂部逐漸由受拉轉變至受壓，在第1 d時壓應力約為50～100 kPa，最大壓應力集中在底部區域，最大拉應力相比初始狀態增大46.7%，達314 kPa；升溫至第3 d時，壓應力逐漸降低，拉應力值亦處于較低水平，僅為7.13 kPa，冰蓋體整體上均處于膨脹狀態，即壓應力分布為主要占比；在第5 d時，冰蓋體頂部亦均為壓應力分布，但在底部還殘留有拉應力，即在溫度傳遞完成之后，整體上拉應力分布亦從頂部至底部轉移。

圖4 冰蓋體壓應力分布云圖

圖5為COMSOL在迭代求解升溫過程中冰蓋體應力特征分布時獲得的最大靜冰壓力曲線。從圖5中可看出，靜冰壓力在升溫過程中約為20～140 kPa，其最大靜冰壓力為142.3 kPa，與文獻[13]計算獲得靜冰壓力范圍具有一致性。另外，亦可看出靜冰壓力變化是分階段，呈先增后減變化。

圖5 靜冰壓力變化曲線

圖6為計算出的蠕變有效位移特征分布云圖。從圖6中可看出，蠕變有效位移從初始狀態至溫度完全傳遞過程中，應變量級增大12個量級，即冰蓋體在升溫過程中會促進蠕變進展，此亦在文獻[14]中得到印證，外界溫度升高，巖石等固體材料流變速率會增大，位移值也會相對提高。另從位移在冰蓋體上分布來看，主要集中在右側靠近渠道邊坡側以及頂面靠近，即與渠道相接觸截面上蠕變位移持續穩定分布， 在冰蓋體頂部右側區域第5 d相比第1 d位移升高1個量級。

圖6 蠕變有效位移特征分布云圖

4.3 渠道結構參數對冰蓋體應力響應影響

由于渠道結構涉及參數眾多，本文以其中渠道側邊坡系數與渠道底寬作為影響冰蓋體應力分析，基于COMSOL數值軟件設定不同的邊坡系數研究工況，對比各工況之間冰蓋體結果差異。

4.3.1 邊坡系數影響

邊坡系數設定為1.5、2.25、3.25、3.5共4個工況，渠道底寬為40 m，其他參數與前述一致，研究升溫過程中冰蓋體應力變化特征，結果見圖7。

圖7 極限冰壓力變化曲線(邊坡系數影響)

從圖7中可看出，各邊坡系數下極限冰壓力變化曲線形態基本一致，均呈先增大后減小，從各邊坡系數下極限冰壓力水平來看，當邊坡系數增大時，極限冰壓力峰值均有所增大。邊坡系數1.5時，極限冰壓力峰值為508.3 kPa；而邊坡系數為2.25時，增大9%，達553.53 kPa；但在邊坡系數3.25后，極限冰壓力峰值逐漸降低，其中邊坡系數3.5時極限冰壓力峰值相比2.25時降低5.3%，為524.3 kPa。由此表明，邊坡系數影響冰蓋體極限冰壓力峰值亦是先增后減，即邊坡系數對冰蓋體極限冰壓力促進作用具有拐點系數，當超過邊坡拐點系數后，冰蓋體極限冰壓力會逐漸降低。從極限冰壓力隨升溫時間變化曲線可看出，各邊坡系數下極限冰壓力在后期流變過程中處于一致性，即極限冰壓力在二次減小階段減小速率為一致，均為每小時減小冰壓力2 kPa。

4.3.2 渠道底寬影響

渠道底寬是灌區輸水渠道設計中重要尺寸，本文以渠道底寬10、20、40和50 m作為研究工況，其他參數均為一致，獲得冰蓋體極限冰壓力變化曲線，見圖8。

圖8 極限冰壓力變化曲線(渠道底寬影響)

從圖8中可看出，隨著渠道底寬增大，冰蓋體極限冰壓力水平逐漸減小，渠道底寬10 m時冰蓋體極限冰壓力峰值為1 220.4 kPa，而底寬40和50 m時分別降低53.7%和62%，且各渠道底寬下冰蓋體極限冰壓力峰值基本均處于同一升溫時刻，即第26 h。在冰蓋體極限冰壓力二次下降階段中，各渠道底寬下冰蓋體的極限冰壓力曲線均為一致性，一次下降階段中，渠道底寬愈小，則下降速率愈大，底寬10 m時冰蓋體在一次下降階段中每小時減小冰壓力15.3 kPa，而底寬40 m時該參數為6.6 kPa。分析表明，渠道底寬尺寸愈小時，溫度熱效應在小斷面內會具有較大變形約束，因而在冰蓋體內部產生較大的冰壓力，即呈現圖8中現象，當完成溫度熱效應傳遞后，即第120 h后，此時溫差傳遞對熱效應減弱，故而冰蓋體極限冰壓力逐漸為一致性變化。

5 結 論

針對吐魯番市典型灌區輸水渠道冬季產生的冰蓋體，利用COMSOL Multiphyisics有限元數值軟件，開展熱-力耦合下冰蓋體力學響應分析，得到以下幾點結論：

1) 研究了多場耦合下冰蓋體升溫過程中即使已完成溫度熱傳遞，溫度場分布并不是均勻狀態，表面層溫度在5 d內完成熱傳遞；冰蓋體初始拉應力分布在底部，壓應力約為50～100 kPa，熱傳遞后拉應力至頂部，壓應力占據冰蓋體大部分區域，但升溫過程中壓應力為先增后減。

2) 獲得了多場耦合下冰蓋體靜冰壓力為20～140 kPa，最大靜冰壓力為142.3 kPa，靜冰壓力變化呈先增后減；流變有效位移在完成熱傳遞后，增大12個量級，且有效位移位于冰蓋體與渠道接觸截面上。

3) 分析了渠道邊坡系數對冰蓋體應力響應影響特征，渠道邊坡系數對冰蓋體極限冰壓力促進作用具有拐點系數2.25，當邊坡系數超過2.25后，極限冰壓力降低，邊坡系數3.5時極限冰壓力峰值相比2.25時降低5.3%；各邊坡系數下極限冰壓力在二次減小階段減小速率為一致，均為每小時減小冰壓力2 kPa。

4) 研究了渠道底寬影響冰蓋體應力特征，各渠道底寬下冰蓋體的極限冰壓力曲線變化均為一致性，隨底寬增大，冰蓋體極限冰壓力水平減小，極限冰壓力峰值均處于第26 h；各渠道底寬的冰壓力在二次下降階段變化為相同狀態。