有限元大作业matlab---课程设计例子

（1）主要变量：

ID： 问题类型码，ID＝1时为平面应力问题，ID=2时为平面应变问题 N_NODE： 节点个数 N_LOAD： 节点载荷个数

N_DOF： 自由度，N_DOF=N_NODE*2（平面问题） N_ELE： 单元个数 N_BAND： 矩阵半带宽

N_BC： 有约束的节点个数 PE： 弹性模量 PR： 泊松比 PT： 厚度

LJK_ELE(I,3)： 单元节点编号数组，LJK_ELE(I,1)，LJK_ELE(I,2)，LJK_ELE(I,3)分别放单元I的三个节点的整体编号

X(N_NODE), Y(N_NODE)：节点坐标数组，X(I),Y(I)分别存放节点I的x，y坐标值

P_LJK(N_BC,3)： 节点载荷数组，P_LJK(I,1)表示第I个作用有节点载荷的节点的编号，P_LJK(I,2)，P_LJK(I,3)分别为该节点沿x,y方向的节点载荷数值 AK(N_DOF,N_BAND)： 整体刚度矩阵 AKE(6,6)： 单元刚度矩阵

BB(3,6)： 位移??应变转换矩阵（三节点单元的几何矩阵） DD(3,3)： 弹性矩阵 SS(3,6)； 应力矩阵

RESULT_N(N_NOF)： 节点载荷数组，存放节点载荷向量，解方程后该矩阵存放节点位移

DISP_E(6):： 单元的节点位移向量 STS_ELE(N_ELE,3)： 单元的应力分量 STS_ND(N_NODE,3)： 节点的应力分量

（2）子程序说明：

READ_IN： 读入数据 BAND_K： 形成半带宽的整体刚度矩阵 FORM_KE： 计算单元刚度矩阵 FORM_P： 计算节点载荷 CAL_AREA：计算单元面积 DO_BC： 处理边界条件 CLA_DD： 计算单元弹性矩阵 SOLVE： 计算节点位移 CLA_BB： 计算单元位移??应变关系矩阵 CAL_STS：计算单元和节点应力

（3）文件管理： 源程序文件：

chengxu.for

程序需读入的数据文件：

BASIC.IN，NODE.IN，ELEMENT.IN（需要手工生成）

程序输出的数据文件：

DATA.OUT

（4）数据文件格式：

需读入的模型 基本信息文件BASIC.IN的格式如下表

实际需输入的数据 问题类型，单元个数，节第1行，每两个数之间用“，”号基本模型数据 点个数，有约束的节点隔开 数，有载何的节点数 材料性质 第2行，每两个数之间用“，”号弹性模量，泊松比，单元隔开 厚度 LJK_U(N_BC,3) 在材料性质输入行之后另起行，每位移约束的节点编号，该节点约束信息 两个数之间用“，”号隔开 节点x方向约束代码，该节点y方向代码， P_IJK(N_LOAD,3) 在节点约束信息输入行之后另起载荷作用的节点编号，该节点荷载信息 行，每两个数之间用“，”号隔开 节点x主向载荷，该节点y方向载荷，?? 栏目 格式说明 需读入的节点信息文件NODE.IN的格式如下表

栏目 节点信息 格式说明 实际需输入的数据 每行为一个节点的信息ND_ANSYS(N_NIDE) （每行三个数，每两个数节点号，该节点的x坐标，之间用空格或“，”分开） 该节点y方向坐标 需读入的单元信息文件ELEMENT.IN的格式如下表

栏目 格式说明 每行为一个单元的信息（每行有14个整型数，前4个为单元节点编号，对于3节点编号，第4个节点编号与第3个节点编号相同，后10个数无用，可输入“0”，每两 个整型数之间用至少一个空格分开） 实际需输入的数据 NE_ANSYS(N_ELE,14) 单元的节点号1（空格）单元的节点号2（空格）单元的节点号3（空格）单元的节点号4（空格） 0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0（空格）0 单元信息 输出结果文件DATA.OUT格式如下表

栏目 节点位移 实际输出的数据 I RESULT_N(2*I_ 1) RESULT_N(2*I) 节点号 x方向位移 y方向位移 IE STE_ELE(IE,1) STE_ELE(IE,2) STE_ELE(IE,3) 单元号 x方向应力 y方向应力 剪切应力 I STS-ND(I,1) STS-ND(I,2) STS-ND(I,3) 节点号 x方向应力 y方向应力 剪切应力 单元应力的三个分量 节点应力的三个分量 [算例原始数据和程序分析]：

（1）模型基本信息文件BASIC.IN的数据为

1,4,6,5,3 1.,0.,1.

1,1,0,2,1,0,4,1,1,5,0,1,6,0,1 1,-0.5,-1.5,3.,-1.,-1,6,-0.5,-0.5

（2）手工准备的节点信息文件NODE.IN的数据为

1 0.0 2.0 2 0.0 1.0 3 1.0 1.0 4 0. 0. 5 1.0 0. 6 2.0 0.

（3）手工准备的单元信息文件ELEMENT.IN的数据为

1 2 3 3 0 0 0 0 1 1 1 1 0 1 2 4 5 5 0 0 0 0 1 1 1 1 0 2 5 3 2 2 0 0 0 0 1 1 1 1 0 3 3 5 6 6 0 0 0 0 1 1 1 1 0 4

（4）源程序文件chengxu.for为：

PROGRAM FEM2D

DIMENSION IJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3), &RESULT_N(500),AK(500,100)

DIMENSION STS_ELE(500,3),STS_ND(500,3) OPEN(4,FILE='BASIC.IN') OPEN(5,FILE='NODE.IN') OPEN(6,FILE='ELEMENT.IN') OPEN(8,FILE='DATA.OUT')

OPEN(9,FILE='FOR_POST.DAT')

READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOAD IF(ID.EQ.1)WRITE(8,20) IF(ID.EQ.2)WRITE(8,25)

20 FORMAT(/5X,'=========PLANE STRESS PROBLEM========') 25 FORMAT(/5X,'=========PLANE STRAIN PROBLEM========') CALL READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT, & IJK_ELE,X,Y,IJK_U,P_IJK)

CALL BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE, & IJK_ELE,X,Y,PE,PR,PT,AK)

CALL FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N)

CALL DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N) CALL SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)

CALL CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, & STS_ELE,STS_ND) c to putout a data file

WRITE(9,70)REAL(N_NODE),REAL(N_ELE) 70 FORMAT(2f9.4)

WRITE(9,71)(X(I),Y(I),RESULT_N(2*I-1),RESULT_N(2*I), & STS_ND(I,1),STS_ND(I,2),STS_ND(I,3),I=1,N_NODE) 71 FORMAT(7F9.4)

WRITE(9,72)(REAL(IJK_ELE(I,1)),REAL(IJK_ELE(I,2)), &REAL(IJK_ELE(I,3)),REAL(IJK_ELE(I,3)),

&STS_ELE(I,1),STS_ELE(I,2),STS_ELE(I,3),I=1, N_ELE) 72 FORMAT(7f9.4) c

CLOSE(4) CLOSE(5) CLOSE(6) CLOSE(8) CLOSE(9) END c

c to get the original data in order to model the problem

SUBROUTINE READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR, &PT,IJK_ELE,X,Y,IJK_U,P_IJK)

DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3), & P_IJK(N_LOAD,3),NE_ANSYS(N_ELE,14) REAL ND_ANSYS(N_NODE,3) READ(4,*)PE,PR,PT