冬季凍區(qū)混凝土輸水渠道凍脹特性的數(shù)值模擬

杜民瑞，宋 玲，陳瑞考，魏 鵬

（石河子大學(xué)水利建筑工程學(xué)院，新疆 石河子 832003）

由于梯形混凝土襯砌渠道斷面形式近似圓形，故其在工程實(shí)際當(dāng)中梯形渠道被廣泛利用，但是在冬季凍區(qū)普遍存在著的渠道襯砌層凍脹破壞現(xiàn)象，又給渠道的建設(shè)，乃至整個(gè)水利工程設(shè)計(jì)帶來(lái)了更大的挑戰(zhàn)。人們對(duì)于普通的襯砌渠道凍脹模型已經(jīng)做出了大量的研究，RL Harlan[1]于1973首次提出了水熱耦合的概念并得出了耦合模型，Loch[2]提出用克勞修斯–可拉貝龍方程來(lái)描述多孔介質(zhì)冰和水相變過(guò)程的熱動(dòng)平衡，Shen Mu和Branko Ladanyi[3]在此基礎(chǔ)上于1990年提出了水、熱、力三場(chǎng)耦合問(wèn)題，得出了簡(jiǎn)化后的耦合模型。Taylor and Luthin[4]基于Harlan模型提出便于求解的有限差分格式，Konrad[5]等通過(guò)實(shí)驗(yàn)數(shù)據(jù)以及理論分析，提出分凝勢(shì)的概念，并把其定義為土體中水分遷移速率與溫度梯度的比值。安維東[6]對(duì)凍土的溫度、水分、應(yīng)力三場(chǎng)耦合從理論上作了論述，并以某些水利工程的問(wèn)題為實(shí)例進(jìn)行數(shù)值模擬。董江偉[7]考慮了渠基土和所接觸復(fù)合土工膜結(jié)構(gòu)面之間的相互作用，對(duì)大斷面剛?cè)峄旌弦r砌渠道凍脹破壞進(jìn)行有限元分析，但都是基于冬季停水狀態(tài)下的研究，還有一類(lèi)冬季運(yùn)行的輸水渠道，如電站的引水渠、給城市供水的防滲渠道以及跨流域調(diào)水工程的輸水渠道，此類(lèi)渠道在運(yùn)行時(shí)不僅存在著水面出現(xiàn)冰蓋問(wèn)題，而且還存在著極大的渠道凍脹破壞的問(wèn)題。宋玲[8-9]從溫度含水率二元復(fù)合場(chǎng)的角度，創(chuàng)建了冬季輸水渠道凍脹數(shù)學(xué)模型。對(duì)于冬季輸水渠道而言，水位線以上是凍結(jié)區(qū)，圖1為新疆瑪納斯河四級(jí)電站引水渠上游段橫向裂縫及隆起破壞，可見(jiàn)邊坡板在渠水位附近產(chǎn)生裂縫。

圖1 引水渠道混凝土板局部破壞Fig.1 Local failure of concrete slab in canal

對(duì)于冬季輸水渠道的凍脹問(wèn)題可供參考的文獻(xiàn)較少。本文利用有限元軟件，以新疆瑪納斯河流域紅山嘴電站四級(jí)電站引水渠道為研究對(duì)象，根據(jù)傳熱學(xué)以及凍土力學(xué)提出了在凍土中采用固液相容且互相傳遞熱量的數(shù)學(xué)控制方程并對(duì)襯砌渠道的溫度場(chǎng)、應(yīng)力場(chǎng)進(jìn)行數(shù)值模擬計(jì)算，分析和預(yù)測(cè)在凍結(jié)期渠基土溫度場(chǎng)、應(yīng)力場(chǎng)以及襯砌材料上凍脹力和凍脹變形的變化規(guī)律，旨在為冬季凍區(qū)輸水混凝土渠道的設(shè)計(jì)、正常輸運(yùn)提出參考。

1 溫度含水率二元場(chǎng)

1.1 停水渠道溫度含水率二元耦合場(chǎng)

季節(jié)性?xún)鐾羺^(qū)的運(yùn)河通道，橫截面是一個(gè)復(fù)合耦合溫度場(chǎng)和水分場(chǎng)。

圖2 渠基土橫截面上溫度、含水率二元耦合場(chǎng)中特殊點(diǎn)的溫度、含水率示意圖Fig.2 Abridged general view of compound temperature and moisture content at special points in the binary variable field

