重庆大学本科学生课程设计任务书
课程设计题目有限元程序设计
学院航空航天学院专业工程力学年级2011级
已知参数和设计要求：
1．独立完成有限元程序设计。
2．独立选择计算算例，并能通过算例判断程序的正确性。
3．独立完成程序设计报告，报告内容包括理论公式、程序框图、程序本体、计算算例，算例结果分析、结论等。
学生应完成的工作：
1.复习掌握有限单元法的基本原理。
2.掌握弹性力学平面问题3节点三角形单元或4节点等参单元有限元方法的计算流程，
以及单元刚度矩阵、等效节点载荷、节点应变、节点应力和高斯积分等的计算公式。
3.用Fortran语言编写弹性力学平面问题3节点三角形单元或4节点等参单元的有限元
程序。
4.在Visual Fortran 程序集成开发环境中完成有限元程序的编辑和调试工作。
5.利用编写的有限元程序，计算算例，分析计算结果。
6.撰写课程设计报告。
目前资料收集情况（含指定参考资料）：
1．王勖成，有限单元法，北京：高等教育出版社，2002。
2．O.C. Zienkiewicz, R. L. Taylor, Finite Element Method, 5th Edition, McGraw-Hall Book Company Limited, 2000。
3．张汝清，董明，结构计算程序设计，重庆：重庆大学出版社，1988。
课程设计的工作计划：
1.第1周星期一上午：教师讲解程序设计方法，程序设计要求和任务安排。
2.第1周星期一至星期二完成程序框图设计。
3.第1周星期三至第2周星期四完成程序设计。
4.第2周星期五完成课程设计报告。
任务下达日期 2014 年 6 月 6 日完成日期 2014 年 6 月30 日
指导教师（签名）学生（签名）
一、前言。
随着科技的发展，理论分析越来越不能满足各行各业的需求，而大量的实
验又耗资巨大，所以依靠计算机的分析方法越来越流行。基于有限单元法的abaqus 等软件就是通过数值方法来解决工程问题的一种高效的方法。
有限元法的基本思想是“化整为零”，把一个复杂的结构分解成相对简单的“单元”，各单元之间通过结点相互连接。单元内的物理量由单元结点上的物理量按一定的假设内插得到，这样就把一个复杂结构从无限多个自由度简化为有限个单元组成的结构。我们只要分析每个单元的力学特性，然后按照有限元法的规则把这些单元“拼装”成整体，就能够得到整体结构的力学特性。
本文要探讨的事有限单元法的核心，利用fortran 编程来编写二维有限单元法的计算过程，并通过abaqus 来验证其准确性，然后在用abaqus 详细分析一些二维的问题，从而较深对有限单元法的理解。
二、理论依据。
虚位移原理和虚应原理分别是平衡方程与力的边界条件和几何方程与位移边界条件的等效积分形式。将物理方程引入虚位移原理和虚应力原理可以分别导出最小位能原理和最小余能原理，它们本质上和等效积分的伽辽金“弱”形式相一致，这是建立弹性力学有限原方程一般表达格式的理论基础。
弹性力学平面问题有限元分析步骤为：
单元定义 单元离散化 单元推导 总装 联立方程组求解
平面4、8 节点有限元公式及计算原理 1、平面问题基本公式。
(1)、有限元平衡方程：
P Ka = (2)、单元刚度矩阵：
ηξ?∑Ω===
e n
i T
i T
e
d d J DBt B H DBtdxdy B K 1||
(3)、单元等效结点载荷：
∑∑??==Ω+=+
=
n
i T n
i T
i S T
T
e
d d J Tt N d d J ft N
H Ttds N ftdxdy N
P e e
1
1
||||η
ξηξσ
采用四边形单元 离散后对单元进行编号，并且给定单元各节点的整体编码及局部坐标系下按逆时针局部编码 单元刚度矩阵，利用等参元性质和雅阁比矩阵进行组集，建立整体刚度矩阵
利用三角分解等方法求解
2.平面问题计算步骤。
（1）通过Serendipity 四边形单元格式构造插值函数。对于4节点单元，插值函数为：
)1)(1(4
1
ηηξξi i i N ++= （i=1,2,3,4）
对于8节点单元，插值函数为： 65112121?N N N N --= 65222121?N N N N --= 7
6332121?N N N N --= 874421?N N N N --= )1)(1(2125ηξ--=N )1)(1(2126
ξη+-=N )1)(1(2127ηξ+-=N )1)(1(2
1
28ξη--=N
)1)(1(4
1?ηηξξi
i i N ++= （i=1,2,3,4） （2）通过?????
?
????
???
?-+=210
00101100
200v v v v E
D ， 平面应力问题：
E E =0 v v =0 平面应变问题：201v E E -= v
v v -=10
式中E 和v 分别为材料Young 氏模量和Poisson 比。
???
??
?
??
?
?????????????????=????????????
??
??????????????????==x N y N y
N x N N N x y y x LN B i
i i i i i
i i 00000
0 [
]
m j i
B B B B =
计算出D 和B 矩阵。
（3）单元矩阵的变换：
?????
?
??????????????=????????????????=??
??????????????????????????????=????????????????-ηη
ξξηξηξy x y x J N N y N x N y N x N J N N i i i i i i i i ，其中：，于是，1J 所以面积微元可表示成ηξηξd d J d d dA ||=?=。
（4）运用高斯积分公式求积分公式近似数值积分解： ∑?==n
i i i b
a
F H d F 1
)()(ξξξ
其中权函数ξξd l H b
a
n i i ?-=)()1(
（5）单元刚度矩阵与载荷列向量可表示为：
η
ξ?∑Ω===e
n
i T i T e
d d J DBt B H DBtdxdy B K 1
||，
∑∑??==Ω+=+=n
i T n
i T
i S T
T
e
d d J Tt N d d J ft N H Ttds N ftdxdy N P e
e
1
1||||η
ξηξσ
将其代入整体刚度矩阵与载荷列向量可得整体刚度矩阵K 和载荷列向量P 。
（6）总势能为：p a Ka a T T p -=2
1
π。
根据最小势能原理
0=??a
p π可得：
P Ka =
运用上式即可求得求得各节点位移。
3.轴对称问题。
对于轴对称问题唯一不同于平面问题的是轴对称问题多了一个旋转分量，从而使问题变得更加复杂。下面简单介绍一下轴对称问题的区别于平面问题的特点。
轴对称应变：[
]
e m j
i
a B B r u r w z u z w r u B εγεεεθrz z r =?????
?
???
???????????????+??????=??????????????=
??????????????
++=
????????
?
??????
?????????????=00021000r z c b r a b c c b A r
N z N r N z N r N B i
i i i i i i i i i i i
i
而平面问题没有旋转项：θε
轴对称问题的应力也不同：0000θrz z r σD ε-εσ)ε-ε(στσσσ+=+=????
??????????=D D
单元刚度阵：???===e
e
e
ΩT ΩT ΩT drdz B π2θdrdz B B DBr DBrd DBdv K e
等效结点载荷：??==e
e
ΩT ΩT N π2N frdrdz fdv P e f
??==e σ
e
σ
π2s T s T e s Trds N Tds N P
??-=-=e
e
v T v T e rdrdz B dv B P 00σσπ2σ
三、有限元分析程序。
具体程序见附录一。
本程序基于Fortran 语言编写，可对二维问题行有限元分析计算。计算采用的单元形式为四边形四节点单元或者四边形八节点单元，对于对称和非对称矩阵采用变带宽存储方法，大量节约内存空间，输入为单元坐标，单元定义等，最终输出结果为单元各节点的位移。
程序简要分析： 1.有限元程序系统的组成及分析过程：
我们本次设计主要针对二维弹性力学静力学问题。
有限元程序系统通常包括处理、有限元程序本体和后处理三部分。
前处理 CAD 建模（几何模型），材料参数，边界条件，分析类型定义，网格划分等。
有限元分析本体程序
核心：各种计算方法的实现
后处理
以图形、曲线和表格等方
式表达、分析计算结果
图 3.1 有限元软件系统的组成及分
2.计算流程：
对于线弹性力学问题，其有限元的计算程序流程如下图：
3.有限元程序主程序流程。
打开输入数据，调用PMESH ，调用PROFIL 计算K 的对角元地址和方程号，分配数组大小，调用ASSEM 组集整体K ，调用PLOAD 组集载荷列向量，调用DATRI 对刚度矩阵进行三角分解，调用DASOL 解方程，调用PRTDIS 输出结点位移，调用TRESSS 计算GAUSS 点和结点应变和应力。 4.变量说明
NUMNP ：节点总数 NUMEL ：单元总数 NUMMAT ：材料种类 NLOAD ：载荷数
输入离散模型数据型数据
计算单元刚度阵
组集结构刚度矩阵
计算单元等效结点载荷
组集结构结点载荷列阵 引入位移边界条件
求解线性方程组 其它辅助计算 输出结果
结 束
单元循环
形成 K
形成 P
消除 K 的奇异
求解 Ka =P 得结点位移 a
计算应力、应变等
图3.2 有限元计算流程
NDM ：每个节点坐标方向数 NDF ：节点自由度 NEN ：每个单元节点数。
ID （K ，J ）：节点J 的第K 个自由度约束编号，=0：无约束；=1：固定约束。 X （K ，J ）：节点J 的第K 个坐标
IX （K ，J ）：J 单元中K 节点的全局节点编号 F （K ，J ）：结点J 第K 方向上的集中载荷
F （NDF ，NLOAD ）：集中载荷；X （NDM ，NUMNP ）：结点坐标 IX （NEN ，NUMEL ）：单元定义；EE ：杨氏弹性模量；XNU ：泊松比 F(K,J)：节点J 在K 方向上集中载荷大小
三、程序验证算例。
通过与abaqus 计算结果比较，验证其准确性。 算例：悬臂梁受均布压缩。
如图所示，左端固定端约束，上面受均布载荷，把它等效到每个节点上。 悬臂梁的尺寸是200*50mm 。均布载荷mm N 2.材料参数：E=200GPa,ν=0.3。
网格划分如下图:
图4.1 悬臂梁的网格点与受力情况
总共划分了40个结点，27个单元，分别采用四、八节点单元进行计算。
1.采用Fortran 程序运行。
INPUT.DAT 文件如下：（四节点）
图4.2 四节点结点及单元定义
说明：由于后面结点数量大，为了避免繁琐的输入结点坐标，单元定义等，选择从abaqus里面导入以上信息，但是再导入过程中会存在格式不匹配，于是把Fortran程序里面的格式化读入部分进行了稍微的修改，以满足可以读入abaqus的数据。
OUTPUT.DAT输出为：(限于篇幅，这里只列出结点的位移)
采用八结点结点及单元定义：
由于输入文件INPUT.DAT 与四节点输入相似，在这里不再累赘列出，只列出 其结点位移。
八结点单元结点位移输出OUTPUT.DAT 。（摘取部分）
图4.3 八结点结点及单元定义
2.abaqus 计算。
Abaqus 同Fortran 计算划分的网格相同，结点及单元定义也相同，分别计算四节点，八结点的input 输出，对比来验证程序的可行性。
首先采用四节点，查看其结点位移：
图4.4 Fortran 八结点结点位移输出
图4.5 abaqus 四节点结点位移输出
采用八结点单元，其结点位移：（部分，为了方便比较，列出相对应的结点位移）
3.理论计算。
应用挠曲线近似微分方程求出其挠度方程：
)24/()64(222EI l lx x qx w +-=
计算得出在40结点处的y 方向位移为：-0.066667mm
与以上结果比较可以看出，百分误差为
()%64%100186962.0-/066667.0--186962.0-=?=）
（σ 分析原因：
照成这样大的误差是因为悬臂梁的长度和宽度之比较大，不能简单地使用上面的近似计算公式。其次，40结点处没有在中性轴上，也会导致误差。由于划分的网格较少，所以计算出来的位移也会有误差。
4.分析比较。
（1）分析采用四节点单元计算结点位移。 随机选取三个结点进行比较，结点4、20和38. Fortran 计算结果：
图4.6 abaqus 八结点单元结点位移
Fortran位移Abaqus位移相对误差4结点X -0.008699 -0.0086987 0.0034%
Y -0.006909 -0.0069092 0.0029% 20结点X -0.009580 -0.0095804 0.0042%
Y -0.1867 -0.18672 0.011% 38结点X 0.02964 0.029644 0.013%
Y -0.1340 -0.13405 0.038%
可以看出，Fortran的计算结果与abaqus计算结果误差非常小，完全可以忽略不计，从而验证了Fortran程序的可行性。
（2）采用八结点单元结点位移。
选取结点2，结点20和结点86.
Fortran位移Abaqus位移相对误差4结点X -0.009762 -0.009735 0.28%
Y -0.007509 -0.007550 0.54% 20结点X -0.01042 -0.010403 0.16%
Y -0.2038 -0.20348 0.15% 86结点X 0.0002880 0.0002879 0.035%
Y -0.2039 -0.2035 0.20%
可以看出，采用八结点单元的相对误差仍然很小，再次说明Fortran程序的可行性。
另外，可以看出，采用的abaqus计算也存在一定的误差，当采用四节点单元时与采用八结点单元（即二次单元）时，以20结点为例，它们之间也存在着（-0.009580+0.0104028）/0.0104028*100%=7.9%，这样的误差也不算小了，所
以，在实际计算时，一定合理的划分网格，采用合适的单元形式，这样才能得到与正解相近的正确答案。
终上所述，Fortran 的程序是可行的。
5.结构变形图。
采用四节点单元的变形图：
采用八结点单元变形图：
说明：颜色越深的地方应力越大。可以看出，在悬臂梁上下两端的应力最大，越远离固定端应力越小，而在中性轴上应力很小，基本为零，验证了材料力学的结论。
图4.6 四节点单元变形图
图4.7 八结点单元变形图
四、算例
考虑如下一个弹性力学问题：
不计体力，长度150mm,记为a,宽度100mm,记为b 。./1,/N 2q 21mm N q mm ==材料参数E=200GPa,v=0.3.
1.先计算其理论解。
采用弹性力学的变分解法。
设 ?
??????+++=?????+++=)
()(23212321y B y B B y v x A x A A x u 里茲法求解，取一阶精度，即???==by
v Ax
u
应变能：()()uAB B A v Eab dxdy uAB B A v E a b
2)1(22)1(2V 222
002
22+++=+++=??ε 变分：
101)(V abq ady q uds f A b
x -=-==????ε
202)(B V abq bdx q vds f a
y -=-==????ε
引入应力边界条件：
??
?==-====-==dx
ds b v q f b y dy
ds a u q f a x y x ,,,,21时，时， 带入变分得到的上式： .,A 1
22121E uq q B E uq q ab
q B
V
ab q A V --=--=??????-=??-=??εε
于是有：
E
u
q q v E x q x q 1221,u --=-=
分别计算11、19和26结点的结点位移，方便后面的比较。 将各个节点的坐标带入以上计算式
11结点：u=-0.0002013mm,v=-0.0001345mm. 19结点：u=-0.0005034mm,v=-0.0001925mm. 26结点：u=-0.0006724mm,v=-0.0002732mm.
2.划分网格，定义单元，计算结点位移。
（1）四节点单元。
首先划分35个结点，24个单元。
图5.1 四节点单元定义
INPUT.DAT 文件与悬臂梁输入方法类似，这里直接给出用Fortran 和abaqus 计算出的OUTPUT.DAT 文件的结点位移是：（部分）
（2）八结点单元。
划分了93个结点，24个单元。
同理其OUTPUT.DAT 结点位移输出为：（部分）
图5.3 八节点单元定义
Fortran 计算结果
Abaqus 计算结果
图5.2 四节点单元结点位移Fortran 和abaqus 比较
3.变形情况。
变形云图：
说明：云图显示的是变形前和变形后的图形。颜色越深，应力越大。
Fortran 计算结果
Abaqus 计算结果
图5.4 Fortran 、abaqus 八结点单元结点位移对比
四节点单元
八节点单元
图5.5 四、八结点单元变形云图
采用减缩积分方式计算得到的云图：
细化网格，划分1500个网格，分别观察其使用四节点单元及其用减缩积分方式、八结点单元及其减缩积分方式变形情况。
4.
分析比较。
随机选取了三个节点进行分析，选取结点11
、19
和26结点。
四节点单元
四节点单元减缩积分
八节点单元
八节点单元减缩积分
图5.7 不同单元，不同计算方式比较
四节点减缩积分
八节点减缩积分
图5.6 四节点、八结点减缩积分云图比较
（1）四节点单元与八结点单元精度比较。
四节点单
元八结点单元精确解相对误差A B
四节点八结点
11结点X -0.0001792 -0.0001818 -0.0002013 10.97% 9.68% 1.43% Y -0.0001168 -0.0001150 -0.0001345 13.16% 14.49% 1.57%
19结点X -0.0004513 -0.0004630 -0.0005034 10.35% 8.02% 2.53% Y -0.0001770 -0.0001698 -0.0001925 8.05% 11.79% 4.24%
26结点X -0.0005966 -0.0006058 -0.0006724 11.27% 9.90% 1.51% Y -0.0002275 -0.0002128 -0.0002732 16.72% 22.10% 6.90%
说明，相对误差A是四节点或八结点与精确解之间的误差，相对误差B是四节点与八结点之间的误差。
分析：
可以看出采用四节点单元和八结点单元之间的误差都在百分之几，误差较小，说明采用四节点单元和八结点单元计算的结果是有区别的，所以在划分网格，选取单元是时候要慎重。
而四节点或八结点与精确解之间的误差达到了百分之二十多，误差相对较大。导致出现这种情况的原因可能是：a.在采用里茲法计算精确解的时候，为了简化计算，只选取了一阶精度的位移函数，得到的精确解与实际的精确解有较大的误差。b.由于电脑的计算都是数值计算，而且网格划分的比较粗糙，计算结果也存在着一定的误差。
在上面的数据中，还能总结出，无论采用四节点单元还是采用八结点单元，计算的结果总是小于精确解，这就验证了用最小位能原理得到的位移总是趋于真解下限的结论。
(2)四节点单元Fortran计算和abaqus计算比较。
Fortran abaqus 相对误差11结点X -0.0001792 -0.000179183 0.0095%
. 一、题目 如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。已知：P=150N/m ，E=200GPa ，=0.25，t=0.1m ，忽略自重。试计算薄板的位移及应力分布。 要求： 1. 编写有限元计算机程序，计算节点位移及单元应力。（划分三角形 单元，单元数不得少于30个）； 2. 采用有限元软件分析该问题（有限元软件网格与程序设计网格必 须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）； 3. 提交程序编写过程的详细报告及计算机程序； 4. 所有同学参加答辩，并演示有限元计算程序。 有限元法中三节点三角形分析结构的步骤如下： 1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。 2）单元分析，建立单元刚度矩阵。 3）整体分析，建立总刚矩阵。 4）建立整体结构的等效节点荷载和总荷载矩阵 5）边界条件处理。 6）解方程，求出节点位移。 7）求出各单元的单元应力。 8）计算结果整理。 一、程序设计 网格划分 如图，将薄板如图划分为6行，并建立坐标系，则
刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种； （2）该单元对应总刚度矩阵的那几行哪几列 （3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列 循环又分为3层循环：（1）最外层：逐行计算 （2）中间层：该行逐个计算 （3）最里层：区分为第 奇/偶 数个计算 单元刚度的集成：[ ][][][][][]' '''''215656665656266256561661e Z e e e Z e Z e e e e k k k K k k k k k k +?++=? =?==?==?=?????? 边界约束的处理：划0置1法 X Y P X Y P
一．问题描述 如图所示的平面矩形结构，设E=1，NU=0.25，h=1，考虑以下约束和外载： 位移边界条件BC（u）：U A=0，V A=0，U D=0， 力边界条件BC（p）：在CD边上有均布载荷q=1， 建模情形：使用四个四节点矩形单元， 试在该建模情形下，求各节点的位移以及各个单元的应力分布。
二．Matlab程序 (1).函数定义： function k= Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID) syms s t; a = (yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s))/4; b = (yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4; c = (xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4; d = (xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s))/4; B1 = [a*(t-1)/4-b*(s-1)/4 0 ; 0 c*(s-1)/4-d*(t-1)/4 ; c*(s-1)/4-d*(t-1)/4 a*(t-1)/4-b*(s-1)/4]; B2 = [a*(1-t)/4-b*(-1-s)/4 0 ; 0 c*(-1-s)/4-d*(1-t)/4 ; c*(-1-s)/4-d*(1-t)/4 a*(1-t)/4-b*(-1-s)/4]; B3 = [a*(t+1)/4-b*(s+1)/4 0 ; 0 c*(s+1)/4-d*(t+1)/4 ; c*(s+1)/4-d*(t+1)/4 a*(t+1)/4-b*(s+1)/4]; B4 = [a*(-1-t)/4-b*(1-s)/4 0 ; 0 c*(1-s)/4-d*(-1-t)/4 ; c*(1-s)/4-d*(-1-t)/4 a*(-1-t)/4-b*(1-s)/4]; Bfirst = [B1 B2 B3 B4]; Jfirst = [0 1-t t-s s-1 ; t-1 0 s+1 -s-t ; s-t -s-1 0 t+1 ; 1-s s+t -t-1 0]; J = [xi xjxmxp]*Jfirst*[yi ;yj ; ym ; yp]/8; B = Bfirst/J; if ID == 1 D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2]; elseif ID == 2 D = (E/(1+NU)/(1-2*NU))*[1-NU NU 0 ; NU 1-NU 0 ; 0 0 (1-2*NU)/2]; end BD = J*transpose(B)*D*B; r = int(int(BD, t, -1, 1), s, -1, 1); z = h*r; k = double(z); end function z = Quad2D4Node_Assembly(KK,k,i,j,m,p) DOF(1)=2*i-1; DOF(2)=2*i; DOF(3)=2*j-1; DOF(4)=2*j; DOF(5)=2*m-1; DOF(6)=2*m; DOF(7)=2*p-1; DOF(8)=2*p; for n1=1:8
结构有限元分析程序设计 绪论 §0.1 开设“有限元程序设计”课程的意义和目的 §0.2 课程特点 §0.3 课程安排 §0.4 课程要求 §0.5 基本方法复习 $0.1 意义和目的 1.有限元数值分析技术本身要求工程设计研究人员掌握 1). 有限元数值分析技术的完善标志着现代计算力学的真正成熟和实用化,已在各种 力学中得到了广泛的应用。比如：，已杨为工程结构分析中最得以收敛的技术手段，现代功用大致有： a). 现代结构论证。对结构设计从内力，位移等方面进行优劣评定，从而进 行结构优化设计。 b)可取代部份实验，局部实验+有限元分析，是现代工程设计研究方法的一大 特点。 c)结构的各种功能分析（疲劳断裂，可靠性分析等）都以有限元分析工具作为 核心的计算工具。 2). 有限元数值分析本身包括着理论+技术实现（本身功用所绝定的） 有限元数值分析本身包括着泛函理论+分片插值函数+程序设计 2. 有限元分析的技术实现（近十佘年的事）更依赖于计算机程序设计 有限元分析的技术取得的巨大的成就，从某种意义上说，得益于计算机硬件技术的发展和程序设计技术的发展，这两者的依赖性在当代表现得更加突出。（如可视化技术） 3.从学习的角度，不仅要学习理论，而且要从程序设计设计角度对这些理论的技术实现有 一个深入的了解，应当致力于掌握这些技术实现能力，从而开发它，发展它。（理论本身还有待于进一步完美相应的程序设计必须去开发） 4.程序设计不仅是实现有限元数值分析的工具和桥梁，而且在以下诸方面也有意义： 1). 精通基本概念，深化理论认识； 2). 锻炼实际工程分析，实际动手的能力； 3). 获得以后工作中必备的工具。（作业+老师给元素库） 目的：通过讲述有限元程序设计的技术与技巧，便能达到自编自读的能力。 §0.2 课程特点 总描述：理论+算法+数据结构（程序设计的意义） 理论：有限元算法，构造，步骤，解的等外性，收敛性，稳定性，误差分析 算法；指求解过程的技术方法，含两方面的含义；a. 有限元数值分析算法，b, 与数据结构有关的算法（总刚稀疏存贮，提取，节点优化编号等） 数据结构：指各向量矩阵存贮管理与实现，辅助管理结构（指针，数据记录等） 具体特点： 理论性强：能量泛函理论+有限元构造算法+数据结构构造算法 内容繁杂：理论方法+技术方法+技术技巧 技巧性强：排序，管理结构（指针生成，整型运算等）
重庆大学本科学生课程设计任务书 课程设计题目有限元程序设计 学院资源及环境科学学院专业工程力学年级2010级 已知参数和设计要求： 1.独立完成有限元程序设计。 2.独立选择计算算例，并能通过算例判断程序的正确性。 3.独立完成程序设计报告，报告内容包括理论公式、程序框图、程序本 体、计算算例，算例结果分析、结论等。 学生应完成的工作： 1.复习掌握有限单元法的基本原理。 2.掌握弹性力学平面问题3节点三角形单元或4节点等参单元有限元方法 的计算流程，以及单元刚度矩阵、等效节点载荷、节点应变、节点应力 和高斯积分等的计算公式。 3.用Fortran语言编写弹性力学平面问题3节点三角形单元或4节点等参 单元的有限元程序。 4.在Visual Fortran 程序集成开发环境中完成有限元程序的编辑和调试 工作。 5.利用编写的有限元程序，计算算例，分析计算结果。 6.撰写课程设计报告。 目前资料收集情况（含指定参考资料）： 1.王勖成，有限单元法，北京：高等教育出版社，2002。 2.O.C. Zienkiewicz, R. L. Taylor, Finite Element Method, 5th Eition, McGraw-Hall Book Company Limited, 2000。 3.张汝清，董明，结构计算程序设计，重庆：重庆大学出版社，1988。 课程设计的工作计划： 1.第1周星期一上午：教师讲解程序设计方法，程序设计要求和任务安 排。 2.第1周星期一至星期二完成程序框图设计。 3.第1周星期三至第2周星期四完成程序设计。 4.第2周星期五完成课程设计报告。 任务下达日期 2013 年 6 月 6 日完成日期 2013 年 07 月 03 日 指导教师（签名） 学生（签名）
湖南大学 《有限单元法》编程大作业 专业：土木工程 姓名： 学号： 2013年12月
目录 程序作业题目： (3) 1、程序编制总说明 (3) 2、Matlab程序编制流程图 (3) 3、程序主要标示符及变量说明 (4) 4、理论基础和求解过程 (5) 4.1、构造插值函数 (5) 4.2位移插值函数及应变应力求解 (5) 5.程序的验证 (6) 附录：程序代码 (15)
程序作业题目： 完成一个包含以下所列部分的完整的有限元程序( Project) 须提供如下内容的文字材料（1500字以上）： ①程序编制说明； ②方法的基本理论和基本公式； ③程序功能说明； ④程序所用主要标识符说明及主要流程框图； ⑤ 1～3 个考题：考题来源、输出结果、与他人成果的对比结果（误差百分比）； ⑥对程序的评价和结论（包括正确性、适用范围、优缺点及其他心得等）。 须提供源程序、可执行程序和算例的电子文档或文字材料。选题可根据各自的论文选题等决定。 1、程序编制总说明 a.该程序采用平面三角形等参单元，能解决弹性力学的平面应力、平面应变问题。 b.能计算单元受集中力的作用。 c.能计算结点的位移和单元应力。 d.考题计算结果与理论计算结果比较，并给出误差分析。 e.程序采用MATLAB R2008a编制而成。 2、Matlab程序编制流程图
图1 整个程序流程图 3、程序主要标示符及变量说明 1、变量说明： Node ------- 节点定义 gElement ---- 单元定义 gMaterial --- 材料定义，包括弹性模量，泊松比和厚度 gBC1 -------- 约束条件 gNF --------- 集中力 gk------------总刚 gDelta-------结点位移 输入结构控制参数 输入其它数据 形成整体刚度阵 引入支承条件 解方程，输出位移 求应力，输出应力 形成节点荷载向量 开始 结束 1 单元面积 求弹性矩阵 单元刚度矩阵 位移-应变矩阵 6 7 8 9 10 2 3 4 5
有限元大作业程序设计 学校：天津大学 院系：建筑工程与力学学院 专业：01级工程力学 姓名：刘秀 学号：\\\\\\\\\\\ 指导老师：
连续体平面问题的有限元程序分析 [题目]： 如图所示的正方形薄板四周受均匀载荷的作用，该结构在边界 上受正向分布压力， m kN p 1=，同时在沿对角线y 轴上受一对集中压 力，载荷为2KN ，若取板厚1=t ，泊松比0=v 。 [分析过程]： 由于连续平板的对称性，只需要取其在第一象限的四分之一部分参加分析，然后人为作出一些辅助线将平板“分割”成若干部分，再为每个部分选择分析单元。采用将此模型化分为4个全等的直角三角型单元。利用其对称性，四分之一部分的边界约束，载荷可等效如图所示。
[程序原理及实现]： 用FORTRAN程序的实现。由节点信息文件NODE.IN和单元信息文件ELEMENT.IN，经过计算分析后输出一个一般性的文件DATA.OUT。模型基本信息由文件为BASIC.IN生成。 该程序的特点如下： 问题类型：可用于计算弹性力学平面问题和平面应变问题 单元类型：采用常应变三角形单元 位移模式：用用线性位移模式 载荷类型：节点载荷，非节点载荷应先换算为等效节点载荷 材料性质：弹性体由单一的均匀材料组成 约束方式：为“0”位移固定约束，为保证无刚体位移，弹性体至少应有对三个自由度的独立约束 方程求解：针对半带宽刚度方程的Gauss消元法
输入文件：由手工生成节点信息文件NODE.IN，和单元信息文件ELEMENT.IN 结果文件：输出一般的结果文件DATA.OUT 程序的原理如框图：
有限元法分析与建模课程设计报告 学院：机械电子工程学院 专业：机械电子工程 指导教师：杜平安 学生：乔林 学号：201221080212 2012-12-10
摘要 摘要 连杆的作用是将活塞的往复运动变成曲轴的旋转运动, 并把活塞上的力传 给曲轴连杆工作的小端做往复运动, 大端作旋转运动, 杆身做复杂的平面运动。本文用Pro/E建立连杆的三维模型，并运用ANSYS强大的有限元分析和优化功能来实现连杆的分析ANSYS 是一款极其强大的有限元分析软件。通过数据接口，ANSYS 可以方便的实现从CAD 软件中导入实体模型。因此，将Pro/E强大的 建模功能与ANSYS 优越的有限元分析功能结合在一起可以极大地满足设计者 在设计过程中对建模与分析的需求。 关键词：连杆，有限元，Pro/E，ANSYS
ABSTRACT ABSTRACT The role oftheconnecting rodisthesmall end ofthereciprocation of the pistonintoarotational movementofthecrankshaft, and to transmittheforceon the pistontothecrankshaft connecting rodreciprocates, thebig endfor pivotal movement, Shaftdo complexplanar motion. The establishment ofalinkageof thethree-dimensionalmodelusingPro / E, thepowerfulANSYSfinite elementanalysis andoptimization capabilitiestoachievetheconnecting rodfatigueanalysisANSYSisan extremelypowerfulfinite element analysis software. Throughthedata interface, ANSYS canfacilitate the realization ofsolid modelsimportedfromCAD software. Therefore,thesuperiorpowerfulmodeling capabilitiesofPro / Eand ANSYSfinite elementanalysis capabilitiestogethercanmeetthedesignersin the design processmodelingand analysis. Keywords:rod, finite element, Pro / E, ANSYS
有限元分析软件 有限元分析是对于结构力学分析迅速发展起来的一种现代计算方法。它是50 年代首先在连续体力学领域--飞机结构静、动态特性分析中应用的一种有效的数值分析方法，随后很快广泛的应用于求解热传导、电磁场、流体力学等连续性问题。 有限元分析软件目前最流行的有：ANSYS、ADINA、ABAQUS、MSC 四个比较知名比较大的公司，其中ADINA、ABAQUS 在非线性分析方面有较强的能力目前是业内最认可的两款有限元分析软件，ANSYS、MSC 进入中国比较早所以在国内知名度高应用广泛。目前在多物理场耦合方面几大公司都可以做到结构、流体、热的耦合分析，但是除ADINA 以外其它三个必须与别的软件搭配进行迭代分析，唯一能做到真正流固耦合的软件只有ADINA。ANSYS是商业化比较早的一个软件，目前公司收购了很多其他软件在旗下。ABAQUS专注结构分析目前没有流体模块。MSC是比较老的一款软件目前更新速度比较慢。ADINA是在同一体系下开发有结构、流体、热分析的一款软件，功能强大但进入中国时间比较晚市场还没有完全铺开。 结构分析能力排名：ABAQUS、ADINA、MSC、ANSYS 流体分析能力排名：ANSYS、ADINA、MSC、ABAQUS 耦合分析能力排名：ADINA、ANSYS、MSC、ABAQUS 性价比排名：最好的是ADINA，其次ABAQUS、再次ANSYS、最后MSC ABAQUS 软件与ANSYS 软件的对比分析： 1．在世界范围内的知名度：两种软件同为国际知名的有限元分析软件，在世界范围内具有各自广泛的用户群。ANSYS 软件在致力于线性分析的用户中具有很好的声誉，它在计算机资源的利用，用户界面开发等方面也做出了较大的贡献。ABAQUS软件则致力于更复杂和深入的工程问题，其强大的非线性分析功能在设计和研究的高端用户群中得到了广泛的认可。由于ANSYS 产品进入中国市场早于ABAQUS，并且在五年前ANSYS 的界面是当时最好的界面之一，所以在中国，ANSYS 软件在用户数量和市场推广度方面要高于ABAQUS。但随着ABAQUS北京办事处的成立，ABAQUS软件的用户数目和市场占有率正在大幅度和稳步提高，并可望在今后的几年内赶上和超过ANSYS。 2．应用领域：ANSYS 软件注重应用领域的拓展，目前已覆盖流体、电磁场和多物理场耦合等十分广泛的研究领域。ABAQUS 则集中于结构力学和相关领域研究，致力于解决该领域的深层次实际问题。 3．性价比：ANSYS 软件由于价格政策灵活，具有多种销售方案，在解决常规的
目录一. 前言 二．有限元设计部分 1 问题阐述 2 解析法求解 3 模型简化 4 ANSYS软件应用说明 5 结果分析 三．机械优化设计部分 1 问题阐述 2 解析算法 3 黄金分割法顺序流程图 4 C语言源程序代码 5 结果分析 四．设计心得 五. 参考文 一．前言
二．有限元设计部分 1、问题阐述 外伸梁上均布载荷的集中度为q=3kN/m，集中力偶矩M e=3kN·m列出剪力方程和弯矩方程，并绘制剪力图。材料力学Ι（刘鸿文第四版）P121
图2-1 外伸梁简化图 2、解析法求解 由梁的平衡方程，求出支反力为 F RA=14.5kN，F RB=3.5kN 梁的C A、ＡＤ、ＤＢ等三段内，剪力和弯矩都不能有同一个方程来表示，所以应分为三段考虑。对每一段都可以用同一个方法计算，列出剪力方程和弯矩方程，方程中ｘ以ｍ为单位，Ｆｓ（ｘ）以ｋＮ为单位，Ｍ（ｘ）以ｋＮ为单位。 在ＣＡ段内： Ｆｓ（ｘ）＝－ｑｘ＝－３ｘ（０＜=ｘ＜２ｍ） （g） Ｍ（ｘ）＝－（3／２）Ｘ２（０＜ｘ＜=2ｍ） （h）
在AD段内： Ｆｓ（ｘ）＝F RA－ｑｘ＝14.5－３ｘ（2m＜ｘ＜=6ｍ） （i） Ｍ（ｘ）＝F RA（x-2）-（1/2）X２=14.5（x-2）- （3/2）X 2 (j) (2m＜ ｘ6m) M（x）是x的二次函数,根据极值条件dM（x）/d(x)=0,得 14.5-3x=0 由此解出x=4.83m，亦即在那这一截面上，弯矩为极值。 代入（j） 式得AD段内的最大弯矩为 M=6.04kN·m 当截面取在DB段，用截面右侧的外力计算剪力和弯矩比较方
有限元分析软件的比较（购买必看）－转贴 随着现代科学技术的发展，人们正在不断建造更为快速的交通工具、更大规模的建筑物、更大跨度的桥梁、更大功率的发电机组和更为精密的机械设备。这一切都要求工程师在设计阶段就能精确地预测出产品和工程的技术性能，需要对结构的静、动力强度以及温度场、流场、电磁场和渗流等技术参数进行分析计算。例如分析计算高层建筑和大跨度桥梁在地震时所受到的影响，看看是否会发生破坏性事故；分析计算核反应堆的温度场，确定传热和冷却系统是否合理；分析涡轮机叶片内的流体动力学参数，以提高其运转效率。这些都可归结为求解物理问题的控制偏微分方程式，这些问题的解析计算往往是不现实的。近年来在计算机技术和数值分析方法支持下发展起来的有限元分析（FEA，Finite Element A nalysis）方法则为解决这些复杂的工程分析计算问题提供了有效的途径。在工程实践中，有限元分析软件与CAD系统的集成应用使设计水平发生了质的飞跃，主要表现在以下几个方面： 增加设计功能，减少设计成本； 缩短设计和分析的循环周期； 增加产品和工程的可靠性； 采用优化设计，降低材料的消耗或成本； 在产品制造或工程施工前预先发现潜在的问题； 模拟各种试验方案，减少试验时间和经费； 进行机械事故分析，查找事故原因。 在大力推广CAD技术的今天，从自行车到航天飞机，所有的设计制造都离不开有限元分析计算，FEA在工程设计和分析中将得到越来越广泛的重视。国际上早20世纪在50年代末、60年代初就投入大量的人力和物力开发具有强大功能的有限元分析程序。其中最为著名的是由美国国家宇航局（NASA）在1965年委托美国计算科学公司和贝尔航空系统公司开发的NASTRAN有限元分析系统。该系统发展至今已有几十个版本，是目前世界上规模最大、功能最强的有限元分析系统。从那时到现在，世界各地的研究机构和大学也发展了一批规模较小但使用灵活、价格较低的专用或通用有限元分析软件，主要有德国的ASKA、英国的PA FEC、法国的SYSTUS、美国的ABQUS、ADINA、ANSYS、BERSAFE、BOSOR、COSMOS、ELAS、MARC和STARDYNE等公司的产品。 以下对一些常用的软件进行一些比较分析： 1. LSTC公司的LS-DYNA系列软件
有限元法分析与建模课程设计报告 学院：机电学院 专业：机械制造及其自动化指导教师：**** 学生：* *** 学号：2012011**** 2015-12-31
摘要 本文通过ANSYS10.0建立了标准光盘的离心力分析模型，采用有限元方法对高速旋转的光盘引起的应力及其应变进行分析，同时运用经典弹性力学知识来介绍ANSYS10.0中关于平面应力问题分析的基本过程和注意事项。力求较为真实地反映光盘在光驱中实际应力和应变分布情况，为人们进行合理的标准光盘结构设计和制造工艺提供理论依据。 关键词：ANSYS10.0；光盘；应力；应变。
目录 第一章引言 (3) 1.1 引言 (3) 第二章问题描述 (4) 2.1有限元法及其基本思想 (4) 2.2 问题描述 (4) 第三章力学模型的建立和求解 (5) 3.1设定分析作业名和标题 (5) 3.2定义单元类型 (6) 3.3定义实常数 (9) 3.4定义材料属性 (12) 3.5建立盘面模型 (14) 3.6对盘面划分网格 (22) 3.7施加位移边界 (27) 3.8施加转速惯性载荷并求解 (30) 第四章结果分析 (32) 4.1 旋转结果坐标系 (32) 4.2查看变形 (33) 4.3查看应力 (35) 总结 (38) 参考文献 (39)
第一章引言 1.1 引言 光盘业是我国信息化建设中发展迅速的产业之一，认真研究光盘产业的规律和发展趋势，是一件非常迫切的工作。光盘产业发展的整体性强，宏观调控要求高，因此，对于光盘产业的总体部署、合理布局和有序发展等问题，包括节目制作、软件开发、硬件制造、节目生产、技术标准等。 在高速光盘驱动器中，光盘片会产生应力和应变，在用ANSYS分析时，要施加盘片高速旋转引起的惯性载荷，即可以施加角速度。需要注意的是，利用ANSYS施加边界条件时，要将内孔边缘节点的周向位移固定，为施加周向位移，而且还需要将节点坐标系旋转到柱坐标系下。 本文通过ANSYS10.0建立了标准光盘的离心力分析模型，采用有限元方法对高速旋转的光盘引起的应力及其应变进行分析，同时运用经典弹性力学知识来介绍ANSYS10.0中关于平面应力问题分析的基本过程和注意事项。
有限元编程的c++实现算例 1. #include<> 2. #include<> 3. 4. 5. #define ne 3 #define nj 4 #define nz 6 #define npj 0 #define npf 1 #define nj3 1 2 #define dd 6 #define e0 #define a0 #define i0 #define pi 16. 17. 18. int jm[ne+1][3]={{0,0,0},{0,1,2},{0,2,3},{0,4,3}}; /*gghjghg*/ 19. double gc[ne+1]={,,,}; 20. double gj[ne+1]={,,,}; 21. double mj[ne+1]={,a0,a0,a0}; 22. double gx[ne+1]={,i0,i0,i0}; 23. int zc[nz+1]={0,1,2,3,10,11,12}; 24. double pj[npj+1][3]={{,,}}; 25. double pf[npf+1][5]={{0,0,0,0,0},{0,-20,,,}}; 26. double kz[nj3+1][dd+1],p[nj3+1]; 27. double pe[7],f[7],f0[7],t[7][7]; 28. double ke[7][7],kd[7][7]; 29. 30. 31. 36. void jdugd(int); 38. void zb(int); 39. void gdnl(int); 40. void dugd(int);
有限元编程的c++实现算例 1.#include
有限元课程设计
目录 0.前言 (3) 1.问题阐述 (4) 2.有限元分析 (5) 2.1.梁的参数设定 (5) 2.2.材料参数 (5) 2.3.单元选择 (5) 2.4.梁的边界条件 (6) 2.5.梁所受的载荷 (6) 2.6.ANSYS软件应用说明 (6) 3.交互式的求解过程 (7) 3.1创建梁的各个节点 (7) 3.2定义单元类型、材料特性和梁的横截面几何参数 (10) 3.3创建单元 (12) 4.施加约束和载荷 (13) 4.1节点自由度约束 (13) 4.2施加节点13处的弯矩m。 (14) 4.3施加单元1到单元12上的的分布载荷q。 (15) 5.求解 (15) 5.1定义分析类型 (15) 5.2求解 (15) 6.后处理 (16) 6.1绘制梁的Y方向变形图 (16) 6.2建立单元结果表 (17) 6.3结果显示 (19) 退出程序 (21) 心得体会 (22) 参考文献 (22)
0.前言 目前，几乎所有高校的力学、土木、机械、航空、航天、船舶、水利、交通、桥梁等理工科专业，都为高年级本科生开设了《有限元方法》基础课程，为研究生开设了《非线性有限元方法》学位课程。学生在学习完有限元课程之后，还必须熟练掌握相关有限元软件的使用，才能将有限元基本理论有效地应用到实际工程问题分析中去。为此，部分有条件的高校也开设了有限元软件应用课程（课程名称可能会因学校及专业的不同而有所差异，但都是以讲解有限元软件ANSYS或其他软件为主）。哈尔滨工业大学航天学院工程力学专业20世纪90年代末即开设了该类课程《应用软件工程--ANSYS》，作者从2003年开始接手讲授该门课程。虽然市面上的ANSYS书籍很多，但却难以找出一本非常适合做教材的书籍，因此作者参考多本书籍自主编写了校内讲稿。经过6年多的试用，目前已基本成型，现将多年的校内讲稿和心得体会完善成书，以期与开设该类课程的兄弟院校分享、共勉，同时也供从事相关科研与工程项目的人员参考阅读。 ANSYS软件是目前国际上最著名的大型通用有限元分析软件，经过三十年的发展，已形成融结构、热、流体、电磁、声学及多物理场耦合为一体的大型通用有限元分析软件，广泛应用于航空航天、石油、化工、汽车、造船、铁道、电子、机械制造、地矿能源、水利、核能、生物、医学、土木工程、轻工、一般工业及科学研究等各个领域，其极强的分析功能覆盖了几乎所有的工程问题。作为世界最具权威的有限元产品和工业化分析标准，目前几乎所有的CAD/CAE/CAM软件都竞相开发了与ANSYS的专用接口，实现数据的共享和交换，如Pro/Engineer、NASTRAN、Alogor、I-DEAS及AutoCAD等。ANSYS软件在Linux 和Windows下均有版本，并同时有32位和64位版本，目前最新的版本为12.0。 本书以ANSYS 12.0版本为依据，以Windows NT为操作平台，将结构有限元分析的基本理论与ANSYS实践操作紧密结合，通过大量精心筛选的具有实际工程应用背景的原创性分析实例，以图形用户界面和命令流两种方式向读者全面介绍了ANSYS结构有限元分析方法。
有限元及程序设计 1. 下列关于高精度单元描述正确的是（ ）。 A.等参元的位移模式和坐标变换采用不同的形函数 B.矩形单元形状规则，因而使用范围较 广 C.6结点三角形单元、10结点三角形单元、8结点矩形单元和 12结点矩形单元的单元 刚度矩阵的建立过程是不一样的 D.6结点三角形单元较容易模拟物体的边界形状 【参考答案】：D 2. 0 =cxy 能解决矩形板（ ）问题。 A.左右均布拉压 B.上下均布拉压 C 纯剪切 D.纯弯曲 【参考答案】：C 3. 下列关于等参元的叙述不正确的是（ A.精度较高 B.能较好的模拟边界条件 【参考答案】：D 4. 薄板的边界不包括（ A.简支边界 B.固定边界 C.自由边界和荷载边界 D.非固定边界 【参考答案】：D ）。 C.输入的信息量较少 D.输入的信息量较多 5. 下列属于平面应力问题的是（ A.平板坝的平板支墩 B 挡土墙 【参考答案】：A 6. 在应力函数上任意增减一个（ A.线性项 B.二次项 C.三次项 【参考答案】：A 7. 下列不属于提高单元精度的方法是 A.增加单元结点数目 B 在单元内增设结点 【参考答案】：C 8. 空间问题的基本平衡微分方程有（ A.2 B.3 C.4 D.5 【参考答案】：C ）。 C.重力水坝 D.受内水压力作用的圆管 ），对应力分量无影响。 D.常数项 ）。 C.减少单元结点数目 D.设等参元 个。 ）问题。 C 纯剪切 D.纯弯曲 A.左右均布拉压 B.上下均布拉压 【参考答案】：B 10. 下列属于不规则单元的有（ A.正四面体单元 B.正三棱体单元 【参考答案】：C 11. 空间问题的基本未知位移分量有（ A.2 B.3 C.4 D.5 【参考答案】：B 1. 薄板小挠度弯曲理论的基本假定是（ ）。 A.直法线假定 B.法向位移假定 C.中面位移假定 D.板内无挤压假定 【参考答案】：A|C|D 2. 弹性力学平面问题 按应力求解具体可分为（ ）两种。 A.逆解法 B.顺解法 C.半逆解法 D.半顺解法 ）。 C.任意四面体单元 D.正六面体单元 ）个。
成绩评定表
课程设计任务书
目录 一、前言 (2) 二、平面梁结构的静力学分析 (3) 1. 问题阐述 (3) 2. 物理参数与几何参数 (3) 三、交互式的求解过程 (3) 1. 创建节点 (3) 2. 定义单元类型和材料特性 (6) 3. 创建单元 (7) 4. 施加约束和载荷 (8) 5. 求解 (10) 6. 后处理 (11) 7. 退出程序 (15) 参考文献 (15)
前言 有限元分析（FEA，Finite Element Analysis）的基本概念是用较简单的问题代替复杂问题后再求解。它将求解域看成是由许多称为有限元的小的互连子域组成，对每一单元假定一个合适的(较简单的）近似解，然后推导求解这个域总的满足条件(如结构的平衡条件），从而得到问题的解。这个解不是准确解，而是近似解，因为实际问题被较简单的问题所代替。由于大多数实际问题难以得到准确解，而有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。 有限元是那些集合在一起能够表示实际连续域的离散单元。有限元的概念早在几个世纪前就已产生并得到了应用，例如用多边形（有限个直线单元）逼近圆来求得圆的周长，但作为一种方法而被提出，则是最近的事。有限元法最初被称为矩阵近似方法，应用于航空器的结构强度计算，并由于其方便性、实用性和有效性而引起从事力学研究的科学家的浓厚兴趣。经过短短数十年的努力，随着计算机技术的快速发展和普及，有限元方法迅速从结构工程强度分析计算扩展到几乎所有的科学技术领域，成为一种丰富多彩、应用广泛并且实用高效的数值分析方法。 有限元方法与其他求解边值问题近似方法的根本区别在于它的近似性仅限于相对小的子域中。20世纪60年代初首次提出结构力学计算有限元概念的克拉夫（Clough）教授形象地将其描绘为：“有限元法=Rayleigh Ritz法＋分片函数”，即有限元法是Rayleigh Ritz法的一种局部化情况。有限元法将函数定义在简单几何形状的单元域上，且不考虑整个定义域的复杂边界条件，这是有限元法优于其他近似方法的原因之一。 随着市场竞争的加剧，产品更新周期愈来愈短，企业对新技术的需求更加迫切，而有限元数值模拟技术是提升产品质量、缩短设计周期、提高产品竞争力的一项有效手段，所以，随着计算机技术和计算方法的发展，有限元法在工程设计和科研领域得到了越来越广泛的重视和应用，已经成为解决复杂工程分析计算问题的有效途径，从汽车到航天飞机几乎所有的设计制造都已离不开有限元分析计算，其在机械制造、材料加工、航空航天、汽车、土木建筑、电子电器、国防军工、船舶、铁道、石化、能源和科学研究等各个领域的广泛使用已使设计水平发生了质的飞跃。
P 9 m 9 m 一、题目 如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。已知：P=150N/m，E=200GPa，=，t=，忽略自重。试计算薄板的位移及应力分布。 要求： 1.编写有限元计算机程序，计算节点位移及单元应力。（划分三角形 单元，单元数不得少于30个）； 2.采用有限元软件分析该问题（有限元软件网格与程序设计网格必 须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）； 3.提交程序编写过程的详细报告及计算机程序； 4.所有同学参加答辩，并演示有限元计算程序。 有限元法中三节点三角形分析结构的步骤如下： 1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。 2）单元分析，建立单元刚度矩阵。 3）整体分析，建立总刚矩阵。 4）建立整体结构的等效节点荷载和总荷载矩阵 5）边界条件处理。 6）解方程，求出节点位移。 7）求出各单元的单元应力。 8）计算结果整理。 一、程序设计
网格划分 如图，将薄板如图划分为6行，并建立坐标系，则 X Y P X Y P 节点编号 单元编号
刚度矩阵的集成 建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。 通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种； （2）该单元对应总刚度矩阵的那几行哪几列 （3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列 循环又分为3层循环：（1）最外层：逐行计算 （2）中间层：该行逐个计算 （3）最里层：区分为第奇/偶数个计算 单元刚度的集成： [][] [][] [][] ' ' ' ' ' ' 2 1 56 56 6 6 56 56 2 6 6 2 56 56 1 6 6 1 e Z e e e Z e Z e e e e k k k K k k k k k k + ? + + = ? = ? = = ? = = ? = ? ? ? ? ? ? 边界约束的处理：划0置1法 适用：这种方法适用于边界节点位移分量为已知(含为0)的各种约束。 做法： （1）将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1，其他元素改为零；同 时将载荷列阵｛R｝中相应元素用已知位移置换。 ◎这样，由该方程求得的此位移值一定等于已知量。 （2）将〔K〕中已知位移相应的列的非主对角成元素也置0，以保持〔K〕的对称性。 ◎当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为 保证方程成立，须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。◎若约束为零位移约束时，此步则可省去。 特点： （1）经以上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。 （2）但这种方法不改变方程阶数，利于存贮。
本科生实验报告书 四节点等参单元有限元分析的FORTRAN程序
目录 1．问题概述 (1) 2．四节点四边形等参单元介绍 (1) 3．单元应力磨平方法介绍 (4) 4．程序流程设计 (6) 程序设计概述 流程图 5．程序结构及程序说明 (8) 6．程序应用及算例分析 (9) 算例概述 算例ANSYS求解 算例程序数值解 算例分析 7. 总结 (15)
1. 问题概述 等参单元是有限元方法中使用最广泛的单元类型。等参单元的位移模式和坐标变换均采用相同的形函数，这种坐标变换叫做等参变换。通过等参变换可以将自然（局部）坐标中几何形状规则的单元转换成总体（笛卡尔）坐标中形状扭曲的单元，因而使得单元有较好的适应性。 本问题首先对平面四节点四边形等参单元的形函数、应力矩阵和等效节点力矩阵、应力磨平公式等的推导和计算求解。并通过设计FORTRAN 求解程序进行编程求解，最后给出算例（受集中荷载的悬臂梁）并进行求解，将解与ANSYS 的解进行比较。在这个过程中，采用了高斯三点积分和高斯两点积分，这种积分方法的求解效率较高而且精度也较好。在问题的最后，尝试去分析引起数值解误差的原因，并分析四节点等参单元的若干特性。 2. 四节点四边形等参单元介绍 边长为2的正方形单元（如下图所示），在其形心处安置一个局部坐标。单元角 结点i 的坐标（,）分别为 ,因此单元四条边界的方程可用简单公式 和 逐一给出。 图2-1 母单元 1 2 3 4 0
图2-2 四边形单元 形函数的表达式： 位移函数： ∑==41 i N u i i u ∑==4 1 i N v i i v 坐标变换式： ∑==41 i x N x i i ，∑==4 1 i y N y i i 单元应变矩阵 {}[]{}[]{}e e B B B B B x v y u y v x u δδε4321==??? ? ??? ?????????????+??????= 式中{ }[]T T T T T 4 321e δδδδδ=——单元节点的位移列阵；
有限元程序设计报告 课程名称:有限元程序设计 指导教师：张亮 学校：重庆大学 专业：工程力学01班 姓名：苏世宏 学号：20126699 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）总结在编写程序和算例中遇到的问题和解决方法，写出自己的心得体会。 二、弹性力学平面问题有限元方法的基本公式