飽和多孔介質(zhì)不同動力耦合形式數(shù)值分析

劉 寶， 蘇 謙,2， 劉 亭， 李 婷

(1. 西南交通大學(xué) 土木工程學(xué)院， 成都 610031；2. 西南交通大學(xué) 高速鐵路線路工程教育部重點(diǎn)實(shí)驗(yàn)室， 成都 610031)

飽和多孔介質(zhì)不同動力耦合形式數(shù)值分析

劉 寶1， 蘇 謙1,2， 劉 亭1， 李 婷1

(1. 西南交通大學(xué) 土木工程學(xué)院， 成都 610031；2. 西南交通大學(xué) 高速鐵路線路工程教育部重點(diǎn)實(shí)驗(yàn)室， 成都 610031)

Biot飽和多孔介質(zhì)波動行為的數(shù)值模擬在眾多工程領(lǐng)域中具有重要的意義和作用，由于固相與液相耦合方程難以解耦，使該問題的數(shù)值模擬難度較大。針對飽和多孔介質(zhì)中部分耦合u-p及全耦合u-p-U方程形式的特征，推導(dǎo)了相應(yīng)動力耦合控制方程的有限元弱形式，并引入不同耦合形式的飽和多孔介質(zhì)時域黏性邊界，綜合利用Comsol Multiphysics提供的偏微分方程應(yīng)用模式進(jìn)行二次開發(fā)求解，通過一維飽和多孔介質(zhì)動力響應(yīng)的解析解和數(shù)值解驗(yàn)證了模型求解技術(shù)的合理性和可行性，基于u-p-U耦合形式探討了沖擊荷載作用下干砂飽和砂地基動力固結(jié)中應(yīng)力波傳播特性。計算結(jié)果表明慢縱波對動力固結(jié)的影響比較顯著，合理的沖擊荷載持續(xù)時間有利于固結(jié)效果的改善。

飽和多孔介質(zhì)；Biot模型；Comsol Multiphysics；u-p方程；u-p-U方程

飽和多孔介質(zhì)的波動問題涉及諸多工程領(lǐng)域，特別在巖土工程、海洋工程、水利水電工程等領(lǐng)域中普遍存在。目前，研究飽和多孔介質(zhì)宏觀動力行為的理論主要有：Biot波動理論[1-2]和多孔介質(zhì)理論[3]，兩種理論在工程中已被廣泛應(yīng)用，比較而言，Biot波動理論在不同工程領(lǐng)域中得到了更為廣泛和深入的應(yīng)用。自Biot飽和多孔介質(zhì)波動理論建立以來，國內(nèi)外學(xué)者采用解析或數(shù)值手段對Biot波動方程及其工程應(yīng)用進(jìn)行了廣泛的研究。

飽和多孔介質(zhì)動力方程常用的解析方法多采用積分變換、級數(shù)展開和勢函數(shù)分解等手段將其偏微分控制方程轉(zhuǎn)化為變換域中的常微分方程，再通過逆變換獲得空間域-時域中的解[4-5]，在一維等簡單情形下可通過逆變換以獲得其解析表達(dá)式[6-7]，對于大多數(shù)的研究問題尚需數(shù)值積分實(shí)現(xiàn)逆變換。時剛等[8]運(yùn)用了薄層法對層狀飽和地基動力響應(yīng)進(jìn)行求解，避免了逆變換過程中的數(shù)值積分，在一定程度上提高了計算效率和計算精度。

