朱媛媛,楊 驍,吳海濤

(1.上海師范大學信息與機電工程學院,上海 200234;2.上海大學力學與工程科學學院,上海 200444)

飽和多孔介質(zhì)熱-流-固耦合系統(tǒng)由于具有特殊的結(jié)構(gòu)和材料性質(zhì),不僅在傳統(tǒng)的應用領域(如土力學、水文學等)得到廣泛應用,而且也逐漸成為許多新興學科和高新技術(shù)(如核廢料污染物處置、地熱資源、熱能貯存、人體關(guān)節(jié)軟骨組織分析等)發(fā)展的關(guān)鍵.因此,研究飽和多孔介質(zhì)熱力學性能的理論、數(shù)值方法及其應用具有重要的意義.

基于Biot[1]的研究工作建立的有關(guān)飽和多孔介質(zhì)的熱彈性理論,目前已有許多研究成果.在連續(xù)介質(zhì)混和物公理體系和體積分數(shù)概念的基礎上,de Boer[2]建立了較完整的多孔介質(zhì)熱彈性理論,該理論可以直接將一些微觀性質(zhì)用來描述宏觀性質(zhì),而且容易反映非線性效應.利用混合物理論,de Boer等[3]為一維流體飽和不可壓縮多孔介質(zhì)的動力學固結(jié)問題提供了解析解.Heider等[4]研究了波在飽和多孔半空間中的傳播.Hu等[5]利用微分求積法對黏彈性流體飽和多孔介質(zhì)的動力學特性進行了研究.劉林超等[6]討論了飽和土中群樁的豎向動力相互作用問題.在de Boer等[7]提出的多孔介質(zhì)熱力學本構(gòu)關(guān)系的基礎上,He等[8]給出了一種在熱局部非平衡條件下的流體多孔介質(zhì)的數(shù)學模型,Yang[9]建立了該模型的相應的Gurtin型變分原理,朱媛媛等[10]對空間軸對稱流體飽和多孔熱彈性柱體動力學特性進行了分析.但是,目前還較少看到有關(guān)飽和多孔介質(zhì)熱-流-固耦合系統(tǒng)非線性理論和應用以及數(shù)值模擬的報道.

本工作基于多孔介質(zhì)理論(porous media theory,PMT)[2,7],研究了在有限變形和熱局部平衡條件下,不可壓流體飽和多孔熱彈性半平面受到表面溫度載荷作用下的動力學特性,并在建立問題非線性數(shù)學模型的基礎上,通過采用微分求積(differential quadrature,DQ)法在空間域內(nèi)離散問題的非線性數(shù)學模型,二階后向差分格式來處理時間導數(shù),Newton-Raphson法求解非線性代數(shù)方程組,從而可得到問題的數(shù)值結(jié)果,分析半平面的動力學特性.為了驗證本方法的正確性,在無熱效應和忽略幾何非線性影響時,計算了不可壓流體飽和多孔彈性介質(zhì)一維動力學問題.與解析解[2]對比發(fā)現(xiàn),本方法是有效可靠的,且具有計算量小、精度高等優(yōu)點.最后,研究了在受到表面溫度載荷作用下,流體飽和多孔熱彈性半平面的動力學特性,考察了達西滲透系數(shù)、熱傳導系數(shù)、熱交換系數(shù)等材料參數(shù)以及幾何非線性的影響.

1 數(shù)學模型

圖1為平面直角坐標系中流體飽和不可壓多孔熱彈性半平面動力學模型.設半平面處于熱局部平衡狀態(tài)(固相介質(zhì)和液相介質(zhì)的變溫相同).不失一般性,假設作用于半平面表面的載荷或者溫度具有某種對稱性,取Oz軸為對稱軸,其正方向為由表面指向半平面內(nèi)部(見圖1(a)),可認為沿Oy軸方向的固相位移uy≈0和流相的相對流速wy≈0,并且各物理量均與坐標y無關(guān).假設平面介質(zhì)為各向同性線性熱彈性多孔固相骨架和理想流體,根據(jù)多孔介質(zhì)理論,在不可壓和小流速的假設下,滿足如下動力學方程.

