多孔介質多相流hm全耦合模型的構建及試驗驗證

1thm耦合模型高廢物地質處置、石油和液化天然氣的地下儲存、熱能源的開發、深層資源的開發、地下空間的開發和利用等建設項目通常與物質、能量和變量交換有關，并表現出復雜的熱耦合過程。thm。經過近30a的發展,特別是在DECOVALEX等國際及歐盟合作項目的推動下,巖土多場耦合研究取得了長足進展,開展了大量室內及現場試驗[3~7],建立了一系列數學模型及數值模擬方法[8~12],開發了數十套數值模擬程序,并在代碼驗證及成果比較等方面取得了豐富成果[13~16]。國內盡管尚未建立大型地下實驗室,但在巖石及緩沖材料的物理力學特性、地下貯庫的概念設計、THM耦合的數學模型及數值模擬等方面也開展了卓有成效的研究工作。理論上,巖土THM耦合數學模型的建立以連續介質力學中的動量、質量和能量三大守恒原理為基礎,其方法包括基于推廣到非等溫過程的Biot固結理論的宏觀唯像方法[23~25]和基于混合體理論的微觀均勻化方法[26~29]。盡管后者在理論上可能突破經驗性本構關系的束縛,但在當前的試驗技術條件下,2種方法都受到有效應力原理、土水特性和熱驅動效應等經驗性物理機制模型的制約。上述經驗性本構模型對于建立THM耦合數學模型是必須的,但由此建立的控制方程組往往出現非對稱、非協調、強非線性及高度病態等特征。特別是將氣體及蒸氣傳輸過程納入THM耦合體系之后,控制方程組的非線性及病態特征進一步惡化,使問題的求解難度增大。在DECOVALEX等國際合作項目的研究成果[13~16]中,氣體及蒸氣傳輸過程要么完全予以忽略,要么僅將蒸氣作為水的組分之一加以考慮,而尚未實現在THM耦合的框架下統一求解氣體及蒸氣傳輸。然而,對于核廢料貯庫等工程而言,氣體及蒸氣傳輸過程對屏障系統的安全性能及設計優化起著重要的作用。它通過水分的傳輸和運移影響屏障系統中的濕度分布及吸濕飽和過程,并通過吸濕膨脹效應、水力特性及熱傳導特性對應力–應變過程、水體傳輸過程和熱傳輸過程產生深刻影響。因此,工程屏障系統中,相對濕度的時空演化是THM耦合作用的綜合反映,同時也是安全性能評價的重要指標。在過去的研究中,相對濕度多通過借鑒Kelvin關系式將其表示為基質吸力和絕對溫度的負指數函數,其實質是在THM耦合體系中忽略氣體及蒸氣傳輸過程的情況下引入的補充關系,本身缺乏明確的物理意義。在物理學上,相對濕度定義為蒸氣壓力與飽和蒸氣壓力的百分比。采用基質吸力(嚴格地講,應該是飽和度或含水量)定義相對濕度顯然混淆了巖土中水分的存在形式?；谶B續介質力學原理和混合體理論,導出了多孔介質多相流THM全耦合數學模型。該模型從固、液、氣三相系統的動量、質量和能量守恒出發,考慮了應力–應變、水體流動、氣體傳輸、蒸氣傳輸、能量傳輸和孔隙率演化等6個過程的耦合作用,實現了對相變、溶解、熱驅動、濕度傳輸和吸濕膨脹等物理現象的模擬,從而力圖全面反映多相流THM耦合理論發展的豐富成果,確保控制方程組的封閉性和協調性,并使相對濕度的計算和評價建立在嚴格的物理學意義之上。其次,通過選取位移、孔隙水壓、氣壓、蒸氣壓、溫度和孔隙率為基本未知量,建立了有限元數值計算格式,自主研發了多孔介質三維八自由度多相流THM全耦合有限元程序THYME3D。該程序通過克服數值求解的困難和對程序架構的合理設計,能夠實現任意有物理意義的單一過程、2個過程直至上述6個過程的全耦合分析。采用法國原子能委員會(CEA)開展的膨潤土THM耦合Mock-up試驗驗證了多相流THM全耦合數值模型及其有限元程序THYME3D的正確性,從而深化了對多相流THM全耦合控制方程組特性以及耦合過程中涉及的大量本構關系及參數的作用效應的認識。2多相流thm全耦合數學模型2.1不含水層存在的不穩定性相流質控制模型將巖土材料視為包含固相(s)、液相(l)和氣相(g)三相系統以及固體(s)、水(w)、干空氣(a)和水蒸氣(v)四種組分的等效多孔連續介質,其在當前構形中的表征單元體積及質量分別為式中:V,M分別為三相混合體的總體積和總質量。定義Vα和Mα=ραVα分別為α相(即s,l和g)的體積和質量,ρα為α相的本征密度。α相占據的體積比例為nα=Vα/V,并定義α相的表觀密度為顯然有式中:ρ為三相混合體的密度,ρ=M/V。三相混合體的孔隙體積為Vv=Vl+Vg,故其孔隙率可表示為定義飽和度為:則根據式(3)~(5),可得如下關系式:為了建立多相流THM全耦合數學模型,引入如下假設:(1)液相流體和氣相流體在介質的骨架空隙中運動,且每一相流體均包含水和干空氣2種組分。其中,氣相被視為由干空氣和水蒸氣2種組分構成的理想氣體,二者完全混合并占據相同的空間;而液相則包含水和溶解氣2種組分。(2)力學響應滿足小變形條件,且假定無論巖土介質處于飽和還是非飽和狀態,修正的Terzaghi有效應力原理都適用。(3)流體運動滿足廣義Darcy定律,氣體擴散符合廣義Fick定律,混合氣體中水蒸氣和干空氣的壓力滿足Dalton分壓定律,且飽和水蒸氣壓力符合Kelvin濕度定律。(4)巖體固相密度是靜水壓力、基質吸力、溫度和平均應力的函數,而液相和氣相密度則是溫度和壓力的函數,且后者滿足理想氣體狀態方程。(5)熱傳導符合廣義Fourier定律,并假定固、液、氣三相處于局部熱平衡狀態,即在空間中的每一點,固、液、氣三相具有相同的溫度,而相間的熱量交換僅存在于相變過程中。2.2種基本理論對于三相混合體中的任一單元,其準靜態過程的動量守恒方程為式中:σ為總應力張量,規定以拉應力為正;g為重力加速度向量;?為梯度算子。將Terzaghi有效應力原理推廣到非飽和情形,并假定孔隙氣和孔隙水的效應可分開考慮,則巖土介質的總應力增量可表示為式中:σ′為有效應力張量;δ為KroneckerDelta張量;pl,pg分別為孔隙水壓力和孔隙氣壓力,規定以壓力為正;αl,αg為增量有效應力系數,一般要求滿足αl+αg=1?？紤]介質在非等溫過程中的熱膨脹效應和吸濕過程中的膨脹效應,則其彈性應力–應變關系可進一步表示為式中:D為介質的彈性張量;ε為應變張量;T為溫度;K為介質的排水體積模量;βT為介質的體積熱膨脹系數;εsw為吸濕膨脹應變,通?？杀硎緸榛|吸力s或飽和度S的函數,其中基質吸力s定義為引入幾何方程則由式(7)~(9),可得動量守恒方程為式中:u為位移向量,t為時間。需要特別指出的是,上述動量守恒方程的建立僅涉及線彈性或非線性彈性本構關系。當考慮塑性或黏性本構關系時,還應使本構關系嚴格滿足熱力學第二定律的約束。2.3混合氣體狀態方程引入參考于固相骨架的物質導數并考慮到式(2),則多相混合體中質量守恒方程的普遍形式可表示為式中:vs為固相骨架的運動速度,vαr≡nαvαr為α相的相對表觀流速,vαr=vα-vs為α相參考于固相骨架的相對本征流速;Qα為相變產生的質量交換項。(1)固相骨架質量守恒方程根據式(6),(14),固相骨架質量守恒方程(即孔隙率演化方程)可表示為式中:ρs為固相的本征密度?？紤]溫度、流體壓力、基質吸力及體積應力對介質骨架的作用,介質的固相密度可表示為式中:εv為體積應變;α1為Biot系數,α1≡1-K/Ks,Ks為固相顆粒的體積模量?？紤]到則由式(15),(16)可得孔隙率演化方程為需要指出的是,式(17)中的近似等式(即前文提到的小變形假設)并非必須。放棄該近似等式,并對式(12)采用Lagrangian描述,即可在完全相同的框架下求解大變形多相流THM全耦合問題。(2)水體質量守恒方程根據式(14),(6),水體質量守恒方程可表示為式中:vlr為液相水體的相對表觀流速;ρl為水的密度;jlg為液相與氣相之間的濕度傳輸速率(汽化為正,液化為負),可表示為式中:psv為飽和蒸氣壓力,pv為水蒸氣壓力,ω為液相遷移系數。根據廣義Darcy定律,并考慮溫度梯度對液相流場的熱驅動效應,則水體的相對表觀流速為式中:k為介質的本征滲透率張量,krl為水體的相對滲透率,klT為水體熱驅動張量,μl為水的動力黏滯系數,vlr為液相參考于固相骨架的相對本征流速。對于正壓流,其質量可表示為壓力和溫度的函數。因此有式中:clp,clT分別為水體的壓縮性系數和熱膨脹系數，。飽和度S通常是基質吸力s、溫度T、孔隙率n等變量的函數,但以基質吸力和溫度的影響最為顯著。因此有其中,將式(21)~(23)代入式(19),并考慮到式(17),有(3)氣體質量守恒方程如前所述,氣體是干空氣(a)和水蒸氣(v)二者的混合體;同時,干空氣既可存在于氣相混合體中,也可溶解于水中。根據式(14),(6),氣相混合氣體的質量守恒方程可表示為式中:vgr為混合氣體的相對表觀流速;ρg為混合氣體的密度。設干空氣和水蒸氣完全混合,則有式中:ρa,ρv分別為干空氣和水蒸氣的密度。根據廣義Darcy定律,并考慮溫度梯度對氣相流場的熱驅動效應,則混合氣體的相對表觀流速為式中:krg為氣體的相對滲透率,kgT為氣體熱驅動張量,μg為氣體的動力黏滯系數,vgr為氣相參考于固相骨架的相對本征流速。視混合氣體為理想氣體,則根據理想氣體狀態方程,有式中:R為普適氣體常數;Mg為混合氣體的分子質量,滿足ρg/Mg=ρa/Ma+ρv/Ml,Ma和M1分別為干空氣和水蒸氣的分子質量;Rg為混合氣體參數,Rg=R/Mg。由式(28)可得式中:cgp,cgT分別為混合氣體的壓縮性系數和熱膨脹系數,cgp=1/pg,cgT=1/T。需要說明的是,氣體狀態方程中涉及的溫度T均應理解為絕對溫度。另一方面,根據式(6),(14)和Henry定律,溶解于水中的干空氣的質量守恒方程可表示為式中:H為Henry溶解系數。干空氣的理想狀態方程為式中:pa為干空氣壓力;Ra為干空氣的氣體常數,Ra=R/Ma。根據Dalton分壓定律,有因此,有式中:cap,caT分別為干空氣的壓縮性系數和熱膨脹系數,cap=1/pa,caT=1/T。將式(23),(27),(29)代入式(25),將式(21),(23),(33)代入式(30),然后對式(25)和(30)求和,并考慮到式(17),可得考慮溶解效應的氣體質量守恒方程:(4)水蒸氣質量守恒方程根據式(6),(14),水蒸氣的質量守恒方程可表示為式中:vvr為水蒸氣的相對表觀流速,可分解為水蒸氣在混合氣體中的擴散流速和水蒸氣隨混合氣體一起運動的對流流速之和,即式中:vvd,vva分別為水蒸氣的擴散速度和對流速度,vvd≡n(1-S)(vvr-vgr),vva≡n(1-S)vgr;vvr為蒸氣參考于固相骨架的相對本征流速。其中,對流速度由式(27)確定。根據廣義Fick定律,水蒸氣的擴散速度與其濃度梯度成正比,即式中:Dv為水蒸氣在孔隙混合氣體中的分子擴散張量。同理,水蒸氣的理想狀態方程為式中:Rv為水蒸氣的氣體常數,Rv=R/Ml。于是,有式中:cvp,cvT分別為水蒸氣的壓縮性系數和熱膨脹系數,cvp=1/pv,cvT=1/T。將式(23),(23),(37),(39)代入式(35),并考慮到式(17),有2.4元熱傳導張量的計算引入質量守恒方程,則多相混合體中能量守恒方程的普遍形式可表示為式中:iα為α相熱傳導通量,σα為α相表觀應力張量,Qα為α相熱源;eα為α相比內能,可表示為式中:Cα為α相比熱,T0為參考溫度。設比熱的導數項可忽略不計,則三相混合體中固、液、氣三相的能量守恒方程可分別表示為式中:L為液相水的汽化潛熱;〈?〉為Macaulay括號函數:當x>0,〈x〉=x;否則〈x〉=0。在局部熱平衡假設條件下,對式(43)~(45)求和,即可得三相混合體的能量守恒方程。其中,根據廣義Fourier定律,三相混合體的熱傳導通量可表示為式中:λ為三相混合體的熱傳導張量,通常是飽和度S的函數。注意到表觀應力的定義式:以及如下恒等關系式:則三相混合體的能量守恒方程可最終表示為3計算thm全耦合值的格式3.1求解目標及邊界條件式(12),(18),(24),(34),(40),(48)構成了巖土介質多相流THM全耦合控制方程組,其基本未知量為位移u(t)、孔隙水壓力1p(t)、孔隙氣壓力pg(t)、水蒸氣壓力pv(t)、溫度T(t)和孔隙率n(t)。對于三維問題,基本未知量的個數為8個。多相流THM全耦合過程的問題提法為:設巖土介質的空間體積為?,邊界為Γ。問題的求解目標是在空間域?及時間域[0,τ]內尋找一組以{u(t),p1(t),pg(t),pv(t),T(t),n(t)}為基本未知量的解答,使之在?內滿足由式(12),(18),(24),(34),(40),(48)組成的控制方程組,同時在t=0時刻滿足如下初始條件:在Γ上滿足如下邊界條件:式(49)~(51)中:Γu,Γσ分別為已知位移和面力邊界;Γl,Γlq分別為已知液相流體壓力和流量邊界;Γg,Γgq分別為已知氣相混合流體壓力和流量邊界;Γv,Γvq分別為已知水蒸氣壓力和流量邊界;ΓT,ΓTq分別為已知溫度和熱通量邊界;分別為邊界上的已知位移和面力;分別為邊界上的已知水壓和水體流量;分別為邊界上的已知氣壓和氣體流量;分別為邊界上的已知蒸氣壓和蒸氣流量;分別為邊界上的已知溫度和熱通量;u0,pl0,pg0,pv0,0T和0n分別為初始位移、初始水壓力、初始氣壓力、初始蒸氣壓、初始溫度和初始孔隙率。需要指出的是,式(50),(51)考慮了各個過程之間的耦合作用對邊界條件的影響,且均可表示為時間的函數。此外,孔隙率不需要給定邊界條件,其演化完全受微分方程及初始條件控制。3.2thm全耦合過程的有限元迭代格式上述多相流THM耦合問題往往具有很強的非線性,因而只能用數值方法求解,這里采用有限單元法。(1)空間域離散引入形函數矩陣N,并對介質位移u(t)、孔隙水壓力p1(t)、氣壓力pg(t)、蒸氣壓力pv(t)、溫度T(t)和孔隙率n(t)等基本未知量采用相同的插值模式,則有式中:和分別為t時刻有限元節點的位移向量、孔隙水壓力向量、氣壓力向量、蒸氣壓力向量、溫度向量和孔隙率向量。顯然,對于u(t),N為矩陣;而對于1p(t),gp(t),pv(t),T(t)和n(t),N為行向量。采用相同的符號可精簡表達,但二者的區別很容易根據具體的表達式判斷。對控制偏微分方程組式(12),(18),(24),(34),(40),(48)應用Galerkin法,并通過引入邊界條件式(50),(51),可得多孔介質多相流THM全耦合過程的有限元支配方程組:其中,式(54a)~(54d)中,系數矩陣C,K以及右端項向量F各元素的計算式詳見附錄。(2)時間域差分現對式(53)進一步應用有限差分法,即可得隨時步遞進的迭代計算格式。令式中:i為時間步,?t為i時步到i+1時步的時間步長,θ為積分常數。當θ=0時,為向前差分格式;當θ=1時,為向后差分格式;當θ=0.5時,為中心差分格式;當θ=2/3時,為Galerkin差分格式。為保證數值計算的穩定性,應選取θ≥0.5。將式(55)代入式(53),得多相流THM全耦合過程的有限元迭代格式(3)求解策略多相流THM全耦合過程涉及8個基本變量,且由于耦合作用關系的強非線性和非對稱性,使有限元迭代方程組式(56)呈現強非線性、高度病態、非對稱、非協調和計算規模隨節點數迅速增長等突出特點。特別是氣體和蒸氣傳輸過程的引入,加劇了方程組的病態及非線性,使問題的求解難度增大。為確保多相流全耦合過程得到順利求解,以下求解策略至關重要:①真實物理力學過程的耦合作用是協調和耗散的,這種協調性和耗散性確保了多相流THM全耦合問題解的存在性。因此,在求解過程中,應盡量減少對數學模型的簡化,并應利用封閉方程組進行聯立統一求解。對數學模型的簡化盡管可以減少基本未知量,但往往是以喪失物理過程的協調性和封閉性為代價的;另一方面,過去常采用交錯法對耦合過程進行逐一迭代,這種方法由于在迭代過程中未能充分利用各過程的最新更新,因而也損失了部分協調性。協調性的喪失往往是導致數值求解收斂緩慢或求解失敗的重要原因。②在求解技術上,THM全耦合過程的求解涉及自動時間步進、耦合迭代和非對稱病態方程組求解三大技術。其中,耦合迭代技術涉及各種耦合變量的更新,是保證耦合關系協調性的關鍵;代數方程組求解技術則直接決定了求解耦合問題的能力和效率。對各項技術進行嚴密設計,并在代數方程組求解層面上施加嚴格的收斂控制,避免其出現非穩定性及這種非穩定性向耦合迭代和時步迭代層面的傳播,是確保耦合過程得到成功求解的關鍵。③不合理的經驗本構關系及其參數的引入是導致耦合過程喪失協調性的又一重要原因。對經驗本構關系的物理意義進行細致分析,正確把握其局限和適用條件,有助于為耦合求解指明正確的方向。④初始條件和邊界條件的協調性也是確保多相流THM全耦合過程得到順利求解的關鍵因素。初始條件的確定應以試驗為依據,結合特定本構關系取值,避免耦合初始參數出現矛盾現象。邊界條件的確定應考慮控制方程組邊界項的物理意義及各過程在邊界處的耦合作用效應。3.3tyne3-4耦合分析軟件的創新點基于前述多相流THM全耦合數學模型、有限元格式及求解策略,自主研發了多孔介質三維八自由度多相流THM全耦合有限元程序THYME3D。其主要特性如下:(1)THYME3D以位移、孔隙水壓、氣壓、蒸氣壓、溫度和孔隙率等8個變量為基本未知量,涉及應力–應變、水體流動、氣體傳輸、蒸氣傳輸、能量傳輸和孔隙率演化等6個過程。THYME3D通過程序架構的合理設計,可求解任意有物理意義的單個過程、2個過程耦合或多個過程耦合問題,這既方便了程序調試,又增強了程序功能。若以應力–應變過程的3個位移自由度為基礎,而將其他耦合過程的自由度數表示為X,此即為周創兵和龔玉鋒提出的3+X維耦合。THYME3D不但實現了3+X維耦合,還賦以了更強的數值模擬能力。(2)THYME3D通過對多相流THM全耦合過程協調性的深入把握和對自動時間步進、耦合迭代及非對稱病態方程組求解等3大核心技術的嚴密設計,確保了數值計算的穩定性,特別是氣體傳輸過程和蒸氣傳輸過程的數值穩定性。這是該程序區別于其他THM耦合分析軟件的重要特色之一。(3)THYME3D包括八節點六面體、六節點棱柱體和四節點棱錐體等單元類型,并實現了THM耦合過程中可能涉及的復雜邊界條件,能夠在任意時刻施加邊界、移除邊界或進行邊界類型轉換。(4)依賴于特定試驗成果的經驗性本構關系完全獨立于程序的其他主要模塊,可通過編譯成動態鏈接庫與主程序鏈接,從而保證了主程序代碼的穩定性和使用的便利性。(5)THYME3D的架構力圖全面反映多相流THM耦合理論發展的最新成果,能夠實現對相變、溶解、熱驅動、濕度傳輸和吸濕膨脹等物理現象的模擬。同時,其良好的程序架構及求解能力也為將來擴展到THMC全耦合分析奠定了堅實的基礎。4動力學過程的測試與分析為了確保耦合求解成果的正確性和可靠性,首先必須確保單一物理或力學過程得到正確求解。采用內壓作用下的厚壁圓筒問題、外水壓力作用下圓柱體的水體傳輸問題、外表面已知溫度下圓柱體的熱傳導問題等對THYME3D的應力–應變過程、流體運動過程、能量傳輸過程及求解器進行了廣泛的測試和考題,計算結果與理論解完全吻合,表明單一物理力學過程的計算成果準確、可靠。但限于篇幅,這方面的成果這里不予介紹。4.1材料和試件加工法國原子能委員會(CEA)針對核廢料處置開展的膨潤土THMMock-up試驗是目前較為先進、且得到廣泛引用的室內試驗之一,其目的是研究核廢料貯庫黏土屏障在高溫、大溫度梯度及吸濕過程中的THM耦合特性。依托歐盟合作項目“應用于貯庫安全性評價的THMC耦合過程”(THERESAProject),瑞典、芬蘭、英國、法國、德國、西班牙、荷蘭等的相關研究團隊應用該試驗成果開展了大量THM耦合分析代碼驗證及成果比較工作。試驗材料采用MX–80膨潤土,經33MPa高壓密實后加工成2塊不同含水量的試件Cell1和Cell2。下面僅考慮Cell1的試驗成果,Cell2完全類似。試件為圓柱形,直徑和高度均為203mm。試件依次置于PTFE套筒、不銹鋼襯及絕熱封套中,以確保試驗過程中體積不變,且熱能損失最小,如圖1所示。試件安裝后施以0.5MPa的軸壓。試驗包括2個階段。第一階段為2003年5月26日15:27~9月16日09:00。在此階段,試樣底部的溫度逐步從室溫(20℃)升高到150℃;第二階段為2003年9月16日14:23~2004年5月25日,此階段在維持試樣底部溫度為150℃的條件下,在試樣頂部施加1MPa的水壓力。在試驗過程中,通過試件內預埋的15個溫度傳感器、7個相對濕度傳感器、4個孔壓傳感器、8個徑向應力傳感器和1個軸壓測力計,對溫度、相對濕度、孔隙水壓力、徑向應力和軸向應力進行了全程監控。但受試驗控制條件和傳感器靈敏度的影響,僅有部分成果相對可靠。4.2維力學模型根據試驗條件,試件內的水、氣、熱傳輸過程為一維問題,但力學過程為軸對稱問題。因此,采用1/4柱體進行研究。數值模擬的有限元網格如圖2所示,其中八節點等參單元432個,節點629個。(2)初始本構關系根據試驗條件,試件的初始軸壓為0.5MPa,初始溫度為20℃,初始含水量為13.66%,初始飽和度為0.755,初始孔隙氣壓力為0.101325MPa,初始相對濕度為60%,初始孔隙率為0.3242,初始孔隙水壓及蒸氣壓則根據上述初始條件由隨后介紹的相關本構關系確定。對于力學過程,模型的上、下表面及側面均取法向支座約束。對于能量傳輸過程,模型的上、下表面均取已知溫度邊界,上表面的溫度值恒為20℃,下表面的溫度值在第一階段由20℃逐步升高到150℃,升溫過程見表1,在第二階段維持150℃不變;側面取絕熱邊界。對于滲流過程,模型的上表面在第一階段取隔水邊界,在第二階段取已知水壓邊界,水壓值為1MPa;下表面及側面均為隔水邊界。對于氣體及蒸氣傳輸過程,模型的所有邊界均為氣密邊界。然而,需要指出的是,試驗條件沒有保證氣密,因此氣體及蒸氣傳輸過程的邊界條件是近似的。(3)水體的本構關系前述多相流THM耦合數學模型涉及大量的本構關系和計算參數。下面結合Mock–up試驗成果,逐一加以介紹。但需要特別指出的是,多相流THM耦合過程涉及的本構關系是相當復雜的,有些經驗性本構關系至今尚未得到學術界的一致認可。限于篇幅,本文不擬對這些本構關系做深入的評述,僅結合試驗驗證對部分關鍵本構關系的作用效應加以討論。①MX–80膨潤土的本構關系及計算參數借鑒對數彈性模型,膨潤土的體積模量K(MPa)采用如下經驗公式:式中:為平均有效應力;κ為柔度系數,取為0.0025;K0為初始體積模量,取為16.7MPa;Kmax為體積模量上限,取為500MPa。根據MX–80膨潤土的吸濕膨脹特性,膨脹應變采用如下經驗公式計算:式中:s為基質吸力(MPa);?為試驗參數,取為0.004。膨潤土的土水特性曲線采用如下經驗關系式:式中:ρd為干密度,取為1.791g/cm3;a=10.96;b=41.89;p0=12.68MPa;m=0.211;η=7.97;n0=0.4;α=0.00647℃-1;T0=20℃。需要指出的是,式(59)是根據含水量–基質吸力關系曲線變換得到的,未必完全滿足s=0時,S=0。因此,使用時應加以修正。有效應力系數αl和αg可分別表示為式中:χ為Bishop系數。為了考慮非飽和介質的進氣延滯效應,并保證飽和條件下αl=α,采用下式計算:式中:sa為非飽和介質的進氣吸力,取為1kPa。膨潤土的熱傳導系數λ(W/(m·K))與飽和度S有關,可表示為式中:A1,A2,x0,dx為參數,分別為0.57,1.28,0.65,0.10。膨潤土本征滲透率k取1×10-21m2;熱膨脹系數βT取為1×10-5℃-1;比熱容Cs取為920J/(kg·℃)。②液相水的本構關系及計算參數借鑒液體–蒸氣界面之間的相變特征,液相遷移系數ω(m·s/kg)可表示為式中:C為調節系數,取為5.5×10-5。式(63)應理解為系數C通過單位體積上的質量遷移能力。水體的密度ρl(kg/m3)采用下式計算:式中:pl0為參考水壓,可取為0或一個大氣壓;ρlT為常壓下水體密度(kg/m3),隨溫度的變化,其表達式為在0.01℃~360℃溫度范圍內,式(65)擬合相關系數為R2=0.9997。水體的熱膨脹系數clT(℃-1)采用如下多項式表達:在0.01℃~260℃范圍內,式(66)擬合相關系數為R2=0.9975。水的動力黏滯系數μl(Pa·s)與絕對溫度T有關,即水體在介質中的相對滲透率krl可表示為飽和度S和孔隙率n的函數:式中:Ss為最大飽和度,取為1.0;Sr為殘余飽和度,取為0.01;0n為初始孔隙率,取為0.3242。考慮飽和度的影響,溫度梯度對水體的熱驅動系數klT(m2/(s·K))可表示為飽和度S的線性函數:式中:klT0為飽和條件下的熱驅動系數,取為0。③氣體和蒸氣的本構關系及計算參數根據Sutherland公式,氣體的動力黏滯系數μg(Pa·s)為式(70)各參數取值為:c=120K,0T=291.15K,μ0=1.827×10-5Pa·s。氣體在介質中的相對滲透率krg可表示為飽和度S的函數:溫度梯度對氣體的熱驅動系數kgT(m2/(s·K))可表示為飽和度S的線性函數:式中:kgT0為完全干燥狀態下的熱驅動系數,取為6.0×10-12m2/(s·K)。由于干空氣主要由氧和氮組成,根據二者在水中的溶解能力,Henry溶解系數H可表示為氣壓pg、蒸氣壓pv和絕對溫度T的函數:氣體的比熱容Cg可表示為pg和pv的函數:借鑒Kelvin關系式,孔隙氣中的飽和蒸氣壓力psv(MPa)可表達為基質吸力s(或含水量)和絕對溫度T的函數:式中:s0為初始基質吸力,根據初始條件由土水特性曲線計算得到;a,b為調節系數,取為a=2.3,b=0.35。式(75)的物理意義應理解為:在特定含水量和溫度條件下,孔隙氣中的水蒸氣在平衡態所具有的蒸氣壓力。由于飽和蒸氣壓可能受土體結構和孔隙特征等因素的影響,因此需要用調節系數a,b加以修正。水蒸氣在孔隙氣中的分子擴散系數Dv取為3.0×10-10m2/s。液–汽蒸發潛熱L取為2.46×106J/kg。4.3成果分析(1)計算結果的分析試件中從底至頂13個溫度傳感器T1~T13的實測溫度及計算溫度變化過程曲線如圖3所示(圖中的數據點為不同位置的實測值,實線為不同位置的計算值)。T0和T14傳感器分別位于試件的底部和頂部邊界,為邊界溫度,予以忽略。從圖中可見,計算溫度過程與實測溫度基本吻合。在第一階段,計算溫度主要受飽和度變化(通過熱傳導系數)和升溫過程的影響。由于在計算過程中,邊界的升溫過程采用每個升溫階段始末的溫度值對時間進行內插確定,而未考慮其穩定時間,因此與實測溫度存在一些偏差。在第二階段,計算溫度主要受飽和度的影響。在施加水壓之后的一段時間內,試件上部的飽和度迅速增大,由式(62)可知熱傳導系數有所增大,因此試件中、上部的計算溫度過程曲線略向下彎曲。隨后,隨著試件從上到下逐步吸濕飽和,熱傳導系數在試件內的差別逐步減小,計算溫度分布也逐漸趨于均勻。但在整個計算過程中,固、液、氣三相的應力功、相變能和對流傳輸對溫度分布的影響很小。此外,埋設在試件頂部197.50mm處T13傳感器的計算溫度與實測溫度差別較大。T13距頂部邊界僅5.5mm,但實測溫度卻明顯大于邊界溫度,且與位于其下6.25mm處T12傳感器的溫度值更為接近,這可能是由于T13的初始溫度率定誤差或埋設位置誤差造成的。(2)試件內部的孔隙率和飽和度分布試件在不同位置處的孔隙率和飽和度變化過程曲線如圖4,5所示。盡管孔隙率和飽和度不是試驗的直接指標,但其變化過程是THM耦合過程的綜合反映。因此,分析孔隙率和飽和度的變化過程有助于理解試件中的THM耦合機制。在第一階段,試件內的THM耦合過程受底部邊界的升溫過程及其內的熱能傳輸過程驅動。隨著試件從下至上逐步升溫,試件中、下部出現升溫膨脹、向上擠壓和孔隙收縮等現象,且水分蒸發逐步加劇,并通過溫度梯度對氣體及蒸氣的熱驅動效應使水分向試件中、上部遷移、凝結。從而表現出整個試件的孔隙率有所減小;同時,試件中、下部的飽和度明顯減小,而頂部的飽和度有所增大。在第二階段,試件頂部的水壓對孔隙率和飽和度的分布起關鍵作用。施加水壓之后,試件頂部迅速趨向飽和,且吸濕飽和過程逐步向試件的中、下部推進。但由于膨潤土的本征滲透率很小,吸濕飽和過程推進較為緩慢。經過250d吸濕過程之后,整個試件還遠未達到完全飽和狀態。與此同時,在試件頂部施加水壓之后,試件頂部附近的孔隙率先因水壓作用而明顯增大,而后則由于試件吸濕膨脹而導致孔隙率逐漸減小;但對于試件中、下部,在頂部水壓作用下,試件的孔隙率先明顯減小,而后隨著水體的入滲,在有效應力、水壓、氣壓及吸濕膨脹的共同制約下,孔隙率有緩慢增大趨勢。顯然,在定容試驗條件下,第二階段試件中、上部孔隙率的先劇增、后漸降以及試件中、下部孔隙率的先劇降、后緩增是固相骨架質量守恒的必然結果,這從一個側面反映了孔隙率演化方程對于維護THM耦合系統協調性所具有的重要作用。(3)濕度傳輸特性如前所述,孔隙介質中的相對濕度Hr常表示為基質吸力s和絕對溫度T的負指數函數,即然而,式(76)混淆了孔隙介質中水分的存在形式,且與濕度傳感器的物理工作機制不符。在物理學上,相對濕度定義為蒸氣壓pv與飽和蒸氣壓psv的百分比,即根據式(77),試件中7個相對濕度傳感器HR1~HR7的實測和計算相對濕度變化過程曲線如圖6所示,相應的計算蒸氣壓變化過程曲線如圖7所示。由圖6可知,計算給出的相對濕度過程曲線總體上較好地反映了試件中的濕度傳輸過程,表明給出的多相流THM耦合模型及相對濕度定義式(77)具備描述多孔介質中濕度傳輸特性的能力。相對濕度及蒸氣壓變化過程的物理機制如下:在試驗的初始階段,試件內的溫度梯度較小,其對氣體及蒸氣的熱驅動效應相對較弱,試件中的THM耦合作用以水分的緩慢蒸發為主,從而表現為相對濕度及蒸氣壓成整體增大趨勢。在第一階段隨后的試驗過程中,試件中的溫度梯度逐步增大。相應地,溫度梯度對氣體及蒸氣的熱驅動效應也逐步增強,促使試件中、下部高溫區的氣體及蒸氣向試件中、上部的低溫區遷移,從而表現出試件中、下部的相對濕度及蒸氣壓減小,而中、上部的相對濕度及蒸氣壓持續增大。在第二階段,頂部邊界的水壓及試件內的水體傳輸導致試件內的飽和度不斷增大,氣體及蒸氣占據的空間減小,從而導致蒸氣壓及相對濕度增大。在此階段,相變和熱驅動將影響曲線的斜率和形態。當蒸氣壓達到飽和蒸氣壓時,液–汽相變由蒸發轉變為凝結,從而使蒸氣壓逐步穩定到飽和蒸氣壓上,相對濕度也相應地穩定到100%。需要說明的是,孔隙介質中濕度傳感器的監測數據可能同時受孔隙氣中的蒸氣濃度和介質中其他水分的影響。因此,在這個意義上,式(77)也是近似的,但其物理意義較式(76)明確,特別是對于含水量較小、基質吸力較大的干土或極干土尤為如此。對于Mock-up試驗,若以式(76)定義相對濕度,則根據圖5給出的飽和度變化過程曲線,不可能得到類似于圖6的相對濕度變化過程曲線,也不可能對控制相對濕度變化過程的物理機制作出合理的解釋。(4)軸壓含量在初始和測定軸壓時軸壓在形態、溫度和模量等基試件的實測及計算軸向壓力變化過程如圖8所示。但需要指出的是,根據試驗報告,試驗的應力成果可靠性可能較差,僅有參考意義。在第一階段,試件所有邊界處于隔水狀態,盡管升溫過程將影響試件內水分的分布,但試件內總的含水量不變。因此,試件的體積熱膨脹是導致軸壓變化的主導因素,而水分遷移導致的體積膨脹及收縮則處于次要地位。在此階段,計算軸壓逐步從初始的0.5MPa增大到1.2MPa;但實測軸壓卻迅速增長到5.9MPa,超過了理論分析可能的軸壓范圍。在第二階段,試件內的溫度分布已基本趨于穩定,而水體傳輸過程導致的吸濕膨脹成為軸壓增大的主要驅動力。頂部邊界施加水壓之后,在水壓、氣壓及飽和度變化的共同作用下,軸壓迅速增大,并通過吸濕膨脹和體積模量調整,使軸壓增長過程與實測值基本吻合。需要指出的是,施加水壓前后實測軸壓的跌落與對試驗設備的操作有關。4.4土水特性分析如前所述,多相流THM耦合過程涉及大量的本構關系及計算參數。為了模擬試驗過