吐魯番市典型灌區(qū)輸水渠道冬季冰蓋多場耦合下力學(xué)響應(yīng)特性分析研究

李 鳴

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

1 概 述

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

2 工程概況

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

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

圖1 冰蓋體微觀示意圖

3 多場耦合力學(xué)響應(yīng)理論

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

(1)

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

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

(2)

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

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

(3)

其中應(yīng)力場變形服從以下方程式：

(4)

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

(5)

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

(6)

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

(7)

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

(8)

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

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

4 力學(xué)特征分析

4.1 研究工況與建模

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

圖2 渠道特征斷面

4.2 多場耦合下冰蓋體應(yīng)力特征

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

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

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

圖4 冰蓋體壓應(yīng)力分布云圖

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

圖5 靜冰壓力變化曲線

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

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

4.3 渠道結(jié)構(gòu)參數(shù)對冰蓋體應(yīng)力響應(yīng)影響

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

4.3.1 邊坡系數(shù)影響

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

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

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

4.3.2 渠道底寬影響

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

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

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

5 結(jié) 論

針對吐魯番市典型灌區(qū)輸水渠道冬季產(chǎn)生的冰蓋體，利用COMSOL Multiphyisics有限元數(shù)值軟件，開展熱-力耦合下冰蓋體力學(xué)響應(yīng)分析，得到以下幾點結(jié)論：

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

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

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

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