基于有限元计算程序的探究(下)—方程求解、单元应变应力计算、后处理
前言
在上一篇文章中，已经介绍了网格的划分、刚度矩阵的求解与载荷的施加；本篇文章是在上一篇文章的基础上，讲解后续操作，主要包括：方程的求解、单元应变应力的计算与后处理（包括变形图与云图）。
(为了上、下两篇文章保持连续性，标题、图号均续上篇)
5．方程的求解
该部分主要是求解节点位移和约束反力，其函数为：
[Ndsp，Rfoc]=SloveS2d3n(KS，SP，cons)
参数意义：
KS：2N(N为节点总数)行2N列整体刚度矩阵；
SP：2N行1列节点载荷矩阵；
cons：RN (位移约束总数) 行3列矩阵，1-3列分别为节点编号、约束位移方向、约束位移值；
Ndsp：N行3列矩阵，1-3列分别为节点编号、x方向位移值、y方向位移值；
Rfoc：RN行3列矩阵，1-3列分别为节点编号、约束反力方向、约束反力数值。
function [Ndsp,Rfoc] = SloveS2d3n(KS, SP, cons)
函数应用实例
运行完程序后，得到所有节点位移值和约束节点的约束反力值，如图5所示。
图5 方程求解结果
6. 单元应变、应力的计算
通过方程求解得到各节点的位移值后，可以得到各单元的应变应力值，其函数为：
[Estrain，Estress]=EstssM2d3n(Nxy，Enod，Emat，Ndsp)
参数意义：
源程序
函数应用实例
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)
Emat(1,1:4) = [1,180e9, 0.35, 5e-2];
for i = 2:M
Emat(i, 1:4) = Emat(1,1:4);
end
KS = SstiffM2d3n(Nxy,Enod, Emat); %%%整体刚度矩阵
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]; %%%载荷施加信息表
[N,M]= size(Nxy)
SP = SloadA2d3n(EP,N); %%%求解载荷列阵
cons = [1,1,0;
1,2,0;
4,2,0;
7,2,0]
[Ndsp, Rfoc] =SloveS2d3n(KS, SP, cons); %%%方程求解
[Estrain,Estress]= EstssM2d3n(Nxy,Enod,Emat,Ndsp); %%% 求解单元应力应变
运行完程序后，得到所有单元的应变、应力值，如图6所示。
图6 单元应变应力7. 后处理
得到节点的位移值后，可以绘制出各单元的变形图，其函数为：
pltdfm2d3n(Enod，Nxy，Ndsp，enlarge)
参数意义：
Enod：M(单元总数)行4列矩阵，1-4列分别为单元编号、单元所含3节点编号；
Nxy： N(节点总数)行3列矩阵， 1-3列分别为节点编号、节点横、纵坐标；
Ndsp：N行3列矩阵，1-3列分别为节点编号、x方向位移值、y方向位移值；
enlarge：变形放大系数。
源程序
function pltdfm2d3n(Enod,Nxy,Ndsp,enlarge)
函数应用实例
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)
Emat(1,1:4) = [1,180e9, 0.35, 5e-2];
for i = 2:M
Emat(i, 1:4) = Emat(1,1:4); %%%为每个单元赋属性
end
KS = SstiffM2d3n(Nxy,Enod, Emat); %%%整体刚度矩阵
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]; %%%求解载荷列阵
[N,M]= size(Nxy)
SP = SloadA2d3n(EP,N); %%%求解载荷列阵
cons = [1,1,0;
1,2,0;
4,2,0;
7,2,0];
[Ndsp, Rfoc] =SloveS2d3n(KS, SP, cons); %%%求解方程
Nvar(:,1:2) =Ndsp(:,1:2); %%%%节点位移
enlarge = 1e5; %%%变形放大系数
pltdfm2d3n(Enod,Nxy,Ndsp,enlarge); %%%画出单元变形图
运行完程序后，得到单元变形前、变形后的图形，如图7所示。
图7 单元变形前后图形对比
除了得到单元变形图之外，还可以得到单元变形云图，其函数为：
pltvc2d3n(Enod，Nxy，Nvar)
参数意义：
Enod：M(单元总数)行4列矩阵，1-4列分别为单元编号、单元所含3节点编号；
Nxy：N(节点总数)行3列矩阵，1-3列分别为节点编号、节点横、纵坐标；
Nvar：N行2列矩阵，1-2列分别为节点编号、节点位移值。
源程序
functionpltvc2d3n(Enod,Nxy,Nvar)
figure
holdon
axisequal
axisoff
m= size(Enod,1)
fori=1:m
k = Enod(i,2:4);
x = Nxy(k,2);
y = Nxy(k,3);
c = Nvar(k,2); %%% 节点位移
fill(x,y,c); %%%填充x和y确定的图形，填充颜色由c来决定
end
colorbar('location','eastoutside')
end
函数应用实例
运行完程序后，得到单元变形后的云图，如图8所示。
图8 单元变形云图
至此，结构有限元所有的步骤及相关程序都已列出，大家可以根据自己兴趣来选择阅读和编写程序。
前面叙述较为零散，为了大家阅读方便，现把问题及程序从头到尾再列一遍（大部分都已在前面列出）。
问题描述如图9(摘自周博老师的《有限元法与MATLAB》)。
图9 问题描述
求解程序：
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);
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); %%%整体刚度矩阵（内嵌了单元刚度矩阵函数）
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]; %%%载荷施加信息矩阵
[N,M]=size(Nxy);
SP= SloadA2d3n(EP,N); %%%节点载荷列阵
cons= [1,1,0;
1,2,0;
4,2,0;
7,2,0]; %%%%约束施加信息矩阵
[Ndsp,Rfoc] = SloveS2d3n(KS, SP, cons); %%%求解节点位移和节点约束反力
Nvar(:,1:2)= Ndsp(:,1:2); %%%把所有节点的编号、x方向位移分量提取出来，为下一步画云图所用
enlarge=1e5;
pltdfm2d3n(Enod,Nxy,Ndsp,enlarge);%%%把离散结构变形图画出
pltvc2d3n(Enod,Nxy,Nvar); %%%把各节点的x方向位移分量用云图表示
运行以上程序后，网格划分见图10、单元变形图见图7、单元变形云图见图8。
图10 网格划分