如圖2所示，耦合場(chǎng)的臨氣面是一彎曲的凹形線，而原狀渠床的臨氣面則是一條直線。當(dāng)其上空溫度降至0℃以下時(shí)，兩者水熱耦合場(chǎng)是不同的。倘若A點(diǎn)處以下的土體含水率等于或者大于土的臨界凍脹含水率時(shí)，該處附近的土層最終會(huì)產(chǎn)生凍脹現(xiàn)象和凍脹效應(yīng)，但與原狀渠床的凍脹不同。

1.2 輸水襯砌渠道溫度含水率二元耦合場(chǎng)

嚴(yán)寒地區(qū)的輸水襯砌渠道，該橫截面仍然是一個(gè)“耦合”場(chǎng)——水分場(chǎng)、溫度場(chǎng)。存在差異的情況是，在冬季渠道供水狀況下，渠堤、渠坡板，底板下的土體含水率梯度走勢(shì)如圖3。

圖3 渠基土橫截面上溫度、含水率二元耦合場(chǎng)中特殊點(diǎn)的溫度、含水率示意圖Fig.3 Abridged general view of compound temperature and moisture content at special points in the binary variable field

渠坡處C點(diǎn)位置以下的土體含水率遠(yuǎn)遠(yuǎn)大于原始狀態(tài)下的土體含水率，初始滲漏階段的滲漏強(qiáng)度較大，滲水隨著入滲時(shí)間的增長(zhǎng)而減小，到最終達(dá)到穩(wěn)定狀態(tài)。這是因?yàn)樵跐B透初期濕潤(rùn)峰的行進(jìn)距離較短，而水力梯度較大，最后隨著水力梯度急劇減小并在后期逐漸趨于穩(wěn)定造成的。

2 力學(xué)模型的建立

2.1 熱傳導(dǎo)分析

在渠基土凍結(jié)的熱傳導(dǎo)過(guò)程中，由于渠基土凍結(jié)過(guò)程緩慢，且歷時(shí)較長(zhǎng)，故可將其視為穩(wěn)態(tài)傳熱的過(guò)程（忽略相變過(guò)程中釋放的大量熱量）。同時(shí)，將渠道凍脹問(wèn)題簡(jiǎn)化為平面問(wèn)題，二維平面應(yīng)變問(wèn)題基于穩(wěn)定狀態(tài)的接觸熱傳導(dǎo)控制方程為：

式中：T為溫度；λx、λy為凍土沿著x、y向的導(dǎo)熱系數(shù)；A為計(jì)算的凍脹區(qū)域。

對(duì)于凍土，其導(dǎo)熱系數(shù)隨溫度變化不大，而隨土體含水量變化較大，含水量越小導(dǎo)熱系數(shù)越小，熱阻越大[11]。

2.2 凍脹機(jī)理及材料對(duì)流傳熱模型

當(dāng)土體中的細(xì)粒土凍結(jié)時(shí)，凍結(jié)速率很小時(shí)，產(chǎn)生的溫度梯度作用下促進(jìn)水分流動(dòng)。在這種力的作用下，水分不停地涌向凍結(jié)鋒面并在該處堆積，體積膨脹。土體受土的自身構(gòu)成、密實(shí)度、水分的影響。本文不考慮水分遷移這一復(fù)雜的過(guò)程，將凍土看作各向同性完全彈性材料，其彈性模量隨著溫度變化而變化[12]，將此問(wèn)題看作流體外略平壁時(shí)產(chǎn)生對(duì)流傳熱問(wèn)題，如圖4所示。

圖4 冬季輸水渠道流體對(duì)流傳熱示意圖Fig.4 diagram of convection heat transfer in canal with water in winter

考慮物性條件，假設(shè)流體速度均勻，且對(duì)流只發(fā)生在平行壁面方向，在黏性流體的對(duì)流傳熱經(jīng)典理論中，在水流與壁面不發(fā)生滑移的前提下，確定水流與固體壁面之間僅憑借熱傳導(dǎo)產(chǎn)生熱量的互換，因此可以得到對(duì)流傳熱系數(shù)表達(dá)式為：

式（2）中：λ為流體的導(dǎo)熱系數(shù)，tw，t∞為來(lái)流的溫度和壁面的溫度。為貼壁處壁面法線方向上的流體溫度變化率。

3 原型概況

