壳单元的径向返回算法应力更新问题

songzijieszj

关于应力更新算法，书上一般只给出了适用于体单元的公式，若是分析壳单元，应力张量中sig33理应恒为0，根据图片上公式，应力更新式中若想sig33恒为0，则n33也需为0，单根据n的计算式，n为试偏应力/试等效应力，由于试应力sigma33=0，试偏应力肯定不为0，则n33肯定不为0 ，那么壳单元到底如何更新sig33呢，求各位大神解答！感激不尽

guojunhang

sigma33必须为零，不用更新。
epsilon33=弹性部分加塑性部分

songzijieszj

guojunhang 发表于 2017-5-23 19:17
sigma33必须为零，不用更新。
epsilon33=弹性部分加塑性部分

谢谢大神的回答，关于epsilon33，如图片2所示对径向回归过程的推导，是否可认为，由于delta sigma33为0，那么塑性应变增量33也为0，也就是epsilon33只有弹性部分？

songzijieszj

guojunhang 发表于 2017-5-23 19:17
sigma33必须为零，不用更新。
epsilon33=弹性部分加塑性部分

如果只考虑弹性部分，eps33=-v/(1-v)*(eps11+eps22),算出来的和软件自带的本构相差甚远，如果忽略弹性，只考虑塑性，也就是eps33=-(eps11+eps22)，算出的结果差不多，但是在做单单元单轴拉伸模拟的时候，eps33和eps22不相等，差了一点，不知道问题出在哪里……

songzijieszj

guojunhang 发表于 2017-5-23 19:17
sigma33必须为零，不用更新。
epsilon33=弹性部分加塑性部分

问题基本解决了，按照大神的说法，将弹性部分和塑性部分相加，其中塑性部分用上一步算出的塑性应变增量11和22求出塑性应变增量33，再与弹性部分相加得应变增量33，结果有细微差别，应该是按上一步算法得到的应变增量来计算的显式算法造成的。

bluerain09

songzijieszj 发表于 2017-5-24 00:31
问题基本解决了，按照大神的说法，将弹性部分和塑性部分相加，其中塑性部分用上一步算出的塑性应变增量11 ...

楼主，你好，我最近也在被这个问题困扰着，有些问题想和你探讨一下。
首先，你帖子里说的“应力更新式中若想sig33恒为0，则n33也需为0，单根据n的计算式，n为试偏应力/试等效应力，由于试应力sigma33=0，试偏应力肯定不为0，则n33肯定不为0”，其实sig33恒为0，n33不需要为0的，因为n33是跟偏应力相关，而不是sig33。
其次，当材料进入塑性阶段，采用径向返回法计算的时候，跟试应力相关的试应变是怎么取的，我觉得如果33方向的试应变能取的比较准确，就什么问题都没有了。而如果33方向的试应变在进入塑性阶段后没发准确给出，这样就会导致33方向的应力无法保证满足平面应力条件。我查阅过一些相关书籍，里面介绍是对33方向的试应变也采用迭代算法得到，但我还没搞明白。
最后，楼主能否将你现在的算法再说得详细一些？其中“再与弹性部分相加得应变增量33”中的弹性应变是怎么得到的？谢谢

songzijieszj

bluerain09 发表于 2017-5-24 08:54
楼主，你好，我最近也在被这个问题困扰着，有些问题想和你探讨一下。
首先，你帖子里说的“应力更新式中 ...

好的，关于你所说的sig33恒为0，与n无关，我是这么想的: 可以看到应力的更新式包含n，虽说n的表达式中只有偏应力，但n决定了delta sigma，关于n的推导，是根据r来的，也就是drucker公设里的加载面f对sigmaij的偏导，含义为梯度方向，如果sigma33恒为0，r必为0，即n为0，从而导致33方向的应力及应变更新都为0。
关于进入塑性阶段的试应变增量，我姑且命名为epsij，eps11和eps22都应该是程序传入的额，eps33的塑性部分eps33p，我采用上一步牛顿迭代算出的各delta eps11p(k)相加得到delta eps11p,及delta eps22p(k)相加得到delta eps22p，之后由于J2流动理论的体积塑性应变为0，所以eps33塑性部分可以用上一步的结果近似表示为-（delta eps11p+delta eps22p）。
关于弹性 部分，我采用的是弹性阶段的elastic eps33 =-v/(1-v)*（elastic eps11+elastic eps22），其中弹性应变是根据总的应变-上一步得到的塑性应变

至于你说的那本书讲到的迭代法，可否告知书名？

songzijieszj

bluerain09 发表于 2017-5-24 08:54
楼主，你好，我最近也在被这个问题困扰着，有些问题想和你探讨一下。
首先，你帖子里说的“应力更新式中 ...

关于你所说的迭代法，我有了想法，待我验证后回复于你

bluerain09

songzijieszj 发表于 2017-5-24 09:38
好的，关于你所说的sig33恒为0，与n无关，我是这么想的: 可以看到应力的更新式包含n，虽说n的表达式中只 ...

非常感谢楼主详细的回复，那本书名是《Computational methods for plasticity: theory and applications》
您提到的上一步牛顿迭代，是指当前UMAT程序内部，还是通过状态变量储存的上一个增量步的结果？
在进行牛顿迭代的时候，试应力是通过ABAQUS传入的初始应力值和由应变增量用纯弹性刚度矩阵计算得到的应力增量相加得到，但是这其中就涉及到33方向的应变增量是不准确的。
所以，我觉得楼主进行牛顿迭代的试应力是不准确的，那么迭代出来的结果也会不准确，这样会不会带来误差的累积？

songzijieszj

bluerain09 发表于 2017-5-24 09:55
非常感谢楼主详细的回复，那本书名是《Computational methods for plasticity: theory and applications ...

我指的是上一个增量步的结果，关于适应力的准确性，由于我们计算得到的eps33包含弹性部分及塑性部分，其中只有弹性部分才能满足弹性本构，使得试应力33为0，所以肯定有误差的，但是具体解决办法我还没有想出来，但根据计算结果跟软件自带的子程序的对比，比较理想，且当应变达到0.46的时候都还是理想的，这时一般材料都失效了，所以误差的累积应该可以忽略。但是如果你有更好的解决办法，可以一同探讨，哈哈

bluerain09

songzijieszj 发表于 2017-5-24 10:35
我指的是上一个增量步的结果，关于适应力的准确性，由于我们计算得到的eps33包含弹性部分及塑性部分，其中 ...

还有一个问题，就是截面的横向剪切刚度K11,K22和K12，这些参数应该只要根据弹性模量跟截面厚度来设置就可以了吧？

songzijieszj

bluerain09 发表于 2017-5-24 11:01
还有一个问题，就是截面的横向剪切刚度K11,K22和K12，这些参数应该只要根据弹性模量跟截面厚度来设置就可 ...

具体你可以参考二楼大神的主题帖的附件，里边有说明，我做的是lsdyna的二次开发，在lsdyna里是直接设置这个值的，之后用在计算sigma23&13的弹性本构中