圖1 流體飽和多孔熱彈性半平面動力學問題模型Fig.1 Model of fluid-saturated porous thermo-elastic half-plane problem

質(zhì)量守恒方程為

動量和動量矩守恒方程為

質(zhì)量守恒方程為

式中:ux,uz為固相位移;為固相有效應力;wx,wz為相對流速,p為有效孔隙壓,θ為變溫;物相φα(α=S,F,分別表示固相及液相)體積分數(shù)nα(nS+nF=1);宏觀質(zhì)量密度ρα(和微觀質(zhì)量密度ραR的關(guān)系為ρα=nαραR);流固耦合系數(shù)Sv=(nF)2γFR/κF,其中γFR為液相有效比重,κF為達西滲透系數(shù);熱交換因子β;熱傳導系數(shù)k;初始絕對溫度Θ0;固相Lame系數(shù)μS,λS,體積模量KSλS+2μS/3;熱膨脹系數(shù)αS;比熱系數(shù)cα和熱源強度rα.記ρc=ρScS+ρFcF,ρr=ρSrS+ρFrF.

對于各向同性線性熱彈性固相,本構(gòu)關(guān)系為

在有限變形的情況下,流體飽和多孔熱彈性介質(zhì)固相材料的幾何關(guān)系可以表示為

假設流體飽和多孔彈性半平面的表面處于理想排水狀態(tài),同時承受變溫φ(x,t)和垂直載荷q(x,t)的作用;底部為不排水剛性絕熱底面;側(cè)面不排水絕熱且在x-方向位移被約束(見圖1(b)).因為φ(x,t)和q(x,t)關(guān)于Oz軸對稱,可以在區(qū)域0≤x≤L0,0≤z≤H0內(nèi)考慮問題(見圖1(b)),這時的邊界條件為

假設t≤0時半平面處于靜止和熱力平衡狀態(tài),初始條件可以寫成

因此,式(1)～(6),邊界條件(7)和初始條件(8)構(gòu)成了在熱局部平衡條件下流體飽和多孔熱彈性半平面動力學問題的一個數(shù)學模型,其中基本未知量有:固相位移ux、uz;固相有效應力相對流速wx、wz;有效孔隙壓p和變溫θ.

可以看出,上述動力學問題的解析或半解析解很難獲得.在本工作中,上述的非線性數(shù)學模型在空間域中使用DQM離散,時間導數(shù)采用二階后向差分格式處理,離散后的非線性代數(shù)方程組采用Newton-Raphson法迭代求解,可得到問題的數(shù)值模擬結(jié)果.

2 DQM和控制方程的離散化

Bellman等[11-12]提出了DQM,該方法具有公式簡單、精度高、計算量少等優(yōu)點,已有很多成功的應用[13-15].DQM的基本思想是將函數(shù)的偏導數(shù)近似為離散點處相應函數(shù)值的加權(quán)和.由于DQM權(quán)函數(shù)的選取與特定問題無關(guān),微分方程能微分求積(differential quadrature,DQ)離散為相應的代數(shù)方程求解.下面運用DQM來對數(shù)學模型進行空間離散.

假設半平面占有區(qū)域Ω={〈x,z〉|0≤x≤L0,0≤z≤H0},分別在x、z方向布置Nx×Nz個節(jié)點(見圖2),節(jié)點坐標為Chebyshev-Lobatto多項式的零點.

圖2 半平面上的節(jié)點布置Fig.2 Nodes collocated on the half-plane

根據(jù)DQM,未知函數(shù)ψ(x,z)對x的n階偏導數(shù)可近似表示為

式中:ψkη是在x=xk,z=zη節(jié)點處的函數(shù)值;是n階偏導數(shù)的權(quán)系數(shù),本工作中由Lagrange插值多項式來決定.同時,可給出ψ(x,z)對z的m階偏導數(shù),權(quán)系數(shù)為

在熱局部平衡條件下,利用式(9),xz平面內(nèi)的基本動力學方程(式(1)～(3))的DQ離散化方程為

式中:ξ=2,3,···,Nx-1;η=2,3,···,Nz-1.

本構(gòu)關(guān)系(式(4))、幾何關(guān)系(式(5))和固相的有效應力分量與總應力分量關(guān)系(式(6))的DQ離散化形式為

