馬懷發(fā)，楊茹

（1.中國水利水電科學研究院流域水循環(huán)模擬與調控國家重點實驗室，北京100038；2.中國水利水電科學研究院工程抗震中心，北京100048）

三維飽和地基人工黏彈性邊界

馬懷發(fā)1，2，楊茹2

（1.中國水利水電科學研究院流域水循環(huán)模擬與調控國家重點實驗室，北京100038；2.中國水利水電科學研究院工程抗震中心，北京100048）

本文針對中低頻率工程振動問題，不計流體加速度，基于球面波動方程導出了無限大三維飽和地基人工黏彈性邊界及其流量邊界條件。在此基礎上給出了具有人工黏彈性邊界的飽和土體動力固結問題虛位移原理。根據有限元語言編寫出算法腳本文件，并基于有限元自動生成系統(tǒng)（FEPG），生成了計算源代碼程序。通過在無限大地基表面施加沖擊荷載和突加荷載的數(shù)值算例，驗證了所建立方程和程序的正確性。本文方法可以應用于三維飽和地基的動力固結問題數(shù)值模擬。

半無限飽和地基；三維人工黏彈性邊界；虛位移原理；孔隙水壓力；有效應力

1 研究背景

根據地震的傳播機制，地震荷載是從遠域地基波及到土體結構，而后產生地震響應。地震波由遠及近，再由近及遠傳播。在采用有限元法進行動力分析時，一般將半無限地基劃分為鄰近結構的近域地基和其外圍的遠域地基。近域地基及土體結構構成了地震響應分析的研究對象，而切割近域和遠域地基的邊界即形成了人工邊界。人工邊界分為全局邊界和局部邊界，全局邊界雖然是對無限介質的精確模擬，但其時空偶聯(lián)的頻域求解方法難以在有限元中實現(xiàn)，非線性問題不好考慮。局部人工邊界是時空解耦的，并且易于實現(xiàn)、計算量小，因而得到了廣泛的研究和應用。目前，用的比較多的是人工透射邊界［1］、黏性邊界［2］和在黏性邊界基礎上發(fā)展的黏彈性邊界［3-7］。黏性邊界因其概念清楚、易于應用的特點，一度得到了廣泛的研究和應用，但其存在低頻穩(wěn)定性問題，而且從一維推廣到多維有很大的誤差。黏彈性邊界是在黏性邊界的基礎上增加了彈簧，即通過沿人工邊界設置一系列線性彈簧和阻尼器組成的簡單力學模型來吸收射向人工邊界的波動能量，從而達到消除反射的效果，并能模擬地基的彈性恢復能力，可以很好地模擬波動在無限地基中的真實傳播過程。Deeks［3］基于柱面波給出了二維黏彈性人工邊界的應力公式，劉晶波［5-6］等人基于球面波推導了三維黏彈性人工邊界，并給出了邊界的等效剛度和阻尼系數(shù)，馬懷發(fā)等［7］在分析黏彈性邊界方法的理論基礎上建立了統(tǒng)一的虛位移方程，郝明輝［8］等將其應用于拱壩計算。

強震作用下土壩、土體邊坡及其地基會產生液化失穩(wěn)問題，因此對其進行抗震安全評價理應將土體看作由土骨架和孔隙水兩相材料組成。動力分析不僅要給出土骨架的位移時程還要計算出孔隙水壓力的變化時程。與單向固體結構的人工黏彈性邊界［3-7］不同，對于飽和土體的人工邊界，除了要處理位移邊界條件外，還要考慮孔隙水壓力在人工邊界的傳播。Modaressi等［9］基于簡化的比奧（Biot）方程，針對實際土木工程的中低頻率振動并且滲透系數(shù)比較小的情況，忽略第二類壓縮波，將傍軸近似應用于兩相介質，提出了飽和介質動力方程u-p形式的黏性邊界。Akiyoshi等［10］進一步給出了u-w和u-U形式的飽和介質時域黏性邊界。王子輝等［11］基于u-U形式分別給出了具有輻射阻尼性質的外行柱面波和球面波在圓柱面和球面人工邊界上引起的法向、切向應力的表達式，同時模擬了二維半空間無限域介質的能量吸收作用。劉光磊［12］基于柱面波給出了二維飽和地基u-p形式黏彈性邊界條件。本文將基于比奧方程的u-p形式，研究三維飽和地基黏彈性邊界條件，基于球面波給出黏彈性邊界的法向和切向的彈簧系數(shù)、阻尼系數(shù)以及流量邊界條件，并建立具有黏彈性邊界的三維飽和地基動力學的弱解形式。

