论文部分内容阅读
摘要： 为了研究具有复杂表面的饱和土体与移动荷载相互作用机理，根据饱和土Biot理论，采用Fourier变换和势函数分解法，推导了饱和土体频域波数内的Green函数；建立了移动荷载作用下，饱和土体中孔洞动力响应的频域波数域内边界元法方程（2.5D boundary element methods）；利用快速Fourier逆变换法，得到了时间空间域内饱和土体的动力响应.研究结果表明：在频域内，利用移动荷载方向的一致性建立的频域波数域内边界积分方程，可将3D空间问题转化为频域波数域的2D平面问题， 3D边界元简化为2D边界元，使得计算面转换为计算线，减小了计算规模.
关键词： 格林函数；移动荷载；饱和土；Fourier积分变换
中图分类号： TU 471文献标志码： AAnalysis of Dynamic Interaction between Hole Embedded in
Saturated Soil and Moving Loads
Using 2.5D Boundary Element MethodXU Bin1，2，GAO Liang2，LEI Xiaoyan3，XU Manqing1，LIU Linya3
（1. Department of Civil Engineering， Nanchang Institute of Technology， Nanchang 330029， China； 2. School of Civil Engineering and Architecture， Beijing Jiaotong University， Beijing 100044， China； 3. Engineering Research Center of Railway Environmental Vibration and Noise， Ministry of Education， East China Jiaotong University， Nanchang 330013， China）
Abstract：To study the mechanism of interaction between the saturated soil with complex surface and moving loads， a specific 2.5D boundary element method （BEM） for the problem of dynamic responses of the hole embedded in the saturated soil under a moving load was derived systematically. Based on Biots theory， the frequency wave number domain Greens function for saturated porous media was developed using the Fourier transform method and the potential decomposition approach. Then， dynamic responses in the timespace domain solution were further obtained by synthesizing the wave number solution via inverse fast Fourier transform （IFFT）. The results show that using the frequency wavenumber boundary integrated equations established by the consistency of the directions of moving load， a 3D spatial model for the dynamic interaction between the holes embedded in the saturated soil with moving loads is convertible to a 2D planar model in the frequency domain， and the 3D complex boundary can be simplified to a problem of 2D boundary element. Thus， the complicated boundary planes can be transformed into computed lines using the 2.5D wave number domain BEM formulation， consequently reducing the amount of computation to a manageable size.
Key words：Greens function； moving loads； saturated soil； Fourier integral transform
工程中穿越软土地基或道路下的市政供水、天然气等管道不可避免受路面交通荷载作用，当埋地管道一旦达到破坏强度就可能产生安全事故，造成水、天然气等泄漏而引发燃烧与爆炸[12].上述工程问题均可归结于移动荷载与软土地基中孔洞结构相互作用的动力响应模型.
对于具有简单几何构形的半空间土体表面受移动载荷问题，可利用解析方法进行计算[35]，但对于移动荷载作用下，地下孔洞结构的动力响应，边界较复杂，其几何轴对称性已破坏，难以采用解析法计算.特别是移动荷载作用下，土体中各点动力响应都随时间和空间不同而改变，采用有限元法对土体内部单元网格划分时，由于时间和空间的变化，不可避免造成单元数量急剧增加，计算困难.边界元方法（BEM）具有自动满足无穷远处辐射条件及降维功能，已被证明是解决具有复杂边界条件问题一种有效手段，并且在工程中广泛应用.文献[6]中以3个土骨架位移和孔隙水压力为基本变量，建立了饱和土体边界积分方程.采用相同西南交通大学学报第48卷第4期徐斌等：移动荷载与土体中孔洞相互作用的2.5D边界元法分析的基本变量，文献[7]中推导了3D饱和土体频域内边界积分方程.文献[8]中利用3D频域边界积分方程，分析了饱和多孔介质中散射问题.文献[9]中利用拉普拉斯域的Green函数和卷积积分法，推导了时域内边界积分方程方法.采用3D边界元法，文献[10]中对饱和土体中群桩土共同作用的动力问题进行了分析. 与传统的有限元法相比，边界元法只对边界进行单元离散，因此，建模单元、计算时间均可节省.然而， 3D时间域内边界元法仍然需占用计算机大量内存和计算时间.文献[1112]在3D边界元的基础上，将结构作为2D分析或采用空间波数变换用来简化问题，提出了波数域的2.5D边界元法.与3D边界元法相比， 2.5D边界元法的边界离散和计算要求大大减少.目前2.5D边界元法已成功应用于弹性介质中的动力响应问题分析.文献[13]根据具有恒定移动速度的点荷载Green函数，推导了任意角入射波下的无限长河谷3D响应的2.5D间接边界积分方程.利用全空间Green函数，文献[14]采用直接2.5D BEM，分析了粘弹性半空间中无限长隧道衬砌的3D地震荷载作用下的动力响应.采用离散波数的有限元与边界元耦合，并结合傅里叶逆变换法.文献[15]对轨道交通荷载引起地面3D空间的动力响应进行了评价.
上述文献的研究表明：土体为单相弹性介质，忽略土体在移动载荷作用下土体中的水相、土颗粒之间的相互耦合作用，则无法考虑土体中孔压的变化，以及孔压变化对周围土体动力响应的影响[3，16].多孔软土中孔隙水压的变化对土体强度有较大的影响，甚至引起震陷、液化等工程现象，这些都使得软土地基具有非线性、非弹性的动力本质特征[1718].传统的线弹性、粘弹塑性土体本构模型难以有效反映土体变形，特别是孔隙水压力变化而导致的孔洞动态受力性状.目前关于饱和土体中孔洞与移动荷载相互作用的2.5D BEM法尚未见报道.本文分别推导了由单位集中力作用在土骨架时所对应Green函数，以及单位散度源作用在孔隙介质流体时所对应Green函数.考虑荷载移动方向与饱和土体中孔洞轴线方向具有一致性，在三维边界元法（3DBEM）的基础上，采用沿饱和土体中孔洞纵向进行空间波数的积分变换，形成2.5D BEM，将三维移动荷载问题转化为平面内边界元积分问题，既可避免移动荷载引起的半无限空间有限元单元网格划分随时间而变化的难题，又可避免3D BEM在移动荷载方向的边界划分问题，能够减少计算工作量.1饱和二相介质土Green函数为推导饱和地基土体模型的Green函数，定义时间t到频域ω，空间z到波数域kz的Fourier变换对为
f-（ω）=∫+∞-∞f（t）e-iωtdt，
f（t）=12π∫+∞-∞f-（ω）eiωtdω，
f^（kz）=∫+∞-∞f（z）eikzzdz，
f（z）=12π∫+∞-∞f^（kz）e-ikzzdkz，（1）
式中：
“-”、“^”分别表示频域、波数域量.
由于饱和多孔介质由土骨架及水相组成， Green函数张量G可分别由土骨架在单位集中力作用下所对应Green函数Gs、Gs以及饱和孔隙介质流体在单位散度源作用时所对应Green函数Gf、Gf 组成，表达式为
G=GsGf
GsGf.（2）
根据Biot理论[1718]，对于饱和二相介质，土骨架与水相分别在单位体积力Fi和fi（i=1，2，3）作用下，土体骨架位移ui（i=1，2，3）及孔隙水压（pf）4个独立变量表示的频域内饱和土体运动方程为
μi，jj+（λ+μ）j，ji+ρgω2i-αgf，i=
-（i-β4f-i），
f，jj+β2ω2f-β3j，j=f-i，i，（3）
式中： m=α∞ρf/；
αg=α-β4；
ρg=ρb-β4ρf；
ρb=（1-）ρs+ρf；
β1=M/[mω2-iω（η/k）]；
β2=1/（β1ω2）；
β3=ρfω2-α[mω2-iω（η/k）]；
β4=ρfω2β1/M，
其中， λ和μ为Lame常数， α、M分别为饱和土体中土骨架、水相的压缩系数， 为饱和多孔介质的孔隙率， wi（i=1，2，3）为饱和土体中流体相对于土骨架的位移， ρs、ρf分别为饱和土体土骨架、孔隙水密度， α∞为孔隙介质弯曲系数， k、η分别为饱和土体水相的渗透系数及孔隙流体的黏质度， ω为移动荷载频率.
饱和多孔介质中土骨架在集中力作用下，对应Green函数的Gs、Gs，可由式（3）求得，
μΔ 2Gs+（λ+μ）Δ （Δ Gs）+
ρgω2Gs-αgGs=-δ（R）I3×3，
Δ 2Gs+β2ω2Gs-β3Δ Gs=0，（4）
式中：
Δ为Laplace算子.
由文献[19]可知，可采用3个标量势函数f，Gs、s，Gs、f，Gs表示Gs、Gs，则有
Gs=Δ Δ f，Gs+Δ Δ s，Gs+
（Δ 2f，GsI3×3-Δ Δ f，Gs），
Gs=AfΔ Δ 2f，Gs+AsΔ Δ 2s，Gs.（5）
利用恒等式
-δ（R）I3×3=Δ Δ 14πR+
Δ 2I3×3-Δ Δ 14πR，（6）
式中：
R为两点的空间距离.
将式（5）代入式（4）中，同时利用式（6），可得
[afΔ 2f，Gs+ρgω2Δ 2f，Gs]+
[asΔ 2s，Gs+ρgω2Δ 2s，Gs]=14πR，
[AfΔ 2f，s+dff，Gs]+
[AsΔ 2s，Gs+dss，Gs]=0，
Δ 2Gs+k2tGs=14πμR，（7）
式中：
af=λ+2μ-αgAf；
as=λ+2μ-αgAs； df=β2Afω2-β3；
ds=β2Asω2-β3；
k2f=ρgω2/（λ+2μ-αgAf）；
k2s=ρgω2/（λ+2μ-αgAs）；
k2t=ρgω2/μ；
Re（kf）≤Re（ks）.
利用式（1），对式（7）进行z→kz的Fourier变换，整理后可得频域、波数域内的势函数为
f，Gs=iAs4k2f（afAs-asAf）[H0，（2）（γfr）-H0，（2）（-ikzr）]，
s，Gs=iAf4k2s（asAf-afAs）[H0，（2）（γsr）-H0，（2）（-ikzr）]，
Gs=i4ρgω2[H0，（2）（γtr）-H0，（2）（-ikzr）]，（8）
式中： “～”表示频域波数内的量；
H0，（2）为第二类的0阶Hankel函数；
r为两点的平面距离.
根据式（1）进行z→kz的Fourier变换，由式（5）可得频域、波数域内的Green函数Gs、Gs，再将式（8）代入后，可得在集中力作用下饱和多孔介质中土骨架所对应的频域波数域内Green函数为
ij，Gs=[δf，（s）H0，（2）（γfr）+δs，（s）H0，（2）（γsr）-δt，（s）H0，（2）（γtr）]，ij-
[δt，（s）H0，（2）（γtr）]δij，
i，Gs=χp[-H0，（2）（γfr）+H0，（2）（γsr）]，ii，j=1，2，3，（9）
式中：
δf，（s）=iAs/[4k2f（afAs-asAf）]；
δs，（s）=iAf/[4k2s（asAf-afAs）]；
δt，（s）=i/（4ρgω2）；
χp=-Ask2fδs，（s）.
利用式（3），可得单位散度源作用在饱和土体中流体介质时所对应Green函数Gf、Gf为
μΔ 2Gf+（λ+μ）Δ （Δ Gf）+
ρgω2Gf-αgGf=0，
Δ 2Gf+β2ω2Gf-β3Δ Gf=-δ（R）.（10）
由于流体介质中不存在剪切波，因此可用2个标量函数f，Gf、s，Gf表示相应的Green函数Gf、Gf，
Gf=Δ f，Gf+Δ s，Gf，
Gf=AfΔ 2f，Gf+AsΔ 2s，Gf.（11）
同单位集中力作用在土骨架时所对应的频域波数域内Green函数类似的步骤，可得单位散度源作用在饱和土体流体介质时，所对应的频域波数域内Green函数为
i，Gf=δ（f）[H0，（2）（γfr）-H0，（2）（γsr）]，i，
Gf=δ（f）[-Afk2fH0，（2）（γfr）+
Ask2sH0，（2）（γsr）]，i=1，2，3，（12）
式中：
δf=ias/[4k2f（asAf-afAs）].2饱和土体中2.5D边界积分方程2.1饱和土体频域内的边界积分方程根据互易定理及饱和土体的本构关系（式（3）），考虑饱和土体任意2个状态（状态1和状态2），可得
ij，（1）ij，（2）+f，（1）ξ（2）=
ij，（2）ij，（1）+f，（2）ξ（1），（13）
式中：下标（1）、（2）分别表示饱和土体的2种运动状态；
σij为饱和土体总应力分量；
ξ为饱和多孔介质中流体的单位体积增量，
ξ=wi，i.
将式（13）沿饱和土体任意区域Ω积分，并结合式（3），可得
∫S[j，（1）j，（2）-j，（2）j，（1）]dS-
∫S[f，（1）n，（2）-f，（2）n，（1）]dS=
∫Ω[f，（1）f，（2）-f，（2）f，（1）]dΩ-
∫Ω[gj，（1）j，（2）-gj，（2）j，（1）]dΩ，（14）
式中：
j为作用在区域Ω边界S上的面力，
j=ijni；
ni为边界S外法线方向余弦；
f=β1i，i/M；
n为沿边界S的外法线方向流体位移，
n=β1（f，j-ρfω2j）nj/M.
将饱和土体的真实解作为状态1，单位集中力或单位散度源所应于Green函数解为状态2，则在不考虑体积力时，由式（14）可得在饱和土体土骨架任意点x上作用一集中力时，饱和土体的Somigliana等式为
i（x，ω）=∫S[ij，Gs（x，ζ，ω）j（ζ，ω）-
ij，Gs（x，ζ，ω）j（ζ，ω）]dS（ζ）-
∫S[ni，Gs（x，ζ，ω）f（ζ，ω）-
i，Gs（x，ζ，ω）n（ζ，ω）]dS（ζ），（15）
式中：
ni，Gs=β1[i，Gs（x，ζ，ω）/ζj-
ρfω2ij，Gs]nj（ζ）/M，
ij，Gs=ijk，Gs（x，ζ，ω）nk（ζ），
其中， Gs、Gs分别为单位集中力作用在土骨架时所对应Green函数.
利用相同方法，饱和土体中流体介质受单位散度源作用时的Somigliana等式为
i（x，ω）=∫S[j，Gf′（x，ζ，ω）j（ζ，ω）-
j，Gf′（x，ζ，ω）j（ζ，ω）]dS（ζ）-
∫S[n，Gf′（x，ζ，ω）f（ζ，ω）- Gf′（x，ζ，ω）n（ζ，ω）]dS（ζ），（16）
式中：
j，Gf′=β1j，Gf/M；
j，Gf′=β1j，Gf/M；
n，Gf′=β1n，Gf/M；
Gf′=β1Gf/M.
根据式（2）及式（15）、（16）可得到饱和土体统一的Somigliana等式为
i（x，ω）=∫S[ij，G（x，ζ，ω）j（ζ，ω）-
ij，G（x，ζ，ω）j（ζ，ω）]dS（ζ）.（17）
利用文献[20]的方法，考虑不同边界的特性，可得频域内饱和土体的边界积分方程为
ciji（x，ω）=∫S[ij，Gs（x，ζ，ω）j（ζ，ω）-
ij，Gs（x，ζ，ω）j（ζ，ω）]dS（ζ），（18）
在光滑边界时， cij=δij/2.2.2 频域波数域的2.5D边界积分方程移动荷载与饱和土体中的孔洞相互作用模型如图1所示.
图1移动点荷载与饱和土体相互作用示意图
Fig.1Schematic diagram of the interaction between
the moving point loads and the saturated soil若荷载沿孔洞轴向以速度v移动，饱和土体及其中孔洞的动力响应均可表示为
ui（x，t）=ui（x⊥，z-vt），
ti（x，t）=ui（x⊥，z-vt）.（19）
对式（19）进行时间、空间Fourier变换，可得
i（x⊥，kz，ω）=
2πδ（ω+ω0-kzv）U^i（x⊥，kz），
i（x⊥，kz，ω）=
2πδ（ω+ω0-kzv）T^i（x⊥，kz），（20）
式中：
ω0为移动荷载初始频率；
U^i（x⊥，kz）=∫∞-∞ui（x⊥，z-vt）eikz（z-vt）d（z-vt），
T^i（x⊥，kz）=∫∞-∞ti（x⊥，z-vt）eikz（z-vt）d（z-vt），
其中， x⊥为平面内坐标点.
将式（20）代入式（18）中，可得到频率波数域的饱和土体中孔洞的2.5D边界元积分方程为
cijU^j（x⊥，kz）=
∫Γij，G（ζ⊥-x⊥，-kz，kzv）T^j（x⊥，kz）dΓ（ζ⊥）-
∫ΓT～ij，G（ζ⊥-x⊥，-kz，kzv）U^j（x⊥，kz）dΓ（ζ⊥），（21）
式中：
Γ 为饱和土体及孔洞在 平面内的边界；
ζ⊥为 平面内任意观察点坐标.
对于饱和土体-孔洞体系内任一点，有
U^j（x⊥，kz）=
∫Γij，G（ζ⊥-x⊥，-kz，kzv）T^j（x⊥，kz）dΓ（ζ⊥）-
∫ΓT～ij，G（ζ⊥-x⊥，-kz，kzv）U^j（x⊥，kz）dΓ（ζ⊥）.（22）
当ζ⊥→x⊥时，式（12）中的Green函数T～ij，G（ζ⊥-x⊥，-kz，kzv）存在奇异性，采用文献[21]的方法，借助辅助问题组合解，消除奇异性.3边界元离散与算例3.1 边界元离散移动荷载与孔洞相互作用的频域波数域的2.5D边界积分方程式（22），沿边界Γ的积分可采用文献[22]中常规的2D边界元法计算.若采用等参单元离散，在局部坐标系-1≤η≤1下，对于节点x有
x（η）=∑Mi=1Nl（η）x（l），
i（η）=∑Mi=1Nl（η）i，（l），
i（η）=∑Mi=1Nl（η）i，（l），（23）
式中：
M为每个等参单元的节点数.
将饱和土体内任一点xm作为场点，若相连单元数为er （ er=1或2）个，且任一单元的节点数为Ne，则将式（23）代入式（22）中，可得离散的积分方程为
1-c∞U^i（xm）=∑Nek=1∑Ml=1T^j，（kl）∫Γk（ij，G-Uij，Ga）Nl（η）J（η）dη+
∑Nek=1∑Ml=1T^j，（kl）∫ΓkUij，GaNl（η）J（η）dη-
∑Nek=1（k≠e1，…，er）∑Ml=1[U^j，（kl）∫ΓkT～ij，GNl（η）J（η）dη]-
∑Nek=e1，…，er∑Ml=1U^j，（kl）∫ΓkT～ij，GNl（η）-δnklTij，Ga]J（η）dη×
U^j（xm）∑Nek=1（k≠e1，…，er）∫ΓkTij，GaJ（η）dη，（24）
式中： c∞=1/2.求解离散后的边界积分方程可得到频域解，采用快速Fourier逆变换法（IFFT），对波数到空间域进行Fourier逆变换，可得饱和土体中孔洞动力响应的时间空间域响应.
采用IFFT进行Fourier逆变换时，样本点数N，空间域样本间距Δz及波数域样本间距Δkz应满足
Δkz=2πNΔz.（25）
当Δz≤λmin/2时（λmin为最小波长），才能保证足够大的截止频率.3.2 与已知文献结果对比为验证本方计算方法正确性，分析不同荷载速度下，半空间饱和土体表面受移动点荷载作用下的动力响应.以图1为例，图中移动点荷载幅值为Fz，以不同速度v沿z轴正方向在饱和土体表面上移动.
饱和土体参数为：
μ=3.0 GPa， a∞=2.0，
λ=3.0 GPa， a=0.95， bp=1.0×1010 （kg/m3）·s，
M=5.0 GPa， =0.3，
ρs=2 500 kg/m3，
ρf=1 000 kg/m3，
v=0.1vSH，0.5vSH，0.9vSH（vSH=μ/ρb）.
采用文中方法，边界离散仅在与荷载运动方向垂直的x方向进行，本文采用2节点等参单元进行x方向的边界离散，单元总数200个， x方向总计算长度100.0 m，观察点（x=1.0 m， y=-1.0 m）的竖向位移在-2.0 m≤z*=z-vt≤2.0 m范围内变化，如图2，图中： u*y=μuyaR/Fz ， aR=1.0 m.
由图2可知，采用文中的2.5D BEM法与文献[3]结果一致.
图2本文解与文献[3]解析解结果比较
Fig.2Comparison of the present results with
the analytical solutions of Ref. [3]
4结论（1） 不同于粗糙的近似线弹性模型的是，采用饱和半无限空间的软土地基模型可以分析交通荷载作用下，土体中的水相、土颗粒之间的相互耦合作用以及考虑土体中孔压的变化对周围土体动力响应的影响.
（2） 对于具有与荷载运动方向具有一致性的饱和土体中孔洞结构的动力响应问题，采用2.5D BEM，可将移动荷载的3D问题转化为平面内2D边界问题，使边界的离散只在与荷载移动方向相垂直的其它2个方向进行，避免了有限元法单元网格随荷载移动而破坏的难题，减少了计算工作量.
（3） 利用文中2.5D BEM，结合有限元法，还可分析饱和土体中，隧道衬砌在地铁交通荷载作用下的动力响应及地震波场散射问题.参考文献：[1]王永强，牛星钢，谭钦文. 重型车辆荷载下埋地天然气管道的安全分析[J]. 中国安全生产科学技术，2011，7（8）： 109114.
WANG Yongqiang， NIU Xinggang， TAN Qinwen. Safety analysis for buried gas pipelines under heavy vehicle loads[J]. Journal of Safety Science and Technology， 2011， 7（8）： 109114.
[2]王文和，易俊，沈士明. 基于风险的城市埋地燃气管道安全评价模型及应用[J]. 中国安全生产科学技术，2010，6（3）： 163166.
WANG Wenhe， YI Jun， SHEN Shiming. Safety evaluation model and application for the urban buried gas pipeline based on risk[J]. Journal of Safety Science and Technology， 2010， 6（3）： 163166.
[3]LU J F， JENG Dongsheng. A halfspace saturated poroelastic medium subjected to a moving point load[J]. International Journal of Solids and Structures， 2007， 44（2）： 573586.
[4]EASON G. The stress produced in a semiinfinite solid by a moving surface force[J]. International Journal of Engineering Science， 1965， 2（6）： 581609.
[5]PAYTON R G. An application of the dynamical Betti Rayleigh reciprocal theorem of moving point loads in elastic medium[J]. The Quarterly of Applied Mathematics， 1964， 20（1）： 299313.
[6]DOMINGUEZ J. Boundary element approach for dynamic poroelastic problems[J]. International Journal for Numerical Methods in Engineering， 1992， 35（2）： 307324.
[7]CHENG A H D， BADMUS D E， BESKOS D E. Integral equations for dynamic poroelasticity in frequency domain with BEM solution[J]. Journal of Engineering Mechanics， ASCE， 1991， 117（5）： 11361157.
[8]ZIMMERMAN C， STERN M. Boundary element solution of 3D wave scatter problems in a poroelastic medium[J]. Engineering Analysis with Boundary Elements， 1993， 12（4）： 223240.
[9]SCHANZ M. Application of 3D time domain boundary element formulation to wave propagation in poroelastic solids[J]. Engineering Analysis with Boundary Elements， 2001， 25（4/5）： 363376. [10]MAESO O， AZNAREZ J J， GARCIA F. Dynamic impedances of piles and groups of piles in saturated soils[J]. Computers & Structures， 2005， 83（10/11）： 769782.
[11]PAPAGEORGIOU A S， PEI D. A discrete wavenumber boundary element method for study of the 3D response of 2D scatterers[J]. Earthquake Engineering & Structural Dynamics， 1998， 27（6）： 619638.
[12]TADEU A J B， KAUSEL E. Greens functions for twoandahalfdimensional elastodynamic problems[J]. Journal of Engineering Mechanics， ASCE， 2000， 126（10）： 10931097.
[13]LUCO J E， WONG H L， BARROS F C P. Threedimensional response of a cylindrical canyon in a layered halfspace[J]. Earthquake Engineering & Structural Dynamics， 1990， 19（6）： 799817.
[14]STAMOS A A， BESKOS D E. 3D seismic response analysis of long lined tunnels in halfspace[J]. Soil Dynamics and Earthquake Engineering， 1996， 15（2）： 111118.
[15]SHENG X， JONES C J C， THOMPSON D J. Prediction of ground vibration from trains using the wavenumber finite and boundary element methods[J]. Journal of Sound and Vibration， 2006， 293（3/5）： 575586.
[16]徐志英，徐新宁. 土的非线性对土地下结构相互作用的影响[J]. 河海大学学报，1993，21（1）： 1115.
XU Zhiying， XU Xinning. The influential factors on the soilunderground structure interaction due to the nonlinear mechanical behavior of soils[J]. Journal of Hohai University： Natural Sciences， 1993， 21（1）： 1115.
[17]BIOT M A. Theory of propagation of elastic waves in a fluidsaturated porous solid， low frequency range[J]. Journal Acoustic Social America， 1956， 28： 168178.
[18]BIOT M A. Generalized theory of acoustic propagation in porous dissipative media[J]. Journal Acoustic Social America， 1962， 34（9A）： 12541264.
[19]ANDERSEN L， NIELSEN S R K. Boundary element analysis of the steadystate response of an elastic halfspace to a moving force on its surface[J]. Engineering Analysis with Boundary Elements， 2003， 27（1）： 2338.
[20]RASMUSSEN K M， NIELSEN S R K， KIRKEGAARD P H. Boundary element method in time domain for a moving timedependent force[J]. Computer & Structure， 2001， 79（7）： 69701.
[21]NORRIS A N. Radiation from a point source and scattering theory in a fluidsaturated porous medium[J]. Journal Acoustic Social America， 1985， 77（2）： 20122023.
[22]姜弘道. 弹性力学问题的边界元法[M]. 北京：中国水利水电出版社，2008： 6070.
关键词： 格林函数；移动荷载；饱和土；Fourier积分变换
中图分类号： TU 471文献标志码： AAnalysis of Dynamic Interaction between Hole Embedded in
Saturated Soil and Moving Loads
Using 2.5D Boundary Element MethodXU Bin1，2，GAO Liang2，LEI Xiaoyan3，XU Manqing1，LIU Linya3
（1. Department of Civil Engineering， Nanchang Institute of Technology， Nanchang 330029， China； 2. School of Civil Engineering and Architecture， Beijing Jiaotong University， Beijing 100044， China； 3. Engineering Research Center of Railway Environmental Vibration and Noise， Ministry of Education， East China Jiaotong University， Nanchang 330013， China）
Abstract：To study the mechanism of interaction between the saturated soil with complex surface and moving loads， a specific 2.5D boundary element method （BEM） for the problem of dynamic responses of the hole embedded in the saturated soil under a moving load was derived systematically. Based on Biots theory， the frequency wave number domain Greens function for saturated porous media was developed using the Fourier transform method and the potential decomposition approach. Then， dynamic responses in the timespace domain solution were further obtained by synthesizing the wave number solution via inverse fast Fourier transform （IFFT）. The results show that using the frequency wavenumber boundary integrated equations established by the consistency of the directions of moving load， a 3D spatial model for the dynamic interaction between the holes embedded in the saturated soil with moving loads is convertible to a 2D planar model in the frequency domain， and the 3D complex boundary can be simplified to a problem of 2D boundary element. Thus， the complicated boundary planes can be transformed into computed lines using the 2.5D wave number domain BEM formulation， consequently reducing the amount of computation to a manageable size.
Key words：Greens function； moving loads； saturated soil； Fourier integral transform
工程中穿越软土地基或道路下的市政供水、天然气等管道不可避免受路面交通荷载作用，当埋地管道一旦达到破坏强度就可能产生安全事故，造成水、天然气等泄漏而引发燃烧与爆炸[12].上述工程问题均可归结于移动荷载与软土地基中孔洞结构相互作用的动力响应模型.
对于具有简单几何构形的半空间土体表面受移动载荷问题，可利用解析方法进行计算[35]，但对于移动荷载作用下，地下孔洞结构的动力响应，边界较复杂，其几何轴对称性已破坏，难以采用解析法计算.特别是移动荷载作用下，土体中各点动力响应都随时间和空间不同而改变，采用有限元法对土体内部单元网格划分时，由于时间和空间的变化，不可避免造成单元数量急剧增加，计算困难.边界元方法（BEM）具有自动满足无穷远处辐射条件及降维功能，已被证明是解决具有复杂边界条件问题一种有效手段，并且在工程中广泛应用.文献[6]中以3个土骨架位移和孔隙水压力为基本变量，建立了饱和土体边界积分方程.采用相同西南交通大学学报第48卷第4期徐斌等：移动荷载与土体中孔洞相互作用的2.5D边界元法分析的基本变量，文献[7]中推导了3D饱和土体频域内边界积分方程.文献[8]中利用3D频域边界积分方程，分析了饱和多孔介质中散射问题.文献[9]中利用拉普拉斯域的Green函数和卷积积分法，推导了时域内边界积分方程方法.采用3D边界元法，文献[10]中对饱和土体中群桩土共同作用的动力问题进行了分析. 与传统的有限元法相比，边界元法只对边界进行单元离散，因此，建模单元、计算时间均可节省.然而， 3D时间域内边界元法仍然需占用计算机大量内存和计算时间.文献[1112]在3D边界元的基础上，将结构作为2D分析或采用空间波数变换用来简化问题，提出了波数域的2.5D边界元法.与3D边界元法相比， 2.5D边界元法的边界离散和计算要求大大减少.目前2.5D边界元法已成功应用于弹性介质中的动力响应问题分析.文献[13]根据具有恒定移动速度的点荷载Green函数，推导了任意角入射波下的无限长河谷3D响应的2.5D间接边界积分方程.利用全空间Green函数，文献[14]采用直接2.5D BEM，分析了粘弹性半空间中无限长隧道衬砌的3D地震荷载作用下的动力响应.采用离散波数的有限元与边界元耦合，并结合傅里叶逆变换法.文献[15]对轨道交通荷载引起地面3D空间的动力响应进行了评价.
上述文献的研究表明：土体为单相弹性介质，忽略土体在移动载荷作用下土体中的水相、土颗粒之间的相互耦合作用，则无法考虑土体中孔压的变化，以及孔压变化对周围土体动力响应的影响[3，16].多孔软土中孔隙水压的变化对土体强度有较大的影响，甚至引起震陷、液化等工程现象，这些都使得软土地基具有非线性、非弹性的动力本质特征[1718].传统的线弹性、粘弹塑性土体本构模型难以有效反映土体变形，特别是孔隙水压力变化而导致的孔洞动态受力性状.目前关于饱和土体中孔洞与移动荷载相互作用的2.5D BEM法尚未见报道.本文分别推导了由单位集中力作用在土骨架时所对应Green函数，以及单位散度源作用在孔隙介质流体时所对应Green函数.考虑荷载移动方向与饱和土体中孔洞轴线方向具有一致性，在三维边界元法（3DBEM）的基础上，采用沿饱和土体中孔洞纵向进行空间波数的积分变换，形成2.5D BEM，将三维移动荷载问题转化为平面内边界元积分问题，既可避免移动荷载引起的半无限空间有限元单元网格划分随时间而变化的难题，又可避免3D BEM在移动荷载方向的边界划分问题，能够减少计算工作量.1饱和二相介质土Green函数为推导饱和地基土体模型的Green函数，定义时间t到频域ω，空间z到波数域kz的Fourier变换对为
f-（ω）=∫+∞-∞f（t）e-iωtdt，
f（t）=12π∫+∞-∞f-（ω）eiωtdω，
f^（kz）=∫+∞-∞f（z）eikzzdz，
f（z）=12π∫+∞-∞f^（kz）e-ikzzdkz，（1）
式中：
“-”、“^”分别表示频域、波数域量.
由于饱和多孔介质由土骨架及水相组成， Green函数张量G可分别由土骨架在单位集中力作用下所对应Green函数Gs、Gs以及饱和孔隙介质流体在单位散度源作用时所对应Green函数Gf、Gf 组成，表达式为
G=GsGf
GsGf.（2）
根据Biot理论[1718]，对于饱和二相介质，土骨架与水相分别在单位体积力Fi和fi（i=1，2，3）作用下，土体骨架位移ui（i=1，2，3）及孔隙水压（pf）4个独立变量表示的频域内饱和土体运动方程为
μi，jj+（λ+μ）j，ji+ρgω2i-αgf，i=
-（i-β4f-i），
f，jj+β2ω2f-β3j，j=f-i，i，（3）
式中： m=α∞ρf/；
αg=α-β4；
ρg=ρb-β4ρf；
ρb=（1-）ρs+ρf；
β1=M/[mω2-iω（η/k）]；
β2=1/（β1ω2）；
β3=ρfω2-α[mω2-iω（η/k）]；
β4=ρfω2β1/M，
其中， λ和μ为Lame常数， α、M分别为饱和土体中土骨架、水相的压缩系数， 为饱和多孔介质的孔隙率， wi（i=1，2，3）为饱和土体中流体相对于土骨架的位移， ρs、ρf分别为饱和土体土骨架、孔隙水密度， α∞为孔隙介质弯曲系数， k、η分别为饱和土体水相的渗透系数及孔隙流体的黏质度， ω为移动荷载频率.
饱和多孔介质中土骨架在集中力作用下，对应Green函数的Gs、Gs，可由式（3）求得，
μΔ 2Gs+（λ+μ）Δ （Δ Gs）+
ρgω2Gs-αgGs=-δ（R）I3×3，
Δ 2Gs+β2ω2Gs-β3Δ Gs=0，（4）
式中：
Δ为Laplace算子.
由文献[19]可知，可采用3个标量势函数f，Gs、s，Gs、f，Gs表示Gs、Gs，则有
Gs=Δ Δ f，Gs+Δ Δ s，Gs+
（Δ 2f，GsI3×3-Δ Δ f，Gs），
Gs=AfΔ Δ 2f，Gs+AsΔ Δ 2s，Gs.（5）
利用恒等式
-δ（R）I3×3=Δ Δ 14πR+
Δ 2I3×3-Δ Δ 14πR，（6）
式中：
R为两点的空间距离.
将式（5）代入式（4）中，同时利用式（6），可得
[afΔ 2f，Gs+ρgω2Δ 2f，Gs]+
[asΔ 2s，Gs+ρgω2Δ 2s，Gs]=14πR，
[AfΔ 2f，s+dff，Gs]+
[AsΔ 2s，Gs+dss，Gs]=0，
Δ 2Gs+k2tGs=14πμR，（7）
式中：
af=λ+2μ-αgAf；
as=λ+2μ-αgAs； df=β2Afω2-β3；
ds=β2Asω2-β3；
k2f=ρgω2/（λ+2μ-αgAf）；
k2s=ρgω2/（λ+2μ-αgAs）；
k2t=ρgω2/μ；
Re（kf）≤Re（ks）.
利用式（1），对式（7）进行z→kz的Fourier变换，整理后可得频域、波数域内的势函数为
f，Gs=iAs4k2f（afAs-asAf）[H0，（2）（γfr）-H0，（2）（-ikzr）]，
s，Gs=iAf4k2s（asAf-afAs）[H0，（2）（γsr）-H0，（2）（-ikzr）]，
Gs=i4ρgω2[H0，（2）（γtr）-H0，（2）（-ikzr）]，（8）
式中： “～”表示频域波数内的量；
H0，（2）为第二类的0阶Hankel函数；
r为两点的平面距离.
根据式（1）进行z→kz的Fourier变换，由式（5）可得频域、波数域内的Green函数Gs、Gs，再将式（8）代入后，可得在集中力作用下饱和多孔介质中土骨架所对应的频域波数域内Green函数为
ij，Gs=[δf，（s）H0，（2）（γfr）+δs，（s）H0，（2）（γsr）-δt，（s）H0，（2）（γtr）]，ij-
[δt，（s）H0，（2）（γtr）]δij，
i，Gs=χp[-H0，（2）（γfr）+H0，（2）（γsr）]，ii，j=1，2，3，（9）
式中：
δf，（s）=iAs/[4k2f（afAs-asAf）]；
δs，（s）=iAf/[4k2s（asAf-afAs）]；
δt，（s）=i/（4ρgω2）；
χp=-Ask2fδs，（s）.
利用式（3），可得单位散度源作用在饱和土体中流体介质时所对应Green函数Gf、Gf为
μΔ 2Gf+（λ+μ）Δ （Δ Gf）+
ρgω2Gf-αgGf=0，
Δ 2Gf+β2ω2Gf-β3Δ Gf=-δ（R）.（10）
由于流体介质中不存在剪切波，因此可用2个标量函数f，Gf、s，Gf表示相应的Green函数Gf、Gf，
Gf=Δ f，Gf+Δ s，Gf，
Gf=AfΔ 2f，Gf+AsΔ 2s，Gf.（11）
同单位集中力作用在土骨架时所对应的频域波数域内Green函数类似的步骤，可得单位散度源作用在饱和土体流体介质时，所对应的频域波数域内Green函数为
i，Gf=δ（f）[H0，（2）（γfr）-H0，（2）（γsr）]，i，
Gf=δ（f）[-Afk2fH0，（2）（γfr）+
Ask2sH0，（2）（γsr）]，i=1，2，3，（12）
式中：
δf=ias/[4k2f（asAf-afAs）].2饱和土体中2.5D边界积分方程2.1饱和土体频域内的边界积分方程根据互易定理及饱和土体的本构关系（式（3）），考虑饱和土体任意2个状态（状态1和状态2），可得
ij，（1）ij，（2）+f，（1）ξ（2）=
ij，（2）ij，（1）+f，（2）ξ（1），（13）
式中：下标（1）、（2）分别表示饱和土体的2种运动状态；
σij为饱和土体总应力分量；
ξ为饱和多孔介质中流体的单位体积增量，
ξ=wi，i.
将式（13）沿饱和土体任意区域Ω积分，并结合式（3），可得
∫S[j，（1）j，（2）-j，（2）j，（1）]dS-
∫S[f，（1）n，（2）-f，（2）n，（1）]dS=
∫Ω[f，（1）f，（2）-f，（2）f，（1）]dΩ-
∫Ω[gj，（1）j，（2）-gj，（2）j，（1）]dΩ，（14）
式中：
j为作用在区域Ω边界S上的面力，
j=ijni；
ni为边界S外法线方向余弦；
f=β1i，i/M；
n为沿边界S的外法线方向流体位移，
n=β1（f，j-ρfω2j）nj/M.
将饱和土体的真实解作为状态1，单位集中力或单位散度源所应于Green函数解为状态2，则在不考虑体积力时，由式（14）可得在饱和土体土骨架任意点x上作用一集中力时，饱和土体的Somigliana等式为
i（x，ω）=∫S[ij，Gs（x，ζ，ω）j（ζ，ω）-
ij，Gs（x，ζ，ω）j（ζ，ω）]dS（ζ）-
∫S[ni，Gs（x，ζ，ω）f（ζ，ω）-
i，Gs（x，ζ，ω）n（ζ，ω）]dS（ζ），（15）
式中：
ni，Gs=β1[i，Gs（x，ζ，ω）/ζj-
ρfω2ij，Gs]nj（ζ）/M，
ij，Gs=ijk，Gs（x，ζ，ω）nk（ζ），
其中， Gs、Gs分别为单位集中力作用在土骨架时所对应Green函数.
利用相同方法，饱和土体中流体介质受单位散度源作用时的Somigliana等式为
i（x，ω）=∫S[j，Gf′（x，ζ，ω）j（ζ，ω）-
j，Gf′（x，ζ，ω）j（ζ，ω）]dS（ζ）-
∫S[n，Gf′（x，ζ，ω）f（ζ，ω）- Gf′（x，ζ，ω）n（ζ，ω）]dS（ζ），（16）
式中：
j，Gf′=β1j，Gf/M；
j，Gf′=β1j，Gf/M；
n，Gf′=β1n，Gf/M；
Gf′=β1Gf/M.
根据式（2）及式（15）、（16）可得到饱和土体统一的Somigliana等式为
i（x，ω）=∫S[ij，G（x，ζ，ω）j（ζ，ω）-
ij，G（x，ζ，ω）j（ζ，ω）]dS（ζ）.（17）
利用文献[20]的方法，考虑不同边界的特性，可得频域内饱和土体的边界积分方程为
ciji（x，ω）=∫S[ij，Gs（x，ζ，ω）j（ζ，ω）-
ij，Gs（x，ζ，ω）j（ζ，ω）]dS（ζ），（18）
在光滑边界时， cij=δij/2.2.2 频域波数域的2.5D边界积分方程移动荷载与饱和土体中的孔洞相互作用模型如图1所示.
图1移动点荷载与饱和土体相互作用示意图
Fig.1Schematic diagram of the interaction between
the moving point loads and the saturated soil若荷载沿孔洞轴向以速度v移动，饱和土体及其中孔洞的动力响应均可表示为
ui（x，t）=ui（x⊥，z-vt），
ti（x，t）=ui（x⊥，z-vt）.（19）
对式（19）进行时间、空间Fourier变换，可得
i（x⊥，kz，ω）=
2πδ（ω+ω0-kzv）U^i（x⊥，kz），
i（x⊥，kz，ω）=
2πδ（ω+ω0-kzv）T^i（x⊥，kz），（20）
式中：
ω0为移动荷载初始频率；
U^i（x⊥，kz）=∫∞-∞ui（x⊥，z-vt）eikz（z-vt）d（z-vt），
T^i（x⊥，kz）=∫∞-∞ti（x⊥，z-vt）eikz（z-vt）d（z-vt），
其中， x⊥为平面内坐标点.
将式（20）代入式（18）中，可得到频率波数域的饱和土体中孔洞的2.5D边界元积分方程为
cijU^j（x⊥，kz）=
∫Γij，G（ζ⊥-x⊥，-kz，kzv）T^j（x⊥，kz）dΓ（ζ⊥）-
∫ΓT～ij，G（ζ⊥-x⊥，-kz，kzv）U^j（x⊥，kz）dΓ（ζ⊥），（21）
式中：
Γ 为饱和土体及孔洞在 平面内的边界；
ζ⊥为 平面内任意观察点坐标.
对于饱和土体-孔洞体系内任一点，有
U^j（x⊥，kz）=
∫Γij，G（ζ⊥-x⊥，-kz，kzv）T^j（x⊥，kz）dΓ（ζ⊥）-
∫ΓT～ij，G（ζ⊥-x⊥，-kz，kzv）U^j（x⊥，kz）dΓ（ζ⊥）.（22）
当ζ⊥→x⊥时，式（12）中的Green函数T～ij，G（ζ⊥-x⊥，-kz，kzv）存在奇异性，采用文献[21]的方法，借助辅助问题组合解，消除奇异性.3边界元离散与算例3.1 边界元离散移动荷载与孔洞相互作用的频域波数域的2.5D边界积分方程式（22），沿边界Γ的积分可采用文献[22]中常规的2D边界元法计算.若采用等参单元离散，在局部坐标系-1≤η≤1下，对于节点x有
x（η）=∑Mi=1Nl（η）x（l），
i（η）=∑Mi=1Nl（η）i，（l），
i（η）=∑Mi=1Nl（η）i，（l），（23）
式中：
M为每个等参单元的节点数.
将饱和土体内任一点xm作为场点，若相连单元数为er （ er=1或2）个，且任一单元的节点数为Ne，则将式（23）代入式（22）中，可得离散的积分方程为
1-c∞U^i（xm）=∑Nek=1∑Ml=1T^j，（kl）∫Γk（ij，G-Uij，Ga）Nl（η）J（η）dη+
∑Nek=1∑Ml=1T^j，（kl）∫ΓkUij，GaNl（η）J（η）dη-
∑Nek=1（k≠e1，…，er）∑Ml=1[U^j，（kl）∫ΓkT～ij，GNl（η）J（η）dη]-
∑Nek=e1，…，er∑Ml=1U^j，（kl）∫ΓkT～ij，GNl（η）-δnklTij，Ga]J（η）dη×
U^j（xm）∑Nek=1（k≠e1，…，er）∫ΓkTij，GaJ（η）dη，（24）
式中： c∞=1/2.求解离散后的边界积分方程可得到频域解，采用快速Fourier逆变换法（IFFT），对波数到空间域进行Fourier逆变换，可得饱和土体中孔洞动力响应的时间空间域响应.
采用IFFT进行Fourier逆变换时，样本点数N，空间域样本间距Δz及波数域样本间距Δkz应满足
Δkz=2πNΔz.（25）
当Δz≤λmin/2时（λmin为最小波长），才能保证足够大的截止频率.3.2 与已知文献结果对比为验证本方计算方法正确性，分析不同荷载速度下，半空间饱和土体表面受移动点荷载作用下的动力响应.以图1为例，图中移动点荷载幅值为Fz，以不同速度v沿z轴正方向在饱和土体表面上移动.
饱和土体参数为：
μ=3.0 GPa， a∞=2.0，
λ=3.0 GPa， a=0.95， bp=1.0×1010 （kg/m3）·s，
M=5.0 GPa， =0.3，
ρs=2 500 kg/m3，
ρf=1 000 kg/m3，
v=0.1vSH，0.5vSH，0.9vSH（vSH=μ/ρb）.
采用文中方法，边界离散仅在与荷载运动方向垂直的x方向进行，本文采用2节点等参单元进行x方向的边界离散，单元总数200个， x方向总计算长度100.0 m，观察点（x=1.0 m， y=-1.0 m）的竖向位移在-2.0 m≤z*=z-vt≤2.0 m范围内变化，如图2，图中： u*y=μuyaR/Fz ， aR=1.0 m.
由图2可知，采用文中的2.5D BEM法与文献[3]结果一致.
图2本文解与文献[3]解析解结果比较
Fig.2Comparison of the present results with
the analytical solutions of Ref. [3]
4结论（1） 不同于粗糙的近似线弹性模型的是，采用饱和半无限空间的软土地基模型可以分析交通荷载作用下，土体中的水相、土颗粒之间的相互耦合作用以及考虑土体中孔压的变化对周围土体动力响应的影响.
（2） 对于具有与荷载运动方向具有一致性的饱和土体中孔洞结构的动力响应问题，采用2.5D BEM，可将移动荷载的3D问题转化为平面内2D边界问题，使边界的离散只在与荷载移动方向相垂直的其它2个方向进行，避免了有限元法单元网格随荷载移动而破坏的难题，减少了计算工作量.
（3） 利用文中2.5D BEM，结合有限元法，还可分析饱和土体中，隧道衬砌在地铁交通荷载作用下的动力响应及地震波场散射问题.参考文献：[1]王永强，牛星钢，谭钦文. 重型车辆荷载下埋地天然气管道的安全分析[J]. 中国安全生产科学技术，2011，7（8）： 109114.
WANG Yongqiang， NIU Xinggang， TAN Qinwen. Safety analysis for buried gas pipelines under heavy vehicle loads[J]. Journal of Safety Science and Technology， 2011， 7（8）： 109114.
[2]王文和，易俊，沈士明. 基于风险的城市埋地燃气管道安全评价模型及应用[J]. 中国安全生产科学技术，2010，6（3）： 163166.
WANG Wenhe， YI Jun， SHEN Shiming. Safety evaluation model and application for the urban buried gas pipeline based on risk[J]. Journal of Safety Science and Technology， 2010， 6（3）： 163166.
[3]LU J F， JENG Dongsheng. A halfspace saturated poroelastic medium subjected to a moving point load[J]. International Journal of Solids and Structures， 2007， 44（2）： 573586.
[4]EASON G. The stress produced in a semiinfinite solid by a moving surface force[J]. International Journal of Engineering Science， 1965， 2（6）： 581609.
[5]PAYTON R G. An application of the dynamical Betti Rayleigh reciprocal theorem of moving point loads in elastic medium[J]. The Quarterly of Applied Mathematics， 1964， 20（1）： 299313.
[6]DOMINGUEZ J. Boundary element approach for dynamic poroelastic problems[J]. International Journal for Numerical Methods in Engineering， 1992， 35（2）： 307324.
[7]CHENG A H D， BADMUS D E， BESKOS D E. Integral equations for dynamic poroelasticity in frequency domain with BEM solution[J]. Journal of Engineering Mechanics， ASCE， 1991， 117（5）： 11361157.
[8]ZIMMERMAN C， STERN M. Boundary element solution of 3D wave scatter problems in a poroelastic medium[J]. Engineering Analysis with Boundary Elements， 1993， 12（4）： 223240.
[9]SCHANZ M. Application of 3D time domain boundary element formulation to wave propagation in poroelastic solids[J]. Engineering Analysis with Boundary Elements， 2001， 25（4/5）： 363376. [10]MAESO O， AZNAREZ J J， GARCIA F. Dynamic impedances of piles and groups of piles in saturated soils[J]. Computers & Structures， 2005， 83（10/11）： 769782.
[11]PAPAGEORGIOU A S， PEI D. A discrete wavenumber boundary element method for study of the 3D response of 2D scatterers[J]. Earthquake Engineering & Structural Dynamics， 1998， 27（6）： 619638.
[12]TADEU A J B， KAUSEL E. Greens functions for twoandahalfdimensional elastodynamic problems[J]. Journal of Engineering Mechanics， ASCE， 2000， 126（10）： 10931097.
[13]LUCO J E， WONG H L， BARROS F C P. Threedimensional response of a cylindrical canyon in a layered halfspace[J]. Earthquake Engineering & Structural Dynamics， 1990， 19（6）： 799817.
[14]STAMOS A A， BESKOS D E. 3D seismic response analysis of long lined tunnels in halfspace[J]. Soil Dynamics and Earthquake Engineering， 1996， 15（2）： 111118.
[15]SHENG X， JONES C J C， THOMPSON D J. Prediction of ground vibration from trains using the wavenumber finite and boundary element methods[J]. Journal of Sound and Vibration， 2006， 293（3/5）： 575586.
[16]徐志英，徐新宁. 土的非线性对土地下结构相互作用的影响[J]. 河海大学学报，1993，21（1）： 1115.
XU Zhiying， XU Xinning. The influential factors on the soilunderground structure interaction due to the nonlinear mechanical behavior of soils[J]. Journal of Hohai University： Natural Sciences， 1993， 21（1）： 1115.
[17]BIOT M A. Theory of propagation of elastic waves in a fluidsaturated porous solid， low frequency range[J]. Journal Acoustic Social America， 1956， 28： 168178.
[18]BIOT M A. Generalized theory of acoustic propagation in porous dissipative media[J]. Journal Acoustic Social America， 1962， 34（9A）： 12541264.
[19]ANDERSEN L， NIELSEN S R K. Boundary element analysis of the steadystate response of an elastic halfspace to a moving force on its surface[J]. Engineering Analysis with Boundary Elements， 2003， 27（1）： 2338.
[20]RASMUSSEN K M， NIELSEN S R K， KIRKEGAARD P H. Boundary element method in time domain for a moving timedependent force[J]. Computer & Structure， 2001， 79（7）： 69701.
[21]NORRIS A N. Radiation from a point source and scattering theory in a fluidsaturated porous medium[J]. Journal Acoustic Social America， 1985， 77（2）： 20122023.
[22]姜弘道. 弹性力学问题的边界元法[M]. 北京：中国水利水电出版社，2008： 6070.