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

有 限 元 大 作 业 程 序 设 计

学校：天津大学

院系：建筑工程与力学学院 专业：01级工程力学 姓名：刘秀 学号：\\\\\\\\\\\ 指导老师：

连续体平面问题的有限元程序分析

[题目]：

如图所示的正方形薄板四周受均匀载荷的作用，该结构在边界

上受正向分布压力，

m kN

p 1=，同时在沿对角线y 轴上受一对集中压

力，载荷为2KN ，若取板厚1=t ，泊松比0=v 。

[分析过程]：

由于连续平板的对称性，

只需要取其在第一象限的四分之一部分

参加分析，然后人为作出一些辅助线将平板“分割”成若干部分，再为每个部分选择分析单元。采用将此模型化分为4个全等的直角三角型单元。利用其对称性，四分之一部分的边界约束，载荷可等效如图所示。

[

用和单元信息文件DATA.OUT。

位移模式：用用线性位移模式

载荷类型：节点载荷，非节点载荷应先换算为等效节点载荷

材料性质：弹性体由单一的均匀材料组成

约束方式：为“0”位移固定约束，为保证无刚体位移，弹性体至少应有对三个自由度的独立约束

方程求解：针对半带宽刚度方程的Gauss消元法

输入文件：由手工生成节点信息文件NODE.IN，和单元信息文件ELEMENT.IN

结果文件：输出一般的结果文件DATA.OUT

程序的原理如框图：

（1）

ID ：

ID=2时为平面应变问题

（平面问题）

，LJK_ELE(I,1)，LJK_ELE(I,2)，

X(I),Y(I)分别存放节点I 的x ，y

表示第I 个作用有节点载荷的节点x,y 方向的节点载荷数值

存放节点载荷向量，解方程后该矩

（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的格式如下表

需读入的节点信息文件NODE.IN的格式如下表

需读入的单元信息文件ELEMENT.IN的格式如下表

输出结果文件DATA.OUT格式如下表

[算例原始数据和程序分析]：

（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)

D IMENSION 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.DA T')

READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOAD

IF(,20)

IF(,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 FORMA T(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)

E ND

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

READ(4,*)((IJK_U(I,J),J=1,3),I=1,N_BC)

READ(4,*)((P_IJK(I,J),J=1,3),I=1,N_LOAD)

READ(5,*)((ND_ANSYS(I,J),J=1,3),I=1,N_NODE)

READ(6,*)((NE_ANSYS(I,J),J=1,14),I=1,N_ELE)

DO 10 I=1,N_NODE

X(I)=ND_ANSYS(I,2)