平面三角形单元常应变单元matlab程序的编制Word文件下载.docx
《平面三角形单元常应变单元matlab程序的编制Word文件下载.docx》由会员分享，可在线阅读，更多相关《平面三角形单元常应变单元matlab程序的编制Word文件下载.docx（15页珍藏版）》请在冰豆网上搜索。
●基本思想：
单元结点按右手法则顺序编号。
●荷载类型：
可计算结点荷载。
●说明：
主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。
%-----------------------------------------------------------------------------------------------------
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中读入结点个数，单元个数，结点荷载个数，受约束边界点数，弹性模量，泊松比，厚度，单元定义数组，结点坐标数组，结点力数组，约束信息数组；
程序中弹性模量仅输入了一个值，表明本程序仅能求解一种材料构成的结构，如：
钢筋混凝土结构、钢结构，不能求解钢筋混凝土－钢组合结构。
infinite
采用了命令fscanf，其中：
’%d’表示读入整数格式，’%f'
’表示读入浮点数；
1表示读取1个数，[A,B]形式表示读A行B列数组，[A,B]’表示将[A,B]转置，inf表示正无穷。
5调用子程生成单刚，组成总刚并加入约束信息
ASSEMBLE（）
6调用子程生成荷载向量
FORMLOAD（）
7计算结点位移向量
ASDISP=ASTIF\ASLOD'
8调用子程计算单元应力
WRITESTRESS（）
9关闭输出数据文件
fclose（FP2）;
%*******************************************************************读取ASSEMBLE子程%*******************************************************************
functionASSEMBLE（）
%所引用的全局变量：
globalNPIONNELEMNVFIXLNODSASTIFTHICK
globalBMATXSMATXAREAFIXED
%----------------------------------------------------------------------------------------------------
%计算单刚并生成总刚
ASTIF（1:
2*NPION,1:
2*NPION）=0;
%张成特定大小总刚矩阵并置0
Arraystiffness
建立单元刚度矩阵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;
%计数器（当前已分析的节点数）
NVFIX：
受约束边界点数
i=0;
%计数器（当前已处理的约束数）
tmp（NVFIX）=0;
%临时存被约束的序号
whilei<
NVFIX
forj=-1:
ifFIXED（NUM,j+3）==1%若发现约束
temporary
i=i+1;
%计数器自增
tmp（i）=FIXED（NUM）*2+j;
%求约束序号
end
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两个方向的结点力。
fori=1:
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;
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。
%----------------------
升级会员 