有限元原理 结构矩阵分析(平面桁架 平面应力) 变分

第二章 结构矩阵分析

由于有限元方法起源于力学中的结构分析，本章的作用是通过三个典型问题说明有限元方法应用于结构分析时的一般步骤，并借此了解有限元方法的一些基本概念。

§2-1平面桁架

（直接法，结构矩阵分析中常用的力法，处理静定问题，位移法，

可处理静定&静不定）

本节讨论的对象是图2-1所示的平面桁架。组成桁架的各杆为等截面直杆，外载荷p 直接作用于杆的铰接点（结点）上。为简单起见不妨设各杆的截面积均为Ａ，材料的弹性模量均为Ｅ。我们可按下述步骤求得桁架的变形和内力。

x

图２－２

图２－１

１、结构的离散化 对结点及单元编号

取组成桁架的每根杆为一个单元(该问题本身为一离散结构的力学问题) ，以①, ②, ③ 加以编号；取杆的铰接点为结点，以１、２、３加以编号（总体结点序号）。如图２－２所示，即：我们所讨论的桁架包括三个单元、三个结点。各单元（杆）仅在结点处连接。

２、建立总体坐标系 并确定结点坐标和自由度

为了描述结构的平衡需要建立一个坐标系，称为总体坐标系，以区别于以后出现的“局部坐标系”。总体坐标系的选择原则上不受限制，但希望

使用方便。本节所选的总体坐标系示于图２－２，坐标原

点与结点１重合。以u, v 分别表示沿 x, y 方向的位移分量， p, q分别表示力沿 x, y 轴的力分量（投影）。 在总体坐标系中各结点的坐标为： （x 1, y 1）=(0, 0 )、（x 2, y 2）=(a, a )、（x 3,

它们将作为程序的输入数据（几何参数）。

每个结点有两个自由度，对结点１、２、３分别为

图２－３

若暂时不考虑支承约束条件，整个结构的结点自由度为 {u 1 v 1 u 2 v 2 u 3 v 3 }T

３、单元分析（建立结点力与结点位移之间的关系）

取一个一般性的单元，设它的两个结点在结构中的编号为i, j （单元内部的结点序号）。由材料力学可知，杆的轴向刚度为EA/L。其中Ｌ为杆的长度：

22 L =x j -x i +y j -y i

（１）单元局部坐标系

现选取一典型单元对其进行单元分析，对所分析的单元按如下方式建立一个坐标系： 原点：与结点i 重合, x ’轴：沿 i ,j方向, y ’轴：与x ’轴垂直。

如图2-3所示。这个坐标系只属于一个单元，故称为单元局部坐标系，不同单元的单元局部坐标系一般是不相同的。在单元局部坐标系中可以规定：

结点自由度 { u ’i v’i }T ， { u ’j v’j }T ； 单元结点自由度 { u ’ }＝{ u ’i v’i u ’j v’j }T 。 （２）局部坐标系中的单元刚度矩阵

在外载荷作用下，结构发生变形，单元必受到来自结点的作用力。桁架中的杆只承受轴向力Ｓ，大小与杆的轴向伸长△Ｌ成正比

EA S =∆L

L

在局部坐标系中这种特性可以得到清楚的表述（这一点也是引入局部坐标系的理由之一）。若以 p ’i , q ’i , p ’j , q ’j 分别表示结点i, j作用于单元的力在 x ’, y ’ 轴上的投影，由①号单元的静力平衡有（图2-3）有

T {u ，

{u 2 ，v 2}T 、 {u 3， v 3}T 1 v 1} 、

用矩阵的形式可以写成

⎧p i ' ⎫⎡10-10⎤⎧u i '⎫

⎪q ' ⎪⎢⎥⎪⎪

EA ⎢0000⎥⎪v i '⎪⎪i ⎪ ⎨⎬=⎨'⎬p ' ⎢⎥-1010L j ⎪⎪⎪u j ⎪

⎢⎥⎪'⎪⎪⎪q ' 0000⎣⎦⎩v j ⎭ ⎩j ⎭

T

若引入单元广义力矢量： ' = p i ' q i ' p 'j q 'j 则上式可缩写为

EA

(u 'j -u i ') L

EA

p i '=-s =(u i '-u 'j )

L

q 'j =q i '=

0p 'j =

(2-1-1)

{}

r '=[k ']u '(2-1-2)

其中 ⎡10-10⎤ ⎢⎥

EA ⎢0000⎥

