三維層狀飽和地基動(dòng)力Green 函數(shù)的薄層法求解

王 鵬，史 吏，王 軍，劉 凱

（1.溫州大學(xué) 建筑與土木工程學(xué)院，浙江 溫州 325035；2.浙江大學(xué) 軟弱土與環(huán)境土工教育部重點(diǎn)實(shí)驗(yàn)室，杭州 310027；3.香港城市大學(xué) 科學(xué)與工程學(xué)院，香港）

1 引 言

Green 函數(shù)在土與結(jié)構(gòu)動(dòng)力相互作用分析中有重要應(yīng)用，Wolf[1]解析推導(dǎo)了單相層狀地基平面內(nèi)、外運(yùn)動(dòng)的動(dòng)力剛度矩陣，并運(yùn)用到土-結(jié)構(gòu)頻域子結(jié)構(gòu)法動(dòng)力相互作用分析中。對(duì)于飽和地基，Halpern[2]采用柱坐標(biāo)勢(shì)函數(shù)分解和Hankel 積分變化求解了飽和半空間表面作用集中簡(jiǎn)諧荷載時(shí)，土骨架和孔隙水的豎向位移響應(yīng)，并將此基本解運(yùn)用于飽和半空間上覆剛性板動(dòng)力阻抗函數(shù)的求解；Philippacopoulos[3]采用類似方法推導(dǎo)了飽和半空間表面施加豎向集中簡(jiǎn)諧荷載至土骨架時(shí)，土骨架的徑向和豎向位移解，但僅給出了變換域中的位移譜圖；陳勝立等[4]采用一組改進(jìn)的飽和彈性土控制方程，通過Hankel 積分變換給出了飽和半空間埋置點(diǎn)荷載的軸對(duì)稱Lamb 問題的積分解答，繼而采用數(shù)值逆變換求解了地表豎向位移；王建華等[5]亦采用Hankel 積分變換方法求解了飽和半空間內(nèi)部作用簡(jiǎn)諧豎向力時(shí)的Green 函數(shù)，飽和土體采用Biot 動(dòng)力控制方程；王小崗[6]基于飽和土體Biot 動(dòng)力控制方程，采用Fourier 級(jí)數(shù)展開和Hankel 變換給出了柱坐標(biāo)下橫觀各向同性飽和半空間波動(dòng)方程的非軸對(duì)稱通解，在此基礎(chǔ)上建立了飽和半空間的精確動(dòng)力剛度矩陣，求解了柱面荷載作用下的層狀橫觀各向同性飽和地基的動(dòng)力Green 函數(shù)；丁伯陽等[7]采用飽和土體的Biot 波動(dòng)方程，推導(dǎo)了固相和液相Green 函數(shù)頻域解析解，并應(yīng)用于飽和土中基樁動(dòng)力特性分析中。

