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

知乎、B站[易木木响叮当]
关注可了解更多的有限元数值仿真技巧。问题或建议,请公众号留言;
如果你觉得木木同学对你有帮助,欢迎赞赏。

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

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

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

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

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

代码获取:


Q4单元网格

单元自动网格划分

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

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

平面板模型被划分为若干个小矩形区域,共有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 + 1n4 = 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
基于Matlab的有限元网格自动生成算法 | Q4、Q8、Abaqus单元网格的图2
Q4单元网格生成(纵向排序)

网格生成算法代码(横向排序)

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

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 + 1n4 = 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的有限元网格自动生成算法 | Q4、Q8、Abaqus单元网格的图3
Q4单元网格生成(横向排序)

绘图修饰

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

基于Matlab的有限元网格自动生成算法 | Q4、Q8、Abaqus单元网格的图4
修改绘图细节

Abaqus-Q4单元

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

Abaqus-Q4单元节点编码顺序为什么是1-2-4-3,而不是按照顺时针或逆时针的顺序编排

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

基于Matlab的有限元网格自动生成算法 | Q4、Q8、Abaqus单元网格的图5
Abaqus-CPS4单元节点顺序

模型节点读取

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

基于Matlab的有限元网格自动生成算法 | Q4、Q8、Abaqus单元网格的图6
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
基于Matlab的有限元网格自动生成算法 | Q4、Q8、Abaqus单元网格的图7
Abaqus网格信息导入进Matlab

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

Q8单元网格

单元自动网格划分

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

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

平面板模型被划分为若干个小矩形区域,共有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 + 1n5 = n3 + 2n6 = n2 + 1n7 = n1 + 2n8 = n1 + 1说明 的基础上,编码加1、加2;
  2. 注意 n1n2n3的编码结构,与上面示意图相符;
  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的有限元网格自动生成算法 | Q4、Q8、Abaqus单元网格的图9
Q8单元网格自动生成

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

往期推荐


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


本次分享仅限于此了,欢迎大家点赞收藏转发!


谢谢你看完木木同学的分享,今日份阅读花费的流量+1M哈哈哈哈哈哈基于Matlab的有限元网格自动生成算法 | Q4、Q8、Abaqus单元网格的图11



-End-


易木木响叮当

想陪你一起度过短暂且漫长的科研生活


有限元分析网格划分

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

基于Matlab的有限元网格自动生成算法 | Q4、Q8、Abaqus单元网格的相关案例教程

知乎、B站:[易木木响叮当] 关注可了解更多的有限元数值仿真技巧。问题或建议,请公众号留言; 如果你觉得木木同学对你有帮助,欢迎赞赏。 今日给大家带来的主要内容是二维问题下有限元网格如何自动生成? 单元网格的形成实际上属于有限元计算中的前处理部分,即确定单元节点信息,当模型较为复杂时,用户可在Abaqus、Ansys等大型商业有限元软件中进行建模,导出网格信息。 当模型较为简单时,如二维平面板模型
1.改写输入数据格式,使之能适应任意几何(可利用节点坐标输 入节点,利用单元-节点关系输入单元); 2. 计算节点应力,给出并实现至少一种应力处理方案,提供处理 前后的应力结果(可用表格和云图表示),可与其它软件对比; 3.提交总结报告(包括方法/方案描述、带详细注释的代码、程序框图、算例描述、结果比较分析等)、可编译源代码、可执行文件、 数据文件、结果文件 program p53 !------
MATLAB启动ABAQUS,提取odb,将ABAQUS模型视作MATLAB的一个函数,用于优化分析 inpName = 'Job-1'; workDir = 'G:\CABLE_CON\TEXT2'; runabaqus(workDir,inpName); % 运行py脚本文件 system(['abaqus viewer noGui=','out_F.py']); % 载入数据 t=load(
MATLAB启动ABAQUS,提取odb,将ABAQUS模型视作MATLAB的一个函数,用于优化分析 inpName = 'Job-1'; workDir = 'G:\CABLE_CON\TEXT2'; runabaqus(workDir,inpName); % 运行py脚本文件 system(['abaqus viewer noGui=','out_F.py']); % 载入数据 t=load(
1、前言 考虑材料参数空间变异性的岩土工程对象的数值分析是岩土工程研究中重要分支。当前,考虑材料参数空间变异性(即参数随机场)的分析手段中,除了极少数非主流的数值分析软件可以实现一键式随机场分析以外,大多数复杂的随机场实现都存在较高的门槛,且难以实现复杂的岩土对象相互作用分析。在主流岩土工程分析软件中,如,flac和abaqus中,前者需要使用fish编程,且基本要借助第三方软件计算随机场才能实现
影响力
粉丝
内容
获赞
收藏
    项目客服
    培训客服
    6 8