还剩4页未读,继续阅读
文本内容:
umat子程序报告这一阶段主要目的是针对于应力——应变曲线问题进行讨论,由于本用户子程序是基于ab叫us有限元软件来模拟混凝土材料的弹塑性损伤行为,而众所周知混凝土材料在单轴压缩的状况下的应力——应变曲线具有以下三个较为明显的阶段(见图1).弹性阶段,应力——应变关系为一条直线.强化阶段,在这一阶段,应力值还在提升,但是其切线刚度值要小于弹性阶段时的弹性刚度.软化阶段,由于在这一阶段,混凝土内部形成了损伤,并且损伤连续进展成微小的裂纹,大大地降低了材料的承载力,应力——应变曲线呈现出一个下降的趋势图
1.应力——应变曲线图
2.弧长法简介作为一个数值模拟程序,如何推断其精确性与否,其中主要的一部分是检验应力——应变曲线是否能够拟合一些试验数据应力——应变曲线本身代表了应用数值模拟形式来表现材料解析的本构关系的这一过程,所以在这一过程中笔者有必要获得一个与以往混凝土压缩试验得到的应力——应变曲线相类似的曲线,这一曲线形式最好能够如图一所示,既有很明显的三个阶段,当然在数值上也必需满足混凝土试验的结果为了达到这一目的,笔者参考了相关的书籍和软件的关心文件,发觉了很多的问题,而主要困扰笔者的问题可以概括为这样几个方面.混凝土软化阶段如何用数值模拟.在ab叫us/standard中平衡方程的迭代是打算于有效应力还是打算于柯西应力.关于子程序的一些疑问
一、混凝土软化阶段的数值模拟采用newton-raphson迭代在很多问题中可以获得令人满足的解决方案但是对于一些具有损伤状态的本构关系时,由于存在承载力的下降,采用newton-raphson迭代在大多数状况下是无法得到下降段曲线在abaqus中自带了弧长法(riksmethod)这一方法可以解决一大部分承载力下降的问题这一方法的独特之处在将每一个增量部的长度也作为一个待求的未知量,同时弧长法还采纳了寻求单一平衡路径的方法,采用在已经求得的平衡点上做一条切线,切线的长度则由abaqus内部自带的求解方法进行求解(此求解方法详情可以参考AbaqusTheoryManual.ModifiedRiksalgorithm)然后如图2所示,在已经获得的Ai点处做这一切线的正交线同时取得Ai点在平衡面上对应的平衡点,并作此平衡面曲线,其与正交线的焦点便可以作为一个近似的平衡点以上便是弧长法的基本简介处理非线性计算的此外一种方法拟牛顿法Quasi-newton在笔者所开发的子程序中这一方法也具有较强的功能,这是由于在笔者的计算中存在一个用来计算应变增量的雅可比矩阵,而在每一次迭代中,雅可比矩阵保持不变,而采纳拟牛顿法可以加快计算运转速度,并解决一些非线性问题
二、平衡方程的迭代问题非线性问题的解是否合理主要是通过验证平衡方程的收敛性与否而定的,而材料非线性问题是无法使得平衡方程达到完全的平衡,在ab叫us中为解决这一问题设置了一个默认的残余应力参数,具体具体状况在此不做介绍笔者求解数值解的过程中,由于材料的非线性状况严峻,损伤与塑形的耦合材料在进入塑形以后,应力——应变曲线会消失明显的软化现象,这使得采用newton-raphson求解平衡方程的过程中遇到更为困难的问题基于以上缘由在abaqus/strandard模块下只能够采纳弧长法进行计算同时还会遇到这样的一个问题,材料存在损伤状况下储存在ab叫us中的应力是柯西应力cauchystress还是有效应力effectivestress这两个应力之间存在着一个与损伤有关的折减系数,分别运用这两种应力作为平衡方程的求解因素,得到的收敛状况是截然不同的,所以假如不认真的分析这两组应力及其应用,就不能够清晰地了解到怎样能够获得平衡方程迭代收敛的状况在ab叫us软件中介绍到应力存储的是cauchy应力,所以在子程序的最终计算结果中应当给出用于求解平衡方程的应力值为cauchy应力,但是笔者采用以cauchy应力作为输出量的umat子程序来分析处理简单结构时,结果总是不怎么抱负消失的问题的状况总是极像类似构件在某一位置消失应力局部化,并发生损伤,而后计算基本陷于停滞状态笔者在原模型基础上尝试过修改残余应力参数;改换应力加载模式为位移加载模式;将20节点完全积分单元转换为20节点减缩积分单元,甚至是8节点一阶单元;加密网格,网格数量达到了20000个;采纳弧长法作为分析部;但是上述的问题仍旧消失,总结其缘由得到如下几条结论.模型本身存在问题.采纳cauchy应力作为迭代求解的参量是不太合理的.子程序上存在着漏洞但是这些问题同样也提示了笔者,将曾经修改过很多因素综合起来进行计算,所以笔者得到了数值结果3见后文,这一结果还是令人满足的子程序是数值计算的基础,一切问题的最根本缘由直接或者间接都是由子程序所导致的,就笔者所采纳的子程序中有几个困扰笔者的问题首先就是笔者对于这一模型的力学概念不是特别的明确,由于笔者所需的学问有限,不能够对于这一模型有一个比较全面的了解,于是这一问题使得笔者不得不重新学习与复习曾经的力学学问,尤其在有限元以及损伤力学这一部分而这一部分中最为紧要的就是对软化段曲线的熟悉依据资料了解到,混凝土材料试件在单轴压缩试验中,应力会消失峰值,而消失峰值后材料便会进入一个软化段这一软化段的成因从力学上来讲可以认为是由损伤所引起的在我们的子程序中是这样认为,这时的应力-应变关系不能够用一个常数的弹性模量来打算而只能采纳增量方法来计算,随着应变的累加、损伤的消失,这些因素使得材料的弹性模量有一个衰减,并最终使弹性模量衰减为0见图1此时的损伤值D为1柯西应力cauchystress为
0.了解到这一事实后笔者修改了应力积分格式,从原来的求取无损的有效应力effectivestress为计算结果的方式,转化到了以求解柯西应力cauchystress为结果的方式这样做的目的很明显是为了能够在画出含下降段的应力-应变曲线,假如采纳求解有效应力effectivestress方式的话是无法得到一个含下降段的曲线,最多能够得到也只不过是一个外形如抱负弹塑性材料的本构关系曲线以上状况都是基于没有循环加载这样的一个假设的前提下然而有一个问题不能够忽视的是假如采纳了cauchy应力作为计算的输出量来求解问题时,当将子程序运用在简单的模型上没有一个有效的方法使得计算可以很好的进行,笔者得到的最为抱负的损伤模型在数值结果中将给出,但是将求解cauchy应力作为计算的输出量这一方法应用在一个20节点完全积分单元时则得到了一个特别好的效果
四、数值试验结果笔者在计算的过程中遇到了很多问题,为了解决这些问题笔者修改了模型使其简化,并在子程序上进行了再编译,然后进行数值计算,笔者得到了一些结果现将其结果列出,列出这些结果的目的在于说明笔者的计算是在正确的路上进行的数值试验1这一部分的数值计算主要计算的是一个20节点的完全积分单元,在这里简化了模型,并只使用一个20节点完全积分单元的目的是为了要去验证子程序所供应的本构关系在数值计算上具有多么抱负的数值特性,是否能够体现出混凝土材料的应力-应变曲线所具有的弹性阶段,塑形损伤阶段,同时依据经典的弹性力学理论来检验这一模型是否正确采用简洁模型得到数值计算的结果在模拟混凝土材料曲线和与理论对比上都很合理,所以笔者的计算在这一目上基本达到了图
3.mises等效应力图图
4.损伤分布图图3所示的20节点立方单元体的一侧受3个方向的位移约束,而此外一个方向受压缩应力其值为3e7从mises等效应力图可知,应力的集中区域在中间部位,这与立方体试块的压缩试验有几分相识之处,但是笔者无法做一个精确的定量分析,这是由于立方体只划分了一个单元,无法得到一个精确的结论在损伤分布图上可以了解到,损伤消失受压的一面上,并向下延长,但是由于底端固定,所以损伤之为零图
5.应力-应变对应曲线图图
6.损伤随弧长变化曲线为了能够了解到材料的每个点的本构关系,笔者选取了如图4所示的一个点绘制应力-应变对应曲线图,从图中可以看到在弹性区域随着应变值的增加应力值程线性增加,而在弹性极限后应力应变关系着转化为曲线,并在峰值后消失下降段,在分析数据时会发觉,应力值的下降是由于塑形的累计以及损伤的消失使得割线模量减小,应力-应变曲线峰值对应的应力值为-
2.01005E+007(负号表示压缩),此时对应的损伤值为
0.195068这一应力应变曲线很好地说明白我们的子程序可以完善地模拟混凝土材料的在压缩状况下的力学行为图6为选取点的损伤值随弧长的变化曲线数值试验2在上一个数值模拟试验中由于本构关系在下降段上的变现不够完全,所以笔者修改了一些参数(边界条件与载荷是相同的),使得材料在软化段时的力学行为体现的更为明显这些修改主要是在子程序中,笔者在子程序中的雅可比矩阵上增加了损伤一项(依据原始的程序所做的修改),使得材料在损伤因子的迭代过程更为稳定图
7.mises等效应力图图
8.损伤分布图数值结果给出了如图7所示的为mises等效应力云图,这时的应力分布比较平缓,从一侧开头向另一个面过度而图8的损伤分布图着与上一个数值模拟试验的结果相类似,只不过损伤值有所增大,扩展的范围也加大了图
9.应力-应变对应曲线图图
10.损伤随弧长变化曲线对于与数值结果1笔者也给出了一点的应力-应变曲线,这一条曲线包括了整个混凝土材料在压缩状况下应力的全部下降段曲线的应力的峰值-
3.34283E+007(负号表示压缩),发生在第41分析部,这时对应的损伤值为
0.2584还要说明的是材料也是在这一位置左右开头消失了塑性损伤而这一点的损伤值曲线也很完善,从开头消失损伤,并快速增加,而到达肯定值,
0.8左右,则越来越趋于平滑同时为了检验数值结果的正确性与合理性,笔者采用经典弹性力学单轴压缩的解析公式计算了应力值如图11所示,蓝线代表了解析解,红线代表了数值解两者在弹性区域时吻合的特别好,但是随着损伤的消失,数值解的值开头增加,这是由于损伤消失导致有效截面积的减小,依据损伤力学公式可知,损伤消失使得应力增加,所以消失了图11所示的状况图
11、理论解-数值解对比图(蓝线代表了解析解,红线代表了数值解)数值试验3这是数值计算模型是将数值试验2的子程序以及参数应用在一个简单的模型上,对于这个计算来看,对于几何模型简单的状况使用虽然在简洁模型上特别胜利的子程序还是无法得到一个比较抱负的结果但是这一模型的数值计算从结果的模拟程度上看还是比较胜利的图
12.mises等效应力图图13损伤分布图模型外形为一有缺口的圆柱,圆柱一面固定没有唯一,而另一端受压力作用单元采纳了8节点一次单元Mises等效云图图12看出,沿着柱子的最细半径处应力有所集中,最大值为
4.38e6而损伤值图13所示的集中区域也消失在这一位置图14-15为笔者所选取的一个点绘制的应力-应变图与损伤演化图,损伤值最大为
0.12左右从图13中可以了解到在损伤位置的点上也存在着软化现象,只不过对比单块模型的数值模拟结果没有那么明显图
14.应力-应变对应曲线图图
15.损伤随弧长变化曲线
五、结论这一时期笔者的主要工作就是调整子程序使得能够获得一个合理的数值解,赐予以上的数值计算,以及其与解析解的对比可以了解到笔者在子程序的调试工作上还是基本正确的,但是仍旧存在的问题是笔者将对于简洁模型成立的子程序应用在简单模型时,无法得到一个特别抱负的应力应变曲线,这还需要连续调试程序和修改参数
六、附录理论来源整看以下内容该理论采纳弹性能等效原理,由SaanouniForster和BenHatira1994提出的,并且由ShenShen和Chen2004进一步完善,塑性损伤模型的基本数学关系如下3%尾=1-府,丫=1-O=碌”1-24/-£』这里的盘和%分别是应力张量和应变张量,上标p代表塑性,而e代表弹性,上〜表示虚JJ拟静材料参数是损伤变量,y是损伤共规力,碎心是无损伤材料的弹性张量定义应力三轴比为力是有效应力主值的和是有效应力偏量的第2不变量,1jJ2=—SijSjjS-=a-——323%虚拟偏应力张量J将应力三轴比引入到损伤塑性加载条件和非弹性势函数中得到如下表示式:f=a1++乩1—”]04即是一个材料常数,用来反映材料压力敏感特性;k是初始剪切强度常数;底是虚拟静材料的应变硬化极限这极限相应于等效塑性应变趋于无穷大;参数b是模型的常数,可以通过拟和试验现象来确定0是在有效应力空间中定义的非弹性势的塑性部分幺为损伤塑性乘子下面给出损伤的演化率V-£=2一1一/=若8SysS是材料参数将式8进行线性化处理,得AD=A⑼对应一个应变增量△为,采纳Newton法迭代求解损伤塑性加载因子AZ时,塑性损伤乘子2可以由下列的全都性条件解析求得7=7=0104T制⑴求得了AA便可进一步求损伤,进而求得数值求解过程中的弹塑性切线的刚度张量:%/=成1—十2最终更新应力和其他数据。