找回密码
 注册
Simdroid-非首页
查看: 430|回复: 21

[基础理论] 主流商业有限元软件计算可压缩材料的错误及弥补方法

[复制链接]
发表于 2014-9-2 19:33:19 | 显示全部楼层 |阅读模式 来自 广东清远
本帖最后由 chunyu 于 2014-9-14 11:02 编辑

     目前我国计算力学所用的工具基本上被国外的主流商业有限元软件(如ABAQUS,ANSYS,NASTRAN等)所垄断。这些软件都是闭源的“黑盒子”,再加上非线性有限元软件本身的复杂性,软件使用者、包括一些专门从事力学工作的研究者通常对这些工具所采用的算法以及实现方法并不是很关注。久而久之,这些工具逐渐成了标准,成了检验其他代码的“裁判”,甚至成了教科书的内容。这些软件真的就“对”吗?
            下面从基础出发来分析一下主流商业软件在处理可压缩材料时的错误,详细的理论推导可参考Z. P. Bazant, M.Gattu和J. Vorel发表在Proc. R. Soc. A上的文章。为便于理解,叙述上可能会采用数学上不太严谨的定义,比如矩阵和二阶张量不做过多区分,函数和映射不做过多区分等,但不会产生误导。
(1)应力与应变的多种形式
           应力是描述物体受力剧烈程度的一个物理量,应变是描述物体变形剧烈程度的一个物理量。理论上,只要能合理地描述受力的剧烈程度和变形的剧烈程度,可以采用不同的方式来定义应力和应变。以单轴拉伸为例,假设截面面积为A0、初始长度为L0的均匀杆件承受大小为P的拉力,变形后的截面面积为A、长度为L。最直观地,应力就有两种定义方式,一是用力除以初始的截面面积,即P/A0;另外一种是力除以变形后的截面面积,即P/A。这两种定义都可以合理地反应杆件中每一点材料受力的剧烈程度,因此都是被接受的,前者称为工程应力,后者成为真实应力。类似地,应变也有两种定义方式,第一种为变形量除以原始长度,即(L-L0)/L0,该种形式的应变称为工程应变。另外一种和定义真实应力的思路一样,即用变形量除以变形后的长度。由于变形过程中长度不断发生变化,因此需要将变形过程分成很多增量步才能计算这种应变。假设每一步的变形增量为一微小量dL,则对应的应变增量为dL/L。利用积分运算对应变增量进行累加,则从L0变形到L产生的应变为int(dl/L),其中int代表积分(integration)运算,积分上下限分别为L0和L。结果为ln(L/L0),其中ln(x)为自然对数函数。这种形式的应变称为对数应变,也称为真实应变。观察工程应变和真实应变,可以发现二者存在确定的关系,即:Ln(L/L0)=Ln(1+(L-L0)/L0)。
           以上只是从直观出发就可以定义两种应力和应变。对更一般的三维结构,还有更多形式、数学上更完备的应力和应变定义,比如基于变形梯度的应变定义。记变形前结构的坐标(矢量)为大X,变形后的坐标(矢量)为小x,x为X的函数。利用x的梯度(即p(x)/p(X),p()表示偏微分算符)即可以全面地描述结构的变形程度和变形特征,该梯度(矩阵)称为变形梯度,通常记为F。为和直觉吻合(即不变形时应变为0,刚体转动时应变也应该为0),人们又利用F定义了一般变形条件下的几种应变,


T表示转置,I为单位二阶张量。m=2,上述应变为Green-Lagrangian应变;m=1,为Biot应变;m=-2,为Alamansi Lagrangian应变;m=0,则趋近于上述对数应变(也叫Hencky应变)。不管哪种应变,都能合理地反映物体变形的剧烈程度,变形程度大,则应变大,变形程度小,则应变小,不变形,则应变为0。
      类似地,应力也存在多种定义方式。最直观的仍然是真实应力,即以变形后的面积来衡量力的剧烈程度。这种应力也叫Cauchy应力。另外还有第一Piola-Kirchhoff应力,第二Piola-Kirchhoff以及Kirchhoff应力等。具体的数学形式可参看:http://en.wikipedia.org/wiki/Stress_measures

     简言之,可以采用多种形式的应力和应变来描述物体的受力和变形。对某一确定的受力和变形状态,这些不同形式的应力和应变的值通常是不同的,见下图所示的同一种材料的单轴拉伸应力-应变关系:




          天哪,这是同一种材料的应力应变关系吗?是的!
          哪一个对啊?都对!定义形式不同而已。但要注意,在不同形式的应力和应变定义下,材料的切线刚度也不一样了!
          选哪种啊?





