等效线性动力模型在ADINA中的实现

1、扬州大学学报(自然科学版)Journal of Yangzhou U niversity (Natural Science Edition)第11卷第4期2008年11月Vol. 11 No. 4Nov. 2008等效线性动力模型在ADINA中的实现费康1收稿日期：2008- 04- 01基金项目：扬州大学高层次人才科研启动基金资助项目 联系人，E-mai l： kfeiyzu. edu. cn,郑憩2（1.扬州大学岩土工程研究所,江苏扬州225009; 2.北京林业大学 园林学院，北京100083）摘要：在有限元软件ADINA的开发平台上开发岀三维等效线性动力模型.应用该模型分别与正弦应力波

2、作用下粘弹性体应变理论解和一维地基场地地震反应分析的SHAKE91计算结果进行对比，结果表明，编制的模型能较好地反映岩土体剪切模量与应力的关系、岩土体的滞 回阻尼特性和剪切 模量、阻尼比随剪应变变化的特性，为岩土地震反应分析提供了一种可供选择的手段.关键词：等效线性模型；三维有限元；地震反应分析；ADINA有限元程序中图分类号：TU 435文献标识码：A文章编号：1007- 824x（2008） 04- 0074- 05岩土地震反应分析一般可分为等效线性和非线性两大类1.非线性方法理论上更符合物理概念,能反映孔压积累、残余变形等现象，非线性分析的准确程度依赖于所选用的本构模型目前分析程序中的非

3、线性模型要么简单，如摩尔库仑模型2;要么复杂，如边界面模型3和多重塑性屈服面模型4.这些模型需要大量的参数和较为复杂的参数确定过程，加上计算费用和计算时间效率的限制，非线性分析往往只用于一维场地反应分析，在大型复杂结构的三维分析中还较少应用而等效线性模型则通过迭代的方法来近似地反映土体的非弹性和非线性这种方法计算效率高，且大多数情况下都能提供较合理的结果，目前等效线性方法仍然是生产应用中的主流.5-6在基于等效线性方法的程序中，目前仅有一些程序提供了前、后处理的功能，比如二维的 QUAKE/ W;而大部分软件都是基于DOS系统下的一维或二维程序，前、后处理工作烦琐，使用还不是很方便另外，一些大

4、型有限元商业软件如ADINA、ABAQU S等具备友好的用户界面和强大的建模功能，并能处理开挖、分期施工、接触等复杂的工程问题，但缺点是并未提供岩土体专用的静、动力本构模型，在岩土地震工程中的应用范围还很有限一些学者结合商业有限元软件进行了岩土体静力本构模型的二次开发，如邓肯非线性模型7-8、沈珠江双屈服面模型等 ，但在ADINA中开发岩土 体动力本构模型的研究还较少为此，笔者拟利用ADINA提供的二次开发接口 ，将三维等效线性模型编入到ADINA中，以扩展ADINA软件的应用范围，可望为岩土地震反应分析提供可供选择的方 便、实用的工具1 ADINA中等效线性模型的实现1.1用户自定义材料接口

5、介绍ADINA 8. 1版本的开发环境为 Compaq Visual Fortran 6. 6A .用户自定义材料的过程包括：1将ADINA提供的ovl3* .f（用于二维单元）或ovl4* . f（用于三维单元）中的相应代码替换为自定义 的本构关系程序语句；o修改ADINA提供的面对不同求解器的Makefile文件，编译后形成动态链接库文件.对本文中三维等效线性模型的开发而言，须令Makefile. adusr文件中的参数 MAT 3D_OBJ= ov|40lt vel. obj，然后在DOS命令窗口下通过 nmake/f makefile. asusr命令编译生成动态链接 第4期费康等：等

