计算结构力学课程讲义.doc

第1章 绪 论
1.1 课程内容
(1) 研究内容
本课程主要研究工程结构计算机分析（数值分析）的常用方法——有限单元法、加权残数（余量）法和边界单元法的基本概念、基本原理及其应用。
(2) 参考书籍
课程的主要参考书籍如下：
唐锦春，孙炳楠，郭鼎康，计算结构力学，浙江大学出版社，1989
丁皓江, 谢贻权, 何福保，弹性和塑性力学中的有限单元法，机械工业出版社，1989
王勖成，有限单元法，清华大学出版社，2003
王勖成，邵敏，有限单元法基本原理与数值方法，第二版，清华大学出版社，1997
徐次达，固体力学加权残数法，同济大学出版社，1987
孙炳楠，项玉寅，张永元，工程中边界单元法及其应用，浙江大学出版社，1991
Bath, K. J. Finite Element Procedures, Prentice-Hall, Inc., 1996.
Zienkiewicz, O. C., The Finite Element Method, 5th Edition, McGraw Hill, 2001.
Brebbia, C.A., The Boundary Element Method for Engineers, Pentech Press, London, 1978.
Chandrupatla, T. R., Belegundu, A.D. Introduction to Finite Elements in Engineering, Prentice-Hall, Inc., 2002.

1.2 结构分析方法概述
一个工程技术问题总可由一组基本方程（通常是微分方程）加一组边界条件描述，即由下式给出：
基本方程：L(u)-p=0， ÎV（域内）
边界条件：B(u)-g=0， ÎS（边界）
式中L、B为算子，p、g为已知函数。
工程技术问题的常用分析方法有：
(1) 解析方法
只适用于少数简单问题，即形状规则且外部作用（如外荷载）简单的结构分析问题。
(2) 数值方法
数值方法可分为区域型方法和边界型方法。常用的区域型方法包括有限差分法、加权残数法、里兹(Ritz)法（变分法）和有限单元法等，其中有限差分法是直接对基本微分方程进行离散，再对离散后的代数方程进行求解；后几种方法则是先建立基本方程（一般是微分方程）的等效积分表达式，再进行离散求解。边界型方法中最典型的是边界单元法。它是先将基本微分方程变换为等效的边界积分方程，再在边界上对其进行离散求解。
例如，图1.1给出了一个受复杂横向荷载（分布荷载、集中力、集中力偶等）作用的两端固定变截面梁。为求梁的挠度和内力，可列出梁的基本方程和边界条件如下：
l
q(x)
y
x
h(x), EI(x)

P
m
A
B
b
h(x)

图1.1 变截面单跨梁受横向荷载作用

基本方程：L(u)-p=0，ÎV（域内）
——EI(x)y ’’= -M(x), 0£x£l.
边界条件：B(u)-g=0， ÎS（边界）
——(y)x=0或x=l=0，(y ’)x=0或x=l=0
以下分别就采用加权残数法、里兹法（位移变分法）和有限单元法的基本原理进行讨论。
(1) 加权残数法
为求近似解，设试探函数代入基本方程和边界条件，得残值：
RL=L(u)-p（域内），RB=B(u)-g（边界）
迫使残值在某种平均意义（加权积分）上等于零，则有

由此可得到关于待定系数ai的代数方程组，解方程可求得待定系数及解答的近似表达式，其中的试函数可以选择多项式、三角函数、样条函数等。

(2) 里兹法（位移变分法）
里兹法的理论依据是最小势能原理。该原理可表述为：给定外力作用下，满足几何条件的各种可能位移中，真实的位移使总势能取极值，据此有
d(U+UR)=0
假设满足位移边界条件的位移函数为：

将其代入方程得到关于待定系数Ai的代数方程，解方程可得Ai。
里兹法需要在整个计算区域上假设近似函数，很难适应形状（边界）较复杂或解答较难预测的问题。

(3) 有限单元法
有限单元法的理论依据是最小势能原理或其他形式的变分原理。该方法与里兹法的主要区别是不在整体计算区域上假设近似函数，而是先将连续的求解区域离散为一个由有限个单元组成并按一定方式相互连接的单元集合体，再以各单元连接结点处的场量（如位移量）作为基本未知量，在各单元内假设近似函数（通过结点未知量插值得到），从而将一个无限自由度问题简化为有限自由度问题。
u
分段假设试函数
x