飽和多孔介質(zhì)波動方程數(shù)值求解常用的方法包括有限元法、有限差分法及邊界元法等[9]，有限元方法由于能夠適應(yīng)任意復(fù)雜幾何形體、邊界條件以及不同的材料模型而被廣泛采用。Schanz[10]對Biot波動方程的有限元數(shù)值分析做了較為深入的研究，根據(jù)不同的簡化假設(shè)，按不同的基本未知量提出了u-p、u-U、u-p-U耦合形式的有限元控制方程，并指出在滲透系數(shù)和荷載頻率較小的情況，可以將流體相對于土骨架運(yùn)動的慣性力忽略不計即u-p簡化形式，但并未明確的給出可忽略流體慣性項(xiàng)的量化條件。既有研究中基于u-p簡化形式采用有限元隱式算法進(jìn)行飽和多孔介質(zhì)波動問題的研究比較多見[11-12]。Zhao等[13]根據(jù)u-U形式采用顯式有限元解法研究了二維飽和地基表面剛性基礎(chǔ)的振動問題，顯式有限元法避免了整體剛度矩陣的組裝，從而提高了計算效率。但u-U形式在固相與液相兩者壓縮性差距較大時易出現(xiàn)體積鎖死問題，而u-p-U形式通過引入壓力場克服了該問題，另外，文獻(xiàn)[14]基于u-p-U控制方程開發(fā)了相應(yīng)的有限元程序，研究了砂土液化的問題，并開展了相關(guān)的室內(nèi)試驗(yàn)進(jìn)行了對比驗(yàn)證[15]，比較而言，目前有關(guān)u-p-U形式的文獻(xiàn)研究報道較少。

綜上可以看出，飽和多孔介質(zhì)不同動力耦合形式的有限元求解多通過自行編制有限元程序?qū)崿F(xiàn)，但有限元程序的編寫涉及飽和多孔介質(zhì)動力單元的開發(fā)、求解算法及合理人工邊界的實(shí)施等關(guān)鍵問題，具有一定的理論深度和難度。當(dāng)前廣泛應(yīng)用的功能強(qiáng)大的商業(yè)有限元軟件Ansys、Abaqus和Comsol Multiphysics等均未提供飽和多孔介質(zhì)不同耦合形式的動力單元，可通過其提供的開放接口編寫內(nèi)嵌語言或子程序?qū)崿F(xiàn)，而目前這方面的研究報道尚不多見。

鑒于Comsol Multiphysics在定義和耦合偏微分方程靈活方便的優(yōu)勢，本文通過積分變換得到飽和多孔介質(zhì)動力控制方程u-p及u-p-U的弱形式，引入Akiyoshi[16]給出的不同耦合形式的飽和多孔介質(zhì)時域黏性邊界，綜合利用該軟件提供的偏微分方程應(yīng)用模式進(jìn)行求解，通過一維飽和多孔介質(zhì)動力響應(yīng)的解析解驗(yàn)證了不同耦合形式模型的正確性，在此基礎(chǔ)上，基于u-p-U耦合形式建立干砂-飽和砂地基耦合模型，并利用該模型著重分析沖擊荷載下該體系的應(yīng)力波傳播特性，為研究飽和多孔介質(zhì)動力響應(yīng)問題提供了一種方便快捷的計算模式。根據(jù)該計算模式，可靈活地對研究中新提出的應(yīng)力場控制方程、滲流場控制方程及固液相耦合關(guān)系發(fā)生變化的飽和多孔介質(zhì)波動方程實(shí)現(xiàn)計算分析，以滿足實(shí)際問題的需要。

1 飽和多孔介質(zhì)動力控制方程

1.1 全耦合u-p-U方程

采用文獻(xiàn)[17]給出的基本方程，考慮土顆粒的壓縮性

土骨架本構(gòu)方程

(1)

整體平衡方程

(2)

流體平衡方程

(3)

滲流連續(xù)方程

(4)

將式(3)代入式(4)中得

(5)

其中流體位移未知量Ui=ui+wi/n，式(2)、式(5)與式(1)構(gòu)成了飽和多孔介質(zhì)u-p-U動力耦合形式的方程。