本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×
 楼主| 发表于 2014-9-2 19:35:29 | 显示全部楼层 来自 广东清远
Simdroid开发平台
本帖最后由 chunyu 于 2014-9-28 10:03 编辑

(2)功共轭
          上次介绍提到,可以定义多种形式的应力和应变来描述结构受力的剧烈程度和变形的剧烈程度。虽然都是描述同一个物体的同一个状态,不同形式的应力和应变的数值并不相同,这种差异在大变形下尤为明显。那么,选用哪一种形式的应力和应变呢?有什么要求吗?我们先看一下功共轭(work conjugate)的概念。
          物体变形必须遵守能量守恒(或功-能转换),即外力做的功一定和物体的变形能及耗散的能量(比如塑性变形会耗散一部分能量)之和相等。对弹性变形,则外力功一定等于物体的变形能,即以弹性势能形式存储起来的能量。为简单起见,仍以单轴拉伸为例,假设截面面积为A0、初始长度为L0的均匀杆件承受大小为P的拉力,变形后的截面面积为A、长度为L。从开始到伸长为(L-L0),外力所做的功为W=1/2*P*(L-L0),乘以1/2是因为外力是从0增加到P的。如果采用工程应力和工程应变(即以变形前的物体为参考),则变形能的密度为1/2*s*e,其中s=P/A0为工程应力,e=(L-L0)/L为工程应变。对整个(变形前的)物体进行积分可以得到总的变形能:1/2*s*e*V0=1/2*P/A0*(L-L0)/L*(A0*L0)=1/2*P*(L-L0),正好等于外力功W。如果采用全真实应力和全真实应变计算变形能:1/2*s*e*V=1/2*P/A*ln(L/L0)*A*L,结果不等于外力功W! 如果根据应力和应变能正确地计算出能量,则应力和应变是功共轭的。可见采用全应力和全应变,工程应力和工程应变是功共轭的一对,但真实应力和真实应变不是。
         如何能根据真实应力和真实应变计算出变形能呢?可以采用增量的方式,即将整个变形过程分成很多小步骤,每一步应力的增量为dP/A,应变的增量为dL/L,则变形能的增量为1/2*dP/A*dl/L*(A*L)=dP*dL。对变形能的增量进行累加,得到最终的变形能为1/2*P*(L-L0),这和外力功W是一致的。
       功共轭也意味着功率共轭,即根据应变率和应力率计算得到的变形能变化率必须要等于外力的功率。
       从上面简单的例子可以看出,并不是任意一对应力和应变都是功共轭的。为使能量守恒和计算结果正确,计算程序必须要确保所采用的应力(率)和应变(率)对是功共轭的。然后再根据编程的方便选择满足功共轭条件的应力形式和应变形式。常见的几种满足功共轭条件的应力-应变对:
     

Strain
Stress
Symmetry
Volume
Orientation
Engineering Strain (based on deformed geometry); True strain; Almansi strain
Cauchy (True stress)
Symmetric
Deformed
Spatial
Engineering Strain (based on deformed geometry); True strain; Almansi strain
Kirchhoff
Symmetric
Original
Spatial
Deformation gradient
First Piola-Kirchhoff (Nominal Stress)
Non-symmetric
Original
Mixed
Green-Lagrange strain
Second Piola-Kirchhoff (Material Stress)
Symmetric
Original
Material


            如果采用增量的形式进行计算,应力率和应变率也必须要满足功共轭条件。定义、推导和证明应力率和应变率是否满足功共轭涉及复杂的张量运算,在此不再深入介绍。目前最常见的客观应力率(即不受刚体运动影响的应力)有如下几种:
    1  Cauchy应力的