图1.2 一维试函数的分段假设

例如图1.2中的曲线是某个一维问题的目标函数曲线，若采用里兹法对整个区段假设一个近似的试函数，显然比较困难。但如果现对整个区域进行分段（如图中短线为分段线），再对各个区段假设试函数，则要简单和准确得多，如可将各区段均假设为二次函数。哟次可见，有限单元法可视为一种分片（或分块、分段）形式的变分法。
虽然有限单元法的理论依据和里兹法是一致的，但采用了分片（或分块、分段）假设试函数的处理方法以后，使得该方法的具体实施变得简便易行，具有了优越的可操作性和更为明确的物理意义，也使得该方法具有了其他方法（如里兹法）所不具备的优点：
1) 概念简单、明确，易为工程人员接受；也可建立严格的数学分析和证明；
2) 适用性十分广泛，适应于各类复杂边界和不同外部作用的问题；
3) 求解过程程序化，易于编程和计算机实现。

1.3 课时安排
课程的总体课时安排如下：
有限单元法部分包括概论、进展；平面三角形、矩形、等参元；杆元、板元等，共约20个课时；加权残数法部分包括基本原理、方法分类，以及伽辽金（Galerkin）方法、最小二乘法的应用，共需约4～6个课时；边界单元法主要包括基本原理（以二维势问题为例）；梁弯曲和板弯曲问题，共需约4～6个课时。

思考题
1.1 区域型分析法和边界型分析法在对问题的基本方程和边界条件的处理上有何不同和相同点？试分别举例说明。
1.2 里兹法和有限单元法的理论依据、基本未知量的选取、试函数的假设等方面有何异同点？
1.3 与里兹法相比，有限单元法在解决复杂问题上的适应性更为广泛，你认为主要的原因在于那些方面？

第2章 有限单元法
2.1 概述
2.1.1 发展概况
有限单元法的发展概况：
1943年R. Courant尝试应用三角形区域上定义的分片连续函数和最小势能原理解决St. Venant扭转问题，是较早的有限元思想的体现：
R. Courant, Variational Methods for the Solution of Problems of Equilibrium and Vibrations, Bulletin of the American Mathematical Society, 49: 1-23, 1943
1956年M.J. Turner，R.W. Clough等将刚架矩阵位移法推广到弹性力学平面问题，开始了有限元的第一个成功尝试和应用；用直接刚度法建立单元刚度特性：
M.J. Turner, R.W. Clough, H.C. Martin and L.T. Topp, Stiffness and deflection analysis of complex structures, J. Aeronaut. Sci., 25: 805-823, 1956
1960年Clough第一次提出“有限单元法（FEM）”的名称，沿用至今。
Zienkiewicz等——编写第一本有限元方面专著：
O.C. Zienkiewicz and Y.K. Cheung, The Finite Element Method in Continuum and Structural Mechanics, McGraw-Hill, New York, 1965
1963-1964年发现该方法是基于变分原理的里兹法的另一种形式，确立其理论基础。我国冯康在同一时期独立提出并证明了该方法：
Melosh证明了位移法就是基于最小势能原理的Rayleigh-Ritz法
冯康，基于变分原理的差分格式，应用数学和计算数学，1965，2(4): 238-262
1960至今：实际工程应用：平面Þ空间Þ板壳；静力Þ动力、波动Þ稳定；弹性Þ塑性Þ粘弹性、复合材料；固体Þ流体、传热等连续介质力学；计算分析Þ优化设计、与CAD技术结合。
E.L. Wilson：编写了第一个公开的有限元软件SAP；
通用有限元软件：SAP、ADINA、NASTRAN、ANSYS、ABAQUS、MIDAS等
从半个多世纪以来有限单元法的萌芽、理论依据的证明和充实及其逐步的广泛应用可以看到，它的发展和计算机软硬件的发展基本上是同步的。如果没有计算机的强大软硬件支撑，有限单元法只有其微不足道的一点理论上的意义，而没有更为重要的实际应用的意义。

