离散元_边界元动力耦合模型

　2001年1月

水 利 学 报SHU IL I XU EBAO 第1期

收稿日期:1999212227

基金项目:“九五”国家攻关项目.

作者简介:金峰(1966-),男,贵州遵义人,教授、博士生导师,主要从事结构动力分析研究.

文章编号:055929350(2001)0120023205离散元2边界元动力耦合模型

金　峰1,贾伟伟1,王光纶1

(11清华大学水利水电工程系,北京　100084)

摘　要:本文提出了一种二维变形体离散元与时域边界元的耦合模型,这一模型可以将非连续体的模拟与无限域的模拟统一在一个模型中,可用于在地震波动输入条件下,考虑辐射阻尼的岩体边坡或地下结构等的动力稳定和变形分析,拓宽了离散元动力分析的领域.算例分析表明本耦合分析模型具有较高的精度.

关键词:离散元;边界元;耦合

中图分类号:O344 文献标识码:A

离散元法是一种模拟离散介质的计算方法,自Cundall 在70年代提出以来,在岩石力学、土力学、结构分析等领域的数值模拟中得到广泛应用,是一种新兴的非连续体分析方法[1,2,3].在动力分析中,许多分析均表明辐射阻尼对分析结果影响很大,应该充分重视[4],动力边界元方法由于在其基本解中包含了无限远处的辐射条件,在处理辐射阻尼的影响时十分方便,是分析地下结构动力响应的一种有力的工具[5].动力边界元法可以分为频域边界元法[5]与时域边界元法[6,7],前者主要适用于线性问题.是发展较早、相对成熟的动力边界元方法.时域边界元方法直接在时域内求解,适合解决非线性问题,根据求解的问题具有非线性的特点,本文选择动力时域边界元方法与可变形体离散元耦合,提出了一个二维时域动力边界元———可变形体离散元耦合模型,充分发挥离散元与边界元的优点

,将非连续体的模拟与无限介质辐射阻尼的模拟统一到一个模型中,为地下结构和岩质边坡的抗震稳定分析提供了全新的手段.

1　离散元2边界元耦合模型

图1　离散体系示意

111　二维可变形体离散元原理　可变形体离散元法将模拟的

计算区域看作是若干可变形块体的组合,见图1.这些块体可

以任意平移、旋转,块体之间的相互作用力,用法向和切向弹

簧表示,称为接触,接触力的大小由块体的相对位置决定.每

个块体又划分为三角形差分网格以模拟变形,每个差分三角形

的顶点称为节点,其应变假定为常数,可以由节点的位置确

定.因为一旦节点位置确定,便可以得到差分三角形乃至整个

块体的变形和应力,从而得到所有的响应历程,因此,所有的

计算将围绕节点进行,边界条件也可以通过给定边界节点位移

或节点力来实现.离散元计算采用显式步进的方法,首先将计算的过程分为若干等长的时步,静力分析作为动力分—

32—

析的特例,可以采用临界阻尼以增加收敛速度,此时静力分析的时步不具有真实的物理意义,可看作迭代步.在每一时步内,对所有块体的所有节点分别进行循环计算,每一节点循环计算的主要步骤为:(1)根据上一循环的结果或边界条件,确定本时步初节点的位置和速度等运动量.(2)根据本节点与相邻节点的位置可确定差分网格的应变,根据本构关系求出应力,再积分得到本节点所受弹性力.若为块体边界上的节点,可根据相邻块体的位置确定接触力,从而得到本节点在本时步所受合力的F ′.以上具体公式和算法可参见文献[1,2,3].(3)由牛顿第二定律或动力平衡方程

m ¨u +αm u =F ′+m g

(1)式中:F ′是接触力F c 、弹性力F e 与刚度阻尼力F d =β(F c +F e )之和;am u 代表质量阻尼力;α、β是Rayleigh 阻尼系数;g 代表重力加速度,最后一项表示重力的影响.由式(1)通过中心差分,可得到本时步(第i 时步)的速度

u i = u i -1(1-015αΔt )+F ′m

+g Δt (1+015αΔt )(2)

进而可以得到本时步末本节点的位置.返回第(1)步,继续下一轮循环计算,最终求得所有时刻所有节点的解.

由于任意块体间在运动过程中都有发生各种接触的可能,在计算中需要对所有的块体之间进行接触判断并计算接触力,这是一项十分耗时的工作,也是离散元法的关键技术之一,为减少计算工作量,针对岩体中构造面的特点,已有一些较好的算法,如Cundall 提出的域算法,充分利用生成离散块体时的信息,将接触检索局限于初始顶点构成的域中,特别适合于由构造面切割形成的离散块体系统,由于篇幅所限,在此不能详细介绍,请参阅有关文献[1]和[2].

112　二维时域边界元基本原理　二维全平面时域动力边界元方程[6,7,8]

c αβu β(S ,t )=∫t 0∫Γu 3αβp β(Q ,τ)-p 3αβu β(Q ,τ)

d Γd τ+u I β(3)

式中:u 3αβ,p 3αβ

为二维全平面时域动力基本解,其具体的表达式可参见文献[6]和[8];u ,p 分别表示位移与面力;u I 表示从无限远处入射的位移场;S ,Q 分别表示源点及场点.α,β=1,2,同时对时间、空间进行离散,可以得到时域边界元方程可写作[6,7,8]:

[H ]L {u}0=[G ]0{p}L +{B }

L (4)式中:{u}L 、{p}L 分别是L 时刻的位移与面力分量,

{B }L =

∑L -1

l =1[G ]L -l {p}l -[ H ]L -l {u}l +{u I }L (5)

式中:[H ]l 、[G ]l 、(l =1,L )共2L 个矩阵均为系数矩阵,他们可由各边界单元的子矩阵集成而得,具体表达式可参见文献[6]和[8].

由于{u}l 、{p}l (l =1,2,…,L -1)是1,2

,…,L -1时刻的量,在求解L 时刻时{B }L 为已知量,

而在{u}L ,{p}L 中有N (总的自由度数)个值已知,另N 个值未知,因此可解此方程组,可求得L 时刻的N 个未知量,依此类推,即可求得全部时程的历程反应.

图2　离散元-边界元耦合模型

113　耦合模型的实现　为了实现时域边界元和离散元的耦合,必须解决以下问题:(1)离散元模拟的离散体与边界元分析的连续体实现耦合;(2)求解方式的统一是实现耦合的关键;(3)离散元计算时步与边界元计算时步的协调及离散元节点与边界元单元

节点的对应可以进一步提高耦合计算的效率.下面,分别

介绍解决的方法.

(1)首先,将计算区域划分为两个区域,见图2,即

主要考虑无限域辐射阻尼和地震输入的连续体区域,称为

边界元域,用时域边界元(B EM )模拟,其模拟的重点是

无限域,地表水平自由边界在离散一段距离后截断,经过

—42—