Truesdell
    2. Cauchy应力的
Green-Naghdi
    3. Cauchy应力的
Jaumann
    4. Kirchhoff应力的Jaumann

          具体的数学形式可见:http://en.wikipedia.org/wiki/Objective_stress_rates
        
  采用某种应力(率)形式,则相应地必须要采用满足功共轭条件的应变(率)形式,否则,计算结果就是错误的。
  主流商业软件采用的是哪种形式的应力(率)和应变(率)呢?满足功共轭条件吗?


回复 不支持

使用道具 举报

 楼主| 发表于 2014-9-2 19:36:31 | 显示全部楼层 来自 广东清远
(3)主流商业有限元软件不满足功共轭
      从上两回的介绍中可以得知,应力(率)和应变(率)有多种定义形式,只要满足功共轭关系,都可以用来定义结构的受力和变形。但要注意,材料的刚度参数也要和应力-应变的形式一致,比如利用工程应力应变测得的切向模量和利用真实应力应变测得的切线模量并不相等,如果计算时采用的应力应变形式和材料测试时采用的应力应变形式不一样,则需要进行折算。
      从ABAQUS的理论手册上可以得知,ABAQUS采用的应变率是对数应变率(Henchy应变率)。在进行隐式计算时,应力率是Jaumann形式的Cauchy应力率,在显式动力学计算中,采用的是Green-Naghdi形式的Cauchy应力率,进行屈曲分析时,采用的是Jaumann形式的Kirchhoff应力率。为什么选择这么多形式的应力率,不知道!ANSYS,Ls-Dyna和Nastran采用的是Jaumann形式的Cauchy应力率。
       不幸是,与Henchy应变率功共轭的应力率是Jaumann形式的Kirchhoff应力率,而不是Jaumann形式的Cauchy应力率!事实上,不存在于Jaumann形式的Cauchy应力率功共轭的应变率!这不是新发现,这个问题早在上世纪70年代就已经被指出,见:
      Baˇ zant, Z. P. 1971 A correlation study of incremental deformations and stability of continuous bodies.ASME J. Appl. Mech.38, 919–928
       但是,不知什么原因,这个问题并没有被认真对待。或者大家觉得即使有误差,也不会太大。
       除此之外,材料的刚度矩阵也直接根据Cauchy应力形式的定义来计算。
       比较Jaumann形式的Kirchhoff应力率和Jaumann形式的Cauchy应力率可以发现,二者差别在于Si,j*Vk,k这一项,其中Si,j为Cauchy分量,Vk,k为体积变化率。如果Vk,k为0或很小,两者的差别确实不大。早期的非线性有限元主要用于金属结构的分析,在塑性变形时,金属材料的体积保持不变,因此没有发现什么问题。但是对泡沫材料、生物材料、土、混凝土等可压缩的材料,丢掉体积项肯定会引起错误。
        那错误有多大呢?


回复 不支持

使用道具 举报

 楼主| 发表于 2014-9-2 19:37:44 | 显示全部楼层 来自 广东清远
本帖最后由 chunyu 于 2014-9-14 11:07 编辑

(4)导致的错误
以飞机和船舶工业中常用的夹心板的压入实验为例,



当压入深度为心材厚度的80%时,最大力的错误为20~30%,能量的错误约为20%。比如对三种不同硬度的心材:




以正交各向异性的柱的屈曲为例,临界载荷的错误达100%!



对其他结构、其他材料的演算还要继续。

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有账号?注册

×
回复 不支持

使用道具 举报

发表于 2014-9-2 21:35:42 来自手机 | 显示全部楼层 来自 新加坡
觉得中国在技术上有实力开发自己的有限元软件的
回复 不支持

使用道具 举报

发表于 2014-9-3 10:02:54 | 显示全部楼层 来自 日本
少见研究基础理论的,这个要点个赞。但是,似乎问题不少,  比如

