饱和土中圆柱形衬砌对瞬态弹性波的散射①

王 滢，高 盟，2，高广运，张先林

（1．同济大学 岩土及地下工程教育部重点实验室，上海 200092；2．山东科技大学 山东省土木工程防灾减灾重点实验室，山东 青岛 266510；3．上海市规划和国土资源局，上海 200003）

0 前言

地下隧道、各类埋置的运输管道、地下石油储运设施等重要工程常常会受到弹性波的作用，波与结构间的动力相互作用将影响到这些地下工程的安全。因此该类地下结构在弹性波作用下的动力响应研究日益受到重视。

通常，这类问题分为两类：稳态问题和瞬态问题。众多学者对稳态弹性波入射下地下结构的动力响应进行了研究。如Pao和 Mow［1］采用波函数展开法开创性地研究了无限空间中洞室在弹性波入射下的动应力集中；Lee和 Trifunac［2－4］，Luco等［5］和Davis等［6］运用复变函数法给出了半空间中无壳洞室对P和SV波散射问题的解析解；Karinski等［7］分析了全空间饱和介质中带有衬砌的洞室在平面波作用下的动应力集中；刘殿魁等［8］研究界面上的圆形衬砌结构对平面SH波散射与动应力集中，给出了关于圆形衬砌结构界面上动应力集中系数的数值结果，并讨论了圆形衬砌结构界面的动应力集中系数的影响因素；周香莲等［9］运用复变函数法研究了饱和土中的圆形衬砌结构对弹性波的散射和动应力集中，得到饱和土的位移、应力和孔压的表达及衬砌结构的位移和应力表达，对波数及衬砌结构厚度等参数进行了分析。刘干斌等［10］给出了全空间中半渗透圆柱形壳结构对平面P波散射的解析解；丁光亚等［11］在Biot波动理论的基础上，引入更符合工程实际的半透水边界条件，给出了半空间均质饱和土中半渗透圆柱形壳结构对平面P波散射的级数解。

然而，关于地下结构对弹性波的散射问题的研究主要集中于稳态弹性波，对于瞬态弹性波则少有涉及。Baron［12］首先研究了弹性介质中圆柱形壳体对瞬态压力波的散射；Yoshihara等［13］和Paul等［14］求解了平面弹性波与圆柱形衬砌的瞬态相互作用问题；Garnet等［15］研究了弹性均匀介质中任意厚度圆形衬砌的瞬态反应；Peralta［16］给出了瞬态弹性波散射问题的解答的近似求解方法；Kobayashi和Nishimura［17］应用边界积分方程法研究了二维瞬态波的散射；Pao等［18］运用广义射线法讨论了柱体的瞬态波散射。上述研究主要限于弹性地基的情况，有关饱和地基中弹性波与衬砌的瞬态相互作用问题，国内外还未见有文献报道。

针对现有研究的不足，本文运用Laplace变换和波函数展开法，求解单位阶跃弹性波入射条件下饱和土中圆柱形衬砌的动力响应解答。利用La－place逆变换的数值方法，给出问题的数值解，并分析衬砌结构的动应力集中系数的波型特性和材料剪切模量、衬砌厚度对动应力集中系数的影响。

1 饱和土中波场的求解

考虑一个沿x正向传播的单位阶跃入射波扰动，如图1所示。入射波可以表示为

式中，φ0为入射波幅值；H［t－（x＋a1）／cp］为单位阶跃函数，其中t＝0是入射波首先到达x＝－a1的时间；cp为入射波波速；a1为圆柱形衬砌外半径；a2为圆柱形衬砌内半径；a2－a1＝h为衬砌厚度。

图1 计算模型Fig．1 Analysis model．

对式（1）进行Laplace变换得

为计算方便，将式（2）改写成级数的形式：

式中，p为拉普拉斯变换参数；In（·）为第一类变形的Bessel函数；β＝p／cp；当n＝0时on＝1，当n≥1时on＝2。

1．1 散射波场的求解

入射波遇到衬砌结构后，在饱和土中产生一个向外传播的散射波场。

在动力荷载作用下，以直角坐标表示，土体的本构关系为［19］

渗流连续方程为

式中，σij为土体的总应力；εij、e分别为土骨架的应变和体积应变；σf为孔压；A、N为土骨架的Lame常数；δij为Kronecker符号；ϑ为单位体积中流体量的改变量；α和M为Biot参数。

土骨架和流体的运动方程为

式中，ui、wi分别为土骨架位移和流体相对于土骨架的位移；ρ、ρf分别为土体总密度及流体密度，ρ＝（1－n）ρs＋ρf，ρs为土骨架固体材料密度；n为土体的孔隙率；η为孔隙流体的动力粘滞系数；kd为动力渗透系数。

式中，eijk为直角坐标系下的置换张量。

将式（9）和式（10）代入式（7）和式（8）可得

式中 ▽2＝∂／∂x2＋∂／∂y2，p为拉普拉斯变换参数。

由式（13）可得

由式（14）可得

式中，

将式（22）进一步分解可得

式（23）和式（24）为修正的 Helmholtz方程，根据分离变量法，可得其通解为

