季凍區渠基凍脹模擬分析

劉 超，郭妍秀

(1.黑龍江省江河流域保護中心，黑龍江 哈爾濱 150069；2.黑龍江省水利科學研究院，黑龍江 哈爾濱 150080)

黑龍江省是我國重要糧食生產基地，2022年產量約占全國總產量的11%，連續七年全國第一。糧食的高產歸因于肥沃的黑土地以及充足的水資源，而水渠作為輸送水分的主要載體成了灌溉農田的關鍵，據統計黑龍江省農業用水約占總用水量的4/5，用水效率為25%～40%[1]。同時，黑龍江省也是典型的季節性凍土區，惡劣的氣候環境導致渠道發生隆起、滑塌、襯砌開裂等凍脹災害，據某大型灌區的早年調查報告顯示，全省約83%渠系工程受到環境影響發生凍脹破壞[2]。因此，渠道凍脹的研究成了我國學者們的熱門選題。

渠道凍脹是渠基土在負溫作用下形成溫度梯度，導致土內水分向凍結鋒面遷移形成冰透鏡體，渠基向上隆起的現象。渠道的研究方法包括現場監測、室內試驗研究和數值模擬分析，傳統的現場監測及室內試驗用時長、資源耗量大，數值模擬可有效降低時間和成本，已成為一種普遍的研究手段。1973年，RL Harlan[3]首次提出水熱耦合的概念，并給出了相應的數學模型，為凍脹的數值模擬分析奠定了基礎。土體凍脹的主要影響因素包含土質、水分、溫度[4-6]，只有當三者滿足一定要求時，才會發生凍脹，因此數值模擬也圍繞著他們展開研究。劉旭東等[7]基于有限元方法探究渠基凍脹敏感性，發現影響凍脹的敏感性因素可依次排序為：渠道斷面形式>渠基土凍脹系數>分縫位置>混凝土襯砌厚度>混凝土襯砌板彈性模量。李爽[8]考慮了混凝土襯砌的彈塑性以及凍土、襯砌接觸的非線性，對襯砌與渠基土之間滑移和脫落的破壞過程進行數值模擬，研究表明：考慮非線性接觸的模擬更加真實可靠，法向凍脹力和切向凍脹力較未考慮非線性接觸的模型減小一半。杜民瑞[9]基于流傳熱公式分別對冬季停水、輸水渠道進行溫度場、應力場分析，發現輸水渠道的邊坡法向應力大于停水渠道。陳瑞考[10]基于水熱力耦合方程以及有限元方法，建立關于THM的三場耦合渠道模型，針對梯形渠道建立水熱耦合模型，分析冬季輸水期間渠基內部的水分對凍脹的影響，探討輸水與停水渠道之間的差異。張海晨[11]在考慮高地下水位的基礎上，利用有限元軟件建立混凝土襯砌模型，分析渠基不同位置處的凍脹變形情況，發現渠坡的凍脹變形較渠底更顯著，且最大凍脹量出現在渠坡1/3處。王羿等[12]考慮氣溫、渠基水分場、地下含水邊界和土體滲流特性，構建了關于渠槽形狀耦合的渠基土不均勻凍脹的水-熱-力耦合數學模型。

目前，有限元法在渠道領域已經得到了普遍應用，但該方法對設備、存儲量要求高、無法解決奇異性問題。有限差分法(FDM)妥善地克服了這一缺點，基于能量守恒方程、傅里葉定律等基本理論，利用網格節點上的函數值差商代替控制方程中的導數進行離散，建立以網格節點上的值為因變量的代數方程組的模型[13-15]。本文以黑龍江省松干渠道為原型，通過室內試驗獲取凍脹模擬的關鍵參數-熱力學參數、凍脹系數、力學參數，采用有限差分方法建立封閉式渠道模型，施加當地冬季溫度曲線模擬季凍區渠基凍脹全過程，分析渠基的凍脹特性及一般規律。

1 渠基模型的建立

基于能量平衡方程和傅里葉熱傳導定律的計算原理，利用有限差分法建立松干渠道渠基模型，模擬冬季停水渠道的凍脹全過程。

1.1 模型尺寸及網格劃分

已知該渠基總尺寸為20.00 m×20.00 m×10.15 m，渠頂寬2 m，渠底寬4 m，橫斷面渠槽深4 m，坡比為1∶1.5，渠底基土厚6 m，縱向總長度為20 m，基土上部為0.15 m厚的襯砌板。將渠基的X軸、Y軸、Z軸分別劃分成100、100、51個網格，共510 000個網格，如圖1所示。

圖1 渠道模型的網格劃分

1.2 模型監測點布置

為了實時監測渠基在凍結過程中的溫度場、位移場的變化，分別在渠底、渠坡中間的不同深度處設置監測點，共設置24個，其中監測點1#、監測點11#分別位于渠底、渠坡表面，監測點的具體分布情況見表1、表2。

表1 渠底監測點布置情況 m

表2 渠坡監測布置情況 m

1.3 模型的力學初始、邊界條件

由于該模型渠道冬季停水，忽略雪荷載等因素的影響，只考慮渠道本身自重，渠基底部和側向受到來自周邊土的擠壓，因此模型的側向及底部分別受到水平約束和豎向約束限制其位移，上部為自由表面。

渠基模型的初始應力場采用彈塑性變形方法，首先將抗拉強度和黏聚力設置為最大值以防止模型發生屈服破壞，待計算完成后將二者恢復成初始數值進行運算，直至達到平衡狀態，渠基模型的初始應力場如圖2所示。