2 基本方程

2.1 動力固結方程基于比奧理論，Zienkiewicz給出了動力固結基本控制方程［13］，若規(guī)定土體中應力以受壓為正，受拉為負，動力固結控制方程為：

式中：σij為飽和土的總應力；為有效應力；p為孔隙水壓力；n為孔隙率；ui為土骨架位移；wi為水相對于土骨架的位移，t定 義為相對于總截面面積的平均滲流速度；εij為土骨架的應變，分別為飽和土、土和水的密度，為動力滲透系數(shù)，與土力學中滲透系數(shù)k的關系為為水的容重；λ，μ為土骨架的拉梅（Lame）常數(shù)α、Q為與固體和流體的壓縮性相關的系數(shù)，，其中Ks、Kf和Kb分別為土顆粒、流體和土骨架的體積模量；θ為混合體體積應變；bi為重力加速度。

2.2 動力固結方程的弱形式Zienkiewicz研究表明［13］，對于包括地震工程在內的大部分中低振動頻率的工程問題，中等速度運動情況下可忽略流體加速度，同時忽略土骨架相對于流體的加速度，動力固結方程變?yōu)槿缦滦问剑?/p>

將上面的動力固結控制方程寫成弱解形式為

式中：Γ為應力邊界，ni、nj為對應邊界的方向余弦。

3 飽和地基黏彈性邊界

3.1 黏彈性邊界條件將結構和近場作為動力平衡系統(tǒng)如圖1所示，由入射波產生的波動場由{} u表示，阻尼[] C和剛度[]

K反映結構和地基的動力特性。這里及下文物理量下腳標“b”表示人工邊界，“s”表示結構。假定該系統(tǒng)中一散射波源以球面波動的形式向人工邊界傳播，基于球面波理論建立黏彈性邊界條件。設球坐標系為（r，φ，θ），波動問題可視為球對稱問題，因此所有力學變量只和r有關，分析問題時可只考慮徑向r和垂直于徑向的兩個切線方向φ和θ。由于問題的對稱性，位移：應變：；應力，式中：ur為徑向位移，εr為徑向正應變，εθ和ε?為切向正應變，σr為球面的徑向正應力，σθ和 σ?為切向正應力。幾何方程：

可忽略流體加速度和土骨架加速度后，在球坐標下將式（2）兩邊同時對r求導得：

圖1 結構-地基動力平衡系統(tǒng)

再將（4）式兩邊同時對t求導，并代入式（7）可得球坐標系下孔壓和位移的關系式：

當滲透系數(shù)比較小時，可假設滲透系數(shù)k=0，因此kf=0。研究表明，當滲透系數(shù)比較小時，此假設所引起的誤差可以忽略［13-14］。則式（8）可表示為：

由此可導出用位移表示的波動方程：

令

式中：Vp為膨脹玻（P波）在飽和介質中的傳播速度。引入位移勢函數(shù)

式（10）可表示為：

則方程（13）的解為：

球面波陣面上法向應力和位移滿足：

在?黏彈性人工邊界的物理模型中［5］，黏彈性人工邊界節(jié)點上的應力與位移滿足的微分方程：

對比式（15）和式（16）可得黏彈性人工邊界法向應力邊界條件

孔壓只對法向應力有影響，剪切方向有效應力與總應力相等，則剪切應力與位移滿足微分方程：

式中Vs為飽和介質中剪切波波速，

由式（18）可得黏彈性人工邊界切向應力邊界條件

由以上參數(shù)可以發(fā)現(xiàn)，三維飽和地基黏彈性邊界的參數(shù)形式和一般的三維黏彈性邊界的參數(shù)形式一樣，不同的是波速考慮了孔隙水的影響。

3.2 黏彈性邊界流量假設動力計算過程飽和土體的滲流符合達西定律，流量邊界條件為：

將式（9）代入得：

將位移勢函數(shù)式（12）代入得邊界法向流量表達式為：

可以看出，雖然二維和三維的控制方程形式不同，波動理論也不同，但是邊界流量公式的形式［12］最終是相同的。

4 人工黏彈性邊界的虛位移原理

4.1 方程的弱解形式根據文獻［7］給出的單相體結構的黏彈性邊界虛位移原理，當求解的數(shù)值模型采用人工黏彈性邊界時，三維飽和地基動力固結控制方程及其邊界條件可以寫成如下統(tǒng)一的形式：