同理，由式（15）～式（17）可得

式中，上标s表示散射波；Kn（·）为第二类变形的Bessel函数。

将式（27）、（28）代入式（13）的第（1）式和第（2）式可得

式中，an、bn、cn、dn、en为待定系数。

1．2 饱和土中的总波场

饱和土中总波场由入射波场和散射波场构成，即

在极坐标系下，饱和土体中势函数与位移、应力和孔压的关系为

2 衬砌结构波场的求解

将衬砌结构视为弹性均匀介质，则衬砌结构存在两种体波：压缩波和剪切波。在直角坐标系下，衬砌结构的控制方程可表示为

式中，ρ1为衬砌结构的密度；A1，N1为衬砌结构的Lamé常数。

引入两个衬砌结构的势函数φ、ψ，则衬砌结构位移ui可表示为

将式（43）代入式（42）可得

极坐标系下，式（45）的解答可表示为

式中fn，gn，mn，nn为待定系数。

根据衬砌结构的几何方程、本构关系及平衡方程可得由式（46）表示的位移、应力表达式：

3 边界条件

在衬砌结构的外边界上，即r＝a1处，根据饱和土与衬砌结构结合部分的连续性条件可得

式中，下标Ⅰ、Ⅱ分别表示衬砌结构与饱和土的变量。

由于衬砌结构不透水，接触面上流体相对于土骨架的法向位移为零：

由上述条件可得

式中系数见附录。至此，待定系数an，bn，cn，dn，en，fn，gn，mn，nn已完全确定。

4 时域解

以上推导了弹性土中圆柱形孔洞内源问题频域解，为在物理空间获得时域解，只须对以上各式施以Laplace逆变换。直接进行Laplace逆变换得到解析解是困难的，而必须借助于数值方法。本文采用Durbin［20］的Laplace数值逆变换的方法。该方法可表示为

NSUM的值由以下的收敛准则确定：

式中：T为所求的时间间隔；通常ε＝10－6～10－10，aT＝5～10，T＝20，NSUM＝50～5 000就可以得到精度较高的结果。本文取aT＝10，NSUM＝50。

5 数值结果与参数分析

定义动应力集中系数σ＊为饱和土体与衬砌结构结合面或衬砌结构内边界上的切应力与入射应力的最大幅值的比值，即

饱和土介质参数取值：土体的孔隙率n＝0．40；土骨架密度ρs＝2 500kg／m3，流体密度ρf＝1 000 kg／m3；土体剪切模量G＝1．0×107Pa；土体泊松比μ＝0．4；流体动力粘滞系数η＝1．0×10－2Pa·s，kd＝1×10－7m2；Biot参数α＝0．999；M＝1．0×108Pa。

衬砌结构的计算参数如下：内半径a2＝3 000 mm，外半径a1＝3 300mm，即衬砌厚度h＝300 mm；衬砌材料密度ρL＝2 800kg／m3；衬砌材料的泊松比μL＝0．2；衬砌材料的剪切模量GL＝2．0×108Pa。

5．1 动应力集中系数的波型和综合特性

由式（46a）、式（46b）及式（50）可得

可以明显的看到，式（62）所表示的瞬态应力可以看作是各个波型，即n＝0，1，2，…时的反应所产生的应力的总和。

图2 波型随阶次n的变化Fig．2 Variation of the wave type with order n．

图2给出了θ＝π／2时，n＝1，2，3，4，5，6，7时应力的波型。由图可知，除n＝1外，所有波型的影响主要都是来自“早期”。如果对时间t进行无量纲化处理，即τ＝cst／a1，则发现所有波型的主要贡献均来自τ＜1。随着n的阶次的增大，波型呈现出明显的衰减趋势，当n＞7时，波型基本趋近于水平状态。

图3 迎波面和背波面波型特性Fig．3 Chracteristics of front and back wave surface．

图3给出了θ＝π／4，θ＝π／2，θ＝3π／4，θ＝π，θ＝5π／4时的波型。显然最先到达衬砌周边的波是在θ＝π处，此时τ＝0。到达θ＝π／2的时间是τ＝1，因此迎波面的动应力集中系数最大。

5．2 剪切模量对动应力集中系数的影响

图4 剪切模量对衬砌动应力集中的影响Fig．4 The influence of shear modulus on dynamic stress concentration factor on the lining．

保持其它参数不变，分别计算：（1）土体剪切模量G＝1．0×107Pa，衬砌剪切模量GL＝2．0×108Pa；（2）土体剪切模量G＝2．0×108Pa，衬砌剪切模量GL＝1．0×108Pa两种情况下衬砌结构的动应力集中系数。图4给出了θ＝π／2时两种情况下动应力集中系数的波型。从图中可以看出，当G＜GL，即土体较衬砌结构软时，动应力集中系数较大；而当G＞GL，即土体较衬砌结构硬时，动应力集中系数显著减小。

5．3 衬砌厚度对动应力集中系数的影响

保持其它参数不变，分别取h＝200mm、250 mm、300mm。图5为不同衬砌厚度动应力集中系数波型。由图可知，衬砌厚度越大，动应力集中系数越小。显然增大衬砌厚度是减小动应力集中的有效手段。

