简洁博客

2026-07-04

第 2 章 刚度（位移）法导论

译稿来自《A First Course in the Finite Element Method》学习整理，保留原章节编号、公式编号和图号；内容仍建议结合原书校对。

本章目标

学完本章后，你应当能够：

引言

本章介绍直接刚度法所依赖的一些基本概念。我们首先讨论线性弹簧，因为它足够简单，同时又能很好地揭示刚度法的核心思想。我们先给出刚度矩阵的一般定义，再推导线弹性弹簧单元的刚度矩阵。随后，我们将用平衡与协调这两个基础概念，说明如何把若干弹簧单元组成的结构装配成整体结构刚度矩阵。接着，还会展示如何以直接叠加各单元刚度矩阵的方式得到组合体的整体刚度矩阵。“直接刚度法”这个名称正是由这种技术而来。

建立整体结构刚度矩阵之后，我们将说明如何施加边界条件，包括齐次边界条件和非齐次边界条件。这样便可以得到包括节点位移和支座反力在内的完整解答。（内力的求解将在第 3 章结合杆单元讨论。）

最后，本章引入最小势能原理：先用它推导弹簧单元方程，再用它求解弹簧组合体问题。我们会先在最简单、自由度数量较少的单元上说明这一原理，这样当后续章节不得不把它用于自由度很多的单元时，概念会更容易理解。

图 2-1 （a）单个弹簧单元；（b）三弹簧组合体

2.1 刚度矩阵的定义

熟悉刚度矩阵是理解刚度法的基础。我们把刚度矩阵定义如下：对一个单元而言，若矩阵 $[k]$ 能够满足

\[\{ f \} = [ k ] \{ d \}\tag{2.1.1}\]

则称 $[k]$ 为该单元的刚度矩阵。其中，$[k]$ 把单个单元的节点位移 $\{d\}$ 与节点力 $\{f\}$ 联系起来。例如图 2-1a 所示的弹簧就是一个单元。

对于由一系列单元组成的连续介质或结构，例如图 2-1b 所示的弹簧组合体，刚度矩阵 $[K]$ 把整体坐标系 $(x,y,z)$ 下的节点位移 $\{d\}$ 与整个介质或结构上的整体力 $\{F\}$ 联系起来，即

\[\{ F \} = [ K ] \{ d \}\tag{2.1.2}\]

其中，$[K]$ 表示整个弹簧组合体的整体刚度矩阵。

2.2 弹簧单元刚度矩阵的推导

下面用直接平衡法推导一维线性弹簧的刚度矩阵。所谓一维线性弹簧，是指遵循胡克定律，并且只在弹簧轴线方向抵抗力的弹簧。考虑图 2-2 所示的线性弹簧单元。参考点 1 和 2 位于单元两端，这些参考点称为弹簧单元的节点。与局部 $x$ 轴相关的局部节点力为 $f_{1x}$ 和 $f_{2x}$。局部轴沿弹簧方向，因此可以沿该轴直接度量位移和力。弹簧单元的局部节点位移为 $u_1$ 和 $u_2$。

这些节点位移称为各节点的自由度。如图所示，各节点力和位移的正方向均取为从节点 1 指向节点 2 的正 $x$ 方向。符号 $k$ 称为弹簧常数，或称弹簧刚度。

图 2-2 线性弹簧单元，以及节点位移和节点力的正号约定

在许多工程问题中，都能找到与实际弹簧常数相对应的量。第 3 章将看到，等截面单轴杆的弹簧常数为 $k=AE/L$，其中 $A$ 为杆的横截面积，$E$ 为弹性模量，$L$ 为杆长。类似地，第 5 章将说明，受扭的等截面圆杆具有弹簧常数 $k=JG/L$，其中 $J$ 为极惯性矩，$G$ 为材料的剪切模量。对一维热传导问题（第 13 章），$k=AK_{xx}/L$，其中 $K_{xx}$ 为材料的导热系数；对多孔介质中的一维流体流动问题（第 14 章），$k=AK_{xx}/L$，其中 $K_{xx}$ 为材料的渗透系数。

由此可以看到，刚度法不仅适用于结构问题，也可用于热传导、流体流动、电网络等非结构问题。关键在于选用合适的本构关系，例如结构问题中的胡克定律、热传导中的傅里叶定律、流体流动中的达西定律、电网络中的欧姆定律；同时配合某个守恒原理，例如节点平衡或能量守恒。

现在，我们希望建立弹簧单元的节点力与节点位移之间的关系。这个关系就是刚度矩阵。因此，我们希望把节点力矩阵与节点位移矩阵联系为

\[{ \left\{ \begin{array} { l } { f _ { 1 x } } \\ { f _ { 2 x } } \end{array} \right\} } = { \left[ \begin{array} { l l } { k _ { 1 1 } } & { k _ { 1 2 } } \\ { k _ { 2 1 } } & { k _ { 2 2 } } \end{array} \right] } { \left[ \begin{array} { l } { u _ { 1 } } \\ { u _ { 2 } } \end{array} \right] }\tag{2.2.1}\]

其中，式（2.2.1）中 $[k]$ 矩阵的单元刚度系数 $k_{ij}$ 是待确定的量。回忆式（1.2.5）和式（1.2.6），$k_{ij}$ 表示在第 $j$ 个自由度施加单位位移 $d_j$、且所有其他位移为零时，在第 $i$ 个自由度产生的力 $F_i$。也就是说，当 $d_j=1$ 且 $d_k=0\ (k\ne j)$ 时，有 $F_i=k_{ij}$。

下面按照第 1.4 节给出的步骤推导弹簧单元的刚度矩阵。不过，对于简单的一维弹簧单元，如果采用直接法，可以跳过步骤 2（选择位移函数）。第 3.2 节会说明选择位移函数的准则；之后在梁单元以及二维、三维单元刚度矩阵的推导中，我们会持续使用步骤 2，因为它能使后续推导更容易组织。由于本书的基本思路是先推导各种单元刚度矩阵，再说明如何用这些单元求解工程问题，所以这里的步骤 1 只涉及单元类型的选择。

步骤 1 选择单元类型

考虑图 2-3 所示的线性弹簧单元。该弹簧可能是弹簧系统中的一个单元，并承受沿弹簧轴线 $x$ 方向的节点拉力 $T$，这些拉力可能来自相邻弹簧的作用，并使单元处于平衡状态。局部 $x$ 轴从节点 1 指向节点 2。我们用两端节点编号和单元编号来表示该弹簧。变形前两节点之间的原始距离记为 $L$。单元的材料性质，即弹簧常数，记为 $k$。

步骤 2 选择位移函数

对于简单的一维弹簧单元，采用直接法时可以跳过步骤 2。第 3.2 节将介绍选择位移函数的准则，并在梁单元、二维单元和三维单元刚度矩阵的推导中使用这一步。

步骤 3 定义应变-位移关系和应力-应变关系

拉力 $T$ 使弹簧产生总伸长量（变形）$\delta$。图 2-4 给出了弹簧的典型总伸长情形。这里 $u_1$ 为负值，因为它的位移方向与正 $x$ 方向相反；而 $u_2$ 为正值。

弹簧的总变形可由节点位移之差表示为

\[\delta = u _ { 2 } - u _ { 1 }\tag{2.2.2}\]

对于弹簧单元，可以直接把弹簧内力与变形联系起来。因此，在这里不需要显式建立应变-位移关系。

应力-应变关系也可以等价地写成力-变形关系：

\[T = k \delta\tag{2.2.3}\]

将式（2.2.2）代入式（2.2.3），得到

\[T = k ( u _ { 2 } - u _ { 1 } )\tag{2.2.4}\]

步骤 4 推导单元刚度矩阵和单元方程

现在推导弹簧单元刚度矩阵。根据节点力的符号约定和平衡条件（见图 2-2 和图 2-3），有

\[f _ { 1 x } = - T f _ { 2 x } = T\tag{2.2.5}\]

由式（2.2.4）和式（2.2.5），可得

\[\begin{array} { l } { { T = - f _ { 1 x } = k ( u _ { 2 } - u _ { 1 } ) } } \\ { { \phantom { T = } } } \\ { { T = f _ { 2 x } = k ( u _ { 2 } - u _ { 1 } ) } } \end{array}\tag{2.2.6}\]

图 2-3 受拉力作用的线性弹簧

图 2-4 变形后的弹簧

将式（2.2.6）改写为

\[\begin{array} { r } { f _ { 1 x } = k ( u _ { 1 } - u _ { 2 } ) } \\ { f _ { 2 x } = k ( u _ { 2 } - u _ { 1 } ) } \end{array}\tag{2.2.7}\]

把式（2.2.7）写成一个矩阵方程，可得

\[{ \left\{ \begin{array} { l } { f _ { 1 x } } \\ { f _ { 2 x } } \end{array} \right\} } = { \left[ \begin{array} { l l } { k } & { - k } \\ { - k } & { k } \end{array} \right] } { \left[ \begin{array} { l } { u _ { 1 } } \\ { u _ { 2 } } \end{array} \right] }\tag{2.2.8}\]

这个关系适用于沿 $x$ 轴布置的弹簧。由刚度矩阵的基本定义，并将式（2.2.1）与式（2.2.8）对照，可得

\[[ k ] = { \left[ \begin{array} { l l } { k } & { - k } \\ { - k } & { k } \end{array} \right] }\tag{2.2.9}\]

这就是线性弹簧单元的刚度矩阵。这里的 $[k]$ 称为该单元的局部刚度矩阵。由式（2.2.9）可以观察到，$[k]$ 具有以下性质：

  1. 它是对称矩阵，即当 $i\ne j$ 时，$k_{ij}=k_{ji}$。这一性质可由 Rayleigh-Betti 互等定理证明 [4]。
  2. 它是方阵，即 $[k]$ 的行数等于列数，因为它把相同数量的节点力与节点位移联系起来。
  3. 它是奇异矩阵。也就是说，$[k]$ 的行列式等于零，因此 $[k]$ 不可逆。

步骤 5 装配单元方程以得到整体方程，并引入边界条件

利用节点力平衡方程、第 2.3 节中的力-变形关系和协调方程，以及第 2.4 节中介绍的直接刚度法，可以装配整体刚度矩阵和整体力矩阵。对于由多个单元构成的结构，有

\[[ K ] = \sum _ { e = 1 } ^ { N } [ k ^ { ( e ) } ] \qquad \mathrm { a n d } \qquad \{ F \} = \sum _ { e = 1 } ^ { N } \left\{ f ^ { ( e ) } \right\}\tag{2.2.10}\]

其中，$[k^{(e)}]$ 和 $\{f^{(e)}\}$ 分别表示已用整体参考坐标系表达的单元刚度矩阵和单元力矩阵。例如第 3 章讨论桁架结构时，这一概念将变得十分重要。（本书中，在这种语境下出现的求和符号 $\sum$ 并不表示对单元矩阵作普通代数加法，而是表示必须按照第 2.4 节所述的直接刚度法，把这些单元矩阵正确装配到整体矩阵中。）

步骤 6 求解节点位移

然后，通过施加边界条件（例如支座条件），并联立求解方程组，便可确定节点位移：

\[\{ F \} = [ K ] \{ d \}\tag{2.2.11}\]

步骤 7 求解单元力

最后，对每个单元进行回代，代入类似式（2.2.7）的单元方程，即可确定单元力。

2.3 弹簧组合体示例

桁架、建筑框架和桥梁等结构，都是由基本结构构件连接形成的整体结构。分析这类结构时，必须先确定由相互连接的单元系统所形成的整体结构刚度矩阵。在讨论桁架和框架之前，我们先利用第 2.2 节已经推导出的弹簧单元力-位移矩阵关系，并结合节点平衡与协调这两个基本概念，求出弹簧组合体的整体结构刚度矩阵。这样也就具体展示了步骤 5 的含义。

