寒旱區梯形渠道襯砌凍脹破壞數值模擬

張海晨　李振　張茂林　姚志鵬　朱興林　王銘巖

摘要：針對剛性襯砌渠道凍脹破壞問題，選取寧夏鹽環定揚黃灌區襯砌渠道為研究對象，在考慮相變及水分遷移模型基礎上，考慮高地下水位對凍脹的影響，建立考慮水分遷移及地下水影響的數學物理模型，并采用Comsol Mutiphysics有限元軟件對梯形渠道混凝土襯砌進行凍脹數值模擬，分別從溫度場、位移場和應力場對剛性襯砌渠道凍脹破壞等方面進行了研究分析。結果表明：陰坡、陽坡凍脹量較大，渠底較小，渠底中部、渠坡1/3坡板長度處凍脹量分別達到最大值。渠底中部的法向凍脹力較小，兩邊逐漸增大，陰坡、陽坡的法向凍脹力分布均勻;渠底的切向凍結力成線性分布。陰坡、陽坡靠近坡腳處切向凍結力較大，左、右兩端坡頂處切向凍結力較小，渠底襯砌板和兩側襯砌板屬于壓彎組合變形構件。

關鍵詞：渠道;襯砌;數值模擬;凍脹破壞;寒旱區

中圖分類號：TV67

文獻標志碼：A

doi：10.3969/j .issn. 1000- 1379.2019.03.032

渠系工程在灌區工程中擔任著主要的輸水配水任務，但由于設計標準低，工程配套差，老化失修嚴重等問題，每年渠道滲漏導致的灌溉水損失占農業用水量的45%，占全國總用水量的33%.嚴重制約渠系工程的安全高效運行[1]。襯砌渠道防滲效果明顯，其整體可減少滲漏損失的70% - 90%.是渠系工程中主要的防滲措施。我國凍土面積廣闊，在北方灌區，每年凍脹融沉導致的渠道襯砌破壞現象屢見不鮮。襯砌渠道凍脹是影響渠系工程安全運行的一大因素，受到國內外相關科研工作者的廣泛關注。Everett[2]通過試驗探索，提出固、液表面的兩相壓力差是造成水分凝冰并引起土體凍脹的主要因素，即毛細理論：Miller[3]最早提出凍結緣的概念，解釋了冰透鏡體的形成過程，即第二凍脹理論：徐學祖等[4]通過研究凍土與融土的熱力學特性及凍脹特性，得出了溫度梯度與水分遷移的關系。這些學者的研究成果為渠道襯砌的凍脹研究提供了有益參考和理論支撐，鑒于試驗統計的差異性及凍土凍脹實際情況的復雜性，使得研究結果存在一定的局限性，因此開展渠道襯砌凍脹的數值模擬分析就顯得十分必要。Harlan[5]基于非飽和凍土的水分遷移理論，最早提出了水一熱耦合模型：Konrad[6]根據凍土中自由水含水量與溫度之間的關系，采用數值方法對凍土中的參數進行預報。然而，這些模型對于高地下水位對凍脹影響的考慮卻較少[7]。鑒于此，本文針對渠道襯砌凍脹破壞現象，選取寧夏鹽環定揚黃灌區襯砌渠道為研究對象，在考慮相變及水分遷移模型的基礎上，考慮高地下水位對凍脹的影響，建立考慮水分遷移以及地下水影響的數學物理模型，并運用Comsol Mu-tiphysics有限元軟件進行模擬，得到溫度場、位移場、應力場的分布，并將數值模擬結果與試驗結果進行對比，驗證數值模擬的合理性與正確性，以期為相似環境的渠道襯砌凍脹破壞治理與防護提供理論依據。

1 數學物理模型的建立

1.1 建立模型的基本假設

渠道防滲襯砌工程是否產生凍脹破壞，與當地的土質、土的含水量、負溫度及工程結構等因素有關[1]。渠基土凍結時，土體、水和冰之間的相互作用的微觀結構及動態過程相當復雜，對于工程設計與優化計算來說，準確模擬整個凍脹過程是非常復雜且沒有必要的[8]。為了便于實際應用，需對其進行適當簡化，主要假設如下[9]。