本文仿真計(jì)算原型為新疆紅山嘴水電站四級(jí)引水渠梯形渠斷面，引水渠全長(zhǎng)4.153 km，底寬3.0m，渠深4.5m，混凝級(jí)土強(qiáng)度為C20的襯砌防滲渠道，渠內(nèi)水深3.5m，邊、底板板厚均為0.2 m，渠道縱坡坡度1/1400，邊坡坡度1/1.75。凍土層的最大厚度1.3 m，渠基土與混凝土板的摩擦系數(shù)為0.35。文獻(xiàn)[14]通過(guò)小水壺多次多點(diǎn)取水測(cè)平均溫度的方法，測(cè)得該渠道在冬季運(yùn)行過(guò)程當(dāng)中，沿流方向水溫都在1.5～3.2℃范圍內(nèi)，因此渠內(nèi)流體在通常情況下不會(huì)結(jié)冰。

3.1 參數(shù)計(jì)算

對(duì)于穩(wěn)態(tài)的熱傳導(dǎo)，不同導(dǎo)熱系數(shù)的材料將影響溫度場(chǎng)的分布?；炷恋膶?dǎo)熱系數(shù)λc=1.65W/(m·℃)[15]。依據(jù)熱傳導(dǎo)方程，土壤凍結(jié)時(shí)導(dǎo)熱系數(shù)λf=1.9870W/(m·℃)，幾何模型需要考慮邊界條件等因素建立，并且能簡(jiǎn)化原渠道模型。

水位面溫度設(shè)為常溫1.5℃，假定渠道模型無(wú)限長(zhǎng)，轉(zhuǎn)化為平面單元。模型的下邊界無(wú)位移，x方向?yàn)樽杂蛇吔纾线吔鐑蓚€(gè)方向均為自由邊界。應(yīng)力邊界條件：左右邊界剪應(yīng)力為0，表1是在不同冰體積含量下，含水量一定時(shí)的凍土彈性模量[16]。線膨脹系數(shù)按照η/Tmin取值，η為凍脹率，Tmin為相應(yīng)部位月平均表面溫度最小值。

相應(yīng)材料力學(xué)參數(shù)[17]的選取見(jiàn)表2。

表1 凍結(jié)砂土的等效模量Tab.1 Equivalent modulus of frozen sand

表2 材料力學(xué)參數(shù)Tab.2 Mechanical parameters of materials

3.2 計(jì)算結(jié)果

3.2.1 溫度場(chǎng)

本文對(duì)整個(gè)模型采用三角形網(wǎng)格劃分，共分成12328個(gè)網(wǎng)格單元，6256個(gè)節(jié)點(diǎn)。邊界條件中，下邊界溫度取10℃，上邊界溫度取為原型渠道相應(yīng)部位大氣溫度最小值，左右邊界絕熱。以溫度場(chǎng)為模擬對(duì)象對(duì)冬季停水渠道和冬季輸水渠道模擬結(jié)果進(jìn)行對(duì)比。

圖5為冬季停水、輸水渠道等溫線分布情況，由圖5可知：

（1）渠道表層溫度梯度大，而隨著深度的增加，渠基土所受溫度的影響也就越來(lái)越?。▓D5a）。

（2）冬季輸水渠道凍結(jié)鋒面深度不一，從水位面以上的渠坡到渠底的土體溫度變化快即溫度梯度大，而在遠(yuǎn)離渠坡處渠堤土的溫度分布呈現(xiàn)為幾組近似平行的直線，并且由于渠堤土與部分渠坡板下土體與大氣負(fù)溫直接接觸，均不同程度產(chǎn)生了不均勻凍脹（圖5b）。而與渠道內(nèi)流體直接接觸的渠坡板（水位線下部）和底板下部土體，因?yàn)榱黧w的介入阻隔了這部分土體與大氣負(fù)溫的接觸，并在熱傳遞和襯砌體滲漏造成土體含水率增加的情況下，該部位土體未產(chǎn)生凍結(jié)（圖5b）。

圖5 冬季停水(a)、輸水(b)渠道等溫線分布圖Fig.5 Isothermal distribution map of temperature canal water cut and with water in winter

3.2.2 應(yīng)力場(chǎng)計(jì)算結(jié)果

根據(jù)數(shù)值計(jì)算得出的冬季停水渠道和輸水渠道法向凍脹力分布圖，繪制出沿襯砌板長(zhǎng)的展開(kāi)圖見(jiàn)圖6a。從圖6a可以看出：

（1）冬季停水渠道底板中部受到的法向凍脹力較小，而在渠底板端部受到的負(fù)法向凍脹力較大。輸水渠道凍脹力沿襯砌板分布方面產(chǎn)生了較為明顯的差異，由于正負(fù)溫的過(guò)渡變化，水位面以下渠坡板、底板的法向凍脹力為0，而水位面以上凍脹力與停水渠道類(lèi)似。在負(fù)溫狀態(tài)下，冬季輸水渠道渠坡板的法向凍脹力峰值要比冬季停水渠道渠坡板法向凍脹力峰值高約18%，切向凍脹力大小峰值兩者相似。