我们考虑图 2-5 所示的二弹簧组合体。这个例子足够一般，可以说明如何用直接平衡法建立弹簧组合体的整体刚度矩阵。这里节点 1 固定，并在节点 3 施加轴向力 $F_{3x}$，在节点 2 施加轴向力 $F_{2x}$。弹簧单元 1 和 2 的刚度分别为 $k_1$ 和 $k_2$。组合体的节点编号为 1、3、2；之所以这样编号，是为了体现更一般的情形，因为在大型问题中，相邻单元之间通常不会总是采用连续编号。

组合体的 $x$ 轴为整体轴。每个单元的局部 $x$ 轴都与组合体的整体轴重合。

对单元 1，由式（2.2.8）有

\[\left\{ \begin{array} { c } { f _ { 1 x } ^ { ( 1 ) } } \\ { f _ { 3 x } ^ { ( 1 ) } } \end{array} \right\} = \left[ \begin{array} { c c } { k _ { 1 } } & { - k _ { 1 } } \\ { - k _ { 1 } } & { k _ { 1 } } \end{array} \right] \left\{ { u _ { 1 } ^ { ( 1 ) } } \atop { u _ { 3 } ^ { ( 1 ) } } \right\}\tag{2.3.1}\]

对单元 2，有

\[\left\{ \begin{array} { c } { { f _ { 3 x } ^ { ( 2 ) } } } \\ { { f _ { 2 x } ^ { ( 2 ) } } } \end{array} \right\} = \left[ \begin{array} { c c } { { k _ { 2 } } } & { { - k _ { 2 } } } \\ { { - k _ { 2 } } } & { { k _ { 2 } } } \end{array} \right] \left[ { u _ { 3 } ^ { ( 2 ) } } \right]\tag{2.3.2}\]

此外，在整个变形过程中，单元 1 和单元 2 必须在公共节点 3 处保持连接。这称为连续性要求或协调条件。由协调条件可得

\[u _ { 3 } ^ { ( 1 ) } = u _ { 3 } ^ { ( 2 ) } = u _ { 3 }\tag{2.3.3}\]

图 2-5 二弹簧组合体

图 2-6 与单元力符号约定一致的节点力

其中，本书中位移 $u$ 上方括号内的上标表示它所对应的单元编号；右下标表示位移所在的节点。因此，$u_3$ 表示整体或全局弹簧组合体中节点 3 的位移。

图 2-6 给出了各单元和各节点的自由体图，且采用图 2-2 中已经建立的单元节点力符号约定。

根据图 2-6 中各节点的自由体图，并考虑每个节点处外力必须等于内力，可以写出节点 3、节点 2 和节点 1 的节点平衡方程：

\[F _ { 3 x } = f _ { 3 x } ^ { ( 1 ) } + f _ { 3 x } ^ { ( 2 ) }\tag{2.3.4}\]
\[F _ { 2 x } = f _ { 2 x } ^ { ( 2 ) }\tag{2.3.5}\]
\[F _ { 1 x } = f _ { 1 x } ^ { ( 1 ) } \tag{2.3.6}\]

其中，$F_{1x}$ 来自固定支座处外加反力。

从节点转到与该节点相连的单元时，需要应用牛顿第三定律：作用力与反作用力大小相等、方向相反。将式（2.3.1）至式（2.3.3）代入式（2.3.4）至式（2.3.6），可得

\[\begin{array} { l } { F _ { 3 x } = ( - k _ { 1 } u _ { 1 } + k _ { 1 } u _ { 3 } ) + ( k _ { 2 } u _ { 3 } - k _ { 2 } u _ { 2 } ) } \\ { F _ { 2 x } = - k _ { 2 } u _ { 3 } + k _ { 2 } u _ { 2 } } \\ { F _ { 1 x } = k _ { 1 } u _ { 1 } - k _ { 1 } u _ { 3 } } \end{array}\tag{2.3.7}\]

把式（2.3.7）写成矩阵形式，可得

\[{ \left\{ \begin{array} { l } { F _ { 3 x } } \\ { F _ { 2 x } } \\ { F _ { 1 x } } \end{array} \right\} } = { \left[ \begin{array} { l l l } { k _ { 1 } + k _ { 2 } } & { - k _ { 2 } } & { - k _ { 1 } } \\ { - k _ { 2 } } & { k _ { 2 } } & { 0 } \\ { - k _ { 1 } } & { 0 } & { k _ { 1 } } \end{array} \right] } { \left[ \begin{array} { l } { u _ { 3 } } \\ { u _ { 2 } } \\ { u _ { 1 } } \end{array} \right] }\tag{2.3.8}\]

若按照节点自由度编号从小到大的顺序重新排列式（2.3.8），则有

\[\left\{ \begin{array} { c } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \end{array} \right\} = \left[ \begin{array} { c c c } { k _ { 1 } } & { 0 } & { - k _ { 1 } } \\ { 0 } & { k _ { 2 } } & { - k _ { 2 } } \\ { - k _ { 1 } } & { - k _ { 2 } } & { k _ { 1 } + k _ { 2 } } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \\ { u _ { 3 } } \end{array} \right\}\tag{2.3.9}\]

式（2.3.9）现在可以写成一个统一的矩阵方程：

\[\{ F \} = [ K ] \{ d \}\tag{2.3.10}\]

其中，$\{F\}=\left\{ \begin{array} { l } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \end{array} \right\}$ 称为整体节点力矩阵；$\{d\}=\left\{ \begin{array} { l } { { u _ { 1 } } } \\ { { u _ { 2 } } } \\ { { u _ { 3 } } } \end{array} \right\}$ 称为整体节点位移矩阵；而

\[[ K ] = \left[ { \begin{array} { c c c } { k _ { 1 } } & { 0 } & { - k _ { 1 } } \\ { 0 } & { k _ { 2 } } & { - k _ { 2 } } \\ { - k _ { 1 } } & { - k _ { 2 } } & { k _ { 1 } + k _ { 2 } } \end{array} } \right]\tag{2.3.11}\]

称为总刚度矩阵、整体刚度矩阵，或系统刚度矩阵。

总之，为了建立弹簧组合体的刚度方程和刚度矩阵，即式（2.3.9）和式（2.3.11），我们使用了力-变形关系式（2.3.1）和式（2.3.2）、协调关系式（2.3.3），以及节点力平衡方程式（2.3.4）至式（2.3.6）。在第 2.4 节介绍一种更实用的总刚度矩阵装配方法、并在第 2.5 节讨论支座边界条件之后，我们将回到这个例题，给出完整求解。

2.4 用叠加法装配总刚度矩阵（直接刚度法）

现在考虑一种更方便的总刚度矩阵构造方法。这个方法的基础，是把构成结构的各个单元刚度矩阵进行正确叠加（也可参见文献 [1] 和 [2]）。

仍以第 2.3 节的二弹簧组合体为例。由式（2.3.1）和式（2.3.2）可得单元刚度矩阵：

\[[ k ^ { ( 1 ) } ] = \left[ \begin{array} { c c } { { u _ { 1 } } } & { { u _ { 3 } } } \\ { { k _ { 1 } } } & { { - k _ { 1 } } } \\ { { - k _ { 1 } } } & { { k _ { 1 } } } \end{array} \right] \ u _ { 1 } \qquad [ k ^ { ( 2 ) } ] = \left[ \begin{array} { c c } { { u _ { 3 } } } & { { u _ { 2 } } } \\ { { k _ { 2 } } } & { { - k _ { 2 } } } \\ { { - k _ { 2 } } } & { { k _ { 2 } } } \end{array} \right] \ u _ { 3 }\tag{2.4.1}\]

这里，写在 $[k]$ 矩阵列上方和行旁边的 $u_i$，表示该单元矩阵每一行、每一列所对应的自由度。

式（2.4.1）中的两个单元刚度矩阵并不对应同一组自由度：单元 1 对应节点 1 和节点 3 的轴向位移，而单元 2 对应节点 2 和节点 3 的轴向位移。因此，这两个单元刚度矩阵不能以当前形式直接相加（叠加）。要叠加单元矩阵，必须先把它们扩展到整体结构（弹簧组合体）刚度矩阵的阶数，使每个单元刚度矩阵都与结构的全部自由度相联系。扩展的方法很简单：对于某个单元不涉及的那些位移，在相应位置添加零行和零列。

对于单元 1，可把刚度矩阵改写成扩展形式，使式（2.3.1）变为

\[k _ { 1 } \left[ \begin{array} { c c c } { u _ { 1 } } & { u _ { 2 } } & { u _ { 3 } } \\ { 1 } & { 0 } & { - 1 } \\ { 0 } & { 0 } & { 0 } \\ { - 1 } & { 0 } & { 1 } \end{array} \right] \quad \left\{ \begin{array} { l } { u _ { 1 } ^ { ( 1 ) } } \\ { u _ { 2 } ^ { ( 1 ) } } \\ { u _ { 3 } ^ { ( 1 ) } } \end{array} \right\} = \left\{ \begin{array} { l } { f _ { 1 x } ^ { ( 1 ) } } \\ { f _ { 2 x } ^ { ( 1 ) } } \\ { f _ { 3 x } ^ { ( 1 ) } } \end{array} \right\}\tag{2.4.2}\]

由式（2.4.2）可以看出，$u_2^{(1)}$ 和 $f_{2x}^{(1)}$ 并不与 $[k^{(1)}]$ 相关。类似地，对单元 2 有

\[k _ { 2 } \left[ \begin{array} { c c c } { { u _ { 1 } } } & { { u _ { 2 } } } & { { u _ { 3 } } } \\ { { 0 } } & { { 0 } } & { { 0 } } \\ { { 0 } } & { { 1 } } & { { - 1 } } \\ { { 0 } } & { { - 1 } } & { { 1 } } \end{array} \right] \quad \left\{ \begin{array} { l } { { u _ { 1 } ^ { ( 2 ) } } } \\ { { u _ { 2 } ^ { ( 2 ) } } } \\ { { u _ { 3 } ^ { ( 2 ) } } } \end{array} \right\} = \left\{ \begin{array} { l } { { f _ { 1 x } ^ { ( 2 ) } } } \\ { { f _ { 2 x } ^ { ( 2 ) } } } \\ { { f _ { 3 x } ^ { ( 2 ) } } } \end{array} \right\}\tag{2.4.3}\]

现在考虑每个节点处的力平衡，可得

\[\left\{ \begin{array} { c } { f _ { 1 x } ^ { ( 1 ) } } \\ { 0 } \\ { f _ { 3 x } ^ { ( 1 ) } } \end{array} \right\} + \left\{ \begin{array} { c } { 0 } \\ { f _ { 2 x } ^ { ( 2 ) } } \\ { f _ { 3 x } ^ { ( 2 ) } } \end{array} \right\} = \left\{ \begin{array} { c } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \end{array} \right\}\tag{2.4.4}\]

式（2.4.4）实际上就是式（2.3.4）至式（2.3.6）的矩阵形式。将式（2.4.2）和式（2.4.3）代入式（2.4.4），得到

\[k _ { 1 } \left[ \begin{array} { c c c } { 1 } & { 0 } & { - 1 } \\ { 0 } & { 0 } & { 0 } \\ { - 1 } & { 0 } & { 1 } \end{array} \right] \left\{ \begin{array} { c } { u_1^{(1)} } \\ { u_2^{(1)} } \\ { u_3^{(1)} } \end{array} \right\} + k _ { 2 } \left[ \begin{array} { c c c } { 0 } & { 0 } & { 0 } \\ { 0 } & { 1 } & { - 1 } \\ { 0 } & { - 1 } & { 1 } \end{array} \right] \left\{ \begin{array} { c } { u_1^{(2)} } \\ { u_2^{(2)} } \\ { u_3^{(2)} } \end{array} \right\} = \left\{ \begin{array} { c } { F_{1x} } \\ { F_{2x} } \\ { F_{3x} } \end{array} \right\}\tag{2.4.5}\]

这里，位移上标仍然表示相应的单元编号。化简式（2.4.5）可得