2.1.2 有限单元法概念
(1) 离散化
离散化的过程是将连续体划分为有限数目、有限大小的单元的集合体。单元与单元之间只在指定点（即结点）连接，其他位置则一般保持连续即可。单元可以具有不同的形状，即单元外形可以不同；单元与单元之间可以有不同的连接方式，即单元的结点数目、位置可以不同。

连续体
rg

图2.1 连续体离散为单元集合体示例

(2) 单元分析
对典型单元假设位移模式（由各结点位移插值），再分析单元的力学特性，建立单元的结点力与结点位移之间的关系，即单元刚度方程：
{F}e=[k]{D}e）
并将各类荷载变换为作用在结点上的等效结点荷载。
(3) 整体分析
将各单元刚度方程集成整体结构的整体刚度方程：
{F}=[K]{D}
根据结点的平衡条件，得最终的有限元方程：
[K]{D}={R}
求解该方程可得到未知的结点位移。
(4) 再次单元分析
求出各单元的应变和应力。

2.2 弹性力学平面问题的矩阵描述
2.2.1两类平面问题
弹性力学的平面问题可分为平面应力问题和平面应变问题两类。实际上，所有的弹性力学问题都是空间问题。所谓平面问题，并不是说这个问题所分析的对象本身（如形状、荷载分布）是平面的，而是指该问题的形状、外部作用以及问题的解答（即由此产生的效应，如位移、应力等）只在平面内有变化，而沿着平面外就保持不变了。因此可以肯定地说，所谓的平面问题就是一个特殊的空间问题。
那么，是不是一个问题的形状和外部作用（即已知的位移和应力边界条件）只在平面内发生变化，而沿着平面外保持不变了，这个问题就是平面问题呢？不是的，还必须附加其他条件，这一结论才能成立。这个附加条件就是该问题沿平面外的尺寸与平面内尺寸相比要么非常小（如无限短），要么非常大（如无限长）。如果符合前者条件，则弹性体内只存在平面内的应力，而平面外的应力均为零，故这类问题称为平面应力问题；如果符合后者条件，则弹性体内只存在平面内的位移或平面内的应变，而平面外的位移及应变均等于零，故这类问题称为平面应变问题。

z
x
t/2
O
y
t/2
y
体积力
表面力

图2.2 平面应力问题示例
x
O
y
z

图2.3 平面应变问题示例

2.2.2 基本量和基本方程的矩阵表示
(1) 基本量
应力：{s}=[sx sy txy]T，
应变：{e}=[ex ey gxy]T，
位移：{f}=[u v]T；
体积力：{p}=[px py]T，
表面力：{q}=[qx qy]T。

(2) 基本方程
几何方程：
物理方程：{s}=[D] {e}，
式中[D]：弹性矩阵。对平面应力问题，

平衡方程的弱形式——能量原理
1) 虚功原理：弹性体在外力作用下处于平衡状态的充分和必要条件是，外力在满足变形协调推进的虚位移上所做的虚功等于其内力在相应的虚应变上所做的虚功，即

2) 最小势能原理：满足几何条件的各种可能位移中，真实位移使P取驻值，即
dP=d (U+UR)=0

(3) 平面问题的常用单元：

三结点三角形单元 六结点三角形单元 四结点矩形单元

四结点任意四边形单元 八结点曲边四边形单元
图2.4 平面问题的常用单元

单元的加密方法可分为p型加密和h型加密。前者是保持单元的大小和形状不变，而提高单元的插值函数的阶数；后者是指采用同一类单元，但加密单元的网格，即减小单元的尺寸，增加离散单元的数目。

2.3 三结点三角形单元分析
2.3.1 离散
将一平面结构离散为ne个单元，设结点综述为n个，则结构的基本未知量取为这n个结点的位移，即
{D}=[u1 v1 u2 v2 …un vn]T

y
ui (Ui)
i
x
vi (Vi)
uj (Uj)
j
vj (Vj)
um (Um)
m
vm (Vm)

图2.5 三结点三角形单元

2.3.2 单元位移模式
单元结点力向量：
{F}e=[Ui Vi Uj Vj Um Vm] T
单元结点位移向量：
{D}e=[ui vi uj vj um vm]T

