混凝土材料umat子程序报告

第 1 页

umat子程序报告

这一阶段主要目的是针对于应力——应变曲线问题进行研究，由于本用户子程序是基于abaqus有限元软件来模拟混凝土材料的弹塑性损伤行为，而众所周知混凝土材料在单轴压缩的情况下的应力——应变曲线具有以下三个较为明显的阶段（见图1）

1．弹性阶段，应力——应变关系为一条直线。

2．强化阶段，在这一阶段，应力值还在上升，但是其切线刚度值要小于弹性阶段时的弹性刚度 3．软化阶段，由于在这一阶段，混凝土内部形成了损伤，并且损伤连续发展成微小的裂纹，大大地降低了材料的承载力，应力——应变曲线呈现出一个下降的趋势。

图1.应力——应变曲线

图2.弧长法简介

作为一个数值模拟程序，如何判断其准确性与否，其中主要的一部分是检验应力——应变曲线是否能够拟合一些实验数据。应力——应变曲线本身代表了应用数值模拟形式来表现材料解析的本构关系的这一过程，所以在这一过程中笔者有必要获得一个与以往混凝土压缩实验得到的应力——应变曲线相类似的曲线，这一曲线形式最好能够如图一所示，既有很明显的三个阶段，当然在数值上也必须满足混凝土实验的结果。

为了达到这一目的，笔者参考了相关的书籍和软件的帮助文件，发现了许多的问题，而主要困扰笔者的问题可以概括为这样几个方面：

1．混凝土软化阶段如何用数值模拟。

2．在abaqus/standard中平衡方程的迭代是决定于有效应力还是决定于柯西应力。

3．关于子程序的一些疑问。

一、 混凝土软化阶段的数值模拟

利用newton-raphson迭代在很多问题中可以获得令人满意的解决方案。但是对于一些具有损伤状态的本构关系时，由于存在承载力的下降，利用

Page 1 of 8

第 2 页

newton-raphson迭代在大多数情况下是无法得到下降段曲线。在abaqus中自带了弧长法（riks method），这一方法可以解决一大部分承载力下降的问题。这一方法的独特之处在将每一个增量部的长度也作为一个待求的未知量，同时弧长法还采用了寻求单一平衡路径的方法，利用在已经求得的平衡点上做一条切线，切线的长度则由abaqus内部自带的求解方法进行求解。（此求解方法详情可以参考Abaqus Theory Manual.2.3.2 Modified Riks algorithm）然后如图2所示，在已经获得的A1点处做这一切线的正交线。同时取得A1点在平衡面上对应的平衡点，并作此平衡面曲线，其与正交线的焦点便可以作为一个近似的平衡点。以上便是弧长法的基本简介。

处理非线性计算的另外一种方法拟牛顿法（Quasi-newton），在笔者所开发的子程序中这一方法也具有较强的功能，这是由于在笔者的计算中存在一个用来计算应变增量的雅可比矩阵，而在每一次迭代中，雅可比矩阵保持不变，而采用拟牛顿法可以加快计算运转速度，并解决一些非线性问题。

二、 平衡方程的迭代问题

非线性问题的解是否合理主要是通过验证平衡方程的收敛性与否而定的，而材料非线性问题是无法使得平衡方程达到完全的平衡，在abaqus中为解决这一问题设置了一个默认的残余应力参数，具体详细情况在此不做介绍。

笔者求解数值解的过程中，由于材料的非线性情况严重，损伤与塑形的耦合，材料在进入塑形以后，应力——应变曲线会出现明显的软化现象，这使得利用newton-raphson求解平衡方程的过程中遇到更为困难的问题。基于以上原因 在abaqus/strandard模块下只能够采用弧长法进行计算。同时还会遇到这样的一个问题，材料存在损伤情况下储存在abaqus中的应力是柯西应力（cauchy stress）还是有效应力（effective stress），这两个应力之间存在着一个与损伤有关的折减系数，分别运用这两种应力作为平衡方程的求解因素，得到的收敛情况是截然不同的，所以如果不仔细的分析这两组应力及其应用，就不能够清楚地了解到怎样能够获得平衡方程迭代收敛的情况。

在abaqus软件中介绍到应力存储的是cauchy应力，所以在子程序的最后计算结果中应该给出用于求解平衡方程的应力值为cauchy应力，但是笔者利用以cauchy应力作为输出量的umat子程序来分析处理复杂结构时，结果总是不怎么理想。出现的问题的状况总是极像类似：构件在某一位置出现应力局部化，并发生损伤，而后计算基本陷于停滞状态。笔者在原模型基础上尝试过修改残余应力参数；改换应力加载模式为位移加载模式；将20节点完全积分单元转换为20节点减缩积分单元，甚至是8节点一阶单元；加密网格，网格数量达到了20000个；采用弧长法作为分析部；但是上述的问题仍然出现，总结其原因得到如下几条结论：

1. 模型本身存在问题。