以上文獻(xiàn)均采用積分變化法求解土體動(dòng)力控制方程，最終所得空間-頻域解答為含Hankel 函數(shù)的無窮積分。文獻(xiàn)中多采用數(shù)值積分來完成上述無窮積分，從而不可避免地引入了截?cái)嗾`差，同時(shí)含Hankel 函數(shù)的被積函數(shù)在無窮積分區(qū)域上具有振蕩特性[8]，數(shù)值積分方法的選擇對(duì)積分結(jié)果穩(wěn)定性和精確性的影響很大。時(shí)剛等[9]應(yīng)用薄層法獲得了層狀飽和地基的Lamb 問題解答，并構(gòu)造了適用于薄層法底層邊界的通用旁軸近似解答。以上解答中，含Hankel 函數(shù)的逆變換是解析給出的，不再需要進(jìn)行數(shù)值逆變換，從而簡(jiǎn)化了計(jì)算過程，提高了計(jì)算效率。

本文建立了三維層狀飽和地基模型并在模型底邊施加Degrande 等[10]提出的頻域無反射邊界條件，采用Fourier 變換和正交變換將飽和地基Biot 動(dòng)力控制方程化為Love 和Rayleigh 模態(tài)方程，將地基沿深度方向劃分為數(shù)個(gè)薄層單元并對(duì)單元位移場(chǎng)采用線性插值函數(shù)，從而將模態(tài)方程離散化，通過求解本征值問題，最終給出了三維層狀飽和地基的頻域動(dòng)力Green 函數(shù)。在此基礎(chǔ)上，分析了均質(zhì)飽和地基的Love 和Rayleigh 模態(tài)的彌散特性，求解了均質(zhì)和成層飽和半空間的 Lamb 問題，并與Philippacopoulos 的數(shù)值積分解對(duì)比，說明了所求Green 函數(shù)的正確性。文中所得Green 函數(shù)無需關(guān)于Hankel 函數(shù)進(jìn)行數(shù)值逆變換且能直接計(jì)入地基的成層性，從而提高了計(jì)算效率。

2 層狀飽和地基Green 函數(shù)

直角坐標(biāo)系中，飽和土體Biot[11]動(dòng)力控制方程為

本構(gòu)方程為

式中：ui、wi(i=x,y,z)分別為土骨架和孔隙流體相對(duì)于土骨架的位移；λ、μ為土骨架Lame 常數(shù)；M、α為Biot 常數(shù)；ρ=nρf+(1-n)ρs，ρs、ρf分別為土骨架和孔隙流體的質(zhì)量密度；m=a∞ρf/n，a∞為孔隙彎曲系數(shù)，n為孔隙度；b=ρfg/kD，kD為土體Darcy 滲透系數(shù)；θ=ui,i、?=-wi,i分別為土骨架和孔隙流體的體積應(yīng)變；σij為土體總應(yīng)力；p為孔壓。采用滯回阻尼表征土骨架的材料阻尼，即土骨架Lame 常數(shù)替換為λ*=λ(1+2iβ)及μ*=μ(1+2iβ)，其中i為單位虛數(shù)，β為滯回阻尼比。定義飽和土體特征頻率為ωc=bρ/(mρ-ρf×ρf)[10]，對(duì)于不同的荷載激振頻率ω，定義無量綱頻率比為χ=ω/ωc。

定義Fourier 變換對(duì)和正交變換對(duì)分別為

式中：上標(biāo)“^”表征波數(shù)-頻率域變量；上標(biāo)“-”表征正交變換域中的變量；kx和ky分別為與空間坐標(biāo)x 和y 對(duì)應(yīng)的波數(shù)，ω為圓頻率；

利用Fourier 正變換將控制方程（1）～（4）變換到波數(shù)-頻率域中，再針對(duì)位移和應(yīng)力引入式（6）～（7）的正交變換，從而可將原本耦合的控制方程化為解耦的Love 和Rayleigh 模態(tài)控制方程，其對(duì)應(yīng)的應(yīng)力表達(dá)式亦為解耦的，見式（8）～（10）。

Love 模態(tài)：

Rayleigh 模態(tài)：

本構(gòu)關(guān)系：

建立如圖1 所示的層狀飽和地基模型，在模型深度L 范圍內(nèi)（沿z 軸）劃分為N個(gè)薄層，薄層厚度為hi（i=1～N）。各薄層土體材料性質(zhì)相同，不同薄層之間的土體材料性質(zhì)和厚度可以不同。模型底部施加Degrande 等[10]提出的無反射邊界條件，從而模擬土體下臥半空間。

圖1 三維層狀飽和地基薄層法分析模型Fig.1 Thin layer element model of 3D layered saturated ground

由式（8）～（10）可知，控制方程均為關(guān)于坐標(biāo)z 的常微分方程。對(duì)于特定薄層i，建立如圖1右下角所示的局部坐標(biāo)體系，并引入如下線性位移插值函數(shù)：

將式（8）～（9）代入式（11），運(yùn)用Galerkin法，可得薄層i 的Love 和Rayleigh 模態(tài)積分弱形式控制方程如下所示。

薄層i 的Love 模態(tài)為

薄層i 的Rayleigh 模態(tài)為

對(duì)于由N個(gè)薄層組成的飽和地基模型，組裝式（12）、（13）可得地表無應(yīng)力時(shí)層狀地基Love 模態(tài)和Rayleigh 模態(tài)的薄層法控制方程（由薄層i 的應(yīng)力矢量表達(dá)式可知，層間應(yīng)力在組裝時(shí)會(huì)相互抵消）如下所示。

飽和層狀地基Love 模態(tài)為

飽和層狀地基Rayleigh 模態(tài)為

式中：

分別為模型的Rayleigh 和Love模態(tài)總體未知節(jié)點(diǎn)位移矢量；分別為模型Rayleigh 和Love 模態(tài)的總體節(jié)點(diǎn)應(yīng)力矢量。

在地基模型底部施加Degrande 等[10]提出的無反射邊界條件，從而模型底部的節(jié)點(diǎn)應(yīng)力可表達(dá)為該節(jié)點(diǎn)未知節(jié)點(diǎn)位移的函數(shù)，如下式所示：

其中，矩陣元素bij(i,j=1,2,3)的表達(dá)式見文獻(xiàn)[10]。若模型底邊為固定邊界，則有位移邊界條件，此時(shí)與這些位移自由度相對(duì)應(yīng)的應(yīng)力將不出現(xiàn)在式（14）、（15）的右端項(xiàng)中。

將式（16）帶入式（14）、（15）并移至方程左端，從而可得飽和層狀地基Rayleigh 模態(tài)和Love模態(tài)的齊次控制方程，即形成了模型的本征值問題。求解Love 模態(tài)本征值問題可得2(N+1)個(gè)共軛特征值對(duì)kLm和對(duì)應(yīng)的特征向量Φym(m=1～2(N+1))；求解Rayleigh 模態(tài)本征值問題可得4(N+1)個(gè)共軛特征值對(duì)kRm和對(duì)應(yīng)的特征向量Φxm、Φz(mì)m(m=1,…,4(N+1))；但由于矩陣Gs和As中與節(jié)點(diǎn)自由度和對(duì)應(yīng)的元素全部為0，從而使得Love 模態(tài)共有(N+1)個(gè)非零共軛特征值對(duì)，而Rayleigh 模態(tài)共有3(N+1)個(gè)非零共軛值對(duì)，且僅虛部為負(fù)的特征值方才滿足無窮處的輻射條件。最終正交變換域中的位移可采用振型疊加法表達(dá)為

其中，

為第m 階特征向量；qRm、qR′m和 qLm為對(duì)應(yīng)的第m階振型參與系數(shù)。

同理，作用

時(shí)，由式（15）和式（17）第1 式可求得以下兩式：

將式（18）～（20）合寫即可得正交變換域內(nèi)的飽和層狀地基Green 函數(shù)：

應(yīng)用式（6）～（7）中的正交逆變換，將式（21）化為波數(shù)-頻率域的Green 函數(shù)可有

式（22）中系數(shù)矩陣各元素表達(dá)式為

應(yīng)用式（5）的Fourier 逆變換將式（22）變換至空間-頻率域，注意到式（22）的系數(shù)矩陣中僅為波數(shù)kx和ky的函數(shù)，記其逆變換為(i=1～9)，可有

其中，函數(shù)F1(x)，F(xiàn)2(x)和F3(x)表達(dá)式為

薄層模型任意節(jié)點(diǎn)（激振點(diǎn)）施加集中點(diǎn)荷載至土骨架或孔隙流體時(shí)，各節(jié)點(diǎn)（觀察點(diǎn)）的土骨架和孔隙流體位移響應(yīng)可由式（25）算得。故式（25）為三維層狀飽和地基空間-頻率域的Green 函數(shù)。但由(i=1～9)的表達(dá)式可知，當(dāng)激振點(diǎn)和觀察點(diǎn)的平面距離r=0時(shí)會(huì)出現(xiàn)奇異性，此時(shí)可用等效均布圓環(huán)或圓盤荷載（分布半徑為r0）代替點(diǎn)荷載來求解觀察點(diǎn)位移響應(yīng)。下文分析和結(jié)果呈現(xiàn)均在空間-頻率域中進(jìn)行，從而均省略變量上標(biāo)“～”。

以上求解過程在Intel Visual Fortran?平臺(tái)下采用Fortran 95 語言編寫，本征值方程采用IMSL?數(shù)學(xué)包提供的子程序求解。

3 算例分析

3.1 均質(zhì)飽和地基Love 及Rayleigh 模態(tài)

為研究不同滲透系數(shù)和激振頻率（ω）下的飽和地基Love 及Rayleigh 模態(tài)的彌散特性，此處僅取均質(zhì)單層飽和地基模型，模型底邊固定。土體材料參數(shù)取文獻(xiàn)[10]中的Molsand 參數(shù)，但滲透系數(shù)此處為可變值，模型參數(shù)如表1 所示。

表1 飽和土體材料及薄層參數(shù)Table1 Parameters of the saturated ground and the thin layer

圖2 kD=1.0×10-4 m/s時(shí)Love 波數(shù)和Rayleigh 波數(shù)隨激振頻率的變化Fig.2 Variations of Love and Rayleigh wave numbers with respect to excitation frequency when kD=1.0×10-4 m/s

圖3 kD=+∞(b=0)時(shí)Love 波數(shù)和Rayleigh 波數(shù)隨激振頻率的變化Fig.3 Variations of Love and Rayleigh wave numbers with respect to excitation frequency when kD=+∞(b=0)

由于模型底邊固定，此時(shí)Love 模態(tài)只有1個(gè)非零特征值（波數(shù)）kL1，Rayleigh 模態(tài)則有3個(gè)非零特征值（波數(shù)）kR1、kR2和kR3。圖2、3 分別呈現(xiàn)了滲透系數(shù)kD為1.0×10-4m/s 和+∞（b=0）時(shí)Love 波數(shù)及Rayleigh 波數(shù)隨激振頻率ω 的變化，波數(shù)實(shí)部（Re）代表了波的傳播特性，虛部（Im）則代表了衰減特性。其中，橫軸激振頻率ω 無量綱化為ωL/cs0，，為飽和地基完全不透水時(shí)的剪切波速；縱軸波數(shù)k(k=kL1,kR1,kR2,kR3)無量綱化為kL。

由圖2(a)可知，Love 波在激振頻率低于S 波的截止頻率時(shí)，kL1實(shí)部為0，即此時(shí)Love 波并不傳播。而當(dāng)激振頻率高于截止頻率時(shí)，其實(shí)部隨ω 線性增加，但虛部接近于0，即此時(shí)Love 波是傳播的；由圖2(b)可知，kR1的實(shí)部和虛部均隨激振頻率明顯增大，即該Rayleigh 波是傳播的，但很快衰減。且該波在激振頻率很小時(shí)就已傳播，無截止頻率；由圖2(c)可知，此Rayleigh 波僅在ω 高于P1 波的截止頻率時(shí)方才傳播，且衰減很??；由圖2(d)可知，該Rayleigh 波存在兩個(gè)截止頻率，分別對(duì)應(yīng)S 和P1的截止頻率。當(dāng)ω 介于二者之間或大于P1 截止頻率時(shí)，該波是傳播的，且衰減很小。當(dāng)前滲透系數(shù)下，對(duì)于表1 中的激振頻率變化范圍，無量綱頻率比χ ＜ 0.01，此時(shí)P2 波傳播速度很小可忽略不計(jì)，且其衰減系數(shù)遠(yuǎn)較P1 波和S 波為大，故而圖2 Rayleigh 波傳播特性中并未體現(xiàn)P2 波截止頻率的影響。

對(duì)比圖2、3 可知，滲透系數(shù)變化對(duì)kL1和kR3影響很小，但對(duì)kR1和kR2影響很大。當(dāng)前滲透系數(shù)下，使得飽和地基的特征頻率ωc為0，即χ為無窮大，此時(shí)P2 波速大于S 波速，故而kR2中體現(xiàn)了P2 波截止頻率的影響。

3.2 均質(zhì)飽和半空間Lamb 問題

為驗(yàn)證所得Green 函數(shù)的正確性，運(yùn)用式（24）求解了均質(zhì)飽和半空間的Lamb 問題并與Philippacopoulos[3]的數(shù)值積分解進(jìn)行對(duì)比。被積函數(shù)包括了Hankel 函數(shù)，為振蕩積分，此處采用Hidenori[8]提供的方法完成數(shù)值積分。飽和土體材料參數(shù)同表1，滲透系數(shù)kD=1.0×10-4m/s。

圖4為無量綱頻率比χ=0.01 和0.10時(shí)，薄層法和數(shù)值積分法所得半空間表面土骨架豎向位移（uz）的對(duì)比圖。圖中橫坐標(biāo)r為激振點(diǎn)與觀察點(diǎn)之間的平面距離，無量綱化為r/λs0，其中λs0=2π cs0/ω；縱坐標(biāo)uz無量綱化為uz/uz0，其中uz0為均質(zhì)單相彈性半空間的Boussinesq 解。單相彈性半空間的材料參數(shù)μ 和ν 同表1。薄層法求解時(shí)，模型深度L=10 λs0，各薄層厚度h=λs0/6，薄層總數(shù)N=60，。經(jīng)試算，數(shù)值積分法在包括Bessel 函數(shù)前2 000個(gè)0 點(diǎn)時(shí)，數(shù)值積分結(jié)果趨于穩(wěn)定，此處取前3000個(gè)0 點(diǎn)。

圖4 均質(zhì)飽和半空間Lamb 問題薄層法求解與數(shù)值積分法求解對(duì)比Fig.4 Comparison between the thin layer element method and numerical integration solutions of Lamb problem of saturated homogeneous ground

由圖4 可知，薄層法與數(shù)值積分法的結(jié)果變化趨勢(shì)完全相同，且兩種方法給出的豎向位移值十分接近，由此說明了所得Green 函數(shù)的正確性。注意r/λs0＜0.5時(shí)，Hidenori 方法引入的截?cái)嗾`差導(dǎo)致了數(shù)值積分解的波動(dòng)。