(1) 位移模式的建立
所谓位移函数或称位移模式，是指利用单元的结点位移将单元内任一点的位移用坐标的函数表示出来。
对于三结点三角形单元，假设其位移模式是坐标的线性函数，则有
u=a1+a2 x+a3 y
v=a4+a5 x+a6y
系数a1~ a6由6个结点位移分量（ui、vi、uj、vj、um、vm）确定。
将位移模式写成结点位移的显式：
u= Niui+ Nj uj+ Nmum
v= Nivi+ Nj vj+ Nmvm
Ni、Nj、Nm称为形状函数（形函数）或插值函数，

其中ai、bi、ci是分母行列式第1行各元素的代数余子式，展开后可表示为
ai=xjym- xmyj，bi=yj-ym，ci= -xj+xm（i, j, m轮换）

j
y
i
m
x
P
j
1
i
m
Ni(x,y)

图2.6 形函数的物理意义

(2) 形函数性质
(a) (Ni)i=1，(Ni)j=0，(Ni)m=0；
(b) 单元内任一点：Ni+Nj+Nm=1。

(3) 位移模式的矩阵表示

其中[N]为形函数矩形，且
。

2.3.3 单元应变矩阵和应力矩阵
将位移模式{f}=[N]{D}e代入几何方程，得单元应变为：
{e}=[B]{D}e
[B]称为单元应变矩阵，[B]= [Bi Bj Bm]，而各子块为。
可见[B]为一常量矩阵，故三结点三角形单元为一常应变单元。
将{e}=[B]{D}e代入物理方程，可得
{s}=[S]{D}e=[D][B]{D}e
[S]称为单元应力矩阵，[S]= [S1 S2 S3]，对平面应力问题
。
可见[S]也是一常量矩阵，故三结点三角形单元为一常应力单元。

2.3.4有限元方程的建立
将连续体近似看作为由一系列只在结点相连的单元组装而成的集合体，并且对各个单元均建立了其位移模式，那么该集合体的最终位移法方程，即最终的有限元方程可由变分原理（此时为最小势能原理）建立。
位移模式确定后，离散结构的可能位移就是由不同结点位移决定的分片函数。各可能位移函数中，真实的位移函数使离散结构的总势能取驻值（处于稳定平衡状态时为最小值），即
dP=d (U+UR)=0

令：
{D}e=[G]{D}

这里[G]是一6´2n的位置矩阵，表示单元结点位移{D}e在整体结点位移向量{D}中的位置；[k]e为单元刚度矩阵；[R]e单元等效结点荷载向量。于是有

令：

这里[K]为整体刚度矩阵；[R]为整体等效结点荷载向量，则有：

由势能驻值原理dP=0可得
或

2.3.5 单元刚度矩阵
(1) 单刚列式

其中，

(2) 单刚性质
单元刚度元素kpq的物理意义为：第q个结点位移分量为单位位移（其它结点位移等于0），所引起的第p个结点力分量。
例如要获得元素k26的直观解释，可先将单元的所有6个结点位移约束住，即在每个结点处沿两个方向同时施加链杆约束，如图2.7所示。
y
m
x
k26
i
j
1
图2.7 单元刚度元素k26的物理意义
单元刚度矩阵[k]的性质：
1) 对称性： kpq = kqp；
2) 奇异性：因具有刚体位移；
3) 每行（列）元素之和为零；
例如，第6列元素的意义：当第六个结点位移等于1（D6=1），其他结点位移等于0时，在各结点位移方向施加的结点力的大小，即图2.7中6个附加链杆约束上的约束力。因单元在这些结点力作用下处于平衡，故有：
k16+ k36+ k56=0；
k26+ k46+ k66=0。
4) 主对角元素恒大于零：kqp>0；
5) [k]取决于单元的形状、方位和弹性常数，与所在位置（即平移或np 转动）无关。

(3) 推导单刚的另一种方法——由单元的平衡条件导出
物理方程
几何方程
位移模式
虚功方程
{f }=[N]{d}e
{e}=[B]{d}e
{s}=[S]{d}e ，[S]= [D][B]
{F}e=[k]{d }e，[k]= [B]T [D] [B]tA
结点位移 位移 应变 应力 结点力
{d}e {f } {e} {s} {F}e

图2.8 单元分析中各物理量的联系图

假设单元发生了虚位移，则结点力在结点虚位移上做的虚功=单元应力在虚应变上的虚功（虚变形能）：

令：