圖2 初始應力場

1.4 模型參數

渠道模擬的關鍵參數包含：熱力學參數、力學參數、凍脹系數，可分別通過導熱系數試驗、三軸試驗、凍脹試驗獲取。上述試驗皆選取含水率為15%的渠基土作為試驗土樣。

熱力學參數依托于導熱系數測定儀，試樣制備完成后，將導熱系數測定儀插入試樣中部并設置不同溫度。力學參數的獲取方法是將含水率為15%的土樣制備成Φ100 mm×H50 mm的試樣，利用GDS三軸儀進行三軸UU試驗(UU表示不固結不排水)。

(1)

式中：η為修正后的凍脹系數，℃-1；η0為初始凍脹系數，℃-1；Ls為試驗測得的凍脹增量，mm；L為模擬所得凍脹量，mm。

查閱文獻得到襯砌板的力學參數、熱力學參數及凍脹系數[16]，結合上述試驗所得數據匯總成如表3所示的模型參數。

表3 模型參數r

1.5 渠基模型的溫度場施加

溫度場演變實質是土質溫度分布不均勻導致的熱量傳輸過程。本次模擬采用顯式的求解算法，待初始靜力達到平衡，對渠基表面施加溫度，運用各向同性傳導方式模擬季凍區冬季渠基土凍脹全過程。應用FLAC3D進行凍脹模擬計算時，作出以下假設：

(1)土體溫度場的熱應力分布僅與溫度有關。

(2)流體的分布不受結構應力的影響。

(3)可將土體凍結過程中的某一單位時間段內的凍脹率、導熱系數視為常量。

1.6 氣溫變化的選取

由于黑龍江省從11月份開始發生凍結，選取大慶市肇源縣2019年11月1日—2020年3月19日的氣溫作為凍脹模擬時間段。對氣溫曲線進行離散化處理，以每連續10 d的平均氣溫作為一組計算氣溫，共14組。該地區日平均氣溫和每組的計算氣溫如圖3所示。

圖3 模擬時間段的每組計算氣溫

1.7 溫度邊界、初始條件

本文所模擬的渠道上邊界采用對流熱邊界條件，即流體溫度隨時間的變化情況。根據黑龍江水利科學研究院現場實際監測數據測得凍結前的渠基表面溫度為3.1 ℃，開始凍結后的溫度按圖3的每組計算氣溫施加；渠基底部采用溫度邊界條件，施加恒定溫度12 ℃保持不變；渠基側面采用熱流邊界條件，熱流密度為零。

2 渠基冬季凍脹模擬結果分析

利用構建的數值模型，模擬季凍區渠基在冬季凍脹的全過程，監測并分析渠底與渠坡的溫度場、位移場變化情況。

2.1 凍結全過程中溫度場變化規律

從渠基在整個冬季的溫度場變化云圖可以發現：渠坡和渠底表層土溫度變化顯著，溫度梯度相對較大，渠基深處溫度基本不受氣溫的影響；渠底下部的土層溫度梯度近似平行。

根據溫度場云圖繪制出如圖4所示的渠基溫度場變化曲線。監測點1#、監測點11#表示渠基底部及坡面的最低溫度，分別為-13.84 ℃、-14.80 ℃。渠基溫度與深度呈負相關，0～60 d期間，監測點2#～4#、12#～14#的溫度隨著氣溫的下降快速降低，其余土層溫度降低速率緩慢；61～100 d期間，基土溫度隨著氣溫的穩定而緩慢降低；100 d后氣溫回升，距離表面近的土層受到影響，溫度有所升高，但依舊處于負溫狀態，遠距離土層不受影響。

圖4 渠基溫度變化曲線

2.2 凍結過程中凍脹量變化規律

圖5為渠底與渠坡的凍脹量隨時間的變化情況，其中監測點1#、監測點11#分別位于渠底和渠坡表面，受溫度影響效果最顯著，監測結果代表該位置處的最終變形量，分別為2.83 cm、3.07 cm，遠離渠基表面的監測點5#～10#、16#～24#，溫度未達到初始凍結溫度，不發生凍脹。

圖5 渠基凍脹量隨時間的變化情況

由圖5可知，渠底與渠坡的凍脹變形趨勢基本一致，前10 d渠基不發生變形，凍脹量為0；11～60 d期間，表面負溫下移形成溫度梯度，導致凍脹量快速增加，此時位于渠基表面的監測點1#和監測點11#的凍脹量增長最迅速，分別為1.94 cm、2.29 cm；2#～4#、12#～15#監測點由于受到土壤熱阻的作用，凍脹量的增長速率減緩。61～100 d期間，隨著溫度梯度的減小，凍脹增長速率逐漸變得緩慢，直至增長至最大凍脹量，渠底與渠坡的最大凍脹量分別為2.83 cm、3.07 cm；之后氣溫回升但最高溫度仍<0 ℃，未達到渠基融化溫度，凍脹量趨于穩定。

3 結 論

本文通過室內試驗獲取凍脹模擬關鍵參數，基于有限差分法建立渠基模型，采用溫度施加方法模擬松干渠道冬季凍脹全過程，探究渠基溫度場、位移場演變規律。研究發現：渠基溫度與渠基深度呈負相關，凍脹量變化包含快速凍結階段、緩慢凍結階段、穩定階段，2月份凍脹量達到最大，渠底與渠坡的最大凍脹量分別為2.83 cm、3.07 cm。該方法極大降低了科研成本，為渠基凍脹研究提供了一定參考。