（2）冬季停水渠道凍脹破壞往往發(fā)生在邊坡板下1/3范圍內(nèi)，而冬季輸水渠道凍脹破壞發(fā)生在邊坡板中部1/3范圍內(nèi)。這是因?yàn)槎疽浪煌陀谇涝O(shè)計(jì)水位的1/5。此水位對(duì)應(yīng)的渠道水位線在邊坡板中部以下1/3處的范圍內(nèi)波動(dòng)，這個(gè)位置以上的渠坡板、渠坡土處于劇烈凍脹區(qū)，并且極易發(fā)生破壞。

（3）兩模型法向凍脹力在水位面處產(chǎn)生了比較大的變化，在負(fù)溫段，冬季輸水混凝土渠道距離水位面1/3處的坡板處兩側(cè)受拉力，距渠頂1/3處兩側(cè)受壓力，由于混凝土襯砌層受力性質(zhì)抗壓強(qiáng)度大抗拉強(qiáng)度低，故距離水位面1/3-1/2的渠坡板處易發(fā)生凍脹破壞，這也與實(shí)際工程情況相類(lèi)似。

混凝土剛性襯砌板下表面受到的切向凍脹力分布如圖6b所示，其中以指向延渠坡板拔出方向?yàn)檎邢騼雒浟?。由圖6b可知：

（1）冬季停水渠道陰、陽(yáng)坡板以距離渠道底部1/3處分界，凍結(jié)力由分界點(diǎn)沿著混凝土襯砌板長(zhǎng)度方向指向拔出方向，因此該部位極易產(chǎn)生破壞。這一特點(diǎn)被工程實(shí)踐證明。

（2）冬季輸水渠道切向凍脹力以渠水位面上部1/3處分界，更好揭示了經(jīng)歷了一個(gè)寒季之后輸水渠道高位坡板破壞的現(xiàn)象。

4 結(jié)論和討論

本文不僅首次從場(chǎng)的角度分析輸水渠道橫斷面土體溫度-含水率兩場(chǎng)的耦合作用，并與冬季停水渠道耦合場(chǎng)進(jìn)行對(duì)比分析，研究了輸水渠道渠基土凍脹后的兩場(chǎng)變化情況，而且首次考慮固液耦合模型對(duì)渠基土溫度和襯砌層所受凍脹力的影響，并結(jié)合有限元軟件對(duì)冬季輸水渠道進(jìn)行分析，明確了冬季輸水渠道溫度分布情況，研究了襯砌板下凍脹力的作用方式與大小，得到了以下結(jié)論：

（1）冬季輸水渠道中與渠道內(nèi)流體直接接觸的渠坡板（水位線下部）和底板下部土體未產(chǎn)生凍結(jié)。

（2）考慮固液相容的模型比不考慮固液相容模型的凍結(jié)鋒面的位置稍高。

（3）經(jīng)過(guò)一個(gè)冬季之后，輸水渠道發(fā)生凍脹破壞的部位要高于停水渠道。

（4）本文采用有限元軟件對(duì)冬季輸水渠道進(jìn)行分析，結(jié)果基本能反映渠道凍脹變形的受力規(guī)律以及凍脹變形規(guī)律，故可為冬季輸水混凝土襯砌渠道的設(shè)計(jì)與施工提供參考與建議。

（5）本文存在缺少準(zhǔn)確的實(shí)測(cè)數(shù)據(jù)、具體的實(shí)驗(yàn)驗(yàn)證等不足，因此有待對(duì)此類(lèi)渠道開(kāi)展更加深入的研究。

[1]Harlan R L.Analysis of Coupled heat-fluid transport in partially frozen soil[J].Water Research，1973，9(5):1314-1323.

[2]Loch J P G.Thermodynamic Equilibrium Between Ice and Water in Porous Media[J].Soil Science，1978，126(2):77-80.

[3]Taylor G S，Luthin J N.A model for coupled heat and moisture transfer during soil freezing[J].Canadian Geotechnical Journal，1978，15(4):548-555.

[4]Shen M，Ladanyi B.Modelling of coupled heat，moisture and stress field in freezing soil[J].Cold Regions Science&Technology，1987，14(3):237-246.

[5]Konrad J M，Morgenstern N R.Prediction of frost heave in the laboratory during transient freezing[J].Canadian Geotechnical Journal，1982，19(3):250-259.