（1）根據現場及室內試驗研究，假設凍土是均勻連續各向同性體。

（2）水分對凍脹的影響表現在兩方面：其一是土壤水分凍結釋放相變潛熱，影響渠基土溫度場及凍深，從而影響凍脹量：其二是水分凍結形成分凝冰，體積增加，特別是地下水不斷遷移結冰引起的膨脹。針對第一方面，建立考慮水分遷移熱效應的凍土凍脹模型[10];對第二方面，參考《渠系工程抗凍脹設計規范》[11]中考慮地下水遷移影響的凍深一凍脹量關系，對凍土本構關系進行修正。

（3）根據試驗研究，假定相變溫度在同一種土中和同種外力條件下為常值，即暫取相變溫度為0℃[12]。

（4）渠道為線狀結構，沿水流方向長度遠大于垂直水流方向寬度，可以作為二維平面應變問題進行建模。

（5）不考慮孔隙水壓的影響。

1.2 考慮相變潛熱的熱傳導方程

大量的試驗結果證明土壤中未凍水在負溫度梯度作用下由未凍區向凍結鋒面遷移、集聚并凍結釋放潛熱[13]。本文將這部分相變潛熱作為渠基土凍土物理場的熱源看待，修正后的熱傳導公式為

1.3 考慮地下水影響的凍土本構方程

凍土在凍結過程中水凍結成冰，除原位水凍結膨脹外，還有從未凍結區向凍結區遷移的水分凍結成冰。工程中一般將原位水以及遷移水凍結成冰的體積膨脹以凍脹率η表示，其取值不僅與土質、凍深相關，也與地下水埋深有關[14]，《渠系工程抗凍脹設計規范》給出了地下水埋深與凍脹量、凍深的關系，見圖1。

渠基土凍脹是由土壤孔隙中水分凍結成冰膨脹引起的，為便于工程應用，將渠基土視為冷脹熱縮材料。平面應變問題中考慮地下水位的凍土本構方程為

以二維非穩態熱傳導方程建立的溫度場模型與以本構方程建立的力學模型耦合，最終建立熱力耦合分析模型。

2 渠道凍脹數值模擬模型及參數選取

2.1 渠道概況

該工程水源為黃河水，工程位于寧夏回族自治區鹽環定揚黃灌區共用工程八干渠上段，為預制混凝土方形板襯砌梯形渠道，渠道斷面見圖2。邊坡系數m=1.5，襯砌板厚度為6 cm，渠基土質為砂壤土，地下水埋深為渠道底板以下20 cm。工程所在區域屬中溫帶干旱區，多年平均降水量290 mm.多年平均蒸發能力1 340 mm;l月（最冷月）平均氣溫-8.9℃，日照時數長，全年日照時數2 867.9 h。對該渠道進行凍脹數值仿真計算，并與現場實測數據對比。

2.2 有限元模型

為了驗證上述數理方程的合理性與實用性，運用數學物理方程并采用Comsol Mutiphysics有限元軟件模擬寧夏鹽環定揚黃灌區渠道的凍脹過程，有限元模型是在原型渠道基礎上的簡化，基礎從底板向下取250 cm.左、右邊界各取75 cm。在邊界條件中，上邊界溫度取原型渠道相應部位表面溫度，陰坡、陽坡、渠底分別取-4.92、-4.56、-3.55℃，下邊界溫度取10℃，左、右邊界近似隔熱，其位移約束條件是x與z方向位移為零;渠底基土下邊界加y向約束。梯形渠道有限元網格見圖3。

2.3 參數選取

將混凝土襯砌板看成各向同性材料，彈性模量取2.4x104 MPa，將凍土與非凍土視為各向同性彈性體。未凍土彈性模量取15 MPa，凍土的彈性模量隨著溫度的改變而改變，取值見表1。其他材料計算參數及凍土的凍脹率與地下水埋深數據見表2、表3。

3 結果分析

3.1 溫度場分析

由已構建模型對襯砌渠道的凍土溫度場進行模擬，結果見圖4。