在式(1)中：εij為固相骨架應(yīng)變；λs、μs為固相骨架Lame常數(shù)；p為孔隙水壓；α為Biot常數(shù)，α=1-Kb/Ks，Kb,Ks分別為固相骨架、顆粒的體積模量，Kb=λs+2μs/3；δij為Kronecker符號。在式(2)中：ρm=nρf+(1-n)ρs為總體有效密度；ρs和ρf分別為固相骨架密度和流體密度；n為孔隙率；u為固相骨架位移。在式(3)中：wi為流體相對于固相骨架的位移，在實(shí)際應(yīng)用中根據(jù)；kf為動力滲透系數(shù)。在式(4)中，Q為流體的壓縮系數(shù)，Q=(α-n)/Ks+n/Kf，Kf為流體體積模量。

1.2 部分耦合u-p方程

Zienkiewicz等通過略去Biot方程中的孔隙流體相對加速度項(xiàng)，推導(dǎo)出以土骨架位移u和孔壓p為變量的簡化控制方程。

平衡方程

(6)

達(dá)西運(yùn)動方程

(7)

滲流連續(xù)方程

(8)

將式(7)代入式(8)中得

(9)

式(6)、式(9)與式(1)構(gòu)成了飽和多孔介質(zhì)u-p動力耦合形式的方程。

2 有限元模型的建立

2.1u-p方程的有限元弱形式

Comsol Multiphysics采用虛位移原理對偏微分方程表達(dá)的問題進(jìn)行描述，因此需根據(jù)虛位移原理推導(dǎo)耦合場控制方程的等效弱形式。在式(2)與式(5)中分別乘以虛變量δu和δp，然后在求解區(qū)域范圍內(nèi)對方程進(jìn)行積分，為使積分方程能夠引入邊界條件，應(yīng)用Green-Gauss公式對積分方程進(jìn)行分部積分，從而可得到耦合方程的弱形式。式(2)的積分形式為

(10)

由Green-Gauss公式得

(11)

式(5)的積分形式為

(12)

同理由Green-Gauss公式得

▽Tp(kij▽p)dΩ-

(13)

利用Galerkin加權(quán)余量法離散得到以土骨架位移u和孔壓p為變量的u-p形式控制方程

u(x,t)=NuU(t);p(x,t)=NpP(t)

(14)

δu(x,t)=Nu;δp(x,t)=Np

(15)

式中，Nu，Np分別為U(t)，P(t)的插值形函數(shù)，以土骨架位移u和孔壓p表達(dá)的u-p形式的有限元控制方程為

(16)

(17)

其中，

2.2u-p-U方程的有限元弱形式

采用相同的方式可得出u-p-U方程的有限元弱形式

(18)

其中，

2.3 黏性人工邊界

飽和多孔介質(zhì)中由于固相和液相間存在相互作用，包含兩類壓縮波和一類剪切波，邊界的處理問題變得較為復(fù)雜[18]，飽和多孔介質(zhì)的人工邊界， 除了要處理位移邊界條件外， 還要考慮孔壓在人工邊界的傳播。目前針對飽和多孔介質(zhì)的人工邊界的建立和研究已取得諸多成果，用的比較多的是透射邊界、 黏性邊界和在黏性邊界基礎(chǔ)上發(fā)展的黏彈性邊界[19]。本文采用Akiyoshi飽和多孔介質(zhì)黏性邊界，通過設(shè)置阻尼器來吸收傳到邊界上的波動能量，從而消除或減小在邊界上的反射。

(1)u-p耦合形式

法向應(yīng)力

(19)

切向應(yīng)力

(20)

邊界壓力

(21)

(2)u-p-U耦合形式

法向應(yīng)力

(22)

切向應(yīng)力

(23)

邊界壓力

(24)

Vp為壓縮波在飽和多孔介質(zhì)中的傳播速度

(25)

Vs為剪切波在飽和多孔介質(zhì)中的傳播速度

(26)

2.4 模型建立的實(shí)施