3.3 成層飽和半空間Lamb 問題

為研究地基成層性的影響，此處考慮上軟、下硬雙層飽和半空間的Lamb 問題。第1 層飽和土體參數(shù)同表1、層厚記為L1、最短波長記為λs1=2π cs1/ω，其中，μ1為土層1 的剪切模量,各薄層厚度h1=λs1/10；第2 層飽和地基土體參數(shù)除剪切模量μ2=2μ1外，其余參數(shù)同表1，厚度為L2、最短波長記為λs2=2π cs2/ω，其中，各薄層厚度h2=λs2/10。兩層地基的滲透系數(shù)kD均為1.0×10-4m/s。計(jì)算時(shí)，L2=4λs2，L1則分別取為0.2λs1、0.5λs1、λs1和4λs1，無量綱頻率比χ=0.01。半空間表面土骨架豎向位移（uz）計(jì)算結(jié)果如圖5所示，其中橫坐標(biāo)為無量綱激振頻率ωr/cs1，縱坐標(biāo)無量綱化為uz/uz0，uz0為均質(zhì)單相彈性半空間的Boussinesq 解，彈性半空間的材料參數(shù)μ 和ν 同土層1。

圖5 上軟下硬雙層飽和地基Lamb 問題Fig.5 Lamb problem of two-layer saturated ground(soft up layer and stiff second layer)