式中：Fˉmb作用在人工邊界上地震波動轉化荷載；為黏彈性人工邊界上的等效剛度系數(shù)；為黏彈性人工邊界上的等效阻尼系數(shù)；Svs為黏彈性邊界，假定散射源到人工邊界的距離rb可近似地取表面中心到人工邊界的距離。上角標：m=±x，±y，-z；下腳標：i，j，k=1，2，3。其中，

這里，αN和αT取值范圍［15］分別為1～2.0和0.5～1。本文αN和αT分別取1.2、0.5。

εt

方程（26）是動力固結方程（5）在引入人工黏彈性邊界后的積分弱解形式。

4.2 算法及FEPG腳本文件利用方程（26）可建立其有限元數(shù)值離散模型，聯(lián)立求解土骨架變形和孔隙水壓力。關于土骨架位移ui的動力學方程是雙曲線型方程，而關于孔隙水壓力p的滲流方程是橢圓型方程。在常規(guī)的邊界上，求解動力學方程需要給出位移邊界和應力（荷載）邊界條件。求解滲流方程需要給出壓力邊界條件。壓力第二類邊界條件由邊界法向加速度表示流量表達式（22）給出。動力學方程中的速度、加速度采用紐馬克（NEWMARK）算法進行插值離散，在空間上進行有限元離散。離散后的式（26）壓力變量僅在剛度矩陣里與位移變量耦合。

本文計算程序是基于有限元程序自動生成系統(tǒng)（Finite Element Program Generator，F(xiàn)EPG）［16］開發(fā)出來的。按照該系統(tǒng)提供有限元語言規(guī)則［17］，根據本文提出的飽和地基人工黏彈性邊界的虛位移原理式（26）編寫腳本。腳本文件包括：PDE文件，F(xiàn)BC文件，算法（NFE，GCN和MDI）文件。體積分項由PDE表達，面積分項由FBC表達，NFE描述NEWMARK算法的廣義剛度矩陣的組合關系及其算法過程，GCN文件給出求解線性方程組的求解器及時間計算流程，MDI文件給出求解未知量及PDE、FBC所對應的單元類型及高斯積分點個數(shù)。然后運行MDI命令即可生成計算源代碼程序。

5 數(shù)值算例

為了驗證本文給出的飽和地基黏彈性邊界及其動力學虛位移方程的可靠度和準確性，將通過兩個數(shù)值算例來驗證。算例通過黏彈性邊界結果和參考邊界結果的對比驗證了黏彈性邊界的可靠性，并和固定邊界結果對比可以看出黏彈性邊界的計算結果更加合理。黏彈性邊界、固定邊界計算模型如圖2所示，模型尺寸為48 m×48 m×24 m。土體表面為自由排水面，在表面中心8 m×8 m的區(qū)域加面荷載。計算邊界采用黏彈性邊界時，土體上表面自由，底面和4個側面加黏彈性邊界。固定邊界計算時，土體上表面自由，前后面y向約束不透水，左右面x向約束不透水，底面z向約束不透水。參考結果是將計算區(qū)域擴大，長寬高分別變?yōu)樵Ｐ偷?倍，即144 m×144 m×72 m，在計算時段內計算結果不受邊界的影響，可以認為該計算結果為真實解。經計算，參考解的邊界采用固定邊界和黏彈性邊界算出的結果是一致的，本文參考解計算采用的是黏彈性邊界。

圖2 飽和地基計算模型

圖3 沖擊荷載

計算參數(shù)與文獻［10］的參數(shù)取值一致，即λ=83.3kPa；μ=125kPa；Ks=1.0×108kPa；Kf=1.0×103kPa；ρ=1.7×103kg/m3；ρf=1.0×103kg/m3；k=1.0×10-2m/s；kf=1.0×10-6m4/(N?s)；n=0.3；α=1.0由計算得Q=3.33 MPa，ν=0.2，E=300 kPa，Kb=166.63kPa。

計算給出內部點A（24，24，22）和黏彈性邊界上點B（48，24，22）的結果，如圖2模型中標注：A點位置為上表面中心正下方垂直距離為2 m處，B點為垂直黏彈性邊界面對稱軸上的點，與A點處于同一高度。計算發(fā)現(xiàn)土體在很大范圍主要受自由表面和荷載的影響，因此內部取最具代表性的點A。B為黏彈性邊界對稱軸上的點，豎向位移最大且具有不為零孔壓。

