SimPC博士：几何非线性有限元基本原理及matlab编程（上）
导读：大家好，我是仿真秀专栏作者SimPC，近日笔者获悉荣获了“仿真秀2022年度优秀讲师”荣誉称号，在此也感谢广大用户对我的仿真秀专栏内容认可。在未来的每一天，我还将继续在仿真秀官网分享自己的原创视频课程、文章、参与用户问答，为用户提供答疑、个性化培训和技术咨询服务等。
全文目录
我的《Matlab有限元编程从入门到精通10讲》视频课程主要围绕几何非线性有限元分析的matlab编程进行讲解。同样也是以具体案例为对象，先从简单的二维梁单元入手，内容涉及的几何非线性有限元分析的基本理论和针对具体几何非线性有限元分析案例的matlab编程实现。课程主要内容：首先Part1介绍基于不同运动描述方法下的（更新拉格朗日法UL和共旋法CR）的切线刚度矩阵的推导（为什么叫切线刚度矩阵，就是因为几何非线性特性最终反应到有限元编程上就是其刚度矩阵是不断变化的，这也是其非线性的体现）;其次Part2就是非线性求解算法的实现（这就是增量迭代求解的过程，每个迭代步都要对刚度矩阵进行更新，所以我们先讲切线刚度矩阵的推导，再讲非线性求解算法，这里的非线性求解算法不仅仅适合几何非线性也同样适合其他非线性，如材料非线性，我们常听说的非线性求解算法有标准牛顿迭代、修正牛顿迭代、弧长法，这只是非线性求解算法里的一小分支，其实完整的非线性系统的求解算法远比我们所了解的这些名称要复杂)；最后Part3我们会针对具体案例的matlab编程实现来进行讲解，打破有限元理论和实际编程之间的壁垒。
Part1：更新拉格朗日法（UL） ·共旋法（CR） ·切线刚度矩阵 几何非线性所涉及的大变形主要包括大挠度、大转动、大曲率等，先通过几个例子理解大变形条件下为何要几何非线性问题。图1所示悬臂梁端部受足够大弯矩发生大变形。不难想象，梁的实际变形应当是左图的情况，由于没有收到轴向力的作用，梁的中性层的长度不会改变。然而，如果按照结构力学分析中几何线性假定，平衡方程是在未变形状态下建立的，并且不会随变形而更新，那样的话计算得到的梁的变形如右图所示，轴向作用和弯曲作用之间没有耦合，梁的长度会明显增加但没有轴力，其变形和内力并不符合实际，所以线性假定在大变形条件下是不适合的。
图1
如下图2所示，在梁的右端，上面的梁可以自由地水平移动，而下面的梁则不能。在线性理论中，如果梁承受垂直载荷，那么这两个端部条件是等效的。轴向作用和弯曲作用之间没有耦合。但是，在几何非线性分析中，不同的端部条件将导致完全不同的结果：当末端自由轴向移动时，梁的垂直位移几乎与几何线性情况相同。当限制轴向位移时，垂直位移将小于线性情况，并且对载荷具有很强的非线性依赖性。当梁偏转时，如果末端不能向内移动，那么它的中心线就会被拉伸。这将引入很大的轴向力，从而使梁在拉力方面的作用类似于张力中的钢丝——拉力越高，它抵抗横向力的能力就越大。图2
可见结构的刚度会由于几何非线性效应而发生显著变化，上述情况会增加结构刚度；而如果上述梁施加轴向受压的预应力，那么其横向刚度实际上会降低。如图4所示的梁结构，在考虑几何非线性条件下，随着荷载的增加，结构的刚度逐渐减小，出现负刚度，然后结构刚度逐渐增加，其荷载比例与位移的关系如图5所示。可见刚度会由于几何非线性效应而发生显著变化。如果不考虑几何非线性，那么结构刚度保持不变，其荷载比例与位移的关系是线性的，显然与实际情况不符。图5
基于小变形的几何线性假定，来分析一个绕原点在 xy 平面上刚性旋转的二维矩形钢板。下图6显示了矩形钢板旋转了 10°之后的mise应力为 572MPa，高于大多数最常见钢材的屈服极限。显然是不合理的，因为刚性位移并不会产生应力。但使用几何非线性公式求解上述问题，可以避免出现此类虚假应力，如图所示，现在，应力水平是纯数字噪声；比屈服极限低12个数量级，显然后者分析得到的结果更符合实际情况。为什么会出现上述情况呢？因为结构力学和材料力学、弹性力学中我们常见的基于小变形的假定应变叫做工程应变，刚体的旋转会导致工程应变张量的非零分量，这将通过本构定律引起应力。通过理论推导可以解释为何旋转导致工程应变张量出现非零分量：旋转角度phi之后，原先位于（X，Y）的点将移至新的位置（x，y），由下式表示
（1）
（2）
（3）对于刚体旋转，所有应变应为零，但显然其中两个应变分量不是零。金属通常会在 0.001 量级的应变下屈服。这种大小的虚拟应变在 2.5° 旋转时已经发生。为了使应变小于 0.0001，刚性旋转不得大于 0.8°。这意味着，即使是在我们经常希望的“小角度近似”充分的角度下，也必须使用几何非线性方法。使用与上述相同的刚体旋转，但使用格林-拉格朗日应变，得出
（4）
从上述案例可知，大变形情况下应力和应变的度量方式与小变形是不同的，应变由格林-拉格朗日Green-Lagrange应变张量代替工程应变来表示，应力由第二类 Piola-Kirchoff 应力张量表示，二者具体形式将在后文给出。物理现象的数学公式的方程通常表示为微分方程。在大多数情况下，这些是非常复杂的方程和/或应用于复杂的域，这使得无法获得解析解。为了克服这些困难，使用数值方法将连续问题离散为有限自由度，因此微分方程可以写成代数方程。离散化代数方程通常写成矩阵系统，其中系数矩阵将状态变量与控制变量联系起来。当应用基于位移的有限元方法来离散结构力学问题时，状态变量是节点的位移和旋转，而控制变量是外部载荷。在这种情况下，系统矩阵是一个刚度矩阵，它提供了关于向结构施加一定位移所需的力的信息。与材料非线性类似，几何非线性特性最终也是体现在结构整体刚度矩阵随加载过程不断变化的过程，结构整体刚度矩阵因其非线性性质而被称为切线刚度矩阵，切线刚度矩阵由一个与材料相关的部分（传统的弹性刚度矩阵）和一个与几何形状相关的部分（称为几何刚度矩阵）组成，该部分用来表示由于大挠度引起的系统刚度变化，因为其系数涉及从位移获得的内力。为了在拉格朗日运动描述（相对于欧拉描述）范围内建立物体的有限元方程，必须建立用于度量应力、应变以及其他物理量的参考平衡构型。原则上，任何先前获得的平衡构型都可以用作参考构型。然而，在实践中，仅使用初始或最后获得的平衡构型来做参考构型。基于参考构型的选择，常用三种运动学描述来建立非线性平衡方程组：总拉格朗日方法Total Lagrangian（TL）、更新拉格朗日法Updated Lagrangian（UL）和共旋法Corotational（CR）。因为UL和CR相对于TL有更强的适应性，因此本文主要讲解UL方法和CR方法。在几何非线性有限元分析中，引入了一个伪时间变量t，该变量与动力分析无关，只是为了实现荷载的逐步加载，假定加载的速度非常慢，因此可以假设为静态分析，也就是有限元软件中经常提到的准静态分析。一个三维物体在任何伪时刻的构型由其几何形状（表面积S和体积V）、密度ρ和位置信息（由每个粒子的笛卡尔坐标 x1、x2、x3 描述）定义。然后将结构响应定义为随着伪时间从零开始以增量Δt变化而获得的一系列平衡构型。假定从分析开始到时间t的所有平衡构型都是已知的，在时间t+Δt建立平衡构型要根据先前t时刻的构型（或初始构型）建立方程组。具体对于本文要分析的梁单元，如图7所示，在达到t时刻的平衡构型之后，如果期望在新的变形构型（t+Δt）上建立平衡，用于建立平衡方程所需的所有力学和运动学变量可根据t时刻的平衡构型来建立。
图7 更新拉格朗日格式的参考构型
对于几何非线性分析中，任意时刻的切线刚度矩阵怎么得到呢？与线性问题的有限元刚度矩阵推导过程类似，也是通过能量原理（最小势能原理、虚功原理）进行推导，只不过计算内力虚功所用的应力应变是按照上述格林应变和PK-II应力来度量。因为采用的是更新拉格朗日方法，所以上述应力应变度量公式应能根据已知的t时刻参考构型来表达位置的t+dt时刻的构型，因此对于t+dt时刻的格林应变张量和PK-II应力张量可表示为：
（5）
（6）
将上述应力应变表达式带入虚功原理的方程中（内力虚功等于外力虚功）进行推导可以得到：
（7）
表示柯西应力张量；左下标t代表的是[t, t + Δt]区间内的物理量增量；左上标代表t时刻的参考构型对应的物理量。
第一个积分项对应于变形产生内力虚功分量，该分量取决于GL应变张量的线性部分，这里称为内力虚功的线性分量， 
; 第二个积分对应内力虚功的非线性分量 
,该分量取决于GL应变的非线性部分。通过对上述虚功方程进行有限元离散化可得到如下离散化表示的线性增量平衡方程组
（8）
内力虚功的线性和非线性部分分别分别对应上式中的弹性刚度矩阵[KE]和几何刚度矩阵[KG]。外力虚功对应节点负载向量{P}。因为虚位移 
具有任意性，可从方程左右两端约掉，所以上式可进一步退化为
（9）
针对欧拉梁单元对应的内力虚功的线性分量 
和非线性分量 
的表达式为
（10）
（11）
上述非线性分量 
可根据需要决定是否舍去。采用Hermitian插值函数对上述内力虚功表达式进行有限元离散之后，可得
（12）
（13）
对应的弹性刚度矩阵[KE]和几何刚度矩阵[KG]分别为
（14）
（15）对于更新的拉格朗日公式，考虑到应变张量的高阶项和形状函数的复杂度，为梁单元提供了其他不同类型的几何刚度矩阵，这里暂不展开讲解。对于共旋公式(CR)，只建立了一个切线刚度矩阵。该矩阵建立基于梁相对于共旋参考构型的变形很小，因此使用线性应变来度量应变。对于共旋（CR）法，如图8所示，参考构型被分为两部分，分别为初始构型和共旋构型，从而将刚体 位移与产生形变的位移分开。初始构型用于衡量刚体运动，而共旋构型用于衡量应变和应力。共旋构型使用与物体一起移动的局部坐标系，因此，相对于该局部坐标系，共旋构型的刚体运动为零。该方法适用于梁、板和壳等结构单元，特别是对于存在刚体运动，但应变较小的分析。图8
全局坐标系XY下，单元自由度相对于初始构型的总位移，包括平移 𝑑𝑥1、𝑑𝑦1、𝑑𝑥2、𝑑𝑦2 和旋转 𝜃1、𝜃2，如式(16)中向量{dg}。局部系统（也叫自然坐标系）具有三个位移分量，在单元中产生应变和应力，分别是是相对于初始构型的伸长率 un，以及相对于共旋构型的两个相对旋转角 θn1 和 θn2。自然系统的位移分量如式(16)矢量 {dn}。
图9
（16）
全局和自然坐标系下，与上述位移分量对应的节点力分量如下图10所示。
图10
（17）
另外根据图9中所示的几何关系，可推导出全局坐标系和自然坐标系下位移的关系可以通过坐标转换矩阵相互实现转换，同样全局坐标系和自然坐标系下节点力向量也可以通过该坐标转换矩阵实现转换，如下式所示，
（18）
（19）
（20）
对上式方程左右两边求变分可得下式，从而建立起fg与fn的关系（ 
的表达式较复杂），这里暂不列出
（21）
基于上述推导出的全局坐标和自然坐标系下各物理量的关系带入上式，可以推导出下图11黑体公式所示的全局坐标系下节点位移和节点力的表达式图11
（22）
（23）
其中等号右侧中括号内的表达式即为切线刚度矩阵，上述刚度矩阵无需经过旋转可直接组装至整体全局刚度矩阵中。表达式中的第一项取决于线弹性刚度系数矩阵 [Cn]，对应于弹性刚度矩阵。其余项取决于单元内力力，对应几何刚度矩阵。考虑与全局坐标系X轴平行(β = 0)的梁单元，，弹性和几何刚度矩阵的结果分别为
（24）
（25）
由于篇幅原因，本文后续内容还将分中下篇分享，希望大家关注笔者仿真秀专栏，另外我为读者提供以下学习资料，请大家在附件自行下载，如果遇到麻烦请在评论留言，我会立即更新，欢迎大家订阅课程加入我的VIP用户交流。
（完）
声明：本文首发仿真秀App，部分图片和内容转自网络，如有不当请联系我们，欢迎分享，禁止私自转载，转载请联系我们。
获赞 10092粉丝 21553文章 3539课程 219
