三角形常应变单元程序的编制与使用共18页文档Word格式文档下载.docx
- 文档编号:17012026
- 上传时间:2022-11-27
- 格式:DOCX
- 页数:14
- 大小:63.35KB
三角形常应变单元程序的编制与使用共18页文档Word格式文档下载.docx
《三角形常应变单元程序的编制与使用共18页文档Word格式文档下载.docx》由会员分享,可在线阅读,更多相关《三角形常应变单元程序的编制与使用共18页文档Word格式文档下载.docx(14页珍藏版)》请在冰豆网上搜索。
杨士勋(唐初学者,四门博士)《春秋谷梁传疏》曰:
“师者教人以不及,故谓师为师资也”。
这儿的“师资”,其实就是先秦而后历代对教师的别称之一。
《韩非子》也有云:
“今有不才之子……师长教之弗为变”其“师长”当然也指教师。
这儿的“师资”和“师长”可称为“教师”概念的雏形,但仍说不上是名副其实的“教师”,因为“教师”必须要有明确的传授知识的对象和本身明确的职责。
对于平面问题,结构离散常用的网格形状有三角形、矩形、任意四边形,以三个顶点为节点的三角形单元是最简单的平面单元,它较矩形或四边形对曲边边界有更好的适应性,而矩形或四边形单元较三节点三角形有更高的计算精度。
Matlab语言是进行矩阵运算的强大工具,因此,用Matlab语言编写有限元中平面问题的程序有优越性。
本章将详细介绍如何利用Matlab语言编制三角形常应变单元的计算程序,程序流程图见图1。
有限元法中三节点三角形分析结构的步骤如下:
1)整理原始数据,如材料性质、荷载条件、约束条件等,离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2)单元分析,建立单元刚度矩阵。
3)整体分析,建立总刚矩阵。
4)建立整体结构的等效节点荷载和总荷载矩阵
5)边界条件处理。
6)解方程,求出节点位移。
7)求出各单元的单元应力。
8)计算结果整理。
计算结果整理包括位移和应力两个方面;
位移计算结果一般不需要特别的处理,利用计算出的节点位移分量,就可画出结构任意方向的位移云图;
而应力解的
图1程序流程图
误差表现在单元内部不满足平衡方程,单元
与单元边界处应力一般不连续,在边界上应
力解一般与力的边界条件不相符合。
1.1程序说明
%三角形常应变单元求解结构主程序
●功能:
运用有限元法中三角形常应变单元解平面问题的计算主程序。
●基本思想:
单元结点按右手法则顺序编号。
●荷载类型:
可计算结点荷载。
●说明:
主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程,即实现程序流程图的程序表达。
1程序准备
formatshorte%设定输出类型
clearall%清除所有已定义变量
clc%清屏
formatshorte-设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法;
clearall-清除所有已定义变量,目的是在本程序的运行过程中,不会发生变量名相同等可能使计算出错的情况;
clc-清屏,使屏幕在本程序运行开始时
2全局变量定义
globalNNODENPIONNELEMNVFIXNFORCECOORDLNODSYOUNGPOISSTHICK
globalFORCEFIXED
globalBMATXDMATXSMATXAREA
globalASTIFASLODASDISP
globalFP1
NNODE—单元结点数,NPION—总结点数,NELEM—单元数,NVFIX—受约束边界点数,NFORCE—结点力数,COORD—结构结点坐标数组,LNODS—单元定义数组,YOUNG—弹性模量,POISS—泊松比,THICK—厚度
FORCE—节点力数组(n,3)n:
受力节点数目,(n,1):
作用点,(n,2):
x方向,(n,3):
y方向;
FIXED—约束信息数组(n,3)n:
受约束节点数目,(n,1):
约束点(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0
BMATX—单元应变矩阵(3*6),DMATX—单元弹性矩阵(3*3),SMATX—单元应力矩阵(3*6),AREA—单元面积
ASTIF—总体刚度矩阵,ASLOD—总体荷载向量,ASDISP—结点位移向量
FP1—数据文件指针
3打开文件
FP1=fopen('
input.txt'
'
rt'
);
%打开输入数据文件存放初始数据
FP1=fopen('
-打开已存在的输入数据文件input.txt,且设置其为只读格式,使程序在执行过程中不能改变输入文件中的数值,并用文件句柄FP1来执行
FP2=fopen('
output.txt'
wt'
-打开输出数据文件,该文件不存在时,通过此命令创建新文件,该文件存在时则将原有内容全部删除。
该文件设置为可写格式,可在程序执行过程中向输出文件写入数据。
4读入程序控制信息
NPION=fscanf(FP1,'
%d'
1)%结点个数(结点编码总数)
NELEM=fscanf(FP1,'
1)%单元个数(单元编码总数)
NFORCE=fscanf(FP1,'
1)%结点荷载个数
NVFIX=fscanf(FP1,'
1)%受约束边界点数
YOUNG=fscanf(FP1,'
%e'
1)%弹性模量
POISS=fscanf(FP1,'
%f'
1)%泊松比
THICK=fscanf(FP1,'
1)%厚度
LNODS=fscanf(FP1,'
[3,NELEM])'
%单元定义数组(单元结点号)
建立LNODS矩阵,该矩阵指出了每一单元的连接信息。
矩阵的每一行针对每一单元,共计NELEM;
每一列相应为单元结点号(编码)、按逆时针顺序输入。
命令中,[3,NELEM]’表示读取NELEM行3列数据赋值给LNODS矩阵。
显然,LNODS(i,1:
3)依次表示i单元的i,j,k结点号。
COORD=fscanf(FP1,'
[2,NPION])'
%结点坐标数组
建立COORD矩阵,该矩阵用来存储各结点x,y方向的坐标值。
从FP1文件中读取全部结点个数NPOIN的坐标值,从1开始按顺序读取。
COORD(i,1:
2)表示第i个结点的x,y坐标。
FORCE=fscanf(FP1,'
[3,NFORCE])'
%结点力数组
(n,3)n:
受力结点数目,(n,1):
y方向
FIXED=fscanf(FP1,'
[3,inf])'
%约束信息数组
约束点(n,2)与(n,3)分别为约束点x方向和y
方向的约束情况,受约束为1否则为0
●总体说明:
从输入文件FP1中读入结点个数,单元个数,结点荷载个数,受约束边界点数,弹性模量,泊松比,厚度,单元定义数组,结点坐标数组,结点力数组,约束信息数组;
程序中弹性模量仅输入了一个值,表明本程序仅能求解一种材料构成的结构,如:
钢筋混凝土结构、钢结构,不能求解钢筋混凝土-钢组合结构。
采用了命令fscanf,其中:
’%d’表示读入整数格式,’%f'
’表示读入浮点数;
1表示读取1个数,[A,B]形式表示读A行B列数组,[A,B]’表示将[A,B]转置,inf表示正无穷。
5调用子程生成单刚,组成总刚并加入约束信息
functionASSEMBLE()
6调用子程生成荷载向量
functionFORMLOAD()
7计算结点位移向量
ASDISP=ASTIF\ASLOD'
8调用子程计算单元应力
functionWRITESTRESS()
9关闭输出数据文件
fclose(FP2);
%*******************************************************************读取ASSEMBLE子程%*******************************************************************
functionASSEMBLE()
%所引用的全局变量:
globalNPIONNELEMNVFIXLNODSASTIFTHICK
globalBMATXSMATXAREAFIXED
%计算单刚并生成总刚
ASTIF(1:
2*NPION,1:
2*NPION)=0;
%张成特定大小总刚矩阵并置0
建立单元刚度矩阵ASTIF,该矩阵的行列数均为2*NPION,NPION表示结点数,每个结点有两个方向的力和位移。
fori=1:
NELEM
FORMSMATX(i)%%%调用应力子程序
ESTIF=BMATX'
*SMATX*THICK*AREA;
%求解单元刚度矩阵
a=LNODS(i,:
%临时向量,用来记录当前单元的节点编号
forj=1:
3
fork=1:
ASTIF((a(j)*2-1):
a(j)*2,(a(k)*2-1):
a(k)*2)=ASTIF((a(j)*2-1):
a(k)*2)+ESTIF(j*2-1:
j*2,k*2-1:
k*2);
%跟据节点编号对应关系将单元刚度分块叠加到总刚矩阵中
end
FORMSMATX(i)调用应力子程序,提取i单元的应力矩阵SMATX;
a=LNODS(i,:
)记录i单元的三个结点编号;
for…end循环语句表示行从1到3循环,列从1到3循环,将单刚中的元素叠加至总刚中:
a(k)*2)表示总刚中第a(j)*2-1到:
a(j)*2行,第a(k)*2-1到a(k)*2列的元素由单刚中第j*2-1到j*2行,第k*2-1到k*2列的元素叠加而得,a(j)*2即将单元中的位移编码对应到整体的位移编码。
%将约束信息加入总刚(置0置1法)
NUM=1;
%计数器(当前已分析的节点数)
i=0;
%计数器(当前已处理的约束数)
tmp(NVFIX)=0;
%临时存被约束的序号
whilei<
NVFIX
forj=-1:
ifFIXED(NUM,j+3)==1%若发现约束
i=i+1;
%计数器自增
tmp(i)=FIXED(NUM)*2+j;
%求约束序号
NUM=NUM+1;
tmp(NVFIX)=0,形成一个元素值均为0的一行NVFIX列的行向量,
执行while语句,首先判断i是否小于控制数据NVFIX,若小于则往下进行,若不小于则退出。
执行for语句,FIXED(NUM,j+3)表示约束信息数组中第NUM行第j+3列的元素,j从-1到0,即j+3表示2到3列,即约束信息数组中描述结点x和y方向受约束的情况,判断FIXED(NUM,j+3)若等于1,则约束数自增,若不等于1,跳出。
FIXED(NUM)表示FIXED(NUM,1),tmp(i)=FIXED(NUM)*2+j计算整体约束序号,将序号放入tmp行向量中的i列。
fori=1:
ASTIF(1:
2*NPION,tmp(i))=0;
%将一约束序号处的总刚列向量清0
ASTIF(tmp(i),1:
%.将一约束序号处的总刚行向量清0
ASTIF(tmp(i),tmp(i))=1;
%将行列交叉处的元素置为1
后处理法中置0置1法,设
(包括
),则将总刚中的主元素
Kjj换为1,j行和j列的其他元素均改为零。
%读取FORMSMATX子程
functionFORMSMATX(ELEMENT)%计算应力矩阵
%引用所需的全局变量
globalNPIONNELEMCOORDLNODSYOUNGPOISS
globalBMATXDMATXSMATXAREA
%生成弹性矩阵D
a=YOUNG/(1-POISS^2);
DMATX(1,1)=1*a;
DMATX(1,2)=POISS*a;
DMATX(2,1)=POISS*a;
DMATX(2,2)=1*a;
DMATX(3,3)=(1-POISS)*a/2;
平面应力问题的弹性矩阵
其中,E为弹性模量,
为泊松比。
i=ELEMENT;
%i为当前所计算的单元号
%计算当前单元的面积
AREA=det([1COORD(LNODS(i,1),1)COORD(LNODS(i,1),2);
...
1COORD(LNODS(i,2),1)COORD(LNODS(i,2),2);
1COORD(LNODS(i,3),1)COORD(LNODS(i,3),2);
])/2;
det表示求矩阵行列式的值,
,其中
分别表示一个三角形单元的三个节点坐标。
MATLAB中若一行中无法写下一个完整的命令,则可以在行尾加入3个连续的英文句号,表示命令余下的部分在下一行出现。
%-----------------------------------------------------------------------------------------------------%生成应变矩阵B
forj=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);
BMATX=[b
(1)0b
(2)0b(3)0;
0c
(1)0c
(2)0c(3);
c
(1)b
(1)c
(2)b
(2)c(3)b(3)]/(2*AREA);
应变矩阵
rem表示求余函数,rem(x,y)命令返回的是x-n.*y,当y
0时,n=fix(x./y),其中fix为最近的整数取整。
SMATX=DMATX*BMATX;
%求应力矩阵S=D*B
%读取FORMLOAD子程
functionFORMLOAD()%本子程生成荷载向量
%所需引用的全局变量
globalASLODNPIONNFORCEFORCE
ASLOD(1:
%张成特定大小的0向量
ASLOD为总体荷载向量,每个结点有x,y两个方向的结点力。
NFORCE
ASLOD((FORCE(i,1)*2-1):
FORCE(i,1)*2)=FORCE(i,2:
3);
end
FORCE(i,1)为作用点,FORCE(i,2:
3)为x,y方向的结点力。
%将有约束处的荷载置为0
NUM=1;
i=0;
tmp(NVFIX)=0;
whilei<
ifFIXED(NUM,j+3)==1
i=i+1;
tmp(i)=FIXED(NUM)*2+j;
end
NUM=NUM+1;
end
ASLOD(tmp(i))=0;
置0置1法,同上。
ASDISP=ASTIF\ASLOD'
%计算结点位移向量
%读取WRITESTRESS子程
functionWRITESTRESS()
%求解内力,即两个方向的正应力与一个剪应力(
)
%所引用的全局变量
globalNELEMNPIONSMATXASDISPLNODS
ELEDISP(1:
6)=0;
%当前单元的结点位移向量
ELEDIS为单元的结点位移
ELEDISP(j*2-1:
j*2)=ASDISP(LNODS(i,j)*2-1:
LNODS(i,j)*2);
FORMSMATX(i)%%%调用子程求应力矩阵
i
STRESS=SMATX*ELEDISP'
%求内力
●说明:
通过循环,依次从结构位移列阵中对号,赋值给第i个单元的单元位移向量ELEDISP。
1.2程序应用举例
【例1】利用三角形三节点常应变单元程序计算图2所示的结构。
设弹性模量
,泊松比为0,厚度为1。
%-----------------------------------------------------------------
输入数据文件input.txt为:
64161.0e00.01
312
524
325
635图2
0.02.0
0.01.0
1.01.0
0.00.0
1.00.0
2.00.0
10-1
110
210
411
501
601
说明:
第一行:
读入程序控制信息
第二、三、四、五行:
读入单元连接信息:
;
%单元定义数组,单元结点号,逆时针输入
第六、七、八、九、十、十一行:
读入结点坐标
[2,NPOIN])'
%结点坐标值,第1列为x坐标值,第2列为y坐标值
第十一行:
读入结点荷载信息
%结点号,X向结点荷载数值,Y向结点荷载数值(与坐标轴方向一致为正)
第十二、十六行:
读入零位移信息
%结点号,X向约束,Y向约束
1.3程序的改进要点
上述三角形三节点程序反映了有限元的基本思路,可以计算简单的平面应力或平面应变问题。
在熟练掌握了程序的编制与使用后,可在以下几方面对程序进行改进,以加深对矩阵位移法及MATLAB语言编程的理解:
1、本程序的弹性模量仅能输入一个数值,意味着程序仅能计算由同种材料构成的结构。
考虑如何改进使程序可计算由不同材料构成的组合结构。
2、本程序仅能计算结点集中荷载类型,考虑如何编制体积力、表面力、跨中集中力等类型的程序。
3、考虑如何编制有支座位移的程序。
4、本程序最后的结果没有生成输出文件,可以编制输出文件。
5、本程序计算的应力没有进行结果处理,可以编制最后结果处理的程序。
6、可以在此程序基础上编制三角形六节点、四边形四节点等程序。
综上所述,本章的三角形三节点常应变程序体现了如何将有限元法的计算方法和过程用MATLAB程序语言表达出来,重点放在程序架构和流程的建立以及算法实现方面,主要依赖手工操作-手工输入初始数据(前处理)、手工绘制计算结果(后处理),目的是使学生能够清晰、明确地掌握有限元法的基本理论和概念,熟练掌握应用MATLAB语言编制、修改和调试简单程序的能力。
希望以上资料对你有所帮助,附励志名言3条:
1、有志者自有千计万计,无志者只感千难万难。
2、实现自己既定的目标,必须能耐得住寂寞单干。
3、世界会向那些有目标和远见的人让路。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 三角形 应变 单元 程序 编制 使用 18 文档