简洁博客

2026-07-03

第 3 章 桁架方程的建立

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

本章目标

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

引言

在第 2 章建立了直接刚度法的基础之后，本章将按照第 1 章给出的一般步骤，推导线弹性杆单元（或桁架单元）的刚度矩阵。推导中会同时引入两类坐标系：一种是根据单元本身选取的局部坐标系，另一种是为了描述整体结构并便于数值计算而选取的整体坐标系或参考坐标系。本章还将讨论如何把向量从局部坐标系变换到整体坐标系，并利用变换矩阵把任意方向杆单元的刚度矩阵表示在整体坐标系中。随后，通过三个平面桁架例题说明如何建立结构的总刚度矩阵和求解方程。图 3-1 给出了典型铁路栈桥平面桁架和伊利诺伊河上的升降桥桁架。

接着，本章把刚度法推广到空间桁架。我们将建立三维空间中的变换矩阵，并分析两个空间桁架问题。之后介绍对称性概念，以及如何利用对称性缩小问题规模、方便求解。本章还会用一个桁架例题说明这一概念，并进一步说明如何处理倾斜支座或斜支座。

随后，我们将用最小势能原理重新推导杆单元方程，并把有限元解与一个承受线性变化分布载荷的杆的精确解进行比较。之后介绍 Galerkin 余量法，并用它推导杆单元方程。最后，还会介绍其他常见余量法，包括配点法、子域法和最小二乘法。本章只把这些方法作为入门介绍，并通过求解承受线性变化载荷的杆问题来展示其应用。

3.1 局部坐标中杆单元刚度矩阵的推导

现在推导图 3-2 所示线弹性、等截面面积（等截面）杆单元的刚度矩阵。这里的推导可直接用于铰接桁架的求解。杆件承受沿杆件局部轴线方向的拉力 $T$，该力作用在节点 1 和节点 2 上。

杆单元假定具有常横截面积 $A$、弹性模量 $E$ 和初始长度 $L$。节点自由度为局部轴向位移，也就是沿杆长方向的纵向位移，分别用单元两端的 $u_1$ 和 $u_2$ 表示，如图 3-2 所示。

由胡克定律式（3.1.1）和应变-位移关系式（3.1.2）或式（1.4.1），可写为

\[\sigma _ { x } = E \varepsilon _ { x }\tag{3.1.1}\]
\[\varepsilon _ { x } = { \frac { d u } { d x } }\tag{3.1.2}\]

由力平衡可得

\[A \sigma _ { x } = T = \mathrm{constant}\tag{3.1.3}\]

这适用于仅在两端施加载荷的杆。（分布载荷将在第 3.10 节讨论。）将式（3.1.2）代入式（3.1.1），再将式（3.1.1）代入式（3.1.3），并对 $x$ 求导，可得到控制线弹性杆行为的微分方程：

（a）

（b） 图 3-1 （a）典型铁路栈桥平面桁架；（b）伊利诺伊河上的升降桥桁架

图 3-2 受拉力 $T$ 作用的杆；正节点位移和正节点力均沿局部 $x$ 方向

\[{ \frac { d } { d x } } \left( A E { \frac { d u } { d x } } \right) = 0\tag{3.1.4}\]

其中，$u$ 是沿单元 $x$ 方向的轴向位移函数。虽然在后续推导中会假定 $A$ 和 $E$ 在整根杆上为常数，但在这个微分方程的一般形式中，仍把 $A$ 和 $E$ 写成似乎可随 $x$ 变化的形式。

推导杆单元刚度矩阵时采用以下假定：

  1. 杆不能承受剪力或弯矩，即
\[f _ { 1 y } = 0, \quad f _ { 2 y } = 0, \quad m _ { 1 } = 0, \quad m _ { 2 } = 0.\]
  1. 忽略横向位移的任何影响。
  2. 胡克定律适用，即轴向应力 $\sigma_x$ 与轴向应变 $\varepsilon_x$ 满足 $\sigma_x=E\varepsilon_x$。
  3. 杆中间没有外加载荷。

下面按照第 1 章给出的步骤推导杆单元刚度矩阵，并说明杆件组合体的完整求解过程。

步骤 1 选择单元类型

用两端节点编号来表示杆单元，通常还要标注单元编号（见图 3-2）。

与推导弹簧单元刚度矩阵类似，在推导一维杆单元刚度矩阵时，暂时可以跳过步骤 2。为简化推导，可直接进入步骤 3。

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

应变-位移关系为

\[\varepsilon _ { x } = { \frac { u _ { 2 } - u _ { 1 } } { L } }\tag{3.1.5}\]

由胡克定律给出的应力-应变关系为

\[\sigma _ { x } = E \varepsilon _ { x }\tag{3.1.6}\]

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

单元刚度矩阵推导如下。由基础力学可知

\[T = A \sigma _ { x }\tag{3.1.7}\]

将式（3.1.5）和式（3.1.6）代入式（3.1.7），得到

\[T = A E \left( { \frac { u _ { 2 } - u _ { 1 } } { L } } \right)\tag{3.1.8}\]

此外，根据图 3-2 中的节点力符号约定，有

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

代入式（3.1.8）后，式（3.1.9）变为

\[f _ { 1 x } = \frac { - A E } { L } ( u _ { 2 } - u _ { 1 } )\tag{3.1.10}\]

类似地，

\[f _ { 2 x } = T\tag{3.1.11}\]

由式（3.1.8）可得

\[f _ { 2 x } = \frac { A E } { L } ( u _ { 2 } - u _ { 1 } )\tag{3.1.12}\]

将式（3.1.10）和式（3.1.12）合写为矩阵形式，有

\[\left\{ \begin{array} { c } { f _ { 1 x } } \\ { f _ { 2 x } } \end{array} \right\} = \frac { A E } { L } \left[ \begin{array} { c c } { 1 } & { -1 } \\ { -1 } & { 1 } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \end{array} \right\}\tag{3.1.13}\]

由于 $\{f\}=[k]\{d\}$，由式（3.1.13）可得

\[[ k ] = { \frac { A E } { L } } \left[ \begin{array} { c c } { 1 } & { -1 } \\ { -1 } & { 1 } \end{array} \right]\tag{3.1.14}\]

式（3.1.14）表示局部坐标中的杆单元刚度矩阵。在式（3.1.14）中，杆单元的 $AE/L$ 类似于弹簧单元的弹簧常数 $k$。

步骤 5 装配单元方程以得到整体方程或总方程

使用第 2 章介绍的直接刚度法，装配整体刚度矩阵、整体力矩阵和整体方程（桁架例题见第 3.6 节）。对于由多个单元组成的结构，这一步可再次写为

\[[ K ] = \sum _ { e = 1 } ^ { N } [ k ^ { ( e ) } ], \qquad \{ F \} = \sum _ { e = 1 } ^ { N } \{ f ^ { ( e ) } \}\tag{3.1.15}\]

其中，在应用式（3.1.15）所示直接刚度法之前，所有局部单元刚度矩阵 $[k^{(e)}]$ 都必须转换成整体单元刚度矩阵，除非局部轴与整体轴重合。坐标变换和刚度矩阵变换的概念将在第 3.3 节和第 3.4 节说明。

步骤 6 求解节点位移

通过施加边界条件，并联立求解方程组 $\{F\}=[K]\{d\}$，即可确定节点位移。

步骤 7 求解单元力

最后，将位移回代入类似式（3.1.5）和式（3.1.6）的方程，求出每个单元中的应变和应力。

下面用一个一维杆问题说明求解过程。

例题 3.1

对图 3-3 所示的三杆组合体，求：（a）整体刚度矩阵；（b）节点 2 和节点 3 的位移；（c）节点 1 和节点 4 处的反力。节点 2 沿 $x$ 方向承受 $15000~\mathrm{N}$ 的外力。每个单元长度均为 $0.6~\mathrm{m}$。对单元 1 和单元 2，取 $E=2.0\times10^{11}~\mathrm{Pa}$、$A=6\times10^{-4}~\mathrm{m^2}$；对单元 3，取 $E=1\times10^{11}~\mathrm{Pa}$、$A=12\times10^{-4}~\mathrm{m^2}$。节点 1 和节点 4 固定。

图 3-3 三杆组合体

解

（a）由式（3.1.14），各单元刚度矩阵为

