ADINA有限元软件中材料本构的二次开发.pdf

第 29 卷第 8 期 岩 土 力 学 Vol.29 No.8 2008 年 8 月 Rock and Soil Mechanics Aug. 2008 收稿日期2007-06-20 作者简介熊玉春，男，1974 年生，博士，工程师，主要从事土动力学与地基基础研究。E-mail xyuchun 文章编号文章编号1000－7598－2008 08－2221－06 ADINA 有限元软件中材料本构的二次开发有限元软件中材料本构的二次开发 熊玉春 1，房营光2 （1. 广东省建筑科学研究院，广州 510500；2. 华南理工大学 土木工程系，广州 510640） 摘摘 要要美国 ADINA R 2. Department of Civil Engineering, South China University of Technology, Guangzhou 510640, China Abstract ADINA software, which is developed by ADINA R for example, there is no Duncan-Chang constitutive model in ADINA commonly used in geotechnical engineering. Taking Duncan-Chang E-B constitutive model as an example, the s of conducting secondary development of material constitutive model and the development process are described. Then the Subroutine Cuser3 is verified by simulating three conventional triaxial compression models. Results show that material constitutive model developed in ADINA has advantages of fast convergence, high precision and convenient preprocessing-postprocessing. Results also show that ADINA has powerful nonlinear solution ability and reliable simulating results. Key words ADINA; nonlinearity; Duncan-Chang E-B constitutive model; secondary development 1 引 言 ADINA[1]（ Automatic Dynamic Incremental Nonlinear Analysis） 软件是在 Bathe K. J.博士的带领 下，由其研究小组共同开发的动力非线性有限元软 件。它具有求解线性问题的能力，还具有分析非线 性问题的强大功能，包括求解结构以及涉及结构场 之外的多场耦合问题。增量方法是数值求解非线性 物理问题本质的方法，对非线性物理问题，计算解 逼近真实解的过程是通过控制增量步逐步实现的， 所谓增量步通常是指载荷增量步或时间增量步。 ADINA 正是基于自动动态增量分析的非线性有限 元软件。所谓自动动态增量，是指当用户给定的时 间步/荷载步太大以致迭代计算出现不收敛时， 它能 自动划分荷载步增量以使计算结果收敛。 ADINA 以其深厚的理论基础与强大的非线性 分析能力和完全的流-固耦合功能，在工程界、科 学研究等领域得到广泛应用。由于 ADINA 在我国 的市场开发起步较晚、可供人们交流学习的资料较 少，这在一定程度上影响了它在我国的普及和发 展。特别是对于处于科学研究前沿的用户来说，由 于科学研究的创新性和前瞻性，几乎所有商业化软 岩 土 力 学 2008 年 件包括 ADINA，都不可能完全满足所有工程问题、 科研课题和所有用户的需求。岩土是一种极其复杂 的工程材料，随土性、排水条件等不同，其力学性 状也千差万别，因此就必须基于该软件的分析平 台，对 ADINA 软件的功能进行扩充以实现不同用 户的不同需求、解决具体的工程问题，这属于基于 ADINA 软件的二次开发问题[2]。 本文以岩土工程中广泛应用的 Duncan-Chang E-B 模型[3]为例，介绍 ADINA 的二次开发过程，并 指出了开发过程中应注意的问题，整个过程力求思 路清晰、步骤明确，以期对从事 ADINA 二次开发 的用户能有所裨益。 2 用户子程序 ADINA 为用户提供了丰富的二次开发功能，允 许用户方便地定义各种用户功能，常见的诸如材料 本构算法、材料破坏准则、接触摩擦形式、断裂力 学判据和开裂扩展规律、用户边界条件/荷载等等。 ADINA8.1 的开发环境为 Compaq Visual Fortran 6.6A，开发时通过 ADINA 提供的 Makefile 文件自 动链接各个*.f 文件生成动态链接库文件，针对不 同求解器的开发，提供不同的 Makefile 文件，包括 ADINA 模块adusr.dll；ADINA-T 模块atusr.dll； ADINA-F 模块afusr.dll；ADINA-FSI 模块 adfusr.dll；ADINA-TMC 模块adtusr.dll。ADINA 二次开发过程如下 （1）备份 C\Program Files\ADINA\ADINA Sys- tem 8.1\bin 目录下的 dll 文件。 （2）修改 C\Program Files\ADINA\ADINA Sys- tem 8.1\usrdll目录下的Fortran源文件为用户材料本 构源文件。ADINA 提供了两种本构示例的 Fortran 源文件，一类是 Ovl3*.f 文件，用于 2D 模型单元 的本构；一类是 Ovl4*.f 文件，用于 3D 模型单元的 本构。 （3）编辑 Makefile 文件，修改其中的参数 MAT2D_OBJ ovl30u_usr.obj MAT2D_OBJ ovl40u_usr.obj 其中usr 为用户选用或定义的文件名。 （4）编译生成新的 dll 文件，编译过程如下 cd c\program files\Microsoft Visual Studio\ df98\bin cd c\program files\ ADINA\ ADINA system 8.1\ usrdll （5）nmake /f Makefile.，其中 替换 为用户要编译的 dll 链接库文件。 （6）将用户生成的 dll 文件拷贝到 bin 目录下。 编译、调试通过后就可以进行计算分析，此时 在 ADINA AUI 界面中采用 User-Supplied 模式，输 入用户模型参数，其中 CTI 数组用于存储温度无关 的材料常数；CTD 数组用于存储温度相关的材料 常数；LGTH1 为用户使用的实型变量 ARRAY 的 长度；LGTH2 为用户使用的整型变量 IARRAY 的 长度。这些参数与 ADINA 子程序中的变量一一对 应，用于材料本构的计算。 ADINA 为方便用户进行材料本构的二次开发， 采用 FORTRAN 语言接口方式，提供了 20 多个用 户子程序即 SUBROUTINE CUSER2（开发二维本 构）和 SUBROUTINE CUSER3（开发三维本构） 。 材料本构关系须按照 FORTRAN语言的算法来编写 并嵌入 SUBROUTINE CUSER2 或 SUBROUTINE CUSER3 子程序中。以本文用到的 SUBROUTINE CUSER3 子程序为例，每一个 3D 单元的应力积分 点通过整型变量 KEY，来控制在不同阶段调用该子 程序完成下面 4 类操作 （1）当 KEY1 时，在进行分析时被调用一次。 这一步主要是对积分点的变量进行初始化，在整个 数据输入阶段只被执行一次。ADINA 自动将数组 STRAIN，STRESS，ARRAY 和 IARRAY 等初始化 为 0。如果 ARRAY 和 IARRAY 数组初始值不为 0， 则用户在编译时须将其设置为所需的初始值。 （2）当 KEY2 时，被每一个荷载步/时间步的 每一个子步调用。这一步主要是在求解平衡方程获 得节点位移以及节点位移求出后的应力计算阶段被 调用来进行单元应力积分点的应力计算。单元每一 应力积分点的应力计算按照嵌入 SUBROUTINE CUSER3 中的用户代码来执行。除进行应力计算 外，如果用户使用 ARRAY 和 IARRAY 来存储历 史相关变量，则 KEY2 时，用户须更新 ARRAY 和 IARRAY 数组。 （3）当 KEY3 时，被调用来计算新的切线系数 矩阵。 这一步是计算本构矩阵 D 用以估算切线系数 矩阵，此时传递给 SUBROUTINE CUSER3 子程序 的变量值是 KEY2 程序段最后一个子步的计算值。 因此这些变量在该阶段不再计算、更新。 （4）当 KEY4 时，被调用来输出应力计算阶段 的单元响应。同样该阶段不进行应力计算，也不对 任何变量进行更新，只按照用户代码制定的格式输 出用户需要的单元结果。 ADINA 中材料的变形描述有 3 种 （1） 材料非 线性（MNO） ，适用于无穷小位移和应变； （2）完 全拉格朗日（TL） ，适用于位移和旋转变形可能比 2222 第 8 期 熊玉春等ADINA 有限元软件中材料本构的二次开发 较大，而位移视材料特性可以是大位移也可以是小 位移的情况； （3）更新拉格朗日（ULH） ，适用于 大应变分析。ADINA 约定剪应变分量 ij ε采用工程 剪应变 ij γ的形式，在材料弹塑性分析中，当需要由 剪应变求剪应力或由剪应力求剪应变时，需要把工 程剪应变转换为 Cauchy 剪应变或把 Cauchy 剪应变 转换为工程剪应变来存储。 另外，在应用用户材料进行有限元分析时，需 要用户设定温度荷载，当用户材料不涉及温度效应 时，可将温度荷载设定为 0。 3 本构模型的数学表达式 在 ADINA 中应力应变的正负号规定与土力学 不同，ADINA 规定以拉应力为正，压应力为负；应 变以拉伸为正，压缩为负。相应地，土力学中大主 应力 1 σ为 ADINA 中最小主应力 3 σ−，而小主应力 3 σ为 ADINA 中最大主应力 1 σ−。因此，编程开发 时应改写本构模型[4]的数学表达式，使之符合 ADINA 约定的数学表达式。 Duncan-Chang E-B 模型中切线弹性模量 t E的 数学表达式为 2 tf 1 i EER S− （1） 式中 i E为初始弹性模量；S为应力水平； f R为破 坏比，其定义分别为 131 a a1 13 f f 13 ult 1sin 2 cos2sin m i EkpS pc R ϕ σσσ ϕσϕ σσ σσ ⎫ ⎡⎤−−− ⎪ ⎢⎥ − ⎪ ⎣⎦ ⎪ ⎬ − ⎪ ⎪ − ⎪ ⎭ ；； （2） 式中k、m为试验确定的参数；ϕ为内摩擦角；c 为土的黏聚力； a p为标准大气压力。 切线体积模量 t B的数学表达式为 1 tba a m Bk p p σ⎡⎤− ⎢⎥ ⎣⎦ （3） 式中 b k为试验确定的参数。 考虑到粗粒料的内摩擦角ϕ会随围压的变化而 变化，摩尔包线往往不是直线，故采用下式来计算 内摩擦角 1 0 a lg p σ ϕϕϕ ⎡⎤− −∆ ⎢⎥ ⎣⎦ （4） 式中 0 ϕ为围压等于一个标准大气压时的内摩擦 角；∆ϕ为反映内摩擦角ϕ随围压 3 σ增大而减小的 参数。 计算中，如果输入参数∆0ϕ，则不考虑围压 对内摩擦角的影响。 在卸荷、再加荷情况下，应力-应变关系接近直 线，这时的弹性模量 ur E仅取决于周围压力，而与 主应力差无关，形式如下 1 urura a n Ek p p σ⎡⎤− ⎢⎥ ⎣⎦ （5） 式中 ur k、n为试验确定的参数。 加载函数的选择是变弹性模型中的一个难题， Duncan等[5 ，6]经过多年探索，于 1984年提出加载 函数 1 4 131 l 13 fa f p σσσ σσ ⎡⎤−− ⎢⎥ − ⎣⎦ （6） 当 l f大于历史上最大值 lmax f时，视为加载，否 则视为卸载或再加载。由于卸荷模量与加荷模量之 比达2个数量级以上，常引起迭代过程的不稳定， 具体计算时采用下面经验方法当 llmax ff≥，切线 弹性模量按式（1）计算；当 llmax 0.75ff≤，切线弹 性模量按式（5）计算；当 lmaxllmax 0.75fff，切 线弹性模量按下式线性内插法计算 lmaxl turt lmax 0.25 ff EEEE f − − （7） 由式（1）（7）即可确定 Duncan-Chang E-B 模型本构矩阵中的所有元素。在有限元分析中，单 元的增量应力-应变关系可用矩阵形式表示为 {}[]{}∆∆σDε （8） 式中[]D为单元切线系数矩阵。 根据 ADINA 应力、应变分量的约定，对于三 维问题，存储应力的数组 STRESS 与 Cauchy 应 力各分量的关系为（存储应变的数组 STRAIN 与 Cauchy 应变的关系与此类似，故不再列出） STRESS1 xx σ；STRESS2 yy σ；STRESS3 zz σ； STRESS4 xy σ； STRESS5 xz σ； STRESS 6 yz σ。 对于 Duncan-Chang E-B 三维模型问题，本构 系数矩阵为 IIIJIJ IJIIIJ IJIJII JJ JJ JJ 000 000 000 [] 00000 00000 00000 ddd ddd ddd d d d ⎡⎤ ⎢⎥ ⎢⎥ ⎢⎥ ⎢ ⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎣⎦ D （9） 对于二维平面应变问题，存储应力的数组 STRESS与Cauchy应力各分量的关系与三维情况不 2223 岩 土 力 学 2008 年 同STRESS1 yy σ；STRESS2 zz σ；STRESS 3 yz τ；STRESS4 xx σ。 对于平面应变问题，本构系数矩阵为 IIIJIJ IJIIIJ JJ IJIJII 0 0 [] 000 0 ddd ddd d ddd ⎡⎤ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎣⎦ D （10） 式（9） ， （10）中系数矩阵的各非 0 元素分别为 IIIJ JJ 3 33 3 99 3 9 BBEBBE dd BEBE BE d BE −⎫ ⎪ ⎪−− ⎬ ⎪ ⎪ −⎭ ；； （11） ADINA 中子步应变增量用数组 DEPS（不包括 子步热应变增量（存储于数组 DEPST 中） 。若为弹 塑性分析，则 DEPS 数组中包括子步塑性应变增 量；若本构系数矩阵采用弹性系数矩阵，则在 KEY 2 进行应力积分时，应从数组 DEPS 中减去所有的 子步非弹性应变增量）表示，若子步应力增量用 DS 数组表示（DS 数组中的应力各分量与数组 STRESS 的各分量对应） ，并把数组 DEPS 和 DS 视 为向量，且向量各分量等于对应的数组元素，则增 量应力-应变关系为 {}[]{}DEPSDSD （12） 即为 ADINA 计算单元应力积分点每一子步应力增 量时的增量应力-应变形式。 4 Duncan-Chang E-B 模型程序概要 一个完整的 ADINA 材料本构用户子程序除了 第一节介绍的由关键字 KEY 控制的 4 个程序模块 外，还包括变量声明（包括数组的维数和大小）和 变量说明以及在 KEY2 时将数组 CTI 和 CTD 中的 值赋值给相关变量。下面针对 Duncan-Chang E-B 模型介绍本文主要的编程过程 （1）对数组 ARRAY 和 IARRAY 初始化，将数 组 CTI 中的材料参数赋值给相关变量； （2）根据当前子步起始时的单元积分点的应 力，计算主应力，并判别大小主应力计算主应力 差。在计算主应力时可调用 ADINA 的子程序 XJ123（三维）或 XJ123（二维）由应力分量计算出 应力张量不变量，然后求解一元三次方程解出主应 力。 （3）计算当前子步的应力水平并进而计算切线 弹性模量、卸荷弹性模量和切线体积模量。 （4）计算加载函数并与最大值进行比较，根据 加卸载准则确定子步的切线弹性模量和体积模量。 （5）由切线弹性模量和体积模量计算当前子步 的系数矩阵。 （6）根据子步应变增量计算当前子步的应力增 量，并更新应力数组 STRESS。 （7）进行单元破坏（拉伸破坏或剪切破坏，取 决于用户确定的破坏准则）判别与应力修正。 以上除（1）中数组 ARRAY 和 IARRAY 的初 始化在 KEY1 程序模块完成外，其余都是在 KEY 2 模块进行的。 （8）根据当前荷载步最后一个子步的应力状态 重新计算系数矩阵用以求解下一荷载步的分析。该 部分内容是在 KEY3 程序模块进行的。 5 程序的验证 本文开发完成了 Duncan-Chang E-B 模型的 SUBROUTINE CUSER3 子程序，经编译调试通过 后计算了若干模型问题，以测试 ADINA 的非线性 计算能力和精度。模型计算参数取自文献[7]，见 表 1，取∆0ϕ。有限元计算模拟常规三轴试验的 加载过程。加载方式、模型尺寸与常规三轴试验一 致。有限元计算模型如图 1 所示，共有 2 640 个单 元和 3 133 个节点，按 Gauss 积分法计算单元积分。 表表 1 Duncan-Chang E-B 模型参数模型参数 Table 1 Parameters for Duncan-Chang E-B model c /kPaϕ0 /（）k /kPakur /kPa n Rf kb /kPam 0.00642.421215.5 231.1 0.92 0.613 151.010.075 表 2、图 2 和图 3 分别给出了围压为 200 kPa、 最大主应力差为 600 kPa、采用不同荷载步数时轴 向应变 a ε的数值解和 Duncan-Chang E-B 模型理论 解及试样轴向位移应力云图。由表 2 和图 3 可知， 在低荷载水平时，ADINA 的数值模拟结果与模型 理论解的误差较小，并且随着荷载步的增加，数值 模拟结果越趋近于理论值，即荷载步数增大，模拟 的精度越高，可见模拟结果与理论值具有很好的一 致性。 表表 2 不同荷载步的计算结果与理论值的比较不同荷载步的计算结果与理论值的比较 Table 2 Results of different load steps and theoretical values 主应力差 /kPa 加载步 /步 0100 200 300 400 500 600 30 0.00.2520.5440.886 1.288 1.7612.304 60 0.00.2600.5660.924 1.345 1.8462.451 120 0.00.2630.5710.933 1.379 1.9222.548 理论值0.00.2650.5750.941 1.390 1.9402.640 2224 第 8 期 熊玉春等ADINA 有限元软件中材料本构的二次开发 图图 1 有限元计算模型及轴向应力（单位有限元计算模型及轴向应力（单位kPa）） Fig.1 Finite element calculation model and axial stress unit kPa 图图 2 三轴试样轴向位移（单位三轴试样轴向位移（单位mm）） Fig.2 Axial displacement of the triaxial Specimen unit mm 图图 3 不同荷载步时的 （不同荷载步时的 （σ σ1-σ σ3）） -ε εa曲线曲线 Fig.3 The σ σ1-σ σ3-ε εa curves of different load steps 图 4 为采用相同模型参数，在围压和最大主应 力差分别等于 100 kPa 和 300 kPa，200 kPa 和 600 kPa，300 kPa 和 900 kPa，且荷载步数分别为 30， 60，90 步条件下，根据二次开发的程序计算的数值 结果即 A1、A2和 A3曲线和 Duncan-Chang E-B 模型 的理论解（即 E1、E2和 E3曲线，其中下标 1∼3 分 别表示围压为 100，200，300 kPa 的固结状态） 。可 见在不同围压下，基于 ADINA 开发的非线性有限 元程序对 Duncan-Chang E-B 模型都有很好的模拟 能力。 图图 4 不同围压时的 （不同围压时的 （σ σ1-σ σ3）） -ε εa曲线曲线 Fig.4 The σ σ1-σ σ3-ε εa curves of different confining pressures 图 5 为由二次开发的程序对加载、卸载的模拟 结果。采用的模型参数与上文相同，围压和最大主 应力差分别为 200 kPa 和 600 kPa，荷载步数为 60 步。由图 5 可知，加载到主应力差为 360 kPa 时沿 直线卸载，卸载后再加载时应力轨迹为曲线，这是 由于模型在卸载时采用常模量，而在加载时根据应 力水平采用变模量的缘故。图 5 亦表明，基于 ADINA 开发的 Duncan-Chang 模型对于加卸载情况 也有很强的模拟能力。 图图 5 加卸载时的 （加卸载时的 （σ σ1-σ σ3）） -ε εa曲线曲线 Fig.5 The σ σ1-σ σ3-ε εa curves of loading and unloading 6 结 语 通过对非线性 Duncan-Chang E-B 模型的开发 与程序测试表明，在 ADINA 中进行用户材料本构 下转第下转第 2240 页页 600 500 400 300 200 100 0 0.0 0.6 1.2 1.8 2.4 εa / （σ1-σ3）/kPa −320.0 −400.0 −480.0 −560.0 −640.0 −720.0 −800.0 −0.133 −0.400 −0.667 −0.933 −0.200 −0.467 −0.733 600 500 400 300 200 100 0 30 荷载步 60 荷载步 120 荷载步 E-B 模型理论解 0.0 0.5 1.0 1.5 2.0 2.5 εa / （σ1-σ3）/kPa 900 800 700 600 500 400 300 200 100 0 A3 A2 A1 E3 E2 E1 0.0 0.6 1.2 1.8 2.4 3.0 εa / （σ1-σ3）/kPa 2225 岩 土 力 学 2008 年 stability of slope[J]. West China Exploration Engineering, 2005, 6 213－214. [2] 栾茂田， 黎勇， 杨庆. 非连续变形计算力学模型在岩体 边坡稳定性分析中的应用[J]. 岩石力学与工程学报， 2000，19（3） 289－294. LUAN Mao-tian, LI Yong, YANG Qing. Discontinuous deation computational mechanics model and its application to stability analysis of rock slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2000, 193 289－294. [3] 陈昌彦， 王思敬， 沈小克. 边坡岩体稳定性的人工神经 网络预测模型[J]. 岩土工程学报，2001，23（2） 157 －161. CHEN Chang-yan, WANG Si-jing, SHEN Xiao-ke. Predicting models to estimate stability of rock slope based on artificial neural network[J]. Chinese Journal of Geotechnical Engineering, 2001, 232 157－161. [4] 卢才金，胡厚田，徐建平，等. 改进的 BP 网络在岩质 边坡稳定性评价中的应用[J]. 岩石力学与工程学报， 1999，18（3） 303－307. LU Cai-jin, HU Hou-tian, XU Jian-ping, et al. Application of improved back-propagation network to the uation of railway rock slope[J]. Chinese Journal of Rock Mechanics and Engineering, 1999, 183 303－307. [5] KOHONEN T. Self-organized ation of topologically correct feature maps[J]. Biological Cybernetics, 1982, 43 39－69. [6] KOHONEN T. Self-organizing Maps2nd edition[M]. Berlin Springer-Verlag, 1977. [7] 吕小平，胡厚田. 铁路边坡稳定性预测系统的研究[J]. 铁道学报，1993，15（1） 94－102. L Xiao-ping, HU Hou-tian. A stability analysis for railway slopes[J]. Journal of the China Railway Society, 1993, 151 94－102. [8] 岩土工程手册编委会. 岩土工程手册[M]. 北京 中国建筑工业出版社，1994. 上接第上接第 2225 页页 的二次开发是可行的，ADINA 具有很强的非线性 求解能力。本文给出的基于 ADINA 有限元软件的 本构二次开发思路和方法，可为在 ADINA 中开发 其他材料本构模型及其他用户子程序提供借鉴和参 考。ADINA 有限元软件的二次开发，有助于将这 一先进的非线性分析软件进行改进和功能扩充，从 而满足不同的用户需求，推动该软件在科学研究和 实际工程中的广泛应用。 参参 考考 文文 献献 [1] BATHE K J. ADINA/ADINAT 使用手册[M]. 赵兴华译. 北京机械工业出版社，1986. [2] 汪仁和，李栋伟，王秀喜. 改进的西原模型及其在 ADINA 程序中的实现[J]. 岩土力学，2006，27（11） 1 954－1 958. WANG Ren-he, LI Dong-wei, WANG Xiu-xi. Improved Nishihara model and realization in ADINA FEM[J]. Rock and Soil Mechanics, 2006, 2711 1 954－1 958. [3] DUNCAN J M, CHANG C Y. Nonlinear analysis of stress and strain in soils[J]. Journal of Soil Mechanics and Foundation Division, 1970, 96SM5 1 629－ 1 653. [4] 钱家欢，殷宗泽. 土工原理与计算[M]. 北京中国水 利电力出版社，1993. [5] 陈火红， 尹伟奇， 薛小香. MSC.Marc 二次开发指南[M]. 北京科学出版社，2004. [6] 朱百里，沈珠江. 计算土力学[M]. 上海上海科学技 术出版社，1990. [7] 冯卫星，常邵东，胡万毅. 北京细砂土邓肯-张模型参 数试验研究[J]. 岩石力学与工程学报，1999，18（3） 327－330. FENG Wei-xing, CHANG Shao-dong, HU Wan-yi. Experimental study on parameters of Duncan-Chang model for Beijing fine sandy soil[J]. Chinese Journal of Rock Mechanics and Engineering, 1999, 183 327－ 330. 2240