Comsol是一款基于求解偏微分方程組的多場耦合計算平臺，具有較強(qiáng)的計算性能及多場直接耦合的能力，其中PDE模塊可以靈活地建立數(shù)學(xué)模型并基于該方程對問題進(jìn)行求解，現(xiàn)將模型建立的關(guān)鍵流程進(jìn)行描述：① 在模型向?qū)е袆?chuàng)建幾何模型維度，通過MathematicsPDE Interfaces添加Weak Form PDE物理場，并將u-p或u-p-U動力耦合方程中對應(yīng)的待求變量輸入，選擇Time Dependent后進(jìn)入主圖形界面；② 在模型樹中，按照u-p或u-p-U中應(yīng)力場、滲流場控制方程形式定義表達(dá)式，根據(jù)軟件約定的弱形式書寫規(guī)范將定義的方程表達(dá)式在”Weak Expressions”欄中輸入；③ 建立幾何模型并添加相應(yīng)的人工邊界，位移邊界可通過PDE模式中內(nèi)置的Dirichlet Boundary Condition(狄式邊界)添加，應(yīng)力邊界通過添加Weak Contribution來實(shí)現(xiàn)；④ 模型網(wǎng)格劃分與求解，網(wǎng)格單元默認(rèn)為拉格朗日二次單元，網(wǎng)格的大小及時間步長選取應(yīng)滿足CFL穩(wěn)定性原則，模型求解由內(nèi)置的直接求解器或迭代求解器自動完成，即可實(shí)現(xiàn)飽和多孔介質(zhì)動力耦合有限元模型的建立。

3 模型的驗(yàn)證及應(yīng)用

Boer等給出了一維不可壓縮飽和多孔介質(zhì)動力響應(yīng)的解析解，本文利用此解析解進(jìn)行對比驗(yàn)證，算例的計算模型，如圖1所示，模型的橫向尺寸取1 m，豎向尺寸取4 m，按照上文中模型實(shí)施流程進(jìn)行建模，邊界條件分兩種工況進(jìn)行設(shè)置；第一類工況為底邊為固定且不透水邊界，左右邊界水平位移固定且不發(fā)生滲流，上表面為自由透水邊界，即孔壓為0；第二類工況為底部為黏性吸收邊界，其余邊界條件同上，位移邊界通過內(nèi)置的狄式邊界添加，應(yīng)力邊界通過在Weak Contribution欄輸入對應(yīng)的表達(dá)式實(shí)現(xiàn)。采用四邊形拉格朗日二次單元劃分網(wǎng)格，每個節(jié)點(diǎn)包含3個自由度。材料參數(shù)取值，如表1所示，土柱承受的動力荷載的時程曲線為函數(shù)f(t)，如圖1所示。計算時間步長取0.001 s，荷載持續(xù)時間為0.3 s。本研究中由于求解問題的規(guī)模不太大，選用精度高、穩(wěn)定性好、高效的稀疏直接求解器MUMPS進(jìn)行求解，時間積分方法采用general-α方法，該方法具有無條件穩(wěn)定、二階精度的優(yōu)點(diǎn)，并且能夠通過控制數(shù)值阻尼，有效濾除高頻數(shù)值頻散響應(yīng)[20]。

(a)幾何模型f(t)=1-cos(20πt)kPa(b)動力荷載

表1 文獻(xiàn)[6]材料參數(shù)

3.1u-p及u-p-U模型的驗(yàn)證

在u-p耦合模型中，基本未知量為固相位移和孔隙水壓，圖2表示數(shù)值解與已知解析解的對比結(jié)果圖，由圖中可知，采用固定邊界時，開始階段的數(shù)值結(jié)果與解析解結(jié)果幾乎一致，隨著時間的推移，應(yīng)力波傳播到底部時，在固定邊界產(chǎn)生反射波，進(jìn)而使應(yīng)力波在上部和下部邊界產(chǎn)生多次反射，造成計算結(jié)果的偏差越來越大。采用Akiyoshi黏性吸收邊界時，數(shù)值計算結(jié)果與解析解的結(jié)果非常吻合，可以看出，黏性邊界比固定邊界的計算結(jié)果更加合理。