\[[ k ^ { ( 1 ) } ] = [ k ^ { ( 2 ) } ] = \frac { (6\times10^{-4})(2\times10^{11}) } { 0.6 } \left[ \begin{array} { c c } { 1 } & { -1 } \\ { -1 } & { 1 } \end{array} \right] = 2\times10^8 \left[ \begin{array} { c c } { 1 } & { -1 } \\ { -1 } & { 1 } \end{array} \right]~\frac{\mathrm{N}}{\mathrm{m}}\]
\[[ k ^ { ( 3 ) } ] = \frac { (12\times10^{-4})(1\times10^{11}) } { 0.6 } \left[ \begin{array} { c c } { 1 } & { -1 } \\ { -1 } & { 1 } \end{array} \right] = 2\times10^8 \left[ \begin{array} { c c } { 1 } & { -1 } \\ { -1 } & { 1 } \end{array} \right]~\frac{\mathrm{N}}{\mathrm{m}}\tag{3.1.16}\]

式（3.1.16）中，每个单元矩阵对应的自由度由该单元的节点编号确定。用直接刚度法装配这些单元刚度矩阵，可得整体刚度矩阵

\[[ K ] = 2\times10^8 \left[ \begin{array} { r r r r } { 1 } & { -1 } & { 0 } & { 0 } \\ { -1 } & { 2 } & { -1 } & { 0 } \\ { 0 } & { -1 } & { 2 } & { -1 } \\ { 0 } & { 0 } & { -1 } & { 1 } \end{array} \right]~\frac{\mathrm{N}}{\mathrm{m}}\tag{3.1.17}\]

（b）式（3.1.17）把整体节点力与整体节点位移联系为

\[\left\{ \begin{array} { c } { F _ { 1 x } } \\ { F _ { 2 x } } \\ { F _ { 3 x } } \\ { F _ { 4 x } } \end{array} \right\} = 2\times10^8 \left[ \begin{array} { r r r r } { 1 } & { -1 } & { 0 } & { 0 } \\ { -1 } & { 2 } & { -1 } & { 0 } \\ { 0 } & { -1 } & { 2 } & { -1 } \\ { 0 } & { 0 } & { -1 } & { 1 } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \\ { u _ { 3 } } \\ { u _ { 4 } } \end{array} \right\}\tag{3.1.18}\]

边界条件为

\[u _ { 1 } = 0, \qquad u _ { 4 } = 0\tag{3.1.19}\]

把边界条件和已知外力代入式（3.1.18），并对式（3.1.18）的第 1 个和第 4 个方程进行分块处理，可由第 2 个和第 3 个方程得到

\[\left\{ \begin{array} { c } { 15000 } \\ { 0 } \end{array} \right\} = 2\times10^8 \left[ \begin{array} { r r } { 2 } & { -1 } \\ { -1 } & { 2 } \end{array} \right] \left\{ \begin{array} { c } { u _ { 2 } } \\ { u _ { 3 } } \end{array} \right\}\tag{3.1.20}\]

联立求解式（3.1.20），得节点位移

\[u _ { 2 } = 5\times10^{-5}~\mathrm{m} = 0.05~\mathrm{mm}, \qquad u _ { 3 } = 2.5\times10^{-5}~\mathrm{m} = 0.025~\mathrm{mm}\tag{3.1.21}\]

（c）把式（3.1.19）和式（3.1.21）回代入式（3.1.18），可得整体节点力，其中包括节点 1 和节点 4 处的反力：

\[F _ { 1 x } = 2\times10^8 (u_1-u_2) = 2\times10^8(0-5\times10^{-5}) = -10000~\mathrm{N}\]
\[F _ { 2 x } = 2\times10^8(-u_1+2u_2-u_3) = 2\times10^8[0+2(5\times10^{-5})-2.5\times10^{-5}] = 15000~\mathrm{N}\]
\[F _ { 3 x } = 2\times10^8(-u_2+2u_3-u_4) = 2\times10^8[-5\times10^{-5}+2(2.5\times10^{-5})-0] = 0\]
\[F _ { 4 x } = 2\times10^8(-u_3+u_4) = 2\times10^8(-2.5\times10^{-5}+0) = -5000~\mathrm{N}\tag{3.1.22}\]

式（3.1.22）表明，反力 $F_{1x}$ 和 $F_{4x}$ 的合力，与节点 2 上的 $15000~\mathrm{N}$ 外加节点力大小相等、方向相反，因此杆件组合体满足整体平衡。此外，式（3.1.22）也显示 $F_{2x}=15000~\mathrm{N}$、$F_{3x}=0$ 只是节点 2 和节点 3 上的外加节点力，这进一步验证了求解结果。

3.2 在一维杆单元刚度矩阵推导的步骤 2 中选择位移函数

为一维杆单元选择位移函数时，可参考以下准则。（关于位移函数和其他近似函数，如温度函数的选择，后续章节还会继续讨论：第 4 章讨论梁单元，第 6 章讨论常应变三角形单元，第 8 章讨论线性应变三角形单元，第 9 章讨论轴对称单元，第 10 章讨论三节点杆单元和平面四边形单元，第 11 章讨论三维应力单元，第 12 章讨论板弯曲单元，第 13 章讨论热传导问题。文献 [1-3] 中也给出了更多说明。）

位移函数选择准则

  1. 必须预先选择一个数学函数，用来表示杆单元受载后的变形形状。由于解析闭式解或精确解有时很难得到，甚至不可能得到，因此通常用合适的数学函数假定单元内部的解形状或位移分布。最常用的函数是多项式。

由于杆单元只抵抗轴向载荷，其局部自由度是沿 $x$ 方向的位移 $u_1$ 和 $u_2$，因此选择位移函数 $u$ 表示单元内任意位置的轴向位移。这里假定杆沿 $x$ 轴的位移呈线性变化（图 3-4b），因为给定两个端点后，线性函数的路径唯一。因此有

\[u = a _ { 1 } + a _ { 2 } x\tag{3.2.1}\]

一般而言，系数 $a$ 的总数等于该单元相关自由度的总数。这里总自由度数为 2，即单元两个节点上各有一个轴向位移。式（3.2.1）写成矩阵形式为

\[u = [ 1 \quad x ] { \left\{ \begin{array} { l } { a _ { 1 } } \\ { a _ { 2 } } \end{array} \right\} }\tag{3.2.2}\]

现在希望把 $u$ 表示为节点位移 $u_1$ 和 $u_2$ 的函数。这样就可以像步骤 3 所要求的那样，直接施加节点位移的物理边界条件，并在步骤 4 中把节点位移与节点力联系起来。具体做法是：在每个节点处计算 $u$，再由式（3.2.1）求解 $a_1$ 和 $a_2$：

（a）

（c）

图 3-4 （a）杆单元；（b）位移函数 $u$；（c）形函数 $N_1$；（d）形函数 $N_2$ 在单元域内的分布

\[u (0) = u _ { 1 } = a _ { 1 }\tag{3.2.3}\]
\[u (L) = u _ { 2 } = a _ { 2 } L + u _ { 1 }\tag{3.2.4}\]

由式（3.2.4）解得 $a_2$：

\[a _ { 2 } = \frac { u _ { 2 } - u _ { 1 } } { L }\tag{3.2.5}\]

将式（3.2.3）和式（3.2.5）代入式（3.2.1），得到

\[u = \left( { \frac { u _ { 2 } - u _ { 1 } } { L } } \right) x + u _ { 1 }\tag{3.2.6}\]

写成矩阵形式为

\[u = \left[ 1 - \frac { x } { L } \quad \frac { x } { L } \right] \left\{ \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \end{array} \right\}\tag{3.2.7}\]

或者

\[u = [ N _ { 1 } \quad N _ { 2 } ] \left\{ \begin{array} { c } { u _ { 1 } } \\ { u _ { 2 } } \end{array} \right\}\tag{3.2.8}\]

其中

\[N _ { 1 } = 1 - { \frac { x } { L } }, \qquad N _ { 2 } = { \frac { x } { L } }\tag{3.2.9}\]

称为形函数。这是因为，当第 $i$ 个单元自由度取单位值、其余自由度为零时，$N_i$ 描述了假定位移函数在单元域（$x$ 坐标）内的形状。本例中，$N_1$ 和 $N_2$ 都是线性函数，且 $N_1$ 在节点 1 处为 1、在节点 2 处为 0；$N_2$ 在节点 2 处为 1、在节点 1 处为 0。图 3-4（c）和（d）给出了这些形函数在杆单元域内的分布。此外，对杆上任意轴向坐标，都有 $N_1+N_2=1$。

形函数之和为 1 的意义会在准则 4 中进一步说明。另外，$N_i$ 也常称为插值函数，因为我们是在已知节点值之间插值得到函数值。除端点或节点外，插值函数可以不同于真实函数；但在端点或节点处，插值函数和真实函数必须等于指定的节点值。

  1. 近似函数在杆单元内部应当连续。式（3.2.1）中的简单线性函数显然在单元内部连续。因此，该线性函数能在单元内部给出连续的 $u$ 值，并由于 $u$ 的连续、光滑变化，避免出现开裂、重叠和跳跃（图 3-5）。