* Cauchy应力与Almansi strain共轭, Kirchhoff应力就不能与Almansi strain共轭, 这里差一个Jacobian。
*  不存在于Jaumann形式的Cauchy应力率功共轭的应变率
   You are definitely wrong! 任何客观应力率都有其共轭量。至于
* Jaumann率导致非现实的结果(如纯剪变形时的应力震动)是公认的事实!研究者们高喊着消灭Jaumann率, 但是
* 为什么选择这么多形式的应力率
   有历史的原因,有计算速度的考虑。更根本的是;没有公认的最佳应力率

评分

1

查看全部评分

回复 不支持

使用道具 举报

 楼主| 发表于 2014-9-3 21:55:24 | 显示全部楼层 来自 广东广州
有讨论是好事。
”Cauchy应力与Almansi strain共轭, Kirchhoff应力就不能与Almansi strain共轭, 这里差一个Jacobian。“您说的是对的,所以在体积那里做了Original和deformed的区分。

至于应力的客观性与功共轭,我的理解是,应力(率)的客观性是指应力不受刚体转动影响的这个性质,而功共轭是指根据应力-应变对得到正确的能量,应力的客观性和功共轭没有必然的联系。
“ 不存在与Jaumann形式的Cauchy应力率功共轭的应变率”这一论断不是我的发现,而是Zdeneˇk P. Bažant等大牛们早在70年代就指出来的问题,可参考:
Baˇ zant, Z. P. 1971 A correlation study of incremental deformations and stability of continuous
bodies.ASME J. Appl. Mech.38, 919–928. (doi:10.1115/1.3408976)
Baˇ zant, Z. P. & Cedolin, L. 1991Stability of structures: elastic, inelastic, fracture and damage
theories. New York, NY: Oxford University Press.

选择Jaumann率估计和您说的类似,有多种原因,历史的原因、方便实现等等。这么多年就这样下来了,对金属材料问题不大,出问题的主要是可压缩材料和正交各向异性材料。


回复 不支持

使用道具 举报

发表于 2014-9-3 23:29:24 | 显示全部楼层 来自 挪威
中国缺的是领袖人物。自己写程序还很远。不是技术问题。
回复 不支持

使用道具 举报

 楼主| 发表于 2014-9-4 07:38:26 | 显示全部楼层 来自 广东广州
个人感觉,既缺领袖人物,又缺深厚的技术实力,更缺有深厚技术实力、具有开放心态、能从全局把握和组织的领袖式人物。就现代的需求而言,对项目总负责人的要求更高,大规模、非线性、多物理、智能化、易扩展、易维护等要求必须需要新的技术才能实现。完成这样的需求需要一个核心人物和一个队伍。
回复 不支持

使用道具 举报

发表于 2014-9-4 21:12:07 | 显示全部楼层 来自 日本
chunyu 发表于 2014-9-3 21:55
有讨论是好事。
”Cauchy应力与Almansi strain共轭, Kirchhoff应力就不能与Almansi strain共轭, 这里差一 ...

1.   你对客观的理解不错,但不准确。这连续体力学中导入客观性条件,Truesdell是鼻祖,可参见他的rational mechanics专著。只是他的这一条件也是有问题的,可参见下面的讨论
http://blog.sciencenet.cn/home.p ... =blog&id=635127

2  在数学中客观性条件倒是一个清晰的概念,即不依赖于坐标系的。翻开一本现代微分几何,它会告诉你covariant derivative, Lie derivative是不依赖于坐标系,是客观的。另外,现在连续体力学中使用的各种客观应力率是Lie derivative的特定形式。

为理解你的回答, 拜读了Bazant教授大作的11.1到11.3节,除了对式(11.3.17)存疑外,没有发现问题,因此他得到的结果公式(11.3.18)-(11.3.20)相信是正确的。他在本书中用自己的方法(即没用到Tuesdell1的定义,也没用到微分几何的概念),得到了Truesdell stress rate, 这个很有意思。但是