由圖5 可知，土體分層對(duì)土骨架豎向位移的影響十分明顯：L1/λs1=0.2時(shí)位移幅值較低；L1/λs1=0.5時(shí)位移幅值增大；L1/λs1增大至1 和4時(shí)則位移幅值減小，但位移變化趨勢(shì)趨于一致。

4 結(jié) 論

（1）滲透系數(shù)變化對(duì)Love 模態(tài)影響很小，對(duì)Rayleigh 模態(tài)有明顯影響。

（2）Love 波和Rayleigh 波在飽和地基中傳播時(shí)存在截止頻率。

（3）本文所得Green 函數(shù)是正確的，且無需關(guān)于Hankel 函數(shù)進(jìn)行數(shù)值逆變換，提高了計(jì)算效率。

（4）土體分層對(duì)地表位移響應(yīng)有明顯影響，本文所得Green 函數(shù)能直接考慮地基土體的成層性。

[1]WOLF J P.Soil-structure-interaction analysis in time domain[M].[S.l.]:Prentice Hall,1988.

[2]HALPERN M R.Response of poroelastic halfspace to steady-state harmonic surface traction[J].International Journal for Numerical and Analytical Methods in Geomechanics,1986,10:609－632.

[3]PHILIPPACOPOULOS A J.Lamb’s problem for fluid-saturated,porous media[J].Bulletin of the Seismological Society of America,1988,78(2):908－923.

