前言
本文是专栏《有限单元法MATLAB程序设计》的第一篇章，将从“程序设计”、“程序改进”、“APP开发”这循序渐进的三部分来详细阐述“平面问题常应变三角形MATLAB分析计算”这一主题。本文的编程思路基于江文强的译著《MATLAB和Abqus有限元分析理论与应用》，笔者也在对该书学习的同时，形成了自己的MATLAB APP开发思路，即：
1.理论的学习与手算的实践
2.利用实时编辑器将“理论与手算”转化为编程语言
3.利用编程语言的优势完善程序（如编写通用函数、丰富结果可视化等）
4.利用MATLAB App Designer进行“UI设计”与“功能迁移”
本篇章以下三节内容详细阐述“平面问题常应变三角形MATLAB分析计算”：
一、平面问题常应变三角形单元(CST)MATLAB编程
二、矩形边界常应变三角形单元(CST)自动划分网格MATLAB编程
三、平面问题常应变三角形单元(CST)MATLAB APP的开发实践
正文
1 例题
假设C=10mm，L=60mm，t=5mm，弹性模量为200000MPa，泊松比为0.3，集中荷载P为1000N。如下图
用24个单元将该区域离散，其中编号为19、20、21的三个结点为固定端。单元划分如下图所示
1 编程思路概述
本程序分为三部分：数据文件、主程序、函数。其中，数据文件用来输入求解随想的基本信息；主程序用来搭建整体有限元编程计算前处理、分析、后处理的框架；函数用来将一些通用求解部分归纳整合起来，方便其他程序调用。
注：笔者的另一篇文章《关于杆系结构MATLAB编程计算的一些思考》将MATLAB有限元所有求解流程均实现了通用函数化， 使得主文件可读性更强，通用效果更好。读者可移步至该文章查看优化思路。
编程思路框架如图所示：
2 程序架构
2.1 数据文件
数据文件输入结点数、单元数、结点自由度数；节点坐标、单元中结点间关系、单元材料和几何参数、边界条件、节点载荷。
%% 本文件为常应变三角形单元数据输入文件（横梁）
global nnd nel nne nodof eldof n geom connec dee nf Nodal_loads
format short e
%%%%%%%%%%%%%%%%%%%%%%% 输入开始 %%%%%%%%%%%%%%%%%%%%%%
%% 节点&单元数
nnd = 21 ; % 节点数
nel = 24 ; % 单元数
nne = 3 ; % 单元节点数
nodof =2; % 节点自由度
eldof = nne*nodof; % 单元自由度
%% 节点坐标
% geom(nnd,2)存储节点的x,y坐标
geom = [ 0, -10 ; % Node 1
0, 0 ; % Node 2
0, 10 ; % Node 3
10, -10 ; % Node 4
10, 0 ; % Node 5
10, 10 ; % Node 6
20, -10 ; % Node 7
20, 0 ; % Node 8
20, 10 ; % Node 9
30, -10 ; % Node 10
30, 0 ; % Node 11
30, 10 ; % Node 12
40, -10 ; % Node 13
40, 0 ; % Node 14
40, 10 ; % Node 15
50, -10 ; % Node 16
50, 0 ; % Node 17
50, 10 ; % Node 18
60, -10 ; % Node 19
60, 0 ; % Node 20
60, 10 ]; % Node 21
%% 节点关联关系
%connec(nel,2)存储节点关联关系。【节点关联关系矩阵】
connec = [ 1, 4, 2; ...% Element 1
4, 5, 2; ...% Element 2
2, 5, 3; ...% Element 3
5, 6, 3; ...% Element 4
4, 7, 5; ...% Element 5
7, 8, 5; ...% Element 6
5, 8, 6; ...% Element 7
8, 9, 6; ...% Element 8
7, 10, 8; ...% Element 9
10, 11, 8; ...% Element 10
8, 11, 9; ...% Element 11
11, 12, 9; ...% Element 12
10, 13, 11; ...% Element 13
13, 14, 11; ...% Element 14
11, 14, 12; ...% Element 15
14, 15, 12; ...% Element 16
13, 16, 14; ...% Element 17
16, 17, 14; ...% Element 18
14, 17, 15; ...% Element 19
17, 18, 15; ...% Element 20
16, 19, 17; ...% Element 21
19, 20, 17; ...% Element 22
17, 20, 18; ...% Element 23
20, 21, 18]; % Element 24
%% 单元材料和几何参数
E = 200000.; % E（Mpa）
vu = 0.3; % 泊松比
thick = 5.; % 横梁厚度（mm）
%% 形成平面应力的弹性矩阵
dee = formdsig(E,vu);
%% 边界条件
%nf(nnd,nodof)储存边界条件。0代表自由度被约束，1（编号）代表没有。
%%1.假设每个节点每个自由度都未被约束。
nf = ones(nnd, nodof); % 将矩阵 nf 初始化为 1
%%2.对实际有约束的自由度设置为0
nf(19,1) = 0; nf(19,2) = 0; % 节点19、20、21为固端约束，自由度0
nf(20,1) = 0; nf(20,2) = 0;
nf(21,1) = 0; nf(21,2) = 0;
%%3.对未被约束的节点进行编号
n=0;
for i=1:nnd
for j=1:nodof
if nf(i,j) ~= 0
n=n+1;
nf(i,j)=n;
end
end
end
%% 节点载荷
%load矩阵储存荷载
Nodal_loads= zeros(nnd, 2);
Nodal_loads(2,1) = 0.; Nodal_loads(2,2) = -1000.; % Node 2
%%%%%%%%%%%%%%%%%%%%%% 输入结束 %%%%%%%%%%%%%%%%%%%%对于该数据输入文件，仍有优化的空间，如：结点数nnd可由结点坐标矩阵geom的行数大小求得；单元数nne可由单元结点关系矩阵connec求得；边界条件的处理可编写成通用函数等。该优化方式可见文章：《关于杆系结构MATLAB编程计算的一些思考》。
2.2 主程序
主程序分为：前处理、分析、后处理三个部分。具体代码如下：
%% 本程序使用3节点常应变三角形单元对二维问题进行线弹性静态分析
clear;clc
global nnd nel eldof n geom connec dee nf Nodal_loads
format long g
%% 前处理
% 1参数的加载
fid =fopen('CST_COARSE_MESH_RESULTS.txt','w');%命名计算结果文件
CST_COARSE_MESH_DATA;%载入数据件
% 2整体力矢的组装
fg=zeros(n,1);
for i=1: nnd
if nf(i,1) ~= 0
fg(nf(i,1))= Nodal_loads(i,1);
end
if nf(i,2) ~= 0
fg(nf(i,2))= Nodal_loads(i,2);
end
end
% 3总刚的组装
kk = zeros(n, n);
for i=1:nel
[bee,g,A] = elem_T3(i); % 形成应变矩阵bee，操作矢量g，单元面积A
ke=thick*A*bee'*dee*bee; % 计算单刚
kk=form_KK(kk,ke, g); % 组装总刚
end
%% 分析
%1 整体有限元方程的求解
delta = kk\fg ; % 求解未知位移
node_disp=zeros(nnd,2);
for i=1: nnd %
if nf(i,1) == 0
x_disp =0.;
else
x_disp = delta(nf(i,1));
end
if nf(i,2) == 0
y_disp = 0.;
else
y_disp = delta(nf(i,2));
end
node_disp(i,:) =[x_disp y_disp];
end
%2 单元应力和应变
%1)检索位于中性轴上的节点的x_coord和y_disp
k = 0;
for i=1:nnd
if geom(i,2)== 0.
k=k+1;
x_coord(k) = geom(i,1);
vertical_disp(k)=node_disp(i,2);%vertical_disp垂直位移
end
end
%2)计算应变应力
for i=1:nel
[bee,g,A] = elem_T3(i); % 形成应变矩阵bee，操作矢量g，单元面积A
eld=zeros(eldof,1); % 将单元位移初始化为零
for m=1:eldof
if g(m)==0
eld(m)=0;
else
eld(m)=delta(g(m)); % 检索单元位移
end
end
eps=bee*eld; % 计算应变（应变矩阵*位移）
EPS(i,:)=eps ; % 存储所有结点的应变
sigma=dee*eps; % 计算应力
SIGMA(i,:)=sigma ; % 存储所有节点的应力
end
%% 后处理
print_CST_results;
% 绘制 x_direction 上的应力（梁的纵向应力云图）
x_stress = SIGMA(:,1);
cmin = min(x_stress);
cmax = max(x_stress);
caxis([cmin cmax]) %设置颜色范围
patch('Faces', connec, 'Vertices', geom, 'FaceVertexCData',x_stress, 'Facecolor','flat','Marker','o')
%patch('Faces',F,'Vertices',V，'FaceVertexCData',x_stress)
% 创建一个或多个多边形，其中 V 指定顶点的值，F 定义要连接的顶点
% 'FaceVertexCData'这个东西是来指定顶点vertex上的颜色的。
% 'Facecolor','flat'设置属性
title('\sigma_{xx}(KN/mm^{2})','fontname','Times New Roman')
colorbar2.3 函数
2.3.1 单元分析elem_T3.m
function[bee,g,A] = elem_T3(i)
%% 该函数返回单元 i 的节点坐标、单元面积、应变矩阵、操作向量
global nne nodof geom connec nf
%% 返回节点坐标
x1 = geom(connec(i,1),1); y1 = geom(connec(i,1),2);
x2 = geom(connec(i,2),1); y2 = geom(connec(i,2),2);
x3 = geom(connec(i,3),1); y3 = geom(connec(i,3),2);
%% 单元面积的计算
A = (0.5)*det([1 x1 y1;
1 x2 y2;
1 x3 y3]);
m11 = (x2*y3 - x3*y2)/(2*A);
m21 = (x3*y1 - x1*y3)/(2*A);
m31 = (x1*y2 - y1*x2)/(2*A);
m12 = (y2 - y3)/(2*A);
m22 = (y3 - y1)/(2*A);
m32 = (y1 - y2)/(2*A);
m13 = (x3 - x2)/(2*A);
m23 = (x1 - x3)/(2*A);
m33 = (x2 -x1)/(2*A);
%% 应变矩阵的计算
bee = [ m12 0 m22 0 m32 0; ...
0 m13 0 m23 0 m33; ...
m13 m12 m23 m22 m33 m32] ;
%% 操作向量
l=0;
for k=1:nne
for j=1:nodof
l=l+1;
g(l)=nf(connec(i,k),j);
end
end2.3.2 总刚的组装函数form_KK.m
function[kk]=form_KK(kk,ke, g)
%% 此函数组装总刚
global eldof
for i=1:eldof
if g(i) ~= 0
for j=1: eldof
if g(j) ~= 0
kk(g(i),g(j))= kk(g(i),g(j)) + ke(i,j);
end
end
end
end
2.3.3 平面应力问题的弹性矩阵的构件函数formdsig.m
function[dee] = formdsig(E,vu)
global vu E
% 此函数形成平面应力问题的弹性矩阵
c=E/(1.-vu*vu);
dee=c*[1 vu 0. ;...
vu 1 0. ;...
0. 0. .5*(1.-vu)];2.3.4 结果打印函数CST_PLANE_STRESS.m
%% 本文件把分析结果输出到fid所代表的文件里。
fprintf(fid, '-------------------------------------------------------- \n');
fprintf(fid, ' \n ********** 分析结果 **********\n');
%% 打印节点位移
fprintf(fid, '------------------------------------------------------ \n');
fprintf(fid, 'Nodal displacements \n');
fprintf(fid, 'Node disp_x disp_y \n');
for i=1:nnd
fprintf(fid,' %g, %8.5e, %8.5e\n', ...
i, node_disp(i,1), node_disp(i,2));
end
fprintf(fid,'\n');
%% 打印节点应力
fprintf(fid, '------------------------------------------------------ \n');
fprintf(fid, ' Element stresses \n');
fprintf(fid, 'element sigma_(xx) sigma_(yy) tau_(xy)\n');
for i=1:nel
fprintf(fid,' %g, %7.4e, %7.4e, %7.4e\n',i, ...
SIGMA(i,1),SIGMA(i,2),SIGMA(i,3));
end
%% 打印节点应变
fprintf(fid, '------------------------------------------------------ \n');
fprintf(fid, ' Element strains \n');
fprintf(fid, 'element epsilon_(xx) epsilon_(yy) gamma_(xy)\n');
for i=1:nel
fprintf(fid,' %g, %7.4e, %7.4e, %7.4e\n',i, ...
EPS(i,1),EPS(i,2),EPS(i,3));
end3 运行结果
3.1 可视化结果
平面主应力可视化结果如下图所示：
可视化结果分析：
常应变三角形单元由于应力及应变矩阵为常量，且单元间应力梯度较大，分析结果较粗糙。
3.2 结果打印
--------------------------------------------------------
********** 分析结果 **********
------------------------------------------------------
Nodal displacements
Node disp_x disp_y
1, 1.45081e-02, -6.49329e-02
2, 3.28049e-04, -6.52078e-02
3, -1.42385e-02, -6.47141e-02
4, 1.42332e-02, -4.97317e-02
5, 1.82950e-04, -4.94530e-02
6, -1.38358e-02, -4.94091e-02
7, 1.29745e-02, -3.50495e-02
8, 1.37982e-04, -3.46630e-02
9, -1.26721e-02, -3.47556e-02
10, 1.09224e-02, -2.19922e-02
11, 8.95233e-05, -2.14870e-02
12, -1.07002e-02, -2.16958e-02
13, 8.08085e-03, -1.13485e-02
14, 2.56420e-05, -1.07261e-02
15, -7.90991e-03, -1.10480e-02
16, 4.46383e-03, -3.88383e-03
17, -6.63586e-05, -3.19069e-03
18, -4.26507e-03, -3.66370e-03
19, 0.00000e+00, 0.00000e+00
20, 0.00000e+00, 0.00000e+00
21, 0.00000e+00, 0.00000e+00
------------------------------------------------------
Element stresses
element sigma_(xx) sigma_(yy) tau_(xy)
1, -7.8546e+00, -7.8546e+00, 7.8546e+00
2, -1.3515e+00, 5.1683e+00, 1.3112e+01
3, 6.6118e-02, 9.8937e+00, 9.1400e+00
4, 9.1400e+00, 3.6192e+00, 9.8937e+00
5, -2.5827e+01, -2.1744e+00, 4.8607e+00
6, 1.5601e+00, 8.1980e+00, 1.5027e+01
7, -6.9913e-01, 6.6741e-01, 5.9323e+00
8, 2.4966e+01, 5.6374e+00, 1.4180e+01
9, -4.2552e+01, -5.0356e+00, 1.6983e+00
10, 2.2662e+00, 1.0785e+01, 1.8024e+01
11, -1.6757e+00, -2.3552e+00, 2.8152e+00
12, 4.1961e+01, 8.4119e+00, 1.7462e+01
13, -5.9121e+01, -7.6315e+00, -1.4550e+00
14, 2.6997e+00, 1.3258e+01, 2.0813e+01
15, -2.7809e+00, -5.0108e+00, -2.2163e-01
16, 5.9202e+01, 1.1322e+01, 2.0864e+01
17, -7.5391e+01, -1.0170e+01, -4.5429e+00
18, 2.5481e+00, 1.4627e+01, 2.3117e+01
19, -4.1445e+00, -7.6816e+00, -3.0783e+00
20, 7.6988e+01, 1.3636e+01, 2.4504e+01
21, -9.3536e+01, -1.4198e+01, -4.9720e+00
22, 1.4584e+00, 4.3753e-01, 2.4544e+01
23, -1.6603e+00, -9.9582e+00, -7.7540e+00
24, 9.3738e+01, 2.8121e+01, 2.8182e+01
------------------------------------------------------
Element strains
element epsilon_(xx) epsilon_(yy) gamma_(xy)
1, -2.7491e-05, -2.7491e-05, 1.0211e-04
2, -1.4510e-05, 2.7869e-05, 1.7045e-04
3, -1.4510e-05, 4.9369e-05, 1.1882e-04
4, 4.0271e-05, 4.3858e-06, 1.2862e-04
5, -1.2587e-04, 2.7869e-05, 6.3189e-05
6, -4.4967e-06, 3.8650e-05, 1.9535e-04
7, -4.4967e-06, 4.3858e-06, 7.7120e-05
8, 1.1637e-04, -9.2623e-06, 1.8434e-04
9, -2.0521e-04, 3.8650e-05, 2.2078e-05
10, -4.8459e-06, 5.0524e-05, 2.3431e-04
11, -4.8459e-06, -9.2623e-06, 3.6597e-05
12, 1.9719e-04, -2.0883e-05, 2.2701e-04
13, -2.8416e-04, 5.0524e-05, -1.8915e-05
14, -6.3881e-06, 6.2239e-05, 2.7057e-04
15, -6.3881e-06, -2.0883e-05, -2.8812e-06
16, 2.7903e-04, -3.2191e-05, 2.7123e-04
17, -3.6170e-04, 6.2239e-05, -5.9058e-05
18, -9.2001e-06, 6.9314e-05, 3.0052e-04
19, -9.2001e-06, -3.2191e-05, -4.0018e-05
20, 3.6448e-04, -4.7301e-05, 3.1856e-04
21, -4.4638e-04, 6.9314e-05, -6.4636e-05
22, 6.6359e-06, 0.0000e+00, 3.1907e-04
23, 6.6359e-06, -4.7301e-05, -1.0080e-04
24, 4.2651e-04, 0.0000e+00, 3.6637e-04
后记
本文是《MATLAB和Abqus有限元分析理论与应用》的学习心得，几乎原样照搬书中提供的代码。这是笔者学习的基础，是笔者打开MATLAB有限元分析的第一步。正是有了“仿写”的基础，才有后来的“应用”、“改进”与“迁移”。
后两期将从“程序改进——矩形边界网格自动化分”、“MATLAB APP的开发实践”两部分阐述该项目的完善，也记录着笔者的学习历程，即从“仿写”到“改进”与“迁移”。
扫一扫
5670
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