\[\left[ \begin{array} { c c c } { k _ { 1 } } & { 0 } & { - k _ { 1 } } \\ { 0 } & { k _ { 2 } } & { - k _ { 2 } } \\ { - k _ { 1 } } & { - k _ { 2 } } & { k _ { 1 } + k _ { 2 } } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \\ { u _ { 3 } } \end{array} \right\} = \left\{ \begin{array} { c } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \end{array} \right\}\tag{2.4.6}\]

在式（2.4.6）中，与节点位移相关的单元上标已经去掉，因为 $u_1^{(1)}$ 实际上就是整体位移 $u_1$，$u_2^{(2)}$ 实际上就是整体位移 $u_2$；并且由式（2.3.3）可知，$u_3^{(1)}=u_3^{(2)}=u_3$，即节点 3 在整体组合体中的位移。通过叠加得到的式（2.4.6）与式（2.3.9）完全相同。

式（2.4.2）和式（2.4.3）中的扩展单元刚度矩阵也可以直接相加，从而得到式（2.4.6）所示的结构总刚度矩阵。把各个单元刚度矩阵直接装配成整体结构刚度矩阵，并得到整体刚度方程组的这种可靠方法，称为直接刚度法。它是有限元法中最重要的步骤。

对于这个简单例子，把单元刚度矩阵先扩展再叠加并不困难。但是，当问题包含大量自由度时，把每个单元刚度矩阵都扩展到整体刚度矩阵的阶数会非常繁琐。为了避免这种扩展，可以采用直接刚度法的直接形式，也就是一种快捷装配方法。对于这里的弹簧组合体，可以先按照各单元相关自由度标记每个单元刚度矩阵的行和列：

\[[ k ^ { ( 1 ) } ] = \left[ \begin{array} { c c } { { u _ { 1 } } } & { { u _ { 3 } } } \\ { { k _ { 1 } } } & { { - k _ { 1 } } } \\ { { - k _ { 1 } } } & { { k _ { 1 } } } \end{array} \right] \ u _ { 1 } \qquad [ k ^ { ( 2 ) } ] = \left[ \begin{array} { c c } { { u _ { 3 } } } & { { u _ { 2 } } } \\ { { k _ { 2 } } } & { { - k _ { 2 } } } \\ { { - k _ { 2 } } } & { { k _ { 2 } } } \end{array} \right] \ u _ { 3 }\tag{2.4.7}\]

随后，只需把 $[k^{(1)}]$ 和 $[k^{(2)}]$ 中与某一自由度相关的项，直接加到 $[K]$ 中相同自由度对应的位置即可。例如，$[K]$ 中 $u_1$ 行、$u_1$ 列的项只由单元 1 贡献，因为只有单元 1 含有自由度 $u_1$，所以 $k_{11}=k_1$。而 $[K]$ 中 $u_3$ 行、$u_3$ 列的项由单元 1 和单元 2 共同贡献，因为自由度 $u_3$ 与两个单元都有关，所以 $k_{33}=k_1+k_2$。用同样的思路，可以得到

\[[ K ] = [ \begin{array} { c c c } { { u _ { 1 } } } & { { u _ { 2 } } } & { { u _ { 3 } } } \\ { { k _ { 1 } } } & { { 0 } } & { { - k _ { 1 } } } \\ { { 0 } } & { { k _ { 2 } } } & { { - k _ { 2 } } } \\ { { - k _ { 1 } } } & { { - k _ { 2 } } } & { { k _ { 1 } + k _ { 2 } } } \end{array} ] \begin{array} { l } { { u _ { 1 } } } \\ { { u _ { 2 } } } \\ { { u _ { 3 } } } \end{array}\tag{2.4.8}\]

这里，$[K]$ 中各元素的位置依据整体结构自由度按节点编号递增排列的顺序确定。第 2.5 节将在讨论支座边界条件的同时，给出二弹簧组合体的完整求解过程。

2.5 边界条件

对于图 2-5 所示的弹簧组合体这类结构模型，必须指定边界条件（或支座条件），否则 $[K]$ 将是奇异矩阵；也就是说，$[K]$ 的行列式为零，其逆矩阵不存在。这意味着结构系统是不稳定的。如果没有指定足够的运动约束或支座条件，结构将能够作为刚体自由移动，不能抵抗任何外加载荷。一般来说，使 $[K]$ 非奇异所需的边界条件数量，等于可能的刚体模态数量。

与弹簧组合体相关的边界条件通常和节点位移有关。这些条件分为两类。齐次边界条件较为常见，它出现在某些位置完全不能移动的情形；非齐次边界条件则出现在指定了非零位移值的地方，例如支座沉降。

从边值问题的数学角度看，在常微分方程、偏微分方程，或由泛函一阶变分得到的方程上施加边界条件时，通常会遇到两大类边界条件（见文献 [4,5,8]）。不过，本书作为基础教材，不深入展开这些数学分类。

第一类是主边界条件、本质边界条件或 Dirichlet 边界条件，它指定解在区域边界上必须满足的数值，例如指定边界位移。

第二类是自然边界条件或 Neumann 边界条件，它指定解的导数在区域边界上必须满足的数值。

为了说明这两类位移边界条件，我们仍考虑图 2-5 中弹簧组合体对应的式（2.4.6）。该组合体沿弹簧组合体运动方向只有一个刚体模态。

齐次边界条件

先考虑齐次边界条件。此时，所有边界条件都使某些节点上的位移为零。由于节点 1 固定，所以 $u_1=0$。因此，式（2.4.6）可写为

\[\left[ \begin{array} { c c c } { k _ { 1 } } & { 0 } & { - k _ { 1 } } \\ { 0 } & { k _ { 2 } } & { - k _ { 2 } } \\ { - k _ { 1 } } & { - k _ { 2 } } & { k _ { 1 } + k _ { 2 } } \end{array} \right] \left\{ \begin{array} { c } { 0 } \\ { u _ { 2 } } \\ { u _ { 3 } } \end{array} \right\} = \left\{ \begin{array} { c } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \end{array} \right\}\tag{2.5.1}\]

把式（2.5.1）展开，可得

\[\begin{array} { r l r } & { } & { k _ { 1 } ( 0 ) + ( 0 ) u _ { 2 } - k _ { 1 } u _ { 3 } = F _ { 1 x } } \\ & { } & { 0 ( 0 ) + k _ { 2 } u _ { 2 } - k _ { 2 } u _ { 3 } = F _ { 2 x } } \\ & { } & { - k _ { 1 } ( 0 ) - k _ { 2 } u _ { 2 } + ( k _ { 1 } + k _ { 2 } ) u _ { 3 } = F _ { 3 x } } \end{array}\tag{2.5.2}\]

其中，$F_{1x}$ 是未知反力，而 $F_{2x}$ 和 $F_{3x}$ 是已知外加载荷。

把式（2.5.2）中的第二、第三个方程写成矩阵形式，有

\[{ \left[ \begin{array} { l l } { k _ { 2 } } & { - k _ { 2 } } \\ { - k _ { 2 } } & { k _ { 1 } + k _ { 2 } } \end{array} \right] } { \left\{ \begin{array} { l } { u _ { 2 } } \\ { u _ { 3 } } \end{array} \right\} } = { \left\{ \begin{array} { l } { F _ { 2 x } } \\ { F _ { 3 x } } \end{array} \right\} }\tag{2.5.3}\]

这样，我们实际上已经从 $[K]$ 中分离掉第一行和第一列，并从 $\{d\}$ 与 $\{F\}$ 中分离掉第一行，从而得到式（2.5.3）。

对于齐次边界条件，式（2.5.3）也可以直接从式（2.5.1）得到：删去与零位移自由度对应的行和列即可。这里删去第 1 行和第 1 列，因为 $[K]$ 的第 1 列实际上要乘以 $u_1=0$。但是，$F_{1x}$ 并不一定为零；一旦求出 $u_2$ 和 $u_3$，就可以继续求得该反力。

由式（2.5.3）求解 $u_2$ 和 $u_3$ 后，有

\[\left\{ \begin{array} { c } { u _ { 2 } } \\ { u _ { 3 } } \end{array} \right\} = \left[ { \begin{array} { c c } { k _ { 2 } } & { - k _ { 2 } } \\ { - k _ { 2 } } & { k _ { 1 } + k _ { 2 } } \end{array} } \right] ^ { - 1 } \left\{ \begin{array} { c } { F _ { 2 x } } \\ { F _ { 3 x } } \end{array} \right\}\tag{2.5.4}\]

现在，$u_2$ 和 $u_3$ 可由式（2.5.4）确定。将它们代入式（2.5.2）的第一个方程，可得反力 $F_{1x}$：

\[F _ { 1 x } = - k _ { 1 } u _ { 3 }\tag{2.5.5}\]

若把式（2.5.4）得到的 $u_3$ 代入式（2.5.5），还可以用外加节点力 $F_{2x}$ 和 $F_{3x}$ 表示节点 1 处的未知节点力（也称反力）：

\[F _ { 1 x } = - F _ { 2 x } - F _ { 3 x }\tag{2.5.6}\]

因此，对于所有齐次边界条件，可以从原始方程组中删去与零位移自由度对应的行和列，然后求解未知位移。这个过程对手算很有用。（不过，附录 B.4 给出了更适合计算机实现的联立方程求解方案。）

非齐次边界条件

现在考虑非齐次边界条件。此时，一个或多个指定位移不为零。为简单起见，在式（2.4.6）中令 $u_1=\delta$，其中 $\delta$ 是已知位移（图 2-7）。此时有

\[\left[ \begin{array} { c c c } { k _ { 1 } } & { 0 } & { - k _ { 1 } } \\ { 0 } & { k _ { 2 } } & { - k _ { 2 } } \\ { - k _ { 1 } } & { - k _ { 2 } } & { k _ { 1 } + k _ { 2 } } \end{array} \right] \left\{ \begin{array} { c } { \delta } \\ { u _ { 2 } } \\ { u _ { 3 } } \end{array} \right\} = \left\{ \begin{array} { c } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \end{array} \right\}\tag{2.5.7}\]

把式（2.5.7）展开，可得

\[\begin{array} { r l r } & { } & { k _ { 1 } \delta + 0 u _ { 2 } - k _ { 1 } u _ { 3 } = F _ { 1 x } } \\ & { } & { 0 \delta + k _ { 2 } u _ { 2 } - k _ { 2 } u _ { 3 } = F _ { 2 x } } \\ & { } & { - k _ { 1 } \delta - k _ { 2 } u _ { 2 } + ( k _ { 1 } + k _ { 2 } ) u _ { 3 } = F _ { 3 x } } \end{array}\tag{2.5.8}\]

图 2-7 节点 1 处具有已知位移 $\delta$ 的二弹簧组合体

其中，$F_{1x}$ 现在是由支座移动 $\delta$ 后产生的反力。考虑式（2.5.8）中的第二、第三个方程，因为它们右端的节点力 $F_{2x}$ 和 $F_{3x}$ 已知，可得

\[\begin{array} { r } { 0 \delta + k _ { 2 } u _ { 2 } - k _ { 2 } u _ { 3 } = F _ { 2 x } } \\ { - k _ { 1 } \delta - k _ { 2 } u _ { 2 } + ( k _ { 1 } + k _ { 2 } ) u _ { 3 } = F _ { 3 x } } \end{array}\tag{2.5.9}\]

把已知的 $\delta$ 项移到式（2.5.9）的右端，得到

\[\begin{array} { c } { { k _ { 2 } u _ { 2 } - k _ { 2 } u _ { 3 } = F _ { 2 x } } } \\ { { - k _ { 2 } u _ { 2 } + ( k _ { 1 } + k _ { 2 } ) u _ { 3 } = k _ { 1 } \delta + F _ { 3 x } } } \end{array}\tag{2.5.10}\]

将式（2.5.10）重新写成矩阵形式，有

\[\left[ \begin{array} { c c } { k _ { 2 } } & { - k _ { 2 } } \\ { - k _ { 2 } } & { k _ { 1 } + k _ { 2 } } \end{array} \right] \left\{ \begin{array} { c } { u _ { 2 } } \\ { u _ { 3 } } \end{array} \right\} = \left\{ \begin{array} { c } { F _ { 2 x } } \\ { k _ { 1 } \delta + F _ { 3 x } } \end{array} \right\}\tag{2.5.11}\]