式中:ξ=1,2,···,Nx;η=1,2,···,Nz.

相應的,邊界條件和對稱性條件的DQ離散化形式為

式中:qξ1,φξ1是半平面表面處給定的外載荷q(x,t)和溫度載荷φ(x,t)的值.

在空間域內(nèi)對半平面問題的數(shù)學模型DQM離散化后,產(chǎn)生了關(guān)于時間t的微分代數(shù)方程組.使用二階向后差分格式來處理關(guān)于微分代數(shù)方程組中的時間導數(shù)

式中:Δt是時間步長;φ(tn,x)是在時刻t=tn的函數(shù)值.

利用式(13),時間離散DQ離散化方程組(10)～(14).最后利用DQ離散的初始條件式(8),對離散化方程組進行Newton-Raphson迭代求解,從而可得到問題的數(shù)值結(jié)果.

3 動力學特性

3.1 驗證方法

忽略體積力的影響,Boer等[3]分析了幾何小變形和常溫條件下不可壓飽和多孔彈性介質(zhì)的一維動力固結(jié)問題,并給出了解析解.為了驗證本方法的正確性、有效性,在給定的數(shù)學模型中設平面表面變溫為φ(x,t)=0,從而可以認為熱局部平衡條件下的解近似等于常溫狀態(tài)時的解.本工作用一個高H0=10 m,寬L0=10 m的流體飽和多孔彈性體來模擬Boer等[6]所考慮的類似問題,并與解析解進行比較.在計算驗證中,分別考慮了階梯載荷q(x,t)=q0h(t)和周期循環(huán)載荷q(x,t)=q0[1-cos(ωt)]的影響,其中q0=3 kN/m2,ω=0.4 s-1,h(t)為Heaviside函數(shù).表1是計算時選取的材料參數(shù).

表1 材料參數(shù)[13]Table 1 Material parameters[13]

圖3分別給出了階梯載荷和周期載荷下對稱軸(x=0)處不同深度的位移uz(沉降),其中時間步長Δt=1 s.從圖3中可以看到:數(shù)值解與解析解吻合良好;當布置結(jié)點數(shù)為7×7時,便能得到令人滿意的結(jié)果.這說明用本方法計算流體飽和多孔彈性介質(zhì)的動力學特性是有效的.

圖3 不同深度的位移時程曲線(x=0)Fig.3 Time-history curves of the displacement in differen depths(x=0)

3.2 參數(shù)和非線性影響

3.2.1 階梯溫度載荷下的動力學特性

考察高H0=10 m,寬L0=10 m的流體飽和多孔熱彈性平面域,邊界條件由式(7)給定.考察表面的階梯溫度載荷φ(x,t)=φ0h(t)對半平面動力學特性的影響,其中φ0=10°C,h(t)是Heaviside函數(shù).計算中布點數(shù)為7×7,時間步長Δt=1 s,平面表面垂直外載q(x,t)=3 kN/m2,材料參數(shù)見表1.

圖4階梯溫度載荷下κF對半平面時-程曲線的影響Fig.4 The effect ofκF on time-history curves of half-plane under the step temperature loading

圖4 為在階梯溫度載荷作用下κF對半平面對稱軸(x=0)處未知量uz,wz,p及θ的影響.圖4(a)表明,隨著時間的逐漸增加,即使在不同的κF下,沉降uz最終也會達到相同的穩(wěn)定狀態(tài).κF較大時,由于在初始階段固結(jié)速度比熱傳導快,導致熱體積膨脹不明顯,所以半平面的沉降uz變大.圖4(b)表明了相對流速wz的變化趨勢.當系統(tǒng)趨于穩(wěn)定時,相對流速wz趨于0;但在初始階段,平面內(nèi)部的流速小于上表面附近,并且κF較大時流速的峰值較大.圖4(c)表明孔隙壓p隨著時間的增加由初始孔壓逐漸至0,并且半平面上表面附近孔隙壓的消散速度較平面內(nèi)部快,同時孔隙壓的消散速度隨κF增大而增大.圖4(d)表明溫度從上表面?zhèn)鞑サ桨肫矫嫔钐?并且逐漸到達至給定的表面溫度,最終達到等溫狀態(tài).同時,κF對變溫θ的影響很小.

