基于Matlab语言的按平面三角形单元划分的结构有限元程序设计

基于Matlab语言的按平面三角形单元

划分的结构有限元程序设计

专 业： 建筑与土木工程 班 级： 建工研12-2 * 名： *** 学 号： *********

基于Matlab语言的按平面三角形单元划分

结构有限元程序设计

一、 有限单元发及Matlab语言概述

1. 有限单元法

随着现代工业、生产技术的发展，不断要求设计高质量、高水平的大型、复杂和精密的机械及工程结构。为此目的，人们必须预先通过有效的计算手段，确切的预测即将诞生的机械和工程结构，在未来工作时所发生的应力、应变和位移因此，需要寻求一种简单而又精确的数值分析方法。有限单元法正是适应这种要求而产生和发展起来的一种十分有效的数值计算方法。

有限元法把一个复杂的结构分解成相对简单的“单元”，各单元之间通过结点相互连接。单元内的物理量由单元结点上的物理量按一定的假设内插得到，这样就把一个复杂结构从无限多个自由度简化为有限个单元组成的结构。我们只要分析每个单元的力学特性，然后按照有限元法的规则把这些单元“拼装”成整体，就能够得到整体结构的力学特性。

有限单元法基本步骤如下：

(1)结构离散：结构离散就是建立结构的有限元模型，又称为网格划分或单元划分，即将结构离散为由有限个单元组成的有限元模型。在该步骤中，需要根据结构的几何特性、载荷情况等确定单元体内任意一点的位移插值函数。

(2)单元分析：根据弹性力学的几何方程以及物理方程确定单元的刚度矩阵。 (3)整体分析:把各个单元按原来的结构重新连接起来，并在单元刚度矩阵的基础上确定结构的总刚度矩阵，形成如下式所示的整体有限元线性方程：

FK ①

式中，F是载荷矩阵, K是整体结构的刚度矩阵, 是节点位移矩阵。

(4)载荷移置:根据静力等效原理，将载荷移置到相应的节点上，形成节点载荷矩阵。

(5)边界条件处理:对式①所示的有限元线性方程进行边界条件处理。 (6)求解线性方程:求解式①所示的有限元线性方程，得到节点的位移。在该步骤中，若有限元模型的节点越多，则线性方程的数量就越多，随之有限元分析的计算量也将越大。

(7)求解单元应力及应变 根据求出的节点位移求解单元的应力和应变。

(8)结果处理与显示 进入有限元分析的后处理部分，对计算出来的结果进行加工处理，并以各种形式将计算结果显示出。

2. Matlab简介

在用有限元法进行结构分析时，将会遇到大量的数值计算，因而在实用上是一定要借助于计算机和有限元程序，才能完成这些复杂而繁重的数值计算工作。而Matlab是当今国际科学界最具影响力和活力的软件。它起源于矩阵运算，并已经发展成一种高度集成的计算机语言。它提供了强大的科学计算，灵活的程序设计流程，高质量的图形可视化与界面设计，便捷的与其他程序和语言接口的功能。Matlab在各国高校与研究单位起着重大的作用。

“工欲善其事，必先利其器”。如果有一种十分有效的工具能解决在教学与研究中遇到的问题，那么Matlab语言正是这样的一种工具。它可以将使用者从繁琐、无谓的底层编程中解放出来，把有限的宝贵时间更多地花在解决问题中，这样无疑会提高工作效率。 目前，Matlab已经成为国际上最流行的科学与工程计算的软件工具，现在的Matlab已经不仅仅是一个“矩阵实验室”了，它已经成为了一种具有广泛应用前景的全新的计算机高级编程语言了，有人称它为“第四代”计算机语言，它在国内外高校和研究部门正扮演着重要的角色。Matlab语言的功能也越来越强大，不断适应新的要求提出新的解决方法。可以预见，在科学运算、自动控制与科学绘图领域Matlab语言将长期保持其独一无二的地位。为此，本例采用Matlab语言编程，以利用其快捷强大的矩阵数值计算功能。

二、 问题描述

一矩形薄板，一边固定，承受150kN集中荷载，结构简图及按平面三角形单元划分的有限元模型图如下所示。

材料参数：弹性模量E2108KN/m2；泊松比：0.2；薄板厚度2mm。在本例中，所取结构模型及数据主要用于程序设计理论分析，与工程实际无关。

F=150kN1m2my6415324x1m11m21m

3

三、 参数输入：

单元个数NELEM=4 节点个数NNODE=6 受约束边界点数NVFIX=2 节点荷载个数NFORCE=1 弹性模量YOUNG=2e8 泊松比POISS=0.2 厚度THICK=0.002