5.1 沖擊荷載計算模型的面荷載為沖擊荷載，荷載形式為半周期正弦波（周期為0.4s），峰值為1000 Pa即P=1000sinωt，如圖3所示。圖4和圖5分別A點的孔壓和豎向位移時程曲線。由圖4可以看出不同的邊界形式對A點的孔壓值影響不大。分析原因可能是因為A點的孔壓主要受表面荷載和自由表面的影響，邊界對其影響很小，而豎向位移圖5中可以看出，黏彈性邊界的結果和參考解一致，比固定邊界結果更合理。圖6和圖7為B點的孔壓和豎向位移時程曲線?？梢悦黠@地看出，黏彈性邊界的結果和參考結果基本一致，而固定邊界在后面則出現(xiàn)出了嚴重的偏離。

圖4 沖擊荷載下A點孔壓時程曲線

圖5 沖擊荷載下A點豎向位移時程曲線

圖6 沖擊荷載下B點孔壓時程曲線

圖7 沖擊荷載下B點豎向位移時程曲線

5.2 突加荷載計算模型的面荷載為突加荷載，在0.1 s前為四分之一周期正弦波（周期為0.4 s），峰值為1 000 Pa，0.1 s后維持峰值不變，如圖8所示。圖9和圖10分別為A點的孔壓和位移時程曲線，圖11和圖12分別為B點的孔壓和豎向位移時程曲線。從這些計算結果可以看出，突加荷載表現(xiàn)出了和沖擊荷載結果趨勢相同，黏彈性邊界結果和參考結果基本一致，而固定邊界結果則發(fā)生了一定程度的偏離，其中B點的響應偏離更為嚴重。

圖8 突加荷載

本文給出兩個數(shù)值算例結果雖然不同，但卻表現(xiàn)出了相同的趨勢，可以看出不同邊界對A點的孔壓值影響不大，是因為A點的孔壓主要受表面荷載和自由表面的影響，邊界對其影響很小。我們還發(fā)現(xiàn)，沖擊荷載和突加荷載的最大值為1 000 Pa，而A點孔壓的最大值都達到了1 200 Pa，出現(xiàn)了曼德爾效應，即上表面自由排水的飽和地基在上表面受到荷載作用時會發(fā)生固結，土中的水排出，孔壓消散，但是孔壓不會立即消散，而是先上升后才慢慢消散。原因是上表面在荷載作用下自由排水，土體體積收縮，體積模量增大，但是內部的水還來不及流出，上部土體對內部土體產生了束縛作用，使得孔壓值有一個上升，隨后才慢慢消散。而B點的孔壓黏彈性邊界與參考邊界的結果不完全一致，主要原因可能是在計算中孔壓和位移的數(shù)量級差別比較大，位移發(fā)生很小的變化，對孔壓的影響都很大，而且由于在推導流量公式的過程中假設k=0，忽略了第二類壓縮波的影響，而計算過程中k=0.01，采用黏彈性邊界的計算結果與參考解相比具有一定誤差，在文獻［12］給出二維模型計算中也得到類似的結果，但這些誤差在工程上是可以接受的，因此該結果適用于中砂和更小顆粒土層。

圖9 突加荷載下A點孔壓時程曲線

圖10 突加荷載下A點豎向位移時程曲線

圖11 突加荷載下B點孔壓時程曲線

圖12 突加荷載下B點豎向位移時程曲線

6 結論

本文基于球面波動方程導出了無限大三維飽和地基人工黏彈性邊界及其流量邊界條件。在此基礎上給出了具有人工黏彈性邊界的飽和土體動力固結方程一般弱解形式。然后基于FEPG有限元自動生成系統(tǒng)編寫了腳本文件，生成了計算源代碼程序。通過數(shù)值算例驗證了黏彈性邊界的可靠性，由計算結果和參考結果、固定邊界結果的對比可以看出，黏彈性邊界作為人工邊界可以很好地模擬原始無限邊界的阻尼效應。進一步分析發(fā)現(xiàn)，當滲透系數(shù)比較小時，結果更接近真實解，但是對于包括地震工程的大部分中低振動頻率的工程問題，該邊界可以很好地模擬中砂和更小顆粒土層的邊界，也驗證了所開發(fā)程序的正確性。所編寫腳本文件，可以便捷的實現(xiàn)黏彈性邊界的應用和施加，從而很方便的用于飽和地基的動力計算。

參考文獻：

［1］廖振鵬.工程波動理論導論［M］.第2版.北京：科學出版社，2002：136-187.

［2］Lysmer J，Kuhlemeyer R L.Finite Dynamic Modelfor Infinite Media［J］.Journal EngineeringMethods，1969，95（EM4）：859-877.