因此，在处理非齐次边界条件时，不能一开始就像齐次条件那样删去式（2.5.7）中与该边界条件对应的第 1 行和第 1 列。式（2.5.11）说明，对应列中的项正在与一个非零数相乘。如果直接删去，就会漏掉式（2.5.11）中的 $k_1\delta$ 项，从而使位移解产生错误。一般而言，对于非齐次边界条件，必须先把与已知位移相关的项移到右端力矩阵中，再求解未知节点位移。这里就是把式（2.5.9）第二个方程中的 $k_1\delta$ 项移到式（2.5.10）第二个方程右端来体现这一点。

接下来也可以像求解式（2.5.3）那样求解式（2.5.11）中的位移。不过，由于这样不会提供新的概念信息，本书不再继续展开。

然而，将位移回代到式（2.5.7）后，反力变为

\[F _ { 1 x } = k _ { 1 } \delta - k _ { 1 } u _ { 3 }\tag{2.5.12}\]

它不同于式（2.5.5）中齐次边界条件下的 $F_{1x}$。

需要注意的是，如果某个节点的位移已知（例如 $u_1=\delta$），那么该节点上与该位移同方向的力 $F_{1x}$ 起初并不知道。必须先求出未知节点位移，再利用式（2.5.7）这样的整体方程来确定它。

到这里，我们可以总结式（2.5.7）中整体刚度矩阵的一些性质。这些性质也适用于有限元法的一般形式。

  1. $[K]$ 是方阵，因为它联系的是相同数量的力和位移。
  2. $[K]$ 是对称矩阵，因为每个单元刚度矩阵都是对称矩阵。如果你熟悉结构力学，就不会觉得这个对称性奇怪。它可以用文献 [3] 和 [4] 中介绍的互等定理来证明。
  3. 在施加足够边界条件之前，$[K]$ 是奇异矩阵，也就是它的行列式为零，因而不存在逆矩阵。必须施加足够的边界条件来消除奇异性，并防止刚体运动。
  4. $[K]$ 的主对角线项总是正的。否则，正的节点力 $F_i$ 可能产生负的位移 $d_i$，这与真实结构的物理行为相反。
  5. $[K]$ 是半正定矩阵，也就是说，对所有非零实向量 $\{x\}$，有 $\{x\}^{T}[K]\{x\}>0$。（关于半正定矩阵的更多内容，可参见附录 A。）

一般而言，指定支座条件可以在数学上通过对整体平衡方程进行分块来处理：

\[\left[ \begin{array} { l l } { [ K _ { 1 1 } ] } & { [ K _ { 1 2 } ] } \\ { [ K _ { 2 1 } ] } & { [ K _ { 2 2 } ] } \end{array} \right] \left\{ \begin{array} { l } { \{ d _ { 1 } \} } \\ { \{ d _ { 2 } \} } \end{array} \right\} = \left\{ \begin{array} { l } { \{ F _ { 1 } \} } \\ { \{ F _ { 2 } \} } \end{array} \right\}\tag{2.5.13}\]

其中，令 $\{d_1\}$ 表示未约束的自由位移，$\{d_2\}$ 表示给定的已知位移。由式（2.5.13）可得

\[[ K _ { 1 1 } ] \{ d _ { 1 } \} = \{ F _ { 1 } \} - [ K _ { 1 2 } ] \{ d _ { 2 } \}\tag{2.5.14}\]

以及

\[\{ F _ { 2 } \} = [ K _ { 2 1 } ] \{ d _ { 1 } \} + [ K _ { 2 2 } ] \{ d _ { 2 } \}\tag{2.5.15}\]

其中，$\{F_1\}$ 是已知节点力，$\{F_2\}$ 是指定移位节点处的未知节点力。先由式（2.5.14）求出 $\{d_1\}$，然后可由式（2.5.15）求出 $\{F_2\}$。在式（2.5.14）中，假定 $[K_{11}]$ 已经不再是奇异矩阵，因此可以用它确定 $\{d_1\}$。

为了说明如何用刚度法求解弹簧组合体，下面给出若干例题。

例题 2.1

对于图 2-8 所示、节点编号任意给定的弹簧组合体，求：（a）整体刚度矩阵；（b）节点 3 和节点 4 的位移；（c）节点 1 和节点 2 的反力；（d）每根弹簧中的力。节点 4 沿 $x$ 方向承受 25 kN 的力。弹簧常数如图中所示。节点 1 和节点 2 固定。

图 2-8 用于求解的弹簧组合体

解

（a）首先利用式（2.2.9）把各单元刚度矩阵表示为

\[\begin{array}{c} [ k ^ { ( 1 ) } ] = \left[ \begin{array} { c c } { 200 } & { -200 } \\ { -200 } & { 200 } \end{array} \right],\qquad [ k ^ { ( 2 ) } ] = \left[ \begin{array} { c c } { 400 } & { -400 } \\ { -400 } & { 400 } \end{array} \right],\qquad [ k ^ { ( 3 ) } ] = \left[ \begin{array} { c c } { 600 } & { -600 } \\ { -600 } & { 600 } \end{array} \right] \end{array}\tag{2.5.16}\]

这里，列上方和行旁边的数字表示每个单元对应的节点自由度。例如，单元 1 与自由度 $u_1$ 和 $u_3$ 相关。同时，每个单元的局部 $x$ 轴都与整体 $x$ 轴重合。

利用叠加思想，也就是直接刚度法，可以得到整体刚度矩阵：

\[[ K ] = [ k ^ { ( 1 ) } ] + [ k ^ { ( 2 ) } ] + [ k ^ { ( 3 ) } ]\]

也就是

\[[ K ] = \left[ \begin{array} { c c c c } { 200 } & { 0 } & { -200 } & { 0 } \\ { 0 } & { 600 } & { 0 } & { -600 } \\ { -200 } & { 0 } & { 200 + 400 } & { -400 } \\ { 0 } & { -600 } & { -400 } & { 400 + 600 } \end{array} \right]\tag{2.5.17}\]

（b）式（2.5.17）给出的整体刚度矩阵把整体力与整体位移联系起来：

\[\left\{ { \begin{array} { c } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \\ { F _ { 4 x } } \end{array} } \right\} = \left[ { \begin{array} { c c c c } { 200 } & { 0 } & { -200 } & { 0 } \\ { 0 } & { 600 } & { 0 } & { -600 } \\ { -200 } & { 0 } & { 600 } & { -400 } \\ { 0 } & { -600 } & { -400 } & { 1000 } \end{array} } \right] \left\{ { \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \\ { u _ { 3 } } \\ { u _ { 4 } } \end{array} } \right\}\tag{2.5.18}\]

对式（2.5.18）施加齐次边界条件 $u_1=0$ 和 $u_2=0$，代入已知外力，并按零位移边界条件删去相应方程与刚度矩阵行列，可得

\[\left\{ \begin{array} { c } { 0 } \\ { 25000 } \end{array} \right\} = \left[ \begin{array} { c c } { 600 } & { -400 } \\ { -400 } & { 1000 } \end{array} \right] \left\{ \begin{array} { c } { u _ { 3 } } \\ { u _ { 4 } } \end{array} \right\}\tag{2.5.19}\]

解式（2.5.19），得到整体节点位移

\[u _ { 3 } = { \frac { 250 } { 11 } }~\mathrm{mm}, \qquad u _ { 4 } = { \frac { 375 } { 11 } }~\mathrm{mm}\tag{2.5.20}\]

（c）为了求整体节点力，也就是包括节点 1 和节点 2 处反力在内的力，把式（2.5.20）以及边界条件 $u_1=0$、$u_2=0$ 回代入式（2.5.18）。得到

\[\left\{ { \begin{array} { c } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \\ { F _ { 4 x } } \end{array} } \right\} = \left[ { \begin{array} { c c c c } { 200 } & { 0 } & { -200 } & { 0 } \\ { 0 } & { 600 } & { 0 } & { -600 } \\ { -200 } & { 0 } & { 600 } & { -400 } \\ { 0 } & { -600 } & { -400 } & { 1000 } \end{array} } \right] \left\{ { \begin{array} { c } { 0 } \\ { 0 } \\ { \frac { 250 } { 11 } } \\ { \frac { 375 } { 11 } } \end{array} } \right\}\tag{2.5.21}\]

矩阵相乘并化简后，可得各节点处的力

\[\begin{array} { r l } F _ { 1 x } = { \frac { -50000 } { 11 } }~\mathrm{N}, \qquad & F _ { 2 x } = { \frac { -225000 } { 11 } }~\mathrm{N}, \\ F _ { 3 x } = 0, \qquad & F _ { 4 x } = { \frac { 275000 } { 11 } }~\mathrm{N} \end{array}\tag{2.5.22}\]

由这些结果可以看出，反力 $F_{1x}$ 与 $F_{2x}$ 的合力，大小等于、方向相反于外加力 $F_{4x}$。这验证了整个弹簧组合体的平衡。

（d）下面利用局部单元方程式（2.2.8）求各单元中的力。

单元 1

\[\left\{ \begin{array} { c } { f _ { 1 x } ^ { ( 1 ) } } \\ { f _ { 3 x } ^ { ( 1 ) } } \end{array} \right\} = \left[ \begin{array} { c c } { 200 } & { -200 } \\ { -200 } & { 200 } \end{array} \right] \left\{ \begin{array} { c } { 0 } \\ { \frac { 250 } { 11 } } \end{array} \right\}\tag{2.5.23}\]

化简式（2.5.23），得到

\[f _ { 1 x } ^ { ( 1 ) } = \frac { -50000 } { 11 }~\mathrm{N}, \qquad f _ { 3 x } ^ { ( 1 ) } = \frac { 50000 } { 11 }~\mathrm{N}\tag{2.5.24}\]

图 2-9（a）给出了弹簧单元 1 的自由体图。该弹簧承受式（2.5.24）给出的拉力。同时，$f_{1x}^{(1)}$ 等于式（2.5.22）给出的反力 $F_{1x}$。图 2-9（b）所示的节点 1 自由体图也显示了这一点。

图 2-9 （a）单元 1 的自由体图；（b）节点 1 的自由体图

单元 2

\[\left\{ \begin{array} { c } { f _ { 3 x } ^ { ( 2 ) } } \\ { f _ { 4 x } ^ { ( 2 ) } } \end{array} \right\} = \left[ \begin{array} { c c } { 400 } & { -400 } \\ { -400 } & { 400 } \end{array} \right] \left\{ \begin{array} { c } { \frac { 250 } { 11 } } \\ { \frac { 375 } { 11 } } \end{array} \right\}\tag{2.5.25}\]

化简式（2.5.25），得到

\[f _ { 3 x } ^ { ( 2 ) } = \frac { -50000 } { 11 }~\mathrm{N}, \qquad f _ { 4 x } ^ { ( 2 ) } = \frac { 50000 } { 11 }~\mathrm{N}\tag{2.5.26}\]

图 2-10 给出了弹簧单元 2 的自由体图。该弹簧承受式（2.5.26）给出的拉力。

图 2-10 单元 2 的自由体图

单元 3

\[\left\{ \begin{array} { c } { f _ { 4 x } ^ { ( 3 ) } } \\ { f _ { 2 x } ^ { ( 3 ) } } \end{array} \right\} = \left[ \begin{array} { c c } { 600 } & { -600 } \\ { -600 } & { 600 } \end{array} \right] \left\{ \begin{array} { c } { \frac { 375 } { 11 } } \\ { 0 } \end{array} \right\}\tag{2.5.27}\]

化简式（2.5.27），得到

\[f _ { 4 x } ^ { ( 3 ) } = \frac { 225000 } { 11 }~\mathrm{N}, \qquad f _ { 2 x } ^ { ( 3 ) } = \frac { -225000 } { 11 }~\mathrm{N}\tag{2.5.28}\]

