中心开孔方板的Matlab有限元编程

---

前言

这篇博文是我初次学习有限元的一次总结，虽然只是对一个很简单的问题的求解，但这个过程已经包含了有限元分析的完整步骤。对于更加复杂的有限元的分析问题，无非是在建模、网格划分和求解大规模矩阵方程的算法等方面进行优化。所以这篇博文对初学者希望对有限元分析的过程以及编程有一个初步了解是有帮助的，更高级、更难的问题的求解这里是没有涉及的，当然希望有这方面经验的大佬能不吝赐教，这里先谢过啦~~

---

一、问题提出

考虑如下图的中心开孔方板，上部受到均匀的剪切荷载，下部固支。圆孔圆心与正方形中心重合。该板处于平面应力状态，弹性模量为E=100，泊松比μ=0.3,厚度h=1。下面将使用有限元程序求解此结构的变形图。

二、求解步骤

1.单元剖分

按照有限元的分析流程，要先对悬臂梁进行单元剖分，确定单元与结点编号、以及单元的自由度编号。因为这里是平面应力问题，所以可以采用常应变三角形单元进行网格划分，并且采用的是非结构化的网格。我这里采用的Matlab的pdetool工具进行划分的，详细步骤可以参照下面。当然这里也可以自己编程进行划分，但限于篇幅和介绍的详细程度，这里不会进行讨论。

• 在matlab中输入pdetool并回车，打开该工具。

>> pdetool

• 在打开的工具中，进行模型的建立，如下图所示。
  

• 网格划分，点击Mesh——Initialize Mesh进行初次网格划分，再次点击Mesh——Refine Mesh加密网格，可得到下图。
  

• 导出网格数据。点击Mesh——Export Mesh…导出网格，重命名变量名为node e element，这里node存储结点的坐标信息，element存储的是单元信息，e我们暂时不会用到。

• 存储网格划分数据。在Matlab的命令窗口输入以下命令，按我们所需要的格式存储node和element到Mesh.mat文件中。

>> element = element';
>> element = element(:,1:3);
>> node = node';
>> save('Mesh.mat','node','element')

2.单元分析

在该结构中采用的是常应变三角形单元，在整体坐标系中单元刚度矩阵均用下式进行计算(i,j,m进行轮换就可以得到$b_j、c_j、b_m、c_m$)：
$$\begin{gathered} b_i=y_j-y_m \\ {c}_i=-x_j+x_m \\ [k{]}^e=[B{]}^T[D][B]h\Delta \\ [B]=\frac{1}{2\Delta }\left(\begin{array}{cccccc}{b}_i& 0& b_j& 0& b_m& 0\\ 0& c_i& 0& c_j& 0& c_m\\ c_i& b_i& c_j& b_j& c_m& b_m\end{array}\right) \end{gathered}$$
因为求解的问题是平面应力问题，这里的本构矩阵$[D]$采用如下形式：
$$[D]=\frac{E}{1-{\mu }^2}\left(\begin{array}{ccc}1& \mu & 0\\ \mu & 1& 0\\ 0& 0& \frac{1-\mu }{2}\end{array}\right)$$

3.单元组装

根据单元组装的对号入座原则，即可得到整体刚度矩阵。另外整体载荷向量中只有悬臂梁上端的结点的向下的自由度上非零，其它位置均为0。

4.引入支座约束求解

在自由度编号为悬臂梁的左侧位置上为固定支座约束，即位移均为零，据此提取出非支座自由度得到刚度方程求解即可得到位移。

大家看到这里对于第3、4步或许有很多疑问，其实这里对于初学者确实如此，但只要入门了这里在编写代码是很容易的。如果大家需要这个项目的源码，可以在链接appendix.rar下载。源码是有详细的注释的！！！

三、结果

变形图如下图所示