Cosserat弹塑性模型在ABAQUS中的数值实施

彭从文，董龙彬，吴 群

(1.长江大学 城市建设学院，湖北 荆州 434023; 2 荆州市城市规划设计研究院，湖北 荆州 434000；3.深圳中广核工程设计有限公司，广东 深圳 518031)

0 引 言

偶应力理论是微极理论的一个特例，Cosserat兄弟最先提出了完整的偶应力理论[1]，Toupin[2]，Mindlin等[3]对该理论作了进一步的发展和完善.Cosserat理论引入了旋转自由度和相应的微曲率，引入了与微曲率能量共轭的偶应力、以及具有“特征长度”意义的尺度参数.该理论可以较好地处理网格敏感性和控制方程失去椭圆性的问题，近年来，由于细观力学、非均质力学的发展，Cosserat理论重新受到关注，逐渐成为研究热点之一.

数值方法是重要的研究手段之一，为了提高计算精度与效率，基于大型通用数值计算平台的二次开发方法得到了广泛的应用.目前,许多数值计算平台没有内嵌Cosserat计算模型，关于Cosserat模型的二次开发还不多见[4-6]. ABAQUS是目前最流行、功能最强的商用有限元软件之一，该软件可以进行结构静、动力分析，具有强大的非线性计算能力、丰富的材料库及良好的扩充功能，自1997年进入我国以来，越来越多的国内企业和研究机构采用ABAQUS作为产品研发和科学研究的工具.本文采用ABAQUS的用户接口程序，研究压力相关弹塑性Cosserat连续 体模型的用户子程序UEL实施方法.

1 Cosserat连续体模型

考虑Cosserat连续体平面问题，每个材料点有三个自由度.

应力、应变分别定义为：

几何方程为：

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)中,

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)

B=

将式(9)、(10)代入Cosserat介质虚功方程(11)进行方程离散.

(11)

对于线性问题，结合材料本构关系，得到式(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)

2.3 计算流程

计算流程如图1所示.

图1 UEL流程图

3 性状分析

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 不同网格密度等效塑性应变云图

图3 不同网格密度等效塑性应变云图(特征长度0.15)

图4 不同网格密度下应力位移曲线

(2)特征长度的影响.Cosserat理论引入特征长度作为正则化机制，特征长度决定Cosserat连续体模型模拟应变局部化问题的能力并影响局部化剪切带宽度大小.图5与图6分别为不同特征长度下采用Cosserat理论计算得到的等效塑性应变云图和应力位移曲线.由图2～6可知，随着特征长度增大，剪切带厚度增大，等效塑性应变峰值减小，材料软化模量降低.

图5 不同特征长度下等效塑性应变云图(网格20×40)

图6 不同特征长度的影响

4 结 语

基于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.