[4]陳勝立,張建民,陳龍珠.飽和地基中埋置點(diǎn)源荷載的動(dòng)力Green 函數(shù)[J].巖土工程學(xué)報(bào),2001,23(4):423－426.CHEN Sheng-li,ZHANG Jian-min,CHEN Long-zhu.Dynamic Green’s functions of saturated soils subjected to the internal excitation[J].Chinese Journal of Geotechnical Engineering,2001,23(4):423－426.

[5]王建華,陸建飛,沈?yàn)槠?半空間飽和土在內(nèi)部簡(jiǎn)諧垂直力作用下的Green 函數(shù)[J].水利學(xué)報(bào),2001,3:54－57.WANG Jian-hua,LU Jian-fei,SHEN Wei-ping.The Green function of harmonic vertical load applied on the interior of the half space saturated soil[J].Journal of Hydraulic Engineering,2001,3:54－57.

[6]王小崗.橫觀各向同性飽和土中柱面荷載的動(dòng)力Green 函數(shù)[J].力學(xué)學(xué)報(bào),2010,42(5):909－918.WANG Xiao-gang.Dynamic Green’s function for internal barrel loads in transversely isotropic saturated soils[J].Chinese Journal of Theoretical and Applied Mechanics,2010,42(5):909－918.

[7]丁伯陽,童曉剛,陳軍,等.排水狀態(tài)下飽和土中基樁動(dòng)力特性分析[J].巖土工程學(xué)報(bào),2013,35(1):85－95.DING Bo-yang,TONG Xiao-gang,CHEN Jun,et al.Dynamic properties of a single pile in saturated soil under drainage[J].Chinese Journal of Geotechnical Engineering,2013,35(1):85－95.

[8]HIDENORI OGATA.A numerical integration formula based on the Bessel functions[J].Research Institure for Mathematical Sciences, 2005,41:949－970.

[9]時(shí)剛,高廣運(yùn),馮世進(jìn).飽和層狀地基的薄層發(fā)基本解及其旁軸邊界[J].巖土工程學(xué)報(bào),2010,32(5):664－671.SHI Gang,GAO Guang-yun,FENG Shi-jin.Basic solution of saturated layered ground by thin layered method and its paraaxial boundary[J].Chinese Journal of Geotechnical Engineering,2010,32(5):664－671.

[10]DEGRANDE G,ROECK G D.An absorbing boundary condition for wave propagation in saturated poroelastic media——Part I:Formulation and efficiency evaluation[J].Soil Dynamics and Earthquake Engineering,1993,12(7):411－421.

[11]BIOT M A.Theory of propagation of elastic waves in afluid saturated porous solid.I:Low frequency range[J].Journal of Acoustic Society of American,1956,28(2):168－178.

[12]蔣通.地基-結(jié)構(gòu)動(dòng)力相互作用分析方法——薄層法原理及應(yīng)用[M].上海:同濟(jì)大學(xué)出版社,2009.