由圖4可知，接近于渠底的溫度分布接近于一組平行的直線：0℃等溫線以上，未凍水凍結導致的劇烈相變使等溫線較密集：0℃等溫線以下，等溫線較稀疏，溫差較大，容易造成渠基土的不均勻凍脹，陰坡凍深為98 cm，渠底凍深為55 cm，陽坡的凍深為91 cm。渠坡及渠底表層的溫度變化快，即溫度梯度大：而在離渠坡較遠的深部，溫度分布呈一組幾乎平行的直線，基本上不受渠坡邊界溫度的影響。從溫度場的分布和變化來看，三維邊界、界面幾何形狀以及基土中熱流對溫度分布均有顯著影響。從整體計算結果與前人的研究結果比較來看，模擬溫度場與實測結果基本一致[15]，即上部凍深較大、底部凍深較小，陰坡凍深較大、陽坡凍深較小。土的最大凍脹率出現在地表下約1/3凍結深度處，而最大凍脹量發生在地表層。一般把凍脹量占總凍脹量的70% - 90%部分土體稱為主凍脹區。主凍脹區一般分布于地表。八干渠渠道不同朝向的渠坡接受太陽輻射的能力不同，導致同一地區渠道陰、陽坡平均地溫有一定的差異，同一區域內相同下墊面條件的太陽福射差異主要是渠道（地形因子）中的坡度和坡向不同引起的，總的來說陽坡的年平均地溫高于陰坡的年平均地溫，所以凍土往往具有明顯的非對稱性。

當冷卻強度很大導致土內溫度梯度很大時，表層凍結面迅速向下推移，土中水分來不及從下臥層向凍結面轉移，就在原地凍結成冰。在這種條件下形成的冰，一般均勻散布于土孔隙中或土粒接觸處。此類凍土一般無明顯凍脹。當冷卻強度小時，凍結面向下推移慢，甚至會因水結冰時放出的潛熱的阻擋而較長時間停留在某一深度。此時，下臥層的遷移水流在克服了沿途的阻力后，還來得及到達凍結面形成冰層。有外部水源補給時，凍結面停留的時間越長，則形成的冰層越厚。這種凍土有明顯的凍脹。

3.2 位移場分析

由已構建模型對襯砌渠道的凍土位移進行計算，并與試驗結果進行比對.結果見圖5。

由圖5可知：

（1）模擬結果與模型試驗結果分布規律一致，陰坡、陽坡凍脹量較大，渠底較小。原因是陰、陽坡凍深較大，發生凍脹的土體較多。因坡腳處襯砌板對基土的約束，故臨近坡腳處的凍脹量很小。溫度逐漸降低，土體中未凍結的水在負溫度梯度的動力驅使下向凍結鋒面遷移并最終凍結成冰，加上土體中原位水的凍結，造成土體的體積膨脹以及凍脹量變大，并最終造成了襯砌板的法向位移逐漸增大。

（2）在渠底中部、渠坡1/3長度處凍脹量分別達到最大值，其中：渠底2.16 cm，陰坡4.76 cm.陽坡5.01cm，與李安國模型計算結果基本一致[15]。陽坡凍脹量稍大于陰坡的，原因是在低溫條件下，溫度越低，凍結向下推進越快，水分越難向凍結邊緣附近遷移，致使凍結區由水分遷移形成的冰的含量較少，因此陽坡凍脹量大于陰坡的。

3.3 應力場分析

3.3.1 法向凍脹力分析

渠基土負溫條件下凍結膨脹會對渠道襯砌產生法向凍脹力以及切向凍結力。取垂直于襯砌板下表面的應力張量可得到法向凍脹力分布情況，見圖6。

渠底兩端由于襯砌板對于基土的約束以及陰坡、陽坡的襯砌板對渠底襯砌板的約束作用，渠道底部兩端的法向凍脹力較大。渠底中部的法向凍脹力較小，兩邊的逐漸增大，其原因是基土與襯砌板底部的約束得到釋放;陰坡、陽坡的法向凍脹力分布均勻，但在渠坡靠近坡腳處法向凍脹力較大。

3.3.2 切向凍結力分布

根據沿襯砌板下表面的應力張量可得到切向凍結力分布，見圖7。

由圖7可知，渠底的切向凍結力分布受兩端約束，兩端切向凍結力較大。渠底左側受拉，右側受壓，因此渠底的切向凍結力成線性分布;陰坡、陽坡靠近坡腳切向凍結力較大，左、右兩端坡頂處切向凍結力較小。

