离散元_无限元非线性耦合法

离散元2无限元非线性耦合法

鲍　鹏1,2　姜忻良1

(1天津大学建工学院,天津300072;2河南大学土木建筑学院,河南开封475001)

摘要:根据变形体动力学方程及模态分解原理,建立了可变形体离散元模型.依照离散元法原理,选用边2边接触模型,推导出相应的理式并实现在边界上与无限元的耦合.采用近场土体由离散元模拟,远场土体由无限元模拟,进行土2地下结构相互作用分析.由静力问题的收敛性,验证了计算程序和计算参数选取的正确性,给出了地下结构在地震波作用下的位移反应.

关　键　词:土2地下结构相互作用;离散元;无限元;非线性;耦合法;可变形体中图分类号:TU431　　文献标识码:A　　文章编号:167124512(2004)1120091203

Nonlinearcouplingmethodfordiscreteelement2infiniteelementBaoPeng　JiangXinliangAbstract:AdiscreteelementmodelforthedeformableforthedeformablebodieswasbuiltaccordingtothedynamicequationandthemodeldecompositionBasedonthefundamentalprincipleofdiscreteelementmethod,thecontactmodelwasselectedandthecorrespondingtheoreticalformuladeduced.Inthebound2ary,thediscreteelement2infiniteelementcouplingwereachieved.Theanalysisofthesoil2undergroundstructureinteractioneffectswascarriedoutwhenthenearbysoilwasdescribedasdiscreteelementsandthesoilbeyondasinfiniteones.Thevalidofthecalculatingprogramandtheselectioncalculatingparametersweretested.Thedynamicreactionofundergroundstructurewasobtainedundertheearthquakewaves.Keywords:soil2undergroundstructureinteractioneffect;discreteelement;infiniteelement;non2linear;

couplingmethod;deformablebodies

BaoPeng　Assoc.Prof.;SchoolofCivilEng.,TianjinUniversity,Tianjin300072,China.

　　随着地下建筑、地铁等许多重要复杂结构的

出现,地震作用下的土2地下结构动力相互作用已成为人们关注的问题[1～5].目前土2结构相互作用的研究数值计算方法很多,但尚不能全面反映地震荷载作用下真实的界面特性[2],即在强震作用下结构与土系统会出现不同程度的非线性现象.为此本文提出用于土2结构相互作用非线性问题的离散元2无限元耦合法.

离散元的基本思想是将结构体处理成单元,使每一个单元满足平衡方程,单元间的力则通过作用力和反作用力来传递[1].离散元具有独特的处理非连续介质的能力,能很好地处理土与结构交界面的性能;而无限元适合于地基土无限远

收稿日期:2004202208.

的实际情况,避免了边界失真现象[5].因此用离散元处理结构物附近的土体(近场),用无限元处理结构物周围伸向无限的土体(远场),发挥各自的优点去解决土-结构相互作用问题是一种行之有效的途径.

1　可变形体运动控制方程及其单元

模型

　　根据力学基本公理,若将非惯性坐标系的原点选择在物质的质心处,同时选取在非惯性坐标系中观察测量到的物体角动量为零的坐标轴.可得物体运动控制方程F=m¨r0;　L=ω󰂻・I+ω×(ω・I)+ω・󰂻I;

作者简介:鲍　鹏(192),男,副教授;天津,天津大学建工学院(300072).

E2mail:baop@eyou.com

基金项目:河南省骨干教师资助计划项目(G2002025).

　　　　　　　　　　　　　　华　中　科　技　大　学　学　报(自然科学版)　　　　　　　　　　　第32卷92

　󰂻r・RdΩ+󰂻r・tdΓ=

Δ

-

∫∫

1ω・󰂻ω+󰂻I・K+2

T・󰂻rdΩ,∫

(1)

触边两部分的接触力.

根据不抗拉条件,若σn<0,则令σn=0,dn=0;σs=0,ds=0.若切向应力超过摩擦力,即σσ|σs|>c+μn,则令σs=c+μn;ds=0,其中μ为摩擦系数,μ=tan<;c为粘聚力系数.后一等式成立是因为摩擦滑动本身就在消耗块体的动能,因此无需再加阻尼力.2.2　力和位移的计算循环

根据块体的几何形状及其与邻近块体的关系,可以计算出作用在某一特定块体上的一组力,进而得出它们的合力和合力矩.根据牛顿第二运动定律确定块体质心的加速度和角加速度,进行差分可确定在时步Δt内的速度和角速度以及位移和转动量.在用动态松驰法时,计算循环是以时步Δt向前差分进行.由于时步选取得非常小,每个单元在一个时步内只能以很小的位移与其相邻接的单元作用,而与较远的单元无关.另外由于是可变形体离散元,单元之间的相对位移增量还和单元变形量有关.

式中,K为非惯性系下的动能;L为对点O的合力矩;I为惯性矢量;t为边界面力;R为体力;T为应力并矢量;ω为惯性坐标系下的角速度.

物体的变形由式(1)确定.通过把变形写成固有振动模态的形式,运用加权余量法、拉格朗日法、哈米顿原理或虚功原理,可对模态方程解耦.

2

得出每个模态的控制方程.即mξjj+mjnjξj=

(t・<)dΓ+(R・<)dΩ,式中,m=∫∫ρξdΩ;n为振动模态的角频率(特征值);ξ为<∫

0

j

0jj

jjjj

广义坐标;假设单元应变方程可写成与刚体模态无关的