图 3-5 二杆结构的单元间连续性

  1. 对离散线单元，近似函数应保证每个节点处所有自由度的单元间连续；对二维和三维单元，则应保证公共边界线和公共边界面上的连续性。对杆单元而言，必须保证两个或多个单元的公共节点在变形后仍为这些单元的公共节点，从而防止单元之间出现重叠或空隙。例如，考虑图 3-5 所示的二杆结构。对该结构，每个单元内部的线性函数 $u$（式 3.2.1）会保证单元 1 和单元 2 仍然相连；单元 1 在节点 2 处的位移等于单元 2 在同一节点 2 处的位移，即 $u_2^{(1)}=u_2^{(2)}$。这个规则也曾由式（2.3.3）说明。因此，该线性函数称为杆单元的协调函数或相容函数，因为它同时保证相邻单元之间的连续性和单元内部的连续性。

一般用符号 $C^m$ 描述分片场变量（例如轴向位移）的连续性，其中上标 $m$ 表示单元间连续的导数阶次。如果函数本身在单元间连续，则该场为 $C^0$ 连续。例如，对图 3-5 所示的轴向位移场，位移在公共节点 2 处连续，因此该位移场称为 $C^0$ 连续。杆单元、平面单元（见第 7 章）和实体单元（见第 11 章）都是 $C^0$ 单元，因为它们强制公共边界上的位移连续。

如果函数的场变量及其一阶导数都在公共边界处连续，则称该场变量为 $C^1$ 连续。后面将看到，梁单元（见第 4 章）和板单元（见第 12 章）是 $C^1$ 连续单元，也就是说，它们会同时强制公共边界上的位移连续和转角/斜率连续。

  1. 近似函数应当允许单元产生刚体位移，并能表示单元内部的常应变状态。一维位移函数式（3.2.1）满足这些要求，因为 $a_1$ 项允许刚体运动，即物体无应变的常量运动；$a_2x$ 项允许常应变，因为 $\varepsilon_x=du/dx=a_2$ 为常数。（实际上，如果单元取得足够小，单元内部就可以出现这种常应变状态。）满足这第四条准则的简单多项式式（3.2.1），称为对杆单元是完备的。

完备性还意味着，通常不能省略低阶项而只保留高阶项。对简单线性函数而言，这意味着不能省略 $a_1$ 而保留 $a_2x$。函数完备性是收敛到精确解的必要条件，例如位移和应力的收敛（图 3-6，见文献 [3]）。图 3-6 展示了随着有限元解中单元数量增加，位移解向精确解的单调收敛。所谓单调收敛，是指一系列近似解（有限元解）一致地向精确解靠近，而不改变符号或方向。

图 3-6 随着有限元解中单元数增加，位移解向精确解收敛

插值函数或近似函数必须允许刚体位移，这意味着函数必须能够给出常数值，例如 $a_1$，因为这种情况确实可能发生。因此考虑

\[u = a _ { 1 }\tag{3.2.10}\]

若 $u=a_1$，则需要节点位移 $u_1=u_2$ 才能得到刚体位移。因此

\[a _ { 1 } = u _ { 1 } = u _ { 2 }\tag{3.2.11}\]

将式（3.2.11）代入式（3.2.8），有

\[u = N _ { 1 } u _ { 1 } + N _ { 2 } u _ { 2 } = ( N _ { 1 } + N _ { 2 } ) a _ { 1 }\tag{3.2.12}\]

由式（3.2.10）和式（3.2.12），得到

\[u = a _ { 1 } = ( N _ { 1 } + N _ { 2 } ) a _ { 1 }\tag{3.2.13}\]

因此，由式（3.2.13）可得

\[N _ { 1 } + N _ { 2 } = 1\tag{3.2.14}\]

所以，式（3.2.14）说明：为了在刚体位移发生时使 $u$ 能够给出常数值，位移插值函数必须在单元内部每一点相加为 1。

3.3 二维向量变换

在许多问题中，同时引入局部坐标系 $(x'-y')$ 和整体坐标系（或参考坐标系）$(x-y)$ 会很方便。局部坐标总是为了便于描述单个单元而选取；整体坐标则为了便于描述整个结构而选取。

给定图 3-7 中用向量 $\mathbf{d}$ 表示的单元节点位移，我们希望建立该向量在一个坐标系中的分量与另一个坐标系中的分量之间的关系。为保持一般性，本节假定 $\mathbf{d}$ 既不与局部轴重合，也不与整体轴重合。在这种情况下，我们希望把整体位移分量与局部位移分量联系起来。这样将得到一个变换矩阵，后续会用它建立杆单元的整体刚度矩阵。定义角度 $\theta$ 从整体 $x$ 轴逆时针转到局部 $x'$ 轴时为正。位移向量 $\mathbf{d}$ 可分别用整体坐标和局部坐标表示为

\[\mathbf { d } = u \mathbf { i } + v \mathbf { j } = u ^ { \prime } \mathbf { i } ^ { \prime } + v ^ { \prime } \mathbf { j } ^ { \prime }\tag{3.3.1}\]

图 3-7 二维一般位移向量 $\mathbf{d}$

其中，$\mathbf{i}$ 和 $\mathbf{j}$ 是整体 $x$ 和 $y$ 方向的单位向量；$\mathbf{i}'$ 和 $\mathbf{j}'$ 是局部 $x'$ 和 $y'$ 方向的单位向量。由图 3-7，可用图中字母表示下列向量：

\[\overline { u } = \overline { OA }, \qquad \overline { v } = \overline { AB }, \qquad \overline { u } ^ { \prime } = \overline { OC }, \qquad \overline { v } ^ { \prime } = \overline { CB }\tag{3.3.2}\]

由图 3-7 沿 $x'$ 轴进行向量相加，可得

\[\overline { OC } = \overline { OD } + \overline { DC }\tag{3.3.3}\]

利用图 3-7 中的三角关系，并结合式（3.3.2），得到

\[\overline { OD } = \overline { OA }\cos\theta = \overline { u }\cos\theta, \qquad \overline { DC } = \overline { AE } = \overline { v }\sin\theta\tag{3.3.4}\]

将式（3.3.2）中 $u'$ 的定义和式（3.3.4）代入式（3.3.3），有

\[u ^ { \prime } = u\cos\theta + v\sin\theta\tag{3.3.5}\]

类似地，在图 3-7 的 $y'$ 方向进行向量相加，可得

\[\overline { CB } = -\overline { AD } + \overline { BE }\tag{3.3.6}\]

再次利用图 3-7 中的三角关系，并结合式（3.3.2），有

\[\overline { AD } = \overline { OA }\sin\theta = \overline { u }\sin\theta, \qquad \overline { BE } = \overline { AB }\cos\theta = \overline { v }\cos\theta\tag{3.3.7}\]

将式（3.3.2）中 $v'$ 的定义和式（3.3.7）代入式（3.3.6），得到

\[v ^ { \prime } = -u\sin\theta + v\cos\theta\tag{3.3.8}\]

把式（3.3.5）和式（3.3.8）合写为矩阵形式，可得

\[\left\{ \begin{array} { c } { u ^ { \prime } } \\ { v ^ { \prime } } \end{array} \right\} = \left[ \begin{array} { c c } { C } & { S } \\ { -S } & { C } \end{array} \right] \left\{ \begin{array} { c } { u } \\ { v } \end{array} \right\}\tag{3.3.9}\]

其中，$C=\cos\theta$，$S=\sin\theta$。

式（3.3.9）把整体位移矩阵 $\{d\}$ 与局部位移矩阵 $\{d'\}$ 联系为

\[\{ d ^ { \prime } \} = [ T ] \{ d \}\tag{3.3.10}\]

其中

\[\{ d \} = \left\{ \begin{array} { c } { u } \\ { v } \end{array} \right\}, \qquad \{ d ^ { \prime } \} = \left\{ \begin{array} { c } { u ^ { \prime } } \\ { v ^ { \prime } } \end{array} \right\}, \qquad [ T ] = \left[ \begin{array} { c c } { C } & { S } \\ { -S } & { C } \end{array} \right]\tag{3.3.11}\]

矩阵 $[T]$ 称为变换矩阵或旋转矩阵。关于该矩阵的更多说明见附录 A。第 3.4 节将用它建立任意方向杆单元的整体刚度矩阵，并把整体节点位移和节点力转换为局部量。

例题 3.2

图 3-8 所示杆单元的节点 2 整体节点位移已知为 $u_2=2.5~\mathrm{mm}$、$v_2=5~\mathrm{mm}$。求节点 2 的局部 $x'$ 方向位移。

图 3-8 局部轴 $x'$ 沿单元方向的杆单元

解

由式（3.3.5），得到

\[u _ { 2 } ^ { \prime } = (\cos60^\circ)(2.5) + (\sin60^\circ)(5) = 5.58~\mathrm{mm}\]

3.4 平面内任意方向杆单元的整体刚度矩阵

现在考虑图 3-9 所示的杆单元。该杆相对于整体 $x$ 轴倾斜，局部轴 $x'$ 沿杆方向由节点 1 指向节点 2。正角度 $\theta$ 取为从整体 $x$ 轴逆时针转到局部 $x'$ 轴。

图 3-9 在整体 $x-y$ 平面内任意方向布置的杆单元

先使用式（3.1.13）。这里用撇号表示局部量：局部单元刚度矩阵 $[k']$ 把局部坐标中的节点力 $\{f'\}$ 与局部节点位移 $\{d'\}$ 联系起来：

\[\left\{ \begin{array} { c } { f _ { 1 x } ^ { \prime } } \\ { f _ { 2 x } ^ { \prime } } \end{array} \right\} = \frac { A E } { L } \left[ \begin{array} { r r } { 1 } & { -1 } \\ { -1 } & { 1 } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } ^ { \prime } } \\ { u _ { 2 } ^ { \prime } } \end{array} \right\}\tag{3.4.1}\]