［3］Deeks A J，Randolph M F.Axisymmetric Time-domain Transmitting Boundary［J］.Journal of Engineering Me?chanics，1994，120（1）：25-42.

［4］劉晶波，呂彥東.結構—地基動力相互作用問題分析的一種直接方法［J］.土木工程學報，1998，31（3）：55-64.

［5］劉晶波，王振宇，杜修力，等.波動問題中的三維時域粘彈性人工邊界［J］.工程力學，2005，22（6）：46-51.

［6］劉晶波，杜義欣，閆秋實.粘彈性人工邊界及地震動輸入在通用有限元軟件中的實現(xiàn)［J］.防災減災工程學報，2007，27：37-42.

［7］馬懷發(fā)，王立濤，陳厚群.粘彈性人工邊界的虛位移原理［J］.工程力學，2013，30（1）：168-174.

［8］郝明輝，張艷紅，陳厚群.基于ABAQUS的黏彈性人工邊界在重力壩分析中的應用［J］.中國水利水電科學研究院學報，2012，10（2）：120-126.

［9］Modaressi H，Benzenati I.An absorbing boundary element for dynamic analysis of two-phase media［C］//Earth?quake Engineering，Tenth Word Conference，Madrid，Spain.Rotterdam：Balkema，1992：1157-1161.

［10］Akiyoshi T，F(xiàn)uchida K，F(xiàn)ang H L.Absorbing boundary conditions for dynamicanalysis of fluid-saturated porous media［J］.Soil Dynamics and Earthquake Engineering，1993，12：387-397.

［11］王子輝，趙成剛，董亮.流體飽和多孔介質黏彈性動力人工邊界［J］.力學學報，2006，38（5）：605-611.

［12］劉光磊，宋二祥.飽和無限地基數(shù)值模擬的粘彈性傳輸邊界［J］.巖土工程學報，2006，28（12）：2128-2133.

［13］Zienkievicz O C，Shiomi T.Dynamic behavior of saturated porous media；The generalized Biot formulation and its numerical solution［J］.International Journal for Numerical and Analytical Methods in Geomechanics，1984，8：71-96.

［14］楊軍，宋二祥，陳肇元.飽和土動力反應中兩類壓縮波的獨立作用分析［J］.巖土力學，2001，22（2）：199-203．

［15］谷音，劉晶波，杜義欣.三維一致粘彈性人工邊界及等效粘彈性邊界單元［J］.工程力學，2007，24（12）：31-37.

［16］梁國平，周永發(fā).有限元語言［M］.北京：科學出版社，2013.

［17］梁國平，梁國平.有限元程序自動生成系統(tǒng)與有限元語言［J］.力學進展，1990，20（2）：199-204.

3D viscous-spring artificial boundaries of infinite saturated foundation

MA Huaifa1，2，YANG Ru2
（1.State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin，China Institute of Water Resources and Hydropower Research，Beijing100038，China；2.Earthquake Engineering Research Center，China Institute of Water Resources and Hydropower Research，Beijing100048，China）

Aiming at the low frequency vibration problems in engineering，regardless of fluid acceleration，the viscous-spring artificial boundary and flux boundary conditions infinite of 3D saturated ground are de?rivedbasedonthe spherical wave equation.Afterwards，the principle of virtual displacements of 3Dvis?cous-spring artificial boundary of infinite saturated foundation is proposed.Then the script files are written according to the finite element language，and based on the Finite Element ProgramGenerator（FEPG），the program source codes were generated.By means of the impact load and sudden load numerical analysis ex?acted on infinite ground surface，the correctness of the established equations and procedures were verified. This method can be applied to the numerical simulation of 3D saturated soil dynamic consolidation.

infinite saturated foundation；3Dviscous-spring artificial boundary；virtual displacement princi?ple；pore water pressure；effective stress

TU470

A

10.13244/j.cnki.jiwhr.2015.02.008

1672-3031（2015）02-0128-08

（責任編輯：韓昆）

2014-11-24

國家自然科學基金項目（51079164）；水利部公益性行業(yè)科研專項（201301057）；中國水利水電科學研究院科研專項（KJ1242）

馬懷發(fā)（1962-），男，山東棗莊人，博士，教授級高級工程師，從事計算力學、水工結構抗震及混凝土細觀力學分析研究。E-mail：mahf@iwhr.com

楊茹（1989-），女，陜西，碩士生，從事工程數(shù)值模擬。E-mail：yr1002151615@163.com