图 2-11 （a）单元 3 的自由体图；（b）节点 2 的自由体图

图 2-11（a）给出了弹簧单元 3 的自由体图。该弹簧承受式（2.5.28）给出的压力。同时，$f_{2x}^{(3)}$ 等于式（2.5.22）给出的反力 $F_{2x}$。图 2-11（b）中的节点 2 自由体图也说明了这一点。

例题 2.2

对于图 2-12 所示的弹簧组合体，求：（a）整体刚度矩阵；（b）节点 2 至节点 4 的位移；（c）整体节点力；（d）局部单元力。节点 1 固定，而节点 5 具有给定的已知位移 $\delta=20.0~\mathrm{mm}$。所有弹簧常数均为 $k=200~\mathrm{kN/m}$。

图 2-12 用于求解的弹簧组合体

解

（a）利用式（2.2.9），可把每个单元刚度矩阵写为

\[[ k ^ { ( 1 ) } ] = [ k ^ { ( 2 ) } ] = [ k ^ { ( 3 ) } ] = [ k ^ { ( 4 ) } ] = \left[ \begin{array} { r r } { { 200 } } & { { -200 } } \\ { { -200 } } & { { 200 } } \end{array} \right]\tag{2.5.29}\]

再次使用叠加法，可得整体刚度矩阵

\[[ K ] = \left[ \begin{array} { c c c c c } { { 200 } } & { { -200 } } & { { 0 } } & { { 0 } } & { { 0 } } \\ { { -200 } } & { { 400 } } & { { -200 } } & { { 0 } } & { { 0 } } \\ { { 0 } } & { { -200 } } & { { 400 } } & { { -200 } } & { { 0 } } \\ { { 0 } } & { { 0 } } & { { -200 } } & { { 400 } } & { { -200 } } \\ { { 0 } } & { { 0 } } & { { 0 } } & { { -200 } } & { { 200 } } \end{array} \right] \frac { \mathrm{kN} } { \mathrm{m} }\tag{2.5.30}\]

（b）式（2.5.30）中的整体刚度矩阵把整体力与整体位移联系为

\[\left\{ { \begin{array} { c } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \\ { F _ { 4 x } } \\ { F _ { 5 x } } \end{array} } \right\} = \left[ { \begin{array} { c c c c c } { 200 } & { -200 } & { 0 } & { 0 } & { 0 } \\ { -200 } & { 400 } & { -200 } & { 0 } & { 0 } \\ { 0 } & { -200 } & { 400 } & { -200 } & { 0 } \\ { 0 } & { 0 } & { -200 } & { 400 } & { -200 } \\ { 0 } & { 0 } & { 0 } & { -200 } & { 200 } \end{array} } \right] \left\{ { \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \\ { u _ { 3 } } \\ { u _ { 4 } } \\ { u _ { 5 } } \end{array} } \right\}\tag{2.5.31}\]

施加边界条件 $u_1=0$ 和 $u_5=20~\mathrm{mm}=0.02~\mathrm{m}$，并代入已知整体力 $F_{2x}=0$、$F_{3x}=0$、$F_{4x}=0$。对式（2.5.31）中与这些边界条件对应的第 1 个和第 5 个方程进行分块，可得

\[\left\{ \begin{array} { c } { { 0 } } \\ { { 0 } } \\ { { 0 } } \end{array} \right\} = \left[ \begin{array} { c c c c c } { { -200 } } & { { 400 } } & { { -200 } } & { { 0 } } & { { 0 } } \\ { { 0 } } & { { -200 } } & { { 400 } } & { { -200 } } & { { 0 } } \\ { { 0 } } & { { 0 } } & { { -200 } } & { { 400 } } & { { -200 } } \end{array} \right] \left\{ \begin{array} { c } { { 0 } } \\ { { u _ { 2 } } } \\ { { u _ { 3 } } } \\ { { u _ { 4 } } } \\ { { 0.02~\mathrm{m} } } \end{array} \right\}\tag{2.5.32}\]

把与已知位移 $0.02~\mathrm{m}$ 相乘的刚度项移到左端，可把式（2.5.32）重写为

\[\left\{ \begin{array} { c } { { 0 } } \\ { { 0 } } \\ { { 4~\mathrm{kN} } } \end{array} \right\} = \left[ \begin{array} { c c c } { { 400 } } & { { -200 } } & { { 0 } } \\ { { -200 } } & { { 400 } } & { { -200 } } \\ { { 0 } } & { { -200 } } & { { 400 } } \end{array} \right] \left\{ \begin{array} { c } { { u _ { 2 } } } \\ { { u _ { 3 } } } \\ { { u _ { 4 } } } \end{array} \right\}\tag{2.5.33}\]

解式（2.5.33），得到

\[u _ { 2 } = 0.005~\mathrm{m}, \qquad u _ { 3 } = 0.01~\mathrm{m}, \qquad u _ { 4 } = 0.015~\mathrm{m}\tag{2.5.34}\]

（c）把边界条件位移和式（2.5.34）回代到式（2.5.31），即可得到整体节点力：

\[\begin{array} { r l } & { F _ { 1 x } = ( -200 ) ( 0.005 ) = -1.0~\mathrm{kN} } \\ & { F _ { 2 x } = ( 400 ) ( 0.005 ) - ( 200 ) ( 0.01 ) = 0 } \\ & { F _ { 3 x } = ( -200 ) ( 0.005 ) + ( 400 ) ( 0.01 ) - ( 200 ) ( 0.015 ) = 0 } \\ & { F _ { 4 x } = ( -200 ) ( 0.01 ) + ( 400 ) ( 0.015 ) - ( 200 ) ( 0.02 ) = 0 } \\ & { F _ { 5 x } = ( -200 ) ( 0.015 ) + ( 200 ) ( 0.02 ) = 1.0~\mathrm{kN} } \end{array}\tag{2.5.35}\]

式（2.5.35）给出的结果表明，反力 $F_{1x}$ 与使节点 5 产生 $\delta=20.0~\mathrm{mm}$ 位移所需的节点力 $F_{5x}$ 大小相等、方向相反。这验证了整个弹簧组合体的平衡。

要记住：如果某个节点在某个方向上的位移已知（本例中为 $u_5=20~\mathrm{mm}$），那么同一节点、同一方向上的力 $F_{5x}$ 起初并不知道。它要在求出未知节点位移之后再确定。

（d）接着，利用局部单元方程式（2.2.10）求各单元中的力。

单元 1

\[\left\{ \begin{array} { c } { f _ { 1 x } ^ { ( 1 ) } } \\ { f _ { 2 x } ^ { ( 1 ) } } \end{array} \right\} = \left[ \begin{array} { c c } { 200 } & { -200 } \\ { -200 } & { 200 } \end{array} \right] \left\{ \begin{array} { c } { 0 } \\ { 0.005 } \end{array} \right\}\tag{2.5.36}\]

化简式（2.5.36），得到

\[f _ { 1 x } ^ { ( 1 ) } = -1.0~\mathrm{kN}, \qquad f _ { 2 x } ^ { ( 1 ) } = 1.0~\mathrm{kN}\tag{2.5.37}\]

单元 2

\[\left\{ \begin{array} { c } { f _ { 2 x } ^ { ( 2 ) } } \\ { f _ { 3 x } ^ { ( 2 ) } } \end{array} \right\} = \left[ \begin{array} { c c } { 200 } & { -200 } \\ { -200 } & { 200 } \end{array} \right] \left\{ \begin{array} { c } { 0.005 } \\ { 0.01 } \end{array} \right\}\tag{2.5.38}\]

化简式（2.5.38），得到

\[f _ { 2 x } ^ { ( 2 ) } = -1.0~\mathrm{kN}, \qquad f _ { 3 x } ^ { ( 2 ) } = 1.0~\mathrm{kN}\tag{2.5.39}\]

单元 3

\[\left\{ \begin{array} { c } { f _ { 3 x } ^ { ( 3 ) } } \\ { f _ { 4 x } ^ { ( 3 ) } } \end{array} \right\} = \left[ \begin{array} { c c } { 200 } & { -200 } \\ { -200 } & { 200 } \end{array} \right] \left\{ \begin{array} { c } { 0.01 } \\ { 0.015 } \end{array} \right\}\tag{2.5.40}\]

化简式（2.5.40），得到

\[f _ { 3 x } ^ { ( 3 ) } = -1.0~\mathrm{kN}, \qquad f _ { 4 x } ^ { ( 3 ) } = 1.0~\mathrm{kN}\tag{2.5.41}\]

单元 4

\[\left\{ \begin{array} { c } { f _ { 4 x } ^ { ( 4 ) } } \\ { f _ { 5 x } ^ { ( 4 ) } } \end{array} \right\} = \left[ \begin{array} { c c } { 200 } & { -200 } \\ { -200 } & { 200 } \end{array} \right] \left\{ \begin{array} { c } { 0.015 } \\ { 0.02 } \end{array} \right\}\tag{2.5.42}\]

化简式（2.5.42），得到

\[f _ { 4 x } ^ { ( 4 ) } = -1.0~\mathrm{kN}, \qquad f _ { 5 x } ^ { ( 4 ) } = 1.0~\mathrm{kN}\tag{2.5.43}\]

你可以画出每个节点和每个单元的自由体图，并利用式（2.5.35）至式（2.5.43）的结果，验证节点平衡和单元平衡。

为了回顾本章介绍的主要概念，下面再求解一个例题。

例题 2.3

（a）对图 2-13 所示的线弹性弹簧系统，利用第 2.3 节中的思想，写出边界条件、类似式（2.3.3）的协调或连续条件，以及类似式（2.3.4）至式（2.3.6）的节点平衡条件。然后建立整体刚度矩阵和整体方程，用来求解未知整体位移和力。各单元的弹簧常数分别为 $k_1$、$k_2$ 和 $k_3$；$P$ 是作用在节点 2 上的外力。

（b）用直接刚度法建立与（a）相同的整体刚度矩阵和整体方程。

图 2-13 用于求解的弹簧组合体

解

（a）边界条件为

\[u _ { 1 } = 0 \qquad u _ { 3 } = 0 \qquad u _ { 4 } = 0\tag{2.5.44}\]

节点 2 处的协调条件为

\[u _ { 2 } ^ { ( 1 ) } = u _ { 2 } ^ { ( 2 ) } = u _ { 2 } ^ { ( 3 ) } = u _ { 2 }\tag{2.5.45}\]

节点平衡条件为

\[\begin{array} { r l } & { F _ { 1 x } = f _ { 1 x } ^ { ( 1 ) } } \\ & { P = f _ { 2 x } ^ { ( 1 ) } + f _ { 2 x } ^ { ( 2 ) } + f _ { 2 x } ^ { ( 3 ) } } \\ & { F _ { 3 x } = f _ { 3 x } ^ { ( 2 ) } } \\ & { F _ { 4 x } = f _ { 4 x } ^ { ( 3 ) } } \end{array}\tag{2.5.46}\]

式（2.5.46）采用了图 2-2 给出的正单元节点力符号约定。图 2-14 给出了各单元和各节点的力自由体图。

图 2-14 图 2-13 中弹簧组合体的单元和节点自由体图

将局部刚度矩阵用于每个单元，并结合式（2.5.45）的协调条件，可得到总平衡方程或整体平衡方程：

\[\begin{array} { r l } & { F _ { 1 x } = k _ { 1 } u _ { 1 } - k _ { 1 } u _ { 2 } } \\ & { P = - k _ { 1 } u _ { 1 } + k _ { 1 } u _ { 2 } + k _ { 2 } u _ { 2 } - k _ { 2 } u _ { 3 } + k _ { 3 } u _ { 2 } - k _ { 3 } u _ { 4 } } \\ & { F _ { 3 x } = - k _ { 2 } u _ { 2 } + k _ { 2 } u _ { 3 } } \\ & { F _ { 4 x } = - k _ { 3 } u _ { 2 } + k _ { 3 } u _ { 4 } } \end{array}\tag{2.5.47}\]

把式（2.5.47）写成矩阵形式，有