[6]安維東，陳肖柏，吳紫汪.渠道凍結(jié)時(shí)熱質(zhì)遷移的數(shù)值模擬[J].冰川凍土，1987，29(1)：35-46．An W D，Chen X B，Wu Z W.Numerical simulation of heat and mass transfer during channel freezing[J].Journal of Glaciology and Geocryology，1987，29(1):35-36

[7]董江偉，湯驊，姜海波，等.大斷面剛?cè)峄旌弦r砌渠道凍脹有限元分析[J].人民黃河，2016，38(11):149-152.Dong J W，Tang h，Jiang H B.et al.Finite Element Analysis of Rigid Mixed Lining Channels’Frost Heaving of Large Section[J].Yellow River，2016，38(11):149-152

[8]宋玲，歐陽(yáng)輝，余書(shū)超.混凝土防滲渠道冬季輸水運(yùn)行中凍脹與抗凍脹力驗(yàn)算[J].農(nóng)業(yè)工程學(xué)報(bào)，2015，31（18）:114-120.Song L，Ou Y H，Yu S C.Frozen heaving and capacity of frozen heaving resistance oftrapezoidalconcrete lining canal with water in winter[J].Transactions of the Chinese Society of Agricultural Engineering，2015，31（18）:114-120.

[9]宋玲.季節(jié)凍土區(qū)剛性襯砌渠道溫度、應(yīng)力及凍脹變形場(chǎng)的計(jì)算[J].石河子大學(xué)學(xué)報(bào)(自然科學(xué)版)，2008，26(6):763-767.Song L.Studies on computational method about temperature，stress and frost deformation of rigid lined canal in seasonal frost region[J].Journal of Shihezi University(Natural Science)，2008，26(6):763-767.

[10]康紹忠，粟曉玲，沈清林，等.石羊河流域水資源利用與節(jié)水農(nóng)業(yè)發(fā)展模式的戰(zhàn)略思考[J].水資源與水工程學(xué)報(bào)，2004，15(4):1-8.Kang S Z，Su X L，Shen Q L，et al.Strategy for water resources utilization and water-saving agriculture development models in Shiyanghe river basin[J].Journal of Water Resources&Water Engineering.2004，15(4):1-8.

[11]徐學(xué)祖，王家澄，張立新.凍土物理學(xué)[M].北京：科學(xué)出版社，2010.

[12]陳肖柏，劉建坤，劉鴻緒，等.土的凍結(jié)作用與地基[M].北京：科學(xué)出版社，2006.

[13]張輝，陳善年，徐益謙.理想流體對(duì)流傳熱問(wèn)題的理論解[J].熱科學(xué)與技術(shù)，2004(3):200-204.Zhang H，Chen S N，Xu Y Q.Analytical solutions to convective heat transfer of ideal fluid[J].Journal of Thermal Science and Technology，2004(3):200-204.

[14]黃酒林，宗全利，劉貞姬，等.高寒區(qū)引水渠道抽水融冰原型試驗(yàn)及分析[J].石河子大學(xué)學(xué)報(bào)(自然科學(xué)版)，2014,32（3）:392-396.Huang J L，Zong Q L，Liu Z J，et al.Mechanism and prototype experiment on pumping well water to melt Ice for diversion channel at high altitude and cold regions[J].Jou-rnal of Shihezi University(Natural Science)，2014,32(3):392-396.

[15]劉雄，寧建國(guó)，馬巍.凍土地區(qū)水渠的溫度場(chǎng)和應(yīng)力場(chǎng)數(shù)值分析[J].冰川凍土，2005，27（6）:932-938.Liu X，Ning J G，Ma W.Numerical analyses of the temperature and stress Fields of channel in frozen soil regions[J].Journal of Glaciology and Geocryology，2005，27（6）:932-938.

[16]李爽，王正中，高蘭蘭，等.考慮混凝土襯砌板與凍土接觸非線性的渠道凍脹數(shù)值模擬[J].水利學(xué)報(bào)，2014，45(4):497-503.Li S，Wang Z Z，Gao L L，et al.Numerical simulation of canal frost heaving considering nonlinear contact between concrete lining board and soil[J].Journal of Hydraulic Engineering，2014，45(4):497-503.

[17]寧建國(guó)，朱志武.含損傷的凍土本構(gòu)模型及耦合問(wèn)題數(shù)值分析[J].力學(xué)學(xué)報(bào)，2007，39(1):70-76Ning J G，Zhu Z W.Constitutive model of frozen soil with damage and numerical simulation of the coupled problem[J].Chinese Journal of Theoretical and Applied Mechanics，2007，39(1):70-76.