1 2 6 2 3 4 单元节点编码数组LNODS =2 4 52 5 60 0 1 0 2 0节点坐标数组COORD =

12 1 10 1节点力数组FORCE =[4 0 -150]

1 1 1 约束信息数组FIXED = 6 1 1以上数值数据为程序运行前输入的初始数据，存为“471220580.txt”文本格式，同时必须放在Matlab工作目录下，路径不对程序不能自动读取指定初始文件，运行出错。初始数据文本文件输入格式如下图：

四、 Matlab语言程序源代码：

1. 程序中变量说明

NNODE 单元节点数 NPION 总结点数 NELEM 单元数 NVFIX 受约束边界点数 FIXED 约束信息数组 NFORCE 节点力数 FORCE 节点力数组 COORD 结构节点坐标数组 LNODS 单元定义数组 YOUNG 弹性模量 POISS 泊松比 THICK 厚度

B 单元应变矩阵(3*6) D 单元弹性矩阵(3*3) S 单元应力矩阵(3*6) A 单元面积 ESTIF 单元刚度矩阵 ASTIF 总体刚度矩阵 ASLOD 总体荷载向量 ASDISP 节点位移向量 ELEDISP 单元节点位移向量 STRESS 单元应力 FP1 数据文件指针

2. Matlab语言程序代码

%******************************************************************************

%初始化及数据调用

format short e %设定输出类型 clear %清除内存变量

FP1=fopen('471220580.txt','rt'); %打开输入数据文件，读入控制数据 NELEM=fscanf(FP1,'%d',1), %单元个数（单元编码总数） NPION=fscanf(FP1,'%d',1), %结点个数（结点编码总数） NVFIX=fscanf(FP1,'%d',1) %受约束边界点数 NFORCE=fscanf(FP1,'%d',1), %结点荷载个数 YOUNG=fscanf(FP1,'%e',1), %弹性模量 POISS=fscanf(FP1,'%f',1), %泊松比 THICK=fscanf(FP1,'%f',1) %厚度

LNODS=fscanf(FP1,'%d',[3,NELEM])' %单元定义数组（单元结点号）

%相应为单元结点号（编码）、按逆时针顺序输入

COORD=fscanf(FP1,'%f',[2,NPION])' %结点坐标数组

%坐标：x，y坐标（共NPOIN组）

FORCE=fscanf(FP1,'%f',[3,NFORCE])' %结点力数组（受力结点编号, x方向,y方向） FIXED=fscanf(FP1,'%d',[3,NVFIX])' %约束信息（约束点，x约束，y约束） %有约束为1，无约束为0

%*****************************************************************************

%生成单元刚度矩阵并组成总体刚度矩阵

ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0

%*****************************************************************************

for i=1:NELEM

%生成弹性矩阵D

D= [1 POISS 0; POISS 1 0;

0 0 (1-POISS)/2]*YOUNG/(1-POISS^2)

%*****************************************************************************

%计算当前单元的面积A

A=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2); 1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2); 1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)])/2

%*****************************************************************************

%生成应变矩阵B

for j=0:2 b(j+1)=

COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);

c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1); end

B=[b(1) 0 b(2) 0 b(3) 0;... 0 c(1) 0 c(2) 0 c(3);...

c(1) b(1) c(2) b(2) c(3) b(3)]/(2*A); B1( :,:,i)=B;

%*****************************************************************************

%求应力矩阵S=D*B

S=D*B;

ESTIF=B'*S*THICK*A; %求解单元刚度矩阵

a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号 for j=1:3 for k=1:3

ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)…

=ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2); %根据节点编号对应关系将单元刚度分块叠加到总刚

%度矩阵中

end end end

%*****************************************************************************

%将约束信息加入总体刚度矩阵（对角元素改一法）

for i=1:NVFIX

if FIXED(i,2)==1

ASTIF(:,(FIXED(i,1)*2-1))=0; %一列为零 ASTIF((FIXED(i,1)*2-1),:)=0; %一行为零 ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1; %对角元素为1 end

if FIXED(i,3)==1

ASTIF( :,FIXED(i,1)*2)=0; %一列为零 ASTIF(FIXED(i,1)*2,:)=0; %一行为零 ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; %对角元素为1 end end

%*****************************************************************************

%生成荷载向量

ASLOD(1:2*NPION)=0; %总体荷载向量置零 for i=1:NFORCE

ASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3); end

%*****************************************************************************

%求解内力

ASDISP=ASTIF\\ASLOD' %计算节点位移向量 ELEDISP(1:6)=0; %当前单元节点位移向量 for i=1:NELEM for j=1:3

ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);