在u-p-U耦合模型中，基本未知量為固相位移、孔隙水壓及液相位移，能夠考慮液相對于固相相對加速度的影響。由圖3可知，采用黏性邊界能夠減少應(yīng)力波在底部邊界的反射，從而使數(shù)值計算與解析解結(jié)果相符合，而采用固定邊界時，由于底部產(chǎn)生的反射波的影響，使得數(shù)值結(jié)果與解析解結(jié)果隨著時間的推移偏差越來越明顯。由圖3(a)、圖3(b)可知，土柱在動荷載作用下固相位移與液相位移的運(yùn)動方向相反，固相向下運(yùn)動，液相向上運(yùn)動從而不斷地擠出孔隙，反映出飽和多孔介質(zhì)中固相與液相由內(nèi)摩擦引起的黏彈性質(zhì)。

(a) z=1 m處固相位移時程曲線

(b) t=0.15 s時超孔隙水壓沿深度分布曲線

(a) z=1 m處固相位移時程曲線

(b) z=1 m處液相位移時程曲線

3.2 模型耦合效應(yīng)討論

在飽和多孔介質(zhì)中，孔隙水與土骨架的相互作用主要包括滲透力和慣性耦合力，滲透系數(shù)、荷載頻率等因素是影響土骨架與孔隙水耦合效應(yīng)的主要參數(shù)，由于在u-p簡化形式中忽略了流體對于固相的相對加速度，因此有必要衡量不同荷載頻率、不同滲透系數(shù)下u-p及u-p-U動力耦合形式的差異性。

圖4為不同滲透系數(shù)下兩種耦合形式的固相位移差變化規(guī)律，可以看出，在滲透系數(shù)較小時，孔隙水和固相骨架之間的相對位移比較小，產(chǎn)生的慣性耦合效應(yīng)比較弱，在此情形下，u-p及u-p-U兩種耦合形式的計算結(jié)果沒有差異，可以認(rèn)為滲透系數(shù)kf<10-3m/s時，慣性耦合效應(yīng)可以忽略。但在滲透系數(shù)較大時，孔隙水和固相骨架之間的相對流動阻力小，兩者之間形成較大的相對位移，從而表現(xiàn)出明顯的慣性耦合效應(yīng)，使得兩者之間的計算結(jié)果出現(xiàn)差異。為考察荷載頻率對慣性耦合效應(yīng)的影響。圖5為不同荷載頻率下兩種耦合形式的固相位移差變化規(guī)律。選取kf=10-6m/s以忽略滲透系數(shù)對慣性效應(yīng)的影響，可以看出，隨著頻率的不斷提高兩種耦合形式的計算結(jié)果差異越來越顯著，頻率對飽和多孔介質(zhì)應(yīng)力波的傳播特性有較大的影響，從上述分析可知，部分耦合形式u-p適用于滲透系數(shù)和荷載頻率較小的情況。

圖4 不同滲透率下兩種耦合形式的固相位移差曲線

Fig.4 Solid displacement of the two dynamic formulation in differernt permeability coefficients

圖5 不同頻率下兩種耦合形式的固相位移差曲線

Fig.5 Solid displacement of the two dynamic formulation in differernt frequencies

3.3 應(yīng)用算例

濱海地區(qū)粉細(xì)砂地基分布廣泛，且地下水位較高，粉細(xì)砂地基受到地下水影響常處于飽水狀態(tài)，在工程實(shí)踐中多應(yīng)用動力固結(jié)法進(jìn)行處理，鑒于目前有關(guān)u-p-U耦合形式的應(yīng)用實(shí)例較少，本節(jié)基于此動力形式建立干砂-飽和砂層狀地基耦合模型，研究沖擊荷載作用下該體系的應(yīng)力波傳播特性。計算模型，如圖6所示，飽和砂層上表面為自由透水邊界，底邊為固定且不透水邊界，左右邊界水平位移固定，材料計算參數(shù)見表2，荷載形式為三角形沖擊荷載，加載周期分別為0.01 s和1 s，峰值都為1 Mpa，時間積分方法采用general-α方法，計算時間步長取0.000 1 s，采用直接求解器MUMPS進(jìn)行求解。

