【客观应力率】Abaqus折叠屏材料弯曲模拟

可折叠显示设备日益走进我们的生活，对此类屏幕分析验证是当今CAE工程师面临的难题之一，因为必须要考虑多层堆叠的复合材料，并进行90度弯曲、展开的大变形模拟；另外，为了预测它的耐用性，需要确定以何种损伤标准进行评估。执行这种高度复杂的显示器分析，先决条件是进行精确的应力、应变计算，在此之前，工程师必须要了解一个基本概念，那就是“客观应力率”。

01

小张的困惑：线弹性材料的“残余应力”！

小张是训练有素的CAE工程师，有一天，他接到一个分析任务：折叠屏材料的弯曲有限元分析，心想，还真是赶时髦呀，来吧。

供应商提供了某一层材料的试验数据曲线，笔直的让人能口算出弹性模量，试验部门也提前告知了弯曲试验完全在此应变范围进行加、卸载。于是，小张确信用线弹性本构无疑，一顿操作，下班前竟完成了弯曲试验对标：仿真得出来的应力、应变、弯矩和试验结果完全一致。

折叠屏某层材料90°弯曲仿真-加载

正要高兴的时候，他看到了卸载的计算结果：

卸载后的应力、应变

线弹性材料加、卸载怎么会出现“残余应力”？于是他又校核了一下模型：线弹性材料模型、静力学分析，几何非线性，ALM接触、沙漏控制，一切都很合理，否则前面的试验对标不会这么顺利，然而并没有定义塑性啊，为什么材料会表现出如此强烈的路径依赖性？

就算是数值误差，也不可能在这个量级的吧？这一度成为小张童鞋的一个困惑。

02

什么是客观应力率？

客观应力率（Objective Stress Rate）是连续介质力学的一个基本概念，为什么会有这么个概念呢？考虑一根受力平衡、两端拉力恒为F的杆，在空间中发生一定角度的刚体转动（地面上的观察者），即使力的大小不变，方向也是不断变化的，所以杆中（柯西）应力变化率不为0；另一方面，杆的长度始终不变，应变和初始状态一样，是恒定的，应变的变化率为0，若按照不变的本构描述（杆上的观察者），可以立刻得出应力变化率为0！换句话说，柯西应力的应力率不是客观的，从下图也可以看出，柯西应力满足张量旋转的协变规则，而应力率不满足。

柯西应力不满足应力率的客观性

客观应力率是应力对时间的导数，不应该依赖于参考系，因为很多本构方程是通过应力（率）和应变（率）之间的关系进行描述的，材料的机械响应和坐标系的选择没有关系，材料本构方程也应该是和框架无关的（客观的），总不能为每一个观察者配一套材料本构方程吧？于是力学家们提出了多种满足客观性的应力率，在各大商业有限元软件中应用比较广泛的有：

• 柯西应力的Truesdell率

• 柯西应力的Green-Naghdi率

• 柯西应力的Jaumann率

Abauqs用的是后两个，通过Nlgeom=ON启用，具体是哪一种和模型的设定有关，取决于单元类型、求解器以及材料本构模型。

Abauqs中的客观应力率

注：

Solid指连续单元

Structural指梁、杆、膜、壳这类结构单元

Nlgeom=ON，Standard+Solid Element，对应就是Jaumann率，这种情况下，应变增量是通过应变率在时间上积分得到的，而应变率的更新又依赖于上一个增量的即时构型坐标系，所以，对于小张的模型，总是存在初始增量应变率为0，最后一个增量应变率不为0的情况，那么闭环积分（加卸载）的时候，肯定会引起所谓的“残余应力”和“残余应变”。

问题找到了，有两种方法可以改善:

1. 使用非常小的增量（Step里通过Maximum Increment Size控制），计算出的“残余应变”基本会和Maximum Increment Size保持一个量级，比如你觉得应变结果在1E-8量级上是可以接受的，那么你就得设置Maximum Increment Size=1E-8，显然，对于很多分析，这样做是不太现实的。
   

   
   

2. 使用超弹性本构，超弹性本构基于应变能势描述材料，自动确保能量守恒，它使用总应变公式，不用更新坐标系，因此从定义上来说就是客观的，这是它的先天优势。

小张尝试把线弹性改为超弹性的Mooney-Rivlin或者Neo-Hookean之后，得到的应力、应变、弯矩结果和线弹性之间的误差只有1%左右，而且“残余应力”和“残余应变”都降低到了可以忽略的程度。

线弹性-超弹性应力结果差异（加载）

线弹性-超弹性应变结果差异（加载）

线弹性-超弹性弯矩结果差异（加载）

线弹性-超弹性“残余应变”结果差异（卸载）

问题终于解决了~

03

如何将线弹性材料“超弹性化”？

看到这里，小张的困惑是解决了，有的朋友可能困惑了：这家伙是怎么把线弹性材料“超弹性化”的？

我们要考虑的应变相对来说比较小，大致在（-0.1~0.1）范围内，通过下面的公式即可完成转换，反过来（超弹性本构线弹性化）也是可以的。

小范围应变下弹性模量与超弹性参数关系式

比如：

线弹性E=200MPa → 超弹性NH，C1=33.33MPa

线弹性E=200MPa → 超弹性MR，(C10 , C01 )=(35MPa, -1.667MPa)

通过绘制这三种本构模型的应力-应变曲线可以看出，在（-0.1~0.1）应变范围内，三者是基本重合的，可以用作互相替换。

三种本构的应力-应变曲线对比

参考：

SIMULIA User Assistance 2022, Dassault Systèmes Simulia Corp., Providence, RI, USA.

Wikipedia-Objective stress rate

Nonlinear Finite Elements for Continua and Structures, Belytschko, T., Liu, W.K. and Moran, B., John Wiley and Sons, p. 99, 2000.

Finite Element Procedures, Bathe, K.-J., Prentice-Hall, p. 612, 1996.

Stress rate objectivity and foldable display simulation, Naver blog 2022, korea

'Finite Deformation Constitutive Equations and a Time Integration Procedure for Isotropic, Hyperelastic-Viscoplastic Solids,' Weber, G. and Anand, L., Computer Methods in Applied Mechanics and Engineering, Vol. 79, pp. 173-202, 1990.