有限元三角形等参单元

北方工业大学

高等有限元课程总结

姓名:韩双鹏学号: 2011303310105 专业班级:结构研-11 系(部、院):建筑工程学院

2012 年5 月25 日

高等有限元学习总结——六节点三角形等参数单元

1 概述

从弹性力学基本方程到有限元原理再到最新进展，经过本课程的学习，比较系统的掌握了有限元相关内容，更学习到了一种方法、一些生活中的哲理。首先从大方向掌握所学内容，避免迷失在局部造成一叶遮目不见泰山之悲剧，比如弹性力学原理从大方向说就是三类方程，以及其在各类问题中的应用；其次了解了科研的相关过程及创新之处，从已知的东西到无知的领域，正如老师所说，能成功地把某一领域的东西搬到相关领域，这就是一大创造，比如有限元中将梁弯曲的理论研究厚板弯曲问题，由有限元标准单元到等参元的研究等；再有，我们生活中的常识、学习中的某些东西值得我们细细品味，也许这就是平时所说的小事反应大道理，老师的理论：“很多想法都是错误的”“很好想到的方法也许很难走通”“有缺陷的东西才更体现出美”“平衡的理论，吃点亏也许是福”等等，受益匪浅。不再一一赘述，本文将取其中的一个知识点，总结六节点三角形等参单元的相关内容。

我们知道，无论三节点或者六节点三角形单元还是四节点或者八节点矩形单元，它们形状简单、规则但计算精度低，且对于复杂边界的适应性差，难以很好的拟合曲边边界，解决这一问题的通用方法是细分边界，以直代曲，利用更多的简单单元去拟合边界复杂的区域。但这样处理仍存在折线代替曲线所带来的误差，且这种误差不能通过提高单元位移函数的精度来补偿。那么能否构造出单元形状任意、边界适应性好、计算精度高的曲边单元，以便在给定的精度下用较少数目的单元去解决实际问题？这就是有限元中一类重要的单元——等参数单元。本文将总结等参数单元的基本概念，并以六节点三角形单元为例讲述等参元实现过程中的三种变换，以及该等参元的收敛性等问题。

2 等参数单元及实现过程

2.1 等参数单元概念

由于实际问题的复杂性，通常需要使用一些形状不规整和形状复杂的单元来离散边界形状复杂的原问题。如下图所示（a）中为常见的几何形状不规整的实际单元，称为实际单元，也称为参数单元。（b）中为对应的形状规整的单元，称为标准单元。对于形状复杂的实际单元的单元分析，若仍采用前面介绍的方法进行，则在单元位移函数的建立和单元刚度矩阵计算方面会遇到许多困难。由此可考虑利用前面介绍过的形状规整的标准单元的单元分析来研究实际单元，几何形状的不同可认为是坐标变换的结果。

下面以六节点三角形单元为例来说明标准单元和实际单元间的坐标变换。首先回顾下六节点三角形单元。如下图所示：

单元节点位移向量为：

T

e