图5 衬砌厚度对衬砌动应力集中的影响Fig．5 Influence of lining＇s thickness on dynamic stress concentratio．

6 结论

采用Laplace变换和波函数展开法，求解了单位阶跃弹性波入射条件下，饱和土中圆柱形衬砌的动力响应解答，研究了衬砌结构的动应力集中系数的波型特性及材料剪切模量和衬砌厚度对动应力集中系数的影响，结论如下：

（1）衬砌结构的瞬态应力可以看作是各个波型，即n＝0，1，2，…时的反应所产生的应力的总和，除n＝1外，所有波型的影响主要都是来自τ＜1。随着n的阶次的增大，波型呈现出明显的衰减趋势。迎波面的动应力集中最显著。

（2）当G＜GL，即土体较衬砌结构软时，动应力集中系数较大；而当G＞GL，即围岩较衬砌结构硬时，动应力集中系数显著减小。

（3）衬砌厚度越大，动应力集中系数越小，增大衬砌厚度是减小动应力集中的有效手段。

［1］ Pao Y H，Mow C C．Diffraction of elastic waves and dynamic stress concentrations［M］．New York：Crane and Russak，1973．

［2］ Lee V W，Trifunac M D．Response of tunnels to incident SH－waves［J］．Journal of Engineering Mechanics，ASCE，1979，105（4）：643－659．

［3］ Cao H，Lee V W．Scattering and diffraction of plane P waves by circular cylindrical canyons with variable depth－to－width ratio［J］．Soil Dynamics and Earthquake Engineering，1990，9（3）：141－150．

［4］ Lee V W，Karl J．On deformations near a circular underground cavity subjected to incident plane P waves［J］．European Earthquake Engineering，1993，7（1）：29－41．

［5］ Luco J E，Debarros F C P．Dynamic displacements and stresses in the vicinity of a cylindrical cavity embedded in a halfspace［J］．Earthquake Engineering and Structural Dynamics，1994，23（33）：321－340．

［6］ Davis C A，Lee V W，Bardet J P．Transverse response of underground cavities and pipes to incident SV waves［J］．Earthquake Engineering and Structural Dynamics，2001，30（3）：383－410．

［7］ Karinski Y S，Shershnev V V，Yankelevsky D Z．Analytical solution of the harmonic waves diffraction by a cylindrical lined cavity in poroelastic saturated medium［J］．International Journal for Numerical and Analytical Methods in Geomechanics，2007，31（5）：667－689．

［8］ 刘殿魁，史守峡．界面上圆形衬砌结构对平面SH波的散射［J］．力学学报，2002，34（5）：796－803．

［9］ 周香莲，周光明，王建华．饱和土中圆形衬砌结构对弹性波的散射［J］．岩石力学与工程学报，2005，24（9）：572－1576．

［10］ Liu Ganbin，Xie Kanghe，Yao Hailin．Scattering of elastic wave by apartially sealed cylindrical shell embedded in poroelastic medium［J］．International Journal for Numerical and Analytical Methods in Geomechanics，2005，29（11）：1109－1126．

［11］ 丁光亚，蔡袁强，徐长节．半空间饱和土中圆形壳结构对平面P波的散射［J］．工程力学，2008，25（12）：35－41．

［12］ Baron M L，Parness R．Displacement and velocities produced by the diffraction of a pressure wave by a cylindrical cavity in an elastic medium［J］．Journal of Applied Mechanics，1962，29（2）：385－395．

［13］ Yoshihara T，Robinson A，Merrit J L．Interaction of plane elastic waves with an elastic cylindrical shell［A］∥structure research series No261［G］．Urbana：University of Illinois，1963．

［14］ Paul S，Robinson A，Ali－Akbarian M．Interaction of plane elastic wave with a thick cylindrical shell，U．S．［R］．New Mexico：Air Force Weapons Laboratory，Kirtland Air Force Base，1965．

［15］ Garnet H，Jacques C P．Transient response of a circular cylinder of arbitrary thickness in an elastic medium［J］．AMSE J．Appl．Mech．，1966，33（3）：521－532．

［16］ Peralta L，Carrier G F，Mow C C．An approximate proce dure for the solution of a class of transient－wave diffraction problems［J］．AMSE J．Appl．Mech．，1966，33（1）：168－183．

［17］ Kobayashi S，Nishimura N．Transient stress analysis of tunnel and caverns of arbitrary shape due to traveling wave in"Development in BEM Ⅱ"ed Banerjee PK and Shaw RD［M］．London：Appl．Sci．，1981：177－210．

［18］ Pao Y H，Ku G C，Ziegler F．Application of the theory of generalized rays to diffraction of transient waves by a cylinder［J］．Wave Motion，1983，5：385－398．

［19］ Biot M A．Mechanics of Deformation and Acoustic Propagation in Porous Media［J］．Journal of Applied Physics，1962，33（4）：1482－1498．

［20］ Durbin F．Numerical inversion of Laplace transforms：an efficient improvement to Dubner and Abate＇s method［J］．The Computer Journal，1974，17（4）：371－376．

附录