\[\left\{ \begin{array} { c } { F _ { 1 x } } \\ { P } \\ { F _ { 3 x } } \\ { F _ { 4 x } } \end{array} \right\} = \left[ \begin{array} { c c c c } { k _ { 1 } } & { - k _ { 1 } } & { 0 } & { 0 } \\ { - k _ { 1 } } & { k _ { 1 } + k _ { 2 } + k _ { 3 } } & { - k _ { 2 } } & { - k _ { 3 } } \\ { 0 } & { - k _ { 2 } } & { k _ { 2 } } & { 0 } \\ { 0 } & { - k _ { 3 } } & { 0 } & { k _ { 3 } } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \\ { u _ { 3 } } \\ { u _ { 4 } } \end{array} \right\}\tag{2.5.48}\]

因此，整体刚度矩阵就是式（2.5.48）右端的方形对称矩阵。利用边界条件式（2.5.44），并考虑式（2.5.47）或式（2.5.48）中的第二个方程，可解得

\[u _ { 2 } = { \frac { P } { k _ { 1 } + k _ { 2 } + k _ { 3 } } }\tag{2.5.49}\]

这个结果也可以按前面介绍的方法得到：删去 $\{F\}$ 和 $\{d\}$ 中与零位移对应的第 1、3、4 行，并删去 $[K]$ 中对应的第 1、3、4 行和列，然后求解 $u_2$。

利用式（2.5.47），可求得整体力

\[F _ { 1 x } = - k _ { 1 } u _ { 2 }, \qquad F _ { 3 x } = - k _ { 2 } u _ { 2 }, \qquad F _ { 4 x } = - k _ { 3 } u _ { 2 }\tag{2.5.50}\]

式（2.5.50）中的力可解释为本例中的整体反力。负号表示这些力指向左侧，也就是与 $x$ 轴正方向相反。

（b）下面用直接刚度法建立整体刚度矩阵。首先，根据弹簧单元刚度矩阵，把每个单元刚度矩阵写为

\[[ k ^ { ( 1 ) } ] = \left[ \begin{array} { c c } { k _ { 1 } } & { - k _ { 1 } } \\ { - k _ { 1 } } & { k _ { 1 } } \end{array} \right]_{(u_1,u_2)}, \quad [ k ^ { ( 2 ) } ] = \left[ \begin{array} { c c } { k _ { 2 } } & { - k _ { 2 } } \\ { - k _ { 2 } } & { k _ { 2 } } \end{array} \right]_{(u_2,u_3)}, \quad [ k ^ { ( 3 ) } ] = \left[ \begin{array} { c c } { k _ { 3 } } & { - k _ { 3 } } \\ { - k _ { 3 } } & { k _ { 3 } } \end{array} \right]_{(u_2,u_4)}\tag{2.5.51}\]

下标中的自由度表示各单元矩阵的行和列对应的整体自由度。按照第 2.4 节介绍的直接刚度法，把各单元刚度矩阵中的项加到整体刚度矩阵中相同自由度对应的位置，可得

\[[ K ] = \left[ \begin{array} { c c c c } { k _ { 1 } } & { - k _ { 1 } } & { 0 } & { 0 } \\ { - k _ { 1 } } & { k _ { 1 } + k _ { 2 } + k _ { 3 } } & { - k _ { 2 } } & { - k _ { 3 } } \\ { 0 } & { - k _ { 2 } } & { k _ { 2 } } & { 0 } \\ { 0 } & { - k _ { 3 } } & { 0 } & { k _ { 3 } } \end{array} \right]\tag{2.5.52}\]

可以看到，每个单元刚度矩阵 $[k]$ 都被加到了整体矩阵 $[K]$ 中与其自由度相同的位置。例如，单元 3 与自由度 $u_2$ 和 $u_4$ 有关，因此它对 $[K]$ 的贡献出现在 $2$-$2$、$2$-$4$、$4$-$2$ 和 $4$-$4$ 位置，这些位置正是式（2.5.52）中的 $k_3$ 项。

用直接刚度法装配出整体 $[K]$ 后，便可按通常方式利用一般方程 $\{F\}=[K]\{d\}$ 建立整体方程。这些方程已经由式（2.5.48）给出，因此不再重复。

罚函数法处理给定位移

另一种处理强制边界条件的方法称为罚函数法。它既可处理齐次（零）指定自由度，也可处理非齐次（非零）指定自由度，而且很容易在计算机程序中实现。

考虑图 2-15 所示的简单弹簧组合体，它受到外力 $F_{1x}$ 和 $F_{2x}$ 作用。假设节点 1 的水平位移被强制规定为 $u_1=\delta$。

如图 2-16 所示，在节点位移 $u_1=\delta$ 的方向上，向组合体中再加入一个具有很大刚度 $k_b$ 的弹簧，这个弹簧常称为边界单元。该弹簧刚度的数量级通常取为最大 $k_{ii}$ 项的大约 $10^6$ 倍。

图 2-15 用于说明罚函数法的弹簧组合体

图 2-16 在节点 1 处加入边界弹簧单元后的弹簧组合体

现在，在 $u_1$ 方向上加入力 $k_b\delta$，然后按通常方法求解该问题。单元刚度矩阵为

\[[ k ^ { ( 1 ) } ] = \left[ \begin{array} { c c } { k _ { 1 } } & { - k _ { 1 } } \\ { - k _ { 1 } } & { k _ { 1 } } \end{array} \right], \qquad [ k ^ { ( 2 ) } ] = \left[ \begin{array} { c c } { k _ { 2 } } & { - k _ { 2 } } \\ { - k _ { 2 } } & { k _ { 2 } } \end{array} \right]\tag{2.5.53}\]

用直接刚度法装配单元刚度矩阵，可得整体刚度矩阵

\[[ K ] = \left[ \begin{array} { c c c } { k _ { 1 } + k _ { b } } & { - k _ { 1 } } & { 0 } \\ { - k _ { 1 } } & { k _ { 1 } + k _ { 2 } } & { - k _ { 2 } } \\ { 0 } & { - k _ { 2 } } & { k _ { 2 } } \end{array} \right]\tag{2.5.54}\]

装配整体方程 $\{F\}=[K]\{d\}$，并施加边界条件 $u_3=0$，有

\[\left\{ \begin{array} { c } { F _ { 1 x } + k _ { b } \delta } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \end{array} \right\} = \left[ \begin{array} { c c c } { k _ { 1 } + k _ { b } } & { - k _ { 1 } } & { 0 } \\ { - k _ { 1 } } & { k _ { 1 } + k _ { 2 } } & { - k _ { 2 } } \\ { 0 } & { - k _ { 2 } } & { k _ { 2 } } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \\ { u _ { 3 } = 0 } \end{array} \right\}\tag{2.5.55}\]

求解式（2.5.55）中的第一个和第二个方程，可得

\[u _ { 1 } = \frac { F _ { 2 x } - ( k _ { 1 } + k _ { 2 } ) u _ { 2 } } { - k _ { 1 } }\tag{2.5.56}\]

以及

\[u _ { 2 } = \frac { ( k _ { 1 } + k _ { b } ) F _ { 2 x } + F _ { 1 x } k _ { 1 } + k _ { b } \delta k _ { 1 } } { k _ { b } k _ { 1 } + k _ { b } k _ { 2 } + k _ { 1 } k _ { 2 } }\tag{2.5.57}\]

当 $k_b$ 趋于无穷大时，式（2.5.57）化简为

\[u _ { 2 } = \frac { F _ { 2 x } + \delta k _ { 1 } } { k _ { 1 } + k _ { 2 } }\tag{2.5.58}\]

而式（2.5.56）化简为

\[u _ { 1 } = \delta\tag{2.5.59}\]

这些结果与一开始直接令 $u_1=\delta$ 所得到的结果一致。

使用罚函数法时，具有极大刚度的边界单元应当与某个自由度平行，就像前面的例子一样。如果 $k_b$ 是倾斜的，或放置在结构内部，它就会同时影响整体刚度矩阵 $[K]$ 的对角项和非对角项。这种情况可能导致求解 $\{F\}=[K]\{d\}$ 时出现数值困难。为避免这种情况，可以把倾斜支座处的位移转换到局部坐标中处理，这将在第 3.9 节说明。

2.6 用势能法推导弹簧单元方程

推导单元方程和单元刚度矩阵时，常用的另一种方法基于最小势能原理。（该原理在结构力学中的应用详见文献 [4]。）与第 2.2 节中使用节点平衡、单元平衡以及单元应力-应变关系的方法相比，这种方法更具一般性。因此，在处理复杂单元时，例如平面应力/应变单元、轴对称应力单元、板弯曲单元和三维实体应力单元，最小势能原理更便于用来确定单元方程。

需要再次说明：虚功原理（附录 E）适用于任意材料行为，而最小势能原理只适用于弹性材料。不过，对于线弹性材料，两种原理会导出相同的单元方程，而本书讨论的也正是线弹性材料。此外，最小势能原理属于变分方法这一大类（虚功原理也是如此）。这类方法还会导向与势能类似的其他变分函数或泛函，可用于其他类型问题，尤其是非结构类问题。这些问题通常称为场问题，包括杆的扭转、热传导（第 13 章）、流体流动（第 14 章）和电势问题（第 14 章）等。

还有一些问题，其变分形式并不容易清楚定义，此时可以用加权余量法来建立方程。第 3.12 节将介绍 Galerkin 法；第 3.13 节将介绍配点法、最小二乘法和子域加权余量法。在第 3.13 节中，还会用这四种余量法分别求解一个一维杆问题，并把各自结果与精确解比较。（关于加权余量法的更多内容，也可参见文献 [5-7]。）

本节介绍如何用最小势能原理推导弹簧单元方程。我们先把这个概念用于最简单的单元，希望这样能让读者在后续章节把它用于更复杂单元时更加自然。

（a）

（b） 图 2-17 （a）受逐渐增大外力 $F$ 作用的弹簧；（b）线性弹簧的力-变形曲线

结构的总势能 $\pi_p$ 用位移表示。在有限元列式中，这些位移通常是节点位移，因此可写为 $\pi_p=\pi_p(d_1,d_2,\ldots,d_n)$。当相对于这些位移使 $\pi_p$ 取最小值时，就会得到平衡方程。对于弹簧单元，我们将说明这样得到的节点平衡方程与第 2.2 节中推导出的 $[k]\{d\}=\{f\}$ 相同。

最小势能原理可以表述如下：

在一个物体所有几何上可能的形状中，真正满足稳定平衡的形状，对应于总势能的最小值。

为了解释这个原理，必须先说明势能和函数驻值这两个概念。下面依次讨论。

总势能定义为内应变能 $U$ 与外力势能 $\Omega$ 之和，即

\[\pi _ { p } = U + \Omega\tag{2.6.1}\]

应变能是结构内部力（或应力）通过变形（应变）做功的能力；$\Omega$ 则是体力、面力和外加节点力等外力通过结构变形做功的能力。

为了理解内应变能，先讨论外力功。本节只考虑外加节点力所做的外功。在第 3 章第 10 节中，将讨论体力（通常为自重）和面力（分布力）所做的功。对于线弹性构件，外功可以这样理解：以图 2-17（a）所示的弹性弹簧为例，将力 $F$ 的大小逐渐增大，直到某个最大值 $F_{\max}$，且该最大值低于会使弹簧产生永久变形的水平。当力达到最大值时，对应的最大变形为 $x_{\max}$，如图 2-17（b）所示。外功等于图 2-17（b）中力-变形曲线下的面积，该直线斜率等于弹簧常数 $k$。根据基本力学，外功 $W_e$ 可写为力向量 $\mathbf{F}$ 与微小位移 $d\mathbf{x}$ 点积的积分：

\[W _ { e } = \int { \boldsymbol { F } \cdot d \boldsymbol { x } } = \int _ { 0 } ^ { x _ { \max } } F _ { \max } \left( \frac { x } { x _ { \max } } \right) d x = F _ { \max } x _ { \max } / 2\tag{2.6.2}\]