6、效线性动力模型在 ADINA中的实现79库 adusr. dll.1.2用户自定义材料子程序的结构三维用户自定义材料子程序(以下简称CUSER3)的主要目的是对单元的每个积分点根据 ADINA传来的参数变量计算出应力增量，并给出模量矩阵.下面介绍几个最重要的参量(数组):1 CTI (NCTI):用户自定义材料与温度无关的本构模型参数数组 ；o ARRAY (LGTH 1):实数型历史 变量数组，常用于塑性应变、应力水平等数据的存储；？ IARRAY (LGTH2):整数型历史变量数组， 常用于程序控制或其他用途 ；？ STRAIN( 6): t+ $t时刻的应变分量数组；？ STRESS(

7、6):上一个 子增量步的应力；？ D(6,6):模量矩阵.用户自定义材料子程序CUSER3的代码结构分为4部分,由参数KEY控制：1 KEY = 1,在分析开始时调用，对每个积分点存储的数据进行初始化，STRAIN ( 6) , ST RESS ( 6) , ARRAY (LGT H1) , IARRAY( LGT H2)数组默认为0,若上列数组不为 0,则须用户自己通过语句控制；0KEY = 2,在每一个子增量步调用 ，根据应变增量求解并更新单元各积分点的应力，同时还须更新ARRAY (LGTH 1)和IARRAY (LGT H2)数组；？ KEY = 3,在求取结构总切线刚度矩阵时调用，此

8、时需要CUSER3给出各单元的切线模量矩阵D(6,6); ? KEY= 4,无需进行任何计算，将所需输出的结果按照一定格式输出.1. 3等效线性模型理论等效线性模型实际上是基于粘弹性理论，即用粘弹性Kelvin模型来反映土体在周期荷载下的滞回性.土体能量损耗(或阻尼比)随剪应变变化的特性和滞回曲线斜率(或剪切模量)随剪应变变化的特性,通过改变模型参数来反映.具体过程为：先假设一个初始值阻尼比和剪切模量，计算过程中记录每个单元经历的最大剪切应变，根据试验求得的剪切模量、阻尼比与周期剪应变之间的关系曲线确 定新的剪切模量和阻尼比，再根据新的材料参数进行计算，整个过程重复数次，直到材料性质不再发生变

9、化为止粘弹性模型的应力应变关系可表示为S= GC+ GC,導效滞回曲线图1滞回曲线Fig. 1 Hysteresisstresstrain curve式中G为剪切模量，G为剪切粘滞系数，S为剪应力，C 为剪应变.模型对应的土体滞回曲线和土体真实的滞回 曲线绘制于图1中.由图可见，土体真实滞回圈最大应力 和最大应变发生的时刻是对应的，而等效线性模型的滞回圈为一个椭圆，激励和反应之间有一定的相位差，与真实情况有所区别，但对于阻尼水平较低、椭圆狭长的情况，等 效线性模型还是适用的.原点和最大剪应变点之间连线 的斜率即代表剪切模量 G,而剪切粘滞系数为(1)G= 2GK/X,式中K为阻尼比，X为圆频率

10、.国内外学者通过大量室内试验10给出的阻尼比、剪切模量和剪切应变的关系曲线常可表示为G/ Gmax =1 + Q Cef , K K°ax = 1- G/ Gmax,(2)fin式中 6为参考应变，Kx为最大阻尼比，Gax为最大剪切模量，可由公式 Gmax = kpa. R/ P确定，式 中k, n分别为剪切模量系数和指数，pa为大气压力，R)为静平均压力.对于三维的情况，可将模量矩阵用体积模量K和剪切模量G表示,体积模量K可由剪切模量 G(3) 由公式(1)和(3)可知G/G=求得：K= 2(1+ T)G/ 3(1- 2T),式中T为泊松比.类似于式(1),材料的体积粘滞系数为G&

