Cosserat弹塑性模型在ABAQUS中的数值实施
基于ABAQUS的内压厚壁圆筒的弹塑性分析
基于ABAQUS的压厚壁圆筒的弹塑性分析学院：航空宇航学院专业：工程力学指导教师：：学号：1. 问题描述一个受压的厚壁圆筒（如图1），半径和外半径分别为mm a 10=和mm b 15=（外径与径的比值2.15.11015b >==a ），受到均匀压p 。
材料为理想弹塑性钢材（如图2），并遵守Mises 屈服准则，屈服强度为MPa Y 380=σ，弹性模量GPa E 200=，泊松比3.0=υ。
图1 压作用下的端部开口厚壁圆筒 图2 钢材的应力-应变行为首先通过理论分析理想弹塑性材料的厚壁圆筒受压作用的变形过程和各阶段的应力分量，确定弹性极限压力e p 和塑性极限压力p p ；其次利用ABAQUS 分析该厚壁圆筒受压的变形过程，以及各个阶段厚壁筒的应力分布，与理论分析的结果进行对比，验证有限元分析的准确性。
2. 理论分析2.1基本方程由于受到压p 的作用，厚壁圆筒壁上受到径向压应力r σ、周向压应力θσ和轴向应力z σ的作用，由开口的条件可推出0=z σ。
因为这是一个轴对称问题，所有的剪应力和剪应变均为零。
平衡方程和几何方程用下式表示：0-=+rd d r r r θσσσ （1）r u dr du r r r ==θεε, （2） 弹性本构关系为：()()r r r E E συσεσυσεθθθ****1,1-=-= （3） 由于此问题为平面应变问题，所以上式中2*1υ-=E E υυυ-=1* 相应的边界条件为：0,=-===b r r a r r p σσ （4）2.2弹性阶段根据弹性力学中的应力解法：取应力分量r σ，θσ为基本未知函数，利用平衡方程和应力表示的协调方程联合求解，可得应力分量的通解⎪⎪⎩⎪⎪⎨⎧=+=221221-r C C r C C r θσσ 将边界条件带入可得应力分量为：⎪⎪⎩⎪⎪⎨⎧⎪⎪⎭⎫ ⎝⎛+-=⎪⎪⎭⎫ ⎝⎛-=11--2222222222r b a b p a r b a b p a r r σσ （5） 因为b r a ≤≤，所以00>≤θσσ且r ，可以观察到：r z σσσθ≥=>0，分析采用Mises 屈服准则，表达为()()()()222222226Y z rz r z z r r στττσσσσσσθθθθ=+++-+-+- （6）该厚壁圆筒是轴对称平面应变问题，即0===θθτττz rz r ，由Mises 屈服条件其表达式可得到：Y Y r σσσσθ155.132==- （7）当压p 较小时，厚壁圆筒处于弹性状态，在a r =处，()r σσθ-有最大值，筒体由壁开始屈服，此时的压为e p ，由式（5）、（7）联立可求得弹性极限压力为()2222155.1b a b p Y e σ-= （8） 代入题目所给数据得到弹性极限压力为：()MPa p e 92.1211521015380155.1222=⨯-⨯= 2.3 弹塑性阶段当e p p <时，圆筒处于弹性状态，当e p p >的情况，在圆筒壁附近出现塑性区，产生塑性变形，随着压的增大，塑性区逐渐向外扩展，而外壁附近仍为弹性区。
Cosserat弹塑性模型在ABAQUS中的数值实施
体模型的用户子程序ｕ Ｌ Ｅ 实施方法 ｌ Ｃ ｓｅａ 连 续 体 模 型 ｏ ｓｒｔ
考虑 Ｃ ｓｅａ 连 续 体 平 面 问 题 ， 个 材 料 点 ｏ ｓｒｔ 每 有 三个 自由度 ．
Ｍ Ｌ “ Ｕ “ Ｊ
应力、 应变分 别 定义 为 ：
Ｃ ｓｅａ 弹塑性模 型在 Ａ ＡＱＵＳ中的数值实施 ｏ ｓｒｔ Ｂ
彭 董龙 吴 从文 ， 彬 ， 群。
（． １ 长江 大学 城 市建设 学 院， 北 荆 州 ４ ４ ２ ； 湖 ３ ０ ３ ２荆 州市城 市规划设 计研 究院 ， 北 荆 州 ４ ４ ０ ； 湖 ３ ０ ０ ３ 深 圳 中广核 工程设 计有 限公 司， 东 深圳 ５ ８ ３ ） ． 广 １ ０ １
有 强大 的非 线 性 计 算 能 力 、 富 的材 料 库 及 良好 丰 的扩充 功能 ， １ ９ 自 ９ ７年 进 入 我 国 以来 ， 来 越 多 越 的国 内企业 和 研 究 机 构 采 用 ＡＢ ＡＱＵＳ作 为 产 品 研 发 和科 学研 究 的工具 ． 文采 用 ＡＢ 本 ＡＱｕｓ的 用 户 接 口程序 ， 究 压 力 相 关 弹 塑 性 Ｃ ｓｅａ 连 续 研 ｏ ｓｒｔ
０ 引 言
偶应力 理 论 是 微 极 理 论 的一个 特 例 ， ｏｓｒｔ Ｃ ｓｅａ
兄弟最 先 提 出 了 完 整 的偶 应 力 理 论［ Ｔ ｕ ｉｌ ， ， ｏ ｐ ｎ Ｍｉｄ ｎ等 对 该 理 论 作 了 进 一 步 的 发 展 和完 善 ． ｎｌ ｉ Ｃ ｓｅａ 理论引 入 了旋 转 自 由度 和 相 应 的微 曲率 ， ｏｓｒｔ 引入 了与微 曲率能量共 轭的偶应力 、 以及具有“ 征 特 长度” 意义 的尺度参数． 该理论 可 以较好地 处理 网格 敏感性 和控 制方程失去椭 圆性 的问题 ， 近年来 ， 由于 细观力学 、 非均质力 学 的发展 ， ｏｓｒｔ Ｃ ｓｅａ 理论 重新受
统一弹塑性本构模型在ABAQUS中的开发与应用_潘晓明
2010年4月 Rock and Soil Mechanics Apr. 2010收稿日期：2009-05-27第一作者简介：潘晓明，男，1979年生，博士研究生，主要从事隧道及地下结构方面。
E-mail: pxm155138@文章编号：1000－7598 (2010) 04－1092－07统一弹塑性本构模型在ABAQUS 中的开发与应用潘晓明1, 2 ，孔 娟1, 2， 杨 钊1, 2， 刘 成1, 2（1.同济大学 地下建筑与工程系，上海 200092；2.同济大学 岩土及地下工程教育部重点实验室，上海 200092）摘 要：基于统一弹塑性本构模型的有限元理论格式，根据ABAQUS 的UMAT 格式要求，编制相应的接口程序，将统一弹塑性本构模型引入ABAQUS 中。
采用退化的统一强度模型（0b =时，为Mohr-Coulomb 模型）与ABAQUS 自带的Mohr-Coulomb 模型，对单轴试验和圆形硐室进行弹塑性分析，验证所开发材料子程序的正确性及高效性。
考虑到统一强度模型的一般情况（0b ≠）和屈服面硬化条件，对圆形硐室进行弹塑性计算，得到应力场的变化规律。
所提出的研究思路具有普遍性，为采用ABAQUS 平台进行本构模型的二次开发提供了借鉴和参考。
关 键 词：统一强度理论；屈服面；流动矢量；奇异点；应力拉回算法；ABAQUS 中图分类号：TD 353.6 文献标识码：ASecondary development and application of unified elastoplasticconstitutive model to ABAQUSPAN Xiao-ming 1, 2, KONG Juan 1, 2, YANG Zhao 1, 2, LIU Cheng 1,2(1.Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China;2.Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University, Shanghai 200092, China)Abstract: Based on finite element theoretical scheme of unified elastoplastic constitutive model, and according to the UMAT interface requirement of ABAQUS, the corresponding UMAT codes are programmed, which will be called by the main analytical module of ABAQUS. Adopting degenerative model of the unified strength (0b =，Mohr-Coulomb model) and the built-in Mohr-Coulomb model of ABAQUS, the uniaxial tests and circular chamber are analyzed to verify the correctness and efficiency of the developed material subroutine. Finally, considering the general situation form of unified elastoplastic constitutive model (0b ≠) and hard condition of yield surface, which are not available in ABAQUS software, circular chamber is simulated and variational discipline of stress field is obtained. The provided basic procedures and programming essentials of the UMAT redefining in ABAQUS are universal and can offer a reference for other developers.Key words: unified strength theory; yield surface; flow vector; singular points; return stress algorithm; ABAQUS1 引 言我国力学专家俞茂宏教授从多滑移单元体力学模型出发，考虑了作用在双剪应力单元体上的所有应力分量对材料屈服或破坏的不同影响，提出了一个能够适用于各种岩土类材料的统一强度理论和统一形式的数学表达式。
ABAQUS后处理二次开发在塑性成形模拟中的应用
的压下挠度矫直后 , 产品残余几何形态合格 。
参考文献 : [ 1 ] 李忠富 , 臧勇 , 王会刚1 H 型钢九辊变辊距矫直力能参数与压
弯挠度关系解析 [J ] 1 北京科技大学学报 , 2004 , (5) : 92 941 [ 2 ] 崔甫1 矫直原理与轿直机械 [M] 1 北京: 冶金工业出版社120021 [ 3 ] 王会刚 1 H 型钢矫直机理及有限元动态仿真 [ D ] 1 北京科技
关键词 : ABAQU S ; Pyt ho n ; 后处理 ; 数值模拟 ; 塑性成形 中图分类号 : TP391 文献标识码 : A 文章编号 : 100023940 ( 2006) 0420111204
Appl ication of second2developed ABAQUS post2process on numerical simulation of plastic forming
模拟的结果 , 本文提出使用 Pyt ho n 语言对 ABAQU S 后处理进行二次开发来达到这一目的 。文中探讨了二次开发 实现的原理 , 以及其中文件的读写与复制 、数据读取与处理 、结果输出与查看等关键技术 , 并以一个厚度处理的实 例来说明 。结果表明使用 Pyt ho n 对 ABAQU S 模拟产生的结果数据库进行处理 , 可以得到所要查看的厚度数据 , 从而为浏览结果以及指导后续的模拟优化提供了便利 。
基于ABAQUS二次开发方法的强化材料梁塑性极限分析
140 科技创新导报 Science and Technology Innovation Herald
工业技术
科技创新导报 2018 NO.25
Science and Technology Innovation Herald
表1 Q390E力学性能要求
屈服点σs，MPa 牌号 质量等级 厚度（直径，边长），mm
对空间问题
σσx' σx'2
' y
σ '2 y
对
[ D ] p ( H ' 9 G 3 G 2 ) σ 2 σσ xx''τ σ x' yy' σσ y'y' τσ x' zy' σ σz 'τ z'2x 'y τ x '2y 称 (4)
分析[J].淮海工学院学报:自然科学版,2004(3):18-21.
对于弹性阶段，塑性阶段分 别建 模 进 行分析，可计 算 出梁的最大位移。
(1)弹性阶段，如图2所示。 (2)塑性阶段，如图3所示。
科技创新导报 Science and Technology Innovation Herald 141
表1为Q39 0 E国家标准（GB / T1591-9 4）值，图2、3为本 次 试 轧的 Q 3 9 0 E 钢 板 拉伸和冲 击力学 性 能。表中 数 据 显 示，成品厚度规格20mm、30mm、40mm Q390E钢板的屈服 强度、抗拉强度、断后延伸率和冷弯等性能均满足国家标 准GB / T1591-9 4中Q39 0 E性能要求，特 别是 低温冲击韧性 值较高，且屈服强度等力学性能指标稳定，不随厚度增加 而明显降低。
弹塑性力学土木工程应用有限元ABAQUS分析课件
q=100Mpa
k
故应力集中因子为：
Kσφmax 279.42.794 q 100
弹塑性力学土木工程应用 有限元ABAQUS分析
误差分析
每边单元数10，最大应力288 每边单元数15，最大应力299
弹塑性力学土木工程应用 有限元ABAQUS分析
对比分析
网格划分的不同，对数据的拟合具有一定的影响， 划分的密集，计算结果更逼近理论值。
验证小孔处的应力集中系数
K σ φmax q
弹塑性力学土木工程应用 有限元ABAQUS分析
验证
基于结构和荷载的对称性，只 取结构的 1/4 进行分析。
弹塑性力学土木工程应用 有限元ABAQUS分析
验证
圆孔边缘应力最大的部位在
90°处，与理论分析的结果
一致，且最大应力279.4Mpa。
右侧施加的均布荷载为
0.09406
380
0.150
437.00
0.13976
0.13831
400
0.200
480.00
0.18232
0.18072
弹性模量E 3.00E+05
弹塑性力学土木工程应用 有限元ABAQUS分析
PART.03
有限元分析验证
弹塑性力学土木工程应用 有限元ABAQUS分析
平板圆孔应力
σρ
q 2
l0d lllnll0
lnl lnl0l
l0
l0
nom
l l0
lnl0 l0lln1nom
弹塑性力学土木工程应用 有限元ABAQUS分析
名义、真实应力（变） 真实应力与名义应力的关系
nom(1nom)
真实应变与名义应变的关系
基于ABAQUS的弹塑性本构关系模型开发
基于ABAQUS的弹塑性本构关系模型开发邱伟民【摘要】根据塑性力学基本公式的谱分解形式,采用回映算法编制用户材料子程序,数值计算结果表明算法有效.【期刊名称】《低温建筑技术》【年(卷),期】2016(038)006【总页数】3页(P47-49)【关键词】塑性力学;线性塑性硬化;谱分解;回映算法;数值一致切线刚度【作者】邱伟民【作者单位】同济大学土木工程学院建筑工程系,上海200092【正文语种】中文【中图分类】TU311.41有限元软件ABAQUS本身自带许多种材料本构关系模型，非线性功能强大，并提供多种具有良好开放性的用户子程序接口，允许用户根据自身面对的实际问题以代码形式来扩展主程序功能。
然而实际问题过于复杂，软件自带的本构关系模型并不能包括所有模型。
故文中决定利用有限元软件ABAQUS用户子程序实现塑性本构关系模型。
目前建立塑性力学模型一般采用回映算法，其基本步骤主要分为两步：(1) 先采用线弹性假定对应力应变进行试算，并将试算结果带入塑性屈服函数。
若材料未进入塑性则更新应力、应变，进入下一增量步；若已进入塑性则将试算结果作为初始值转入塑性修正步骤。
(2) 进行塑性修正，通常采用Newton-Raphson算法对事先已转化为非线性代数方程组的本构方程进行迭代直至，在误差范围内，塑性屈服函数得到满足。
更新有效应力、塑性应变、塑性硬化参数并计算弹塑性切线刚度。
回映算法是无条件稳定并且收敛于精确解的本构方程数值算法。
首先，这种方法将本构方程，根据不同的微分方程数值方法，转化为非线性代数方程的形式；然后，采用迭代方法求解非线性代数方程。
一般而言，可采用向后Euler方法将本构方程改写为非线性方程组的形式并采用经典的Newton-Raphson方法迭代求解。
对任意非线性方程f(x)=0，采用上述迭代方法可有如下迭代格式：f(x)=0⟹⟹⟹x(k+1)=x(k)+δx文中约定，上标k、k+1表示第k、k+1次迭代值；下标n、n+1表示第n、n+1次增量步对应值。
横观各向同性岩土材料应变局部化现象的有限元模拟
横观各向同性岩土材料应变局部化现象的有限元模拟常江芳;徐远杰;楚锡华【摘要】The anisotropic properties of geomaterials are significantly related to their inherent microstructures.In this paper,a modified Drucker-Prager yield criterion for transversely isotropic materials is developed by evaluating the internal friction angle with the stress state,the microstructure tensor and the material principal direction.Then,based on the Cosserat continuum theory,a consistent return mapping algorithm for the modified criterion is formulated,and a consistent tangent modulus matrix is achieved.Moreover,the codes are implemented through the user defined element subroutine (UEL) in the finite element software Abaqus,and the correctness of the program is verified by comparing the theoretical results with the relationships of the material strength to the principal direction and anisotropic degree of the material in the integration points.Finally,the influences of the principal direction and anisotropic degree of the material on the bearing capacity and the failure mode of the structure are emphatically analyzed by numerical examples,which are then compared with the results based on the classical continuum theory.It is found that the above-mentioned method is effective in simulating the strain localization of transversely isotropic geomaterials.%岩土材料的各向异性性质与其微结构紧密相关,文中考虑到内摩擦角是应力状态、组构张量及材料主方向的函数,发展了适用于横观各向同性岩土材料的修正的Drucker-Prager屈服准则.基于Cosserat连续体理论,推导了该准则的一致性映射返回算法,形成了一致性切线模量矩阵,并利用有限元软件Abaqus的用户单元子程序(UEL)进行了数值实现.通过将积分点材料强度随材料主方向及各向异性程度的变化关系与理论结果进行比较,验证了程序开发的正确性.数值算例重点分析了材料主方向和各向异性参数对结构极限承载力及破坏模式的影响,且与经典连续体的结果进行了比较.结果表明,文中方法能够较好地模拟具有横观各向同性的岩土材料的应变局部化现象.【期刊名称】《华南理工大学学报（自然科学版）》【年(卷),期】2016(044)010【总页数】11页(P70-80)【关键词】岩土材料;岩土力学;横观各向同性;Cosserat连续体;屈服准则;应变局部化【作者】常江芳;徐远杰;楚锡华【作者单位】武汉大学工程力学系,湖北武汉430072;石家庄铁道大学工二程力学系,河北石家庄050043;武汉大学工程力学系,湖北武汉430072;武汉大学工程力学系,湖北武汉430072【正文语种】中文【中图分类】O34自然界的岩土材料一般均呈现很强的各向异性.各向异性又分为固有各向异性和诱导各向异性.由于其沉积过程、孔隙、裂纹等固有属性而使其在原位状态下表现出的各向异性称为固有各向异性；岩土材料在非比例加载或塑性变形的影响下，由于主应力方向不断变化，导致其细观结构发生演化而表现出的各向异性称为诱导各向异性.关于此两种各向异性，国内外学者做了相关的理论及实验研究[1- 3]，提出了多种各向异性屈服准则[4- 7].姚仰平等[4]以横观各向同性土为例，阐述了同时考虑材料外部应力分布与材料内部强度分布下各向异性土的破坏机制.贾乃文[5]以Hill正交异性屈服准则为基础，讨论了横观各向同性条件下，岩土中出现塑性屈服的Drucker准则修正表达式.Duveau等[6] 对一些强各向异性材料的屈服准则做了一个综合评价.Gao等[7]在屈服函数的摩擦系数中引入了一个各向异性变量，考虑它为应力不变量及微结构张量的联合不变量的函数，提出了一个适用于横观各向同性材料的屈服准则.然而基于唯象方法提出的各向异性屈服准则一般包含较多的材料参数，且这些参数与材料微观结构之间的关系不明确.Pietruszczak等[8- 9]提出了一个各向异性屈服准则，基于物质点给出了各向异性程度及加载方向对材料强度影响的理论分析，其材料参数是应力和微结构张量联合不变量的显式函数，使得实验验证更易实现.余天堂等[10- 11]基于Pietruszczak等[8- 9]提出的各向异性屈服准则，考虑各向异性参数和单轴抗压强度是一个由微结构张量和加载方向表示的分布函数，给出了沉积岩的一种各向异性模型.Zhong等[12]对Pietruszczak-Morz 屈服准则中的材料参数进行了分析.在Pietruszczak等[8- 9]的研究基础上，Lade[13- 15]推导了一个基于 Lade 强度准则的各向异性强度准则，并展开了相应的理论及实验研究.然而上述研究多关注的是对各向异性破坏准则的描述，是各向异性对一个物质点破坏的影响.关于各向异性对结构承载力及破坏模式影响的数值模拟仍相对较少.Chu 等[16]在Pietruszczak等[8- 9]提出的各向异性屈服准则的基础上，考虑内摩擦角为组构张量、应力状态及材料主方向的函数，对Drucker-Prager屈服准则进行了修正，以横观各向同性材料为例研究了各向异性对应变局部化的影响，然而不足之处是其基于经典连续体理论，在伴随应变软化出现应变局部化现象的同时会产生网格依赖性.为克服这一问题，文中基于Cosserat连续体理论[17- 19]，推导了修正的Drucker-Prager屈服准则的一致性映射返回算法及切线模量矩阵，重点分析了材料主方向和各向异性程度对横观各向同性材料极限承载力及破坏模式的影响，最后通过与经典连续体计算结果进行比较来验证解决网格依赖性问题的有效性.1 基于Cosserat连续体的横观各向同性弹塑性模型1.1 Cosserat连续体的基本方程与经典连续体相比，Cosserat连续体中引入了独立的转动自由度，平面应变情况下，其位移矢量u可表示为[17]u=[ux uy ωz]T其中，ux、uy为平动位移，ωz为独立的转动位移.应变矢量ε除了与经典意义上所对应的正应变和剪应变有关以外，还包含两个微曲率，如下所示：ε=[εx εy εz εxy εyx lcκxz lcκyz]T其中：可以看到这里引入了一个内部长度参数lc，保证了Cosserat连续体的控制方程在计算过程中的正定性，可有效克服网格依赖性问题.几何方程可表示如下：ε=Lu其中，Cosserat连续体的应力矢量σ引入了两个偶应力mxz和myz(如图1所示)：弹性应力-应变关系可以表示为σ=Deε其中本构矩阵为这里,和G是传统意义上的拉梅常数，Gc则为Cosserat剪切模量.准静态条件下，若忽略体积力及体力偶的影响，则平衡方程可表示为LTσ=01.2 修正的Drucker-Prager准则适用于各向同性材料的Drucker-Prager准则以应力不变量的形式表示为F=q+Aφp+Bφ=0其中,p和q可分别用应力第一不变量及偏应力第二不变量表示如下:其中，P矩阵可表示为式(10)中参数Aφ、Bφ为内摩擦角φ和黏聚力c的函数，若Drucker-Prager准则为Mohr-Coulomb准则的内接圆，则可将其分别表示为通常对黏聚力考虑线性硬化或软化，即其中，c0为初始黏聚力，为硬化(软化)参数，p为等效塑性应变.基于应力不变量与组构张量，Pietruszczak与Mroz给出了一个能够描述各向异性材料的屈服准则[8]：η=η0(1+Ωijlilj)其中：η为描述屈服函数的材料参数，与应力状态有关；η0为η的平均值；Ωij 为描述η偏离η0程度的偏量.对正交各向异性材料Ωij有两个不等的特征值，对横观各向同性材料，其可用一个标量描述，对各向同性材料，Ωij为零，此时η=η0.li、lj为相对于材料主轴的加载方向.式(16)也可用Ωij的主值表示为对横观各向同性材料按如下方式建立坐标系，以竖直方向为y轴，以水平面为x-z 平面，建立整体坐标系Oxyz.令局部坐标系与材料主方向共轴，即各向同性面为1-3面，其法向为e(2)方向，其中方向e(3)与z轴重合，如图2所示.此时，Ω1=Ω3，且Ω1+Ω2+Ω3=0以及=1，式(17)可进一步表示为加载向量按如下方式计算：式中,表示应力张量，N=e(2)⊗e(2).为了简化，仅考虑内摩擦角的各向异性，即令则有式(10)、式(3)、式(1)构成了描述横观各向同性材料的屈服函数.对于岩土材料，通常采用非关联流动法则，塑性势可取为G=q+Aψp+Bψ类似地：其中，ψ为膨胀角.由一致性算法最终得到更新后的应力[17- 18]，σ=CαsE+(pE-KAψΔ)m其中：Δ其中FE=qE+AφpE+BE.以及一致性弹塑性切线模量矩阵：其中:可以验证，当采用关联流动法则时，该矩阵具有主对称性.以Abaqus软件的UEL子程序接口为平台，对文中发展的单元类型和弹塑性本构模型进行了数值实现[20].以平板压缩模型为例，单元类型为平面应变8节点减缩单元，模型尺寸为0.6 m×0.8 m，上下边界约束x方向自由度，左右边界自由，上下边界施加对称均布位移荷载uy，如图3(a)所示.计算中所采用的材料参数如表1所示.1)E为弹性模量；υ为泊松比；β为材料主方向角.Pietruszczak等[8]基于他们提出的各向异性屈服准则对横观各向同性材料物质点的理论分析显示，对于不同的加载方向，材料轴向强度随Ω1的演化规律不同，存在一个转折点，该点材料的轴向强度将独立于各向异性程度，且Zhong等[12]证明了此转折点所对应的材料主方向与材料固有属性有关.这里，为与已有理论结果进行对比，选取了结构的两个代表性单元，标号为Ⅰ和Ⅱ，分别位于板上边缘中间位置和剪切带附近，如图4所示，给出了两个单元1号积分点(见图3(b))上等效应力随加载位移的变化曲线，其中参数取β=30°，Ω1=-0.1.易知峰值qmax代表了该点的强度，图5(a)和5(b)则给出了积分点(物质点)上归一化的材料强度随各向异性程度及材料主方向的变化趋势，可见其趋势与已有理论结果[8]较为吻合(见图5(c),其中法向强度ζ=Δσ1/(2σ0))，从而验证了本程序开发的正确性.下面主要就材料主方向β及各向异性程度Ω1对结构极限承载力及破坏模式的影响做详细分析. 2.1 材料主方向β的影响由图2可知，β为各向同性面相对x轴逆时针转动的角度，描述了材料主轴相对整体坐标系的偏移，β的改变也可理解为加载方向的改变.图6为Ω1=-0.10，加载位移为0.05 m时等效塑性应变分布云图.值得注意的是，Ω1取负值意味着各向同性面法向(e(2)方向)有较大的摩擦强度[14].可以看出剪切带呈现了X型分布，当β=0°和β=90°时，由于加载方向和材料主方向保持一致或垂直，剪切带呈对称分布，如图6(a)和6(e)所示.当β≠0°和β≠90°时，剪切带两个分支呈现不对称性，沿各向同性面方向的剪切带较宽，等效塑性应变值较高，文中称之为强剪切带，与之共轭的另一条剪切带称为弱剪切带，如图6(b)和6(h)所示，与姚仰平等[4]的理论分析相比，其破坏面的位置可能就是强剪切带出现的位置，数值模拟显示还存在一个弱剪切带，只是在破坏时强剪切带占据了主导地位.由图6(b)-6(d)可见，当0°<β<90°时剪切带呈“/”型分布，且随着β的增大“/”方向剪切带逐渐变宽.反之，当90°<β<180°时，剪切带呈“\”型分布，且随着β的增大“\”方向剪切带逐渐变窄.对于一对互补的β角，如β=30°和β=150°，等效塑性应变峰值相等，但剪切带分布模式相反.随着β从0°演化到180°，强剪切带上等效塑性应变峰值出现了增大-降低-增大-降低的对称过程，弱剪切带则具有互补的变化过程，如图7所示.图8给出了随β变化材料的承载力-位移曲线,可以清楚地看到当各向同性面位于水平位置，即β=0°时，结构承载力最大，随着β的增大承载力逐渐降低，90°时最低，这与横观各向同性材料屈服准则中摩擦强度参数各向异性的分布规律一致. 2.2 各向异性参数Ω1的影响如前所述，Ωij为描述η偏离η0程度的偏张量，对于横观各向同性材料可以用一个标量Ω1来描述.需指出文中各向异性参数取值来自于文献[8]，依次取Ω1=0，-0.05，-0.1,-0.15,β=30°,加载位移为0.05 m，其他参数同表1.图9给出了等效塑性应变随各向异性程度变化的结果.可以看出，各向异性程度越大，X型剪切带的非对称性越强，强剪切带等效塑性应变的峰值越大，弱剪切带峰值则越小.图10为承载力-位移曲线，随着Ω1绝对值的增大结构承载力逐渐增大.这里同时给出了归一化的承载力随各向异性程度及材料主方向的综合变化规律，如图11所示，其中，纵轴表示各承载力峰值除以Ω1=0时的承载力峰值，可见结构整体响应亦呈现出了与图5类似的规律，即在β约为47°的位置出现转折点，β小于该值时，结构承载力随Ω1绝对值的增大而增大，大于该值时则相反.2.3 网格依赖性调查基于经典连续体理论，采用有限单元法计算弹塑性问题，当应变软化和局部化现象发生时，不可避免地会出现网格依赖性问题.本节通过对经典连续体和Cosserat连续体两种理论框架下的计算结果进行比较，调查了有限元计算的网格依赖性问题. 取材料主方向β=90°，Ω1=-0.1，其他材料参数同表1，网格密度分别取6×8，12×16，18×24，24×32，30×40，60×80进行比较.图12给出了基于经典连续体等效塑性应变分布图，可以看到随着网格加密，等效塑性应变的峰值逐渐增大，且剪切带宽度明显变窄.图13显示结构极限承载力随网格密度增大而逐渐降低，且模拟软化段的能力亦逐渐降低.图14则为基于Cosserat连续体得到的等效塑性应变分布图，随着网格加密，等效塑性应变峰值有所增大但剪切带宽度基本保持不变.图15显示虽然结构承载力有些许降低，但其变化越来越小，计算结果最终收敛于较为接近的值.需说明的是，图中网格为6×8时计算结果与其他结果相比差异较大，这是由于网格稀疏到一定程度，有限元计算本身的误差所致.经过比较可见引入Cosserat连续体后，网格依赖性问题得到了有效的解决.基于Pietruszczak 和 Mroz提出的各向异性屈服准则，考虑摩擦角为材料主方向及应力张量与微结构张量联合不变量的函数，提出了一个适用于横观各向同性岩土材料的修正的Drucker-Prager屈服准则，并基于Cosserat连续体推导了其一致性映射返回算法和切线模量矩阵.数值算例重点分析了材料主方向及各向异性程度对结构极限承载力及破坏模式的影响，结果表明：(1)材料主方向β=0°或β=90°时剪切带呈对称分布，前者结构承载力最大，后者最低.0°<β<90°时，随β增大，剪切带逐渐变宽，结构承载力逐渐降低，90°<β<180°则相反.β从0°到180°，强剪切带等效塑性应变峰值经历了增大-降低-增大-降低的对称过程，弱剪切带则相反.β互补时剪切带呈现相反的分布模式. (2)Ω1绝对值越大，材料各向同性面的强度越低；对于β≠0°、β≠90°的非对称情况，剪切带呈现的X型分布的非对称性更强.承载力变化趋势存在一个转折点，约为β=47°，β小于该值时，结构承载力随Ω1绝对值的增大而增大，β大于该值时则相反.(3)与经典连续体的模拟结果比较表明，文中基于Cosserat连续体的数值方法较好地解决了网格依赖性问题.† 通信作者: 徐远杰(1956-)，男，博士，教授，主要从事工程结构破坏数值仿真研究.E-mail:***************【相关文献】[1] ARTHUR J R F,Menzies B.Inherent anisotropy in a sand[J].Geotechnique,1972,22(1):115- 128.[2] ARTHUR J R F,CHUA K S,DUNSTAN T.Induced anisotropy in a sand[J].Geotechnique,1979,27(1):13- 30.[3] ODA M.Inherent and induced anisotropy in plasticity theory of granular soils[J].Mechanics of Materials,1993,16(1/2):35- 45.[4] 姚仰平,祝恩阳.横观各向同性土的简明破坏机制解释 [J].岩土力学,2014,35(2):328- 333. YAO Yang-ping,ZHU En-yang.Concise interpretation of damage mechanism for cross-anisotropic soil [J].Journal of Geotechnical Mechanics,2014,35(2):328- 333.[5] 贾乃文.岩土工程中横观各向同性Drucker屈服准则 [J].华南理工大学学报(自然科学版)，1993,21(2):49- 54. JIA Nai-wen.Drucker yield criterion of horizontally isotropic rock and soil [J].Journal of South China University of Technology(Natural ScienceEdition),1993,21(2):49- 54.[6] DUVEAU G,SHAO J F,HENRY J P.Assessment of some failure criteria for strongly anisotropic materials [J].Mechanics of Cohesive-Frictional Materials,1998,3(1):1- 26.[7] GAO Z W,ZHAO J D,YAO Y P.A generalized anisotropic failure criterion for geomaterials [J].International Journal of Solids and Structures,2010,47(22/23):3166- 3185.[8] PIETRUSZCZAK S,MROZ Z.Formulation of anisotropic failure criteria incorporating a microstructure tensor [J].Computer and Geotechnics,2000,26(2):105- 112.[9] PIETRUSZCZAK S,MROZ Z.On failure criteria for anisotropic cohesive-frictional materials [J].International Journal for Numerical and Analytical Methods in Geomechanics,2001,25(5):509- 524.[10] 余天堂,卢应发,SHAO J F,等.沉积岩的一种各向异性模型 [J].岩土力学,2002,23(1):47- 50. YU Tian-tang,LU Ying-fa,SHAO J F,et al.An anisotropic model of sedimentary rocks [J].Journal of Geotechnical Mechanics,2002,23(1):47- 50.[11] 余天堂.岩土材料固有各向异性的模拟 [J].岩石力学与工程学报,2004,23(10):1604- 1607. YU Tian-tang.Modeling of inherent anisotropy for geotechnical material [J].Chinese Journal of Rock Mechanics and Engineering,2004,23(10):1604- 1607.[12] ZHONG S Y,XU W Y,LING D S.Influence of the parameters in the Pietruszczak-Mroz anisotropic failure criterion [J].International Journal of Rock Mechanics and Mining Sciences,2011,48(6):1034- 1037.[13] LADE P V.Modeling failure in cross-anisotropic frictional materials [J].International Journal of Solids and Structures,2007,44(16):5146- 5162.[14] LADE P V.Failure criterion for cross-anisotropic soils [J].Journal of Geotechnical and Geoenvironmental Engineering,2008,134(1):117- 124.[15] LADE P V.Three-dimensional failure in geomaterials:experimentation and modeling [M]∥Constitutive Modeling of Geomaterials.[S.l.]:Springer Berl in Heidelberg,2013:47- 58.[16] CHU X H,CHANG J F,JIANG Q H,et al.Numerical simulation of progressive failure in transversely isotropic geomaterials[C]∥Advance in Heterogeneous Material ncaster:Destech Publication Inc.,2011:965- 968.[17] DE BORST R.Simulation of strain localization:a reappraisal of the Cosserat continuum [J].Engineering Computations,1991,8(4):317- 332.[18] LI X K,TANG H X.A consistent return mapping algorithm for pressure-dependent elastoplastic Cosserat continua and modeling of strain localization [J].Computers & Structures,2005,83(1):1- 10.[19] 余村,楚锡华,唐洪祥,等.基于Cosserat连续体的颗粒破碎影响研究 [J].岩土力学,2013,34(增刊1):67- 79. YU Cun,CHU Xi-hua,TANG Hong-xiang,et al.Study of effect of particle breakage based on Cosserat continuum [J].Journal of Geotechnical Mechanics,2013,34(Sup 1):67- 79.[20] ARSLAN H,STURE S.Finite element analysis of localization and micro-macro structure relation in granular materials(Part I):Formulation [J].Acta Mechanica,2008,197(3/4):135- 152.。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
Cosserat弹塑性模型在ABAQUS中的数值实施彭从文;董龙彬;吴群【摘要】The subroutine UEL embedded in commercial FEM software ABAQUS is utilized to develop a user element based on the material of pressure-dependent elasto-plastic Cosserat continuum model. The user element adopt plane eight-node isoparametric unit which include three degree of freedom (two translation and one rotation) at each node. The Drucker-Prager material model with associate flow rule is used, and the stress integration algorithm is four-order Runge-Kutta method. The user element is used to analyze the influence of mesh density and the material characteristic length in the material localization problem. The results show that the thickness of shear band and the equivalent plastic strain are independent of mesh density, with the increase of the material characteristic length, the thickness of shear band and the equivalent plastic strain decrease constantly.%利用大型有限元软件ABAQUS提供的接口程序UEL,开发了压力相关弹塑性Cosserat连续体材料的用户单元.采用平面八节点等参单元,包括平动与转动三个自由度,考虑相关联流动的Drycjer-Prager材料模型,应力计算采用四阶龙格-库塔法.利用该单元分析了材料局部化问题中网格密度与材料特征长度的影响.结果表明,网格密度对材料剪切带厚度与等效塑性应影响很小；随着特征长度增大,剪切带厚度增大,等效塑性应变峰值减小.【期刊名称】《武汉工程大学学报》【年(卷),期】2011(033)006【总页数】5页(P102-106)【关键词】Cosserat模型；FEM；ABAQUS；UEL；局部化【作者】彭从文;董龙彬;吴群【作者单位】长江大学城市建设学院,湖北荆州434023;荆州市城市规划设计研究院,湖北荆州434000;深圳中广核工程设计有限公司,广东深圳518031【正文语种】中文【中图分类】TB120 引言偶应力理论是微极理论的一个特例，Cosserat兄弟最先提出了完整的偶应力理论[1]，Toupin[2]，Mindlin等[3]对该理论作了进一步的发展和完善.Cosserat理论引入了旋转自由度和相应的微曲率，引入了与微曲率能量共轭的偶应力、以及具有“特征长度”意义的尺度参数.该理论可以较好地处理网格敏感性和控制方程失去椭圆性的问题，近年来，由于细观力学、非均质力学的发展，Cosserat理论重新受到关注，逐渐成为研究热点之一.数值方法是重要的研究手段之一，为了提高计算精度与效率，基于大型通用数值计算平台的二次开发方法得到了广泛的应用.目前,许多数值计算平台没有内嵌Cosserat计算模型，关于Cosserat模型的二次开发还不多见[4-6]. ABAQUS是目前最流行、功能最强的商用有限元软件之一，该软件可以进行结构静、动力分析，具有强大的非线性计算能力、丰富的材料库及良好的扩充功能，自1997年进入我国以来，越来越多的国内企业和研究机构采用ABAQUS作为产品研发和科学研究的工具.本文采用ABAQUS的用户接口程序，研究压力相关弹塑性Cosserat连续体模型的用户子程序UEL实施方法.1 Cosserat连续体模型考虑Cosserat连续体平面问题，每个材料点有三个自由度.u=应力、应变分别定义为：S=E=几何方程为：ui,j=uj,i-eijkωk(1)κij=ωj,i(2)静力平衡方程为：σij,j+fi=0(3)mij,i+eijkσik+qj=0(4)式(1)～(4)中，fi、qj分别为体积力与体积力偶；eijk为排列算子；ux,uy,ωz分别是平面内平移与转动自由度；mxz,myz偶应力；κxz,κyz为微曲率.对于弹性材料，其本构关系为[7]S=DeE(5)式(5)中，G,v,a,l分别是材料的剪切模量、泊松比、COSSERAT材料参数及特征长度.对于弹塑性Cosserat材料，采用基于Drucker-Prager屈服准则的弹塑性Cosserate连续体模型，其屈服函数与流动势函数分别为[8](6)(7)式(7)中,其中，α,β,A,B为材料参数，σ,s分别为应力与偏应力张量.Cosserat连续体弹塑性本构关系推导方法同经典连续介质力学，应力应变增量关系为(8)式(8)中，内变量采用等效塑性应变,为硬化项定义为2 UEL实现方法及计算流程2.1 子程序编写注意事项(1)ABAQUS提拱了两类单元自定义方法.一类是线性单元，可以通过结果文件或INP文件直接给出，不需要编写UEL；另一类就是通用单元，通过UEL子程序定义；(2)UEL子程序中更新变量与分析问题类别有关.同一个模型中可能遇到不同的分析步，如地应力平衡、静力分析、摄动步分析等，因此，编写UEL时要区别处理；(3)UEL允许自定义荷载.包括集中荷载、均布荷载及弯矩等.其中，对于均布荷载，须定义荷载标志号；(4)自定义单元在ABAQUS/CAE中不可见.若想在ABAQUS/CAE中显示自定义单元变形图，可以将ABAQUS标准单元与自定义单元绑定，同时将标准单元材料参数设为小值.UEL子程序中所有输出变量均通过SDV写入结果文件(.fil、.dat)，其分量在ABAQUS/CAE中不可见.2.2 AMATRX与RHS计算UEL界面与ABAQUS内核主要通过AMATRX与RHS等变量进行数值传递.本文采用八节点等参单元，设单元节点位移、插值函数与位移应变转换矩阵分别为d、N与B，单元内任一点位移u及应变ε为u=∑Nidi=Nd(9)ε=Bd(10)式(9～10)中,B=将式(9)、(10)代入Cosserat介质虚功方程(11)进行方程离散.(11)式(11)中，f=，t=对于线性问题，结合材料本构关系，得到式(12).(12)式(12)中，将式(12)简记为Kd=f′(13)对于非线性问题，采用Newton-Raphson方法，将式(13)改写为ψ(d)=K(d)d-f′(14)设ψ(d)为具有一阶导数的连续导数，初始近似值为d(0)，第i次迭代的近似值为d(i).将函数ψ(d)在d(i)处展开，保留线性项，忽略高阶项得:(15)d(i+1)=d(i)+Δd(i)(16)式(15)中，与Ψ(d)分别为UEL中需更新变量AMATRX与RHS.2.3 计算流程计算流程如图1所示.图1 UEL流程图Fig.1 Flow of UEL3 性状分析3.1 有限元模型模型几何尺寸10×5 m2，采用三种不同网格密度，单元数分别是15×30、20×40、25×50.边界条件为底端竖向固定，左侧水平向固定.材料弹性模量25 GPa，泊松比0.3，内摩擦角35°，粘聚力1.5 MPa，软化模量15 MPa.采用相关联流动法则.顶部采用位移加载方式，加载量为20 mm，加载方向向下.为了触发局剪切带，对左下角单元弱化处理.采用高斯完全积分，四阶龙格-库塔显式应力积分方法.3.2 计算结果计算模型分析了不同网格密度及特征长度的影响，计算结果如图2～6所示.(1)局部化带的客观性.图2为经典连续介质理论得到的等效塑性应变云图，图3与图4分别是采用Cosserat理论计算得到的等效塑性应变云图与应力应变曲线.由图可知，采用Cosserat理论计算时，随着网格密度增加，剪切带厚度与等效塑性应变峰值基本不变.当采用经典连续介质理论计算时，计算结果有明显的网格依赖性，随着网格密度增加，软化带逐渐变窄，等效塑性应变峰值也不断增大,计算收敛趋于弱化.图2 不同网格密度等效塑性应变云图Fig.2 Contour of equivalent plastic strain with different mesh density注：(a)15×30, (b)20×40,(c)25×50图3 不同网格密度等效塑性应变云图(特征长度0.15)Fig.3 Contour of equivalent plastic strain with different mesh density (characteristic length:0.15)注：(a)15×30, (b)20×40,(c)25×50图4 不同网格密度下应力位移曲线Fig.4 Stress-displacement curve with different mesh density(2)特征长度的影响.Cosserat理论引入特征长度作为正则化机制，特征长度决定Cosserat连续体模型模拟应变局部化问题的能力并影响局部化剪切带宽度大小.图5与图6分别为不同特征长度下采用Cosserat理论计算得到的等效塑性应变云图和应力位移曲线.由图2～6可知，随着特征长度增大，剪切带厚度增大，等效塑性应变峰值减小，材料软化模量降低.图5 不同特征长度下等效塑性应变云图(网格20×40)Fig.5 contour of equivalent plastic strain with different characteristic length (mesh density: 20×40)注：(a)0.2 (b)0.15 (c)0.05图6 不同特征长度的影响Fig.6 Effect of different characteristic length4 结语基于ABAQUS接口程序UEL，开发了压力相关弹塑性Cosserat连续体材料的用户单元，并采用该单元分析了有限元网格密度及材料特征长度对材料局部化的影响.结果表明，采用Cosserat理论计算时网格密度对材料剪切带厚度、等效塑性应变影响很小，这也在一定程度上说明本文方法的正确性.要特别说明的是，基于ABAQUS平台进行二次开发能有效地利用现有程序代码，减小开发工作量，缩短有限元程序开发周期，极大地提高科研工作效率.参考文献：[1] Cosserat E, Cosserat F. Theorie des Corps Deformables [M].Paris: Herman et Files, 1909.[2] Toupin R A. Elastic materials with couple stresses [J].Archive Rational Mechanics and analysis, 1962(11): 385-414.[3] Mindlin R D, Tiersten H F. Effects of couple stresses in linear elasticity [J]. Archive Rational Mechanics and analysis, 1962(11): 415-448.[4] 杨乐, 吴德伦, 许年春. 偶应力理论的层状岩体洞室数值模拟[J]. 重庆建筑大学学报, 2008, 30(3):73-77.[5] 尹雪英, 杨春和, 李银平. 层状盐岩体三维Cosserat介质扩展本构模型的程序实现[J]. 岩土力学, 2007,28(7):1415-1420,1426.[6] 朱珍德, 秦天昊, 王士宏, 等. 基于Cosserat理论的柱状节理岩体各向异性本构模型研究[J]. 岩石力学与工程学报,2010,29(增2): 4068-4076.[7] Sharbati E, Naghdabadi R. computational aspects of the cosserat finite element analysis of localization phenomenon[J]. Computational materials science,2006(38):303-315.[8] HAYDAR ARSLAN. Localization analysis of granular materials in cosserat elastoplasticity -formulation and finite element Implementation [D]. Colorado: University of Colorado, 2006, 85-88.。