基于Matlab的有限元网格自动生成算法 | Q4、Q8、Abaqus单元网格

今日给大家带来的主要内容是二维问题下四边形单元有限元网格如何自动生成？

单元网格的形成实际上属于有限元计算中的前处理部分，即确定单元节点信息，当模型较为复杂时，用户可在Abaqus、Ansys等大型商业有限元软件中进行建模，导出网格信息。

当模型较为简单时，如二维平面板模型，用户可基于一些较为基础的网格生成算法，在自己的程序中通过控制模型长、宽等信息，即可生成有限元网格。

看似应用有限，但是在一些比较复杂的领域内，往往需要先在简单的模型中得到理论验证，如此以来，有利于自编程代码的完整性，即前处理、内核计算、后处理于一体。

本篇推文，木木就带着大家学习一下Q4、Q8单元网格的自动生成以及Abaqus网格节点顺序解读。

代码获取：

Q4单元网格

单元自动网格划分

如下图所示，为4节点四边形单元网格生成示意图，图中NXE和NYE分别是模型横向和纵向单元个数，dhx和dhy分别是单元的横向、纵向长度。

平面板模型被划分为若干个小矩形区域，共有4个节点，分别是 ，一个矩形中是1个四边形单元。

该模型总的单元数目和节点数目分别为 ， 。

网格生成算法代码（纵向排序）

global nnd nel nne  nodof eldof  n 
global geom connec dee nf Nodal_loads
global Length Width NXE NYE X_origin Y_origin dhx dhy 
%
%
nnd = 0;
k = 0;
% 纵向排序
for i = 1:NXE
   for j=1:NYE
        k = k + 1;
        n1 = j + (i-1)*(NYE + 1);
        geom(n1,:) = [(i-1)*dhx - X_origin    (j-1)*dhy - Y_origin ];
        n2 = j + i*(NYE+1);
        geom(n2,:) = [i*dhx - X_origin       (j-1)*dhy - Y_origin  ];
        n3 = n1 + 1;
        geom(n3,:) = [(i-1)*dhx - X_origin       j*dhy - Y_origin  ];
        n4 = n2 + 1;
        geom(n4,:) = [i*dhx- X_origin       j*dhy - Y_origin       ];
        nel = k;
        connec(nel,:) = [n1  n2  n4  n3];
        nnd = n4;
        end
end

代码解读

1. 采用全局变量 global的形式，进行变量的传递；
2. 从两层循环结构上看，最外层是 for i = 1:NXE...end说明网格划分的过程中，x不动，遍历每一个y，节点 纵向排序；
3. 先确定节点号 和 ， n3 = n1 + 1、 n4 = n2 + 1说明 和 在 和 的基础上，编码加1；
4. n1 = j + (i-1)*(NYE + 1)行不动，每次按照列增加，说明 按照纵向排序；
5. n2 = j + i*(NYE+1)比 多了一列的节点，说明 与 同行；
6. 的坐标 geom由 相对位置-坐标轴原点 (X_origin,Y_origin)，该数值由主程序中给出；
7. nel = k指的是每次循环中矩形个数，当两层循环结束时， nel指的是全部四边形单元的个数；
8. connec(nel,:)指的是四边形单元节点编码顺序；
9. connec(nel,:) = [n1  n2  n4 n3]指的是该单元由 节点组成， 并不是固定的，按照一定的顺序即可（顺时针或者逆时针）；
10. nnd = n4当循环结束时， 的数值就是节点的最大值，也就是节点的个数。

网格绘制

figure('Name','Q4单元有限元网格模型');
patch('Faces', connec, 'Vertices', geom,  'Facecolor','none','EdgeColor','b',...   'LineWidth',2,'Marker','pentagram','MarkerSize',8,'MarkerFaceColor','r',...
    'MarkerEdgeColor','r');
axis off
% 节点编号显示
for i=1:nnd
    txt =num2str(i);
    text(geom(i,1)+dhx/10,geom(i,2)+dhy/10,txt);
