matlab自带有限元工具分析圆孔应力集中问题

有限元分析结构的力学行为是力学人必备的知识技能，一般分析的步骤为几何建模、划分网格、施加荷载和边界条件、求解，下面我将利用matlab自带的有限元求解器来解决圆孔应力集中问题，本文的matlab版本为2020a及以上。

1. 几何模型的建立

matlab有自带的简单几何模型的建立，这里只需要一个方形的板左下角去掉一个1/4的圆形图形。给定方形板的几何特征是1×1，圆孔的半径为0.1，那么可以首先利用matlab创建pde模型：

model = createpde('structural','static-planestress');

这里指定分析类型为结构分析下面的静力问题-平面应力。

然后用decsg(gm,sf,ns)函数创建几何模型，这个函数有三个输入参数，第一个是描述几何特征的多个列向量gm(意思是geometry matrix)，第二个是描述几何相互之间操作关系（比如差集、并集等）的一个符号表达式sf(意思是string function)，第三个参数是给你的每个描述的几何命名ns(意思是name space)。

首先是描述几何特征的列向量：这个列向量第一个数字是1-4的整数，1表示圆，2表示多边形，3表示矩形，4表示椭圆。对于圆来说，1后面跟的数字是圆心的x坐标、y坐标、半径值；对于多边形来说，2后面跟的数字是多边形的顶点个数n、n个顶点的一次x坐标、n个顶点的依次y坐标；对于矩形来讲，3后面跟的是4、4个点的x坐标、4个点的y坐标；对于椭圆的话，4后面跟着的是椭圆中心的x坐标、y坐标、x半轴、y半轴。

本例子中描述几何的列向量有两个，分别是方形（矩形）列向量和圆形列向量，注意因为描述几何的列向量长度并不一样，所以短的要在后面补0对齐长的，然后按两列的形式放在decsg函数的第一个输入参数的位置。我们gm定义为

square_x = 0; square_y=0; square_sidelength = 1;
square = [3, 4, 0, 1, 1, 0, 0, 0, 1, 1]';
circle_x = 0; circle_y = 0; circle_radius = 0.1;
circle = [1,0,0,circle_radius]';
circle = [circle; zeros(numel(square)-numel(circle),1)]; % 补齐
gm = [square, circle];

其次我们首先确定decsg的第三个输入参数ns，两列代表两个几何分别是正方形和圆，我们用字符串'sq'代表正方形，'cir'代表圆形，那么ns可以表达为：

ns = char('sq','cir')';

最后再确定几何相互之间操作关系sf，这里是正方形减去圆，也即是：

sf = 'sq-cir'; % 用正方形的名字减去圆的名字

最后就可以使用decsg函数构建圆孔应力集中的模型的edge了

g = decsg(gm, sf, ns); % 注意模型g只是edge

利用geometryFromEdges函数可以得到模型的几何

geometryFromEdges(model, g); % 从边界构建几何

我们用pdegplot函数可以画出模型的几何信息，包括edge、face，也即是边的编号、面的编号。

pdegplot(model,'EdgeLabels','on','facelabels','on')

效果展示如下

2. 网格的划分

matlab网格划分很有意思，对于二维只有三角形网格，对于三维只有四面体网格，但是有线性和二次的差别。这里模型是二维的，为了较高的精度，我们用二次单元。网格划分函数generateMesh函数，里面hmax参数是最大单元尺寸，这里给出为0.02。GeometricOrder是单元阶数。

generateMesh(model,'GeometricOrder','quadratic','hmax',0.02);

同时可以用pdeplot显示划分好网格的模型

pdeplot(model)

最后划分的网格如下图

3. 给定模型的物理参数

物理参数其实就是弹性模量泊松比，如果是动力学问题，还需要给定密度之类的。matlab里面一句话解决：

structuralProperties(model,'YoungsModulus', 210e9,'PoissonsRatio',0.2);

4. 荷载和边界条件施加

之前的pdegplot已经显示出了模型的边界编号，这样就方便施加荷载了，我们这个例子，下边界需要给定y方形位移约束，左边界需要给定x方向位移约束，上边界需要一个y方向的应力，方便起见，上边界y方向的应力设置为1，边界约束条件以及荷载施加的代码如下：

structuralBC(model,'Edge',3,'Ydis',0) 
structuralBC(model,'Edge',4,'Xdis',0)
structuralBoundaryLoad(model,'Edge',2,'surfaceTraction',[0,1]);

5. 求解

matlab求解的函数很傻瓜，就一个solve就行！！

structuralresults = solve(model);

6. 后处理

我们可以发现，matlab有自带的后处理函数，很方便，这里就直接贴代码了

figure(3)
clf
pdeplot(model,'XYData',structuralresults.Displacement.x, ...
                          'Deformation',structuralresults.Displacement, ...
                          'DeformationScaleFactor',0, 'colormap','jet')
axis equal
figure(4)
pdeplot(model,'XYData',structuralresults.Displacement.y, ...
                          'Deformation',structuralresults.Displacement, ...
                          'DeformationScaleFactor',0, 'colormap','jet')
axis equal
figure(5)
clf
pdeplot(model,'XYData',structuralresults.Stress.yy, ...
                          'Deformation',structuralresults.Displacement, ...
                          'DeformationScaleFactor',0, 'colormap','jet')
axis equal

7. 画图出的结果

y方向位移云图如下：

x方向位移云图如下：

y方向应力云图sigma_yy如下

可见应力集中系数约为3，可以用以下命令得到具体的值：

>> max(structuralresults.Stress.yy)

ans =

    3.0537

以上就是我的全部分享，对于三维结构、动力学、热传导、模态分析等都可以套这个模式计算呢！！