波动问题中的三维时域粘弹性人工边界

1、第22卷第6期 工 程 力 学 V ol.22 No.6 2005年 12 月ENGINEERING MECHANICSDec. 2005收稿日期:2003-12-03;修改日期:2004-05-24基金项目:国家973项目(2002CB412706;国家自然科学基金(50078028;北京市自然科学基金重点项目(8011002资助杜修力(1963,男,四川广安人,教授,博士,博士生导师,从事结构工程和防灾减灾工程研究; 杜义欣(1979,男,河北保定人,博士生,从事结构工程和人防工程研究.文章编号:1000-4750(200506-0046-06波动问题中的三维时域粘弹性人工边界*刘晶波1,

2、王振宇1,杜修力2,杜义欣1(1. 清华大学土木工程系,北京 100084;2. 北京工业大学,北京 100022摘 要:应用弹性波动理论,发展了实现波动直接模拟的三维时域粘弹性人工边界。首先基于三维波动方程推导了三维粘弹性人工边界的法向与切向边界方程,然后研究了时域粘弹性人工边界的数值模拟技术,最后通过典型的波动问题算例证明了给出的三维人工边界具有较高精度,可以方便地应用于三维波动问题的模拟分析。 关键词:波动;粘弹性人工边界;三维介质;时域;无限域 中图分类号:O347.4, P315.3 文献标识码:ATHREE-DIMENSIONAL VISCO-ELASTIC ARTIFICIAL

3、BOUNDARIESIN TIME DOMAIN FOR WA VE MOTION PROBLEMS*LIU Jing-bo 1 , WANG Zhen-yu 1 , DU Xiu-li 2 , DU Yi-xin 1(1. Department of Civil Engineering, Tsinghua University, Beijing 100084, China; 2. Beijing University of Technology, Beijing 100022, ChinaAbstract: Three-dimensional visco-elastic artificial

4、 boundaries in time domain are developed by means of theory of elastic wave motion in this paper. Three-dimensional visco-elastic artificial boundary equations along normal and tangent directions are derived based on three-dimensional wave motion equations. Numerical simulation techniques of visco-e

5、lastic artificial boundaries in time domain are studied. Numerical examples of classic wave motion problems demonstrate that high precision is achieved by use of three-dimensional visco-elastic boundaries, and that the boundaries can be used in analysis of three-dimensional wave motion problems easi

6、ly.Key words: wave motion; visco-elastic artificial boundary; three-dimensional media; time domain; infinitedomain二十世纪七十年代以来,与土-结构动力相互作用问题相关的动力基础振动、地震波散射等近场波动问题的数值模拟方法得到了很大的发展13。在无限域介质的瞬态波动问题模拟中,一种广泛采用的数值模拟技术是引入虚拟的人工边界从无限介质中切取出有限尺寸的近场计算区;然后对计算区采用有限元(或有限差分技术完成运动微分方程和物理边界条件的时空离散化,使对波动的模拟归结为代数和算术运算,从而实

7、现对真实波动的直接模拟。显然,构造合适的人工边界是实现上述方法的关键,人工边界的本质在于允许来自广义结构的外行散射波穿过人工边界进入无限域。以此为出发点,近三十年来国内外对人工边界进行了广泛而深入的研究,并发表了一系列综述35。在人工边界波动问题中的三维时域粘弹性人工边界 47两个主要方向:全局与局部人工边界的研究中,后者由于其时空解耦特性和广泛适用性,在无限域波动模拟问题中受到重视。在局部人工边界的发展历程中,较早出现的粘性边界因其概念清楚、应用方便,在相当长的时间内得到了广泛应用6,7。但粘性边界仅考虑了对散射波能量的吸收,从物理概念上理解,施加粘性边界后的力学模型为悬浮在空中的脱离体,在

8、低频力作用下可能发生整体漂移;此外,粘性边界是基于一维波动理论提出的,简单地将其推广到多维情况将导致相当误差8。为了克服上述缺点,Deeks 9、刘晶波8,10基于柱面波动方程建立了二维粘弹性人工边界,并指出该边界相对于粘性边界的优点是能模拟人工边界外半无限介质的弹性恢复性能,具有良好的高频和低频稳定性。实际的波动问题在很多情况下是三维的,一方面我们缺乏三维形式的粘弹性人工边界,另一方面简单地将低维边界应用于多维情况将不可避免地导致理论的合理性和应用的可靠性问题。因此,发展实用的三维粘弹性人工边界及其数值模拟技术具有重要意义,本文将基于近场弹性波动理论推导三维粘弹性人工边界方程,并讨论相应的数

9、值模拟技术,最后通过经典波动问题数值模拟结果与解析解的对比,证明三维粘弹性人工边界的可靠性与实用性。1 三维粘弹性人工边界法向条件1.1 法向边界方程球坐标系中球面膨胀波(P 波的波动方程为22222(1t R c R R p =(1 其中为位移势函数,c p 为介质的P 波波速,R 为径向坐标。方程(1的通解可表示为 (1(1,(t c R g Rt c R f R t R P P += (2式中f (和g (为任意函数,分别表示外行扩散波和内行会聚波。考虑外行扩散波,垂直于波阵面的位移可写为(1(12t c R f Rt c R f R R u p p = (3而法向应力由下式计算R uR

10、u 22(+= (4式中、为拉梅常数。由式(3可得(2(2(132t c R f R t c R f R t c R f R R u p p p += (5(1(132t c R f Rt c R f R R u p p = (6 式(5、(6代入式(4可得波阵面上用函数f (表示的法向应力为(4(4(232t c R f R t c R f R t c R f R p p p += (7 为建立法向应力与位移u 之间的关系式,引入以下方程(2t c R f R c t c R f R c t uup p p p +=& (8 (22222t c R f R c t c R f Rc t u

11、up p p p =& (9(4(4(232t c R f R c t c R f R c t c R f R c t p p p p p p+=(10 由式(3、式(7至(10可得在波阵面上法向应力和位移满足+=+uG R u c R u R G c R pp &442 (11 方程(11即为三维法向人工边界方程,该方程给出了波阵面上法向应力与位移的关系式,式中=G 为介质剪切模量,为介质密度。注意推导方程(11时利用了关系式:22p c =+ (12 式(12中、为拉梅常数,为介质密度。当时间趋于无穷大时,三维法向人工边界方程(11退化为u RG4= (13 式(13即为静力边界方程11,

12、该方程表示弹性体中球形空腔受均匀压力时,空腔附近位移与应力的关系。1.2 法向边界数值模拟技术为建立人工边界,将无限连续介质截断,在截断处,即人工边界上施加连续的弹簧-阻尼器-集中质量系统,如图1所示。图1给出的物理系统的运动方程为48 工 程 力 学=+(M R R u uC Ku & (14 0(=+R M M u u C uM & (15 式(14与(15中R u 、M u 分别表示人工边界节点与集中质量沿荷载作用方向位移。 图1 法向边界上施加的物理系统Fig.1 Physical system imposed on vertical boundary由式(14可得(1+=R R M

13、uC Ku C u & (16 (1&+=R R M u C u K Cu(17 式(16、(17代入式(15,可得到关于施加物理系统的人工边界节点应力与位移满足的微分方程:+=+22t u K M t u C M u K t C M R R R (18 将式(11与(18进行对比,可以发现当物理元件参数为RGK 4=;p c C =;R M = (19 时,人工边界上力和位移的条件与原连续介质的完全相同,即只须采用相应参数的弹簧、阻尼器和集中质量单元即可实现三维粘弹性人工边界。应该说明,在实际的有限元实现过程中,连续分布的物理常数可采用一定规则进行离散化,例如,法向边界参数采用式(19的量值

14、与对应有限元节点所代表网格面积的乘积。2 三维粘弹性人工边界切向条件2.1 切向边界方程球坐标系中球面剪切波动(S 波位移的近似解为12(1(1,(t c R g Rt c R f R t R u s s += (20其中s c 为介质剪切波(S 波波速,等式右边第一和第二项分别表示外行扩散波和内行会聚波。考虑扩散的球面剪切波,波阵面的切向位移为(1,(t c R f Rt R u s = (21由式(21所确定的剪应变和剪应力分别为(1(2,(2t c R f R t c R f RR uR u t R s s +=(22+=(1(2,(2t c R f R t c R f R G G t

15、R s s (23在坐标R 处质点的运动速度为(,(t c R f R c t t R u u s s =& (24由式(22、(23、(24可得在波阵面上: ,(,(2,(t R uc t R u RGt R s &= (25 式(25即为三维切向人工边界方程,该方程给出了波阵面上切向应力与位移的关系式,式中=G 为介质剪切模量,为介质密度。注意推导方程(25时利用了关系式2s c G =。2.2 切向边界数值模拟技术对于如式(25所示的三维切向边界方程,容易证明该方程等价于并联的弹簧-阻尼器系统,对应物理元件的参数为RGK 2=;s c C = (26 同样,在实际的有限元实现过程中,切向

16、边界物理元件参数的实际取值等于式(26给出的量值与对应有限元节点所代表网格面积的乘积。3 粘弹性人工边界为了模拟实际连续介质条件,方程(19、(26分别给出了人工边界法向和切向应施加的物理元件形式,其中法向边界物理元件中的质量M 与阻尼器相连(如图1所示,显然这是一个不稳定的系统。为了克服实际处理和计算时可能引起的不便,我们将质量M 忽略,并将与质量M 相连的阻尼器的一端固定,从而形成粘性阻尼器+弹簧的人工边界,与切向的人工边界一起,统称为粘(阻尼器弹(弹簧性人工边界(Viscous-Spring Boundary,关于这一处理的有效性和精度将通过算例证明。应当注意,粘弹性人工边界所模拟的是人

17、工边界上的应力条件,因此,这是一种连续分布的人工边界条件。当采用有限元法或其它离散化方法将人工边界所包围的计算区离散化时,人工边界面也将随之离散化,此时,可以采用有限元法的形函数将连续分布的人工波动问题中的三维时域粘弹性人工边界 49 边界物理元件化为耦联的人工边界,可称为一致粘弹性人工边界;也可以简单地采用集中处理方法,形成解耦的人工边界,称为集中粘弹性人工边界。一般情况下,可以采用相对简单的集中粘弹性人工边界条件,其具体实施方法如图2所示,图中坐标X 、Y 沿人工边界的切向,Z 为法向,图中粘弹性人工边界节点上物理元件的参数为=Ii ip Ii i Ii i s Ii i A c C A

18、R G K A c C C A R G K K 1313121121,4,2 (27其中,A i 为人工边界上节点所代表的面积,对图2所示情况I =4。 图2 三维粘弹性人工边界示意图 Fig.2 Three dimensional viscous-spring boundary4 数值算例为了验证本文提出三维粘弹性人工边界的精度和可靠性,考察两个具有解析解的半无限介质波动问题,通过不同计算条件下结果的对比,分析三维粘弹性人工边界的精度和适用性。算例中采用的人工边界为集中粘弹性边界。 4.1 Lamb 问题算例考察Lamb 表面源问题:一均匀、各向同性半空间,在自由面上受加载函数的集中垂直荷载

19、作用,如图3所示。图3中加载函数的表达式为 1(43(421(641(4(16(44444+=G G G G G (28Tt H G =,(34 其中T 为持续时间,(H 为Heaviside 阶梯函数。图3 Lamb 问题算例示意图 Fig.3 Lambs problem算例中对应于加载函数的解析解由Lamb 问题基本解13,14积分求得。而对应的数值解则分别采用不同的人工边界获得,这些边界为:三维粘弹性边界(3D viscous-spring boundary、粘性边界(viscous boundary、固定边界(fixed boundary。此外为了与带有悬挂质量块的粘弹性人工边界(ma

20、ss-viscous-spring boundary进行对比,也同时计算了有悬挂质量块的人工边界工况,算例中称为三维带质量块粘弹性边界。从理论上讲,人工边界的形状采用图3所示的半球形将提供更好的模拟精度。但考虑到实际有限元计算中模型选取和对人工边界处理的简便性,算例均取平直的人工边界,即取人工边界方向与直角坐标轴垂直,空间上形成矩形计算区。此外,人工边界物理元件参数中的距离R 并不取坐标原点到人工边界各个节点处的距离,而是统一取为坐标原点到相应人工边界面的最短距离,即每一个平直人工边界上物理元件的参数相同。数值计算对象范围选取5.0=b b x x ,b z =1.0;有限元网格1.0=z y

21、 x ;介质剪切波速C s =4;密度=1.0;泊松比=0.25。观测量为自由面上点A (r =0.2与B (r =0.4的垂直位移反应。图4和图5分别为观测点A 和B 的位移反应时程曲线,可以清楚看到,粘性边界导致了明显的漂移,而本文给出的三维粘弹性人工边界具有较高的精度,其模拟效果明显优于粘性边界与固定边界。从图中还可以发现,采用三维粘弹性边界和带悬挂质量的三维粘弹性边界给出的位移时程曲线符合较好,这说明对于本问题,三维粘弹性边界与带悬挂质量的三维粘弹性边界的模拟精度相同。50工 程 力 学 -0.005时间/s 图4观测点A(r =0.2位移时程比较Fig.4 Displacement

22、time-history of observation point A (r =0.2 -0.0050.0000.0250.030位移/m时间/s图5 观测点B(r =0.4位移时程比较Fig.5 Displacement time-history of observation point B (r =0.44.2 内源问题算例考察在均匀、各向同性线弹性无限介质中集中力源产生的波动。采用直角坐标系,集中力作用于坐标原点,如图6所示。 图6 内源问题示意图 Fig.6 Inner source problem数值计算对象范围选取5.0=b b b z x x |,观测点的坐标为x =y =z =

23、0.1,其它计算参数与前述表面源Lamb 问题相同。图7给出了采用不同人工边界条件下观测点z 向位移时程与解析解的对比情况。由图7可以看出,同Lamb 问题算例一样,粘性边界导致了明显的漂移,而此时三维粘弹性边界仍然具有良好精度。同时也可以看到,当参数R 不是取坐标原点至人工边界的最短距离0.5,而是统一取为R=0.707时,粘弹性人工边界的模拟精度进一步提高。位移/m时间/s图7 内源问题加载函数作用下观测点z 向位移时程 Fig.7 Displacement time-history in z direction of observationpoint under -function in

24、 inner source problem5 结语建立合理有效、有足够模拟精度且简单易用的三维人工边界对于无限域介质中的波动模拟具有重要意义。本文基于弹性波动理论推导了三维时域粘弹性人工边界,研究了相应的数值模拟技术。数值算例表明,三维时域粘弹性人工边界具有较高精度;等效离散后的人工边界条件方程与有限元方法相结合,可以方便地求解无限域介质的瞬态波动问题。另外,数值算例还表明三维时域粘弹性人工边界的精度与带悬挂质量的粘弹性人工边界相当,但三维时域粘弹性人工边界的应用更方便。时域粘弹性人工边界中,弹簧刚度对人工边界的模拟精度有重要影响,它的选取与波源到人工边界点的距离有关,同时与波动传播方向与人工

25、边界法向之间的夹角有关。本文算例仅简单地取波源到人工边界面的距离代替波源到人工边界上各节点的距离,并假设散射波的入射方向与人工边界法向之间的夹角为零,因此,弹簧刚度的取值总体偏大。波动问题中的三维时域粘弹性人工边界 51 另外，本文推导粘弹性人工边界时，应用的是弹性 全空间理论，全空间中波动的位移场和应力场与半 空间中的有所不同，基于全空间理论得到的粘弹性 边界用于处理实际的半空间问题时，为提高模拟精 度，相关参数需予以调整和修正。粘弹性人工边界 弹簧刚度的优化选取方案及其鲁棒性是两个值得 进一步研究的问题。 参考文献： 1 廖振鹏. 工程波动理论导论(第二版M. 北京: 科学 出版社, 20

26、02. Liao Zhenpeng. Introduction to wave motion theories in engineering M. Beijing: Science Press, 2002. (in Chinese 刘晶波, 王振宇. 考虑土-结构相互作用大型动力机器 基 础 三 维 有 限 元 分 析 J. 工 程 力 学 , 2002, 19(3: 3438. Liu Jingbo, Wang Zhenyu. 3D finite element analysis of large dynamic machine foundation considering soil-str

27、ucture interaction J. Engineering Mechanics, 2002, 19(3: 3438. (in Chinese 廖振鹏. 近场波动数值模拟J. 力学进展, 1997, 27(2: 193212. Liao Zhenpeng. Numerical simulation of near-field wave motion J. Advances in Mechanics, 1997, 27(2: 193212. (in Chinese Kausel E. Local transmitting boundaries J. Journal of Engineeri

28、ng Mechanics, ASCE, 1988,114(6:10111027. Wolf P J. A comparison of time-domain transmitting boundaries J. Earthquake Engineering and Structure Dynamics, 1986, 14: 655673. Lysmer J, Kulemeyer R L. Finite dynamic model for 9 8 7 10 2 11 12 3 13 14 4 5 infinite media J. Journal of Engineering Mechanics

29、. ASCE, 1969, 95: 759877. 中华人民共和国(GB 50267-97. 核电厂抗震设计规范 S. 北京: 中国计划出版社, 1998. The Peoples Republic of China (GB 50267-97. Seismic Design Code of Nuclear Power Plants S. Beijing: China Planning Press, 1998. (in Chinese Liu Jingbo, Lu Yandong. A direct method for analysis of dynamic soil-structure in

30、teraction based on interface idea M. In: Zhang Chuhan, Wolf J P, edited, Dynamic Soil-Structure Interaction. Academic Press, 1997. 258273. Deeks A J, Randolph M F. Axisymmetric time-domain transmitting boundaries J. Journal of Engineering Mechanics, ASCE, 1994, 120(1: 2542. 刘晶波, 吕彦东. 结构-地基动力相互作用问题分析

31、的 一种直接方法J. 土木工程学报, 1998, 31(3: 5564. Liu Jingbo, Lu Yandong. A direct method for analysis of dynamic soil-structure interaction J. China Civil Engineering Journal, 1998, 31(3: 5564. (in Chinese 徐芝伦. 弹性力学M. 北京: 高等教育出版社, 1990. Xu Zhilun. Elasticity mechanics M. Beijing: Higher Education Press, 1990. (

32、in Chinese Aki K, P. G. Richards. Quantitative Seismology Theory and Methods, Vol. I M. W. H. Freeman and Company, 1980. Pekeris C L. Proceedings of the national academy of sciences C. USA, No.41, 1955. 469. 王振宇. 大型结构地基系统动力反应的计算理论及 其应用研究D. 北京: 清华大学, 2002. Wang Zhenyu. Computational theory of dynamic

33、 response of large structure-soil systems and its application D. Beijing: Tsinghua University, 2002. (in Chinese 6 (上接第 45 页 3 M Ozturk, F. Erdogan. Axisymmetric crack problem in bonded materials with a graded interfacial region J. Int. J. Solids Structures 1996, 33: 193219. P Gu, R J Asaro. Cracks

34、in functionally graded materialsJ. Int. J. Solids Structures. 1997, 34: 117. H J Choi. The problem for bonded half-planes containing a crack at an arbitrary angle to the graded interfacial zone J. Int. J. Solids Structures, 2001, 38: 65596588. 周维垣, 肖洪天, 吴劲松. 三维间断位移法及其强奇 异和超奇异积分的处理方法J. 力学学报, 2002, 34(4: 645651. Zhou Weiyuan, Xiao Hongtian, Wu Jingsong. Three-dimensional discontinuous displacement method and the strongly singular and hypersingular integrals J. Acta Mechanica Sinica, 2002, 34: 645651. (in Chinese Z Q Yue, H T Xiao. Generalized K