滲流、溫度場、應(yīng)力場的耦合作用

滲流、溫度場、應(yīng)力場的耦合作用

1巖石體多孔介質(zhì)的耦合地下水和地應(yīng)力是地質(zhì)環(huán)境的主要因素。它們之間的聯(lián)系、相互作用和制約性是相互聯(lián)系、相互制約的，因此在地下水、應(yīng)力場和溫度場之間形成了多階段耦合效應(yīng)。該多場耦合問題廣泛存在于水利水電工程、環(huán)境工程、石油開發(fā)、地?zé)衢_發(fā)、地震預(yù)報(bào)與控制、地下工程以及核廢料處置等諸多領(lǐng)域,是人類活動和地質(zhì)環(huán)境相互作用研究的重要課題。在水利水電工程和巖土工程實(shí)際應(yīng)用中,考慮巖土類多孔介質(zhì)的多場耦合作用,可以更加準(zhǔn)確地評價(jià)高壩復(fù)雜基礎(chǔ)、高地應(yīng)力和高外水壓力下的深埋洞室、高邊坡巖體的穩(wěn)定性,并為采取有效、合理的滲控措施提供科學(xué)依據(jù)。巖土類多孔介質(zhì)的多場耦合問題實(shí)際上是由固體力學(xué)、流體力學(xué)、傳熱學(xué)、傳質(zhì)學(xué)、物理化學(xué)與反應(yīng)理論等基礎(chǔ)學(xué)科與眾多工程學(xué)科相互交叉而形成的新興邊緣學(xué)科。研究巖土體多場耦合問題的關(guān)鍵在于耦合數(shù)學(xué)模型的建立,目前國內(nèi)外研究者已經(jīng)提出了很多數(shù)學(xué)模型,經(jīng)歷了由初期的研究多孔介質(zhì)飽和狀態(tài)下的多場耦合情形,發(fā)展到現(xiàn)在主要研究多孔介質(zhì)在非飽和-飽和狀態(tài)下的多場耦合情形,國外如Nguyen等在Biot固結(jié)理論的基礎(chǔ)上提出了適用于飽和孔隙介質(zhì)的非等溫固結(jié)的數(shù)學(xué)模型;Bower等針對裂隙巖體提出了飽和雙重介質(zhì)的多場耦合數(shù)學(xué)模型;Rutqvist等提出了一種非飽和介質(zhì)的多場耦合數(shù)學(xué)模型,該模型中較為充分地考慮了多孔介質(zhì)的多相效應(yīng)以及各場間的相互作用關(guān)系。國內(nèi)方面,如賴遠(yuǎn)明、李寧等針對凍土問題提出了帶相變的多場耦合數(shù)學(xué)模型;劉亞晨等對核廢料貯存庫圍巖裂隙巖體介質(zhì)多場耦合問題進(jìn)行了探討,等等。這些數(shù)學(xué)模型假設(shè)有些是結(jié)合具體的問題提出,這樣得到的數(shù)學(xué)模型應(yīng)用范圍較窄;有些提出的假設(shè)雖然可以使模型大大簡化,便于數(shù)值計(jì)算,但是與多場的實(shí)際耦合過程并不一定完全符合。本文在總結(jié)國內(nèi)外學(xué)者對多孔介質(zhì)多場耦合的機(jī)理、數(shù)學(xué)模型建立及求解方法的研究成果基礎(chǔ)上,以巖土工程中的巖土類多孔介質(zhì)為研究對象,綜合運(yùn)用固體力學(xué)、水力學(xué)、熱學(xué)等基本理論,結(jié)合多孔介質(zhì)的熱本構(gòu)關(guān)系以及孔隙流體的熱運(yùn)動規(guī)律,建立了多孔介質(zhì)多相、多場全耦合數(shù)學(xué)模型方程組,在此基礎(chǔ)上采用加權(quán)余量法和有限差分法對其分別進(jìn)行了空間和時(shí)域上的離散,并討論了數(shù)值求解方法的穩(wěn)定性,最終得到了多孔介質(zhì)的多相、多場耦合的有限元求解模型。2密度vl將非飽和狀態(tài)下巖土工程中的巖土類多孔介質(zhì)視為連續(xù)性介質(zhì),其概念模型為:假設(shè)一個(gè)非飽和多孔介質(zhì)的特征單元體主要由三相構(gòu)成:固相、液相和氣相,其體積V可以用下式表示:其中Vs、Vl、Vg分別是固體介質(zhì)、液體介質(zhì)和氣體介質(zhì)的體積。假設(shè)ρs、ρl、ρg分別代表固體介質(zhì)、液體介質(zhì)和氣體介質(zhì)的密度,ms、ml、mg分別表示固體介質(zhì)、液體介質(zhì)和氣體介質(zhì)的質(zhì)量,且滿足m=ρV,ms=ρsVs,ml=ρlVl,mg=ρgVg。定義ns=Vs/V,nl=Vl/V,ng=Vg/V分別代表固相、液相、氣相的體積分?jǐn)?shù);定義表征單元體的孔隙體積為Vv,則有Vv=Vl+Vg,相應(yīng)的,定義n為孔隙率,有:定義單元體中液體的飽和度為S,且滿足S=Vl/Vv。3多孔熱耦合結(jié)構(gòu)的關(guān)系模型3.1根據(jù)本構(gòu)關(guān)系求解熱應(yīng)力應(yīng)變關(guān)系當(dāng)溫度發(fā)生變化時(shí),多孔介質(zhì)將由于受熱發(fā)生膨脹、受冷而發(fā)生收縮,從而產(chǎn)生應(yīng)變。如果多孔介質(zhì)熱變形受到邊界的約束,就會產(chǎn)生應(yīng)力。這種由于溫度變化而引起的應(yīng)力可以稱為“熱應(yīng)力”或“溫度應(yīng)力”?？紤]等溫條件下的巖土體本構(gòu)關(guān)系比較常用,如反映非飽和土中力學(xué)行為的Camclay模型、Alonso-Gens模型和Bolzon-Schrefler模型等;反映材料塑性屈服行為的Mohr-Coulomb模型、Drucker-Prager模型等,為了利用這些成熟的本構(gòu)關(guān)系,同時(shí)考慮溫度變化對多孔介質(zhì)的影響,可以給出以增量形式表述的熱本構(gòu)關(guān)系。另外,在非飽和多孔介質(zhì)中,對變形起主要作用的是有效應(yīng)力,于是可以將非飽和多孔介質(zhì)的熱應(yīng)力應(yīng)變關(guān)系寫為:式中:σ′為有效應(yīng)力張量;D為剛度張量,可以為彈性張量,可以為彈塑性張量;ε為總應(yīng)變張量;εT=βT(T-T0)δ可以為受溫度影響而產(chǎn)生的應(yīng)變張量;βT為固體骨架的熱膨脹系數(shù)。3.2相對滲透率系數(shù)式中:k為介質(zhì)的本征飽和滲透張量;krl為液體的相對滲透性系數(shù),非飽和區(qū)域krl∈(0,1),飽和區(qū)域krl=1.0;klT為液相的熱耦合張量。3.3相熱運(yùn)動的廣義darcy規(guī)則式中:krg為氣體的相對滲透性系數(shù);kgT為氣相的熱耦合張量;μg為氣相的動力粘滯系數(shù);vgr為氣體相對于固體骨架的速度向量。3.4熱通量矢量t式中:q為熱流密度矢量即單位時(shí)間內(nèi)通過單位面積的熱通量矢量;T表示溫度分布矢量;λ為導(dǎo)熱系數(shù),其影響因素包括介質(zhì)的種類、材料組分、濕度、壓力等。4廣義darcy定律推導(dǎo)多孔介質(zhì)非飽和狀態(tài)下多場耦合的控制方程將主要基于上述非飽和多孔介質(zhì)的概念模型,同時(shí)將遵循以下基本假設(shè):(1)假設(shè)多孔介質(zhì)為連續(xù)、均勻介質(zhì),多孔介質(zhì)中的液體和氣體只在連通孔隙中流動,氣體為理想氣體,且液體和氣體的運(yùn)動滿足式(4)和式(5)所示的廣義Darcy定律;(2)假設(shè)多孔介質(zhì)中的固相介質(zhì)的密度是應(yīng)力、吸力、溫度以及水壓力的函數(shù);流體介質(zhì)(包括液體和氣體)的密度是溫度和氣壓的函數(shù);(3)假設(shè)多孔介質(zhì)的變形滿足小變形假設(shè),而且其應(yīng)力狀態(tài)滿足Biot有效應(yīng)力原理;無論是飽和狀態(tài),還是非飽和狀態(tài),應(yīng)力應(yīng)變關(guān)系滿足式(3)所示的熱應(yīng)力應(yīng)變關(guān)系;應(yīng)力以拉應(yīng)力為正,壓應(yīng)力為負(fù);(4)假設(shè)多孔介質(zhì)滿足局部熱平衡,即同一局部區(qū)域內(nèi)固、液、氣三相的溫度相同;液相與氣相不相混溶,且不考慮液相與氣相的相互轉(zhuǎn)化,熱量在多孔介質(zhì)中的傳導(dǎo)過程滿足式(6)所示的Fourier定律。多孔介質(zhì)多場耦合的控制方程主要包括變形場控制方程、連續(xù)性方程以及能量守恒方程構(gòu)成。4.1應(yīng)力增量算子pg式中σ為總應(yīng)力張量。對于多孔介質(zhì)非飽和狀態(tài)下,總應(yīng)力張量的增量形式可以用Biot有效應(yīng)力原理來表示,即:式中σ′為有效應(yīng)力張量;δ為Kronecker張量算子;αl、αg為有效應(yīng)力增量系數(shù);pl,pg分別為孔隙水壓力和孔隙氣壓力。將介質(zhì)熱應(yīng)力應(yīng)變關(guān)系代入上式,并將總應(yīng)力表示的平衡方程替換成有效應(yīng)力表示的平衡方程,于是得到增量形式的變形場控制方程有:4.2物質(zhì)導(dǎo)數(shù)的及其應(yīng)用考慮到液體的密度是壓力和溫度的函數(shù)1/ρldρl=clpdpl+clTdT,其中clp為液體的壓縮性系數(shù)、clT為液體的熱膨脹系數(shù),二者滿足clp=?ρl/ρl?pl,clT=-?ρl/ρl?pl?？紤]到飽和度為溫度、孔隙率以及吸力等的函數(shù),而其中吸力和溫度是影響飽和度的最重要因素,因此有:dS=cSpd(pg-pl)+cSTdT,其中cSp=?S/?s為吸力對飽和度的影響系數(shù);cST=?S/?T為溫度對飽和度的影響系數(shù)。引入物質(zhì)導(dǎo)數(shù),式中vs為固相介質(zhì)的絕對速度向量。對于物性導(dǎo)數(shù)來說,在小變形條件下vs一般較小,可以認(rèn)為,因此物性導(dǎo)數(shù)可以認(rèn)為是D(×)/Dt≈?(×)/?t,于是。由式(10)可以得到液相介質(zhì)的連續(xù)性方程:式中us為固體骨架的位移向量;εv為固體骨架的體積應(yīng)變。同樣地,多孔介質(zhì)中的氣相介質(zhì)的連續(xù)性方程可以寫成:式中vgr為氣體平均相對速度向量?？紤]到氣相為理想狀態(tài)下的氣體,因此其狀態(tài)可以用以下方程來描述:式中:mg為氣體的分子量;R為常數(shù)。一般地,R/mg為常數(shù),于是可以得到:1/ρgdρg=1/pgdpg-1/TdT,且定義cgp=1/pg,cgT=1/T分別為氣體的壓縮系數(shù)和溫度膨脹系數(shù),結(jié)合式(10)—式(13),可以得到多孔介質(zhì)中理想氣體的連續(xù)性方程:同樣地,將固相介質(zhì)的連續(xù)性方程寫為:固相介質(zhì)的密度可以寫為:式中Ks為固體顆粒的壓縮模量。如果將多孔介質(zhì)視為彈塑性材料,近似地認(rèn)為,其中K1、K2為系數(shù),如果與多孔介質(zhì)的體積變形模量K結(jié)合起來,有:,式中λ為Lame常數(shù)。于是可以得到考慮材料彈塑性的固相介質(zhì)連續(xù)性方程:4.3多孔介質(zhì)的熱重分析式中:qs為流入單元體中固相的熱通量;Cs是固相介質(zhì)的比熱;ρs是固相介質(zhì)的密度;t是時(shí)間;σ′s為作用于單元體中固相的應(yīng)力;Qs為單元體內(nèi)部熱源,如水泥水化熱等,表示在單位時(shí)間內(nèi)單位體積中發(fā)出的熱量。非飽和介質(zhì)單元體中液相的能量守恒方程:非飽和多孔介質(zhì)單元體中氣相的能量守恒方程:可以看出,σ′s、σ′l、σ′g滿足下式:且有:根據(jù)廣義Fourier定律,多孔介質(zhì)的總熱通量可以寫成:式中:q、qs、ql、qg分別多孔介質(zhì)總熱通量、固相的熱通量、液相的熱通量和氣相的熱通量;λ為多孔介質(zhì)的導(dǎo)熱向量,是由固相、液相、氣相的導(dǎo)熱通量集合而成的,通常認(rèn)為是多孔介質(zhì)飽和度S的函數(shù)。聯(lián)合式(18)—式(21)可以得到非飽和多孔介質(zhì)的能量守恒方程:綜上所述,非飽和多孔介質(zhì)THM耦合的控制方程主要由式(9)、式(11)、式(12)、式(17)和式(22)構(gòu)成,其中包含位移向量u、孔隙水壓力pl、孔隙氣壓力pg、溫度T和孔隙率n等七個(gè)基本未知量。4.4面力向量及邊界液體流量(1)非穩(wěn)定非飽和狀態(tài)下的初值條件,即t=0時(shí)刻在Ω內(nèi)有:式中,t、分別為面力向量和已知面力向量;ql、分別為邊界液體流量和已知邊界液體流量;qg、分別為邊界氣體流量和已知邊界氣體流量;qT、分別為邊界熱通量和已知邊界熱通量;n為邊界面的外法向向量。5多孔介質(zhì)多孔耦合模型的有限離散5.1求解方程對非飽和多孔介質(zhì)水-熱-力全耦合數(shù)學(xué)模型采用加權(quán)余量法進(jìn)行有限元計(jì)算公式的推導(dǎo):將計(jì)算域Ω離散成有限個(gè)單元(如NE個(gè)),并將基本未知量寫成節(jié)點(diǎn)上的插值形式:,記應(yīng)變矩陣B=LNu,并經(jīng)過適當(dāng)?shù)恼?可以得到非飽和多孔介質(zhì)水-熱-力耦合問題的有限元求解方程組:其中:5.2耦合方程差分法非飽和多孔介質(zhì)的三相多場耦合方程組(36)是一個(gè)以時(shí)間t為變量的常微分線性代數(shù)方程組,有限元分析中求解該方程組通常有差分法:將求解的時(shí)間域劃分為若干個(gè)時(shí)間步:t0、t1、t2、…、tn、…,假設(shè)在步長Dt時(shí)間內(nèi)未知量列陣X是線性變化的,將未知量列陣X按泰勒級數(shù)在ti附近展開并取線性項(xiàng),并取X、F的差分格式,于是將耦合方程組(36)寫成如下增量形式:式中含上標(biāo)i、i+1量分別代表了時(shí)步i、i+1對應(yīng)的量;α為差分系數(shù),在0到1之間變化,α=0表示前差分格式(Euler差分格式),α=1表示后差分格式,α=0.5表示中間差分格式(Crank-Nicholson差分格式),α=0.667表示Galerkin差分格式;Dt為迭代的時(shí)間步長。5.3差分解穩(wěn)定收斂的條件考慮到微分方程的穩(wěn)定性與非齊次項(xiàng)無關(guān)上式成為:易知該方程的解析解為X=ae-ωt,其中α為任意常數(shù);ω=K/C且總是正值,相當(dāng)于上式的特征值。容易看出該解析解隨時(shí)間單調(diào)趨于零。假設(shè)α為正的常數(shù),則該解析解單調(diào)恒大于零。將式(39)寫成如下差分格式為:將ΔXi+1=Xi+1-Xi代入可以得到:可以看出,為了使差分解穩(wěn)定收斂,則λ必須滿足兩個(gè)基本要求:(1)|λ|<1,這是因?yàn)槿绻鹼λ|>1,則解是發(fā)散的;(2)λ>0,這是因?yàn)槿绻?lt;0,則解會在正負(fù)之間來回震蕩,不滿足解解析解單調(diào)恒正或恒負(fù)的要求。由于ω恒正且0ue08dαue08d1,則λ的表達(dá)式中1-ωΔti+1(1-α)/(1+αωΔti+1)總是小于1,因此為了滿足|λ|<1,必須有:可以看出,當(dāng)αue08e0.5時(shí),解是無條件穩(wěn)定的;當(dāng)0<α<0.5時(shí),解是有條件穩(wěn)定的,此時(shí)要求時(shí)間步長Δti+1滿足Δti+1<2/ω(1-2α)。另一方面,還需要保證解不發(fā)生震蕩即λ>0,此時(shí)時(shí)間步長Dt還需要滿足Δti+1<1/(ω(1-α))。可以看出當(dāng)α=1時(shí),Δti+1<∞。因此,為了保證耦合有限元方程解的穩(wěn)定收斂,將取α=1即后差分格式,此時(shí)解是無條件穩(wěn)定的,時(shí)間步長Δti+1的選擇從理論上講無任何限制條件此時(shí)耦合方程組(37)的實(shí)際迭代格式為:式中6基于廣義darcy定律的多相耦合非飽和巖土類多孔介質(zhì)多場耦合作用的數(shù)學(xué)模型,研究的主要是某一物理場方程中因變量或源匯項(xiàng)受其他物理場作用的變化的數(shù)學(xué)描述,也包括本構(gòu)規(guī)律的影響在控制方程中的反映,因此該類數(shù)學(xué)模型較為復(fù)雜,必定包含多個(gè)控制方程:(1)本文根據(jù)連續(xù)介質(zhì)方法給出了非飽和多孔介質(zhì)的概念模型,歸納總結(jié)了非飽和多孔介質(zhì)的熱本構(gòu)關(guān)系模型,包括變形場的熱本構(gòu)關(guān)系,流體運(yùn)動的廣義Darcy定律以及溫度場熱傳導(dǎo)的廣義Fourier;在給出基本假設(shè)的基礎(chǔ)上,基于質(zhì)量守恒定律、能量守恒定律詳細(xì)推導(dǎo)了非飽和多孔介質(zhì)多相多場耦合的主要控制方程,包括變形場控制方程;固相、液相、氣相的連續(xù)性方程;能量守恒方程等,建立了以位移、孔隙水壓力、孔隙氣壓力、溫度、孔隙率為未知量的非飽和多孔介質(zhì)多相多場全耦合研究的數(shù)學(xué)模型。(2)本文給出了非飽和多孔介質(zhì)多相多場耦合數(shù)學(xué)模型的相關(guān)定解條件,包括初值條件和水力、應(yīng)力、溫度邊值條件,采用加權(quán)余量法對非飽和多孔介質(zhì)多場耦合的數(shù)學(xué)模型進(jìn)行了空間和時(shí)間域的離散,獲得了有限元格式的多場耦合求解方程組,并詳細(xì)討論了該方程組解的穩(wěn)定性,為進(jìn)一步開展工程應(yīng)用分析奠定了基礎(chǔ)。在考慮溫度變化的條件下,溫度梯度同樣是引起滲流運(yùn)動的因素,為此可以將Dar