基于有限元计算程序的探究(上)—网格划分、刚度矩阵求解、载荷施加

前言  

在我们从学校踏入职场，从事有限元计算的初期，绝大多数人都是一拿到模型就导入商业软件中，一顿操作，洋洋洒洒从几何处理、网格划分、加载计算，最后看到五颜六色的云图，颔首微笑，自信满满，以为完全掌握了有限元计算这项技能，打算在该领域大显身手、扬名立万。随着从业时间的增长，开始的几年里进步神速，真正实现了“从入门到精通”，但是在以后几年中，还是在“原地踏步”，成长性不够。这可能就是太依赖于商业软件导致的，商业软件是一把双刃剑，给予我们方便的同时，也扼杀了我们探求其内在算法的欲望，没有理论上升的空间。在疲于应付项目的过程中，不经意间可能已经忘记了有限元计算的理论基础：什么是形函数？如何用节点位移来表示单元内任意一点的位移？单元刚度矩阵如何“组装”成整体刚度矩阵？等等。因此在闲暇时刻，有必要翻一翻理论教材，细细品味一下其中的奥秘，对我们真正理解“有限元”是有帮助的。这里推荐周博老师的《有限元法与MATLAB》著作，它是一本非常好的书，在讲解基础知识的同时，也用Matlab程序实现了各知识点的贯穿，帮助我们更好的理解一些“生涩”的词汇。

1. 网格划分    

有限元计算的第一步，就是划分网格，其函数为：  

[Nxy, Enod]=mesh2d3n(x12，y12，m，n)  

源程序：    

function[Nxy,Enod] = mesh2d3n(x12,y12, m, n)  

x =linspace(x12(1),x12(2),m);

y =linspace(y12(1),y12(2),n);

[X,Y] =meshgrid(x,y);

Enod =delaunay(X,Y); %%%生成3结点单元结点信息矩阵

h1 =triplot(Enod,X,Y); %%%绘制三角形单元的离散网格

set(h1,'color','k')

axisequal

[en,ny]= size(Enod); %%%en为单元数，ny为每个单元节点数3

for i =1:en

     x_sum = 0;

     y_sum = 0;

     xy_sum = 0;

           for j = 1:ny

               ID = Enod(i,j);

               x_sum = x_sum + X(ID);

           end

       xc = x_sum/ny;    %%%计算三角形单元形心横坐标

       yc = y_sum/ny;    %%%计算三角形单元形心纵坐标

       h2 = text(xc,yc,num2str(i)); %%%在三角形形心处标注单元编号

       set(h2,'color','b')

       set(h2,'fontsize',10)

end

for k =1:m*n

     Nxy(k,1:3) = [k, X(k), Y(k)];

     h3 = text(X(k), Y(k),num2str(k)); %%%在节点坐标处标注节点编号

    set(h3,'color','r')

    set(h3,'fontsize',10)

end

Enod =[(1:en)', Enod]; %%%把单元编号放在第1列，2-4列为节点编号

end

clear；clc

x12=[0,3];

y12=[0,3];

[Nxy,Enod] = mesh2d3n(x12,y12, m, n);

Axis([-0.5,3.5,-0.5,3.5]);

图1 网格划分结果

A =0.5*det(A); %%%三角形单元面积

求得单元应变矩阵后，再求单元刚度矩阵，其函数为：

KE=EstiffM2d3n(xy，mat)  

参数意义：  

xy：3行2列矩阵， 1-2列分别是单元节点的横、纵坐标；  

mat：1行3列矩阵，1-3列分别是材料弹性模量、泊松比、单元厚度；  

KE：6行6列单元刚度矩阵  

源程序：  

function KE = EstiffM2d3n(xy, mat)  

函数应用实例

运行程序后，可以得到这两个单元的刚度矩阵，见图2。

图2 计算得到的单元刚度矩阵

3.  整体刚度矩阵  

KS=SstiffM2d3n(Nxy，Enod，Emat)  

参数意义：  

Nxy： N(节点总数)行3列矩阵， 1-3列分别为节点编号、节点横、纵坐标；  

Enod：M(单元总数)行4列矩阵，1-4列分别为单元编号、单元所含3个节点编号。

Emat：M(单元总数)行4列矩阵，1-4列分别为单元编号、弹性模量、泊松比、单元厚度

源程序：  

functionKS = SstiffM2d3n(Nxy, Enod, Emat)  

M =size(Enod, 1);%%%单元总数量

KS =zeros(2*N, 2*N); %%% 整体刚度矩阵初始化

for j =1:M

     ii = Enod(j, 2);  %%% 依次提取各单元的节点编号

     jj = Enod(j, 3);

    mm = Enod(j, 4);

    sn = [ii, jj, mm];  %%% 节点编号矩阵

    xy = Nxy(sn,2:3);  %%% 提取节点坐标

    mat = Emat(j, 2:4);  %%% 提取材料属性

    KE = EstiffM2d3n(xy, mat);   %%% 调用EstiffM2d3n生成单元刚度矩阵

    sn = [2*ii-1, 2*ii, 2*jj-1, 2*jj, 2*mm-1,2*mm];   %%% 对号入座数组

    KS(sn, sn) = KS(sn, sn) + KE;   %%%将单元刚度矩阵累加到整体刚度矩阵

end

end  

函数应用实例

clear;clc;

x12= [0,2];

y12= [0,1.5];

[Nxy,Enod]= mesh2d3n(x12,y12, 3, 3); %%%结构的三角形网络离散

axis([-0.2,2.2,-0.2,1.7])

[M,N] = size(Enod);  %%%M为单元个数

Emat(1,1:4)= [1, 180e9, 0.35,  5e-2];  %%%材料属性

fori = 2:M

     Emat(i, 1:4) = Emat(1,1:4); %%%对每一个单元赋材料属性

end  

KS =SstiffM2d3n(Nxy, Enod, Emat) %%%生成整体刚度矩阵

运行后，整体刚度矩阵如图3。

图3 整体刚度矩阵

4. 载荷的施加  

需要根据载荷施加情况，需要计算等效节点载荷，其函数为：  

SP=SloadA2d3n(EP，N)  

EP：4列矩阵，1-4列分别为单元编号、节点编号、载荷方向、载荷大小；  

N：节点总数  

SP：2N行1列节点载荷矩阵  

源程序：    

functionSP = SloadA2d3n(EP,N)

[n,m]=size(EP);%%% 提取载荷数量n

SP(1:2*N,1)= 0;%%% 载荷矩阵初始化

fori=1:n

      sn = 2*EP(i,2)+EP(i,3)-2;  %%%施加载荷的位置与方向

      SP(sn) = SP(sn) + EP(i,4);  %%%施加载荷的大小

End  

函数应用实例  

clear;clc;  

N= 9;

EP= [7, 8, 1, 3.5e3；

7, 6, 2, -(2+2/3)*1e3；

7, 9, 2, -(2+4/3)*1e3；

5, 3, 2, -2/3*1e3；

5, 6, 2, -4/3*1e3]；%%%受载荷信息：1-4列分别为单元编号、节点编号、载荷方向、 载荷大小；

SP= SloadA2d3n(EP,N)；%%%载荷列阵  

运行后，得到载荷列阵，结果如图4。  

图4 载荷列阵