损伤在vumat中的实现

ljhzyq · 发表于 2011-1-4 22:59:12

浏览了版中相关的帖子尤其是版主敦诚在下面的帖子中的回复
http://forum.simwe.com/viewthread.php?tid=863730&highlight=%CB%F0%C9%CB
受益很多

我现在想实现的东西很简单在abaqus中已经完全包括了
但是为了以后有更多的修改空间还是决定用vumat实现这个最简单的情况而且方便和abaqus作比较

Abaqus自带: von-Mises+isotropic 硬化失效准则使用damage for ductile metals -> ductile damage 失效过程选择能量控制的线性的方程
Vumat: 在vumat中前面的都简单都有相关的实例可以参考唯独不确定的就是失效后
由于损伤引入的材料软化应力也开始软化
这里我就用了简单的应力更新公式有效应力=应力*(1-D)
应力就是假设损伤不存在下的应力
D按照abaqus帮助中给的和peeq线性相关的公式

在用单胞拉伸测试中用vumat得到的结果和abaqus计算完全吻合
当我应用到一个简单拉伸实验时问题产生了
在第一个单元满足实效准则之后整个试样的力和位移的曲线开始不稳定出现很大的波浪(图1)
提取其中的一个单元看它本身的应力和应变的曲线同样是当实效后波浪发生 (图2)
然而同样的模拟用abaqus本身自带的损伤去做不管力和位移或应力和应变都是非常平滑的曲线

这个问题困扰了很久期望各位高手能够指点一二
有直接解决问题的方法当然好哪怕是一些相关的信息的帮助也非常感谢

敦诚 · 发表于 2011-7-30 23:21:06

损伤和颈缩到底是怎么样对于结构的反映产生影响的，不知道你思考过没有。
我说说我的想法。首先损伤和颈缩出现导致了控制方程的退化，比如从椭圆形方程退化为双曲型方程，而在方程退化过程中就出存在间断解（你可以发现解得震颤是从下降段开始的）。而间断解导致结果发生不稳定效果，这个不是错的，就应该震颤，这正说明了abaqus在解时间差分方程上面具有强大的鲁棒性（这个解的合理性你们可以查阅一下差分格式的数值耗散和数值色斑的概念）。但是我们不希望这样的效果，我们希望的是一天平滑的曲线，所以我们把原本作为错误存在的数值耗散引入到explicit的差分方程中（因为explicit在时间上是差分的），怎么引入数值耗散方法很多，其核心就行让刚度矩阵鳖小于0，具体怎么做，研究的太多了，说不过来，可以去看看相关的文章。
还有我要强调，人家算出来的震颤的解是对的，不要认为人家不对，我们让它变平滑了才是不对的。

敦诚 · 发表于 2011-3-24 00:27:14

单元为1的时候边界都是施加在这个单元上的，这样就相当于增加了单元的约束，提高了刚度，所以衰减行为可以比较缓慢的释放，而多个单元模型中很多单元没有这个优势，所以还是快速衰减，并发生扭曲和失效，建议你定义一个单元删除变量，让单元在damage达到一定数值的时候立刻删除。同时在单元类型里面设置一下比较大的体积粘滞系数，看看还有没有这种情况。
最后减小增量步！做两个分析部一个让材料基本接近损伤的起始，第二个分析部设置步长很小，比如1e-9，然后让他每一步一存储。来看看损伤行为的发展是不是被explicit跳跃了。这样做可能odb文件会很大，但是值得，我以前做过，odb文件甚至达到过20G。

敦诚 · 发表于 2011-1-5 00:48:40

第一点：你的应力更新公式中是否存在迭代求解，如果存在和那么出现震颤可能性非常大。
第二点：是否材料的下降段非常剧烈，一般设置一个粘性因子来缓慢衰减，否则会有震颤，尤其是一个单元损伤很大，刚度很低，这个时候波速变慢，而变形过快，就会导致explicit中的中心差分方法的震颤因素。
第三点：设置体积阻尼等参数来衰减结果的震颤，abaqus中的损伤模型可能都有！

jerson1020 · 发表于 2011-1-10 09:42:32

敦诚版主：
我采用VUMAT做一根悬臂柱的轴压准静态分析，当材料出现刚度退化时，单元由受压变成受拉，应变增量由负值变为正值，请问这个是什么原因呢？我一直调试都不对，请指教！

ljhzyq · 发表于 2011-1-11 00:03:17

第一点：你的应力更新公式中是否存在迭代求解，如果存在和那么出现震颤可能性非常大。
第二点：是否材料的下降段非常剧烈，一般设置一个粘性因子来缓慢衰减，否则会有震颤，尤其是一个单元损伤很大，刚度很低，这个 ...
敦诚发表于 2011-1-5 00:48

首先多谢版主的关注和快速的帮忙可惜前几天没时间上来
关于你给的建议
首先应力更新公式中并无迭代求解稍后我会上传程序和相关的excel文件及图片方便大家帮忙分析
后面两条建议我的水平有限还不能完全理解希望版主能稍微多阐述点或者给出相应的关键词我能在帮助中读一读
不过我想补充一点附加的信息可能也比较重要
首先是这个子程序在于单胞运算所有曲线都很平滑
问题只是在当我运行上千elment的运算时这些震颤出现
下面是两幅von Mises应力分布的云图前一张是震颤出现之前后一张是震颤出现的时候

