凍結(jié)條件下凍土物理力學(xué)特性的數(shù)值模擬

凍結(jié)條件下凍土物理力學(xué)特性的數(shù)值模擬

中國(guó)是世界上最大的冷凍土壤分布區(qū)域。多年凍土和季節(jié)凍土的面積分別占全國(guó)面積的21.5%和53.5%。53.5%.中國(guó)位于中國(guó)東部30n的以北，川陜地區(qū)33n的以北。土壤的凍結(jié)是一個(gè)十分復(fù)雜的過(guò)程,在凍結(jié)過(guò)程中,土壤中水分狀況、熱狀況和溶質(zhì)狀況的變化并不是獨(dú)立的,而是相互制約、相互影響的,是一個(gè)耦合變化過(guò)程。綜合性的描述和研究土壤中水、熱、溶質(zhì)耦合運(yùn)移的文獻(xiàn)較為少見(jiàn)。20世紀(jì)90年代有些學(xué)者從不同側(cè)面對(duì)土壤中水、熱、鹽運(yùn)移問(wèn)題進(jìn)行了理論及試驗(yàn)研究,Nassar利用土壤中水、熱、鹽運(yùn)動(dòng)方程和連續(xù)方程,提出了描述土壤中水、熱鹽運(yùn)動(dòng)的耦合方程。本文通過(guò)對(duì)棕壤土在凍結(jié)條件下水、熱、溶質(zhì)耦合運(yùn)移的模擬,得出水、熱、溶質(zhì)遷移規(guī)律,為后續(xù)研究含多種溶質(zhì)的水、熱、溶質(zhì)耦合運(yùn)移以及田間水、熱、溶質(zhì)耦合運(yùn)移奠定基礎(chǔ)。1凍結(jié)時(shí)間和溫度室內(nèi)試驗(yàn)在中科院寒區(qū)旱區(qū)環(huán)境與工程研究所進(jìn)行,所采用的土壤是棕壤土,其主要物理指標(biāo)見(jiàn)表1,加入5%氯化鉀溶液,土樣高度和直徑均為10cm,從上端降溫進(jìn)行凍結(jié),降溫速率為-0.1℃/h,凍結(jié)時(shí)間為120h,底板控制溫度在1℃,土樣上、下兩端密封。試驗(yàn)結(jié)束后分層取樣測(cè)定含水率的分布狀況。實(shí)驗(yàn)土柱初始含水率為40%,初始溶液濃度為5%,初始溫度分布如表2,土壤凍結(jié)特征曲線如圖1和圖2。2數(shù)學(xué)模型2.1一維和一維熱流運(yùn)移基本方程假定土壤介質(zhì)不可壓縮、均質(zhì)且各向同性;凍結(jié)過(guò)程中,水分以液態(tài)形式遷移;忽略水汽遷移和熱量對(duì)流作用的影響及溶質(zhì)勢(shì)梯度對(duì)水分遷移的影響;凍土中未凍水含量與土壤負(fù)溫處于動(dòng)態(tài)平衡中;水、熱、溶質(zhì)的運(yùn)動(dòng)主要發(fā)生在垂向,可近似為一維問(wèn)題。一維情況下,根據(jù)Darcy定律和質(zhì)量守恒原理,得非飽和土壤水分運(yùn)動(dòng)基本方程:?θ?t=??z[D(θ)?θ?z]-?Κ(θ)?z-ρΙρw?θΙ?t(1)式中θ為液態(tài)水含水率;θI為含冰率;ρI為冰密度;ρw為水密度;D為土壤水分?jǐn)U散系數(shù);K為土壤導(dǎo)水系數(shù)。一維情況下,根據(jù)Fourier定律和質(zhì)量守恒原理,得熱流基本方程:cv?Τ?t=??z[Κh?Τ?z]+LΙρΙ?θΙ?t(2)式中cv為土壤的熱容量;T為溫度;Kh為土壤熱傳導(dǎo)系數(shù);LI為冰的溶解潛熱。一維情況下,根據(jù)Fick定律和質(zhì)量守恒原理,得溶質(zhì)運(yùn)移基本方程:?(θc)?t=??z[Dsh(ν,θ)?c?z]-?(qc)?z(3)式中c為溶質(zhì)濃度;Dsh為水動(dòng)力彌散系數(shù)。2.2土體含水率的確定凍土中,未凍水含量與土體負(fù)溫保持動(dòng)態(tài)平衡,即土體中未凍水含水率和溫度的關(guān)系:θ≤θmax(Τ)(4)式中θmax為相應(yīng)土體負(fù)溫下可能的最大未凍水含量(土壤凍結(jié)特征曲線)。3熱流差分方程凍結(jié)土柱的水、熱、溶質(zhì)運(yùn)動(dòng)狀況是由土柱兩端的土壤水、熱、溶質(zhì)條件決定的。由實(shí)驗(yàn)可知,熱流差分方程上、下邊界條件均為已知,屬于第1類邊界條件。水流差分方程和溶質(zhì)差分方程上下邊界為零通量面,屬于第2類邊界條件。3.1土柱冷端差分方程水流差分方程計(jì)算時(shí),采用半控制容積法處理,即原差分網(wǎng)格不變,對(duì)邊界處的半個(gè)單元,直接根據(jù)水量平衡列出差分方程。如圖3,對(duì)土柱冷端i為0～0.5處列平衡方程:(R-q1+1/2)△t=(ρΙρw△θΙ+△θu)△z12(5)式中表面通量R=0。q1+1/2=-D(θu)?θu?z+Κ(θu)(6)對(duì)土柱暖端采取同樣的方法處理。3.2b土壤蒸發(fā)通量溶質(zhì)差分方程邊界采用與水流差分方程邊界相似的處理,如圖4,在i為0～0.5處列平衡方程:(Es-J1+1/2)△t=△c△z12(7)式中表面蒸發(fā)通量Es=0。J1+1/2=-Dsh?c?z+qc(8)下邊界采取同樣處理方法。4土壤特征參數(shù)4.1抗凍土法未凍土的導(dǎo)水率和擴(kuò)散率可由實(shí)驗(yàn)測(cè)得,而凍土的導(dǎo)水率和擴(kuò)散率測(cè)定非常困難,目前國(guó)內(nèi)尚無(wú)測(cè)定先例。確定凍土導(dǎo)水率和擴(kuò)散率的方法有很多種,本文采用阻抗系數(shù)法,即用凍土區(qū)的導(dǎo)水率和擴(kuò)散率除以阻抗系數(shù)I來(lái)表示,I的大小取決于土壤中含冰量的多少。Taycor和Luthin提出I由下式確定:Ι=1010θΙ(9)未凍區(qū)的導(dǎo)水率和擴(kuò)散率采用實(shí)驗(yàn)值:Κ(θ)=1.5089×10-10e38.666θ(10)D(θ)=0.0017e17.596θ(11)凍結(jié)區(qū):Κ(θ)=1.5089×10-10e38.666θ/Ι(12)D(θ)=0.0017e17.596θ/Ι(13)4.2土壤中各組成物的熱容量的計(jì)算估算土壤熱導(dǎo)率是十分困難的,它不但取決于土壤中各組成物的比例,還和各組成物的形狀等因素有關(guān)。DeVries提出了一種半理論半經(jīng)驗(yàn)的計(jì)算方法,并且福田提出,土壤熱導(dǎo)率應(yīng)該包括4相,土顆粒、水、空氣和冰。Κh=1σΚ1+1-σΚ2(14)式中K1、K2分別為串、并聯(lián)土層的熱導(dǎo)率;σ為串、并聯(lián)土層的體積比,一般取0.5。上式中K1、K2的計(jì)算,關(guān)鍵在于確定土壤中各組成物的熱導(dǎo)率:{Κ1=1Xλs+θλw+θΙλΙ+AλaΚ2=Xλs+θΙλΙ+θλw+Aλaθ+θΙ＜θs(15){Κ1=1Xλs+θλw+θΙλΙΚ2=Xλs+θΙλΙ+θλwθ+θΙ≥θs(16)式中X為土顆粒所占的體積百分?jǐn)?shù),%;θ為未凍水所占的體積百分?jǐn)?shù),%;θI為冰所占的體積百分?jǐn)?shù),%;A為空氣所占的體積百分?jǐn)?shù),%;λs為土顆粒的熱傳導(dǎo)系數(shù);λw為水分的熱傳導(dǎo)系數(shù);λI為冰的熱傳導(dǎo)系數(shù);λa為空氣的熱傳導(dǎo)系數(shù)。土壤體積比熱容具有可加性,是土壤中各組成物熱容量之和,本文采用經(jīng)驗(yàn)公式。cv=cvmxm+cvomxom+cvwxw+cvΙxΙ(17)式中xm、xom、xw、xI分別表示單位體積土壤中礦物質(zhì)、有機(jī)質(zhì)、水、冰占有的體積。cvm、cvom、cvw、cvI分別為相應(yīng)組成物的體積比熱容J/(cm3·K),其值分別為1.926、2.512、4.186、1.73。4.3u3000動(dòng)力擴(kuò)散系數(shù)目前,對(duì)水動(dòng)力彌散系數(shù)的結(jié)構(gòu)形式的認(rèn)識(shí)尚不統(tǒng)一。從理論上講,水動(dòng)力彌散系數(shù)Dsh為分子擴(kuò)散系數(shù)Ds與機(jī)械彌散系數(shù)Dh之和。一般將溶質(zhì)在土壤中的分子擴(kuò)散系數(shù)僅表示為含水率的函數(shù),而與溶質(zhì)的濃度無(wú)關(guān)。本文采用實(shí)驗(yàn)值:Dsh(θ)=0.418θ0.5796(18)5凍結(jié)特征曲線的建立考慮到非粘性土壤水中溶質(zhì)對(duì)水熱運(yùn)動(dòng)參數(shù)的影響和溶質(zhì)梯度在水熱運(yùn)動(dòng)中所起的作用一般可以忽略,可認(rèn)為土壤中的水熱運(yùn)動(dòng)和其中溶質(zhì)濃度的大小無(wú)關(guān)。因而,土壤水熱運(yùn)動(dòng)可先獨(dú)立求解,在此基礎(chǔ)上再求解溶質(zhì)的運(yùn)移。水、熱、溶質(zhì)運(yùn)移基本方程取垂直向下為正方向。鉀鹽溶液的濃度對(duì)凍結(jié)過(guò)程中未凍水含量有較大影響,濃度越大,影響越大,不同濃度的試樣在凍結(jié)過(guò)程中具有不同的凍結(jié)特征方程和決定系數(shù)。本文采用5%和7%的兩種不同濃度的氯化鉀凍結(jié)特征曲線進(jìn)行計(jì)算。差分方程與其邊界條件一起組成三對(duì)角方程組,可通過(guò)MATLAB語(yǔ)言編程求解。方程組是非線性的,本文中采用迭代法求解。在計(jì)算過(guò)程中采用變時(shí)間步長(zhǎng),以提高迭代計(jì)算速度。6凍結(jié)后土體含水率的變化在上述的邊界條件下對(duì)土壤剖面的溫度、含水率和溶質(zhì)變化進(jìn)行了仿真分析,如圖5、6、7。從圖5可以看出,土柱各層溫度模擬值和實(shí)驗(yàn)值吻合較好。3h以后,可以從圖中清楚地看出溫度隨時(shí)間的線性變化規(guī)律非常明顯。在127h,隨土柱高度的升高,溫度呈線性降低。圖6給出的是凍結(jié)后含水率實(shí)驗(yàn)值與模擬值的比較。從圖中可以看出,除4～6cm模擬誤差稍大外,其它各層溫度模擬值與實(shí)驗(yàn)值吻合較好。在4cm高度以上出現(xiàn)了水分增加,高于初始含水率,增加幅度較均勻,在5～7cm和9～10cm處出現(xiàn)了峰值,在4cm高度以下水分減少,含水率低于初始含水率,且隨著土柱高度的降低,含水率不斷降低,但在0～2cm處土柱含水率數(shù)值趨于穩(wěn)定,這主要原因在于隨著凍結(jié)端不斷結(jié)冰,未凍水含量減小,導(dǎo)致凍結(jié)端水勢(shì)梯度減小,水分向凍結(jié)端運(yùn)移,而不斷結(jié)冰累積的結(jié)果。圖7給出的是凍結(jié)后各土層中氯離子含量實(shí)驗(yàn)值與模擬值,模擬效果較好,說(shuō)明土體對(duì)氯離子的吸附量很小。從圖中可以看出,隨著土柱高度的降低,氯離子的濃度不斷降低,且降低幅度較均勻,水分從土柱暖端向冷端遷移,溶質(zhì)隨水分的運(yùn)動(dòng)產(chǎn)生對(duì)流運(yùn)動(dòng),在土柱冷端集聚,氯離子濃度在7cm高度以上高于初始值,在7cm高度以下低于初始值。7未考慮到正反合物-水溶液-液(1)本文通過(guò)建立一維水、熱、溶質(zhì)耦合運(yùn)移數(shù)學(xué)模型,利用Crank-Nicolsen差分格式及線性化迭代法對(duì)凍結(jié)條件下土壤水、熱、溶質(zhì)耦合運(yùn)移進(jìn)行模擬,含水率誤差為10%、溫度誤差為15%、氯