3  Bazant认为Jaumann stress rate不是客观应力率,原因是11.3.20中不包含这一形式。这一结论是非常不谨慎的,因为要得出这一结论,你必须证明11.3.20是客观应力率的充分且必要条件!更何况微分几何已经给出了这个问题反例,Jaumann stress rate 是客观应力率!
  另外,注意Bazant的原文是Jaumann stress rate cannot be obtained by substituting the most general second-order finite strain approximation指的是Jaumann stress rate 不能表述为11.3.18-20,而不是你认为的功共轭问题。

4  (11.3.23)的上端Jaumman stress rate may be remedied by applying the Jaumman rate of the Kirchhoff stress. 这个有点不可思议,因为前面他已经说了Jaumann stress rate cannot be obtained by substituting the most general second-order finite strain approximation.这一事实并未因其作用到 Kirchhoff stress而改变。也许我的理解有误,这样做的原因是为了导入体积变化项—凑答案?

俗话说尽信书不如不读书,读一读其他的书对比一下如何?如,偏应用的话,看看Bathe的Finite element procedure, 偏理论的可以看看Mardsen的Mathematical foundation of elasticity. J.C. Simon, Gurtin等都是连续体力学的大家,他们的论著也是很好的。
回复 不支持

使用道具 举报

 楼主| 发表于 2014-9-4 22:38:04 | 显示全部楼层 来自 广东广州
太复杂了:)
我把这个问题转述一下,更多的是抛砖引玉,这个问题很值得深入探讨,因为牵涉到工具到底可不可靠的问题。诚挚邀请hillyuan深入地把这个问题探究一下。
回复 不支持

使用道具 举报

 楼主| 发表于 2014-9-11 11:00:21 | 显示全部楼层 来自 广东清远
(5)改正办法

方法1:等软件开发者修改
方法2:通过编写UMAT或vumat。
回复 不支持

使用道具 举报

发表于 2014-9-28 03:15:52 | 显示全部楼层 来自 美国
“天哪,这是同一种材料的应力应变关系吗?是的!
哪一个对啊?都对!定义形式不同而已。但要注意,在不同形式的应力和应变定义下,材料的切线刚度也不一样了!
选哪种啊?”

这个工程上都有相关规定。一般有效塑性应变超过2%,应力-应变曲线应该选择真实应力-应变曲线。
回复 不支持

使用道具 举报

发表于 2014-9-28 05:30:31 | 显示全部楼层 来自 山东济宁
本帖最后由 jishuzu 于 2014-9-28 05:58 编辑

只谈一下我对楼主第一条的一点点看法,如有错误请各位不吝指正!!

楼主在第一条最后写道:
   “天哪,这是同一种材料的应力应变关系吗?是的!
          哪一个对啊?都对!定义形式不同而已。但要注意,在不同形式的应力和应变定义下,材料的切线刚度也不一样了!
          选哪种啊?”
      
        这段话我觉得会对刚刚接触有限元软件的人造成困惑,所以发表一下我的一点点理解。
       名义(工程)应力应变测试的数据没有考虑拉伸或者压缩时试件截面积的变化,在计算名义(工程)应力应变时,截面积A是不变的,是个常数。但是我们在实验能看到明显的截面积变化(塑性材料),所以我们需要通过换算把名义(工程)应力应变转换为真实应力应变。       可以看看William Hosford 的《Solid Mechanics》(2010版)第36、37页,有解释和计算公式,也可在网上搜到!       在abaqus,ANSYS中选哪种关系在帮助文件里都已经有了详细的描述,就是选用真实应力和真实应变,在涉及塑性变形时要将材料力学性能测试后的名义(工程)应力应变换算为真实应力应变。
      最后提一句,真实应力应变不是真实应力和真实塑性应变,还需要换算!网上可查到换算公式就不写了!软件输入参数时注意了!!!!
      多看理论教材,多看多学帮助文件!
      以上是我的一点看法,还请大家指正!谢谢!


回复 不支持

使用道具 举报

 楼主| 发表于 2014-9-28 09:53:52 | 显示全部楼层 来自 广东广州
本帖最后由 chunyu 于 2014-9-28 10:01 编辑