4 結論

（1）在考慮相變及水分遷移模型基礎上，考慮高地下水位對凍脹的影響，利用Comsol Mutiphysics軟件將凍土與渠道襯砌體作為一個整體進行熱力耦合數值仿真計算，計算結果與試驗結果基本一致，揭示了剛性襯砌渠道凍脹變形的基本規律及受力規律，可為渠道凍脹破壞研究與治理提供參考。

（2）溫度場數值分析表明渠底的溫度分布接近于一組平行的直線。0℃等溫線以上，等溫線較密集;0℃等溫線以下，等溫線較稀疏：渠坡及渠底表層的溫度變化大，溫度梯度大：深部溫度分布不受渠坡邊界溫度的影響。上部凍深較大，底部凍深較小;陰坡凍深較大，陽坡凍深較小。

（3）位移場數值分析表明陰坡、陽坡凍脹量較大，渠底較小。渠底中部、渠坡1/3坡板長度處凍脹量分別達到最大值，陽坡凍脹量稍大于陰坡凍脹量，襯砌板凍脹量總體分布不均勻。

（4）應力場數值分析表明渠道底部兩端的法向凍脹力較大。渠底中部的法向凍脹力較小，向兩邊逐漸增大，陰坡、陽坡的法向凍脹力分布均勻，都表現在渠坡靠近坡腳處法向凍脹力較大。渠底的切向凍結力成線性分布，左側受拉，右側受壓，坡腳兩端切向凍結力較大，陰坡、陽坡靠近坡腳處切向凍結力較大，左、右兩端坡頂處切向凍結力較小。

參考文獻：

[1] 王文杰，凍土水熱力三場耦合的襯砌渠道凍脹數值模擬研究[D].楊凌：西北農林科技大學，2013：1-5.

[2] EVERETT D H.The Thermodynamics of Frost Damage toPorous Solids[J].Transactions of Faraday Society， 1961， 57（5）：1541-1551.

[3]MILLER R D.Freezing and Heaving of Saturated and Unsat-urated Soils[C]//Highway Research Record. Washington：Highway Research Board， 1972：1-11.

[4]徐學祖，邱維林，溫度梯度誘導薄膜水遷移的凍脹機理[J].科學通報，1997，42（9）：956-959.

[5]HARLAN R L.Analysis of Coupled Heat-Fluid Transport inPartially Frozen Soil[J].Water Resources Research， 1973，9（5）：1314-1323.

[6]KONRAD J M.The Influence of Heat Extraction Rate in Freez-ingsoils[J].Cold Regions Science&Technology，1987， 14（2）：129-137.

[7] 劉月，王正中，王羿，等，考慮水分遷移及相變對溫度場影響的渠道凍脹模型[J].農業工程學報，2016，32（17）：83-88.

[8] 朱志武，寧建國，馬巍，土體凍融過程中水、熱、力三場耦合本構問題及數值分析[J].工程力學，2007，24（5）：138-144.

[9] 陳立杰，王正中，劉旭東，等，高地下水位灌排渠道襯砌結構抗凍脹數值模擬[J].長江科學院院報，2009，26（9）：66-70.

[10] 劉月，考慮遷移水熱效應的土體凍脹模型及雙膜渠道襯砌防凍脹研究[D].楊凌：西北農林科技大學，2016：21-26.

[11] 中華人民共和國水利部，渠系工程抗凍脹設計規范：SL23-2006 [S].北京：中國水利水電出版社，2007：5-13.

[12]柯潔銘，楊平，凍土凍脹融沉的研究進展[J].南京林業大學學報（自然科學版），2004，28（4）：105-108.

[13]曾桂軍，正凍土水分遷移及凍脹模型研究[D].北京：中國科學院大學，2015：28-31.

[14] 郭瑞，考慮分凝冰效應的飽和土凍脹模型及“適應型”渠道襯砌結構研究[D].楊凌：西北農林科技大學，2016：17-23.

[15] 李安國，陳瑞杰，渠道凍脹模擬試驗及襯砌結構受力分析[J].水利與建筑工程學報，2000，6（1）：5-16.