(a)幾何模型示意圖(b)T/s=0.01s(c)T/s=1.0s

表2 模型材料參數(shù)

在干砂中沖擊荷載激發(fā)的應(yīng)力縱波傳播速度為197 m/s，荷載應(yīng)力峰值在12 ms時到達(dá)飽和砂區(qū)域，圖7表示此時飽和砂中有效應(yīng)力與超孔隙水壓力的云圖，由圖7可知，在干砂中產(chǎn)生的應(yīng)力波在飽和砂上部形成了超孔隙水壓較高的區(qū)域，并表現(xiàn)出了曼德爾效應(yīng)，即自由排水的上表面受到荷載作用時，會發(fā)生固結(jié)，土中的水排出，孔壓消散，但是孔壓不會立即消散，而是先上升后才慢慢消散。

由圖8可知，沖擊荷載產(chǎn)生的能量以應(yīng)力波的形式向體系內(nèi)部傳播，通過干砂層與飽和砂層分界面時由于介質(zhì)阻抗的不同發(fā)生了明顯的反射和透射現(xiàn)象，在飽和中快縱波和慢縱波波速大致為1 580 m/s和1 m/s，超孔隙水壓的波動特性與快縱波的傳播規(guī)律一致，而慢縱波波速較小且急劇衰減。由此可知，在該情形下，沖擊荷載產(chǎn)生的應(yīng)力波到達(dá)飽和砂層時有效應(yīng)力迅速減少，轉(zhuǎn)換為超孔隙水壓以快縱波的形式進(jìn)行傳播，由較高的區(qū)域向深部和水平向發(fā)展，圖9表示相同滲透率下荷載持續(xù)時間為1 s時，豎向有效應(yīng)力與超孔隙水壓力的分布圖。由圖可知，該情形下飽和砂中形成了明顯的慢縱波，超孔隙水壓峰值的傳播速度大約為8 m/s，與飽和砂慢縱波的傳播速度一致。

(a) 豎向有效應(yīng)力

(b) 超孔隙水壓

(a) 豎向有效應(yīng)力

(b) 超孔隙水壓

Fig.8 Normalized vertical effective stress and excess pore pressure profiles for loading duration of 0.01 s atx=1 m

由上述分析可知，不同荷載加載速率下，應(yīng)力波在該體系中的傳播規(guī)律表現(xiàn)出顯著的差異性。荷載加載速率過快，荷載在干砂中產(chǎn)生的縱波到達(dá)飽和砂層時轉(zhuǎn)化為超孔隙水壓以快縱波形式進(jìn)行傳播，慢縱波發(fā)揮的作用不大，由于流體的壓縮模量遠(yuǎn)大于土骨架的壓縮模量，從而導(dǎo)致飽和砂中有效應(yīng)力的增加幅度甚微，即快縱波對土體固結(jié)的影響不明顯。當(dāng)荷載加載速率減少時，形成了明顯的慢縱波使得土骨架壓縮從而引起孔隙水排出，但如果滲透系數(shù)較小，土骨架與孔隙水的耦合效應(yīng)使得兩者同時被壓縮，由于流體的壓縮模量遠(yuǎn)大于土骨架的壓縮模量使得荷載轉(zhuǎn)換為超孔隙水壓的比例升高，因此，為促使飽和砂中有效應(yīng)力的增加，宜提高沖擊荷載的周期加載時間和土體滲透性，例如在強(qiáng)夯法加固地基工程實(shí)踐中，采用中粗砂作為強(qiáng)夯墊層。目前，動力排水固結(jié)法加固軟基的加固效果及其影響因素的研究仍待進(jìn)一步深入開展，本文建立的數(shù)值研究方法可為動力排水固結(jié)法加固機(jī)制的研究提供一種方便快捷的手段。