又因{dD}e是任意的，故有
{F}e=[k]e{D}e

2.3.6 单元等效结点荷载
在将连续体用近似的单元集合体代替以后，单元集合体中的每个结点可认为一方面受到外荷载的作用，另一方面受到连接该结点的各个单元的对其产生的约束力（即结点力）的作用。各结点在这两组力的作用下总是处于平衡的。但是，如果外荷载不是直接作用在结点上，而是作用在单元内部，那么怎么建立结点的平衡关系呢？一种有效方法是先建立单元的等效结点荷载，再对结点建立平衡关系。
y
Yi
i
x
Xi
Yj
j
Xj
Ym
m
Xm
Px
Py
建立单元的等效结点荷载，或称将非结点荷载需等效移置到结点上，该等效移置可利用虚功原理进行。

图2.9 单元的等效结点荷载分量

2.3.7 整体刚度矩阵
(1) 建立有限元方程的方法：
1) 由最小势能原理建立：
[K]=å([G]e)T[k] e [G]e
2) 由结点平衡条件建立

(2) 整体刚度矩阵[K]的集成方法：对号入座。
[K]=å([G]e)T[k] e [G]e，
其中([G]e)T决定[kpq]e在[K]中的行，[G]e决定在[K]中的列位置。

(3) 整体刚度矩阵的性质
[K]元素的物理意义
Kpq：结构第q个结点位移分量为单位位移（Dq=1，其它结点位移=0），在第p个结点位移方向所施加的结点力的大小。
如 K25为结构第5个结点位移分量为单位位移，即D5=1，其它结点位移均为零时，在第2个结点位移方向所需施加的结点力的大小。
[K]的性质：
1) 对称性： Kpq = Kqp；主对角元恒大于零；
2) 稀疏矩阵，且一般呈带状分布；
例如，平面问题最大半带宽=2´(单元结点号之差最大值+1)。
3) 引入位移边界条件前为奇异矩阵，引入后为正定矩阵。

2.3.8 位移边界条件的引入
以上集成的总体整体刚度矩阵并不包含位移边界条件的信息，这种刚度矩阵常称为原始刚度矩阵，它是一个奇异矩阵。要使得最终的有限元方程可解，必须引入唯一边界条件，排除结构的刚体位移。引入位移边界条件的常用方法有以下几种：
(1) 直接引入法（矩阵缩小法）
将[K]{D}={R}改写为如下子块形式：

其中{Da}、{Db}分别为未知结点位移向量和已知结点位移向量，[Kaa]、[Kab]、[Kba]、[Kbb]、{Ra}、{Rb}分别为与之对应的刚度矩阵和荷载向量子块。
将上述方程的上半部分（对应已知位移）展开，得
[Kaa] {Da}={Ra}-[Kab] {Db}
由该方程可解出未知结点位移。该方法将改变原始刚度矩阵的阶数及元素的顺序，不利于程序实现，因此在计算机编程中一般很少采用。

(2) 对角元素改1法（零位移边界）
设结构的总自由度数（即结点位移总数）为N，若第i个结点位移为零（即Di=0），则将[K]中对角元素Kii改为1，而第i行和第i列的其他元素改为0，荷载向量[R]中的第i个元素也改为0。其实质是将原有限方程中的第i个方程用方程Di=0代替，而其他方程中与Di对应的系数也改为0，表明Di对其他方程没有影响，同时保证了修改后的刚度矩阵仍具有对称性。

(3) 乘大数法
若第i个结点位移已知，即，则将整体刚度矩阵[K]中的对角元素Kii改为a Kii，其中a为一大数（如1020），而荷载向量[R]中的第i个元素Ri改为，原方程成为以下形式

其实质是将原有限方程中的第i个方程用的以下近似方程代替：

将上述方程各项同除以大数a，除第i项及方程右端项外，其他各项均趋于0，故等价于。

