有限元基础编程 | 常应变单元（静力分析全面讲述）
本文摘要：（由ai生成）
本文深入介绍了有限元方法在平面问题中的应用，着重讲解了常应变三角形单元的基础概念、计算方法和技术。包括有限元基本格式、等效节点荷载处理、共节点应力平滑处理技术等。通过数值案例展示了考虑节点等效荷载时的位移和应力场结果，并提供了相关代码获取方式。旨在帮助读者理解有限元计算原理，设计程序，强调编程语言的重要性。
恭喜读者来到了第三章—平面单元，有限元的基本列式保持不变，只是对于不同的单元类型，对应的单元刚度有所不同，当然还有一些场量的不同，在程序的设计上是这样的。
本章将从新手的角度尽可能全面的讲述平面问题的有限元描述，对于高级的有限元操作有待以后的更新，目前的愿景是努力使读者在阅读完这个文档后能够明白有限元计算的原理，快速设计出属于自己的一套程序，编程语言不是最主要的（但也是要学的）。闲话不多说，来看看本章中要分享哪些知识点吧：
常应变单元介绍 单元刚度矩阵计算方法：直接坐标法、高斯积分 节点等效荷载类型：表面力、体力（重力） 应力均匀化技术：绕节点直接平均、绕节点加权平均 应力类型：正应力、剪应力、等效Mises应力、最大/最小主应力
数据前处理、后处理、边界条件、整体刚度组装在前文均有讲到，本章不在过多赘述。
单元介绍
常应变三角形单元具有三个节点，每个节点包含两个平动自由度，至于为什么叫"常应变"，咱们待会儿讨论。
大家有没有发现在几乎所有的国内有限元教材中（国外的我看的少），三节点三角形单元都作为平面问题的开端，都喜欢拿它作为第一个案例来讲，我想除了单元结构较为简单之外，还有一个很重要的因素就是：
Clough在1960年发表的用有限元法求解平面应力问题文章作为有限元法发明（工程应用角度）的一个机极其重要的标志。
在学习之余，逐渐养成对学科的兴趣是保持对此持续热爱的前提，了解一些故事背景可以很容易让人产生兴趣，这也就是本文档在分享知识的同时会穿插一些故事背景的原因。好了，让我们来揭开该单元神秘的面纱吧！
单元描述如下图所示，节点序号按照逆时针排列（顺序无所谓，但是要一致），单元共有6个自由度，记作 ，相应的有6个节点力，记作 ：
有限元格式
接下来会带着大家再次熟悉有限元基本格式，在以后的有限元问题描述上我都会不断的串讲有限元基本格式，加深印象。因为格式真的很固定，特别适合计算机处理，当然这不代表有限元简单，只不过很多内容不需要我们去推倒，直接拿来用就行（应用型教学）。
以下的理论读者不要觉得枯燥，是和编程直接相关的理论基础，本书将与编程无关的理论均剔除，在后面的程序中，以下的理论公式均有所体现。
形函数
其中，
其中， 表示单元面积：
位移场
单元内某一点的位移，可通过形函数与节点位移描述。
单元内位移场可描述为：
矩阵表示法：
简写：
几何方程
单元内某一点的应变，可通过几何矩阵 与节点位移描述。
其中，
再次展开：
将式三角形单元形函数代入式三角形单元几何矩阵，可得：
从上式中可以看出，当节点坐标确定时，几何矩阵B中的元素只与具体的坐标位置有关，不随 、 变化，是个常数，可推至单元内任意一点的应变数值都为一个常数，故三节点三角形单元又被称为常应变单元，在实际网格划分过程中，当出现应变梯度较大的区域应适当细化网格，否则将不能反映应变(应力)的真实变化情况，从而导致较大的误差。
本构方程
针对线弹性各向同性材料，平面应力状态下：
平面应变状态：将上式中 的换为 即可。
单元刚度方程
单元刚度方程表示如下式：
等效节点荷载
注意到上式中的右项 是节点力列阵，不止包括读取的前处理数据文件中已给的节点力，还包括一些非节点荷载转化为节点荷载的力，即：
其中， 表示已知的节点力， 表示体力， 表示表面力。
这就牵扯到节点等效荷载的概念了，之前在梁单元的介绍中我们提过一次这个概念，现在来看一下在平面问题中非节点荷载如何等效到节点上吧！
体力
体力，本文就只考虑重力因素了，在平面问题上，将重力只分配到 方向上，进行平均分配，上式中的 可表示为：
其中， 表示单元厚度， 表示密度， 表示重力加速度， 表示单元面积。
三角形单元重力节点等效过程如下图所示：
表面力
体力相对容易处理，表面力值得商榷，本小节好好谈一下表面力是如何等效的吧！
以均布荷载为例，如下图所示，该单元作用有大小为 的均布载荷，该均布载荷作用在单元的边2-3上，并且与 轴的夹角为 ，表面力矢量可表示为 ，式中的 可表示为:
单个单元均布荷载节点等效如下图所示：
当均布荷载作用于多个单元边时，很多人都以为每个节点的值等于荷载除以节点数，都是相同的，如下图右所示，这样的处理方式是不对的，我们应该按照"叠加原则"，在共节点处的等效节点应力大小等于单个单元情况下的首尾等效节点力之和，如下图左所示。
共节点应力平滑处理
由刚度方程求解得到的节点位移后，如果我们还想获得节点的应力云图，我们需要根据几何方程将节点位移转化为单元中心的应变，然后根据本构方程将单元中心的应变转化为应力，接下来就是将单元中心的应力转化为单元节点的应力，由于三节点三角形单元为常应变（应力单元），所以直接外推到单元节点时的应力值保持一样。
给读者抛出一个问题，如果遇到四边形单元，积分点的应力值如何外插到节点处呢？
现在我们可以获得每个单元的节点应力了，那么新的问题又来了：如下图所示，红色节点为8个三角形单元的共节点，这8个单元在共节点处的应力不尽相同，那红色共节点的具体以什么为准呢？
目前有两个方法：
绕节点直接平均（应力值累加除以共节点单元数） 绕节点面积加权平均，带上"加权"两个字是不是瞬间感觉高大上起来了，读者不要怕，我们一层一层揭开这神秘的面纱！
绕节点直接平均
假设这8个单元计算所得在红色共节点上的应力为 ，对其进行平均处理，最终该节点的应力为 ：
该方法适合共节点连接的单元形状和大小一致或相近，否则将引起较大的误差。
绕节点加权平均
当共节点连接的单元形状和大小相差较大时，就是将单元的面积（三维情况下为体积）因素也加入计算，共节点的应力为 ：
以上的处理只是计算结果后处理的一种局部改善，并不能从根本上解决节点应力（应变）精度差的问题。
应力分量
本小节再给读者科普一下弹性力学中常用到的一些应力的种类，总结了一些用到比较多且广为人知的名字：Von Mises应力、最大主应力、最小主应力。
在三维应力状态下，Von Mises应力可表示为：
平面应力状态：
平面应变状态下：
值得注意的是，平面应变状态下， 方向的应力不为0，在编写程序时一定要注意。
平面状态下的最大和最小主应力可表示为：
数值案例
不考虑节点等效荷载
先考虑如下的单侧受拉模型，如下图所示，模型左侧节点的自由度被固定（UX=UY=0），右侧的节点承受大小为100向右的拉力，材料属性见前处理文件：
*Material
1e+06,0.3,1.0,1,1.0
弹性模量，泊松比，厚度，平面应力/应变标识，密度
场输出变量有：
节点位移：Umag、U1、U3 节点直接应力: Sx、Sy、Sxy 节点直接应力: Von Mises、最大主应力、最小主应力
考虑节点等效荷载：体力+表面力
逐步施加非节点荷载，首先施加重力荷载，如下图左所示，重力加速度取9.8；然后在模型上表面施加向下的压强荷载，大小为1，如下图右所示，为了节省篇幅，只对照位移场和VonMises应力场。
位移、应力场精度对比
不考虑节点等效荷载
考虑重力
考虑重力+压强
相关代码
本篇推文所讲理论涉及到的全部代码均上传至《有限元基础编程百科全书》中，感兴趣的读者可加入知识星球进行获取，可在后台回复星球，或扫描下方二维码，即可加入，加入知识星球后，查看置顶文章即可获得《有限元基础编程百科全书》最新更新版本。