谢谢tonnyw和jishuzu的评论。
一惊一诧是故意夸张了,毕竟这种枯燥的讨论可能不吸引人。目的仍然是想强调:应力和应变有多种度量,在计算中,一定要注意定义的一致性。包括在定义本构关系时,这一点也要注意。对小变形来讲,这不是问题。但对大变形来讲,特别是对程序开发者而言,一定得搞清。如果选用不同的应力和应变度量,切线刚度矩阵可能是不一样的(但可以换算),比如采用代码中如果使用Green–Naghdi应力率(如abaqus 的vumat),而定义本构关系时采用Cauchy应力,则在计算中需要对切线刚度矩阵(由输入的本构参数求得)进行修正(与当前应力有关!),否则会出现计算错误,且这个错误很难发现。
回复 不支持

使用道具 举报

发表于 2014-9-28 11:11:14 | 显示全部楼层 来自 美国
chunyu 发表于 2014-9-2 19:35
(2)功共轭
          上次介绍提到,可以定义多种形式的应力和应变来描述结构受力的剧烈程度和变形的剧烈 ...
所提到的一维例子有待商榷,非线性问题,外力P和位移并不是线性的关系,应变能的计算公式可能存在问题。

回复 不支持

使用道具 举报

发表于 2014-9-28 11:17:42 | 显示全部楼层 来自 美国
Jaumann formulation 一般用于显式有限元计算而且要求时间步长极小才行,隐式运算采用Hencky formulation。

回复 不支持

使用道具 举报

发表于 2014-9-28 11:23:28 | 显示全部楼层 来自 美国
chunyu 发表于 2014-9-2 19:37
(4)导致的错误
以飞机和船舶工业中常用的夹心板的压入实验为例,

需要对你所提到的两个算例做一下法医学分析,然后再做出结论。由于是非线性大变形,问题可能出在很多地方。最好列出所有操作细节,例如单元类型,应力积分算法,T.L.还是U.L.。
回复 不支持

使用道具 举报

发表于 2014-10-7 16:09:03 | 显示全部楼层 来自 日本
本帖最后由 hillyuan 于 2014-10-7 16:35 编辑
chunyu 发表于 2014-9-28 09:53
谢谢tonnyw和jishuzu的评论。
一惊一诧是故意夸张了,毕竟这种枯燥的讨论可能不吸引人。目的仍然是想强调: ...

我倒是同意你的一惊一诧,这是一个很多人都未必理解的事实。楼上
〉一般有效塑性应变超过2%,应力-应变曲线应该选择真实应力-应变曲线。
〉所以我们需要通过换算把名义(工程)应力应变转换为真实应力应变。
是不正确的或至少是不准确的。选择何种应力-应变曲线取决于你采用何种本构方程,如你采用超弹性本构模型,一般要选用PK应力-Green-Lagrangian应变关系,而不是真实应力-应变曲线。

你的
〉比如采用代码中如果使用Green–Naghdi应力率(如abaqus 的vumat),而定义本构关系时采用Cauchy应力,则在计算中需要对切线刚度矩阵(由
〉输入的本构参数求得)进行修正(与当前应力有关!),否则会出现计算错误,且这个错误很难发现。

这个不会,如果你的应力,等效节点力的计算没问题的话。但是会影响收敛速度。
回复 不支持

使用道具 举报

发表于 2014-10-7 16:20:55 | 显示全部楼层 来自 日本
本帖最后由 hillyuan 于 2014-10-7 16:23 编辑
tonnyw 发表于 2014-9-28 03:15
“天哪,这是同一种材料的应力应变关系吗?是的!
哪一个对啊?都对!定义形式不同而已。但要注意,在不同 ...

从理论上讲,使用真实应力-应变不是必须的。从这一点来看,我想楼主是对的。
另外采用不同的应力率得到的是不同的真实应力,楼主强调这一点也是对的。
回复 不支持

使用道具 举报

您需要登录后才可以回帖 登录 | 注册

本版积分规则

Archiver|小黑屋|联系我们|仿真互动网 ( 京ICP备15048925号-7 )

GMT+8, 2024-3-28 20:27 , Processed in 0.042162 second(s), 12 queries , Gzip On, MemCache On.

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表