2.3.9 位移模式与解答的收敛准则
(1) 有限元解答的收敛准则
为使有限元解答能够收敛于精确解，单元位移模式应满足以下条件：
1) 位移模式必须包含单元的刚体位移；
对弹性力学平面问题，其刚体位移表达式为：
u=u0-wy，v=v0+wx
因此，平面问题的位移模式必要包含上述刚体位移表达式中的各项，才能保证最终解答的收敛性。
2) 位移模式必须包含单元的常量应变；
条件(1)、(2)合并起来可称为完备性要求。对平面问题来说，就是要求位移模式必须包含常数项和完整的一次项。完备性条件是解答收敛的必要条件。
3) 位移模式应尽可能保证位移的连续性。
该条件实际上就是协调性条件，但一般情况下并不是一个必要性条件。
如果位移模式同时满足上述完备性和协调性条件，那么就组成了解答收敛的充分条件。
对一般的弹性力学平面或空间问题，只需要求单元内部以及相邻单元的公共边界上的位移本身连续，即位移模式具有C0连续性。对有些问题，可以放松对协调性的要求，只要通过分片试验，那么也能保证解答的收敛性。
三结点三角形单元的位移模式既满足完备性，又满足协调性的要求（在单元边界上呈线性分布，可由两个结点位移唯一确定），是一种协调单元。证明如下：
单元内：单值连续；
y
x
(1)
i
j
m
p
(2)
相邻单元之间：uij(1)= uij(2)？ vij(1)= vij(2)？
ij边的方程：y=ax+b，则
uij=a1+ a2x+a3(ax+b)= cx+d
uij(1)、uij(2)均为坐标的线性函数，
故可由i、j两点的结点位移唯一确定。

图2.10 两相邻单元

(2) 多项式位移模式的一般选择规则
位移模式应与单元局部坐标的选取无关，即满足所谓的几何各向同性。对于平面问题，位移模式中的x和y的各阶项应保持对称，即有了xmyn项，则应同时具有xnym项。为保证位移模式关于x和y坐标的对称性，通常从以下的Pascal三角形中选取多项式位移模式的各项：
1
x y
x2 xy y2
x3 x2y x y2 y3
x4 x3y x2y2 xy3 y4
…………
图2.11 Pascal三角形

根据最小势能原理建立的有限元，是以结点位移作为基本未知量，这种有限单元称为位移元。由位移元得到的近似解答总体上是精确解的一个下限。

2.4 高精度的三角形单元、矩形单元
2.4.1 高精度三角形单元
(1) 六结点三角形单元（二次单元T6）
位移模式取坐标的完整二次式：
u=a1+a2 x+a3 y+a4 x2+a5 xy+a6 y2
v=b1+b 2 x+b 3 y+b 4 x2+b5 xy+b6 y2
该位移模式包含了常数项和完整的一次项，满足完备性要求；在单元的边界上位移呈二次抛物线分布，可由三个结点位移唯一确定，故又满足协调性的要求，是一种协调单元。

(2) 十结点三角形单元（三次单元T10）
位移模式取坐标的完整三次式：
u=a1+a2 x+a3 y+a4 x2+a5 xy+a6 y2+a7x3+a8x2y+a9xy2+a10y3
v=b1+b 2 x+b 3 y+b 4 x2+b5 xy+b6 y2+b7x3+b8x2y+b9xy2+b10y3
该位移模式包含了坐标的完整一次式（常数项和纯一次项），满足完备性要求；在单元的边界上位移呈三次曲线分布，可由4个结点位移唯一确定，故又满足协调性的要求，是一种协调单元。

i
j
m
i
j
m
3
1
2

2.12 高精度三角形单元

(3) 面积坐标表示的三角形单元形函数
在推导三角形单元的列式时，直接利用整体坐标系（为直角坐标）下的位移模式将使得运算十分繁琐和复杂。如果采用三角形单元内的一种局部坐标——面积坐标作为自然坐标，则可以使列式推导大为简化。
定义单元内任一点P的无量纲面积坐标（Li, Lj, Lm）为
i(1,0,0)
P
j(0,1,0)
m(0,0,1)
Ai
Aj
Am
Li= Ai/ A （i, j, m）
各种单元的形函数：
3结点三角形单元（线性单元T3）：
Ni=Li （i, j, m）
6结点三角形单元（二次单元T6）：
Ni=(2Li-1) Li （i, j, m）
N1=4Lj Lm, （1, 2, 3；i, j, m）
2.13 三个结点单元示意图
面积坐标的特点：
1) 三角形三个角点的面积坐标：i(1,0,0)、j(0,1,0)、m(0,0,1)
三条边用面积坐标表示的方程为：
jm边 Li=0
mi边 Lj=0
ij边 Lm=0
2) 三个面积坐标不独立，其相互关系为
Li+ Lj+ Lm=1
面积坐标与直角坐标之间的关系：

