基于Matlab的有限元网格自动生成算法 | CST、LST单元网格

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

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

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

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

本篇推文，木木就带着大家，学习一下CST、LST单元网格的自动生成。

CST单元网格

单元自动网格划分

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

 

平面板模型被划分为若干个小矩形区域，共有4个节点，分别是    ，一个矩形中又进一步划分两个三角形单元，第一个单元的节点为    ，第二个单元的节点为    。

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

网格生成算法代码

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 = 2*k;
        m = nel -1;
        connec(m,:) = [n1  n2  n3];
        connec(nel,:) = [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 = 2*k指的是每次循环中矩形个数的2倍，当两层循环结束时，nel指的是全部三角形单元的个数；
8. m = nel -1指的是每次循环中矩形个数的2倍-1；
9. connec(m,:)指的是矩形区域下面那个三角形单元节点编码顺序，connec(nel,:)指的是矩形区域下面那个三角形单元节点编码顺序；
10. connec(m,:) = [n1 n2 n3]指的是该单元由      节点组成，并不是固定的，按照一定的顺序即可（顺时针或者逆时针）；
11. nnd = n4当循环结束时，      的数值就是节点的最大值，也就是节点的个数。

网格绘制

figure('Name','LST单元有限元网格模型');
patch('Faces', connec, 'Vertices', geom,  'Facecolor','c','Marker','o','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:2:nel
    txt1 = num2str(j);
    txt2 = num2str(j+1);   text(geom(connec(j,1),1)+dhx/3,geom(connec(j,1),2)+dhy/4,txt1,'Color','red','FontWeight', 'bold');
    text(geom(connec(j+1,3),1)-dhx/3,geom(connec(j+1,3),2)-dhy/4,txt2,'Color','red','FontWeight', 'bold');
end

 

LST单元网格

单元自动网格划分

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

 

平面板模型被划分为若干个小矩形区域，共有9个节点，分别是    ，一个矩形中又进一步划分两个三角形单元，第一个单元的节点为    ，第二个单元的节点为    。

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

网格生成算法代码

global nnd nel geom connec XIG YIG
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 = (2*j-1) + (2*i-2)*(2*NYE+1) ;
            n2 = (2*j-1) + (2*i-1)*(2*NYE+1);
            n3 = (2*j-1) + (2*i)*(2*NYE+1);
            n4 = n1 + 1;
            n5 = n2 + 1;
            n6 = n3 + 1 ;
            n7 = n1 + 2;
            n8 = n2 + 2;
            n9 = n3 + 2;
            %
            geom(n1,:) = [(i-1)*dhx - X_origin             (j-1)*dhy - Y_origin];
            geom(n2,:) = [((2*i-1)/2)*dhx - X_origin       (j-1)*dhy - Y_origin  ];
            geom(n3,:) = [i*dhx - X_origin                 (j-1)*dhy - Y_origin  ];
            geom(n4,:) = [(i-1)*dhx - X_origin           ((2*j-1)/2)*dhy - Y_origin ];
            geom(n5,:) = [((2*i-1)/2)*dhx - X_origin     ((2*j-1)/2)*dhy - Y_origin ];                         
            geom(n6,:) = [i*dhx - X_origin               ((2*j-1)/2)*dhy - Y_origin ];
            geom(n7,:) = [(i-1)*dhx - X_origin              j*dhy - Y_origin];
            geom(n8,:) = [((2*i-1)/2)*dhx - X_origin        j*dhy - Y_origin];
            geom(n9,:) = [i*dhx - X_origin                  j*dhy - Y_origin];
            %
            nel = 2*k;
            m = nel -1;
            connec(m,:) = [n1  n2  n3    n5    n7   n4];
            connec(nel,:) = [n3   n6    n9    n8    n7   n5];
            max_n = max([n1  n2  n3   n4  n5  n6  n7   n8    n9]);
            if(nnd <= max_n); nnd = max_n; end;
    end
end

代码解读

1. 先确定节点号      、      、      ，n4 = n1 + 1、n5 = n2 + 1、n6 = n3 + 1、n7 = n1 + 2、n8 = n2 + 2、n9 = n3 + 2说明      、      、      、      、      、      在      、      、      的基础上，编码加1、加2；
2. n1 = (2*j-1) + (2*i-2)*(2*NYE+1)行不动，每次按照列增加，说明      按照纵向排序，注意的是，6节点单元边的中心也有节点，所以是(2*j-1)；
3. n2 = (2*j-1) + (2*i-1)*(2*NYE+1)比      多了一列的节点，说明      与      同行；
4.      的坐标geom由相对位置-坐标轴原点(X_origin,Y_origin)，该数值由主程序中给出；
5. nel = 2*k指的是每次循环中矩形个数的2倍，当两层循环结束时，nel指的是全部三角形单元的个数；
6. m = nel -1指的是每次循环中矩形个数的2倍-1；
7. connec(m,:)指的是矩形区域下面那个三角形单元节点编码顺序，connec(nel,:)指的是矩形区域下面那个三角形单元节点编码顺序；
8. connec(m,:) = [n1 n2 n3 n5 n7 n4]指的是该单元由      节点组成，并不是固定的，按照一定的顺序即可（顺时针或者逆时针）；
9. if(nnd <= max_n); nnd = max_n; end;当循环结束时，max([n1 n2 n3 n4 n5 n6 n7 n8 n9])的数值就是节点的最大值，也就是节点的个数。

网格绘制

figure('Name','LST单元有限元网格模型');
patch('Faces', connec, 'Vertices', geom,  'Facecolor','c','Marker','o','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:2:nel
    txt1 = num2str(j);
    txt2 = num2str(j+1);
    text(geom(connec(j,1),1)+dhx/3,geom(connec(j,1),2)+dhy/4,txt1,'Color','red','FontWeight', 'bold');
    text(geom(connec(j+1,3),1)-dhx/3,geom(connec(j+1,3),2)-dhy/4,txt2,'Color','red','FontWeight', 'bold');
end

 

代码解读

1. 对于二维网格绘制过程中patch函数中的Faces对应于单元节点组成数组connec，Vertices对应于节点信息数组geom;
2. Facecolor决定面的颜色，c表示的是颜色代码，Matlab支持的颜色代码如下：        
      
3. Marker是Matlab标记符号的命令，可选择多种标记，Matlab支持的标记如下：    
4. MarkerFaceColor可指定标记的面颜色，相应的颜色代码，上图已写明；
5. 编写了显示节点、单元编码程序段，遍历每个节点，num2str(i)将数值      转换为字符串类型；
6. 使用text命令将节点号标记在图中指定位置(geom(i,1)+dhx/10,geom(i,2)+dhy/10);
7. text命令中的Color、FontWeight分别指定文本的颜色和字体粗细。