(2-1-3) [k ']=

⎢⎥L -1010

⎢⎥

0000⎣⎦

称为局部坐标系中的单元刚度矩阵, 它只与杆的几个参数E 、A 、L 有关，与杆的方位无关。

（３）坐标变换

局部坐标系中的单元刚度矩阵公式简捷。但不同单元的局部坐标系一般不同，为了研究结构整体的平衡，必须将结点给单元的力以及相应的单元刚度矩阵转换到统一的坐标系──总体坐标系。在总体坐标系中

单元结点自由度 { u }＝{ u i v i u j v j }T 结点给单元的力 ＝{ p i q i p j q j }T 在图2-3中，x ’ 轴与x 轴的夹角为 α x j -x i y j -y i 　　 cos α = , 　　 sin α =

L L

结点的位移分量的坐标变换为 ⎧u i '⎫⎡cos αsin α⎤⎧u i ⎫⎧u i ⎫

⎨⎬=⎢ ⎥⎨v ⎬=[t ]⎨v ⎬'v -sin αcos α⎦⎩i ⎭⎩i ⎭⎣⎩i ⎭

⎧u 'j ⎫⎡cos αsin α⎤⎧u j ⎫⎧u j ⎫

⎨'⎬=⎢⎥⎨v ⎬=[t ]⎨v ⎬v -sin αcos α⎦⎩j ⎭⎩j ⎭⎣⎩j ⎭

单元的位移分量的坐标变换为

⎧u i '⎫⎡cos αsin α00⎤⎧u i ⎫⎧u i ⎫

⎪v '⎪⎢⎪v ⎪⎪v ⎪

00⎥⎪i ⎪⎢-sin αcos α⎪⎪i ⎪i ⎪⎥[]==T ⎨'⎬⎨⎬⎨⎬ u u ⎢⎥00cos αsin αj j ⎪⎪⎪⎪⎪u j ⎪

⎢⎥⎪⎪ ⎪⎪⎪'v 00-sin αcos α⎦⎩v j ⎭⎩j ⎭⎣⎩v j ⎪⎭

或缩写为 = [T ] '} (2-1-４ )

类似，{r ’} 与 {之间的转换关系为 (2-1-５) r '=[T ]r 由于

⎡cos αsin α⎤

[t ]=⎢⎥(2-1-6) ⎣-sin αcos α⎦

是正交矩阵，因此

⎡t 0⎤[T ]=⎢⎥

⎣0t ⎦

(2-1-7)

也是正交矩阵。所以有 [T ]-1=[T ]T

将（2-1-4）、（2-1-5）代入（2-1-2）有

[T ]r }=[k '][T ]u }

从上式可得到 =T T k 'T =k }(2-1-8) 其中 {k }=T T k 'T (2-1-9)

称为单元在总体坐标系中的单元刚度矩阵。 以后将会看到，（2-1-9）是一个具有普遍意义的公式。它表明，当单元的自由度由一种形式换成另一种形式时，单元刚度矩阵只需进行一次相似变换。对于平面桁架单元，将（2-1-3）、（2-1-6）、（2-1-7）代入（2-1-9）可得到更便于应用的单元刚度矩阵公式

22

⎡cos αcos αsin α-cos α-cos αsin α⎤

⎢⎥22cos αsin αsin α-cos αsin α-sin αEA ⎥[k ]=⎢22

L ⎢-cos α-cos αsin αcos αcos αsin α⎥(2-1-10)

⎢⎥22 -cos αsin α-sin αcos αsin αsin α⎢⎥⎣⎦

（４）具体结果

由（2-1-10）可求得各单元的刚度矩阵的具体形式如下：

单元①：单元自由度 {u 1 v 1 u 2 v 2 }T ， α 单元刚度矩阵为 =45

[][][][]

[][][]

⎡2⎢⎢⎢24[k ]1=EA ⎢⎢a -2⎢⎢4⎢-2⎢⎣4

224-24-24-24-2224

-2⎤

⎥4⎥-2⎥⎥2⎥⎥⎥2⎥⎥4⎦

(2-1-11)

单元②：单元自由度 {u ２ v ２ u ３ v ３ }T ， α 单元刚度矩阵为 = -90

[k ]2

⎡00⎢

EA ⎢01=

a ⎢00⎢

⎣0-100⎤0-1⎥⎥00⎥

⎥01⎦

单元刚度矩阵为

(2-1-12)

单元③：单元自由度 {u 1 v 1 u ３ v ３ }T ， α = 0

[k ]3

⎡1⎢EA ⎢0=

a ⎢-1⎢⎣0

0-1000100

0⎤0⎥⎥0⎥⎥0⎦

(2-1-13)

请注意，单元刚度矩阵与单元自由度中位移分量的排列次序有关。如果改动这种排列次序，例如对①号单元，将单元自由度次序由{u 1 v 1 u 2 v 2 }T 改为{ u2 v 2 u1 v 1 }T ，必然导致刚度矩阵（2-1-11）元素位置的变动。

（５）单元刚度矩阵的物理意义和特点

设平面桁架单元在总体坐标系中刚度矩阵的一般形式为

⎡k 11k 12k 13k 14⎤

⎢k ⎥k k k 21222324⎥ [k ]=⎢

⎢k 31k 32k 33k 34⎥

⎢⎥

k k k k 424344⎦⎣41

T

由（2-1-8），当单元结点位移为{1 0 0 0 }时，在单元各结点上施加的力刚好为单元刚度

T

矩阵中的第一列：{k11 k 21 k31 k41 }。对[k]的其他各列也可做出类似的解释。即单元刚度矩阵的每一列相当于一组特定位移下的结点力，如表2-1所示。由图2-4可以获得更为直观的理解。

图２－４

对图2-4中的各种情况，据平面力系的平衡条件应有

k 1s +k 3s =0 k 2s +k 4s =0

x -x i k 4s -y j -y i k 3s =0　　（¨s =1~4）

j

这三个关系说明，[k]的四个行向量中只有一个线性独立（四个元素有三个约束方程）。

从以上分析可以看出，一般的单元刚度矩阵均具备以下两个特征。（对平面桁架单元而言，从（2-1-10）也可以得出这些结论）

（i ）单元刚度矩阵是对称矩阵，这是线性系统互易定理的具体体现。由于对称性，对行向量或列向量两者之一得到的结论，对另一个也适用。

(ii)单元刚度矩阵是奇异矩阵。它的行向量（或列向量）线性相关，具有零特征值，det[k]=0。对平面桁架的单元刚度矩阵而言，它的四个行向量（或列向量）中只有一个线性独立，而[k]有三个零特征值。这三个零特征值对应的特征向量相当于三种独立的刚体位移模式：两个平移，一个旋转。这是我们在单元分析中不考虑位移约束条件的自然结果。

４、总体刚度矩阵的组装 总体平衡方程

将图2-1所示的桁架中的支承约束以约束反力代替，如图2-5所示。下面来建立平衡问题的有限元方程。

（１） 结点平衡条件

作用于图2-5每个结点上的外载荷、支座反力以及来自单元的力应处于平衡。

以p i (m)、q i (m) 示结点i 作用于单元m 的力在x, y 轴上的投影，则单元m 给结点 i 的力在x, y 轴上的投影应为 －p i (m) 、－q i (m)。 对结点１：

(1) (3)

R 1X -p 1-p 1=0

(1) (3) R 1Y -q 1-q 1=0

对结点２：

(1) (2) P -p -p =0

()()

对结点３；

22

-q 2-p 3

(1)

-q 2-p 3

(2)

(2)

=0=0

(3)

图２－５

(2) (3)

R 3Y

-q 3-q 3

=0

可以合并成

⎧p 1(1) ⎫⎧0⎫⎧p 1(3) ⎫⎧R 1X ⎫

⎪(1) ⎪⎪⎪⎪(3) ⎪⎪R ⎪0q ⎪q 1⎪⎪1Y ⎪⎪1⎪⎪⎪(2) (1) ⎪⎪p 2⎪⎪⎪⎪p 2⎪⎪⎪⎪0⎪⎪⎪⎪P ⎪⎪

⎨(1) ⎬+⎨(2) ⎬+⎨⎬=⎨⎬

q 0⎪q 2⎪⎪2⎪⎪0⎪⎪⎪(2-1-1４) (2) (3) ⎪0⎪⎪p 3⎪⎪p ⎪⎪0⎪

3

⎪⎪⎪⎪⎪⎪⎪(2) (3) ⎪⎪⎪⎪⎪⎭⎩0⎪⎭⎩q 3⎭⎪⎩q 3⎪⎭⎩R 3Y ⎪

式（2-1-14）的右边为外载荷和支反力。左边则为单元给结点的力，它们是未知的，但可以借助单元刚度矩阵以结点位移来表示。

（２） 单元刚度矩阵的扩充

为了表示（2-1-14）左边的各个列向量，设想将每个单元的自由度扩充到与结构总体自由度相同（本例为６），并在单元刚度矩阵中补充零元素，由（2-1-11）、（2-1-12）、（2-1-13）和（2-1-8）可以用结点位移表示（2-1-14）左边的各列向量。 由单元①

⎡2⎤2-2-2 00⎥⎢

44⎢⎥⎧u 1⎫⎧p 1(1) ⎫⎧u 1⎫

22-2-2⎢⎥⎪⎪⎪(1) ⎪⎪v ⎪ 00v q 1⎢4⎥⎪⎪⎪1⎪⎪1⎪444⎢⎥⎪⎪⎪(1) ⎪⎪⎪

⎪p 2⎪EA -2⎢⎨(1) ⎬=

a ⎢4⎪q 2⎪

⎢-2⎪0⎪

⎢⎪⎪

⎪⎢4⎩0⎪⎭

⎢0⎢⎣0

-2

4-[1**********]400

0000

⎪u 2⎪⎪u 2⎪

[]=k ⎥0⎨⎬⎬1⎨v v ⎥⎪2⎪⎪2⎪(2-1-15) ⎥⎪⎪⎪u 3⎪0⎥⎪u 3⎪⎪⎪

⎪⎥⎪v ⎪⎩v 3⎪⎭ 0⎥⎩3⎭

0⎥⎦

由单元②

⎧0⎫⎡0

⎪0⎪⎢0

⎪(2) ⎪⎢

⎪ ⎪p 2⎪⎪EA ⎢0⎨(2) ⎬=⎢ q a ⎪2⎪⎢0

⎪p 3(2) ⎪⎢0

⎪(2) ⎪⎢

⎪⎪⎢⎣0⎩q 3⎭

由单元③

⎧p 1(3) ⎫⎡

⎪(3) ⎪⎢ q ⎪1⎪⎢

⎪⎪⎪0⎪EA ⎢

000000

0000

00010

00000

00-10

0⎤⎧u 1⎫⎧u 1⎫⎪v ⎪⎪v ⎪0⎥⎪1⎪⎥⎪1⎪

⎪0⎥⎪⎪u 2⎪⎪⎪u 2⎪⎪

⎥⎨⎬=]2⎨⎬-1⎥⎪v 2⎪⎪v 2⎪

⎪u 3⎪0⎥⎪u 3⎪

⎥⎪⎪⎪⎪⎪⎪⎪1⎥⎦⎩v 3⎭⎩v 3⎪⎭

(2-1-16)

1

00

⎨⎬=⎢

a 0⎪⎪⎢0

⎪p (3) ⎪⎢-13

⎪(3) ⎪⎢⎪⎢⎣0⎩q 3⎪⎭

[1**********]0-[1**********]

0⎤⎧u 1⎫⎧u 1⎫⎪v ⎪⎪v ⎪0⎥⎪1⎪⎥⎪1⎪

⎪0⎥⎪⎪u 2⎪⎪⎪u 2⎪⎪

⎥⎨⎬=[]3⎨⎬0⎥⎪v 2⎪⎪v 2⎪

⎪u 3⎪0⎥⎪u 3⎪

⎥⎪⎪⎪⎪⎪⎪⎪0⎥⎦⎩v 3⎭⎩v 3⎪⎭

(2-1-17)

（３） 组装总体刚度矩阵 将（2-1-15）、（2-1-16）、（2-1`-17）代入（2-1-14）得到 或

(]1+[]]2

⎧u 1⎫⎧R 1X ⎫⎪v ⎪⎪R ⎪⎪1⎪⎪1Y ⎪⎪⎪u 2⎪⎪⎪⎪P ⎪⎪

+[]]3)⎨⎬=⎨⎬

v 02⎪⎪⎪⎪⎪u 3⎪⎪0⎪⎪⎪⎪⎪⎪⎪⎪v R ⎩3⎭⎩3Y ⎪⎭

⎧u 1⎫⎧R 1X ⎫

⎪v ⎪⎪R ⎪⎪1⎪⎪1Y ⎪⎪u 2⎪⎪⎪⎪P ⎪⎪[]⎪=⎨⎬⎨⎬

v 0⎪2⎪⎪⎪⎪u 3⎪⎪0⎪⎪⎪⎪⎪⎪⎪⎪v R ⎩3⎭⎩3Y ⎪⎭

(2-1-18)

其中

]=]1+]2+]3=∑[]m

(2-1-19)

上式称为没有考虑位移约束条件情况下的总体刚度矩阵（求和对所有单元进行）。对本

节所分析的平面桁架有

⎡⎢⎢⎢⎢EA ⎢⎢a ⎢⎢⎢⎢⎢⎢⎣

22

+142244-22

-4-22

-44-1000-2

42-22400

⎤-2

-10⎥4⎥⎧u 1⎫⎧R 1X ⎫2-00⎥⎪v 1⎪⎪R 1Y ⎪

⎥⎪⎪⎪⎪⎥⎪⎪u 2⎪⎪⎪⎪P ⎪⎪2

=⎥00⎨⎬⎨⎬v 0⎥⎪2⎪⎪⎪⎥⎪⎪⎪2⎪

+10-1⎥⎪u 3⎪⎪0⎪4⎥⎪v ⎪⎪R ⎪010⎥⎩3⎭⎩3Y ⎭-101⎥⎦

(2-1-20)

要形成总体平衡方程（2-1-18），只需组装出它的右端项和总体刚度矩阵[k]就足够了，

这只是同一件事的两种不同的提法而已。

（2-1-19）用来“书写“组装总刚度矩阵的过程是简单而明了的。但事实上不能这样做，将所有[k]m 补充零元素扩充为[k]m 极大地浪费了宝贵的存贮空间，这些零元素仅起到使[k]m 的元素在总刚度矩阵中就位的作用。实际上采用的是另一种方法。如果已选定各结点位移在结构总体自由度中的排列次序为{u1 v 1 u 2 v 2 u 3 v 3 }T ，对每个单元在形成单元刚度矩阵的同时还形成了个定位数组LM ，它将指出单元自由度的各分量在总体自由度中的序号，如下表：

单元刚度矩阵中第s 行第t 列的元素k st 加到总刚度矩阵的第LM(s)行LM(t)列即可。这一组装总体刚度矩阵的方法被形象的称为“对号入座”。

从（2-1-20）可以看出总刚阵[k]是奇异阵，它的六个行向量（或列向量）中只有三个线性独立。这是尚未考虑位移约束条件，结构的刚体位移未受到限制的必然结果。

（４） 引入位移约束条件

由图2-1，平面桁架的位移约束条件为

u 1=v1=v3=0 (2-1-21) 代入方程（2-1-18），得到

⎡

⎢⎢⎢⎢EA ⎢⎢a ⎢⎢⎢⎢⎢⎢⎣

22

+142244-22

-4-22

-44-1000-2

42-422400

⎤-2

-10⎥4⎥⎧0⎫⎧R 1X ⎫2-00⎥⎪0⎪⎪R 1Y ⎪

⎥⎪⎪⎪⎪4

⎥⎪⎪u 2⎪⎪⎪⎪P ⎪⎪2

=⎥00⎨⎬⎨⎬

0v ⎥⎪2⎪⎪⎪

⎥⎪u ⎪⎪0⎪2

+10-1⎥⎪3⎪⎪⎪4⎪⎪⎪⎪R 0⎥010⎥⎩⎭⎩3Y ⎭-101⎥⎦

它显然可以分成两个方程组

⎡

⎢⎢EA ⎢a ⎢⎢⎢⎢⎣

2240

⎤2

0⎥⎥⎧u 2⎫⎧P ⎫2⎪⎪⎪⎪+10⎥⎨v 2⎬=⎨0⎬

⎥4⎪⎪⎪⎪

01⎥⎩u 3⎭⎩0⎭

⎥⎥⎦-242-4-1

⎤-1⎥

⎥⎧u 2⎫⎧R 1X ⎫⎪⎪⎪⎪0⎥⎨v 2⎬=⎨R 1Y ⎬⎥

⎪⎪⎪⎪0⎥⎩u 3⎭⎩R 3Y ⎭⎥⎥⎦

(2-1-22)

和

⎡-2

⎢⎢4EA ⎢2

-a ⎢4⎢0⎢⎢⎣

(2-1-23)

（2-1-22）就是平常所说的平衡问题的有限元方程，一般把它写成 (2-1-24) [K ]{U }={F }

[K]为考虑了约束条件（2-1-21）后的总体刚度矩阵。{U}为非约束自由度位移向量，{F}为作用于这些自由度上的载荷向量。矩阵[K]相当从矩阵[k]中划去了第１、２、６行和第１、２、６列。这说明（2-1-21）这种零位移约束条件，可以通过对总刚度矩阵划行划列的方法来实现。至于其他形式的位移约束条件的实现方法，后面将有专门的一章加以讨论。实际上，在组装总刚度矩阵的过程中，只要不组装应划去的行和列即可直接得到[K]和{F}，即方程（2-1-24）。

至于方程（2-1-23），尽管求得结点位移后可以由它求得支承反力，但是有限元分析中一般不形成这个方程。如确有必要求支座反力，将另找其他途径。

５、解有限元方程

对于平衡问题，归结为解一个代数方程。本例的自由度较少，从（2-1-22）可以方便地求得非约束自由度的位移

aP ⎫⎧

(22+1) ⎪EA ⎪⎧u 2⎫⎪⎪ aP ⎪⎪⎪⎪v 2⎬=⎨-⎨⎬ EA ⎪u ⎪⎪⎪

0⎩3⎭⎪ ⎪⎪⎪ ⎩⎭

被约束自由度根据（2-1-21）可直接赋零

⎧u 1⎫⎧0⎫

⎪⎪⎪⎪

⎨v 1⎬=⎨0⎬

⎪v ⎪⎪0⎪ ⎩3⎭⎩⎭

６、求单元内力

桁架单元的内力只有一个轴力Ｓ。由（2-1-1）S 的大小和正负与 p j ’相同。由（2-1-2）和（2-1-4）

'=[k ']'}=[k '][T ]}

⎧u i ⎫⎪v ⎪⎪i ⎪sin α]⎨⎬

⎪u j ⎪⎪⎩v j ⎪⎭

取出它的第三行即得到

EA

[-cos αS =

L

-sin αcos α

(2-1-24)

对单元①、②、③ 求得内力Ｓ于下表

由于本节所讨论的桁架是一个十分简单的静定桁架，用理论力学的知识即可得到各杆的

内力，结果相同。但是，若桁架为静不定桁架且杆的数目有上千个，那么本节所讨论的方法原则上不会遇到任何困难，我们要解决的课题将转为如何管理有关的大量数据和如何解一个数千阶的代数方程组这样一些技术问题。

§2-3 平面应力问题 常应变三角形

图2-8为一边长为a 、厚度为t 的正方形薄板。其中AB 边固定，BC 、CD 边自由，AD 边作用均布压力q 。对这一问

题，有限元分析的步骤是：

１、将ABCD 划分（离散）为８个

三角形（单元），编号①—⑧。各单元仅在顶点（结点）铰接，结点编号1—9。建立坐标系后，不难定出各结点的坐标（x i , yi ）。 ２、单元分析 任取一个一般性的单元，如图２－９所示。三个结点的编号为i, j, k。结点位移为( ui , vi

)

、( uj , vj ) 、( uk , vk ) 单元

结点位移为

图２－８ T

u }=u i v i u j v j u k v k

（１） 假定单元内位移场u, v 是x, y的一次函数

u =α1+α2x +α3y

(2-3-1)

{}

v =α4+α5x +α6y

图（２－９）

⎧u i ⎫⎡1x i ⎪⎪⎢

⎨u j ⎬=⎢1x j ⎪u ⎪⎢1x

k ⎩k ⎭⎣

y i ⎤⎧α1⎫⎧α1⎫α1~α6为待定常数，在结点处应有 ⎪⎪⎪⎪

[]y j ⎥α=DA ⎨⎬⎨α2⎬可解出：⎥2

⎪α⎪⎪α⎪⎡a i ⎧α1⎫⎧u i ⎫y k ⎥⎦⎩3⎭⎩3⎭

1⎢⎪⎪⎪-1⎪

⎨α2⎬=[DA ]⎨u j ⎬=⎢b i

2∆⎪α⎪⎪u ⎪⎢c i

⎩3⎭⎩k ⎭⎣

a j b j c j

a j b j c j

a k ⎤⎧u i ⎫

⎥⎪⎪b k ⎥⎨u j ⎬

⎪⎪c k ⎥⎦⎩u k ⎭

⎡a i ⎧α4⎫⎧v i ⎫

1⎢⎪⎪⎪-1⎪

b i ⎨α5⎬=[DA ]⎨v j ⎬=⎢

⎪α⎪⎪v ⎪2∆⎢c ⎩6⎭⎩k ⎭⎣i

其中

⎧v ⎡1x i i ⎫a k ⎤⎧v i ⎫⎪⎥⎪⎪⎪⎢

v ⎨⎬b k ⎥j ⎨v =⎢1x j j ⎬⎪v ⎪⎪⎪⎢1x ⎥k ⎭k ⎩⎣c k ⎦⎩v k ⎭

y i ⎤⎧α4⎫⎧α4⎫

⎪⎪⎪⎪

[]y j ⎥α=DA ⎨⎬⎨α5⎬⎥5

⎪⎪⎪α⎪y k ⎥⎦⎩α6⎭⎩6⎭

a i =

x j x k

y j y k

=x j y k -x k y j a j =-

x i x k

y i y k

=x k y i -x i y k a k =

x i x j

y i y j

=x i y j -x j y i

b i =-

y j

=y j -y k i

i y k

b j =

y i

=y k -y i i

y k

b k =-

y i

y j

=y i -y j

i

x j

c i ==x k -x j

x k x i

c j =-=x i -x k

x k

x i

c k ==x j -x i

x j

2∆=det DA =a i +a j +a k

1

(a i +b i x +c i y )u i +1a j +b j x +c j y u j +1(a k +b k x +c k y )u k u =

2∆2∆2∆ =N i (x , y ) u i +N j (x , y ) u j +N k (x , y ) u k

当 i, j, k 的位置为逆时针排列时，2Δ恒正，且等于三角形单元面积的两倍。将这些结果代入（2-3-1）有

类似可得到

v =N i (x , y ) v i +N j (x , y ) v j +N k (x , y ) v k

可以合并成

()

⎧u ⎫⎡N i ⎨⎬=⎢⎩v ⎭⎣0

0N i

N j 0

0N j

N k 0

⎧u i ⎫⎪v ⎪⎪i ⎪0⎤⎪⎪u j ⎪⎪⎡N i

⎨⎬=⎢N k ⎥⎦⎪v j ⎪⎣0

⎪u k ⎪⎪⎪⎪⎩v k ⎪⎭

0N i

N j 0

0N j

N k 0

0⎤

}N k ⎥⎦

（２）单元的应变、应力 利用（2-3-1）不难求得

⎡∂⎤⎡∂⎤

0⎥0⎥⎢⎢

⎧εx ⎫⎢∂x ⎥⎢∂x ⎥

∂⎥⎧u ⎫⎢∂⎥⎡N i 0N j ⎪⎪⎢

⎨εy ⎬=⎢0⎬=⎢0⎥⎨⎥⎢0v ∂y ∂y ⎩⎭⎣0N i ⎪ε⎪⎢⎢∂⎩xy ⎭∂∂⎥∂⎥

⎢⎥⎢⎥⎢⎢⎣∂y ∂x ⎥⎦⎣∂y ∂x ⎥⎦⎡b i 0b j 0b k 0⎤1⎢⎥=0c 0c 0c =[B ]}i j k ⎥2∆⎢

⎢c i b i c j b j c k b k ⎥⎣⎦

0N j

N k 0

0⎤

u }N k ⎥⎦

(2-3-2)

⎧σx ⎫⎧εx ⎫

⎪⎪⎪⎪(2-3-3) σ=E ⎨y ⎬⎨εy ⎬=E B ⎪σ⎪⎪ε⎪

⎩xy ⎭⎩xy ⎭

其中 ⎡b i 0b j 0b k 0⎤

⎥ [B ]=1⎢0c 0c 0c i j k ⎥

2∆⎢

⎢c i b i c j b j c k b k ⎥(2-3-4) ⎣⎦

在假定单元内位移场u 、v 是 x, y的一次函数的前提下，单元内的应变和应力将是常数，故这种单元又称为常应变三角元。

[][][]}

0(a)

F 3x )

(b)

图（２－10）

（３）为了在单元内构成均匀应力场，必须在单元的各边施加均布载荷，它们的合力一定

作用在各边的中点，如图2-10(a) 所示。再将各边上的合力平分到这边的两点结点，由图2-10 (b)不难得出

1t p i =(F 1x +F 3x )=y j -y i σx +x i -x j τxy -y k -y j σx -(x i -x k )τxy 22

t =y j -y k σx +x k -x j τxy 2

t =b i σx +c i τxy 2

类似可求得 q i 、p j 、q j 、p k 、q k 并可合并写成 ⎧p i ⎫⎡b i 0c i ⎤ ⎪q ⎪⎢0c b ⎥

i i i ⎥⎪⎪⎢⎧σx ⎫⎧σx ⎫

⎪⎪⎢⎥⎪⎪p j ⎪t b j 0c j ⎪⎪T ⎪T

[]=σ=t ∆B σ⎢⎥⎨y ⎬⎨⎬⎨y ⎬=t ∆[B ][E ][B ]⎪q j ⎪2⎢0c j b j ⎥⎪τ⎪⎪τ⎪

xy ⎭⎩⎩xy ⎭⎢b k 0c k ⎥⎪p k ⎪

⎢⎥⎪⎪

⎪⎪q 0c b ⎢k k ⎥⎩k ⎭⎣⎦

根据单元刚度矩阵[k]的直观意义，常应变三角元的单元刚度矩阵即为 [k ]=

[B ]T [E ][B ]⋅∆⋅t (2-3-5)

（４）单元①和②的边上作用着均布的外载荷，可以把它们的合力平分到两结点，如图2-11所示。 t t

()f =q x -x f =q (x 1-x 4)1472 22

q f [()()()]

[(

(

)()]

)

f

图 （２－１１）

（５）为了组装总体刚度矩阵，每个单元还应形成一个数组LM 。元素LM(1)~LM(6) 分别为 u i 、v i 、u j 、v j 、u k 、v k 在总体自由度中的序号。

由于单元的结点位移和结点力都是用总体坐标描述的，单元刚度矩阵不必再进行坐标变换。

３、组装总体刚度矩阵和载荷向量。为了实现位移约束条件

u 7=v 7=u 8=v 8=u 9=v 9=0

这些自由度对应的行和列可以不必组装。总体平衡方程为

[K ]{U }={F }(2-3-6)

其中, 非约束自由度位移为 {U }={u 1v 1u 2v 2u 3v 3u 4v 4u 5v 5u 6v 6}T

相应自由度的载荷向量为 T

{F }={0-f 200000-(f 1+f 2)0000}

解方程（2-3-6）可得到各非约束自由度的位移，再由（2-3-3）可求得各单元的应力。

§2-4结束语

１、本章所讨论的三具体问题尽管力学背景不同，但分析步骤大体相同。区别仅在于结点参数不同，单元分析公式不同。而组装总体矩阵的方法、处理约束条件的方法以及代数方程组的求解方法完全通用。

２、从建立有限元方程的方法看，本章属于直接方法一类。这种方法只在简单情况下有效，但是物理意义清晰。尽管今天有限元方法已发展成为一种抽象的数学分析方法，本章提供的直观解释仍然常常被用来说明一些问题。

３、对于第一、二两个问题（平面桁架和平面框架）本章求的将是精确解（真实解）。而对第三个问题（平面应力问题），本章求得的只是近似解（位移场是假设的），这样就自然引出以下两个问题：

（１） 所求的近似解与真实解的接近程度如何？或者说，怎样才能接近得更好一些？ （２） 在进行单元分析时，曾作了两条关键性的假设：

（i ）单元内位移u 、v 是坐标x, y的一次函数；

(ii)各单元仅在结点处铰接。

这样的假定有没有依据？是否必需？这些问题将在后续的章节中讨论。

变分原理简介

一、引言

变分涉及求以函数为变量的函数——泛函的极值的问题。

我们首先以三个非常简单的极值问题为例，引出泛函和变分的概念。这三个问题是：最短路径问题；最速下落问题；最小旋转曲面问题。

（a ）最短路径问题 求连接给定两点的最短平面曲线。坐标系取为直角坐标系，其中的一个点位于原点，另一个点的坐标为(x 1, y 1) 。如果的曲线，I 是两点间曲线的长度，显然

或

（1）

是一条通过（0,0）和（x 1,y 1）

所求问题为，求函数y ，使上述积分值为最小。

（b ）最速下落问题 求铅直面内一曲线，当粒子在重力作用下从曲线的一点由静止状态滑落到曲线上另一点时所用时间最少。建立如下图所示的坐标，起点位于原点，终点坐标为

（x 1,y 1）。任一时刻粒子运动的速度

为

。粒子从原点下落到（x 1,y 1）所需时间为

所求问题为，求函数y ，使上述积分值为最小。

（c ）最小旋转曲面问题 给定两点和一条直线，点和直线在一个平面内，求连接两点的一条曲线，使曲线绕直线旋转生成的曲面面积为最小。取直线为x 轴，其中的一个点在y 轴上，曲线端点的坐标为（0，y 0）和（x 1,y 1）。用S 表示旋转曲面的面积，则

记

所求问题为，求函数y ，使上述积分值为最小。

这三个问题的共同特点是，在给定条件下，求函数y=f(x)使积分

为最小。

我们先来求解第一个问题。最短路径对应的积分式可表示为

设所求曲线方程为y=f(x)，而y=F(x)是不同于f(x)的任意一条曲线的方程。两者之差为

，位于F(x)和f(x)之间的曲线可表示为

不依赖于x 的参数。如果α足够小，则

，其中α是

位于y=f(x)的邻域内。如果曲线

y=f(x)和y=F(x)是给定的，则积分的条件为

=0，即

。

仅是α的函数。α=0处积分取最小值

在α=0处

对上式进行分部积分

在端点处差值是

=0。由于

是任意的，故上式积分为零的的条件

求解该微分方程，得

代入边界条件

，得

，即所求曲线为过给定点的直线。

下面进一步求解第二个问题，导出求解泛函极值更一般的表达式。

以上三个问题的一般表达式为，参照问题一的求解过程，对给定的曲线

F(x)和f(x)，位于两者之间的其它曲线方程为，积分

仅是α

的函数。积分

对求α导数，得

对上式进行分部积分，得

该偏微分方程称作拉格朗日方程。

二、变分 符号 术语 实例

当x 固定不变时，差值

称作y 的变分, 记作δy 。变分和导数的区别：

'

y 的变分是在x 比变的情况下函数曲线发生变化引起的y 的变化，而y 的导数y =

'

dy

是对dx

'

固定的函数曲线，当x 变化时，y 对x 的变化率。同样地y 的变分表示函数曲线变化时y 的

' ' ' '

改变量，δy =δη=F (x ) -f (x ) 。δy =δ

'

dy d

=δy ，即变分和求导运算可交换顺序。dx dx

使用变分记号，泛函极值可表示为：，即

最速降落问题的解为