有限元编程之杆单元

0.前言

      最近一段时间重新学习了下有限元分析，果然温故而知新，主要是加深了对有限元概念的理解。接下来我跟大家分享下近期用Mathematical编写的关于杆单元的有限元编程，主要包括形函数定义，单元刚度矩阵的求解与组装，最后求解节点位移等过程。最后分享下学习所得。

​01.代码讲解

      可能部分读者对Mathematical不太熟悉，过段时间我出个这个软件的简单教程，入门还是很快的。不过这并不妨碍我今天想要分享的主题，我分享学习经验的核心是理解有限元概念，至于代码等内容，仅仅是供有需要的人使用，而且并不是每个人都有有限元编程的需要。下面感兴趣的朋友结合注释看看吧。

      首先，研究的问题主要是对应下图的工况，我们通过有限元编程将下面的模型用杆单元来代替，通过选取不同网格尺寸来模拟该模型。

Clear["Global`*"]
(*Parameters*)EA = 1;  定义一些几何尺寸和材料参数L = 1;p = 1*x1;nElements = 2;​​(*Preliminary calculations*)nNodes = nElements + 1;               (*节点个数*)​​eL = L/nElements;                        (*单元长度*)l = Table[i*eL, {i, 0, nElements}] ;     (*节点坐标列表*)Kg = Table[0, {i, 1, nNodes}, {j, 1, nNodes}];  (*初始化刚度矩阵*)Kg // MatrixForm;                         (*初将Kg改为矩阵形式*)​​fn = Table[0, {i, 1, nNodes}];        (*初始化节点荷载列阵*)​​(*通过循环求解每个单元的刚度矩阵，节点荷载，并组装为整体刚度矩阵k，全局节点荷载数组，最后求解kx=F*)For[j = 1, j <= nElements, j++, (*displacement approximation*) Disp = a0 + a1*x1;​ (*Swtch to Nodal Displacements*) Eq1 = Disp /. {x1 -> l[[j]]};将坐标点代入位移插值函数​ Eq2 = Disp /. {x1 -> l[[j + 1]]}; sol = Solve[{Eq1 == u1, Eq2 == u2}, {a0, a1}]; (*这里u1和u2不用随循环改变，因为改成u3，u4没区别，因为我们要求知道形函数和单元区域范围即可求刚度矩阵，形函数与单元节点坐标和插值\节点坐标有关*) {a0, a1} = {a0, a1} /. sol[[1]]   ;         ​​ (*Shape Functions*) Disp; N1 = Coefficient[Disp, u1];      (*求Disp表达式中的u1变量的系数*) N2 = Coefficient[Disp, u2]; Ni = {N1, N2}; (*Derivative of the Shape Function*) dNi = D[Ni, x1]; (*Local stiffness Matrix*) Ke1 = Integrate[EA*dNi*dNi[[1]], {x1, l[[j]], l[[j + 1]]}]; Ke2 = Integrate[EA*dNi*dNi[[2]], {x1, l[[j]], l[[j + 1]]}];  Ke = Table[   Integrate[EA*dNi*dNi[[i]], {x1, l[[j]], l[[j + 1]]}], {i, 1,     2}];(*这里借助Table的隐式循环直接求出单元刚度矩阵，用Ke取代Ke1,Ke2*) Ke // MatrixForm; (*Nodal Forces Vector*) fe = Table[Integrate[p*Ni[[i]], {x1, l[[j]], l[[j + 1]]}], {i, 1, 2}]; (*这里误将x1写成x1，很久都没发现原因，以后命名小写为主*)​ (*Add Element Stiffness to Global Stiffness*) (*这里全局阵的每个元素等于上次迭代的值，加上当前单元的分量值*) Kg[[j, j]] = Kg[[j, j]] + Ke[[1, 1]];   Kg[[j, j + 1]] = Kg[[j, j + 1]] + Ke[[1, 2]]; Kg[[j + 1, j]] = Kg[[j + 1, j]] + Ke[[2, 1]]; Kg[[j + 1, j + 1]] = Kg[[j + 1, j + 1]] + Ke[[2, 2]]; Kg // MatrixForm;​ (*Add Nodal Forces to the Global Forces*) fn[[j]] = fn[[j]] + fe[[1]];   fn[[j + 1]] = fn[[j + 1]] + fe[[2]]; Clear[a0, a1]; ]fn // MatrixForm;Kg // MatrixForm;      写成矩阵形式(*Boundary Condition*)KgBC = Drop[  Kg, {1}, {1}];               (*节点1处位移为0，因此划去全局刚度矩阵的1行1列；节点荷载向量的1行*)\​​fnBC = Drop[fn, {1}];                       (*对应的划去节点荷载的对应行*)(*Solve our System*)uBC = LinearSolve[KgBC, fnBC];         (*解KX=F*)u = Insert[uBC, 0, 1];                (*将节点1的位移值0插入uBC中的位置1处*)​(*Determine the Approximation Functions of each elements*)​uApprox = Table[0, {i, 1, nElements}];​(*根据以上求出的节点位移，代入插值函数，求出每个单元的位移场函数的参数a0，a1*)For[i = 1, i <= nElements, i++, Disp = a0 + a1*x1; Eq1 = Disp /. {x1 -> l[[i]]}; Eq2 = Disp /. {x1 -> l[[i + 1]]}; sol = Solve[{Eq1 == u[[i]], Eq2 == u[[i + 1]]}, {a0, a1}]; {a0, a1} = {a0, a1} /. sol[[1]]; uApprox[[i]] = Disp; Clear[a0, a1] ]

 02.结果分析

     由于在这个案例中用到的是线性单元（例如，杆单元的每个单元只有两个节点），有限元理论告诉我们线性杆单元的应变，应力是不连续的，这似乎与我们的是冲突的，那么线性单元模拟得到的结果精度如何呢？下面我们来看看。

2个单元的情况

▲ 位移插值函数

▲ 位移的有限元结果与理论解

▲应力的有限元结果与理论解

5个单元的情况

▲ 位移插值函数

▲ 位移的有限元结果与理论解

▲ 应力的有限元结果与理论解

50个单元的情况

      位移函数太长了，我就不列啦

▲ 位移的有限元结果与理论解

▲ 应力的有限元结果与理论解

       好啦，大家看到结果的对比图，想必也都发现了尽管线性杆单元的应力不连续，但是当网格划分的足够细时，数值结果是收敛于解析解的。可能有人会说这个还用你告诉我吗？这个是个常识啊。可是，问题就在这里，大家做数值模拟的时候，你的网格尺寸满足收敛了吗？如果没有，后续的结果是否有意义呢？好啦，后续我会分享更多学习经验，未来我可能打算做视频讲解这些内容，喜欢的朋友可以关注我的公众号 （获取更多资料和交流欢迎大家关注公众号冬生亦东生），获取完整代码回复杆单元即可，拜拜。