%取出当前单元的节点位移向量 end i

STRESS=D*B1(:, :, i)*ELEDISP' %求内力 end

fclose(FP1); %关闭数据文件

五、 程序运行结果

NELEM = 4 NPION = 6 NVFIX = 2 NFORCE = 1 YOUNG =

200000000 POISS =

2.0000e-001 THICK =

2.0000e-003 LNODS =

1 2 6 2 3 4

2 4 5 2 5 6 COORD =

0 0 1 0 2 0 2 1 1 1 0 1 FORCE =

4 0 -150 FIXED =

1 1 1 6 1 1 D =

2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007 A =

5.0000e-001 D =

2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007 A =

5.0000e-001 D =

2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007 A =

5.0000e-001 D =

2.0833e+008 4.1667e+007 0 4.1667e+007 2.0833e+008 0 0 0 8.3333e+007 A =

5.0000e-001 ASDISP =

0 0 -8.0607e-004 -1.5848e-003 -1.0281e-003

-4.4727e-003 1.1937e-003 -4.6947e-003 8.6670e-004 -1.8880e-003 0 0 （说明：以上12个ASDISP输出结果数据表示节点位移向量，即各节点分别在X，Y方向上产生的位移。） i = 1 STRESS =

-1.6793e+005 -3.3586e+004 -1.3207e+005 i = 2 STRESS =

-5.5503e+004 -5.5503e+004 -5.5503e+004 i = 3 STRESS =

5.5503e+004 -4.9526e+004 -9.4497e+004 i = 4 STRESS =

1.6793e+005 -2.7040e+004 -1.7932e+004 （说明：以上四组STRESS输出结果数据表示单元应力SX,SY,SXY，i为单元号。）

六、 ANSYS建模比较分析

利用ANSYS有限元分析软件，完全按照前面Matlab程序设计的各项条件参数以及单元划分方式建立ANSYS模型，其求解结果与以上程序计算结果比较，验证程序是否正确。

ANSYS模型节点单元如下图所示：

ANSYS模型求解变形图如下所示：

ANSYS求解节点位移结果列表显示如下：（单位：m）

PRINT U NODAL SOLUTION PER NODE

***** POST1 NODAL DEGREE OF FREEDOM LISTING *****

THE FOLLOWING DEGREE OF FREEDOM RESULTS ARE IN THE GLOBAL COORDINATE SYSTEM

NODE UX UY UZ USUM 1 0.0000 0.0000 0.0000 0.0000 2 -0.80607E-03-0.15848E-02 0.0000 0.17780E-02 3 -0.10281E-02-0.44727E-02 0.0000 0.45893E-02 4 0.11937E-02-0.46947E-02 0.0000 0.48441E-02 5 0.86670E-03-0.18880E-02 0.0000 0.20774E-02 6 0.0000 0.0000 0.0000 0.0000 MAXIMUM ABSOLUTE VALUES

NODE 4 4 0 4 VALUE 0.11937E-02-0.46947E-02 0.0000 0.48441E-02

ANSYS求解单元应力结果列表显示如下：（单位：KN/m2）

PRINT S ELEMENT SOLUTION PER ELEMENT

***** POST1 ELEMENT NODAL STRESS LISTING *****

THE FOLLOWING X,Y,Z VALUES ARE IN GLOBAL COORDINATES ELEMENT= 1 PLANE182

NODE SX SY SZ SXY SYZ 1 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 2 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 6 -0.16793E+06 -33586. 0.0000 -0.13207E+06 0.0000 ELEMENT= 2 PLANE182

NODE SX SY SZ SXY SYZ 2 -55503. -55503. 0.0000 -55503. 0.0000 3 -55503. -55503. 0.0000 -55503. 0.0000 4 -55503. -55503. 0.0000 -55503. 0.0000 ELEMENT= 3 PLANE182

NODE SX SY SZ SXY SYZ 2 55503. -49526. 0.0000 -94497. 0.0000 4 55503. -49526. 0.0000 -94497. 0.0000 5 55503. -49526. 0.0000 -94497. 0.0000 ELEMENT= 4 PLANE182

NODE SX SY SZ SXY SYZ 2 0.16793E+06 -27040. 0.0000 -17932. 0.0000

5 0.16793E+06 -27040. 0.0000 -17932. 0.0000 6 0.16793E+06 -27040. 0.0000 -17932. 0.0000

结论

通过比较Matlab语言设计程序运行结果和ANSYS建模分析结果可知，两种方式算出的结果完全一致，说程序设计正确。所以本程序适用于按三角形单元划分的平面结构有限元分析。