或

2.4.2 四结点矩形单元（R4单元）
(1) 采用双线性的位移模式：
u=a1+a2 x+a3 y+a4xy
1
y
x
2
3
4
h
x
o
v=b1+b 2 x+b 3 y+b 4xy
可保证位移模式的完备性和协调性。
若写成形函数形式，则为

2.14 四结点矩形单元

这里的Ni（i=1, 2, 3, 4）可以先求出8个待定系数在获得，也可以根据形函数的性质直接构造，例如对N1，可设

该表达式满足在结点2、3、4取值均为零的性质；再令Ni(1)=1，可得待定系数

这里设矩形单元的边长各为2a、2b。

(2) 局部坐标下的形函数
设矩形长和宽各为2a和2b，在矩形形心为原点建立局部坐标（自然坐标）系xoh，它与直角坐标的变换式：
x
4
h
1
2
3
x= -1
x=1
h= -1
h=1
o

2.15 局部坐标系下的矩形单元

则四个角点的局部坐标分别为1(-1, -1)，2(1, -1)，3(1,1)，4(-1,1)。
位移模式为：
Ni（i=1, 2, 3, 4）可根据形函数的性质直接构造出来：
或

(3) 单元应变矩阵和应力矩阵
单元应变为：
{e}=[B]{D}e
其中
[B]=[B1 B2 B3 B4]，。
单元应力为：
{s}=[S]{D}e=[D][B]{D}e
其中
[S]= [S1 S2 S3 S4]，
[Si]=[D] [Bi]

2.5 C0连续型单元形函数的构造
在有限单元法中，根据结点布置方式的不同可将矩形单元分为两类，一类称为Lagrange矩形单元，另一类称为Serendipity矩形单元。前者在单元纵横网格线的交点上均布置结点，而后者仅在单元的边界线上布置结点，如图2.17所示。

(a) Lagrange矩形单元 (b) Serendipity矩形单元
2.16 两类矩形单元

2.5.1 Lagrange矩形单元
(1) 一维Lagrange插值多项式
过n个结点（坐标分别为x1, x2, …, xn）的一维Lagrange插值多项式为

x
x1
x2
x3
xn
xi
1

2.17 一维Lagrange插值多项式l1(x)

例如n=2，则有：
。
若令x1=0，x2= l，则
。
用该多项式作为形函数，可满足形函数的性质。

(2) Lagrange矩形单元的形函数
在矩形的各个网格交点上均布置结点，
如水平方向r+1个，竖向p+1个。于是
第I列J行结点i的相应形函数：
Ni= NIJ= lI(r)(x) lJ(p)(h )
其中，lI(r)(x)、lJ(p)(h )均为一维Lagrange插值多项式，

Ni在结点i为1，在其他结点处为零；单元边界上的结点数=形函数阶数+1，能够保证边界位移的唯一性和协调性。
这类单元包含较多的内部结点，增加了单元的自由度，而实践证明这些自由度的增加通常并不能有效提高单元的精度。对n´n结点的Lagrange单元，其Pascal三角形中包含的项数如图2.19所示。
(I,J)
(r, p)
(0, 0)
1
1
lI (x)
lJ (h)
1
xn
yn
xnyn

x
x2
x3
y
y2
y3

图2.18 Lagrange矩形单元及插值模式 图2.19 Pascal三角形包含项数

2.5.2. Serendipity矩形单元
(1) 结点布置特点
这类单元只在其边界上布置结点，但不同边界上可布置有不同数目的结点。

(2) 形函数的构造方法
如果只有四个角点上布置结点（R4单元），则

如果增加结点5使之成为5结点矩形单元（或称R5单元），则

满足(N5) j=d5j （j=1, 2, 3, 4, 5）。原不满足，需对其进行修正才能满足：(N1)5=0，(N2)5=0，且保证N1，N2在除5结点外的其他4个结点取值不变。
x
4
h
1
2
3
5

1
0
1

1/2
1
N5

图2.20 五结点矩形单元 图2.21 形函数N1的修正

