△ SyracuseU ? 作者简介： 雪城第一98k， 美国雪城大学PhD在读；高级摸鱼运动员；有限元科研工作者；PUBG端游爱好者。 经过上一篇我潜水在群里看大家的反响，感觉似乎我写的东西好像对大家的睡眠质量起到了至关重要的作用，为了广大设计师的睡眠质量，我下决心这一期把Beam element（梁单元）的东西更新完。
（也有大白的原因，最近赶项目在研究楼盖舒适度和TMD，郑博的投稿被大白拖了好一阵子，在此say sorry）
关于上一次我夹带私货的提出杆端弯矩是否会稳定性特征值的问题，我特意的写了一篇小作文请教了王博士，然后王博士也回复给我了一篇辨析文，我只能说给我这个半瓶水数学水平上了一课，但是我这里在最后给出了我的理解和我的分析，王博士给我的回复我也就先不放了，我怕误解了别人的表达，正好结合梁单元放出来一起讨论讨论看看。如果有问题可以直接联系白总。
关于这一篇的某些概念，比如型函数，基函数，试探函数等已经在上一篇讨论过了，如果实在回忆不起来建议先翻翻上篇内容。
△ 点击图片即可阅读
1.基本假设 (Basic assumptions) 和控制方程(Govern equation)
假设我们的梁单元在x方向的位移为u，在y方向的位移为v如图一(a)所示：
根据上一篇的内容，正应变和x方向的变形关系可以表达为:
这里强调一下梁单元，梁单元虽然忽略了在轴力作用下的轴向变形，但是并不代表在弯曲作用下没有正应力产生，集中力，均布荷载以及附加弯矩都可以让梁产生弯曲，然后由于弯曲的作用下产生轴向正应力。
这里我想读过上篇内容的应该会有印象，在伽辽金法中，另变分法中的权函数(weight function)等与方程的近似解，那么弱形式就可以表达为：
当然这只是数学上的一种表达，换一种说法你可以从物理方面去理解，做功(或者虚功方程)说白了就是杆件内部的应变能（由于形变产生的能量）与外力做功（外力与作用点的位移的乘积）要相等，而做功与有限元的核心联系就是把一个整体（杆件）离散化，再通过计算每个小微元体（杆件单元）的内力做功（应变能），与外力做功相等，就可以求出近似解（解的形式一般是位移）。把上述思想体现到(4) 上就可以表达为：
注：(4)还停留在最初正应力的范畴，因此将(2)带入(4)，再通过化简就可以得到上面的表达式。
再从无穷小（微元体）的角度出发，把所有的单元都看作一个弹簧，
再 把每个弹簧的位移做功累加：
通过对比(7)和弱形式的公式，是不是从侧面验证了弱形式的物理正确性。
举个简单的例子：假设一根简支梁如图2所示，被划分为3个单元，在均布荷载作用下的位移图可以表达为v
这里如果实在不理解可以接着往下看，仅仅给一个例子帮助理解，后面我会给一个详细的算例。
2.型函数 (Shape function) 和刚度矩阵（Stiffness Matrix）
上一篇中，已经详细的介绍了什么是型函数，这一篇主要是讲型函数怎么确定，首先说一下梁单元有4个型函数(Shape function)，原因是梁的自由度有4个（2个竖向位移d1,d2，2个转角位移r1,r2）如图3所示：
4个自由度各自分离又相互联系，假设四个型函数分别表达为: 
，杆件内部的任意一点的位移值v就可以表达为：
注：具体型函数的定义，上一篇中讲得
很详细，如果有什么疑惑可以 回头去翻一翻。
对于梁单元来说，实际的位移由两部分组成，第一部分是直接的竖直位移产生的竖直位移 (d1,d2)，第二部分是因为转角(r1,r2)产生的竖直位移。
这里给出一个 型函数的推导，剩余d2, r2与d1, r2非常类似，会直接给出结论。
首先我们假设d1的型函数是一个3次函数，并且带有4个未知量：
边界条件1：x=0时，也就是微元体的最左边时位移值等于杆件的位移值也就是说型函数的值在x=0的时候，y=1：
边界条件2：x=0时，微元体最左边的转角值等于0：
边界条件3：x=1时，即微元体最右边的位移值与杆件左端的位移值无关：
边界条件4：x=1时，即微元体最右边的转角值等于0：
我们假设d2的型函数是一个3次函数，并且带有4个未知量：
边界条件1：x=0时，也就是微元体的最左边
的 转角值等于杆件的转角值也就是说型函数的值在x=0的时候，y=1：
边界条件2：x=0时，微元体最左边的转角值等于1：
边界条件3：x=1时，即微元体最右边的位移值与杆件左端的转角值无关：
边界条件4：x=1时，即微元体最右边的转角值等于0：
(9)中d1, d2, r1, r2都是结点位移值，可以看作一个常数，那么对于[v]的导数就可以近似看作对[B]的导数，x1, x2分别是单元的起始点，那么 [10] 就可以看作是[B]的导数的平方在区间 [x1, x2] 上面的积分乘以材料和截面性质EI的积分。
注：因为大部分的有限元计算程序都要引入高斯积分，这里比较原始的方法我计算过程就写的相当于简略一点。
对于 :
注：上述的型函数中自变量x一般用x与L的比值，单元长度为L，实际的所在坐标是x，对于局部坐标系中的自变量一般用x/L
如果再将剩下的全部积分，那么弱形式化简后的左边也就是(10)就会得到结果(14)：
(14) 是不是看着有那么点眼熟，没错就是结构力学中梁单元的刚度矩阵[K]。
假设一根1m长的悬臂梁悬臂梁，截面弹性模量为200MPa，惯性矩为1x10^7 mm^4，划分为2个单元，末端受一个1kN的力，左端为固定支座，如图5所示：
将L=0.5m，I= 1*10^7 mm^4 E=200MPa带入单元①②的刚度矩阵（单位同一为kN，m）：
已知悬臂梁前两个自由度是0，所以求解可以忽略掉这两个已知量：
注：例子只是随手编的，然后随便算了一算，结果如果有什么错误的话，也麻烦理解一下，主要是看一下给出的计算步骤。
首先，我想先大体说一下如何在线弹性的情况下计算失稳的临界荷载（Critical buckling load）
假设由一
根 悬臂梁（Cantilever Beam）仅受一个集中轴力P，两端端点分别为A和B，弹性模量为E，截面的惯性矩为I，EI均为常数，悬臂端有一个未知位移Δ，杆件长为L，失稳模态未知，用y方向的函数值y(x)来表示。 具体如图所示：
因为结构是对称的，向上失稳和向下失稳都有可能，在这里为了方便后面的计算，先假设失稳模态是向下的，如果我们在杆件的固定端切断，取Free body（隔离体？隐约记得是这么叫）做受力分析，因为考虑了大变形，那么P和位移会产生一个制作反力作用在固定端，那么A点的两个支座反力分别可以计算出来，即: 和 ，如图所示
支座反力外力做功都有了，那么杆件在距离坐标原点x的位置的内力也可以根据已知的情况用y， ， 等表达出来：
因为在最初的模型中不考虑存在水平向的剪力，所以杆件内力中也不含有剪力（v=0）:
通过材料力学（如果没记错的话），梁元模型内部的弯矩与y向位移二次导数以及弹性模量和截面惯性矩相关：
一般书上都会写带上负号，此处不带负号是因为我们定义y轴向下为正方向，与一般的默认坐标系不同。
在此处，我们令两种不同方法得到的弯矩内力相等，再将其移项：
为了方便计算，整个等式同时除以EI，再将P/EI用k^2来表示：
如果用待定系数法来求解微分方程（如果觉得引起不适可以自动跳过只看结论）：
微分方程的解由两部分组成，分别是方程的 通解 和方程的特解 组成，用待定系数法，假设特解也是一个常数，再将特解带入原微分方程，可以得到：
此处有3个未知量： ，需要三个方程来求解。
可以解得
边界条件2：在x=0处，固定端转角θ=0,即
可以解得：
最后将边界条件3：在x=L时，位移= 代 入：
因为 =1
再将 带入：
从这个公式看，确实杆端弯矩似乎对最特征值不产生影响，因为通过欧拉公式，M再增大，还是不影响失稳临界值，但是如果返回到受力平衡部分：
假设你 在 杆端加一个大小为 得杆端弯矩
应该不是这个方程 的 解，等式右边带入值=1，但是左端因为 ，所以左边 的 结果不等于1。
如果从应变能
的 方向来论述这个问题，那就是作用一个杆端弯矩在杆件内，因为杆件要失稳，所以一定会有弯矩在杆件末尾，那么外力做功就可以表达为：
一般来讲，有限元计算的时候用到的一个概念叫Stiffness Stiffened Matrix（刚度折减矩阵？），换通俗易懂的话来说就是上式中 的部分，这一部分写入到弱形式的Govern equation中，就表达为刚度折减矩阵。
但是一般的刚度折减矩阵的假定是仅有轴力作用，如上面算例中自由端位置布置一个外力弯矩M，那么在自由端只要存在转角的情况下，外力做功就不仅仅只有 这一项。
但是神奇的是我检查过几个小软件，在计算失稳临界荷载的时候似乎并不能考虑杆端弯矩的作用，因此我个人的观点还是说对于电算结果要有一定的批判思维。