或者

\[\{ f ^ { \prime } \} = [ k ^ { \prime } ] \{ d ^ { \prime } \}\tag{3.4.2}\]

现在，希望对图 3-9 中相对于整体坐标轴任意方向布置的杆单元，建立整体单元节点力 $\{f\}$ 与整体节点位移 $\{d\}$ 之间的关系。这个关系将给出该单元的整体刚度矩阵 $[k]$。也就是说，我们要找到矩阵 $[k]$，使得

\[\left\{ \begin{array} { c } { f _ { 1 x } } \\ { f _ { 1 y } } \\ { f _ { 2 x } } \\ { f _ { 2 y } } \end{array} \right\} = [k] \left\{ \begin{array} { c } { u _ { 1 } } \\ { v _ { 1 } } \\ { u _ { 2 } } \\ { v _ { 2 } } \end{array} \right\}\tag{3.4.3}\]

简写为

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

由式（3.4.3）可见，在整体坐标中，每个平面杆单元具有 4 个力分量和 4 个位移分量。然而，对弹簧或杆单元的局部坐标表示而言，只有 2 个力分量和 2 个位移分量，如式（3.4.1）所示。通过建立局部力分量与整体力分量之间、局部位移分量与整体位移分量之间的关系，就可以得到整体刚度矩阵。由变换关系式（3.3.5）可知

\[\begin{array} { r } { u _ { 1 } ^ { \prime } = u _ { 1 } \cos\theta + v _ { 1 } \sin\theta } \\ { u _ { 2 } ^ { \prime } = u _ { 2 } \cos\theta + v _ { 2 } \sin\theta } \end{array}\tag{3.4.5}\]

写成矩阵形式为

\[\left\{ \begin{array} { c } { u _ { 1 } ^ { \prime } } \\ { u _ { 2 } ^ { \prime } } \end{array} \right\} = \left[ \begin{array} { c c c c } { C } & { S } & { 0 } & { 0 } \\ { 0 } & { 0 } & { C } & { S } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } } \\ { v _ { 1 } } \\ { u _ { 2 } } \\ { v _ { 2 } } \end{array} \right\}\tag{3.4.6}\]

或者

\[\{ d ^ { \prime } \} = [ T ^ { * } ] \{ d \}\tag{3.4.7}\]

其中

\[[ T ^ { * } ] = \left[ \begin{array} { c c c c } { C } & { S } & { 0 } & { 0 } \\ { 0 } & { 0 } & { C } & { S } \end{array} \right]\tag{3.4.8}\]

类似地，由于力与位移一样都是向量，它们按相同方式变换。因此，将式（3.4.6）中的位移替换为力，可得

\[\left\{ \begin{array} { c } { f _ { 1 x } ^ { \prime } } \\ { f _ { 2 x } ^ { \prime } } \end{array} \right\} = \left[ \begin{array} { c c c c } { C } & { S } & { 0 } & { 0 } \\ { 0 } & { 0 } & { C } & { S } \end{array} \right] \left\{ \begin{array} { c } { f _ { 1 x } } \\ { f _ { 1 y } } \\ { f _ { 2 x } } \\ { f _ { 2 y } } \end{array} \right\}\tag{3.4.9}\]

类似式（3.4.7），式（3.4.9）可写为

\[\{ f ^ { \prime } \} = [ T ^ { * } ] \{ f \}\tag{3.4.10}\]

将式（3.4.7）代入式（3.4.2），得到

\[\{ f ^ { \prime } \} = [ k ^ { \prime } ] [ T ^ { * } ] \{ d \}\tag{3.4.11}\]

再把式（3.4.10）代入式（3.4.11），有

\[[ T ^ { * } ] \{ f \} = [ k ^ { \prime } ] [ T ^ { * } ] \{ d \}\tag{3.4.12}\]

不过，如果要写出单元整体节点力与整体节点位移之间的最终表达式，就需要对式（3.4.12）中的 $[T^]$ 求逆。但 $[T^]$ 不是方阵，因此不能直接求逆。为此，即使 $f_{1y}'$ 和 $f_{2y}'$ 为零，也必须把 $\{d'\}$、$\{f'\}$ 和 $[k']$ 扩展到与整体坐标形式一致的阶数。对每个节点位移使用式（3.3.9），可得

\[\left\{ \begin{array} { c } { u _ { 1 } ^ { \prime } } \\ { v _ { 1 } ^ { \prime } } \\ { u _ { 2 } ^ { \prime } } \\ { v _ { 2 } ^ { \prime } } \end{array} \right\} = \left[ \begin{array} { c c c c } { C } & { S } & { 0 } & { 0 } \\ { -S } & { C } & { 0 } & { 0 } \\ { 0 } & { 0 } & { C } & { S } \\ { 0 } & { 0 } & { -S } & { C } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } } \\ { v _ { 1 } } \\ { u _ { 2 } } \\ { v _ { 2 } } \end{array} \right\}\tag{3.4.13}\]

即

\[\{ d ^ { \prime } \} = [ T ] \{ d \}\tag{3.4.14}\]

其中

\[[ T ] = \left[ \begin{array} { c c c c } { C } & { S } & { 0 } & { 0 } \\ { -S } & { C } & { 0 } & { 0 } \\ { 0 } & { 0 } & { C } & { S } \\ { 0 } & { 0 } & { -S } & { C } \end{array} \right]\tag{3.4.15}\]

同理，因为力与位移一样都是向量，可写为

\[\{ f ^ { \prime } \} = [ T ] \{ f \}\tag{3.4.16}\]

同时，$[k']$ 也必须扩展为 $4\times4$ 矩阵。因此，式（3.4.1）的扩展形式为

\[\left\{ \begin{array} { c } { f _ { 1 x } ^ { \prime } } \\ { f _ { 1 y } ^ { \prime } } \\ { f _ { 2 x } ^ { \prime } } \\ { f _ { 2 y } ^ { \prime } } \end{array} \right\} = \frac { A E } { L } \left[ \begin{array} { c c c c } { 1 } & { 0 } & { -1 } & { 0 } \\ { 0 } & { 0 } & { 0 } & { 0 } \\ { -1 } & { 0 } & { 1 } & { 0 } \\ { 0 } & { 0 } & { 0 } & { 0 } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } ^ { \prime } } \\ { v _ { 1 } ^ { \prime } } \\ { u _ { 2 } ^ { \prime } } \\ { v _ { 2 } ^ { \prime } } \end{array} \right\}\tag{3.4.17}\]

