能量平衡、沙漏及结果评估.doc

1、 【子程序】vumat有沙漏问题么？沙漏问题和VUMAT无关，跟你选择的单元有关系，如果你采用减缩积分单元，则会存在沙漏。有限元的一个核心就是单元模型，其思想是采用单元近似连续体，单元内采用形函数进行插值。采用全积分的话，可以精确地积出刚度矩阵，但是采用全积分会导致有限元过刚，例如体积锁死和剪切锁死等，因此很多力学及提出了各种各样的单元模型来解决这些问题。现在用的较多的低阶单元就是一点积分，一点积分的单元由于积分点过少而存在零能模式（沙漏），即在某些变形模式下会出现零应变，这个可以从形函数的公式中推导出来。 所以，沙漏模式是否存在取决你选用的单元，但是你采用ABAQUS的默认设置基本上就可以解决这个问题。不知道我有没有说清楚2、 【基础理论】【概念】剪切锁死、体积锁死、沙漏、零能模式1. 剪切锁死（shear locking）简单地说就是在理论上没有剪切变形的单元中发生了剪切变形。该剪切变形也常称伴生剪切（ parasitic shear）。发生的条件：1.一阶、全积分单元；2.受纯弯状态；产生的结果：使得弯曲变形偏小，即弯曲刚度太刚。解决方法：1.采用减缩积分；2.细化网格；3.非协调单元；4.假定剪切应变法；2.体积锁死（volumetric locking）简单地说就是应该有单元的体积变化的时候体积却没发生变化。该原因是受到了伪围压应力（Spurious pressure stresses ）。发生的条件：1.全积分单元；2.材性几乎不可压缩；二阶单元：对于弹塑性材料（塑性部分几乎属于不可压缩），二阶全积分四边形和六面体单元在塑性应变和弹性应变在一个数量级时会发生体积锁死。二次减缩积分单元发生大应变时体积锁死也伴随出现。但值得注意的是，一阶全积分单元当采用选择性减缩积分（selectively reduced integration）时可以避免出现体积锁死。产生的结果：使得体积不变，即体积模量太大，刚度太刚。解决方法：1.将大应变区域网格细化；2.mixed formulation法；检查方法：输出积分点的围压应力，分析围压应力是否在相邻积分点存在突变，是否显棋格式分布，是的话就说明出现体积锁死。3.沙漏（hourglassing）简单地说就是单元只有一个积分点，周边的节点可以随意变形。发生的对象：1.一阶、减缩积分单元；产生的结果：单元太柔；解决方法：1.对一阶减缩单元，合理细化网格；荷载避免使用点荷载； 2.在大应变区或大应变梯度区使用一阶单元，而不是使用二阶单元。4.零能模式（zero-energy mode）采用一阶减缩积分时会出现零能模式。即单元只有一个积分点，在受弯时该积分点没有任何的应变能，此时此单元没有任何刚度，就无法抵抗变形。解决方法：1.提供人工的“沙漏刚度”；2.细化网格（一般在高度方向至少要有4个单元）3、 请问在abaqus里面如何查看模型是否出现沙漏？是查看伪能ALLAE么？伪能为零就是没有放生沙漏么？如果单纯想观察单元的变形，是在自动状态下查看，还是设置放大系数？答：是的,看伪应变能,如果约占内能的1%左右是没关系的,如果超10%就基本是错的了,里说的还有可以看,单元的形状,如果单元变成了交替出现的梯形形状的,基本也沙漏的比较厉害了你可以查看一下伪应变能与弹性应变能的比值，一般希望控制在5%以内，你也可以看看伪应变能随时间有没有变小，如果变小表明沙漏现象没有恶化。4、 “沙漏”问题是指什么？在非线形动力分析中，采用高斯单点积分的单元计算可以极大的节省运算时间，但单点积分可能引起零能模式，即沙漏模态，如采用多点积分虽不会产生沙漏模态，但计算时间三次方递增。当单元变形的沙漏模态丢失，它对单元应变能的计算没有影响，在动力响应计算时，出现数值动荡，所以要对沙漏模态进行控制。在程序中有相应的沙漏控制项，取默认值一般就可以对沙漏得到较好的控制。总之，当采用单元单点高斯积分进行动力计算时，必须对沙漏模态进行控制。五、A1:有限元方法一般以节点的位移作为基本变量，单元内各点的位移以及应变均采用形函数对各节点的位移进行插值计算而得，应力根据本构方程由应变计算得到，然后就可以计算单元的内能了。如果采用单点积分（积分点在等参元中心），在某些情况下节点位移不为零（即单元有形变），但插值计算得到的应变却为零（譬如一个正方形单元变形为一个等腰梯形，节点位移相等但符号相反，各形函数相同，所以插值结果为0），这样内能计算出来为零（单元没变形！）。这种情况下，一对单元叠在一起有点像沙漏，所以这种模式称之为沙漏模式或沙漏。现在有很多 控制沙漏的专门程序，如控制基于单元边界的相对转动。但这些方法不能保持完备性。: 我主要讲一下物理的稳定性，在假设应变方法的基础上，建立沙漏稳定性的过程。在这些过程中，稳定性参数基于材料的性能。这类稳定性也称为物理沙漏控制。对于不可压缩材料，即使当稳定性参数是一阶的时候，这些稳定性方法也将没有自锁。在建立物理沙漏控制中，必须做出两个假设：1.在单元内旋转是常数。2.在单元内材料响应是均匀的。A2:沙漏(hourglass)模式是一种非物理的零能变形模式，产生零应变和应力。沙漏模式仅发生在减缩积分(单积分点)体、壳和厚壳单元上。LS-DYNA里面有多种算法用于抑制沙漏模式。缺省的算法(type 1)通常不是最有效的算法，但却是最经济的。一种完全消除沙漏的方法是转换到全积分或者选择减缩积分(S/R)方程的单元。但这种方法是一种下策。例如，第一，类型2体单元比缺省的单点积分体单元计算开消大; 其二，在大变形应用时更不稳定(更容易出现负体积);其三，类型2体单元当单元形状比较差时在一些应用中会趋向于剪切锁死(shear-lock)，因而表现得过于刚硬。三角形壳和四面体单元没有沙漏模式，但缺点是在许多应用中被认为过于刚硬。减小沙漏的一个好的方法是细化网格，但这当然并不总是现实的。加载方式会影响沙漏程度。施加压力载荷优于在单点上加载，因为后者更容易激起沙漏模式。为了评估沙漏能，在*control_energy卡片中设置HGEN2，而且用*database_glstat和*database_matsum卡分别输出系统和每一个部件的沙漏能。这一点是要确认非物理的沙漏能相对于每一个part的峰值内能要小(经验上来说Mischourglass energy。对于流体部件，缺省的沙漏系数通常是不合适的(太高)。因此对于流体，沙漏系数通常要缩小一到两个数量级。对流体用基于粘性的沙漏控制。缺省的沙漏方程(type 1)对流体通常是可以的。对于结构部件一般来说基于刚性的沙漏控制(type 4,5)比粘性沙漏控制更有效。通常，当使用刚性沙漏控制时，习惯于减小沙漏系数到0.030.05的范围，这样最小化非物理的硬化响应同时又有效抑制沙漏模式。对于高速冲击，即使对于固体结构部件，推荐采用基于粘性的沙漏控制(type 1,2,3)。粘性沙漏控制仅仅是抑制沙漏模式的进一步发展，刚性沙漏控制将使单元朝未变形的方向变形。类型8沙漏控制仅用于单元类型16的壳。这种沙漏类型激活了16号壳的翘曲刚度，因此单元的翘曲不会使解退化。如果使用沙漏控制8，16号壳单元可以用于解被称为扭曲梁(Twisted Beam)问题。对于单元类型1的体和减缩积分2D体(shell types 13 & 15)类型6沙漏控制调用了一种假设应变协同转动方程。使用沙漏控制类型6和系数1.0，一个弹性部件在厚度方向仅仅需要划分一层类型1的体单元就可以获得正确的弯曲刚度。在隐式计算里面，对于类型1的体单元应该总是使用类型6的沙漏控制(实际上，在V970里面这是自动设置的)。(More on type 6 HG control from Lee Bindeman)类型6的沙漏控制与类型4,5不在于它用了一个假设应变场和材料属性来估算出假设应力场。这个应力在单元封闭域内进行积分得到沙漏力，因此单元表现的像一个有同样假设应变场的全积分单元。这种假设应变场设计成用来阻止纯弯曲中不真实的剪切变形和近似不可压材料中的体积锁死。类型4和5的沙漏控制基于单元体积，波速和密度像在LS-DYNA理论手册中方程3.21那样来计算沙漏刚度。沙漏类型6主要的改进是应力场在单元域内积分。这使得当使用大的长细比或者歪斜形状的体单元时沙漏控制非常鲁棒。类型4和5的沙漏控制对大长细比和歪斜形状单元反应变不好，它趋向于对某些沙漏模式反应的过于刚硬而对其它模式反应得过弱。沙漏控制类型6另一个理论上的优点是对在厚度方向只有一个单元的梁可以在弹性弯曲问题中得到准确的解。要做到这一点，设置沙漏刚度参数为1.0。同样，对弹性材料方形截面杆的扭曲问题，当沙漏系数设为1.0时可以用很少的单元来解。然而，对于非线性材料，用粗糙的网格得到好的结果是不可能的，因为应力场不是像沙漏类型6假设的那样线性变化的。在梁厚度方向上如果没有更多积分点的话，没有办法捕获应力场的非线性状态。对于选择沙漏控制，下面几个问题要考虑。对于单元有大的长细比或者明显歪斜(不管是初始还是变形过程中)，推荐采用类型6的沙漏控制。类型6的沙漏控制通常对软的材料更好，像泡沫或蜂窝材料在计算中会有非常明显的变形。在材料不是特别软或者单元有合理的形状且网格不是太粗糙时，类型4,5和6沙漏控制似乎都能得到同样的结果。这种情况推荐用类型4的沙漏控制，因为它比其它的更快A3:总能量=内能+动能+滑移界面能能量之间是可以相互转化的，但是对于动力学问题，总能量一般是不太变的。也就是能量守恒原理。 沙漏模式也就零能模式，他在理论上是一种存在的一种变形模式，但是在实际模型中是不可能存上的。零能模式就是指有变形，但是不消耗能量。显然是一种伪变形模式，若不加以控制，计算模型会变得不稳定，并且计算出来的结果也是没有多大意义的。要加抵制这种变形模式就得相应的消耗一定的能量，也就是沙漏能，如果这个比值太多，就说明模型和实际的变形有很大的差别，当然是不正确的。这也是缩减积分所付出的代价。用全积分单元可以解决这个问题，但是效率不高，有可能导致体积锁死，过刚的一些问题。沙漏要控制的，沙漏能一般不大于总能量的10%。如果大于这个值，说明你的计算结果不可靠A4:用全积分单元可以很大程度上减少沙漏。全积分在计算动力问题的时候还是有很多的缺点的，所以一般是选用单点积分的方式，因为此而引起的沙漏问题，dyna提供了多种算法可以减少到内能的5以下，）控制沙漏：首先：您的模型如果是很多PART组成的装配体，那么您需要找出沙漏能最大的PART！ 然后：才能进行一些适当的控制！第一：还得从载荷谈起，避免集中载荷； 第二：从网格谈起，尽量做到网格协调；第三：从沙漏控制谈起，采用不同的控制方法；第四：从单元算法谈起，不行就将一个PART化分为多个PART，除了单元算法不同以外，其他相同，采用全积分！ 有限元方法一般以节点的位移作为基本变量，单元内各点的位移以及应变均采用形函数对各节点的位移进行插值计算而得，应力根据本构方程由应变计算得到，然后就可以计算单元的内能了。如果采用单点积分（积分点在等参元中心），在某些情况下节点位移不为零（即单元有形变），但插值计算得到的应变却为零（譬如一个正方形单元变形为一个等腰梯形，节点位移相等但符号相反，各形函数相同，所以插值结果为0），这样内能计算出来为零（单元没变形！）。显然，麻烦来了。所以必须避免这种情况的出现。判别出现0能模式的方法最简单的是察看单元变形情况，如果单元变成交替出现的梯形形状（两的在一起有点像沙漏，windows里面那个动画光标），就得小心了。另外就是察看hourglass energy，最好不要超过总能量的5%.）对PART进行沙漏控制，使用use *hourglass card 和*PART card。感觉与*CONTROL_HOURGLASS card的功能相同，前者控制具体part，后者控制整个模型）书上说，通过使用好的模型方式可以减少沙漏的产生，如网格的细化、避免施加单点载荷、在易产生沙漏模式的部件中分散一些全积分的“种子”单元，从而减少沙漏。我的理解是在那些地方放一些全积分的单元。其他的都用单点积分单元 。我的理解是增加积分点关于沙漏问题，建议看看abaqus的帮助文档，感觉讲的非常好，由浅入深，把深奥的东西讲的很容易理解。沙漏的产生是一种数值问题，单元自身存在的一种数值问题，举个例子，对于单积分点线性单元，单元受力变形没有产生应变能-也叫0能量模式，在这种情况下，单元没有刚度，所以不能抵抗变形，不合理，所以必须避免这种情况的出现，需要加以控制，既然没有刚度，就要施加虚拟的刚度以限制沙漏模式的扩展-人为加的 沙漏刚度就是这么来的。关于沙漏现象的判别，也就是出现0能模式的方法最简单的是察看单元变形情况，就像刚才所说的单点积分单元，如果如果单元变成交替出现的梯形形状，如果多个这样的单元叠加起来，是不是象我们windows中的沙漏图标呢？)畸变：应该可以用自适应网格划分，要注意的是：*CONTROL_ADAPTIVE 和*PART中的相应的设置）164单元不能自动分裂单元。大变形时确实存在网格畸形的问题，一般可以如下解决：1、当畸形单元少，位置不重要时，删除畸形单元，继续计算；2、适当调整参数，减少畸形的出现，或者使畸形单元破坏（ADD_EROSION）。3、使用ALE。版主所说的第一种解决办法，是通过重启动来实现的么？ 、使用ALE。23)今天突然注意到dyna关键字中的*section_shell_ale中的第八项setype中有3个选项，1-lagrange;2-euler,3-ale，而在*control_ale中的第一项dct有4项，1-lagrange;2-euler,3-ale，4-ale ambient。请高手指教：这两个关键字中的选项是不是要对应，那一个优先级更高？还有，如果都选1，就是lagrange单元，那还定义ale干什么？我一直困惑，希望给出详细的解释，谢谢。根据经验判断，*section_shell_ale的定义应该更高一些，因为在dyna中，普通的关键字总是小于特定part的关键字定义的。在dyna中，默认Lagrange单元、结构动力学分析、显示。因此在设置ALE（*control_ale）、结构热力学分析（*CONTROL_SOLUTION）、隐式（*CONTROL_IMPLICIT_GENERAL）的时候，都会发现系统默认的都是前者的现象，这是DYNA的风格，也是为了避免出现意外。)完全重启动是一个全新的分析，只不过需要考虑前面分析的相关的part的变形和应力情况，这个传递需要由关键字*STRESS_INITIALIZATION_OPTION来实现。在dos命令下输入：LS-DYNA I=restartinput.k R=D3DUMPnn与小型重启动不同的是restartinput.k文件种具有完整的关键字输入，节点、单元、增加的PART等，对于需要初始化的part，在该输入文件种的节点、单元数及排列与拓扑关系都应该与上次求解的输入文件中的一样，但part号可以不同。动力显示算法指南我刚才看了一下dyna的范例手册，也是用part来定义，简单方便。我是用hypermesh来作为前处理的，很方便。具体和大家说一下，在hypermesh中 BCs-entity sets-选择comps，然后输入要定义的set的名称，选择要包括的component，最后再选择create就行了。 hypermesh作为dyna的前处理非常好，推荐使用！Q1:我以前做过一个类似轧钢的问题，使用的单位制为国际单位制，当采用实际的密度（7800kg/m3）时，网格畸变得特别厉害，和你所说的沙漏变形有点象。可是无意间将物质的密度修改大（7.8e7）以后，网格就没有畸变那么厉害了，并且计算速度大幅度提高。不知这里面有什么说道，请明白人给指点一下，谢谢！Q2:我做了一个爆炸分析的流固耦合，但是流体和固体没作用上？流固耦合关键字里面有一个mcoup选项，你选的是1，表示只与密度最大的多物质材料耦合，这样的话就只和炸药耦合了，空气对他的作用就没有了六、能量平衡、沙漏及结果评估1.能量ALLIE Internal energy, ALLIE = ALLSE + ALLPD + ALLCD + ALLAE. 内能增量，简称内能ALLKE Kinetic energy, 动能ALLVD Viscous dissipated energy, 粘性耗散能ALLFD Frictional dissipated energy, 摩擦耗散能ALLCD Energy dissipated by viscoelasticity, 粘弹性耗散能ALLWK Work of the external forces, 外力功ALLPK Work done by contact penalties,ALLCK Work done by constraint penalties,ALLMK Work done by propelling added mass (due to mass scaling),ALLEE electrostatic energyALLSD Total static dissipation energy (due to stabilization）ALLDMD Damage dissipation energyALLCD Creep dissipation energyALLSE Elastic strain energy, 弹性应变能（可恢复）ALLPD Inelastic dissipated energy, plastic dissipation 非弹性耗散能，塑性能。ALLAE Artificial strain energy, 伪应变能（也称沙漏能）ETOTAL Energy balance: 能量平衡（有翻译为“总能量”，我觉得不太确切。）注：1 内能是物体内部分子动能和分子势能之总和，是一个不可计算的量，abq里面的内能，我认为指的是系统内能增量，简称内能。塑性能也是增量。 2 Abq中凡是ALL打头的都是对整体（比如一个set，一个part，甚至整个系统）而言，EL或者E打头的就是对单元而言。如何输出能量见图 3 内能ALLIE=ALLSE + ALLPD + ALLCD + ALLAE。ALLPK Work done by contact penalties,ALLCK Work done by constraint penalties,ALLMK Work done by propelling added mass (due to mass scaling),等翻译留给高手来补充了。2 能量平衡根据能量守恒定律，应该有： ETOTAL=ALLKE+ALLIE+ALLVD+ALLSD+ALLKL+ALLFD+ALLJDALLWK =常量（constant）。由于数值上的原因，这个ETOTAL并非恒为常数。aba的manual指出，ETOTAL变化1%就视为常数了。我的经验，对于静态问题这个数值在0时刻为0，然后开始增大，但是一般不会很大，然后微小变化，这样即为合理。如果拿一个线性静力分析的例子来算一算的话，结果中这个值一直保持为0. 对于能量损失小、无系统外力作用的动态问题，ETOTAL则基本保持为一个常数（比如弹性正碰分析）。少数情况下，存在外力做功，ETOTAL随着外力做功而变化。 ETOTAL变化很大的结果应该慎用。到底变化多大就不可接受呢？这个没有固定的界限，Manual原文：“The sum of these energy components is , which should be constant. In the numerical model is only approximately constant, generally with an error of less than 1%”详见Doc文件：Getting Started with Abaqus: Interactive Edition-9. Nonlinear Explicit Dynamics-9.6 Energy balance Getting Started with Abaqus: Interactive Edition-13. Quasi-Static Analysis with Abaqus/Explicit-13.4 Energy balance也可参见庄茁书之P228：9.6能量平衡。3 准静态与能量平衡 准静态分析（Quasi-Static ）中允许存在一定速度和大位移，只是惯性力在此过程中可以忽略（比如加载速率远小于材料波速，加载时间远大于结构的第一阶固有周期的十倍），即可用准静态来模拟。定性上说，缓慢加载的问题可以用准静态来模拟。abq对于这类问题，最好还是用implicit来算，但是如果模型太复杂不好收敛，准静态就是个不错的选择了。 所得结果是否贴近实际工况呢？ 这是个大问题。Manual中指出，动能必须足够小才能认为是准静态：As a general rule the kinetic energy of the deforming material should not exceed a small fraction (typically 5% to 10%) of its internal energy throughout most of the process。也就是说动能与内能的比值：ALLKE/ALLIE不大于5-10%就可以接受。 定性的讲，ALLKE标定动态特征，准静态采把动态问题用静态方法来算，所以这个量就不应该太大，ALLKE太大就说明这个前提已经不存在了，结果也就自然是不准确的了。 另外，如果用线性减缩积分计算，还有可能出现沙漏（后有详述），还应该看ALLAE/ALLIE，这个又有什么要求呢？看Manual： ALLIE, is a summation of all internal energy quantities. The artificial strain energy is approximately 2% of the total internal energy, indicating that hourglassing is not a problem。Manual说比值小于2%就没有问题，但是没有说大于多少就有问题啊！论坛上一些高手一般也认同ALLAE/ALLIE不大于5-10%就可以接受。但是显然这个也不是绝对的界限，只能是大概的定量化而已，好让大家有个参照。 这方面的讨论论坛很多，比如：用能量来考察quasi-static分析结果的疑问/viewthread.php?tid=895257&highlight=%C4%DC%C1%BF 如果发现ALLKE/ALLIE比值太大，该如何更改模型呢？ 一般就是加大加载时间，以使模型更接近静态；用了质量放大系数的话，这个也可能有影响的，可以减小质量放大系数看看。动态分析，尤其注意输出相应能量，以考核结果是否精确。以下摘录Manual里面一段供大家参考：Energy output is particularly important in checking the accuracy of the solution in an explicit dynamic analysis. In general, the total energy (ETOTAL) should be a constant or close to a constant; the “artificial” energies, such as the artificial strain energy (ALLAE), the damping dissipation (ALLVD), and the mass scaling work (ALLMW) should be negligible compared to “real” energies such as the strain energy (ALLSE) and the kinetic energy (ALLKE).In a quasi-static analysis the value of the kinetic energy (ALLKE) should not exceed a small fraction of the value of the strain energy (ALLIE).It is a good practice to output the constraint penalty work (ALLCW) and the contact penalty work (ALLPW) in analyses involving constraints (such as ties and fasteners) and contact. The value of these energies should be close to zero.来自.au:2080/v6.10/books/usb/default.htm?startat=pt03ch06s03at08.html4 质量放大 为加快运算速度，常常采用质量放大，详见大米斑竹的强帖：关于explicit质量放大的一点小心得/viewthread.php?tid=828207&highlight=%D6%CA%C1%BF%B7%C5%B4%F3。对于这类模型也需要考察结果中的ALLAE/ALLSEA和ALLKE/ALLSE，而且这两个比值应该比较小才适宜。这和1楼所述区别不大。 从ALLKE/ALLSE这个比值来看，宜尽量使用局部质量放大（就是那些稳定时间增量比较小的那些单元）；质量放大系数也不宜太大。在一般非成型的动态分析分析中,mass change percentage不宜超过5%；准静态分析的mass scale factor可以大一些，几十甚至几百都可以尝试下。5沙漏hourglass 1）什么是沙漏？-沙漏(hourglass)模式是发生在减缩积分(线性减缩积分的四边形、六面体单元)体、壳和厚壳单元上的一种非物理的零能变形模式，产生零应变和应力。 2）如何判定沙漏？-判别沙漏最简单的方法是察看单元变形情况，（a）肉眼观察：如果如果单元变成交替出现的梯形形状，一般是沙漏；应力应变明显异常；单元变形异常等（b）定量分析：查看ALLAE/ALLIE的比值以及Etotal来判断。Doc文件-Getting Started with Abaqus: Keywords Edition-10.5.7 Postprocessing有句话：Such a variable is the total internal energy, ALLIE, which is a summation of all internal energy quantities. The artificial strain energy is approximately 2% of the total internal energy, indicating that hourglassing is not a problem。注意manual里面的标准（2%）都比较严，一般工程问题都可以放宽到5%，甚至5%-10%。 用的比较多的就是1）查看单元是否出现交替梯形 2）ALLAE/ALLIE 3）如何避免沙漏？a）三角形壳和四面体单元没有沙漏模式，但缺点是在许多应用中被认为过于刚硬b）采用完全积分c）细化网格d）用光滑分析步的幅值曲线加载，amplitude采用smooth模式。e）减缓加载速度f） 避免单点加载。事实上，点载荷在实际中是不存在的，加载的面积总是或大或小的存在的，只是有时候面积特别小才被抽象为点载荷而已。点载荷最好改为等效的压力载荷。 经验表明：点载荷更易激起沙漏模式，也容易导致单元扭曲、不收敛以及应力奇异。 4）定量评判：ALLAE/ALLIE5%-10%. 5）如何控制沙漏？-除了3）中所述几点以外，可以在前处理中添加沙漏控制措施。也可以采用添加inp的方式：*SOLID SECTION, CONTROLS=name, ELSET=elset*SECTION CONTROLS, NAME=name, HOURGLASS=ENHANCED如果采用ST