(a) 豎向有效應(yīng)力

Fig.9 Normalized vertical effective stress and excess pore pressure profiles for loading duration of 1 s atx=1 m

4 結(jié) 論

(1)針對u-p及u-p-U飽和多孔介質(zhì)動力方程，借助COMSOL Multiphysics偏微分方程應(yīng)用模式，實(shí)現(xiàn)了飽和多孔介質(zhì)兩種動力耦合方程的有限元求解，為飽和多孔介質(zhì)動力響應(yīng)問題提供了一種方便快捷的計算模式。

(2)通過一維飽和多孔介質(zhì)動力響應(yīng)解析解和本文數(shù)值解的對比分析，驗(yàn)證了模型求解技術(shù)的合理性和可行性，并便捷地實(shí)現(xiàn)了飽和多孔介質(zhì)黏性邊界的應(yīng)用和施加，可方便地用于飽和多孔介質(zhì)無限域的動力計算，兩種動力耦合形式的計算結(jié)果表明部分耦合形式u-p適用于滲透系數(shù)和荷載頻率較小的情況。

(3)基于u-p-U耦合形式探討了沖擊荷載作用下干砂-飽和砂地基動力固結(jié)中應(yīng)力波的傳播特性，計算結(jié)果表明慢縱波對動力固結(jié)的影響比較顯著，合理的沖擊荷載持續(xù)時間有利于固結(jié)效果的改善。

[1] BIOT M A. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range[J]. Journal of the Acoustical Society of America, 1956, 28(2): 168-178.

[2] BIOT M A. Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range[J]. Journal of the Acoustical Society of America, 1956, 28(2): 179-191.

[3] BOWEN R M. Compressible porous media models by use of the theory of mixtures[J]. International Journal of Engineering Science, 1982, 20(6): 697-735.

[4] THEODORAKOPOULOS D D. Dynamic analysis of a poroelastic half-plane soil medium under moving loads[J]. Soil Dynamics and Earthquake Engineering, 2003, 23(7): 521-533.

[5] 胡安峰, 孫波, 謝康和. 下臥基巖飽和地基在移動荷載作用下的動力響應(yīng)[J]. 振動與沖擊, 2012, 31(4): 151-156.

HU Anfeng, SUN Bo, XIE Kanghe. Dynamic response of a saturated subgrade with rock substratum to moving loads[J]. Journal of Vibration and Shock, 2012, 31(4): 151-156.

[6] BOER D I R D, EHLERS D I W, LIU Z. One-dimensional transient wave propagation in fluid-saturated incompressible porous media[J]. Archive of Applied Mechanics, 1993, 63(1): 59-72.

