劉 超，郭妍秀

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

黑龍江省是我國(guó)重要糧食生產(chǎn)基地，2022年產(chǎn)量約占全國(guó)總產(chǎn)量的11%，連續(xù)七年全國(guó)第一。糧食的高產(chǎn)歸因于肥沃的黑土地以及充足的水資源，而水渠作為輸送水分的主要載體成了灌溉農(nóng)田的關(guān)鍵，據(jù)統(tǒng)計(jì)黑龍江省農(nóng)業(yè)用水約占總用水量的4/5，用水效率為25%～40%[1]。同時(shí)，黑龍江省也是典型的季節(jié)性?xún)鐾羺^(qū)，惡劣的氣候環(huán)境導(dǎo)致渠道發(fā)生隆起、滑塌、襯砌開(kāi)裂等凍脹災(zāi)害，據(jù)某大型灌區(qū)的早年調(diào)查報(bào)告顯示，全省約83%渠系工程受到環(huán)境影響發(fā)生凍脹破壞[2]。因此，渠道凍脹的研究成了我國(guó)學(xué)者們的熱門(mén)選題。

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

目前，有限元法在渠道領(lǐng)域已經(jīng)得到了普遍應(yīng)用，但該方法對(duì)設(shè)備、存儲(chǔ)量要求高、無(wú)法解決奇異性問(wèn)題。有限差分法(FDM)妥善地克服了這一缺點(diǎn)，基于能量守恒方程、傅里葉定律等基本理論，利用網(wǎng)格節(jié)點(diǎn)上的函數(shù)值差商代替控制方程中的導(dǎo)數(shù)進(jìn)行離散，建立以網(wǎng)格節(jié)點(diǎn)上的值為因變量的代數(shù)方程組的模型[13-15]。本文以黑龍江省松干渠道為原型，通過(guò)室內(nèi)試驗(yàn)獲取凍脹模擬的關(guān)鍵參數(shù)-熱力學(xué)參數(shù)、凍脹系數(shù)、力學(xué)參數(shù)，采用有限差分方法建立封閉式渠道模型，施加當(dāng)?shù)囟緶囟惹€(xiàn)模擬季凍區(qū)渠基凍脹全過(guò)程，分析渠基的凍脹特性及一般規(guī)律。

1 渠基模型的建立

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

1.1 模型尺寸及網(wǎng)格劃分

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

圖1 渠道模型的網(wǎng)格劃分

1.2 模型監(jiān)測(cè)點(diǎn)布置

為了實(shí)時(shí)監(jiān)測(cè)渠基在凍結(jié)過(guò)程中的溫度場(chǎng)、位移場(chǎng)的變化，分別在渠底、渠坡中間的不同深度處設(shè)置監(jiān)測(cè)點(diǎn)，共設(shè)置24個(gè)，其中監(jiān)測(cè)點(diǎn)1#、監(jiān)測(cè)點(diǎn)11#分別位于渠底、渠坡表面，監(jiān)測(cè)點(diǎn)的具體分布情況見(jiàn)表1、表2。

表1 渠底監(jiān)測(cè)點(diǎn)布置情況 m

表2 渠坡監(jiān)測(cè)布置情況 m

1.3 模型的力學(xué)初始、邊界條件

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

渠基模型的初始應(yīng)力場(chǎng)采用彈塑性變形方法，首先將抗拉強(qiáng)度和黏聚力設(shè)置為最大值以防止模型發(fā)生屈服破壞，待計(jì)算完成后將二者恢復(fù)成初始數(shù)值進(jìn)行運(yùn)算，直至達(dá)到平衡狀態(tài)，渠基模型的初始應(yīng)力場(chǎng)如圖2所示。

圖2 初始應(yīng)力場(chǎng)

1.4 模型參數(shù)