其中，式（2.6.2）中的 $F$ 为

\[F = F _ { \max } ( x / x _ { \max } )\tag{2.6.3}\]

在式（2.6.2）中，写出右端第二个积分时，$F$ 与 $dx$ 的方向相同。

根据机械能守恒原理，由外力 $F$ 所做的外功会转化为弹簧的内应变能 $U$。因此，弹簧中的应变能为

\[W _ { e } = U = F _ { \max } x _ { \max } / 2\tag{2.6.4}\]

当力逐渐减小到零时，弹簧恢复到原来的未变形状态。储存在变形弹性弹簧中、并在卸载时释放出来的这部分能量，称为内应变能，简称应变能。此外，

\[F _ { \max } = k x _ { \max }\tag{2.6.5}\]

将式（2.6.5）代入式（2.6.4），可把应变能表示为

\[U = k x _ { \max } ^ { 2 } / 2\tag{2.6.6}\]

外力势能与外功符号相反，因为外力做功时外力势能减少。因此有

\[\Omega = - F _ { \max } x _ { \max }\tag{2.6.7}\]

将式（2.6.6）和式（2.6.7）代入式（2.6.1），得到总势能

\[\pi _ { p } = \frac { 1 } { 2 } k x _ { \max } ^ { 2 } - F _ { \max } x _ { \max }\tag{2.6.8}\]

一般地，对弹簧任意变形 $x$ 及其对应力 $F$，可用 $x$ 替代 $x_{\max}$，用 $F$ 替代 $F_{\max}$，于是

\[U ( x ) = k x ^ { 2 } / 2\tag{2.6.8a}\]
\[\Omega ( x ) = - F x\tag{2.6.8b}\]

把式（2.6.8a）和式（2.6.8b）代入式（2.6.1），可得总势能表达式

图 2-18 函数的驻值

\[\pi _ { p } ( x ) = { \frac { 1 } { 2 } } k x ^ { 2 } - F x\tag{2.6.9}\]

图 2-18 展示了函数 $G$ 的驻值概念，这个概念用于最小势能原理的定义。这里 $G$ 表示为变量 $x$ 的函数。驻值可以是最大值、最小值，也可以是中性点。为了找到使 $G(x)$ 取得驻值的 $x$，可对 $G$ 关于 $x$ 求导并令其等于零：

\[{ \frac { d G } { d x } } = 0\tag{2.6.10}\]

后面会用类似过程，把 $G$ 换成 $\pi_p$，把 $x$ 换成离散值，也就是节点位移 $d_i$。如果掌握变分法（见文献 [8]），可用 $\pi_p$ 的一阶变分来使 $\pi_p$ 最小，其中 $\delta\pi_p$ 表示势能的任意变化或变分。不过，本书不展开变分法细节，而是说明可以用熟悉的微积分来完成 $\pi_p$ 的最小化。为了应用最小势能原理，也就是使 $\pi_p$ 最小，需要对作为节点位移 $d_i$ 函数的 $\pi_p$ 取变分，一般写为

\[\delta \pi _ { p } = \frac { \partial \pi _ { p } } { \partial d _ { 1 } } \delta d _ { 1 } + \frac { \partial \pi _ { p } } { \partial d _ { 2 } } \delta d _ { 2 } + \dots + \frac { \partial \pi _ { p } } { \partial d _ { n } } \delta d _ { n }\tag{2.6.11}\]

该原理指出：当 $d_i$ 描述的结构状态满足 $\delta\pi_p=0$，并且这种条件对从平衡状态出发的任意容许位移变分 $\delta d_i$ 都成立时，结构处于平衡状态。所谓容许变分，是指变分后的位移场仍然满足边界条件和单元间连续性。图 2-19（a）展示了一个具有给定位移 $u_1$ 和 $u_2$ 的弹簧，其假想真实轴向位移函数以及一个容许位移函数。图 2-19（b）展示了不可容许函数：一种在端点 1 和端点 2 之间斜率不连续，另一种不满足右端边界条件 $u(L)=u_2$。这里 $\delta u$ 表示 $u$ 的变分。在一般有限元列式中，$\delta u$ 将由 $\delta d_i$ 取代。这意味着任意一个 $\delta d_i$ 都可能非零。因此，为了满足 $\delta\pi_p=0$，所有与 $\delta d_i$ 相乘的系数必须分别为零。于是有

图 2-19 （a）真实位移函数和容许位移函数；（b）不可容许位移函数

\[\frac { \partial \pi _ { p } } { \partial d _ { i } } = 0 \qquad ( i = 1 , 2 , 3 , \ldots , n ) \qquad \mathrm { or } \qquad \frac { \partial \pi _ { p } } { \partial \{ d \} } = 0\tag{2.6.12}\]

其中，需要求解 $n$ 个方程以得到定义结构静力平衡状态的 $n$ 个 $d_i$。式（2.6.12）表明，在本书范围内，可以把 $\pi_p$ 的变分理解为一种紧凑记号，它等价于对 $\pi_p$ 关于未知节点位移求导。对处于平衡状态的线弹性材料而言，$\pi_p$ 确实取最小值，这一点可参见文献 [4]。

在讨论弹簧单元方程的建立之前，先通过例题 2.4 分析一个受外力作用的单自由度弹簧，以说明最小势能原理的含义。该例将表明，弹簧的平衡位置对应于最小势能。

图 2-20 受力弹簧及其载荷-位移曲线

例题 2.4

对于图 2-20 所示、承受 $5000~\mathrm{N}$ 外力的线弹性弹簧，计算不同位移值下的势能，并说明最小势能也对应于弹簧的平衡位置。

解

总势能为

\[\pi _ { p } = U + \Omega\]

其中

\[U = \frac 1 2 ( k x ) x, \qquad \Omega = - F x\]

下面用常规数学方法说明 $\pi_p$ 的最小化。根据式（2.6.11）和式（2.6.12），对 $\pi_p$ 关于 $x$ 取变分；等价地，也可以对 $\pi_p$ 关于 $x$ 求导，因为这里 $\pi_p$ 只是单个位移 $x$ 的函数：

\[\delta \pi _ { p } = \frac { \partial \pi _ { p } } { \partial x } \delta x = 0\]

由于 $\delta x$ 是任意的，并且不一定为零，因此

\[\frac { \partial \pi _ { p } } { \partial x } = 0\]

代入前面得到的 $\pi_p$ 表达式，有

\[\frac { \partial \pi _ { p } } { \partial x } = 125 x - 5000 = 0\]

因此

\[x = 40~\mathrm{mm}\]

把这个 $x$ 值回代入 $\pi_p$，得到

\[\pi _ { p } = 62.5 ( 40 ) ^ { 2 } - 5000 ( 40 ) = -100000~\mathrm{N\cdot mm}\]

这正对应于表 2-1 中用搜索法得到的最小势能。这里 $U=\frac12(kx)x$ 是应变能，也就是图 2-20 所示载荷-位移曲线下的面积；$\Omega=-Fx$ 是载荷 $F$ 的势能。对给定的 $F$ 和 $k$，有

\[\pi _ { p } = { \frac { 1 } { 2 } } (125)x^2 - 5000x = 62.5x^2 - 5000x\]

现在对若干弹簧变形 $x$ 搜索 $\pi_p$ 的最小值，结果列于表 2-1。图 2-21 给出了 $\pi_p$ 随 $x$ 的变化曲线。

表 2-1 不同弹簧变形下的总势能

变形 $x$，mm总势能 $\pi_p$，N·m
-80800
-60525
-40300
-20125
00
20-75
40-100
60-75
800
100125

图 2-21 势能随弹簧变形的变化

从图和表可以看出，当 $x=40~\mathrm{mm}$ 时，$\pi_p$ 取得最小值。这个变形位置也对应于平衡位置，因为

\[\frac { \partial \pi_p } { \partial x } = 125(40)-5000=0\]

下面用最小势能原理推导弹簧单元方程和刚度矩阵。考虑图 2-22 所示、承受节点力的线性弹簧。由式（2.6.9）可知，总势能为

\[\pi _ { p } = \frac { 1 } { 2 } k ( u _ { 2 } - u _ { 1 } ) ^ { 2 } - f _ { 1 x } u _ { 1 } - f _ { 2 x } u _ { 2 }\tag{2.6.13}\]

其中，$u_2-u_1$ 是式（2.6.9）中的弹簧变形。式（2.6.13）右端第一项是弹簧中的应变能。展开式（2.6.13），得到

\[\pi _ { p } = \frac { 1 } { 2 } k ( u _ { 2 } ^ { 2 } - 2 u _ { 2 } u _ { 1 } + u _ { 1 } ^ { 2 } ) - f _ { 1 x } u _ { 1 } - f _ { 2 x } u _ { 2 }\tag{2.6.14}\]

为了使 $\pi_p$ 对每个节点位移取最小值，需要分别对每个节点位移求偏导，即

\[\begin{array} { l } { \displaystyle { \frac { \partial \pi _ { p } } { \partial u _ { 1 } } } = \frac { 1 } { 2 } k ( - 2 u _ { 2 } + 2 u _ { 1 } ) - f _ { 1 x } = 0 } \\ { \displaystyle { \frac { \partial \pi _ { p } } { \partial u _ { 2 } } } = \frac { 1 } { 2 } k ( 2 u _ { 2 } - 2 u _ { 1 } ) - f _ { 2 x } = 0 } \end{array}\tag{2.6.15}\]

图 2-22 承受节点力的线性弹簧

化简式（2.6.15），得到

\[\begin{array} { r } { k ( - u _ { 2 } + u _ { 1 } ) = f _ { 1 x } } \\ { k ( u _ { 2 } - u _ { 1 } ) = f _ { 2 x } } \end{array}\tag{2.6.16}\]

将式（2.6.16）写成矩阵形式，有

\[\left[ \begin{array} { c c } { k } & { -k } \\ { -k } & { k } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \end{array} \right\} = \left\{ \begin{array} { c } { f _ { 1 x } } \\ { f _ { 2 x } } \end{array} \right\}\tag{2.6.17}\]

由于 $\{f\}=[k]\{d\}$，由式（2.6.17）可得弹簧单元刚度矩阵

\[[ k ] = { \left[ \begin{array} { l l } { k } & { - k } \\ { - k } & { k } \end{array} \right] }\tag{2.6.18}\]

正如预期，式（2.6.18）与第 2.2 节中由式（2.2.9）得到的刚度矩阵完全相同。

上面通过对节点位移最小化总势能，研究了单个弹簧单元的平衡（见例题 2.4），也由此建立了有限元弹簧单元方程。下面进一步说明：整个结构（这里是弹簧单元组合体）的总势能也可以对每个节点自由度取最小值；这种最小化会导出与直接刚度法相同的有限元求解方程。

例题 2.5

求例题 2.1 中弹簧组合体（图 2-23）的总势能，并求其最小值。由此可以看出，单元方程的装配过程也可以由总势能最小化得到。

图 2-23 弹簧组合体

解

由式（2.6.8a），弹簧 1 中储存的应变能为

\[U ^ { ( 1 ) } = k _ { 1 } ( u _ { 3 } - u _ { 1 } ) ^ { 2 } / 2\tag{2.6.19}\]

其中，节点位移差 $u_3-u_1$ 就是弹簧 1 的变形 $x$。式（2.6.19）可写成矩阵形式：

\[U ^ { ( 1 ) } = { \frac { 1 } { 2 } } [ u _ { 3 } \quad u _ { 1 } ] { \left[ \begin{array} { l l } { k _ { 1 } } & { - k _ { 1 } } \\ { - k _ { 1 } } & { k _ { 1 } } \end{array} \right] } { \left\{ \begin{array} { l } { u _ { 3 } } \\ { u _ { 1 } } \end{array} \right\} } = { \frac { 1 } { 2 } } \{ d \} ^ { T } [ K ] \{ d \}\tag{2.6.20}\]

由式（2.6.20）可以看出，应变能 $U$ 是节点位移的二次函数。