qihu · 发表于 2011-1-12 22:40:41

本构模型有较为强烈的下降段时，计算容易表现出部化效应（损伤、塑性变形、或者应变在局部单元集中）导致计算不稳定，出现所谓网格依赖特性，就是网格划分少时计算结果好，划分多时表现出局部化效应计算不稳定，打个比方来说，对于一个truss单元，承受受压荷载，如果其材料本构有明显下降段，在explicit分析中如果划分1－3个单元，计算能稳定性好，能计算出比较好的承载力曲线，如果划分10个以上单元，在计算时构件中的一个单元会率先进入下降段，其压缩应变迅速增长，其它单元会出现拉伸（虽然此时构件总体还是处于压缩状态），计算根本得不到下降段曲线，一旦越过最高点，构件的承载力就迅速降为0，理论上要通过所谓的梯度理论才能解决这个问题，但梯度理论目前很不成熟，需要颠覆传统有限元的单元构建方式，使用困难，目前能较好解决方法就是如敦诚版主版主所说的引入粘滞。对于你的情况，建议你将模型中的损伤变量粘滞化使之能考虑应变率效应，看效果是不是好一些。给一篇文献希望对你有帮助。

qihu · 发表于 2011-1-12 22:46:20

粘滞化不仅仅是为了考虑应变加载速率效应对本构行为的影响，它能提高模型的稳定性，在一定程度上消除网格依赖，在考虑应变率的本构模型中如果将加载时间取的足够长，则可认为没有考虑应变率效应，但这样处理却能提高模型的稳定性，部分的消除网格依赖和局部化效应。

ljhzyq · 发表于 2011-1-13 07:35:44

本构模型有较为强烈的下降段时，计算容易表现出部化效应（损伤、塑性变形、或者应变在局部单元集中）导致计算不稳定，出现所谓网格依赖特性，就是网格划分少时计算结果好，划分多时表现出局部化效应计算不稳定，打个 ...
qihu 发表于 2011-1-12 22:40

多谢qihu兄弟的详细耐心的解释让我更深入的理解了数值化方面的不少概念
题外话: 受了版主和qihu的帮忙不得不说在这个冬天我感到了simwe的温暖
受qihu老兄的启发既然问题是出现在应变速率的问题上因为我现在只是von Mises Plasticity + isotropic hardening 本构中不考虑任何应变速率的效应现在的explicit时间是0.01, 我会加长这个时间到1或者更长来比较一下
不过对于小的模型这么做还好大的模型计算就太耗时了无奈目前的基础有限对于粘滞规则化还没有任何概念要加载到vumat损伤变量中更是需要时间如果能有相关进一步的资料麻烦还能告知

jerson1020 · 发表于 2011-1-13 14:34:32

请问一下版主，怎么使用第二种方法呢？是在本构里面设置，还是在分析步中调整某个参数？ 2# 敦诚

ljhzyq · 发表于 2011-1-20 21:10:38

版主和qihu以及jerson，
经过我这几天连续的测试和改进子程序
现在之前的现象倒是不见了但是有了些别的问题比如说单元变形扭曲当损伤变大时
先把我的算法写下来大家帮忙分析一下看看对不对
damage=stateOld(i,15)
e=(1-damage) *props(1)
xnu = props(2)
twomu = e / ( one + xnu )
alamda = xnu * twomu / ( one - two * xnu )
.....
eps=(1-damg)*eps
......
stateNew(k,15)=...

ljhzyq · 发表于 2011-1-20 21:17:26

有趣的是用这个算法对于单个的element结果又是完全相同和abaqus自带的ductile断裂模型
即便损伤达到几乎快1
但是对于多个element的模型时候
单元开始大量的扭曲当损伤稍微变大
实在让人非常困惑
先提前谢谢各位的帮忙

birdva · 发表于 2011-4-4 16:11:13

11# 敦诚

兄台，请问你是怎么施加粘性因子或体积阻尼的？

surfin · 发表于 2011-4-4 21:22:58

学习一下。。。

heraculia · 发表于 2011-7-27 15:30:04

楼主可否留个联系方式，相关内容想跟您探讨一下

jdbear · 发表于 2011-7-27 17:05:40

这个粘滞阻尼的损伤算法，郭版主能否提一个大概的思路，供大家参考，谢谢

chaff · 发表于 2012-2-17 12:37:55

谢谢分享

magicalarm · 发表于 2012-2-17 12:54:02

ljhzyq 发表于 2011-1-20 21:10
版主和qihu以及jerson，
经过我这几天连续的测试和改进子程序
现在之前的现象倒是不见了但是有了些别的问 ...

请问楼主
eps代表什么含义？误差么？

请问你的损伤最大值设为多少？1吗还是0.99

chejialing · 发表于 2012-5-2 14:59:05

讨论的太高深了，看来我还真的加倍努力了，谢谢上面各位

sr楚风 · 发表于 2015-4-16 11:19:27

挖坟……失效阶段，应力如何更新的？

账号		自动登录	找回密码
密码			注册

账号		自动登录	找回密码
密码			立即注册

[子程序] 损伤在vumat中的实现

本帖子中包含更多资源

相关帖子

评分

评分

本帖子中包含更多资源

评分

本帖子中包含更多资源