渠道模擬的關(guān)鍵參數(shù)包含：熱力學(xué)參數(shù)、力學(xué)參數(shù)、凍脹系數(shù)，可分別通過(guò)導(dǎo)熱系數(shù)試驗(yàn)、三軸試驗(yàn)、凍脹試驗(yàn)獲取。上述試驗(yàn)皆選取含水率為15%的渠基土作為試驗(yàn)土樣。

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

(1)

式中：η為修正后的凍脹系數(shù)，℃-1；η0為初始凍脹系數(shù)，℃-1；Ls為試驗(yàn)測(cè)得的凍脹增量，mm；L為模擬所得凍脹量，mm。

查閱文獻(xiàn)得到襯砌板的力學(xué)參數(shù)、熱力學(xué)參數(shù)及凍脹系數(shù)[16]，結(jié)合上述試驗(yàn)所得數(shù)據(jù)匯總成如表3所示的模型參數(shù)。

表3 模型參數(shù)r

1.5 渠基模型的溫度場(chǎng)施加

溫度場(chǎng)演變實(shí)質(zhì)是土質(zhì)溫度分布不均勻?qū)е碌臒崃總鬏斶^(guò)程。本次模擬采用顯式的求解算法，待初始靜力達(dá)到平衡，對(duì)渠基表面施加溫度，運(yùn)用各向同性傳導(dǎo)方式模擬季凍區(qū)冬季渠基土凍脹全過(guò)程。應(yīng)用FLAC3D進(jìn)行凍脹模擬計(jì)算時(shí)，作出以下假設(shè)：

(1)土體溫度場(chǎng)的熱應(yīng)力分布僅與溫度有關(guān)。

(2)流體的分布不受結(jié)構(gòu)應(yīng)力的影響。

(3)可將土體凍結(jié)過(guò)程中的某一單位時(shí)間段內(nèi)的凍脹率、導(dǎo)熱系數(shù)視為常量。

1.6 氣溫變化的選取

由于黑龍江省從11月份開(kāi)始發(fā)生凍結(jié)，選取大慶市肇源縣2019年11月1日—2020年3月19日的氣溫作為凍脹模擬時(shí)間段。對(duì)氣溫曲線(xiàn)進(jìn)行離散化處理，以每連續(xù)10 d的平均氣溫作為一組計(jì)算氣溫，共14組。該地區(qū)日平均氣溫和每組的計(jì)算氣溫如圖3所示。

圖3 模擬時(shí)間段的每組計(jì)算氣溫

1.7 溫度邊界、初始條件

本文所模擬的渠道上邊界采用對(duì)流熱邊界條件，即流體溫度隨時(shí)間的變化情況。根據(jù)黑龍江水利科學(xué)研究院現(xiàn)場(chǎng)實(shí)際監(jiān)測(cè)數(shù)據(jù)測(cè)得凍結(jié)前的渠基表面溫度為3.1 ℃，開(kāi)始凍結(jié)后的溫度按圖3的每組計(jì)算氣溫施加；渠基底部采用溫度邊界條件，施加恒定溫度12 ℃保持不變；渠基側(cè)面采用熱流邊界條件，熱流密度為零。

2 渠基冬季凍脹模擬結(jié)果分析

利用構(gòu)建的數(shù)值模型，模擬季凍區(qū)渠基在冬季凍脹的全過(guò)程，監(jiān)測(cè)并分析渠底與渠坡的溫度場(chǎng)、位移場(chǎng)變化情況。

2.1 凍結(jié)全過(guò)程中溫度場(chǎng)變化規(guī)律

從渠基在整個(gè)冬季的溫度場(chǎng)變化云圖可以發(fā)現(xiàn)：渠坡和渠底表層土溫度變化顯著，溫度梯度相對(duì)較大，渠基深處溫度基本不受氣溫的影響；渠底下部的土層溫度梯度近似平行。

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

圖4 渠基溫度變化曲線(xiàn)

2.2 凍結(jié)過(guò)程中凍脹量變化規(guī)律

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

圖5 渠基凍脹量隨時(shí)間的變化情況

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

3 結(jié) 論

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