有限元程序设计报告报告

有限元程序设计报告

课程名称:有限元程序设计

指导教师：张*

学校：重庆大学

专业：工程力学01班

******

学号：********

2015年7月8日

有限元程序设计报告

一、前言

有限元方法（the Finite Element Method）是起源与上个世纪50、60 年代，基于弹性力学变分原理的一种近似计算方法，也是当今工程分析中获得最广泛应用的数值计算方法。由于它的通用性和有效性，受到工程技术界的高度重视。伴随着计算机科学和技术的快速发展，现已成为计算机辅助设计（CAD）和计算机辅助制造（CAM）的重要组成部分。有限元程序系统通常包括前处理、有限元程序本体和后处理三部分。前处理包括几何实体模型的建立、材料参数的赋值、位移边界条件的定义、载荷的定义、分析问题类型的定义、单元类型的选择和网格的划分等。（分析问题类型如静力分析、动力特性分析、动力响应、温度场分析、电磁场分析、流体动力学分析等）有限元程序本体是有限元程序系统的核心部分，其功能是实现各种问题的计算。后处理则是将计算结果用图形、曲线和表格的形式表达。（通常包括结构的变形图、应力、应变分布云图等）

本课程设计则是针对有限元程序本体，参照教学程序（FEATP），编写简单的有限元程序以计算简单的平面应力、平面应变和轴对称问题，并将其结果与有限元商用软件（ANSYS）的计算结果，以及问题的理论值进行比较，从而验证程序以及问题模型建立的正确性。

1．设计目的

1）通过编写简单的有限元程序熟悉用有限元方法解决实际问题的基本步骤和过程，体会这种方法的处理手段。

2）在V isual Fortran 中编写程序，熟悉并巩固F ortran 语言的语法、算法，学习程序的调试方法，并体会其在执行某个具体算例时，文件的输入、输出以及程序的执行过程。

2．设计内容

1）以教学程序（FEATP）为参照，编写程序，计算简单的平面应力（Plane Stress），平面应变(Plane Strain)问题，验证程序的正确性。

2）在具体的算例中，对同一问题，在程序和ANSYS 中采用不同的单元和网格划分方式，将其结果与理论值进行对比，体会不同的单元和网格划分对问题解的影响，从而判断模型的正确性和合理性。

3）总结在编写程序和算例中遇到的问题和解决方法，写出自己的心得体会。

二、弹性力学平面问题有限元方法的基本公式

⎢ ⎢ ⎥ i

m 平面问题

1．三角形单元

(1) 位移

⎡ u i ⎤ ⎢ ⎥ ⎢ v i ⎥ ⎡ ⎤ ⎡

⎤⎢ ⎥ u

N i u = ⎢ ⎥ = ⎢

0 N j 0 N m 0 u j

⎥⎢ ⎥ ⎣v ⎦ ⎣ 0

N i

N j

N m ⎦⎢ v j ⎥

⎢ ⎥

⎡ a i ⎤ u m ⎥ ⎢⎣v m ⎥⎦

= [

IN IN j IN m ]

⎢ a j ⎥

= [

N N ⎢⎣a m ⎥⎦

N ]

a e = Na e

其中，形函数

a i ,

b i ,

c i 是取决于节点坐标的常数

应变矩阵 B 的分块子矩阵是

⎥

⎥⎥

⎦⎤⎢⎢

⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎥

⎥⎥⎥⎥⎥⎦

⎤

⎢⎢⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂==i i i i i i i i b c c b A N N x y y x LN B 00210000

3 节点单元的应变矩阵是

1

⎥

(2) 形成单元的刚度矩阵和等效节点载荷列阵

i j

()

()m j i y i c x i b A N ,,i

a 21

i ++=[]

⎥⎥⎥⎦

⎤

⎢⎢⎢

⎣⎡==m m j j i i m j i m j i m j i b c b c b c c c c b b b A B B B B 00000021

苏世宏《有限元课程设计》

⎰

B T ⎢ [

i ⎥

i

m m j m i j m

K e = B T DBdV = B T DBtA

V e

P e = P e + P e + P e + P e

f

s

σ 0 ε0

(3) 集成结构的刚度矩阵和等效节点载荷列阵

K = ∑ K e = ∑⎰V DBdV

e

P = P f e

+ P s P σ e

+ P 0

+ P F ∑

f S σ 0 ε 0

F

= (P e + P e

P e

e

+ P e

) + P

(4) 引入强制边界条件（消除 K 的奇异性） (5) 求解有限元解方程，得到节点位移 a

Ka = P

(6) 计算单元应变和应力

⎡ ε x ⎤

ε = ⎢ ε ⎢γ ⎥

e ⎥ = Lu = LNa = L N i N j

⎥ N ]

a e

⎣ = [

B xy ⎦

B j B ]

a e = Ba e

⎡σ x ⎤

⎢ ⎢σ y ⎢τ ⎥ = D ε = DBa e = Sa e

⎥ ⎣ xy ⎦

应力矩阵 S = DB = D [

B B B ]= [

S S S ]

S 的分块子矩阵为

⎡ ⎢ b i

⎤

ν 0 c i ⎥

S i = DB i =

E 0

⎢

2

⎢ ν 0b i

c ⎥

（i, j, m ）

i

⎥

2(1 -ν 0 ) A ⎢1 -ν 0 1 -ν 0 ⎥

⎣⎢ 2 c i 2

b i

⎥⎦

对于平面应力问题

σ = 0 ε y