format short e clear

FP1=fopen('LinearTriangleElement of Thin plate.txt','rt'); NELEM=fscanf(FP1,'%d',1) NPION=fscanf(FP1,'%d',1) NVFIX=fscanf(FP1,'%d',1) NFORCE=fscanf(FP1,'%d',1) YOUNG=fscanf(FP1,'%e',1) POISS=fscanf(FP1,'%f',1) THICK=fscanf(FP1,'%f',1)

LNODS=fscanf(FP1,'%d',[3,NELEM]) COORD=fscanf(FP1,'%f',[2,NPION]) FORCE=fscanf(FP1,'%f',[3,NFORCE]) FIXED=fscanf(FP1,'%d',[3,NVFIX])

ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0 for i=1:NELEM %生成弹性矩阵D

D= YOUNG/(1-POISS^2)*[1 POISS 0; POISS 1 0;

0 0 (1-POISS)/2] %计算当前单元的面积A

A=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2); 1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2); 1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)])/2 %生成应变矩阵B for j=0:2 b(j+1)=

COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);

c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1); end

B=[b(1) 0 b(2) 0 b(3) 0;... 0 c(1) 0 c(2) 0 c(3);...

c(1) b(1) c(2) b(2) c(3) b(3)]/(2*A); B1( :,:,i)=B; %求应力矩阵S=D*B S=D*B;

ESTIF=B'*S*THICK*A; a=LNODS(i,:); for j=1:3 for k=1:3

ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)…

=ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2

); end end end

%将约束信息加入总体刚度矩阵 for i=1:NVFIX if FIXED(i,2)==1

ASTIF(:,(FIXED(i,1)*2-1))=0; ASTIF((FIXED(i,1)*2-1),:)=0;

ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1; end

if FIXED(i,3)==1

ASTIF( :,FIXED(i,1)*2)=0; ASTIF(FIXED(i,1)*2,:)=0;

ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; end end

%生成荷载向量

ASLOD(1:2*NPION)=0; for i=1:NFORCE

ASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3) end

%求解内力

ASDISP=ASTIF\\ASLOD' ELEDISP(1:6)=0; for i=1:NELEM for j=1:3

ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2); end i

STRESS=D*B1(:, :, i)*ELEDISP' end

fclose(FP1);

format short e clear

FP1=fopen('LinearTriangleElement of Thin plate.txt','rt'); NELEM=fscanf(FP1,'%d',1) NPION=fscanf(FP1,'%d',1) NVFIX=fscanf(FP1,'%d',1) NFORCE=fscanf(FP1,'%d',1) YOUNG=fscanf(FP1,'%e',1) POISS=fscanf(FP1,'%f',1) THICK=fscanf(FP1,'%f',1)

LNODS=fscanf(FP1,'%d',[ NELEM,3]) COORD=fscanf(FP1,'%f',[ NPION,2]) FORCE=fscanf(FP1,'%f',[ NFORCE,3]) FIXED=fscanf(FP1,'%d',[ NVFIX,3])

ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0 for i=1:NELEM %生成弹性矩阵D

D= YOUNG/(1-POISS^2)*[1 POISS 0; POISS 1 0;

0 0 (1-POISS)/2] %计算当前单元的面积A

A=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2); 1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2); 1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)])/2 %生成应变矩阵B for j=0:2 b(j+1)=

COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);

c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1); end

B=[b(1) 0 b(2) 0 b(3) 0;... 0 c(1) 0 c(2) 0 c(3);...

c(1) b(1) c(2) b(2) c(3) b(3)]/(2*A); B1( :,:,i)=B; %求应力矩阵S=D*B S=D*B;

ESTIF=B'*S*THICK*A; a=LNODS(i,:); for j=1:3 for k=1:3

ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)…

=ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2

); end end end

%将约束信息加入总体刚度矩阵 for i=1:NVFIX if FIXED(i,2)==1

ASTIF(:,(FIXED(i,1)*2-1))=0; ASTIF((FIXED(i,1)*2-1),:)=0;

ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1; end

if FIXED(i,3)==1

ASTIF( :,FIXED(i,1)*2)=0; ASTIF(FIXED(i,1)*2,:)=0;

ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; end end

%生成荷载向量

ASLOD(1:2*NPION)=0; for i=1:NFORCE

ASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3) end

%求解内力

ASDISP=ASTIF\\ASLOD' ELEDISP(1:6)=0; for i=1:NELEM for j=1:3

ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2); end i

STRESS=D*B1(:, :, i)*ELEDISP' end

fclose(FP1);