v u v u v u v u v u v u } {}{665544332211=δ 单元内任意一点的位移是单元节点位移的插值函数,已知u(x,y)和v(x,y)

在六个顶点的函数值，可分别设u(x,y)和v(x,y)为具有六个待定系数的插值多项式

22

12345622

123456u A A x A y A x A xy A y

v B B x B y B x B xy B y

?=+++++??=+++++?? 将已知条件代入，可解得其中的待定系数。

u(x,y)和v(x,y)也可用Lagrenge 插值公式表示为：

∑=+++++==6

1665544332211k k k u N u N u N u N u N u N u N u

∑=+++++==61

665544332211k k k v N v N v N v N v N v N v N v

写成矩阵形式：

????????

?

???????????????????????????????????=???

???66554433221165

4

3

2

1

654321

v u v u v u v u v u v u 0

000000N N N N N N N N N N N N v u 缩写为： {}e ][δδN v u =?

??

???=

单元内任意一点的应变

x y xy u x v y u v y x εεγ??=?????=???

???=+????

矩阵形式： {}{}e N H v u H x y y x δε]][[][v u 00=??????=???

????????

??

?

????????????????= e {}[B]{}εδ=缩写为：

[]12

365

43

2

1

65

4

3

2

1

654321

0000000 0

0]][[][?=??

????????????

????????????????==B B B B B B N N N N N N N N N N N N x y y x

N H B []??

??

???

?????????????????=x N y N y N x

N B k k k k

k 0

0 k = 1,2,3…6 N k (k=1,2,3…6)为六节点三角形单元的形函数，N k 的表达式：

2

6524321),(y a xy a x a y a x a a y x N k k k k k k k +++++= k =1,2,3 (6)

(,)1

(,)0

k k k k i i N x y i k N x y =?≠?

=? 故可得：

245356 2 2k k k k k

k k k N a a x a y x

N a a x a y y

?=++??=++? 单元内任意一点的应力：}]{][[}]{[}{e

B D D δεσ== 单元刚度矩阵： ?

=A

T

e

dxdy B D B t k ]][[][][

有了以上六节点三角形单元的知识，下面我们建立其等参元的概念。如下图所示，在母元中采用面积坐标（123,,L L L ），在子元中采用直角坐标（x,y ）。设两个坐标系间的变换关系为

123123(,,)

(,,)x x L L L y y L L L =??

=?

利用六个节点在两个坐标系间的一一对应关系，即：

123123(,,)

(1,2...6)(,,)i i i

i x x L L L i y y L L L =?=?

=?

很容易建立两个坐标系间的坐标变换关系。

对照六节点三角形单元的位移函数可知，坐标变换式与其具有完全相同的形式，故可将坐标变换式整理成插值函数的形式：

6

1231

6

1231(,,)(,,)i i i i i

i x N L L L x y N L L L y

==?

=????=??

∑∑ 其中(1,2,3)i i N L i == 4124N L L = 5324N L L = 6314N L L =与六节点三角形单元的位移插值函数完全相同。

上述坐标变换式建立了两个坐标系之间的坐标变换关系，从而可将局部坐标系下的形状规则的标准单元变换为整体坐标系下的形状复杂的实际单元。上述变换关系中插值基函数

i

N 和插值节点数目n 是两个关键的参数。若变换函数中的

插值基函数（即形函数）以及插值节点与描述单元位移函数的插值基函数及插值节点完全相同，则这种变换称为等参数变换，这种单元就称为等参单元。也就是说，等参单元的位移函数和坐标变换式具有相同的插值基函数（即形函数），且他们分别用同一节点的位移值和坐标值进行函数插值来表示单元内任意点的位移和几何坐标。

2.2 导数的坐标变换和微元变换

在进行等参单元的单元分析时，需要将对整体坐标的求导和积分运算变换为对局部坐标的求导和积分运算。下面讨论其变换关系。

2.2.1 导数的坐标变换

根据偏导数变换式：

ξξξ????+????=??y

y x x ()()() η

ηη????+????=??y

y x x ()()() 以上两式的矩阵形式：

?????

?????????=??

?

?????????????

??????????????????=??????????????????y x J y x y x y x ()()][()()()()ηηξξηξ

????

?

?

???

???????????=ηηξξy x y x J ][ 称为Jacobi 矩阵，由局部坐标与整体坐标间的变换函数确定。 显然若J 的逆存在，则任意函数对整体坐标的求导可转化为对局部坐标的求

导，即

???????

?

??????????????????????????-??-

??=??

??

??????????????=?????

?????????-ηξξη

ξηηξ()()1()()][()()1x x y y J J y x

在六节点三角形单元用面积坐标表示时，可令：

123,,1L L L ξηξη==--=

得偏导数之间的变换矩阵，进而可求得应变。

2.2.2微元变换

图中面积元的面积：

ηξ

d d d ?=Ω

j d y

d y i d x d x j dy i dx d

)()(ηη

ξξηηξξξ??+??+??+??=?+?=

在d ξ方向，d η=0, ∴ j d y i d x

d

???+???=ξξ

ξξξ 同理可得： j d y i d x

d

???+???=ηη

ηηη ()()x y x y d d i d j d i d j ξξηηξξηη

????Ω=?+???+?????故

注意到：1

0 0=?=?=?j i j j i i

(

)x y y x

d d d i j J d d ξηξηξηξη

????Ω=?-???=????故 将其用123,,1L L L ξηξη==--=替换，并带入刚度矩阵公式，可得该等参元的刚度矩阵为：

2

111200

e

L e T T K B DBtdA B DB J tdL dL -Ω==??

?

3 六节点三角形等参单元的收敛性

利用该等参单元求解问题是否收敛，还需看位移函数是否满足三个条件:（1）位移函数中包含有刚体位移（2）位移函数中包含有常应变（3）位移函数内部连续和边界协调。

由于刚体位移和常应变状态可写成下列形式:

012u=a +a x+a y 012v b b x b y

=++ 但上二式可化成下式:

i i u N u =∑ i i

v N v =∑

即：()0

12012i

i i i i i i i u N a

a x a y a N a N x a N y =

++=++∑∑∑∑

()012012i i i i i i i i v N b b x b y b N b N x b N y =++=++∑∑∑∑ 根据形函数的性质: ΣNi=1 则上二式都化成： i i

u N u =

∑ i i

v N v =∑

所以构成的位移函数包含刚体位移和常应变。

多项式在单元内部无奇异点，所以连续。相邻单元的边界上可由线上节点座标唯一确定，而母单元相邻边界都是直线，且形函数相同，所以交界线上的位移由节点位移唯一确定，相邻单元连续。

4等参单元的几点说明

（1）等参单元为协调元，满足有限元解收敛的充要条件。 （2）等参单元存在的充要条件是：为了保证能进行等参变换(即总体坐标与局部坐标一一对应)，通常要求总体坐标系下的单元为凸，即不能有内角大于或等于或接近180度情况。

（3）等参单元的优点是当单元边界呈二次以上的曲线时，容易用很少的单

元去逼近曲线边界。

（4）坐标变换时也可取描述单元几何形状的几何模式不是形函数的，如p-element.

5 总结

通过有限元的学习，为以后的相关软件的学习打下理论基础；这门课程的学习过程、学习方法可以应用于其他课程的学习，掌握了科研的思路和过程；日常生活、学习中养成善于思考的习惯，不断开阔新的思路。

平面三角形单元有限元程序设计

. 一、题目 如图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

有限元法的基本思想及计算 步骤

有限元法的基本思想及计算步骤 有限元法是把要分析的连续体假想地分割成有限个单元所组成的组合体，简称离散化。这些单元仅在顶角处相互联接，称这些联接点为结点。离散化的组合体与真实弹性体的区别在于：组合体中单元与单元之间的联接除了结点之外再无任何关联。但是这种联接要满足变形协调条件，即不能出现裂缝，也不允许发生重叠。显然，单元之间只能通过结点来传递内力。通过结点来传递的内力称为结点力，作用在结点上的荷载称为结点荷载。当连续体受到外力作用发生变形时，组成它的各个单元也将发生变形，因而各个结点要产生不同程度的位移，这种位移称为结点位移。在有限元中，常以结点位移作为基本未知量。并对每个单元根据分块近似的思想，假设一个简单的函数近似地表示单元内位移的分布规律，再利用力学理论中的变分原理或其他方法，建立结点力与位移之间的力学特性关系，得到一组以结点位移为未知量的代数方程，从而求解结点的位移分量。然后利用插值函数确定单元集合体上的场函数。显然，如果单元满足问题的收敛性要求，那么随着缩小单元的尺寸，增加求解区域内单元的数目，解的近似程度将不断改进，近似解最终将收敛于精确解。 用有限元法求解问题的计算步骤比较繁多，其中最主要的计算步骤为： 1）连续体离散化。首先，应根据连续体的形状选择最能完满地描述连续体形状的单元。常见的单元有：杆单元，梁单元，三角形单元，矩形单元，四边形单元，曲边四边形单元，四面体单元，六面体单元以及曲面六面体单元等等。其次，进行单元划分，单元划分完毕后，要将全部单元和结点按一定顺序编号，每个单元所受的荷载均按静力等效原理移植到结点上，并在位移受约束的结点上根据实际情况设置约束条件。 2）单元分析。所谓单元分析，就是建立各个单元的结点位移和结点力之间的关系式。现以三角形单元为例说明单元分析的过程。如图1所示，三角形有三个结点i，j，m。在平面问题中每个结点有两个位移分量u，v和两个结点力分量F x，F y。三个结点共六个结点位移分量可用列

有限单元法与有限元分析

有限单元法与有限元分析 1.有限单元法 在数学中，有限元法（FEM，Finite Element Method）是一种为求解偏微分方程边值问题近似解的数值技术。求解时对整个问题区域进行分解，每个子区域都成为简单的部分，这种简单部分就称作有限元。它通过变分方法，使得误差函数达到最小值并产生稳定解。类比于连接多段微小直线逼近圆的思想，有限元法包含了一切可能的方法，这些方法将许多被称为有限元的小区域上的简单方程联系起来，并用其去估计更大区域上的复杂方程。它将求解域看成是由许多称为有限元的小的互连子域组成，对每一单元假定一个合适的(较简单的）近似解，然后推导求解这个域总的满足条件(如结构的平衡条件），从而得到问题的解。这个解不是准确解，而是近似解，因为实际问题被较简单的问题所代替。由于大多数实际问题难以得到准确解，而有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。 随着电子计算机的发展，有限单元法是迅速发展成一种现代计算方法。它是50年代首先在连续体力学领域--飞机结构静、动态特性分析中应用的一种有效的数值分析方法，随后很快广泛的应用于求解热传导、电磁场、流体力学等连续性问题。 1.1.有限元法分析本质 有限元法分析计算的本质是将物体离散化。即将某个工程结构离散为由各种单元组成的计算模型，这一步称作单元剖分。离散后单元与单元之间利用单元的节点相互连接起来；单元节点的设置、性质、数目等应视问题的性质，描述变形形态的需要和计算精度而定（一般情况单元划分越细则描述变形情况越精确，即越接近实际变形，但计算量越大）。所以有限元中分析的结构已不是原有的物体或结构物，而是同新材料的由众多单元以一定方式连接成的离散物体。这样，用有限元分析计算所获得的结果只是近似的。如果划分单元数目非常多而又合理，则所获得的结果就与实际情况相符合。 1.2.特性分析 1）选择位移模式： 在有限单元法中，选择节点位移作为基本未知量时称为位移法；选择节点力作为基本未知量时称为力法；取一部分节点力和一部分节点位移作为基本未知量时称为混合法。位移法易于实现计算自动化，所以，在有限单元法中位移法应用范围最广。 当采用位移法时，物体或结构物离散化之后，就可把单元总的一些物理量如

有限元单元介绍

第二章单元 在显式动态分析中可以使用下列单元： ·LINK160杆 ·BEAM161梁 ·PLANE162平面 ·SHELL163壳 ·SOLID164实体 ·COMBI165弹簧阻尼 ·MASS166质量 ·LINK167仅拉伸杆 本章将概括介绍各种单元特性,并列出各种单元能够使用的材料类型。 除了PLANE162之外,以上讲述的显式动态单元都是三维的,缺省时为缩减积分(注意：对于质量单元或杆单元缩减积分不是缺省值)缩减积分意味着单元计算过程中积分点数比精确积分所要求的积分点数少。因此,实体单元和壳体单元的缺省算法采用单点积分。当然,这两种单元也可以采用全积分算法。详细信息参见第九章沙漏,也可参见《LS-DYNA Theoretical Manual》。 这些单元采用线性位移函数；不能使用二次位移函数的高阶单元。因此，显式动态单元中不能使用附加形状函数，中节点或P-单元。线位移函数和单积分点的显式动态单元能很好地用于大变形和材料失效等非线性问题。 值得注意的是，显单元不直接和材料性能相联系。例如，SOLID164单元可支持20多种材料模型，其中包括弹性，塑性，橡胶，泡沫模型等。如果没有特别指出的话（参见第六章，接触表面），所有单元所需的最少材料参数为密度，泊松比，弹性模量。参看第七章材料模型，可以得到显式动态分析中所用材料特性的详细资料。也可参看《ANSYS Element Reference》,它对每种单元作了详细的描述，包括单元的输入输出特性。 2.1实体单元和壳单元 2.1.1 SOLID164 SOLID164单元是一种8节点实体单元。缺省时，它应用缩减（单点）积分和粘性沙漏控制以得到较快的单元算法。单点积分的优点是省时，并且适用于大变形的情况下。当然，也可以用多点积分实体单元算法（KEYOPT（1）=2）；关于

有限元2-弹性力学平面问题有限单元法(2.1三角形单元,2.2几个问题的讨论)综述

第2章 弹性力学平面问题有限单元法 2.1 三角形单元(triangular Element) 三角形单元是有限元分析中的常见单元形式之一，它的优点是: ①对边界形状的适应性较好，②单刚形式及其推导比较简单,故首先介绍之。 一、结点位移和结点力列阵 设右图为从某一结构中取出的一典型三角形单元。 在平面应力问题中，单元的每个结点上有沿x 、y 两个方向的力和位移，单元的结点位移列阵规定为： 相应结点力列阵为: (式2-1-1) 二、单元位移函数和形状函数 前已述及，有限单元法是一种近似方法，在单元分析中，首先要求假定（构 造）一组在单元内有定义的位移函数作为近似计算的基础。即以结点位移为已知量，假定一个能表示单元内部（包括边界）任意点位移变化规律的函数。 构造位移函数的方法是：以结点(i,j,m)为定点。以位移(u i ,v i ,…u m v m )为定点上的函数值，利用普通的函数插值法构造出一个单元位移函数。 在平面应力问题中，有u,v 两个方向的位移,若假定单元位移函数是线性的,则可表示成: (,)123 u u x y x y ααα==++ 546(,)v v x y x y ααα==++ (2-1-2)a 式中的6个待定常数α1 ,…, α6 可由已知的6个结点位移分量(3个结点的坐标) {}??? ?? ?????=????? ???? ?????????????=m j i m e d d d d m j j i v u v u v u i {} i i j j m X Y X (2-1-1)Y X Y i e j m m F F F F ?? ?? ???? ???? ??==??????????????????

有限元法基础试题

有限元法基础试题（A ） 一、填空题（5×2分） 1.1单元刚度矩阵e T k B DBd Ω = Ω? 中，矩阵B 为__________，矩阵D 为___________。 1.2边界条件通常有两类。通常发生在位置完全固定不能转动的情况为_______边界，具体指定有限的非零值位移的情况，如支撑的下沉，称为_______边界。 1.3内部微元体上外力总虚功： ()(),,,,e x x xy y bx xy x y y by d W F u F v dxdy δστδτσδ??=+++++??+(),,,,x x y y xy y x u v u u dxdy σδσδτδδ??+++??的表达式中，第一项为____________________的虚功，第二项为____________________的虚功。 1.4弹簧单元的位移函数1N +2N =_________。 1.5 ij k 数学表达式：令j d =_____，k d =_____，k j ≠，则力i ij F k =。 二、判断题（5×2分） 2.1位移函数的假设合理与否将直接影响到有限元分析的计算精度、效率和可靠性。（ ） 2.2变形体虚功原理适用于一切结构（一维杆系、二维板、三位块体）、适用于任何力学行为的材料（线性和非线性），是变形体力学的普遍原理。 （ ） 2.3变形体虚功原理要求力系平衡，要求虚位移协调，是在“平衡、协调”前提下功的恒等关系。 （ ） 2.4常应变三角单元中变形矩阵是x 或y 的函数。 （ ） 2.5 对称单元中变形矩阵是x 或y 的函数。 （ ） 三、简答题（26分） 3.1列举有限元法的优点。（8分） 3.2写出有限单元法的分析过程。（8分） 3.3列出3种普通的有限元单元类型。（6分） 3.4简要阐述变形体虚位移原理。（4分） 四、计算题（54分） 4.1对于下图所示的弹簧组合，单元①的弹簧常数为10000N/m ，单元②的弹簧常数为20000N/m ，单元③的弹簧常数为10000N/m ，确定各节点位移、反力以及单元②的单元力。（10分） 4.2对于如图所示的杆组装，弹性模量E 为10GPa ，杆单元长L 均为2m ，横截面面积A 均为2×10-4m 2，弹簧常数为2000kN/m ，所受荷载如图。采用直接刚度法确定节点位移、作用力和单元②的应力。（10分）

有限元法中的几个基本概念

诚信·公平·开放·共赢 Loyalty Fair Opening Win-win 有限元法中的几个基本概念 有限元法是把要分析的连续体假想地分割成有限个单元所组成的组合体，简称离散化。 这些单元仅在顶角处相互联接，称这些联接点为结点。 离散化的组合体与真实弹性体的区别在于：组合体中单元与单元之间的联接除了结点之外再无任何关联。但是这种联接要满足变形协调条件，即不能出现裂缝，也不允许发生重叠。显然，单元之间只能通过结点来传递内力。 通过结点来传递的内力称为结点力，作用在结点上的荷载称为结点荷载。当连续体受到外力作用发生变形时，组成它的各个单元也将发生变形，因而各个结点要产生不同程度的位移，这种位移称为结点位移。 在有限元中，常以结点位移作为基本未知量。并对每个单元根据分块近似的思想，假设一个简单的函数近似地表示单元内位移的分布规律，再利用力学理论中的变分原理或其他方法，建立结点力与位移之间的力学特性关系，得到一组以结点位移为未知量的代数方程，从而求解结点的位移分量。然后利用插值函数确定单元集合体上的场函数。显然，如果单元满足问题的收敛性要求，那么随着缩小单元的尺寸，增加求解区域内单元的数目，解的近似程度将不断改进，近似解最终将收敛于精确解。 附：FELAC 2.0软件简介 FELAC 2.0采用自定义的有限元语言作为脚本代码语言，它可以使用户以一种类似于数学公式书写和推导的方式，非常自然和简单的表达待解问题的微分方程表达式和算法表达式，并由生成器解释产生完整的并行有限元计算C程序。 FELAC 2.0的目标是通过输入微分方程表达式和算法之后，就可以得到所有有限元计算的程序代码，包含串行程序和并行程序。该系统采用一种语言(有限元语言)和四种技术(对象技术、组件技术、公式库技术生成器技术)开发而成。并且基于FELAC 1.0的用户界面，新版本扩充了工作目录中右键编译功能、命令终端输入功能，并且丰富了文本编辑功能，改善了用户的视觉体验，方便用户快速便捷的对脚本或程序进行编辑、编译与调试。其中并行版在前后处理上进行了相应的改进。

基于Matlab语言的按平面三角形单元划分的结构有限元程序设计模板

基于Matlab语言的按平面三角形单元划分的结构有限元程序设计 专业：建筑与土木工程 班级：建工研12-2 姓名：韩志强 学号： 471220580

基于Matlab语言的按平面三角形单元划分 结构有限元程序设计 一、有限单元发及Matlab语言概述 1. 有限单元法 随着现代工业、生产技术的发展，不断要求设计高质量、高水平的大型、复杂和精密的机械及工程结构。为此目的，人们必须预先通过有效的计算手段，确切的预测即将诞生的机械和工程结构，在未来工作时所发生的应力、应变和位移因此，需要寻求一种简单而又精确的数值分析方法。有限单元法正是适应这种要求而产生和发展起来的一种十分有效的数值计算方法。 有限元法把一个复杂的结构分解成相对简单的“单元”，各单元之间通过结点相互连接。单元内的物理量由单元结点上的物理量按一定的假设内插得到，这样就把一个复杂结构从无限多个自由度简化为有限个单元组成的结构。我们只要分析每个单元的力学特性，然后按照有限元法的规则把这些单元“拼装”成整体，就能够得到整体结构的力学特性。 有限单元法基本步骤如下： (1)结构离散：结构离散就是建立结构的有限元模型，又称为网格划分或单元划分，即将结构离散为由有限个单元组成的有限元模型。在该步骤中，需要根据结构的几何特性、载荷情况等确定单元体内任意一点的位移插值函数。 (2)单元分析：根据弹性力学的几何方程以及物理方程确定单元的刚度矩阵。 (3)整体分析:把各个单元按原来的结构重新连接起来，并在单元刚度矩阵的基础上确定结构的总刚度矩阵，形成如下式所示的整体有限元线性方程： {}[]{}δ F=① K 式中，{}F是载荷矩阵,[]K是整体结构的刚度矩阵,{}δ是节点位移矩阵。 (4)载荷移置:根据静力等效原理，将载荷移置到相应的节点上，形成节点载荷矩阵。 (5)边界条件处理:对式①所示的有限元线性方程进行边界条件处理。 (6)求解线性方程:求解式①所示的有限元线性方程，得到节点的位移。在该步骤中，若有限元模型的节点越多，则线性方程的数量就越多，随之有限元分析的计算量也将越大。 (7)求解单元应力及应变根据求出的节点位移求解单元的应力和应变。

平面三角形单元有限元程序的设计说明

. . P 9 m 9 m 一、题目 如图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 节点编号 单元编号

有限元法基础重点归纳(精)

1、有限元这种数值计算方法起源于20世纪50年代中期航空工程中飞机结构的矩阵分析。 2、有限单元法的基本思想:在力学模型上将一个原来连续的物体离散成为有限个具有一定 大小的单元,这些单元仅在有限个节点上相连接,并在节点上引进等效力以代替实际作用于单元上的外力。 3、节点:网格间相互连接的点。 4、边界:网格与网格的交界线。 5、有限元的优点:①理论基础简明,物理概念清晰,且可在不同的水平上建立起对该法的 理解②具有灵活性和适用性,应用范围极为广泛③该法在具体推导运算中,广泛采用了矩阵方法。 6、有限单元法分类(从选择基本未知量的角度:位移法(以节点位移为基本未知量,通用 性广、力法(以节点力、混合法(一部分以节点位移,另一部分以节点力 7、有限元法分析计算的基本步骤:①结构的离散化②单元分析(选择位移模式,建立单元 刚度方程,计算等效节点力③整体分析④求解方程,得出节点位移⑤由节点位移计算单元的应变与应力。 8、单元划分:将某个机械结构划分为由各种单元组成的计算模型。 9、有限元法基本近似性------几何近似。

10、弹性力学的任务:分析弹性体在受外力作用并处于平衡状态下产生的应力、应变和位移状态及其相互关系等。 11、弹性力学假设所研究的物体是连续的、完全弹性的、均匀的、各向同性的、微小变形的和无初应力的 12、外力:体力(分布在物体体积内的力---重力、惯性力、电磁力面力(分布在物体表面上的力---流体压力、接触力、风力 13、应力:物体受外力作用,或由于温度有所改变,其内部发生的内力。σ={ σx σy σz τx τy τz } = [σx σy σz τx τy τz ]T 14、应变:物体受到外力作用时,其形状发生改变时的形变。---长度和角度。 ε={ εx εy εz γx γy γz } = [εx εy εz γx γy γz ]T 15、位移:弹性体在载荷作用下,不仅会发生形变,还将产生位移,即弹性体位置 的移动。 δ={u v w }=[u v w ]T 16:、变形协调条件:设想在变形前,把弹性体分为许多微小立方单元体。变形后,每个单元体都产生任意变形而变成一些六面体。可能发生这样的情况,这些六面体

平面三角形单元有限元程序设计

平面三角形单元有限元程序设计 P 9 m 9 m 一、题目 如图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)最里层:区分为第 奇/偶 数个计算 X Y P X Y P 节点编号 单元编号

有限元单元法复习资料

1.1有限单元法中“离散”的含义是什么？有限单元法是如何将具 有无限自由度的连续介质的问题转变为有限自由度问题的？位移有限单元法的标准化程式是怎样的？（1）离散的含义即将结构离散化，即用假想的线或面将连续体分割成数目有限的单元，并在其上设定有限个节点；用这些单元组成的单元集合体代替原来的连续体，而场函数的节点值将成为问题的基本未知量。（2）给每个单元选择合适的位移函数或称位移模式来近似地表示单元内位移分布规律，即通过插值以单元节点位移表示单元内任意点的位移。因节点位移个数是有限的，故无限自由度问题被转变成了有限自由度问题。（3）有限元法的标准化程式：结构或区域离散，单元分析，整体分析，数值求解。 1.2单元刚度矩阵和整体刚度矩阵各有哪些性质？各自的物理意义是什么?两者有何区别? 单元刚度矩阵的性质：对称性、奇异性（单元刚度矩阵的行列式为零）。整体刚度矩阵的性质：对称性、奇异性、稀疏性。 单元刚度矩阵Kij物理意义Kij即单元节点位移向量中第j个自由度发生单位位移而其他位移分量为零时，在第i个自由度方向引起的节点力。 整体刚度矩阵K中每一列元素的物理意义是：要迫使结构的某节点位移自由度发生单位位移，而其他节点位移都保持为零的变形状态，在所有个节点上需要施加的节点荷载。 2.1 为了使计算结果能够收敛于精确解，位移函数需要满足什么条件？为什么？ 满足完备性和协调性。 原因：完备性包括两个条件：即刚体位移条件与常应变条件。首先，位移函数必须包含单元的刚体位移。结构中的单元不仅产生与该单元本身变形相应的位移，还可能因其他单元变形而通过节点位移产生单元刚体位移。为了正确反映单元的实际位移形态，位移函数必须具有反映刚体位移的能力。 其次，由于单元位移函数采用多项式，故在单元内部协调条件总能满足，要求反映在相邻单元之间。实质上来说，要求相邻单元间协调是为了保证单元交界面上应变有限。 3.1构造单元形函数有那些基本原则？试采用构造单元几何方法，构造T10单元的形函数，并对其收敛性进行讨论。 答：形函数是定义于单元内坐标的连续函数。通常单元位移函数采用多项式，其中的待定常数由节点位移参数确定，因此其个数应与单元节点自由度数相等。根据实体结构的几何方程，单元的应变是位移的一次导数。为了反映单元刚体位移和常应变即满足完备性要求，位移函数中必须包含常数项和一次项，即完全一次多项式。 3.3 何谓面积坐标？其特点是什么？ 答：三角形单元中，任一点P(x,y)与其3个角点相连形成3个子三角形,其位置可以用下述称为面积坐标的三个比值来确定： L1=A1/A L2=A2/A L3=A3/A 其中A1,A2,A3分别为P23,P31,P12的面积。 各三角形面积为：Ai=1/2* =(ai+bi+ci)/2 由于A1+A2+A3=A,所以有L1+L2+L3=1，Li=(ai+bi+ci)/(2A) 特点：①T3单元的形函数Ni就是面积坐标Li ②面积坐标与三角形在整体坐标系中的位置无关，故称为局部坐标。 ③三个节点的面积坐标分别为节点1（1,0,0），节点2（0,1,0），节点3（0,0,1），形心的面积坐标（1/3，1/3 ，1/3）。④单元边界方程为Li=0 (i=1,2,3)； ⑤在平行于2,3边的一条直线上，所有点都要相同的面积坐标。⑥面积坐标与直角坐标互为线性关系。 体积坐标：P点与四面体四个面围成的四个子四面体的体积与原来四面体体积的比值。即 剪切闭锁现象：当梁的高度与梁的长度之比t/l趋于零时，这种单元将出现这种现象，算得的挠度趋于零。 为克服剪切闭锁，使C0型单元适用于各种高度的梁。采用减缩积分方案与假设减应变法。 零能模式：对应于某种非刚体位移模式，减缩积分时高斯点上的应变正好等于零，此时的应变能当然也为零，这种非刚体位移模式称为零能模式。采用减缩积分时会发生零能模式。 5.1、等参单元：将整体坐标系中xy中形状中较复杂的真实单元变换成局部坐标系xy中规则的标准单元，然后在标准单元中构造形函数。由于坐标变换式与单元位移函数中用了相同的形函数N i（ξ，η），故称这种变换为等参变换，相应的单元称为等参单元。 2、等参单元的优越性：①有些工程较复杂，用直边单元离散这些结构需要大量的单元才能得到较好的近似，而曲边的等参单元可方便地离散复杂结构。②如在单元内多取些节点，单元便具有较多的位移自由度，从而就能够插值表示较复杂的单元内部位移场，这样就提高了单元本身的精度。③等参单元刚度矩阵、荷载矩阵的计算是在规则单元域内进行的因此不管被积函数多么复杂，都可以方便地采用标准化数值积分。 3、数值积分的阶次：对于N点积分，当被积函数为m次多项式且m<=2N-1时，可得精确积分值。反之，对于m次多项式的被积函数，精确积分要求的积分点数N>=（m+1）/2。 6.1、工程梁和剪切梁的基本假设？有两种梁弯曲理论①工程梁理论基本假设：平截面假设与横向纤维无挤压假设。前者认为梁横截面变形后仍为平面，且垂直于变形后的中性轴。该假设意味着横向剪切应变γxy =0，后者认为梁的横向纤维无挤压，即εy=0。②剪切梁理论基本假设：横向纤维元无挤压与另一假设认为法平面变形后仍为平面，但不再垂直于变形后的中性轴。 6.2. 剪切梁怎么考虑剪切影响：在结构单元分析中，可在工程梁单元的基础上考虑剪切变形的影响，也可通过挠度与转角各自独立插值直接构造剪切梁单元。 6.3对于杆系结构单元，为什么要在局部坐标系内建立单元刚度矩阵？为什么还要坐标变换？（1）在局部坐标系内可以更方便的建立单 元刚度矩阵。（2）在整体分析中，对所有单元都应采用同一个坐标系即整体坐标系X Y，否则围绕同一节点的不同单元对节点施加的节点力不能直接相加。因此，在进行整体分析之前，还需要进行坐标转换工作，把局部坐标系中得出的单元刚度方程转换成整体坐标系中的单元刚度方程，从而得出整体坐标系中的单元刚度矩阵。 7.1. 薄板弯曲理论基本假定：第一条：板厚方向的挤压变形可忽略不计，即εz=0,。这项假设类似于梁的横向纤维间无挤压假设。第二条：在板弯曲变形中，中面法线保持为直线且仍为弹性曲面（挠度曲面）的法线。第三条：薄板中面只发生弯曲变形，没有面内的伸缩变形，即中面水平位移（u）z=0=(v) z=0=0. 7.2. 厚板理论基本假设：板的中面法线变形后基本保持为直线，但因横向变形的缘故，该直线不在垂直于变形后的中面。因此，法线绕坐标轴的转角θx、θy不再是挠度的导数，而是独立变量；中面内的线位移和板厚方向的挤压变形也可忽略。 7.3. 薄板、厚板基本假定的不同：薄板：板弯曲变形中，中面法线保持为直线且仍为弹性曲面法线。厚板：板中面法线变形后仍基本保持为直线，但该直线不再垂直于变形后的中面。 7.4. DKT单元：离散Kirchhoff理论的基本思想是在若干离散点上满足Kirchhoff直法线假设。基于这种理论构造薄板单元时，w,θx,,θy 也各自独立插值；然后在若干离散点上引入直法线假设。这样构造的单元叫做DKT单元 8.1. 薄壳单元基本假设：薄壳理论假设：薄壳发生微笑变形时，忽略沿壳体厚度方向的挤压变形；且认为直法线假设成立，即变形后中面法线保持为直线且仍为中面的法线；壳体变形时中面不但发生弯曲，而且面内也将产生面内伸缩变形；折板假设；非耦合假设。 薄壳与薄板理论的假设的异同点：相同点：直法线假设和法向（板厚度方向）的纤维无挤压假设均成立。不同点：薄板中面只发生弯曲变形，没有面内的伸缩变形，即中面水平位移为零，而壳体变形时中面不但发生弯曲，而且也将产生面内伸缩变形。 厚壳分析的假设：变形前后的中曲面法线变形后仍基本保持为直线，但因横向剪切变形的缘故，该直线不再垂直于变形后的中曲面，此外，壳体厚度方向的挤压变形可以忽略。 与厚板理论的假设的 相同点：中面法线变形后仍基本保持为直线，但因横向剪切变形的缘故，该直线不在垂直于变形后的中面。厚度方向的挤压变形忽略不计。不同点：厚板理论的假设中，中面内的线位移可以忽略，而厚壳理论的假设中，中面内的位移不可忽略，并且厚壳的位移场可用中面位移表示。 8.2. 平板型单元：组成的折板系统去代替原来的壳体，由平面应力状态与平板弯曲应力状态加以组合而得壳体的应力单元。 分析这种单元时所提出的假设：理论假设：薄壳发生微笑变形时，忽略沿壳体厚度方向的挤压变形，且认为直法线假设成立，即变形后中面法线保持为直线且仍为中面的法线。，折板假设，非耦合假设。 应用平板型壳单元可能会出现的问题，如何解决：1.单元共面问题，解法：引入唯一边界条件可解方程Ka=P 。2.虚拟旋转刚度，解法：在特殊节点上给以任意的虚拟刚度系数。Kθzθzθzi=0，经坐标变换，整体坐标系中该节点平衡方程将满足唯一解条件。赋予Kθzθz任何值。3.新型平面膜元：在平面膜元角点上增加旋转自由度θz，使其有对应的刚度。 8.3. 面内变形与弯曲变形之间非耦合的假设是针对什么提出来的？试说明单元组装时，面内效应与弯曲效应的耦合将会出现。 答：面内变形与弯曲变形之间非耦合的假设是针对局部坐标系下的单元提出的。 9.1. 减少自由度的措施有哪些？各自基本概念如何？ 答：1.恰当利用结构对称性。基本思想：利用结构的对称性，取结构一部分建立有限元模型。根据荷载对称性，分析对称面上的位移状态，以确定对称面上节点的位移边界条件。2.采用子结构技术。基本思想：在大型复杂结构的有限元分析中，可将原结构分成若干区域，每个区域作为一个子结构，这些子结构在其公共边界上互相连接起来。 2. 为什么说位移法中应力解的精度低于位移解？ 答：在位移有限单元法中，沿单元边界是连续的，而位移的导数通常不连续，因此，在单元边界上应力是不连续的；基本未知量是位移，而单元应变和应力是由位移求导得到的，因此应力精度低于位移精度。 3. 在无法获得精确解的条件下，如何进行误差估计？ 答：有限元解法的误差估计有：残值法，后处理法。后处理法：由于无法获得精确解，一般以修匀后的改进值σ*作为“精确解”进行误差估计，通过与精确值误差范数对比，这样做非常有效。

有限元d 分析与介绍

有限元分析（FEA，Finite Element Analysis）利用数学近似的方法对真实物理系统（几何和载荷工况）进行模拟。还利用简单而又相互作用的元素，即单元，就可以用有限数量的未知量去逼近无限未知量的真实系统。 有限元分析是用较简单的问题代替复杂问题后再求解。它将求解域看成是由许多称为有限元的小的互连子域组成，对每一单元假定一个合适的(较简单的）近似解，然后推导求解这个域总的满足条件(如结构的平衡条件），从而得到问题的解。这个解不是准确解，而是近似解，因为实际问题被较简单的问题所代替。由于大多数实际问题难以得到准确解，而有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。 有限元是那些集合在一起能够表示实际连续域的离散单元。有限元的概念早在几个世纪前就已产生并得到了应用，例如用多边形（有限个直线单元）逼近圆来求得圆的周长，但作为一种方法而被提出，则是最近的事。有限元法最初被称为矩阵近似方法，应用于航空器的结构强度计算，并由于其方便性、实用性和有效性而引起从事力学研究的科学家的浓厚兴趣。经过短短数十年的努力，随着计算机技术的快速发展和普及，有限元方法迅速从结构工程强度分析计算扩展到几乎所有的科学技术领域，成为一种丰富多彩、应用广泛并且实用高效的数值分析方法。 有限元方法与其他求解边值问题近似方法的根本区别在于它的近似性仅限于相对小的子域中。20世纪60年代初首次提出结构力学计算有限元概念的克拉夫（Clough）教授形象地将其描绘为：“有限元法=Rayleigh Ritz法+分片函数”，即有限元法是Rayleigh Ritz法的一种局部化情况。不同于求解（往往是困难的）满足整个定义域边界条件的允许函数的Rayleigh Ritz法，有限元法将函数定义在简单几何形状（如二维问题中的三角形或任意四边形）的单元域上（分片函数），且不考虑整个定义域的复杂边界条件，这是有限元法优于其他近似方法的原因之一。

平面三角形单元有限元程序设计

平面三角形单元有限元程序设计

P 9 m 9 m 一、题目 如图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行，并建立坐标系，则 X Y P X Y P

边界约束的处理：划0置1法 适用：这种方法适用于边界节点位移分量为已知(含为0)的各种约束。 做法： （1）将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1，其他元素改为零；同 时将载荷列阵｛R｝中相应元素用已知位移置换。 ◎这样，由该方程求得的此位移值一定等于已知量。 （2）将〔K〕中已知位移相应的列的非主对角成元素也置0，以保持〔K〕的对称性。 ◎当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为 保证方程成立，须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。◎若约束为零位移约束时，此步则可省去。 特点： （1）经以上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。 （2）但这种方法不改变方程阶数，利于存贮。 （3）不过，若是要求出约束反力，仍要重新计算各个划去的总刚元素。 程序如下： 变量说明 NNODE 单元节点数 NPION 总结点数 NELEM 单元数 NVFIX 受约束边界点数 FIXED 约束信息数组 NFORCE 节点力数 FORCE 节点力数组

有限元法与有限差分法的主要区别

有限元法与有限差分法的主要区别 有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑，可分为中心格式和逆风格式。考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。在河道数值模拟中，常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等。根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法，从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形网格，从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。不同的组合同样构成不同的有限元计算格式。对于权函数，伽辽金(Galerkin)法是将权函数取为逼近函数中的基函数；最小二乘法是令权函数等于余量本身，而内积的极小值则为对代求系数的平方误差最小；在配置法中，先在计算域内选取N个配置点。令近似解在选定的N个配置点上严格满足微分方程，即在配置点上令方程余量为0。插值函数一般由不同次幂的多项式组成，但也有采用三角函数或指数函数组成的乘积表示，但最常用的多项式插值函数。有限元插值函数分为两大类，一类只要求插值多项式本身在插值点取已知值，称为拉格朗日(Lagrange)多项式插值；另一种不仅要求插值多项式本身，还要求它的导数值在插值点取已知值，称为哈密特(Hermite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标，有对称和不对称等。常采用的无因次坐标是一种局部坐标系，它的定义取决于单元的几何形状，一维看作长度比，二维看作面积比，三维看作体积比。在二维有限元中，三角形单元应用的最早，近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元，常采用的插值函数为有La grange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。对于有限元方法，其基本思路和解题步骤可归纳为(1)建立积分方程，根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。(2)区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。(3)确定单元基函数，根据单元中节点数目及对近似解精度的要求，选择满足一定插

三角形单元有限元程序设计

三角形单元有限元程序设计 %******************************************************************* % NNODE 单元节点数NPION 总结点数 % NELEM 单元数NVFIX 受约束边界点数 % FIXED 约束信息数组NFORCE 节点力数 % FORCE 节点力数组COORD 结构节点坐标数组 % LNODS 单元定义数组 % YOUNG 弹性模量POISS 泊松比 % THICK 厚度 B 单元应变矩阵(3*6) % D 单元弹性矩阵(3*3) S 单元应力矩阵(3*6) % A 单元面积ESTIF 单元刚度矩阵 % ASTIF 总体刚度矩阵ASLOD 总体荷载向量 % ASDISP 节点位移向量ELEDISP 单元节点位移向量 % STRESS 单元应力FP1 数据文件指针 %******************************************************************* format short e %设定输出类型 clear %清除内存变量 FP1=fopen('1-1.txt','rt'); %打开输入数据文件存放初始数据 %读入控制数据 NELEM=fscanf(FP1,'%d',1), %单元个数（单元编码总数） NPION=fscanf(FP1,'%d',1), %结点个数（结点编码总数） NVFIX=fscanf(FP1,'%d',1) %受约束边界点数 NFORCE=fscanf(FP1,'%d',1), %结点荷载个数 YOUNG=fscanf(FP1,'%e',1), %弹性模量 POISS=fscanf(FP1,'%f',1), %泊松比 THICK=fscanf(FP1,'%f',1) %厚度 LNODS=fscanf(FP1,'%d',[3,NELEM])' %单元定义数组（单元结点号） %相应为单元结点号（编码）、按逆时针顺序输入 COORD=fscanf(FP1,'%f',[2,NPION])' %结点坐标数组 %坐标：x，y坐标（共NPOIN组） FORCE=fscanf(FP1,'%f',[3,NFORCE])' %结点力数组 %(n,3) n:受力结点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向 FIXED=fscanf(FP1,'%d',[3,NVFIX])' %约束信息数组 % (n,1):约束点(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0 %******************************************************************* %生成单元刚度矩阵并组成总体刚度矩阵 ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0 %*******************************************************************