[7] GAJO A, MONGIOVI L. An analytical solution for the transient response of saturated linear elastic porous media[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1995, 19(6): 399-413.

[8] 時剛, 高廣運(yùn), 馮世進(jìn). 二維及軸對稱條件下飽和層狀地基的Lamb問題解答[J]. 土木工程學(xué)報,2010,43(8):125-131.

SHI Gang, GAO Guangyun, FENG Shijin. The Green’s function of 2D and axisymmetric saturated layered half-space soil[J]. China Civil Engineering Journal,2010,43(8):125-131.

[9] 黃茂松, 李進(jìn)軍. 飽和多孔介質(zhì)土動力學(xué)理論與數(shù)值解法[J]. 同濟(jì)大學(xué)學(xué)報(自然科學(xué)版), 2004, 32(7): 851-856.

HUANG Maosong, LI Jinjun. Dynamics of fluidsaturated porous media and its numerical solution[J]. Journal of Tongji University (Natural Science), 2004,32(7): 851-856.

[10] SCHANZ I M, CHENG A H D. Transient wave propagation in a one-dimensional poroelastic column[J]. Acta Mechanica, 2000, 145(1): 1-18.

[11] RAHMANI A, PAK A. Dynamic behavior of pile foundations under cyclic loading in liquefiable soils[J]. Computers and Geotechnics, 2012, 40(3): 114-126.

[12] KHOSROJERDI M, PAK A. Numerical investigation on the behavior of the gravity waterfront structures under earthquake loading[J]. Ocean Engineering, 2015, 106: 152-160.

[13] ZHAO C, LI W, WANG J. An explicit finite element method for Biot dynamic formulation in fluid-saturated porous media and its application to a rigid foundation[J]. Journal of Sound and Vibration, 2005, 282(3): 1169-1181.

[15] TASIOPOULOULA P, TAIEBAT M, TAFAZZOLI N, et al. On validation of fully coupled behavior of porous media using centrifuge test results[J]. Coupled Systems Mechanics, 2015,4(1):37-65.

[16] AKIYOSHI T, FUCHIDA K, FANG H L. Absorbing boundary conditions for dynamic analysis of fluid-saturated porous media[J]. Soil Dynamics and Earthquake Engineering, 1994, 13(6): 387-397.

[17] ZIENKIEWICZ O C, SHIOMI T. Dynamic behaviour of saturated porous media; the generalized Biot formulation and its numerical solution[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1984, 8(1): 71-96.

[18] 劉光磊, 宋二祥. 飽和無限地基數(shù)值模擬的黏彈性傳輸邊界[J]. 巖土工程學(xué)報, 2006, 28(12): 2128-2133.

LIU Guanglei, SONG Erxiang. Visco-elastic transmitting boundary for numerical analysis of infinite saturated soil foundation[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(12): 2128-2133.

[19] LI P, SONG E. A high-order time-domain transmitting boundary for cylindrical wave propagation problems in unbounded saturated poroelastic media[J]. Soil Dynamics and Earthquake Engineering, 2013, 48(5): 48-62.

[20] HAN B, ZDRAVKOVIC L, KONTOE S. Stability investigation of the Generalised-α time integration method for dynamic coupled consolidation analysis[J]. Computers and Geotechnics, 2015, 64: 83-95.

Numerical simulation for different dynamic coupling forms of saturated porous media

LIU Bao1, SU Qian1,2, LIU Ting1, LI Ting1

(1. School of Civil Engineering, Southwest Jiaotong University, Chengdu 610031, China;2. MOE Key Laboratory of High Speed Railway Engineering, Southwest Jiaotong Univeristy, Chengdu 610031, China)

Numerical simulation of dynamic behavior of saturated porous media is of great importance in many engineering problems. As pore fluid and solid skeleton interaction is hard to be decoupled, there are a lot of difficulties in numerical simulation. Aiming at characteristics of the partly coupled dynamic field equationu-pformulation and the fully coupled dynamic field equationu-p-Uformulation, the FE weak forms of the corresponding dynamic coupling control equations were derived. These FE weak forms were successfully implemented in the finite element software Comsol Multiphysics and the time domain viscous boundarys with different dynamic coupling forms were introduced in simulating unbounded saturated porous media domain. The reasonability and feassibility of the model solving technique were verified by using the analytical solution and numerical one to dynamic response of one-dimensional saturated porous meida. Finally, based onu-p-Ucoapled form, the propagating features of stress waves in the dry sandsaturated sand foundation dynamic consolidation under impacting loads were investigated. The computation results showed that the slow longitudinal were has an obvious effect on the dynamic consolidation; the reasonable time duration of impacting load is beneficial to the improvement of consolidation effects.

saturated porous media; Biot model; Comsol Multiphysics;u-pequation;u-p-Uequation

國家自然基金項(xiàng)目(51378441；51578467)

2015-12-30 修改稿收到日期：2016-03-01

劉寶 男，博士生，1988年生

蘇謙 男，教授，博士生導(dǎo)師，1972年生 E-mail:suqian@126.com

TU443

A

10.13465/j.cnki.jvs.2017.09.022