凍脹土水熱耦合模型的構(gòu)建

凍脹土水熱耦合模型的構(gòu)建

0凍結(jié)土水熱耦合模型在這項(xiàng)工作中，土壤充血破壞頻繁發(fā)生。影響凍脹的主要因素是溫度和水分,溫度負(fù)梯度引起水分遷移,水分凍結(jié)引起體積膨脹,而凍脹力是造成凍脹破壞的主要原因。Harlan.R.L.1973年提出土凍結(jié)過程中熱質(zhì)遷移的數(shù)學(xué)模型。此水熱耦合模型,描述了部分凍結(jié)土的水熱遷移問題。文對凍土水熱耦合問題進(jìn)行了總結(jié),并對凍土應(yīng)力場進(jìn)行了計(jì)算。文在水熱耦合模型的基礎(chǔ)上,考慮了凍土應(yīng)力,提出水熱力相互耦合的數(shù)學(xué)模型。本文根據(jù)此模型,進(jìn)行了水熱力的耦合迭代計(jì)算,并得出輸冷管道附近土體凍結(jié)過程的溫度場、水分場和應(yīng)力場及其凍脹位移。1土體凍脹位移計(jì)算本模型在Harlan.R.L.的水、熱耦合模型的基礎(chǔ)上考慮了冰壓力及水壓力,并將土體作為空間彈性體,進(jìn)行與溫度場、水分場相耦合的應(yīng)力場分析,最終可獲得溫度場、水分場、應(yīng)力場及土體凍脹位移。1.1經(jīng)驗(yàn)系數(shù)的確定ˉC?Τ?t=divCˉˉˉ?T?t=div(λgradΤ)+Lρw?θw?tT)+Lρw?θw?t,其中ˉC=C+Lρw?θu?Τ(1)Cˉˉˉ=C+Lρw?θu?T(1)溫度邊界條件:T|s=Tb或?Τ?n|s=Τb?T?n∣∣s=Tb或[?Τ?n+Τ]|s=Τb(1b)[?T?n+T]∣∣s=Tb(1b)溫度初始條件:T|t=0=T0(1i)方程中的變量分別為絕對溫度T,含水量θw,未凍水含量θu。其中未凍水含量是溫度的函數(shù),可表示為θu=A|T|-B,A、B為與土質(zhì)有關(guān)的經(jīng)驗(yàn)系數(shù)。參量分別為土質(zhì)熱容量C,土質(zhì)導(dǎo)熱系數(shù)λ,水相變潛熱L,水密度ρw。1.2w為水壓項(xiàng)?θw?t=div?θw?t=div(KgradPw)(2)水分初始條件:θw|t=0=θ0(2i)式中Pw為水壓項(xiàng),可由Clapeyron方程得:Ρw=(1-χ)ρwρiLln(Τ/Τ0)+ρwˉσχρw+(1-χ)ρi(3)Pw=(1?χ)ρwρiLln(T/T0)+ρwσˉχρw+(1?χ)ρi(3)參量分別為土質(zhì)導(dǎo)水系數(shù)K,冰密度ρi,權(quán)值χ=(θu/n)1.5,水的冰點(diǎn)T0,靜水壓力ˉσσˉ,土質(zhì)空隙率n,α為系數(shù)(1/9.81×10-3m/Pa)。1.3凍土溫度場、水分場、應(yīng)力場相互耦合的數(shù)學(xué)模型[L]{σ}+[G]=0,{σ}=[D]({ε}-{ε0}),{ε}=[L]{u}(4)應(yīng)力邊界條件:[M]{σ}=[F],位移邊界條件:{u}|s={ˉu}(4b){u}|s={uˉ}(4b)初應(yīng)變{ε0}=13{εV,εV,εV?0?0?0}Τ方程中{u}為位移量,{ˉu}為邊界位移,[G]為體積力,[F]為邊界應(yīng)力,[D]由楊氏模量E和泊松比μ表示,[M]由邊界外法線的方向余弦表示。εV為土體體積變形量,可由水相變成冰引起的體積增量和水分遷移引起的體積增量來表示,即ΔεV=0.09(Δθw-Δθu)+Δθw。由此,方程(1)、(1b)、(1i)、(2)、(2i)、(3)、(4)、(4b)構(gòu)成了凍土溫度場、水分場、應(yīng)力場相互耦合的數(shù)學(xué)模型。方程中土質(zhì)的熱容量C、導(dǎo)熱系數(shù)λ、導(dǎo)水系數(shù)K都依賴于溫度,而在實(shí)際工程中給不出這些參量與溫度的函數(shù)關(guān)系。因此,在計(jì)算中將土體以溫度為界分成3個區(qū)域,即凍結(jié)區(qū)、劇烈相變區(qū)、未凍區(qū),在這3個區(qū)域分別給出上述參量,即如下所示(其中T0為水的冰點(diǎn),Ts為分凝溫度):C={C+(C++C-)/2C-λ={λ+(λ++λ-)/2λ-Κ={Κ+(Κ++Κ-)/2Κ-Τ＞Τ0Τ0≥Τ≥ΤsΤ＜Τs未凍區(qū)劇烈相變區(qū)凍結(jié)區(qū)2水熱耦合模型的數(shù)值模擬2.1溫度和水分生成論對方程(1)、(2)用迦遼金法進(jìn)行數(shù)值離散,其形式如下:∫V[ˉCΝi?Τ?t+λ(?ΤΝi)?(?Τ)-LρwΝi?θw?t]dV-∫sλΝi(?Τ)dS=0∫V[Νi?θw?t+αΚ(?ΤΝi)?(?Ρw)]dV-∫sαΚΝi(?Ρw)dS=0其中{Nj}為形函數(shù),令T=Nj(T)j,θw=Nj(θw)j,Pw=Nj(Pw)j,則有[A]{˙Τ}+[B]{Τ}+[C]{˙θw}=0[D]{˙θw}+[E]{Ρw}=0}(5)其中[A]=(aij)=∫VˉCΝiΝjdV,[B]=(bij)=∫Vλ(?ΤΝi)?(?Νj)dV-∫sλΝi?(?Νj)dS[C]=(cij)=-∫VLρwΝiΝjdV,[D]=(dij)=∫VΝiΝjdV[E]=(eij)=∫VαΚ(?ΤΝi)?(?Νj)dV-∫sαΚΝi?(?Νj)dS{˙Τ}={˙Τj}?{Τ}={Τj}?{θw}={(θw)j}?{Ρw}={(Ρw)j}(i=1,??n；j=1,??n)式中,字母上面的符號‘·’表示對時間的導(dǎo)數(shù)。從方程(5)中可看出溫度{T}、水分{θw}、與應(yīng)力相關(guān)的水壓{Pw}是相互聯(lián)系的。再將{T}和{θw}的變化速率{˙Τ}和{˙θw}用后差分法進(jìn)行離散化,則方程(5)變?yōu)?[Η]{{Τ}t+1{θw}t+1}=[W]{{Τ}t{θw}t}+{R}(6)其中[Η]=[[A]t+1+Δt[B]t+1[C]t+10?[D]t+1]?[W]=[[A]t+1[C]t+10[D]t+1]?{R}=[0-Δt[E]t+1{Ρw}t+1]2.2項(xiàng)耦合的實(shí)現(xiàn)方程(6)、(3)及(4)是溫度場、水分場、應(yīng)力場的3項(xiàng)耦合方程組。但從方程(6)中可看出,其3項(xiàng)耦合可分為兩種耦合:一種是溫度{T}與水分{θw}的耦合;另一種是溫度{T}、水分{θw}與水壓{Pw}的耦合。因此,在這里用兩層迭代來實(shí)現(xiàn)這兩種耦合,即內(nèi)層迭代實(shí)現(xiàn)溫度{T}與水分{θw}的耦合,外層迭代實(shí)現(xiàn)溫度{T}、水分{θw}與水壓{Pw}的耦合。下面分別寫出這兩層迭代的計(jì)算步驟。2.2.1凍土的物理參量和t.t+1的確定在這一層迭代中假設(shè)水壓項(xiàng){Pw}為不變的已知量,則從方程(6)可寫出它的迭代步驟。步驟1:設(shè)允許的收斂誤差為ε、時間步長為Δt,已知t時刻的溫度{T}t、水分{θw}t、水壓{Pw}t,則可確定凍土在此刻的一些物理參數(shù),如熱容量ˉC、導(dǎo)熱系數(shù)λ、導(dǎo)水系數(shù)K等,然后計(jì)算t+1時刻的溫度{T}t+1、水分{θw}t+1。步驟2:再次確定凍土在已求出的溫度{T}t+1、水分{θw}t+1狀態(tài)下的物理參量ˉC、λ、K,并重新計(jì)算溫度{T}′t+1和水分{θw}′t+1。步驟3:判斷是否滿足|{T}′t+1-{T}t+1|<ε、|{θw}′t+1-{θw}t+1|<ε。若滿足,則令{T}t+1={T}′t+1、{θw}t+1={θw}′t+1,然后退出迭代循環(huán);否則,令{T}t+1={T}′t+1、{θw}t+1={θw}′t+1,返回步驟2。2.2.2應(yīng)力和位移層狀的見表1在這一層迭代中,根據(jù)已知t+1時刻的溫度{T}t+1、水分{θw}t+1,計(jì)算出水壓{Pw}t+1,并把它代回到內(nèi)層迭代中,從而重新求出溫度{T}t+1和水分{θw}t+1。具體計(jì)算步驟如下。步驟1:設(shè)允許的收斂誤差為ε,已知t時刻、t+1時刻的溫度、水分,即{T}t、{θw}t、{T}t+1、{θw}t+1和t時刻的應(yīng)力{σ}t、位移{u}t。步驟2:計(jì)算t時刻、t+1時刻的未凍水含量{θu}t、{θu}t+1和水分{θw}t、{θw}t+1,并求得t+1時刻的初應(yīng)變的增量{ΔεV}t+1。步驟3:在已知外荷載和體積初應(yīng)變增量{ΔεV}t+1的情況下,計(jì)算t+1時刻的應(yīng)力增量{Δσ}t+1和位移增量{Δu}t+1,由此可求出t+1時刻的應(yīng)力{σ}t+1={Δσ}t+1+{σ}t和位移{u}t+1={Δu}t+1+{u}t以及靜水壓力(ˉσ)t+1。并確定t+1時刻的水壓{Pw}t+1。步驟4:將已知的水壓{Pw}t+1代回到內(nèi)層迭代,從而求出新的t+1時刻的溫度{T}′t+1和水分{θw}′t+1。步驟5:判斷是否滿足|{T}′t+1-{T}t+1|<ε、|{θw}′t+1-{θw}t+1|<ε。若滿足,則令{T}t+1={T}′t+1、{θw}t+1={θw}′t+1,然后退出迭代循環(huán);否則,令{T}t+1={T}′t+1、{θw}t+1={θw}′t+1,返回步驟2。3為解決關(guān)閉管道附近土壤的凍結(jié)過程，數(shù)值模擬3.1輸入數(shù)據(jù)管道區(qū)域的截面圖如圖1所示。(1)初始條件溫度T|t=0=4.0℃含水量θ|t=0=40%(2)應(yīng)力值溫度T|S1=-0.75℃T|S2=-2.0℃T|S5=4.0℃、?Τ?n|S3+S4+S6=0位移:在S3、S4上ux=0,在S5上ux=0、uy=0應(yīng)力:在S3、S4上τxy=0,在S6上σx=0、τxy=0(3)k-、t015.λ-=7.236kJ/h·m·℃;K+=2.8979×10-9m/h,K-=5.95×10-14m/h;L=334.7kJ/kg;T0=0℃,Ts=-0.1℃;E=11.2MPa、μ=0.3;A=0.068988;B=0.4735;t∈[0,400h]3.2計(jì)算(1)圖2顯示了管道附近土壤區(qū)域的有限網(wǎng)格(2)溫度場圖3分別給出了管道附近土體區(qū)域凍結(jié)40h、200h、400h時候的凍結(jié)鋒面位置。圖4是凍結(jié)400h的等溫線分布圖。(3)b凍結(jié)鋒面周邊的含水量圖5是管道底部(x=0(m)土體的含水量隨深度的分布圖。因?yàn)楣艿赖撞客馏w的溫度梯度比其他部分大,所以發(fā)生水分遷移最多,凍結(jié)鋒面附近的含水量達(dá)到最大值。圖6是管道近處的、不同深度的含水量隨時間的變化圖。(4)土體凍脹的影響圖7的虛線是凍結(jié)400h后的凍脹位移。冷管的低溫引起管道附近土體凍脹,這使管道往上升,而且管道的上部位移比底部位移略大。在地表面上,管道近處的位移比其他地方大得多,其最大位移約為0.88×10-2m。(5)管道的深度分布圖8是管道近處土體的σx,σy,τxy隨深度的分布圖。圖9和圖10是在管道及其上部和底部位置土體的σx,σy隨深度的分布圖。由于應(yīng)力集中,σx在管道的上部和底部大,而在管道中部小,這是因?yàn)楣艿雷冃?在中部釋放了一部分應(yīng)力的原因。σy則相反地在管道上部和底部小,而在管道中部大,這是因?yàn)楣艿郎喜亢拖虏康拇怪眱雒浳灰撇煌鸬?。另外在各?yīng)力圖中虛線‘—·—·—·’處應(yīng)力發(fā)生變化,對照管道底部的含水量分布圖5,我們可看出,此處正好是含水量最大處。3.3凍土應(yīng)力場驗(yàn)證上面計(jì)算出的輸冷管道附近土體溫度場、水分場與文獻(xiàn)的實(shí)測情況相符合。而對于應(yīng)力場,因目前實(shí)測技術(shù)受限制,無法得到凍土應(yīng)力場的精確值,也就無法精確驗(yàn)證凍土應(yīng)力場。但可以定性地知道,即使在無約束條件下,由于凍土的不均勻凍脹,凍土內(nèi)部產(chǎn)生較大的應(yīng)力,且發(fā)生應(yīng)力集中。由凍土水熱力耦合方程和計(jì)算模擬結(jié)果可知,對于凍結(jié)土體溫度場、水分場、應(yīng)力場的三場耦合關(guān)系上,可以得到如下的結(jié)論:水分→應(yīng)力:強(qiáng)耦合,水分遷移多,應(yīng)力變化大;溫