end
% 单元编号显示(与单元节点编码顺序有关)
for j=1:nel
    txt1 = num2str(j);
    text(geom(connec(j,1),1)+dhx/2,geom(connec(j,1),2)+dhy/2,txt1,'Color','red','FontWeight', 'bold');
end

网格生成算法代码（横向排序）

书中所给的代码一律纵向排序，我稍微修改了一下，使得单元按照 方向横向排序，具体如下：

global nnd nel nne  nodof eldof  n 
global geom connec dee nf Nodal_loads
global Length Width NXE NYE X_origin Y_origin dhx dhy 
%
%
nnd = 0;
k = 0;
% 横向排序
for i = 1:NYE
   for j=1:NXE
        k = k + 1;
        n1 = j + (i-1)*(NXE + 1);
        geom(n1,:) = [(j-1)*dhx - X_origin    (i-1)*dhy - Y_origin ];
        n2 = n1 + 1;
        geom(n2,:) = [j*dhx - X_origin       (i-1)*dhy - Y_origin  ];        
        n3 = j + i*(NXE+1);
        geom(n3,:) = [(j-1)*dhx - X_origin       i*dhy - Y_origin  ];
        n4 = n3 + 1;
        geom(n4,:) = [j*dhx- X_origin       i*dhy - Y_origin       ];
        nel = k;
        connec(nel,:) = [n1  n2  n4  n3];
        nnd = n4;
        end
end

代码解读

1. 从两层循环结构上看，最外层是 for i = 1:NYE...end说明网格划分的过程中，y不动，遍历每一个x，节点 横向排序；
2. 先确定节点号 和 ， n2 = n1 + 1、 n4 = n3 + 1说明 和 在 和 的基础上，编码加1；
3. n1 = j + (i-1)*(NXE + 1)列不动，每次按照行增加，说明 按照横向排序；
4. n3 = j + i*(NXE+1)比 多了一行的节点，说明 与 同列；
5. 的坐标 geom由 相对位置-坐标轴原点 (X_origin,Y_origin)，该数值由主程序中给出；
6. nel = k指的是每次循环中矩形个数，当两层循环结束时， nel指的是全部四边形单元的个数；
7. connec(nel,:)指的是四边形单元节点编码顺序；
8. connec(nel,:) = [n1  n2  n4 n3]指的是该单元由 节点组成， 并不是固定的，按照一定的顺序即可（顺时针或者逆时针）；
9. nnd = n4当循环结束时， 的数值就是节点的最大值，也就是节点的个数。

网格绘制

绘图修饰

众所周知，Matlab的可视化能力强的一批，接下来木木稍微修改一下patch函数里面的参数，即可更改填充面的颜色以及标记的形状：

Abaqus-Q4单元

在这个小节，木木想要依据自己的理解，解释一下:

如下图所示，是在Abaqus中建立一个CPS4(Q4)单元的节点顺序：

模型节点读取

在Abaqus建立平面板模型，划分为 个单元，节点、单元编码如下图所示：

由上图可知，Abaqus在进行CPS4单元节点编码时，时一排一排地排序，所以说，在一个单元中节点编码先是同一方向上的1-2，再是第二行的3-4。

在Matlab中使用Readmesh.m函数，将.inp文件的节点、单元信息导入进Matlab中，命令如下：

[geom,connec] = Readmesh('Q4_abaqus.inp');

网格绘制

figure('Name','Q4单元有限元网格模型');
patch('Faces', connec, 'Vertices', geom,  'Facecolor','c','Marker','o','MarkerSize',4,'MarkerFaceColor','k');
axis off
% 节点编号显示
[nnd,ntemp] = size(geom);
[nel,netemp] = size(connec);
for i=1:nnd
    txt =num2str(i);
    text(geom(i,1)+0.5,geom(i,2)+0.5,txt);
end
% 单元编号显示(与单元节点编码顺序有关)
for j=1:nel
    txt1 = num2str(j);
    text(geom(connec(j,1),1)-2.5,geom(connec(j,1),2)-2,txt1,'Color','red','FontWeight', 'bold');
end

由上图可看到，节点、单元顺序与Abaqus显示一致。