在式（3.4.17）中，由于 $f_{1y}'$ 和 $f_{2y}'$ 为零，$[k']$ 中对应这些力分量的行都是零行。现在将式（3.4.14）和式（3.4.16）代入式（3.4.2），得到

\[[ T ] \{ f \} = [ k ^ { \prime } ] [ T ] \{ d \}\tag{3.4.18}\]

式（3.4.18）就是扩展后的式（3.4.12）。等式两端左乘 $[T]^{-1}$，有

\[\{ f \} = [ T ] ^ { -1 } [ k ^ { \prime } ] [ T ] \{ d \}\tag{3.4.19}\]

其中 $[T]^{-1}$ 是 $[T]$ 的逆矩阵。可以证明（见习题 3.28）

\[[ T ] ^ { -1 } = [ T ] ^ { T }\tag{3.4.20}\]

其中 $[T]^T$ 是 $[T]$ 的转置。满足式（3.4.20）这种性质的方阵称为正交矩阵。关于正交矩阵的更多内容见附录 A。矩形坐标系之间的变换矩阵 $[T]$ 是正交矩阵，这一性质将在本书中反复使用。将式（3.4.20）代入式（3.4.19），得到

\[\{ f \} = [ T ] ^ { T } [ k ^ { \prime } ] [ T ] \{ d \}\tag{3.4.21}\]

比较式（3.4.4）和式（3.4.21），可得单元的整体刚度矩阵

\[[ k ] = [ T ] ^ { T } [ k ^ { \prime } ] [ T ]\tag{3.4.22}\]

将式（3.4.15）中的 $[T]$ 和式（3.4.17）中的扩展 $[k']$ 代入式（3.4.22），可得显式形式

\[[ k ] = \frac { A E } { L } \left[ \begin{array} { c c c c } { C ^2 } & { CS } & { -C ^2 } & { -CS } \\ { CS } & { S ^2 } & { -CS } & { -S ^2 } \\ { -C ^2 } & { -CS } & { C ^2 } & { CS } \\ { -CS } & { -S ^2 } & { CS } & { S ^2 } \end{array} \right]\tag{3.4.23}\]

式（3.4.23）就是在 $x-y$ 平面内任意方向布置的杆单元显式刚度矩阵。

由于试探位移函数式（3.2.6）以及图 3-5 所示假定是逐单元分片连续的，因此可用直接刚度法把各单元刚度矩阵相加，得到

\[\sum _ { e = 1 } ^ { N } [ k ^ { ( e ) } ] = [ K ]\tag{3.4.24}\]

其中 $[K]$ 为总刚度矩阵，$N$ 为单元总数。类似地，各单元整体节点力矩阵也可以相加：

\[\sum _ { e = 1 } ^ { N } \{ f ^ { ( e ) } \} = \{ F \}\tag{3.4.25}\]

于是，$[K]$ 把整个结构的整体节点力 $\{F\}$ 与整体节点位移 $\{d\}$ 联系为

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

例题 3.3

对图 3-10 所示杆单元，求其相对于 $x-y$ 坐标系的整体刚度矩阵。设杆横截面积为 $6\times10^{-4}~\mathrm{m^2}$，长度为 $1.2~\mathrm{m}$，弹性模量为 $2\times10^{11}~\mathrm{Pa}$。杆与 $x$ 轴的夹角为 $30^\circ$。

图 3-10 用于刚度矩阵计算的杆单元

解

为了求杆单元的整体刚度矩阵 $[k]$，使用式（3.4.23）。角度 $\theta$ 定义为从整体 $x$ 轴逆时针转到局部 $x'$ 轴时为正。因此

\[\theta = 30^\circ, \qquad C=\cos30^\circ=\frac{\sqrt{3}}{2}, \qquad S=\sin30^\circ=\frac12\]

于是

\[[ k ] = \frac { (6\times10^{-4})(2\times10^{11}) } { 1.2 } \left[ \begin{array} { c c c c } { \frac34 } & { \frac{\sqrt3}{4} } & { -\frac34 } & { -\frac{\sqrt3}{4} } \\ { \frac{\sqrt3}{4} } & { \frac14 } & { -\frac{\sqrt3}{4} } & { -\frac14 } \\ { -\frac34 } & { -\frac{\sqrt3}{4} } & { \frac34 } & { \frac{\sqrt3}{4} } \\ { -\frac{\sqrt3}{4} } & { -\frac14 } & { \frac{\sqrt3}{4} } & { \frac14 } \end{array} \right]~\frac{\mathrm{N}}{\mathrm{m}}\tag{3.4.27}\]

化简式（3.4.27），得

\[[ k ] = 10^8 \left[ \begin{array} { r r r r } { 0.75 } & { 0.433 } & { -0.75 } & { -0.433 } \\ { 0.433 } & { 0.25 } & { -0.433 } & { -0.25 } \\ { -0.75 } & { -0.433 } & { 0.75 } & { 0.433 } \\ { -0.433 } & { -0.25 } & { 0.433 } & { 0.25 } \end{array} \right]~\frac{\mathrm{N}}{\mathrm{m}}\tag{3.4.28}\]

3.5 $x-y$ 平面内杆单元的应力计算

现在讨论如何确定杆单元中的应力。对杆单元而言，局部力与局部位移之间的关系由式（3.4.1）或式（3.4.17）给出。为方便起见，将该式重复如下：

\[\left\{ \begin{array} { c } { f _ { 1 x } ^ { \prime } } \\ { f _ { 2 x } ^ { \prime } } \end{array} \right\} = \frac { A E } { L } \left[ \begin{array} { r r } { 1 } & { -1 } \\ { -1 } & { 1 } \end{array} \right] \left\{ \begin{array} { c } { u _ { 1 } ^ { \prime } } \\ { u _ { 2 } ^ { \prime } } \end{array} \right\}\tag{3.5.1}\]

轴向拉应力的通常定义是轴向力除以横截面积。因此轴向应力为

\[\sigma = { \frac { f _ { 2 x } ^ { \prime } } { A } }\tag{3.5.2}\]

这里使用 $f_{2x}'$，是因为它是图 3-11 所示拉伸杆件的轴向拉力。由式（3.5.1），有

\[f _ { 2 x } ^ { \prime } = \frac { A E } { L } \left[ -1 \quad 1 \right] \left\{ \begin{array} { c } { u _ { 1 } ^ { \prime } } \\ { u _ { 2 } ^ { \prime } } \end{array} \right\}\tag{3.5.3}\]

因此，将式（3.5.2）和式（3.5.3）合并，可得

图 3-11 具有正节点力的基本杆单元

\[\{ \sigma \} = \frac { E } { L } [ -1 \quad 1 ] \{ d ^ { \prime } \}\tag{3.5.4}\]

再利用式（3.4.7），可得

\[\{ \sigma \} = \frac { E } { L } [ -1 \quad 1 ] [ T ^ { * } ] \{ d \}\tag{3.5.5}\]

式（3.5.5）可进一步简写为

\[\{ \sigma \} = [ C ^ { \prime } ] \{ d \}\tag{3.5.6}\]

其中，若使用式（3.4.8）中的 $[T^*]$，有

\[[ C ^ { \prime } ] = { \frac { E } { L } } [ -1 \quad 1 ] \left[ \begin{array} { c c c c } { C } & { S } & { 0 } & { 0 } \\ { 0 } & { 0 } & { C } & { S } \end{array} \right]\tag{3.5.7}\]

完成式（3.5.7）中的矩阵乘法后，得到

\[[ C ^ { \prime } ] = \frac { E } { L } [ -C \quad -S \quad C \quad S ]\tag{3.5.8}\]

例题 3.4

对图 3-12 所示杆单元，求轴向应力。设 $A=4\times10^{-4}~\mathrm{m^2}$，$E=210~\mathrm{GPa}$，$L=2~\mathrm{m}$，且 $x$ 轴与 $x'$ 轴之间的夹角为 $60^\circ$。假定整体位移已经求得为 $u_1=0.25~\mathrm{mm}$、$v_1=0$、$u_2=0.50~\mathrm{mm}$、$v_2=0.75~\mathrm{mm}$。

解

可用式（3.5.6）计算轴向应力。首先由式（3.5.8）计算 $[C']$：

\[[ C ^ { \prime } ] = { \frac { 210\times10^6~\mathrm{kN/m^2} } { 2~\mathrm{m} } } \left[ -\frac12 \quad -\frac{\sqrt3}{2} \quad \frac12 \quad \frac{\sqrt3}{2} \right]\tag{3.5.9}\]

图 3-12 用于应力计算的杆单元

这里在式（3.5.9）中使用了 $C=\cos60^\circ=1/2$ 和 $S=\sin60^\circ=\sqrt3/2$。现在整体位移向量为

\[\{ d \} = \left\{ \begin{array} { c } { u _ { 1 } } \\ { v _ { 1 } } \\ { u _ { 2 } } \\ { v _ { 2 } } \end{array} \right\} = \left\{ \begin{array} { c } { 0.25\times10^{-3}~\mathrm{m} } \\ { 0 } \\ { 0.50\times10^{-3}~\mathrm{m} } \\ { 0.75\times10^{-3}~\mathrm{m} } \end{array} \right\}\tag{3.5.10}\]

将式（3.5.9）和式（3.5.10）代入式（3.5.6），得到杆的轴向应力

\[\sigma_x = \frac { 210\times10^6 } { 2 } \left[ -\frac12 \quad -\frac{\sqrt3}{2} \quad \frac12 \quad \frac{\sqrt3}{2} \right] \left\{ \begin{array} { c } { 0.25 } \\ { 0 } \\ { 0.50 } \\ { 0.75 } \end{array} \right\}\times10^{-3} = 81.32\times10^3~\mathrm{kN/m^2}=81.32~\mathrm{MPa}\]

3.6 平面桁架的求解

现在用第 3.4 节和第 3.5 节建立的方程，并结合直接刚度法装配总刚度矩阵和方程，来求解下面的平面桁架例题。平面桁架是由同一平面内的杆单元组成、并由无摩擦铰连接的结构。平面桁架的载荷也必须只作用在同一平面内，并且所有载荷都必须施加在节点或铰接点上。

例题 3.5

对图 3-13 所示、由三个单元组成的平面桁架，节点 1 受到向下的 $50~\mathrm{kN}$ 外力。求节点 1 的 $x$、$y$ 方向位移，以及每个单元中的应力。所有单元取 $E=200~\mathrm{GPa}$，$A=6\times10^{-4}~\mathrm{m^2}$。各单元长度见图。

图 3-13 平面桁架

解

首先用式（3.4.23）确定每个单元的整体刚度矩阵。这需要确定每个单元的局部 $x'$ 轴与整体 $x$ 轴之间的角度 $\theta$。本例中，各单元的局部 $x'$ 轴均取为从节点 1 指向另一节点，如图 3-13 所示。每个单元的节点编号可以任意选择；但一旦方向选定，角度 $\theta$ 就按从正 $x$ 轴逆时针转到 $x'$ 轴为正来确定。于是，单元 1 的局部 $x_1'$ 轴从节点 1 指向节点 2，故 $\theta^{(1)}=90^\circ$；单元 2 的局部 $x_2'$ 轴从节点 1 指向节点 3，故 $\theta^{(2)}=45^\circ$；单元 3 的局部 $x_3'$ 轴从节点 1 指向节点 4，故 $\theta^{(3)}=0^\circ$。为便于计算，可建立表 3-1。

表 3-1 图 3-13 桁架的方向余弦数据

单元$\theta$$C$$S$$C^2$$S^2$$CS$
1$90^\circ$01010
2$45^\circ$$\sqrt2/2$$\sqrt2/2$$1/2$$1/2$$1/2$
3$0^\circ$10100

在施加边界约束之前，该桁架共有 8 个节点位移分量，也就是 8 个自由度。因此总刚度矩阵的阶数为 $8\times8$。可以像第 2.4 节前半部分那样，把每个单元的 $[k]$ 矩阵通过添加零行零列扩展为 $8\times8$；也可以像第 2.4 节后半部分那样，按各单元相关的位移分量标注单元刚度矩阵的行和列。本例以及本书后续内容采用后一种方法：把各单元刚度矩阵中的项直接加到 $[K]$ 中对应自由度的位置。

对单元 1，利用式（3.4.23）并结合表 3-1 中的方向余弦，有

\[[ k ^ { (1) } ] = \frac{(2\times10^{11})(6\times10^{-4})}{3} \left[\begin{array}{c c c c} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & -1 \\ 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 1 \end{array}\right]_{(u_1,v_1,u_2,v_2)}\tag{3.6.1}\]

类似地，对单元 2，有

\[[ k ^ { (2) } ] = \frac{(2\times10^{11})(6\times10^{-4})}{3\sqrt2} \left[\begin{array}{r r r r} 0.5 & 0.5 & -0.5 & -0.5 \\ 0.5 & 0.5 & -0.5 & -0.5 \\ -0.5 & -0.5 & 0.5 & 0.5 \\ -0.5 & -0.5 & 0.5 & 0.5 \end{array}\right]_{(u_1,v_1,u_3,v_3)}\tag{3.6.2}\]

对单元 3，有

\[[ k ^ { (3) } ] = \frac{(2\times10^{11})(6\times10^{-4})}{3} \left[\begin{array}{c c c c} 1 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{array}\right]_{(u_1,v_1,u_4,v_4)}\tag{3.6.3}\]

式（3.6.1）至式（3.6.3）中都有公共因子

\[\frac{(2\times10^{11})(6\times10^{-4})}{3}=4\times10^7.\]

在式（3.6.2）中，方括号内各项还要乘以 $1/\sqrt2$。把各单元刚度矩阵的项加到 $[K]$ 中对应位置，可得总刚度矩阵

\[[ K ] = 4\times10^7 \left[\begin{array}{r r r r r r r r} 1.354 & 0.354 & 0 & 0 & -0.354 & -0.354 & -1 & 0 \\ 0.354 & 1.354 & 0 & -1 & -0.354 & -0.354 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 1 & 0 & 0 & 0 & 0 \\ -0.354 & -0.354 & 0 & 0 & 0.354 & 0.354 & 0 & 0 \\ -0.354 & -0.354 & 0 & 0 & 0.354 & 0.354 & 0 & 0 \\ -1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{array}\right]\tag{3.6.4}\]

整体矩阵 $[K]$ 把整体力与整体位移联系起来。考虑节点 1 上的外力以及节点 2、3、4 的边界约束，总结构刚度方程可写为

\[\left\{\begin{array}{c} 0 \\ -50000 \\ F_{2x} \\ F_{2y} \\ F_{3x} \\ F_{3y} \\ F_{4x} \\ F_{4y} \end{array}\right\} = [K] \left\{\begin{array}{c} u_1 \\ v_1 \\ u_2=0 \\ v_2=0 \\ u_3=0 \\ v_3=0 \\ u_4=0 \\ v_4=0 \end{array}\right\}\tag{3.6.5}\]

现在可以使用第 2.5 节前半部分介绍的分块方法，把式（3.6.5）中的前两个方程与第 3 至第 8 个方程分开，以求得未知位移 $u_1$ 和 $v_1$。也可以像第 2.5 节后半部分那样，直接删去总刚度矩阵中与零位移相对应的行和列。这里采用后一种方法：由于第 3 至第 8 行/列对应零位移，删去这些行和列。（注意，对于非齐次边界条件，这种直接删行删列方法必须按第 2.5 节所述进行修正。）于是得到

\[\left\{\begin{array}{c}0\\-50000\end{array}\right\} =4\times10^7 \left[\begin{array}{c c} 1.354 & 0.354\\ 0.354 & 1.354 \end{array}\right] \left\{\begin{array}{c}u_1\\v_1\end{array}\right\}\tag{3.6.6}\]

通过左乘 $2\times2$ 刚度矩阵的逆，或直接联立求解两个方程，可得位移

\[u_1=2.59075\times10^{-4}~\mathrm{m}, \qquad v_1=-9.90925\times10^{-4}~\mathrm{m}\]

$v_1$ 为负，表示节点 1 的 $y$ 向位移与整体坐标中正 $y$ 方向相反，也就是节点 1 向下位移。

利用式（3.5.6）和表 3-1，可求得各单元应力：

\[\sigma^{(1)}=\frac{2\times10^{11}}{3}[0\quad -1\quad 0\quad 1] \left\{\begin{array}{c} u_1=2.59075\times10^{-4}\\v_1=-9.90925\times10^{-4}\\u_2=0\\v_2=0\end{array}\right\} =66.06~\mathrm{MPa}\]
\[\sigma^{(2)}=\frac{2\times10^{11}}{3\sqrt2} \left[-\frac{\sqrt2}{2}\quad -\frac{\sqrt2}{2}\quad \frac{\sqrt2}{2}\quad \frac{\sqrt2}{2}\right] \left\{\begin{array}{c} u_1=2.59075\times10^{-4}\\v_1=-9.90925\times10^{-4}\\u_3=0\\v_3=0\end{array}\right\} =24.40~\mathrm{MPa}\]
\[\sigma^{(3)}=\frac{2\times10^{11}}{3}[-1\quad0\quad1\quad0] \left\{\begin{array}{c} u_1=2.59075\times10^{-4}\\v_1=-9.90925\times10^{-4}\\u_4=0\\v_4=0\end{array}\right\} =-17.27~\mathrm{MPa}\]

最后，可通过检查节点 1 处的整体平衡来验证结果。沿整体 $x$ 方向求和，有

\[\sum F_x=0\]
\[(24.40~\mathrm{MPa})(6\times10^{-4}~\mathrm{m^2})\frac{\sqrt2}{2}-(17.27~\mathrm{MPa})(6\times10^{-4}~\mathrm{m^2})=0\]

沿整体 $y$ 方向求和，有

\[\sum F_y=0\]
\[(66.06~\mathrm{MPa})(6\times10^{-4}~\mathrm{m^2})+(24.40~\mathrm{MPa})(6\times10^{-4}~\mathrm{m^2})\frac{\sqrt2}{2}-50000=0\]

例题 3.6

对于图 3-14 所示的两杆桁架，求节点 1 沿 $y$ 方向的位移，以及各单元中的轴向力。节点 1 受到沿正 $y$ 方向的力 $P=1000~\mathrm{kN}$，同时节点 1 沿负 $x$ 方向发生沉降，沉降量为 $\delta=50~\mathrm{mm}$。各单元取 $E=210~\mathrm{GPa}$、$A=6.00\times10^{-4}~\mathrm{m^2}$。单元长度见图。

图 3-14 两杆桁架

图 3-14 两杆桁架

解

仍从式（3.4.23）出发，分别求各单元刚度矩阵。

单元 1

\[\cos\theta^{(1)}=\frac{3}{5}=0.60\]
\[\sin\theta^{(1)}=\frac{4}{5}=0.80\]

于是

\[[k^{(1)}]=\frac{(6.0\times10^{-4}~\mathrm{m^2})(210\times10^6~\mathrm{kN/m^2})}{5~\mathrm{m}} \left[ \begin{array}{c c c c} 0.36&0.48&-0.36&-0.48\\ 0.48&0.64&-0.48&-0.64\\ -0.36&-0.48&0.36&0.48\\ -0.48&-0.64&0.48&0.64 \end{array} \right]\tag{3.6.7}\]

化简式（3.6.7），得

\[[k^{(1)}]=25200 \left[ \begin{array}{c c c c} 0.36&0.48&-0.36&-0.48\\ 0.48&0.64&-0.48&-0.64\\ -0.36&-0.48&0.36&0.48\\ -0.48&-0.64&0.48&0.64 \end{array} \right] \quad (u_1,v_1,u_2,v_2)\tag{3.6.8}\]

单元 2

\[\cos\theta^{(2)}=0.0, \qquad \sin\theta^{(2)}=1.0\]
\[[k^{(2)}]=\frac{(6.0\times10^{-4})(210\times10^6)}{4} \left[ \begin{array}{c c c c} 0&0&0&0\\ 0&1&0&-1\\ 0&0&0&0\\ 0&-1&0&1 \end{array} \right]\tag{3.6.9}\]

为便于与式（3.6.8）使用相同的矩阵前因子，将式（3.6.9）写成

\[[k^{(2)}]=25200 \left[ \begin{array}{c c c c} 0&0&0&0\\ 0&1.25&0&-1.25\\ 0&0&0&0\\ 0&-1.25&0&1.25 \end{array} \right] \quad (u_1,v_1,u_3,v_3)\tag{3.6.10}\]

将式（3.6.8）和式（3.6.10）中的单元刚度矩阵叠加，可得整体刚度矩阵，并有整体力与整体位移关系

\[\left\{ \begin{array}{c} F_{1x}\\F_{1y}\\F_{2x}\\F_{2y}\\F_{3x}\\F_{3y} \end{array} \right\} =25200 \left[ \begin{array}{c c c c c c} 0.36&0.48&-0.36&-0.48&0&0\\ 0.48&1.89&-0.48&-0.64&0&-1.25\\ -0.36&-0.48&0.36&0.48&0&0\\ -0.48&-0.64&0.48&0.64&0&0\\ 0&0&0&0&0&0\\ 0&-1.25&0&0&0&1.25 \end{array} \right] \left\{ \begin{array}{c} u_1\\v_1\\u_2\\v_2\\u_3\\v_3 \end{array} \right\}\tag{3.6.11}\]

我们仍然可以把与已知位移对应的方程分离出来，只对未知位移对应的方程联立求解。边界条件为

\[u_1=\delta, \qquad u_2=0, \qquad v_2=0, \qquad u_3=0, \qquad v_3=0\tag{3.6.12}\]

因此，将式（3.6.11）中的第 2 个方程从第 1、3、4、5、6 个方程中分离出来，并代入边界条件，可得

\[P=25200(0.48\delta+1.89v_1)\tag{3.6.13}\]

其中 $F_{1y}=P$，且 $u_1=\delta$。将式（3.6.13）写成 $P$ 与 $\delta$ 的形式，可以清楚看出外力和支座沉降对 $v_1$ 的贡献。解得

\[v_1=0.000021P-0.254\delta\tag{3.6.14}\]

代入 $P=1000~\mathrm{kN}$、$\delta=-0.05~\mathrm{m}$，得

\[v_1=0.0337~\mathrm{m}\tag{3.6.15}\]

正值表示节点 1 的位移沿整体坐标的正 $y$ 方向。

局部单元力可由式（3.4.11）求得。具体如下。

单元 1

\[\left\{ \begin{array}{c} f'_{1x}\\f'_{2x} \end{array} \right\} =25200 \left[ \begin{array}{c c} 1&-1\\-1&1 \end{array} \right] \left[ \begin{array}{c c c c} 0.60&0.80&0&0\\ 0&0&0.60&0.80 \end{array} \right] \left\{ \begin{array}{c} u_1=-0.05\\v_1=0.0337\\u_2=0\\v_2=0 \end{array} \right\}\tag{3.6.16}\]

完成式（3.6.16）中的矩阵三重乘积，得到

\[f'_{1x}=-76.6~\mathrm{kN}, \qquad f'_{2x}=76.6~\mathrm{kN}\tag{3.6.17}\]

单元 2

\[\left\{ \begin{array}{c} f'_{1x}\\f'_{3x} \end{array} \right\} =31500 \left[ \begin{array}{c c} 1&-1\\-1&1 \end{array} \right] \left[ \begin{array}{c c c c} 0&1&0&0\\ 0&0&0&1 \end{array} \right] \left\{ \begin{array}{c} u_1=-0.05\\v_1=0.0337\\u_3=0\\v_3=0 \end{array} \right\}\tag{3.6.18}\]

完成式（3.6.18）中的矩阵三重乘积，得到

\[f'_{1x}=1061~\mathrm{kN}, \qquad f'_{3x}=-1061~\mathrm{kN}\tag{3.6.19}\]

可通过检查节点 1 处是否满足平衡来验证上述计算。

例题 3.7

为说明如何在同一结构中组合弹簧单元和杆单元，下面求解图 3-15 所示、由弹簧支承的两杆桁架。两根杆件均取 $E=210~\mathrm{GPa}$、$A=5.0\times10^{-4}~\mathrm{m^2}$。杆 1 的长度为 $5~\mathrm{m}$，杆 2 的长度为 $10~\mathrm{m}$。弹簧刚度为 $k=2000~\mathrm{kN/m}$。

图 3-15 带弹簧支承的两杆桁架

图 3-15 带弹簧支承的两杆桁架

解

首先使用式（3.4.23）分别确定各单元的刚度矩阵。

单元 1

\[\theta^{(1)}=135^\circ, \qquad \cos\theta^{(1)}=-\frac{\sqrt2}{2}, \qquad \sin\theta^{(1)}=\frac{\sqrt2}{2}\]
\[\left[k^{(1)}\right] =\frac{(5.0\times10^{-4}~\mathrm{m^2})(210\times10^{6}~\mathrm{kN/m^2})}{5~\mathrm{m}} \left[ \begin{array}{r r r r} 0.5&-0.5&-0.5&0.5\\ -0.5&0.5&0.5&-0.5\\ -0.5&0.5&0.5&-0.5\\ 0.5&-0.5&-0.5&0.5 \end{array} \right] \tag{3.6.20}\]

化简式（3.6.20），得

\[\left[k^{(1)}\right] =105\times10^{2} \left[ \begin{array}{c|r r r r} &u_1&v_1&u_2&v_2\\ \hline u_1&1&-1&-1&1\\ v_1&-1&1&1&-1\\ u_2&-1&1&1&-1\\ v_2&1&-1&-1&1 \end{array} \right] \tag{3.6.21}\]

单元 2

\[\theta^{(2)}=180^\circ, \qquad \cos\theta^{(2)}=-1.0, \qquad \sin\theta^{(2)}=0\]
\[\left[k^{(2)}\right] =\frac{(5.0\times10^{-4}~\mathrm{m^2})(210\times10^{6}~\mathrm{kN/m^2})}{10~\mathrm{m}} \left[ \begin{array}{r r r r} 1&0&-1&0\\ 0&0&0&0\\ -1&0&1&0\\ 0&0&0&0 \end{array} \right] \tag{3.6.22}\]

化简式（3.6.22），得

\[\left[k^{(2)}\right] =105\times10^{2} \left[ \begin{array}{c|r r r r} &u_1&v_1&u_3&v_3\\ \hline u_1&1&0&-1&0\\ v_1&0&0&0&0\\ u_3&-1&0&1&0\\ v_3&0&0&0&0 \end{array} \right] \tag{3.6.23}\]

单元 3

\[\theta^{(3)}=270^\circ, \qquad \cos\theta^{(3)}=0, \qquad \sin\theta^{(3)}=-1.0\]

使用式（3.4.23），但将 $AE/L$ 替换为弹簧常数 $k$，可得弹簧的刚度矩阵为

\[\left[k^{(3)}\right] =20\times10^{2} \left[ \begin{array}{c|r r r r} &u_1&v_1&u_4&v_4\\ \hline u_1&0&0&0&0\\ v_1&0&1&0&-1\\ u_4&0&0&0&0\\ v_4&0&-1&0&1 \end{array} \right] \tag{3.6.24}\]

施加边界条件，有

\[u_2=v_2=u_3=v_3=u_4=v_4=0\tag{3.6.25}$$ 在组装后的整体方程中代入式（3.6.25）的边界条件，可得到缩减方程 \]

\left\{ \begin{array}{c} F_{1x}\\ F_{1y}=-25~\mathrm{kN} \end{array} \right\} =10^{2} \left[ \begin{array}{r r} 210&-105\\ -105&125 \end{array} \right] \left\{ \begin{array}{c} u_1\\v_1 \end{array} \right\} \tag{3.6.26}

\[ 求解式（3.6.26），得到整体位移 \]

u_1=-1.724\times10^{-3}~\mathrm{m}, \qquad v_1=-3.448\times10^{-3}~\mathrm{m}\tag{3.6.27}$$

可用式（3.5.6）求杆单元中的应力。对单元 1，有

\[\sigma^{(1)} =\frac{210\times10^{3}~\mathrm{MN/m^2}}{5~\mathrm{m}} \left[0.707\quad -0.707\quad -0.707\quad 0.707\right] \left\{ \begin{array}{c} -1.724\times10^{-3}\\ -3.448\times10^{-3}\\ 0\\ 0 \end{array} \right\}\]

化简得

\[\sigma^{(1)}=51.2~\mathrm{MPa}\quad(\text{拉})\]

同理，可得单元 2 中的应力为

\[\sigma^{(2)} =\frac{210\times10^{3}~\mathrm{MN/m^2}}{10~\mathrm{m}} \left[1.0\quad 0\quad -1.0\quad 0\right] \left\{ \begin{array}{c} -1.724\times10^{-3}\\ -3.448\times10^{-3}\\ 0\\ 0 \end{array} \right\}\]

化简得

\[\sigma^{(2)}=-36.2~\mathrm{MPa}\quad(\text{压})\]

3.7 三维空间中杆单元的变换矩阵和刚度矩阵

现在推导一个在三维空间中任意取向的杆单元所需的变换矩阵，并由此得到其一般刚度矩阵，如图 3-16 所示。设节点 1 的坐标为 $x_1,y_1,z_1$，节点 2 的坐标为 $x_2,y_2,z_2$。再设 $\theta_x,\theta_y,\theta_z$ 分别为整体 $x,y,z$ 轴到局部 $x'$ 轴之间的夹角。这里，$x'$ 轴沿单元从节点 1 指向节点 2。现在需要确定 $[T^*]$，使得

\[\{d'\}=[T^*]\{d\}\]

推导 $[T^*]$ 时，先考虑三维中同一位移向量 $\mathbf d'=\mathbf d$ 的表达式：

\[u'\mathbf i'+v'\mathbf j'+w'\mathbf k'=u\mathbf i+v\mathbf j+w\mathbf k\tag{3.7.1}\]

其中 $\mathbf i',\mathbf j',\mathbf k'$ 分别为局部 $x',y',z'$ 轴对应的单位向量，$\mathbf i,\mathbf j,\mathbf k$ 分别为整体 $x,y,z$ 轴对应的单位向量。同时，$w$ 和 $w'$ 分别表示整体 $z$ 方向和局部 $z'$ 方向的位移。对式（3.7.1）两边与 $\mathbf i'$ 作点积，有

\[u'+0+0=u(\mathbf i'\cdot\mathbf i)+v(\mathbf i'\cdot\mathbf j)+w(\mathbf i'\cdot\mathbf k)\tag{3.7.2}\]

根据点积的定义，有

图 3-16 三维空间中的杆件及其局部节点位移

图 3-16 三维空间中的杆件及其局部节点位移

\[\begin{array}{l} \mathbf i'\cdot\mathbf i=\dfrac{x_2-x_1}{L}=C_x\\[6pt] \mathbf i'\cdot\mathbf j=\dfrac{y_2-y_1}{L}=C_y\\[6pt] \mathbf i'\cdot\mathbf k=\dfrac{z_2-z_1}{L}=C_z \end{array} \tag{3.7.3}\]

其中

\[L=\left[(x_2-x_1)^2+(y_2-y_1)^2+(z_2-z_1)^2\right]^{1/2}\]

且

\[C_x=\cos\theta_x, \qquad C_y=\cos\theta_y, \qquad C_z=\cos\theta_z\tag{3.7.4}\]

这里，$C_x,C_y,C_z$ 分别是局部轴单位向量 $\mathbf i'$ 在 $\mathbf i,\mathbf j,\mathbf k$ 方向上的投影。因此，将式（3.7.3）代入式（3.7.2），可得

\[u'=C_x u+C_y v+C_z w\tag{3.7.5}\]

对于沿 $x'$ 轴方向的空间向量，式（3.7.5）给出了该向量在整体 $x,y,z$ 方向上的分量关系。利用式（3.7.5），可将节点 1 和节点 2 的局部轴向位移显式写成

\[\left\{ \begin{array}{c} u_1'\\u_2' \end{array} \right\} = \left[ \begin{array}{c c c c c c} C_x&C_y&C_z&0&0&0\\ 0&0&0&C_x&C_y&C_z \end{array} \right] \left\{ \begin{array}{c} u_1\\v_1\\w_1\\u_2\\v_2\\w_2 \end{array} \right\} \tag{3.7.6}\]

现定义

\[\{d'\}=\left\{\begin{array}{c}u_1'\\u_2'\end{array}\right\},\qquad \{d\}=\left\{\begin{array}{c}u_1\\v_1\\w_1\\u_2\\v_2\\w_2\end{array}\right\}\]

并定义

\[[T^*]= \left[ \begin{array}{c c c c c c} C_x&C_y&C_z&0&0&0\\ 0&0&0&C_x&C_y&C_z \end{array} \right] \tag{3.7.7}\]

利用式（3.7.7），可将式（3.7.6）写成矩阵形式：

\[\{d'\}=[T^*]\{d\}\tag{a}\]

其中，$[T^*]$ 是变换矩阵，它使局部位移矩阵 $\{d'\}$ 能够用整体坐标系中的位移矩阵 $\{d\}$ 的分量来表示。

根据式（a），可方便地用 $[T^*]$ 将整体力矩阵表示为局部力矩阵的形式：

\[\{f\}=[T^*]^T\{f'\}\tag{b}\]

在局部坐标中，局部力与局部位移之间满足

\[\{f'\}=[k']\{d'\}\tag{c}\]

将式（a）中的 $\{d'\}$ 代入式（c），并在等式两边左乘 $[T^*]^T$，得

\[[T^*]^T\{f'\}=[T^*]^T[k'][T^*]\{d\}\tag{d}\]

再利用式（b），可得

\[\{f\}=[T^*]^T[k'][T^*]\{d\}\tag{e}\]

整体力与整体位移之间的关系为

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

比较式（e）和式（f）右端，可知空间中任意取向杆单元的整体刚度矩阵为

\[[k]=[T^*]^T[k'][T^*]\tag{g}\]

将式（3.7.7）中的 $[T^*]$，以及式（3.4.1）和式（3.4.2）给出的 $[k']$ 代入式（g），得

\[[k]= \left[ \begin{array}{c c} C_x&0\\ C_y&0\\ C_z&0\\ 0&C_x\\ 0&C_y\\ 0&C_z \end{array} \right] \frac{AE}{L} \left[ \begin{array}{r r} 1&-1\\ -1&1 \end{array} \right] \left[ \begin{array}{c c c c c c} C_x&C_y&C_z&0&0&0\\ 0&0&0&C_x&C_y&C_z \end{array} \right] \tag{3.7.8}\]

化简式（3.7.8），可得 $[k]$ 的显式形式：

\[[k]=\frac{AE}{L} \left[ \begin{array}{c c c c c c} C_x^2&C_xC_y&C_xC_z&-C_x^2&-C_xC_y&-C_xC_z\\ C_yC_x&C_y^2&C_yC_z&-C_yC_x&-C_y^2&-C_yC_z\\ C_zC_x&C_zC_y&C_z^2&-C_zC_x&-C_zC_y&-C_z^2\\ -C_x^2&-C_xC_y&-C_xC_z&C_x^2&C_xC_y&C_xC_z\\ -C_yC_x&-C_y^2&-C_yC_z&C_yC_x&C_y^2&C_yC_z\\ -C_zC_x&-C_zC_y&-C_z^2&C_zC_x&C_zC_y&C_z^2 \end{array} \right] \tag{3.7.9}\]

式（3.7.9）是三维空间中任意取向杆单元刚度矩阵的基本形式。下面将分析一个简单的空间桁架，以说明本节所建立的概念。可以看到，直接刚度法为求解空间桁架问题提供了清晰而直接的步骤。

翻译进度说明：3.7 节已补齐；后续将继续翻译例题 3.8。