11、lt; = 2KK/X.这里的体积粘滞系数并不反映由孔压引起的体积变形随时间的变化1.4等效线性模型子程序的编写1）如何反映模量与应力之间的关系 一般认为，岩土体在地震过程中的动模量取决于震前岩土体的静应力即初始应力的大小但目前ADINA中的用户自定义材料无法应用其自带的初始应力选项，因而必须从外部导入,具体操作为：在KEY= 1时，通过open语句建立与静应力结果文件的通 道,读入各单元的静应力大小静应力结果文件可以是 ADINA的计算结果，也可以由其他程序结果 经简单处理后形成2）如何反映阻尼项 ADINA中不允许用户直接修改总阻尼矩阵，因此无法根据刚度矩阵直接给出总阻尼矩阵另外，若采用R

12、ayleigh阻尼，则一种材料只允许设定一组 Rayleigh阻尼参数，这就 意味着尽管同一种材料中不同位置的岩土体单元在动力过程中所经历的应变水平（最大剪切应变）不同，但它们的阻尼比是相同的，这不符合岩土体的基本动力性质 .因此，要正确反映阻尼项，只有在模 型内部解决，即根据等效线性模型的应力应变关系来反映滞回特性和能量损失3）如何反映模量、阻尼比随应变水平的变化.由前所述,等效线性分析中每一次迭代过程的模量和阻尼是保持不变的，在每个迭代过程中记录各单元经历的最大剪应变，然后乘以一等效系数，如0. 65,再按公式（2）确定出与应变水平相协调的模量和阻尼比后重新计算，直到材料参数不发生变化为止

13、.理论上，这个过程可以内嵌入到用户子程序中让ADINA自动运行，但对于复杂结构（如三维土石坝等）的动力分析，这样会占用大量的计算资源和时间，普通PC用户难以承受，且不利于对中间结果的监控和调试.为了避免这一问题，笔者采用了如下的做法:1 KEY= 2过程中储存并更新各 单元的最大剪应 变；o KEY= 4时,将历史中最大剪应 变结果输出到 一外部文件；？另外编一 FORTRAN程序,根据各单元的应变水平确定G/Gmax和K/K-;?保持ADINA模型和数据文件不变,在读入初始应力的同时读入G/Gmax和KK-,重新计算，直到前后两次分析得到的应变水平误差在允许范围之内2模型的验证为验证编制的模

14、型是否正确，笔者将本文方法与SHAKE91程序的计算结果进行了对比.计算工况为一坐落在基岩上的水平地基，厚45 m,地基由沙土和黏土组成，具体分层见图2,33021 m内土体容重为 19.66 kN m , 2145 m内土体容重为 20. 45 kN m .图2同时给出了 由现场剪切波速试验推导的最大剪切模量随深度的分布.图3给出了黏土和沙土的剪切模量比、 阻尼比与剪应变的关系曲线，图4是基岩的输入加速度时程曲线，最大峰值为0. 1gn, gn为标准重力加图2 土层分布Fig. 2 Soil profile速度.h黏1：羽切HUE-口沙土唧切模就比 粘土阻尼比.o沙土阻足岀图3剪切模量比、阻

15、尼比与剪应变的关系曲线F ig. 3 Shear modulus ratio and damping ratio versus shear stra in图4输入加速度时程曲线Fig. 4 Tine history of input ground accelerationFig. 5 Time history of acceleration at soil surface图5 土层表面加速度反应时程曲线在ADINA三维分析中，取一长1 m （x方向）、宽1 m （y方向）、高45 m （z方向）的土柱进行分 析.土柱底面固定，顶面自由，整个土柱中只允许发生x方向变形，以此来防止单元的弯曲和模拟

16、水平地基的特性.地震波在土柱底部沿x方向输入.土柱采用三维八节点块体单元划分，考虑到土层和最大剪切模量的变化，将地基沿高度划分为30个单元.土体材料为用户自定义的等效线性模型，最大剪切模量直接采用图2中的值根据图3的曲线拟合，黏土和沙土的参考应变分别取0.234%和0. 06% ,最大阻尼比分别取0. 3和0.2.为方便比较，将SHAKE 91中基岩的剪切模量取为2 000GPa来模拟底部的刚性约束利用SHAKE 91和本文方法求解的土层顶面的绝对加速度的最大值分 别为0.29gn和0. 31gn，误差约为6%.另外，两者所得的地震反应加速度时程曲线比较吻合，为清晰起见，图5中只给出本文方法的

