加筋土挡土墙动力特性分析

1、第17卷第2期2004年4月中国公路学报China Journal of Highway and TransportVol.17No.2Apr.2004文章编号:100127372(2004022*收稿日期:2003208203作者简介:李海深(19642,男,湖南宁乡人,湖南大学副教授,工学博士研究生.E 2mail :xiangtanlhs 加筋土挡土墙动力特性分析李海深1,杨果林1,2,邹银生1(1.湖南大学土木工程学院,湖南长沙410082;2.中南大学土木建筑学院,湖南长沙410075摘要:通过对常用数值计算方法的分析、对比及评价,运用有限元计算方法,建立了加筋土弹塑性本构模型,编制

2、了加筋土挡土墙在动荷载作用下的通用数值计算程序。筋材与其上下表层一定厚度的填土层共同工作,视为筋土单元,成功地解决了筋材与土介质的接触问题,较好地解决了墙面板与筋材连接处两个接触面单元的特殊问题,不仅使计算方便,而且也符合实际情况。关键词:道路工程;加筋土挡土墙;有限元分析;动力特性;弹塑性本构模型中图分类号:U417.115文献标识码:AAnalysis of dynamic character of reinforced earth retaining w allL I Hai 2shen 1,YAN G Guo 2lin 1,2,ZOU Y in 2sheng 1(1.School of

3、 Civil Engineering ,Hunan University ,Changsha 410082,China ;2.School of Civil Engineering and Architecture ,Central S outh University ,Changsha 410075,China Abstract :By analyzing ,comparing and appraising numerical analysis methods in common used ,a com 2mon analytical programme has been worked ou

4、t with elastic 2plastic stress 2strain relationship and finite element method is used to analyze reinforced earth retaining wall under dynamic load.A satisfying so 2lution has been put forward to the contact boundary ,which element of reinforcement 2soil will be seen a total.The especial problem of

5、the joints of contact face element between plate and reinforcement is solved.It not only calculates convenient but also tallies with the fact.K ey w ords :road engineering ;reinforced earth retaining wall ;finite element analysis ;dynamic charac 2ter ;elastic 2plastic stress 2strain relationship0引言求

6、解加筋土挡土结构弹塑性动力问题,实际上可归结为按初始条件和边界条件求解偏微分方程M U +C U t +KU t =R t 的初值边值问题。当只有初始条件而没有边界条件就成为初值问题,反之则为边值问题。区域内的偏微分方程称为基本方程,初始条件是表示初始状态的条件,边界条件是表示边界约束情况的条件。对于工程上提出的问题能采用解析法按照边值条件求解偏微分方程的仅限于极少数情况。所以,一般只能采用近似方法求解。随着计算机的广泛应用,数值解法逐渐地成为解边值问题的一种有效的方法。数值解法分为区域型和边界型两大类。区域型数值解法主要是有限元法(Finite Element Method ,FEM 和有限

7、差分法(Finite Difference Method ,FDM 。边界型数值解法主要是边界元法(Bound 2ary Element Method ,B EM 。采用差分法时,将所考虑的区域织成网格,用差分近似微分,把差分方程变成微分方程。通过数学上的近似,把求解微分方程的问题变换成求解关于结点未知量的代数方程的问题。采用有限元法时,将所考虑的区域分割成有限小区域,称为有限单元。这些单元仅在有限结点上相连接,根据变分原理把微分方程变换成变分方程,它是通过物理上的近似,把求解微分方程的问题变换成求解关于结点未知量的代数方程的问题。边界元法是继有限差分法、有限元法之后发展起来的又一数值计算方法

8、。采用边界元法求解时,根据积分定理,将区域内微分方程变换成边界上的积分方程,然后,将边界分割成有限大小的边界元,把边界积分方程离散成代数方程,把求解微分方程的问题变换成求解关于结点未知量的代数方程的问题。它具有降一维的特性,所占的计算机时小、计算时间省,与有限元相耦合能较好地解决工程实际问题。由于边界元法本身适用于无限域和半无限域,所以这一方法已在各个工程领域获得了广泛应用。边界单元法只有在边界上剖分单元,通过基本解把域内未知量化为边界未知量来求解,这就使自由度数目大大减少,而且由于基本解本身的奇异性特点,使得边界元法在解决奇异问题时精度较高,另外,基本解可以根据实际问题的特点适当选择,以达到

9、最大限度地节约之功效,甚至可以避免直接处理无限边界问题。再者,边界元法的降维作用使得问题简化许多,特别是三维问题。边界元法特别适用于大区域、无限域和断裂、耦合问题1。结合的目的在于利用边界元法处理无限域边界所具有的优点以及有限元法求解非线性问题的灵活性,用以消除有限元法上的“边界效用”及边界元法域内剖分的不便,在减小解题规模方面也有明显的效果,这种方法最先被用于弹塑性分析,以后又发展到用于求解动力问题。文献2中,用有限元边界元耦合法进行动态响应分析,将求解区域分割为二,分别采用有限元法和边界元法,在其交界面上通过迭代法满足界面条件,然后利用WIL SON2法求解动力微分方程,此法兼有有限元法和

10、边界元法的特点,具有效率高、输入数据少、节省机时等优点。文献3中从分区势能和分区混合能原理出发,讨论了弹性、弹塑性力学问题的有限元边界元对称耦合法,它能够用于常规有限元和杂交有限元等不同的有限元模型和边界元的耦合分析。无限域单元或称无限元是有限元中专门模拟无限域边界的特殊单元,可视为另一种耦合方法。无限元的特点是采用一种特殊的形函数及位移插值函数使其能反映在无限元处的边界条件。无限元应用的主要优点在于:有效地解决了有限元分析中的“边界效应”及人工边界的缺点,在动力问题中尤为突出;提高了求解精度及计算效率,对三维问题尤为显著;显著减小了解题规模,为微机应用提供了十分有利的条件。另外,马立明等人4

11、,从运动平衡方程的形式入手,应用变分原理,在由空间和时间张成的四维广义空间上,如同建立常规的空间有限元模型一样,对时间域进行离散,建立时间有限元模型。导出了动力分析的单步及两步时间元法。在对两种时间元法的稳定性进行分析的基础上,构造出相应的无条件稳定格式。有限元法、有限差分法和边界元法这些有力的数值技术建立在连续性假设的基础上,离散元法主要处理物体间具有不连续性问题的数值方法,重点是求解多个物体间的接触和冲击问题。一个物体与另一个物体是通过边界接触联系的,边界接触可以随时间变化,一个单元与另一个单元在任何处接触没有限制,可以是结点与结点的接触,也可以是结点与边界接触,接触单元之间产生的接触力遵

12、循不同的接触原则。综上所述,当今有限元法被公认为是一种用数值方法求解工程中遇到的各种问题的最有效、最通用的方法。有限元法目前已成为求解具有已知边界和初始条件或两者之一的偏微分方程组的一种通用的数值解法。鉴于此,笔者用有限元法求解加筋土挡土结构的动力问题。1有限元动力分析方法有限元系统的运动平衡方程M d+C d+K d=R(1式中:M、C、K分别为质量矩阵、阻尼矩阵和刚度矩阵;R是外荷载向量;d、d和d分别为有限元集合体的位移、速度和加速度向量。式(1是通过考虑在时刻t的静力平衡而推导出来的,即式(1可写成F1(t+F D(t+F E(t=R(t(2式中:F1(t为惯性力,F1(t=M d;F

13、 D(t=C d为阻尼力;F E(t为弹性力,F E(t=K d。它们均与时间t有关,因此,在动力分析中,原则上可认为是考虑与加速度有关的惯性力和与速度有关的阻尼力作用在时刻t的静力平衡。在数学上,式(1是一个二阶线性微分方程组,原则上可用于求解常系数微分方程组的标准过程来求得方程的解。但是,如果矩阵的阶数很高,除非特别利用系数矩阵M、C和K的特殊性质,否则,采用求解一般微分方程组的过程可能要付出很高的代92第2期李海深,等:加筋土挡土墙动力特性分析价。在实际有限元分析中有两种基本的求解方法:直接积分法和振型叠加法。1.1直接积分法在直接积分法中对式(1逐步地进行数值积分,进行数值积分前没有把

14、方程进行另一种形式的变换,实质上,直接积分是基于下面的两种想法,第一个想法是只在相隔t的一些离散的时间区间上而不是试图在任意时刻t满足式(1。即包含有惯性力和阻尼力作用的(静力平衡是在求解区间上的一些离散时刻点上获得的,因此,似乎在静力分析中使用过的所有求解方法,在直接积分法中可有效地使用。第二个想法是假定位移、速度和加速度在每一时间区间t内变化,位移、速度和加速度在每一个时间区间内变化,决定着求解的精度、稳定性和求解过程的费用。假设分别用d0、d0、d0来表示初始时刻(t=0的位移、速度和加速度向量为已知且要求出式(1从t=0到t=T的解,在求解时把时间全程划分为n 个相等的时间区间t(t=

15、T/n,所用的积分格式是在时刻0,t,2t,3t,t,t+t,T上解方程的近似解。由于计算下一个时刻的解的算法要考虑到前面各个时刻的解,因此假定在时刻0,t,2t,3t,t的解为已知,来推导出求解时刻t+t的解的算法。计算时刻t+t的解对于计算自此以后的时刻t上的解是具有代表意义的,这样可建立用来计算所有离散时间点上解的一般算法。1.2振型叠加法直接积分法所需的运算次数直接正比于分析中的时间步数,因此,一般来说,当要求较短的时间(n 个时间步数的响应时,可以预料,使用直接积分法是有效的,但是,如果积分必须对许多时间步进行,则先把平衡方程式(1变换,使之能以较少的费用进行逐步求解就可能会更有效。

16、由于所需的运算次数直接正比刚度矩阵的半带宽m k,因而m k的减小会按比例地降低逐步解题的费用。振型叠加法和直接积分法之间的惟一区别是振型叠加法在时间积分之前进行了基的变换,即从有限元坐标基变换为广义特征向量问题K=2M的特征向量基。从数学观点上看,由于n个特征向量所张成的空间和有限元的结点位移所张成的空间是相同的,所以两种分析所得到的解必定是相同的,因而选择直接积分还是选择振型叠加仅需根据效果来考虑。2隐式一显式瞬变动态分析在非线性瞬变动态应力分析中,可用振型叠加法,但实际上,通常都采用时间步进法,这种直接积分法大体上可分为显式法和隐式法两类5。采用非常通用且易于执行的显式中心差分法,在计算

17、每一时间步时,由于没有必要进行形式上的矩阵分解,因此只需要较小的计算量。可惜这种方法的稳定是有条件的,而且常常只有当时间步取得很小时才能稳定。采用隐式方法,需要对矩阵进行因式分解,然而可以选择一种无条件稳定的隐式算法,其时间步长只取决于精度的要求。本程序中可任意选择下列算法之一:隐式解法、显式解法、隐式显式组合解法。在结构有限元分析时,结构中总有一部分区域由于变形大,而通过减小单元尺寸加密单元网格方法来精确模拟其变形和应力分布,其余区域部分则由于变形和应力梯度相对较小而使用较大尺寸单元和粗疏网格剖分形式。在这种情况下,通常一部分区域作隐式积分计算,而后一部分区域作显式积分计算,这样可充分应用显

18、式、隐式算法的各自优点。B EL YSTCH KO、HU GHES分别基于结点分离法和单元分离法提出了两种不同的显式隐式混合算法。在此基础上CHIAN G通过设计方程组并行求解器来实现显式隐式混合并行的功能。冠哲君等人6基于区域并行法,提出了一种新的较前者具有更高并行粒度的显式隐式混合并行法,给出了其算法执行过程。在物理直观意义上解释了该算法的设计思想,其本质是一种单元弱耦合的混合积分算法,并在网络机群上得以实现。隐式、显式时间积分方法结合起来,这时有限元网格包括两组单元:隐式组和显式组,用上标I和E 分别表示隐式组和显式组。隐式显式法中,为满足方程M a n+1+P I(d n+1,V n+

19、1+P E( d n+1, V n+1=f n+1式中:f n+1=f I n+1+f E n+1;M=M I+M E,并假定M E为对角线矩阵,就需要在每个时间步进行迭代。本程序参考文献5的程序,在这些程序基础上作了补充、修改、改进工作,使之适用于加筋土挡土结构动力分析,兼顾通用性,采用4、8和9结点等参四边形单元处理平面应力、平面应变和轴对称问题,并用总的LA GRAN GE公式来处理几何非线性问题,假定材料特性是弹塑性材料的模型。03中国公路学报2004年3动态分析程序3.1基本假设(1假设加筋土挡墙足够长,按平面应变问题进行计算,即取单位长度的墙长进行有限元分析;(2采用摩尔库仑屈服准

20、则;(3假定为相关联的流动法则,关系式为fQ,所以有:(dijp=d5f 5ij;(4弹塑性变形在应力和应变之间的完全递增关系为dij=dij2+1-2vEijdkk+d5f5ij(5假设材料屈服后,屈服面是初始屈服曲线的无平移的均匀扩展,即假定应变强化模型为各向同性强化模型。3.2程序的分析模型目前,构成加筋土结构的有限元模型主要有三种形式:整体式、分离式和组合式。在整体式有限元模型中,将筋带离散于整个单元中,并把单元视为连续均匀的材料。筋带对整个结构的贡献,可以通过调整材料的力学参数来体现,例如提高材料的屈服强度、材料的弹性模量。另一种处理方法是一次求得综合单元刚度矩阵,把钢度矩阵改为由筋

21、带和土体两部分组成,李广信等人提出了一种新的处理方法,把筋带的作用当成一个附加周围压力施加于土体。整体式模型可采用各种平面单元,如三角形单元、平面矩形单元、四结点单元或八结点等参单元,也可根据需要采用三维单元。整体式模型明显缺点是无法揭示筋带与土之间相互作用的微观机理,因而计算结果与工程实际有偏差。在分离式模型中,将土体和筋带各自划分为足够小的单元,按照土体和筋带不同的力学性能,选择多种不同的单元型式。将有限元网格划分为:填土单元、筋带单元、面板单元,筋带单元用一维杆单元,面板单元用梁单元,没有考虑筋带与填土之间可能发生的滑移、嵌入情况。事实上,受到外力作用后,筋带与土体之间在相互约定作用下会

22、产生相对滑移,为模拟两者之间这种粘结约束作用和相对滑移,可插入四结点或六结点节理单元,分离式模型可揭示筋带和土体之间相互作用的微观机理。但对加筋土挡墙来说,在墙面板与筋带连接处出现了两个接触面单元交接的特殊单元(无面积且无厚度,给计算带来了相当的不便。组合式分析模型介于整体式与分离式模型之间。组合式模型假定筋带与土体之间的相互粘结很好,不会产生滑移,这类模型有交替正交层系模型,这种模型用于只有单向加筋的情况;弹塑性层板模型,这种模型很适合土工织物加筋土的应力应变特性的分析。3.2.1有限元网格的划分在参考以上三种有限元模型的基础上,结合本研究课题试验观测结果,提出了有别于以上三种模型的新型模型

23、。在试验中发现:由于筋材为土工格栅,具有网孔且表面高低不平,所以筋材在填料之间的滑移面并不在土工格栅表面,而是在距土工格栅表面一定距离的填料中,即土工格栅与周围一定厚度的填土共同工作,这样可视土工格栅与其上下表面一定厚度的土层为筋土单元,从而解决了筋材无厚度问题,又不需要设置接触面单元,不仅使解题方便,而且也符合实际情况。计算模型采用5层筋带,共15个面板单元,70个筋土单元,140个填土单元,共225个单元,256个结点。均采用4结点单元。3.2.2参数的确定面板单元采用混凝土参数,填土单元、筋土单元各参数如表1所示。3.2.3激励力激励力采用模型试验中施加的正弦激励力f(t=a i+b i

24、 sin(w i t=A max i+A min i2+A max i-A min i2sin(w i t式中:a1=50kPa,a2=60kPa,a3=70kPa;b1=10 kPa,b2=20kPa,b3=30kPa;w1=2f1=4,w2 =2f2=8,w3=2f3=12。取a2=60kPa,b2 =20kPa,w1=2f1=4进行计算。表1单元材料性能参数T ab.1P arameters of element material performance单元材料性能参数填土单元筋土单元面板单元弹性模量E/MPa2152508150泊松比v0.30.30.3单元体积质量密度/kNm-318

25、1823强化参数H/MPa100150340内摩擦角/°304560内聚力c/MPa2.02.53.0参考屈服值F01.7321.7701.500 3.3计算结果本次动力分析计算结果数据很多,仅输出t=0(下转第49页13第2期李海深,等:加筋土挡土墙动力特性分析的承压板4,但没能给出量化的指标。从笔者采用的非线性接触分析结果看,采用较厚的垫板和较薄的承压板的组合,既避免了厚钢板的焊接问题,也解决了承压板抗弯不足的问题,受力合理,是一个可行的办法。(4对比两种方法计算的结果发现,采用两种方法计算对于承压板、与承压板相连的板在连接处的应力影响较大,但对于远离连接处的板的应力影响却很小,

26、说明垫板解决的仅仅是承压板及周围的局部受力问题,而对整体受力影响不大。(5腹板对于承压板相当于弹性固定端的作用,所以索力对于腹板的偏心作用使腹板在腹板与承压板相接处的应力较大,但由于腹板幅面较大,应力大的区域衰减较快。(6承压板与腹板的连接的角隅处应力较为集中,应倒圆角。(7 主梁腹板加劲为充分受力受件,对抵抗腹板由于偏心索力产生的面外变形起重要作用,应仔细设计。(8从腹板加劲的受力可以看出,腹板面外变形的影响区域沿纵向分布不宽,因此,主梁腹板为抵抗偏心索力而设的加劲主要布置在锚箱所在区段。(9笔者采用了非线性分析等较为复杂的计算方法,由于受目前计算理论的限制,计算结果应有试验验证,因此应考虑进行钢锚箱的实测试验或缩微模型试验57。参考文献:1林元培.斜拉桥M .北京:人民交通出版社,1994.2严国敏.现代斜拉桥M .成都:西南交通大学出版社,1996.3藤原亨,森山彰,川西直树.日本多多罗大桥上部结构实施设计J .日本本四技报,1998,22(88:1227.4远藤武夫,松本毅,豚原弘光.大规模斜张桥定着构造的疲劳A .日本