Q8单元网格

单元自动网格划分

如下图所示，为8节点三角形单元网格生成示意图，图中NXE和NYE分别是模型横向和纵向单元个数，dhx和dhy分别是单元的横向、纵向长度。

平面板模型被划分为若干个小矩形区域，共有8个节点，分别是 ，一个矩形就是一个四边形单元。

网格生成算法代码

global nnd nel nne  nodof eldof  n 
global geom connec dee nf Nodal_loads
global Length Width NXE NYE X_origin Y_origin dhx dhy 
%
%
nnd = 0;
k = 0;
for i = 1:NXE
   for j=1:NYE
    k = k + 1;
%
     n1 = (i-1)*(3*NYE+2)+2*j - 1;
     n2 = i*(3*NYE+2)+j - NYE - 1;
     n3 = i*(3*NYE+2)+2*j-1; 
     n4 = n3 + 1;
     n5 = n3 + 2;
     n6 = n2 + 1;
     n7 = n1 + 2;
     n8 = n1 + 1;
 %
    geom(n1,:) = [(i-1)*dhx - X_origin    (j-1)*dhy - Y_origin ];
    geom(n3,:) = [i*dhx - X_origin       (j-1)*dhy - Y_origin  ]; 
    geom(n2,:) = [(geom(n1,1)+geom(n3,1))/2  (geom(n1,2)+geom(n3,2))/2];
    geom(n5,:) = [i*dhx- X_origin       j*dhy - Y_origin ]; 
    geom(n4,:) = [(geom(n3,1)+ geom(n5,1))/2  (geom(n3,2)+ geom(n5,2))/2];
    geom(n7,:) = [(i-1)*dhx - X_origin       j*dhy - Y_origin  ];
    geom(n6,:) = [(geom(n5,1)+ geom(n7,1))/2  (geom(n5,2)+ geom(n7,2))/2];    
    geom(n8,:) = [(geom(n1,1)+ geom(n7,1))/2  (geom(n1,2)+ geom(n7,2))/2];  
 %    
    nel = k;       
    nnd = n5;
    connec(k,:) = [n1 n2 n3 n4 n5 n6 n7 n8];
   end
end

代码解读

1. 先确定节点号 、 、 ， n4 = n3 + 1、 n5 = n3 + 2、 n6 = n2 + 1、 n7 = n1 + 2、 n8 = n1 + 1说明 、 、 、 、 在 、 、 的基础上，编码加1、加2；
2. 注意 n1、 n2、 n3的编码结构，与上面示意图相符；
3. 的坐标 geom由 相对位置-坐标轴原点 (X_origin,Y_origin)，该数值由主程序中给出；
4. nel = k指的是每次循环中单元个数，当两层循环结束时， nel指的是全部四边形单元的个数；
5. connec(k,:)指的是8节点四边形单元节点编码顺序；
6. connec(k,:) = [n1 n2 n3 n4 n5 n6 n7 n8]指的是该单元由 节点组成， 并不是固定的，按照一定的顺序即可（顺时针或者逆时针）；
7. nnd = n5循环结束后，节点总个数等于最后一次循环的 n5编码号。

网格绘制

figure('Name','Q8单元有限元网格模型');
patch('Faces', connec, 'Vertices', geom,  'Facecolor','c','Marker','o','MarkerSize',4,'MarkerFaceColor','k');
axis off
% 节点编号显示
for i=1:nnd
    txt =num2str(i);
    text(geom(i,1)+dhx/10,geom(i,2)+dhy/10,txt);
end
% 单元编号显示(与单元节点编码顺序有关)
for j=1:nel
    txt1 = num2str(j);
    text(geom(connec(j,1),1)+dhx/2,geom(connec(j,1),2)+dhy/2,txt1,'Color','red','FontWeight', 'bold');
end

本文的主要参考内容及Matlab代码均来自Amar Khennane编写的《Introduction to Finite Element Analysis Using MATLAB and Abaqus》，想要进一步了解有限元编程的小伙伴可以入手，强烈推荐！

往期推荐