弹簧 2 和弹簧 3 的应变能表达式类似，为

\[U ^ { ( 2 ) } = k _ { 2 } ( u _ { 4 } - u _ { 3 } ) ^ { 2 } / 2, \qquad U ^ { ( 3 ) } = k _ { 3 } ( u _ { 2 } - u _ { 4 } ) ^ { 2 } / 2\tag{2.6.21}\]

它们也可以写成与式（2.6.20）类似的矩阵形式。

由于应变能是标量，可以把各弹簧中的能量相加，得到系统总应变能：

\[U = \sum _ { e = 1 } ^ { 3 } U ^ { ( e ) }\tag{2.6.22}\]

按该弹簧组合体的节点编号顺序，外加节点力的势能为

\[\Omega = - ( F _ { 1 x } u _ { 1 } + F _ { 3 x } u _ { 3 } + F _ { 4 x } u _ { 4 } + F _ { 2 x } u _ { 2 } )\tag{2.6.23}\]

式（2.6.23）也可写成矩阵形式：

\[\Omega = - [ u _ { 1 } \quad u _ { 2 } \quad u _ { 3 } \quad u _ { 4 } ] \left\{ \begin{array} { l } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \\ { F _ { 4 x } } \end{array} \right\}\tag{2.6.24}\]

组合体的总势能等于应变能与外力势能之和。把式（2.6.19）、式（2.6.21）和式（2.6.23）相加，得到

\[\Pi _ { p } = U + \Omega = \frac { 1 } { 2 } k _ { 1 } ( u _ { 3 } - u _ { 1 } ) ^ { 2 } + \frac { 1 } { 2 } k _ { 2 } ( u _ { 4 } - u _ { 3 } ) ^ { 2 } + \frac { 1 } { 2 } k _ { 3 } ( u _ { 2 } - u _ { 4 } ) ^ { 2 } - ( F _ { 1x }u_1 + F _ { 2x }u_2 + F _ { 3x }u_3 + F _ { 4x }u_4 )\tag{2.6.25}\]

对每个节点位移最小化 $\pi_p$，得到

\[\begin{array} { l } { \displaystyle \frac { \partial \pi _ { p } } { \partial u _ { 1 } } = - k _ { 1 } u _ { 3 } + k _ { 1 } u _ { 1 } - F _ { 1 x } = 0 } \\ { \displaystyle \frac { \partial \pi _ { p } } { \partial u _ { 2 } } = k _ { 3 } u _ { 2 } - k _ { 3 } u _ { 4 } - F _ { 2 x } = 0 } \\ { \displaystyle \frac { \partial \pi _ { p } } { \partial u _ { 3 } } = k _ { 1 } u _ { 3 } - k _ { 1 } u _ { 1 } - k _ { 2 } u _ { 4 } + k _ { 2 } u _ { 3 } - F _ { 3 x } = 0 } \\ { \displaystyle \frac { \partial \pi _ { p } } { \partial u _ { 4 } } = k _ { 2 } u _ { 4 } - k _ { 2 } u _ { 3 } - k _ { 3 } u _ { 2 } + k _ { 3 } u _ { 4 } - F _ { 4 x } = 0 } \end{array}\tag{2.6.26}\]

把式（2.6.26）写成矩阵形式，可得

\[\left[ \begin{array} { c c c c } { k _ { 1 } } & { 0 } & { - k _ { 1 } } & { 0 } \\ { 0 } & { k _ { 3 } } & { 0 } & { - k _ { 3 } } \\ { - k _ { 1 } } & { 0 } & { k _ { 1 } + k _ { 2 } } & { - k _ { 2 } } \\ { 0 } & { - k _ { 3 } } & { - k _ { 2 } } & { k _ { 2 } + k _ { 3 } } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \\ { u _ { 3 } } \\ { u _ { 4 } } \end{array} \right\} = \left\{ \begin{array} { c } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \\ { F _ { 4 x } } \end{array} \right\}\tag{2.6.27}\]

把 $k_1$、$k_2$ 和 $k_3$ 的数值代入式（2.6.27），得到

\[\left[ \begin{array} { c c c c } { 200 } & { 0 } & { -200 } & { 0 } \\ { 0 } & { 600 } & { 0 } & { -600 } \\ { -200 } & { 0 } & { 600 } & { -400 } \\ { 0 } & { -600 } & { -400 } & { 1000 } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \\ { u _ { 3 } } \\ { u _ { 4 } } \end{array} \right\} = \left\{ \begin{array} { c } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \\ { F _ { 4 x } } \end{array} \right\}\tag{2.6.28}\]

式（2.6.28）与式（2.5.18）完全相同，而式（2.5.18）是第 2.4 节用直接刚度法得到的。因此，用最小势能原理装配得到的方程，与直接刚度装配法得到的方程一致。

总结公式

单元刚度矩阵的定义：

\[\{ f \} = [ k ] \{ d \}\tag{2.1.1}\]

结构整体刚度矩阵或总刚度矩阵的定义：

\[\{ F \} = [ K ] \{ d \}\tag{2.1.2}\]

弹簧单元节点力与节点位移之间的基本矩阵方程：

\[{ \left\{ \begin{array} { l } { f _ { 1 x } } \\ { f _ { 2 x } } \end{array} \right\} } = { \left[ \begin{array} { l l } { k } & { - k } \\ { - k } & { k } \end{array} \right] } { \left\{ \begin{array} { l } { u _ { 1 } } \\ { u _ { 2 } } \end{array} \right\} }\tag{2.2.10}\]

线性弹簧单元刚度矩阵：

\[[ k ] = { \left[ \begin{array} { l l } { k } & { - k } \\ { - k } & { k } \end{array} \right] }\tag{2.2.11}\]

弹簧组合体的整体方程：

\[\{ F \} = [ K ] \{ d \}\tag{2.2.13}\]

总势能：

\[\pi _ { p } = U + \Omega\tag{2.6.1}\]

对于弹簧系统：

\[U = \frac { 1 } { 2 } \{ d \} ^ { T } [ K ] \{ d \}\tag{2.6.20}\]

参考文献

[1] Turner, M. J., Clough, R. W., Martin, H. C., and Topp, L. J., “Stiffness and Deflection Analysis of Complex Structures,” Journal of the Aeronautical Sciences, Vol. 23, No. 9, pp. 805-824, Sept. 1956.

[2] Martin, H. C., Introduction to Matrix Methods of Structural Analysis, McGraw-Hill, New York, 1966.

[3] Hsieh, Y. Y., Elementary Theory of Structures, 2nd ed., Prentice-Hall, Englewood Cliffs, NJ, 1982.

[4] Oden, J. T., and Ripperger, E. A., Mechanics of Elastic Structures, 2nd ed., McGraw-Hill, New York, 1981.

[5] Finlayson, B. A., The Method of Weighted Residuals and Variational Principles, Academic Press, New York, 1972.

[6] Zienkiewicz, O. C., The Finite Element Method, 3rd ed., McGraw-Hill, London, 1977.

[7] Cook, R. D., Malkus, D. S., Plesha, M. E., and Witt, R. J. Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, New York, 2002.

[8] Forray, M. J., Variational Calculus in Science and Engineering, McGraw-Hill, New York, 1968.

习题

2.1 （a）对图 P2-1 所示的弹簧组合体，通过叠加各单个弹簧的刚度矩阵，求整体刚度矩阵 $[K]$。其中 $k_1$、$k_2$ 和 $k_3$ 为图中所示各弹簧的刚度。

（b）若节点 1 和节点 2 固定，且力 $P$ 沿正 $x$ 方向作用于节点 4，求节点 3 和节点 4 位移的表达式。

（c）求节点 1 和节点 2 处的反力。

提示：先仿照第 2.4 节开头部分，写出节点平衡方程，并利用每个单元的力-位移关系求解；然后再用直接刚度法求解同一问题。

图 P2-1

2.2 对图 P2-2 所示的弹簧组合体，求节点 2 的位移和每个弹簧单元中的力。同时求力 $F_3$。已知：由于力 $F_3$ 的作用，节点 3 沿正 $x$ 方向位移 $d=20~\mathrm{mm}$；且 $k_1=k_2=100~\mathrm{N/mm}$。

图 P2-2

2.3 （a）对图 P2-3 所示的弹簧组合体，用直接叠加法求整体刚度矩阵。

（b）若节点 1 和节点 5 固定，并在节点 3 施加力 $P$，求各节点位移。

（c）求固定节点 1 和节点 5 处的反力。

图 P2-3

2.4 求解习题 2.3 的变体：令 $P=0$（即节点 3 不施加外力），并使节点 5 具有图 P2-4 所示的给定已知位移 $d$。

图 P2-4

2.5 对图 P2-5 所示的弹簧组合体，用直接刚度法求整体刚度矩阵。令

\[k^{(1)}=200~\mathrm{N/mm},\quad k^{(2)}=400~\mathrm{N/mm},\quad k^{(3)}=600~\mathrm{N/mm},\quad k^{(4)}=800~\mathrm{N/mm},\quad k^{(5)}=1000~\mathrm{N/mm}.\]

2.6 对图 P2-5 中的弹簧组合体，在节点 2 沿正 $x$ 方向施加 $10000~\mathrm{N}$ 的集中力，求节点 2 和节点 4 的位移。

2.7 不要像图 P2-3 那样假定拉伸单元，而是假定压缩单元。也就是说，对弹簧单元施加压力，并推导其刚度矩阵。

2.8-2.16 对图 P2-8 至图 P2-16 所示的弹簧组合体，求节点位移、每个单元中的力以及反力。所有问题均使用直接刚度法。

图 P2-8

图 P2-9

2.17 对图 P2-17 所示的五弹簧组合体，求节点 2 和节点 3 的位移，以及节点 1 和节点 4 的反力。假定连接节点 2 和节点 3 处弹簧的刚性竖杆始终保持水平，但可以向左或向右滑动或位移。节点 3 上有一个向右的 $1000~\mathrm{N}$ 外力。

图 P2-17

\[\mathrm{Let}\quad k^{(1)}=500~\mathrm{N/mm},\quad k^{(2)}=k^{(3)}=300~\mathrm{N/mm},\quad k^{(4)}=k^{(5)}=400~\mathrm{N/mm}.\]

2.18 使用第 2.6 节建立的最小势能原理求解图 P2-18 所示的弹簧问题。也就是说，针对弹簧自由端位移的不同取值，绘制总势能曲线，以确定最小势能。注意：使总势能取最小值的位移，也对应于稳定平衡位置。

2.19 把例题 2.4 中载荷的方向反过来，重新计算总势能。然后利用该总势能求平衡位移值。

2.20 图 P2-20 中的非线性弹簧满足力-变形关系 $f=k\delta^2$。写出该弹簧的总势能，并用这个势能求平衡位移值。

2.21-2.22 用势能法求解习题 2.10 和习题 2.15（参见例题 2.5）。

2.23 电阻型单元常用于电路分析。考虑图 P2-23 所示的典型电阻单元，它具有节点 1 和节点 2。欧姆定律的一种表述为：两点之间的电势差等于通过导体的电流 $I$ 乘以这两点之间的电阻 $R$。用方程表示为 $V=IR$，其中 $I$ 表示电流，单位为安培（A）；$V$ 表示电势或电压降，单位为伏特（V）；$R$ 表示导体电阻，单位为欧姆（$\Omega$）。使用第 2.2 节的方法，推导把节点处电压降与电流联系起来的“刚度”矩阵：

\[{ \left\{ \begin{array} { l } { V _ { 1 } } \\ { V _ { 2 } } \end{array} \right\} } = R { \left[ \begin{array} { l l } { 1 } & { -1 } \\ { -1 } & { 1 } \end{array} \right] } { \left\{ \begin{array} { l } { I _ { 1 } } \\ { I _ { 2 } } \end{array} \right\} } \qquad \mathrm{or} \qquad \{ V \} = [ K ] \{ I \}\]

图 P2-23

翻译进度说明：第 2 章正文、例题、总结公式、参考文献和习题已补齐，待最终全章校对。