从ABAQUS中索单元几何刚度的测试到几何非线性.doc

从ABAQUS中索单元几何刚度的测试到几何非线性曲哲2008-5-61.测试1所示，索长为2.0，设测试模型为一根两端固定的索，在中点施加横向集中力，如图 索的轴向刚度 EA=1.0，截面积A=1.0。设索的初始线应变为，受到中点为F的横向力的作用后产生的附加应变为gable,则有:几何方程：gable =1门(J + ?2)本构方程：FCable = gEA+ gableEAK =2EA平衡方程:% + ln其中？为索中点挠度，Fcable为索的轴力，EA为索的轴向刚度。由此式可见，随着初始 应变的增大，索的横向刚度将线性增大，而随着挠度的增大，索的横向刚度将以更快的速度增大。在ABAQUS中用Truss单元模拟索，并采用温度应力施加初始应变。在分析中打开几 何非线性开关，得到图 2、3所示的结果。eO=O.O1O图2挠度较小时ABAQUS分析结果与理论计算结果的比较图3挠度较大时ABAQUS分析结果与理论计算结果的比较在索中间点添加质量为1.0的集中质量，在索中间点加载的轴向变形的刚度4EA/l=2.0,横向变形的初始刚度（即 ？=0时）为Ko =2EAeo由此可以计算出具有不同初始应变的索的两阶圆频率如表1所示，ABAQUS计算的索的圆频率与理论计算的结果基本一致。另外在ABAQUS的计算中，施加索的初始应变之前，索横向振动的频率非常小，仅为1e-8量级，施加初始应变后则与理论计算一致，符合实际情况。表1索的圆频率（括号中数值为ABAQUS计算结果）©©=0.001©=0.005©=0.010横向振动0.045(0.045) 0.100(0.100)0.141(0.141)轴向振动1.414(1.407)1.414(1.411) 1.414(1.414)2对于几何非线性的认识上面这个简单的几何刚度的测试算例促使我对几何非线性再次学习并有了一些新的认 识，但对于我这样数学差劲的工科脑袋来说，要看懂那一大堆几何非线性方面的公式推导还真不容易。所以下面仅仅总结一下目前的认识，以免过几天就忘了。在力学三大方程中都可以有几何非线性的身影，归纳如下：（i）平衡方程中的几何非线性一一平衡参考位形（几何刚度）上面的算例是典型的在平衡方程中考虑几何非线性的例子，算例中索的几何刚度产生的本质是平衡参考位形的变化。如果始终在初始位形（即索是直的）建立平衡方程，则索中不 会产生任何索力，而上述推导正是以变形后的位形（即索有了挠度）建立平衡方程，更直观 的就体现在下面这个式子中：同样著名的例子就是高层建筑中的所谓P-?效应，如果始终以结构的初始构形（即没有发生侧移）为平衡参考位形，则其竖向刚度将保持不变。但若以结构发生侧移后（有时可能 是较在的侧移）后构形为平衡参考位形，则其竖向刚度将不断变小。由此可以简单认为，几何刚度就是由模型变形引起的刚度，可以通过在变形后的构形 上建立平衡方程加以考虑。在离散的有限元计算中，变形后的构形是未知的，待求解的，因此以变形后的构形为参考构形在求解上不太可能。所以在增量求解中将增量步设得多一些，并用上一步求解结束时的构形作为参考构形，以近似考虑参考构形变化对模型刚度的影响，这就是所谓的UpdatedLagrange方法。而始终以初始构形为参考构形的则称为Total Lagrange方法。（2）几何方程中的几何非线性一一应变中的二次项通过直接研究线段长度的变化，可以定义以下两种应变张量:Gree n应变:1 /t2他+0tuj,i + ；Uk,i ；Uk,j )(又称为 Lagrange 应变)Almansi 应变:t tUj,i-:Uk,i；Uk,j)(又称为Euler应变)可见，Lagrange应变是以其中左上标表示当前时刻， 左下标表示参考位形对应的时刻。初始位形为参考的，而 Almansi应变是以变形后的位形为参考的。二者都包括了一个位移导数的二次项。 对于小应变的情况，这一项可以忽略，如果再同 时忽略参考位移之间的差别，则二者是相等的。固体力学的有限元程序，比如ABAQUS都采用Lagrange应变。是否考虑几何非线性，对于几何方程的影响就在于是否考虑Green应变中的二次项。对于上文中的算例，其 Green应变应该是GEable但在上文的推导中采用了对数应变 即（Logarithmic Strain ,又称"汉奇应变”或"真应变”）,生able = ln ( 1+?2 )图4比较了采用不同的应变定义进行上文推导得到的不同结果。总的来说，如果变形前后的长度一定，则对数应变最小，小应变次之，格林应变最大。其中小应变即指不考虑二阶 项的格林应变。0.12理论计算（小应变）出ABAQI计算/f/0.10.080.060.040.020.10.40.50.20.3挠度0.14理论计算（格林应变）- abaquM1 /一10.120.10.080.060.040.020.10.40.50.20.3挠度图4采用不同应变定义得到的理论结果与ABAQUS计算结果的比较另外对于上文算例，不论采用哪种应变，当挠度较大时，理论结果与ABAQUS计算结果之间都有或多或少的差异，但以对数应变相差最小。 ABAQUS在它的Theory Manual中讨论了一维情况下的对数应变与格林应变的定义，并把格林应变推广到三维情况，同时称在处理大转动、大位移、小应变的问题时采用格林应变，但在程序中却声称对于几何非线性问题采用对数应变，令人无所适从。（3）本构方程中的几何非线性一一真应变与真应力量测材料的单轴应力-应变关系时，经常采用的是名义应力和名义应变（即小应变），名义应力是指力除以初始面积。而对于大应变情况，应该采用真应力（即力除以变形后的面积）与真应变（即对数应变）来定义材料的本构关系。上文算例中当索中央的挠度为0.5时，对数应变达0.1116,应属于大应变问题，但上述推导的本构方程中虽然采用了对数应变，但应力仍然是名义应力，即FcableA= E( +%able )A不会随着拉伸而变小，这这是引起图4中差异的原因吗？对于索单元，它的截面积是索单元的形函数决定的。所以即使ABAQUS采用真应力，其数值与名义应力也应该没有任何差别。至此，已经分析了上文算例中各个方程中可能包含的几何非线性因素，并认为在推导过程中都已经考虑到了，但理论结果和 ABAQUS计算结果之间还是有图 4中看到的差异，如 何解释？图5中对于初始应变为 0的情况，比如了三种应变定义得到的理论结果与ABAQUS、MARC两种软件计算得到的结果。可见，两种软件之间本身就有较大的差异。这下没招儿 了。挠度图5软件差异