17、计算结果 .图6给出了最大地震反应加速度沿深度的分布 ，两种计算方法呈现出相同的变化趋势 ，即土层底 部以上近1/3范围加速度变化较小，在此范围之外，加速度放大较明显，最大加速度出现在土层顶面， 放大系数分别为2. 9和3.1.图7是将输入地震波的峰值加速度 ah调整到0. 001,0. 01,0. 1, 1gn情况 下计算得到的加速度放大系数 amax/ah随峰值加速度的变化曲线，整体呈现出下降的趋势.这是因为随着峰值加速度的增加，地震过程中剪切应变增加 减小,这体现了岩土体材料的非线性.,相应的剪切模量减小，阻尼比增加，因而放大系数图&加速度沿深度前分布 Fig, 6 Distri

18、bmicms of areel eration along the depth图7 土层表面加速度竝大乘数 Fiji- 7 Acceleration amplifications ut thv soil sxirfiice参考文献1 OKAMOTO T . Recent trend for earthquake induced residual settlement of rockfill dam and someconsideration onaffecting fact ors Cl/ Proceedings of the 4th Inter national Conference on

19、Dam Engineering： New Developmentsin Dam Engineering. Rotterdam: BalkemaA A, 2004: 705-716.2 FENG 乙 TSAI PH, CHANG Y H, et al. Dynamic responseof Li-yu-tan earth dam subjected to the 1999 ChiChi earthquake in T aiwan C/ Proceedings of the 4th I nter national FLAC Symposium. Madrid: Itasca, 2006: 271-

20、289.3 WANG Z L, MAKDISI F I, EGAN J. Practical applications of a nonlinear approach to analysis of earthquake- induced liquefaction and deformation of earth structures J. Soil Dyn & Earthquake Engin, 2006, 26( 2): 231252.4 刘汉龙，丰土根，高玉峰，等.砂土多机构边界面塑性模型及其试验验证J.岩土力学，2003, 24(5) : 696700.5 KRAMER S L

21、, PAULSEN S B. P ractical use of geotechnical site response models C / / Proceedings of I nternational Workshop on Uncertainties in Nonlinear Soil Properties and their Impact on M odeling Dynamic Soil Response Berkeley: University of California, 2004: 162-165.6 SINGH R, ROY D, DAS D. A correlation f

22、or permanent earthquake-induced deformation of earth embankments J. Engin Geol, 2007, 90(3): 174-185.7 徐远杰，王观琪，李 健，等.在ABAQUS中开发实现Duncan-Chang本构模型J.岩土力学，2004, 25(7): 1032-1036.8 张 欣，丁秀丽，李术才.ABAQUS有限元分析软件中 Duncan-Chang模型的二次开发J.长江科学院院报， 2005, 22(4) : 45-47.9 岑威钧，朱岳明.基于ABAQUS的土石料本构模型二次开发及其应用J.水利水电科技进展，2

23、005, 25(6): 78-81.10 ZHANG J F, ANDR US R D , JUANG C H. Nor malized shear modulus and material damping ratio relationships J . J Geotech & Geoenviron Engin, 2005, 131(4): 453-464.Impleme ntati onof equivale nt lin ear model in ADINA1* 2FEI Kang , ZHENG Qi(1. Geotech Res Inst, Yangzhou Univ, Yangzhou 225009, China;2. Coll of Landscape Arch it, Beijing For U niv, Beijing 100083, China)Abstract: A three dime nsional model based o n the con cept of equivale nt lin ear earthquake an alysis is built into ADIN