对于8结点Serendipity矩形单元（R8单元），利用上述方法，可得到前4个结点对应的形函数分别为

或写成统一表达式：

后4个形函数分别为

对于p次单元的边界内结点，其相应的形函数为x（或h）的p次和h（或x）的一次Lagrange多项式的乘积；而角结点的形函数为四结点单元的相应函数与各内结点形函数乘以一分数的差。
利用该方法可以很方便地构造一些过渡单元的形函数。

2.6 平面等参单元
直接用形状规则的单元对复杂形状的结构进行离散绘遇到边界难于处理的问题，因此应设法采用适当的方法对规则单元进行坐标变换，使之转化为斜边或曲边的单元。有限元中最常用的变换方法是等参变换，即坐标变换采用与单元位移插值一致的形函数。通过等参变换的单元称为等参单元。

2.6.1 四结点四边形等参元
对于4结点任意四边形单元，若采用双线性的位移模式，则位移在单元斜边界上呈二次抛物线分布，由两个结点位移不能唯一确定。为此在单元内设法建立一个局部坐标（自然坐标），使得单元各边界的局部坐标分别为一常量（或±1），这样，如果位移模式采用局部坐标的双线性函数，则能够满足公共边界位移的唯一性。
如果单元边界的局部坐标分别为±1，则在局部坐标下，单元的形态就是图示的4结点正方形单元，该单元称为基本单元或母单元，其形函数为
。
位移模式为：

1
y
x
2
3
4
x
h
x
4
h
1
2
3
x= -1
x=1
h= -1
h=1
o

(a) 基本单元 (b) 实际单元
2.22 四结点等参数单元

如果同时利用上述形函数作为局部（自然）坐标(x, h)向整体（直角）坐标(x, y)的变换函数，则可以建立两种坐标之间的映射关系如下：

根据形函数的性质：(Ni) j=dij，SNi=1，容易验证上述变换式在四个结点处给出了节点的整体坐标，而在单元的四边上，其中一个局部坐标为±1，另一局部坐标按线性变化，因而正好给出了整体坐标下四边的直线方程（由两个结点整体坐标唯一确定）。例如对12边，
因此上述变换式正确地反映了局部（自然）与整体（直角）坐标之间的映射关系。经变换后的实际单元称为子单元。子单元的位移模式仍为：

利用等参变换可以构造8结点、12结点、20结点等更高次的四边形曲边等参单元（参见图2.23）。
h
x
o
1
5
2
4
7
3
8
6
y
x
1
5
2
4
7
3
8
6

图2.23 8结点曲边四边形等参单元

2.6.2 局部与整体坐标的微分和积分变换式
根据复合函数的求导法则，

写成矩阵形式，

式中J称为雅可比（Jacobi）矩阵，对于具有m个结点的平面等参单元，

若反过来用局部坐标表示整体坐标，则可对上式作反变换，

式中
，
Jij*是J各元素Jji的代数余子式；½J½是J的行列式，称为雅可比（Jacobi）行列式。
面积微分的变换：
dA=dxdy=½J½dxdh。

2.6.3 单元刚度矩阵和单元等效结点荷载向量
对于一具有m个结点的平面等参单元，其应变向量可写为
{e}=[B]{D}e
其中，单元应变矩阵：
[B]=[B1 B2 … Bm]
。
单元应力向量：
{s}=[S]{D}e=[D][B]{D}e
单元刚度矩阵：

单元等效结点荷载列阵：

式中{Rp}e、{Rq}e分别为体积力和表面力引起的等效结点荷载，且

对其中一个局部坐标（x或h）为常数的边界，线积分ds为

或

上述积分式通常采用Newton-Cotes或高斯（Gauss）数值积分求得。
可见，等参单元的刚度矩阵、等效结点荷载列阵等的计算都在规则的母单元中进行，因此各种形式的积分运算都可以采用标准的数值积分方法进行，使得不同工程问题的有限元分析能够采用统一的通用化程序实现。

2.6.4 等参变换的应用条件
等参单元的应用条件是母单元与实际单元之间存在一一对应关系。具体到局部与整体坐标的变换式，要求变换矩阵即雅可比[J]非奇异，雅可比（Jacobi）行列式½J½¹0。
为确保坐标变换一一对应，在单元划分时应避免：
(1) 使任意