kkkk)=σ形式:mk(εA-σI,式中,m为相对于应变

k

模态k的单元广义(有效)质量;ε为应变模态k

k

的内应变率;σA为相对于应变模态的单元施加应k力;σI为相对于应变模态的单元内应力.施加应力σij由作用在块体周边上的接触力表示,即ACCσ表示对所有接触点求和;Vij=∑FiXj/V,∑A

3　离散元与无限元的耦合3.1　无限元单元模型

为块体的体积(对于二维问题是面积S);FCi为作用于接触点C的集中力,XCj为接触点C相对于块体形心的坐标.单元内应力σIij由弹性关系获得:

IσIij=2Gεεδij+λkkij,式中,σij为单元内应力张量;

无限元基本原理是通过拉格朗日插值函数和

衰减函数的乘积来构造形函数.对于位移无限元,一定要选择能反映位移衰减特征的衰减函数,以反映在介质中由近场至远场的位移分布规律,及满足在无穷远处位移为零的条件.本文运用六节点超参无限元,易于和离散元的位移边界条件实

η坐标之间的变换现耦合.图1中,x2y坐标和ξ2

εij为单元应变张量;εkk为单元主应变张量;G为

)・材料剪切模量;λ为拉梅常数,λ=E/[(1+ν

(1-2ν)];δij为Kronekerdelta函数.

2　离散元法的基本方程

2.1　力和位移的关系

对于岩土工程中的许多问题,采用弹塑性无拉力边-边接触模型较为合理.由于接触边作用的是应力,确定应力在接触边的分布情况,最为简便的方法是把该应力近似地表示成两个均匀分布的应力.引入刚度系数Kn和Ks,于是法向应力和切向应力为

σΔun;　σΔus.n=Kns=Ks

考虑接触点的阻尼力,并假定阻尼力与位移增量成正比,则有

Δun;　ds=-βKsΔus,dn=-βKn

式中β2块体的刚度阻尼比例系数.由此可得到接

关系是:

4

图1　无限元单元

4

X=

′

i=1

∑Nx

′

i

i

;　Y=

′

i=1

y∑N′

ii

,

′

式中,N′1=(1/2)(ξ-1)(η-1);N2=-(1/2)・′(ξ-1)(η+1);N′3=(1/2)ξ(η+1);N4=-(1/

2)ξ(η-1).

单元位移模式取为

2

2

U=

′

i=1

∑Nu;　V

″

ii

′

=

i=1

v∑N″

i

i

,

);N″式中,N″1=-(1/2)(η-1)f(ξ2=(1/2)・

第11期　　　　　　　　　　　　鲍　鹏等:离散元2无限元非线性耦合法　　　　　　　　　　　　　　93

(η+1)f(ξ);f(ξ)=e-ξ为单元的位移衰减函数.

注意到位移模式中只选用了节点1和节点2两个位移分量,实际上通过衰减函数,位移模式仍反映了中节点(节点3,4)以及其余节点的位移分布.将形函数写成矩阵形式,可组成单元刚度矩阵和单元质量矩阵,按对号入座的办法,形成整体刚度矩阵和质量矩阵.3.2　与离散元的耦合

图4　人工波(中震)加速度时程

在离散元区域(近场)和无限元区域(远场)的交界处,在t时刻离散元对边界的作用,以力边界条件的形式施加到无限元1,2节点上,由无限元整体平衡方程解出1,2节点的位移,作为t+Δt时刻离散元的位移边界条件.

4　地下结构箱体的地震反应分析及

结论

　　根据上述原理,编制了相应的计算程序.首先计算了某地铁车站箱体结构(见图2)在地面竖向

图5　人工波(小震)作用下的箱体水平位移图6　人工波(中震)作用下的箱体水平位移

图2　地下结构计算简图

荷载作用下的平面静力问题.验证计算程序和计算参数选取的正确性,然后对某地铁车站箱体结构进行地震反应分析.箱体结构取为的刚性离散单元,土体近场取为可变形离散单元,远场取为无限元.土体材料参数:弹性模量E=1×108N/m2,泊松比ν=0.3,摩擦系数μ=0.3,粘聚力系数c=15kPa,重度γ=18kN/m3.离散元法计算参数:Kn=2.0×108N/m2,Ks=7.5×107N/m2,α=1.0×102,β=3.0×10-5,Δt=5.0×10-3s.地震波从底部边界输入.求得人工

(小震)作用下的箱体水平方向位移反波(中震)、

应(其结果见图3～6).

　　结论:a.将基于可变形体离散元模型的离散元2无限元非线性耦合法用于地下结构2土体系相互作用分析,结果证明是可行的.b.在近场采用可变形体模型的离散元,既可描述土体的变形性质,又可通过调整块体接触边的弹塑性模型,模拟块体间线性、非线性以及连续、非连续的接触关系.而在远场采用无限元,可较好地模拟无限边界的性质.c.离散元2无限元非线性耦合法为解决桩2土相互作用以及其他存在不连续界面的岩土工程问题提供了新的途径.

参

考

文

献

[1]王泳嘉,邢纪波.离散单元法及其在岩土力学中的应

用.沈阳:东北工学院出版社,1991.

[2]雷晓燕.岩土工程数值计算.北京:中国铁道出版社,

1999.

[3]瓦尔夫JP.土2结构动力相互作用.吴世明译.北京:

地震出版社,19.

[4]鲍　鹏,姜忻良.适用于可变形体的离散元法.河南大

学学报(自然科学版),2004(1):100～103

图3　人工波(小震)加速度时程

[5]赵崇斌.用无穷元模拟半无限平面弹性地基.清华大

学学报,1986(2):51～