圖5為熱傳導系數(shù)k對半平面對稱軸(x=0)處未知量uz,wz,p及θ的影響.由圖5可以看到,在給定的參數(shù)下,當k較大時,初始階段時熱傳導比固結(jié)過程快,半平面的沉降uz較小,溫度場達到等溫狀態(tài)所需的時間更短.

圖5階梯溫度載荷下k對半平面時-程曲線的影響Fig.5 The effect of k on time-history curves of half-plane under the step temperature loading

圖6 為附加熱交換系數(shù)β對半平面對稱軸(x=0)處未知量uz,wz,p及θ的影響.當β較小時,流相和固相之間的耦合作用減弱,初始階段熱體積膨脹不明顯,半平面的沉降uz較大.

圖6 階梯溫度載荷下β對半平面時-程曲線的影響Fig.6 The effect ofβon time-history curves of half-plane under the step temperature loading

圖7為半平面對稱軸不同深度處沉降uz的幾何非線性結(jié)果與相應的線性解的比較.可以看到:當載荷較小時(q(x,t)=3 kN/m2),非線性解和線性解趨于一致;當載荷較大時(q(x,t)=30 kN/m2),非線性解略大于線性解.

圖7 半平面的非線性解與線性解的對比Fig.7 The comparison between nonlinear solutions and linear solutions of half-plane

3.2.2 循環(huán)溫度載荷下的動力學特性

限于篇幅,本工作只給出了循環(huán)溫度載荷下達西滲透系數(shù)κF對系統(tǒng)的影響(見圖8).由圖8可以看出:在循環(huán)溫度載荷和階梯溫度載荷作用下,半平面整體熱動力學性能相似;各物理量呈周期性變化且周期相同,但在深度方向上存在相位差.

圖8 循環(huán)溫度載荷下κF對半平面時-程曲線的影響Fig.8 The effect ofκF on time-history curves of half-plane under the cycle temperature loading

4 結(jié)束語

在幾何非線性條件和熱局部平衡條件下,本工作對流體飽和不可壓多孔熱彈性半平面在受到表面溫度載荷作用下的動力學特性進行了研究.首先基于PMT并考慮幾何非線性的影響,給出了問題的數(shù)學模型;其次提出了一個系統(tǒng)的綜合數(shù)值計算方法來模擬問題的數(shù)值結(jié)果,通過DQM和二階后向差分格式分別在空間域和時間域離散數(shù)學模型,利用Newton-Raphson法求解非線性代數(shù)方程組,從而可得到問題的數(shù)值結(jié)果,分析半平面的動力學特性.

與現(xiàn)有解析解的對比性研究驗證了本方法是正確可靠的,且具有精度高、計算量小、數(shù)值穩(wěn)定等優(yōu)點.在此基礎上,分析了在表面溫度載荷作用下非線性流體飽和多孔熱彈性半平面的熱力學特性,考察了系統(tǒng)的材料參數(shù)達西滲透系數(shù)κF、附加熱交換系數(shù)β、熱傳導系數(shù)k以及幾何非線性對半平面動力學特性的影響,獲得了一些有益的結(jié)論.

(1)在階梯溫度載荷和循環(huán)溫度載荷下,半平面整體熱動力學性能相似:隨著時間增加,半平面沉降趨于穩(wěn)定,流速逐漸趨于0,孔隙壓由初始值逐漸至零,溫度逐漸上升并由上表面向縱深處傳導和擴散.在循環(huán)溫度載荷作用下,各物理量呈周期性變化且周期相同,但在深度方向上存在相位差.

(2)達西滲透系數(shù)較大或附加熱交換系數(shù)較小時,固結(jié)作用過程快于熱傳導過程,初始階半平面體沉降較大.

(3)熱傳導系數(shù)較大時,初始階段時熱傳導比固結(jié)過程快,半平面的沉降較小,溫度場達到等溫狀態(tài)所需的時間更短.

(4)當載荷較小時半平面的非線性解和線性解趨于一致,當載荷較大時,非線性解略大于線性解,因此當載荷較大時,有必要考慮非線性的影響.