2. 采用cauchy应力作为迭代求解的参量是不太合理的。 3. 子程序上存在着漏洞。

Page 2 of 8

第 3 页

但是这些问题同样也提示了笔者，将曾经修改过很多因素综合起来进行计算，所以笔者得到了数值结果3（见后文），这一结果还是令人满意的。

三、 关于子程序的一些问题

子程序是数值计算的基础，一切问题的最根本原因直接或者间接都是由子程序所导致的，就笔者所采用的子程序中有几个困扰笔者的问题。首先就是笔者对于这一模型的力学概念不是非常的明确，由于笔者所需的知识有限，不能够对于这一模型有一个比较全面的了解，于是这一问题使得笔者不得不重新学习与复习曾经的力学知识，尤其在有限元以及损伤力学这一部分。而这一部分中最为紧要的就是对软化段曲线的认识。

根据资料了解到，混凝土材料试件在单轴压缩实验中，应力会出现峰值，而出现峰值后材料便会进入一个软化段。这一软化段的成因从力学上来讲可以认为是由损伤所引起的（在我们的子程序中是这样认为），这时的应力-应变关系不能够用一个常数的弹性模量来决定，而只能采用增量方法来计算，随着应变的累加、损伤的出现，这些因素使得材料的弹性模量有一个衰减，并最终使弹性模量衰减为0（见图1），此时的损伤值D为1，柯西应力（cauchy stress）为0.了解到这一事实后笔者修改了应力积分格式，从原来的求取无损的有效应力（effective stress）为计算结果的方式，转化到了以求解柯西应力（cauchy stress）为结果的方式。这样做的目的很显然是为了能够在画出含下降段的应力-应变曲线，如果采用求解有效应力（effective stress）方式的话是无法得到一个含下降段的曲线，最多能够得到也只不过是一个形状如理想弹塑性材料的本构关系曲线（以上情况都是基于没有循环加载这样的一个假设的前提下）。然而有一个问题不能够忽视的是如果采用了cauchy应力作为计算的输出量来求解问题时，当将子程序运用在复杂的模型上没有一个有效的方法使得计算可以很好的进行，笔者得到的最为理想的损伤模型在数值结果中将给出，但是将求解cauchy应力作为计算的输出量这一方法应用在一个20 节点完全积分单元时则得到了一个非常好的效果。

四、 数值实验结果

笔者在计算的过程中遇到了很多问题，为了解决这些问题笔者修改了模型使其简化，并在子程序上进行了再编译，然后进行数值计算，笔者得到了一些结果，现将其结果列出，列出这些结果的目的在于说明笔者的计算是在正确的路上进行的。

4.1数值实验1

这一部分的数值计算主要计算的是一个20 节点的完全积分单元，在这里简化了模型，并只使用一个20 节点完全积分单元的目的是为了要去验证子程序所提供的本构关系在数值计算上具有多么理想的数值特性，是否能够体现出混凝土材料的应力-应变曲线所具有的弹性阶段，塑形损伤阶段，同时根据经典的弹性力学理论来检验这一模型是否正确。

利用简单模型得到数值计算的结果在模拟混凝土材料曲线和与理论对比上

Page 3 of 8

第 4 页

都很合理，所以笔者的计算在这一目上基本达到了。

图3.mises等效应力图

图4. 损伤分布图

图3所示的20节点立方单元体的一侧受3个方向的位移约束，而另外一个方向受压缩应力其值为3e7，从mises等效应力图可知，应力的集中区域在中间部位，这与立方体试块的压缩实验有几分相识之处，但是笔者无法做一个精确的定量分析，这是由于立方体只划分了一个单元，无法得到一个精确的结论。在损伤分布图上可以了解到，损伤出现受压的一面上，并向下延伸，但是由于底端固定，所以损伤之为零。

图5.应力-应变对应曲线图

图6. 损伤随弧长变化曲线

为了能够了解到材料的每个点的本构关系，笔者选取了如图4所示的一个点绘制应力-应变对应曲线图，从图中可以看到在弹性区域随着应变值的增加应力值程线性增加，而在弹性极限后应力应变关系着转化为曲线，并在峰值后出现下降段，在分析数据时会发现，应力值的下降是由于塑形的累计以及损伤的出现使得割线模量减小，应力-应变曲线峰值对应的应力值为-2.01005E+007（负号表示压缩），此时对应的损伤值为0.195068。这一应力应变曲线很好地说明了我们的子程序可以完善地模拟混凝土材料的在压缩情况下的力学行为。图6为选取点的损伤值随弧长的变化曲线。

4.2数值实验2

在上一个数值模拟实验中由于本构关系在下降段上的变现不够完全，所以笔者修改了一些参数（边界条件与载荷是相同的），使得材料在软化段时的力学行为体现的更为明显。这些修改主要是在子程序中，笔者在子程序中的雅可比矩阵上增加了损伤一项（根据原始的程序所做的修改），使得材料在损伤因子的迭代过程更为稳定。

Page 4 of 8