任意方向桁架单元的刚度矩阵

简单的 $\frac{E A}{L} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}]$ 刚度矩阵（假设横截面积和杨氏模量恒定）仅适用于与坐标轴对齐的桁架单元。在实际结构中，桁架单元可以任意角度定向。为处理此情况，我们必须推导一个能在共同全局坐标系（x，y）中关联力与位移的刚度矩阵。

过程如下：1. 从与单元对齐的局部坐标系（x’，y’）中的简单刚度矩阵开始。2. 建立一个关联局部和全局坐标系的变换矩阵 T。3. 利用此变换矩阵将局部刚度矩阵转换为所需的全局刚度矩阵。

局部刚度矩阵（K^local）

对于沿其局部 x’ 轴对齐的桁架单元，它只能抵抗沿该轴的力并产生沿该轴的位移。垂直于单元（y’ 方向）的任何力或位移都与零刚度相关，因为它不会引起长度变化。

因此，我们可以将该单元在局部二维坐标系中的刚度矩阵写为（四个自由度：u₁’、v₁’、u₂’、v₂’）：

这个 4×4 矩阵就是局部坐标系下的刚度矩阵，即 K^local。

坐标变换

设全局坐标系（x，y）相对于局部坐标系（x’，y’）旋转一个角度 θ。

前面我们学过，任意向量 $𝐕$ 都可以借助标准的二维旋转矩阵从一个坐标系变换到另一个坐标系。局部向量分量 [V_x’， V_y’] 与全局分量 [V_x， V_y] 之间的关系为：对两边取转置，得到令 c = cos(θ)，s = sin(θ)。旋转矩阵为 $𝐑 = [\begin{matrix} c & s \\ - s & c \end{matrix}]$ 我们可以将其应用于各节点的位移向量。设全局位移为 q，局部位移为 u’：

在节点 1：
在节点 2：

我们可以将它们组合为整个单元的单一变换。这就定义了 变换矩阵 T：

力向量也存在类似关系：其中 F 是全局 xy 坐标系中的节点力向量，F’ 是局部 x’y’ 坐标系中的节点力向量。

由于 T 是正交矩阵，其逆矩阵等于其转置 T^-1 = TT。因此，由上式可得：其中 #### 推导全局刚度矩阵（K^global）

现在我们可以通过组合这些关系推导出全局坐标系下的刚度矩阵：

从局部系统的基本关系开始：
代入位移变换 u’ = Tq：
将两边左乘 T^T 以获得全局力 F：
注意左侧是全局力向量，即 F = **T^TF’： $𝐅 = (𝐓^{𝖳} 𝐊^{local} 𝐓) 𝐪$

此方程现在为全局形式 F = K^globalq。因此，我们已经证明了刚度矩阵的主变换方程：

通过完成该矩阵乘法，我们得到任意定向桁架单元的全局刚度矩阵的显式形式：

$𝐊^{global} = \frac{E A}{L} [\begin{matrix} c^{2} & c s & - c^{2} & - c s \\ c s & s^{2} & - c s & - s^{2} \\ - c^{2} & - c s & c^{2} & c s \\ - c s & - s^{2} & c s & s^{2} \end{matrix}]$

现在，该 4×4 矩阵可用于第一节中描述的直接叠加（组装）过程，以构建整个桁架结构的全局刚度矩阵。

以下示例演示了坐标变换和刚度组装的完整过程。

示例：三杆桁架结构分析

考虑以下承受一个水平力和一个垂直力的桁架（见图）。假设所有杆的杨氏模量和横截面积相同：E = 100 GPa，A = 200 mm²。

全局坐标系与节点： 我们以节点 1 为原点建立全局 $(x, y)$ 坐标系。

节点 1： (0, 0)
节点 2： (0, 2)
节点 3： (2, 2)

全局自由度（DOF）： 该结构有 3 个节点，每个节点 2 个自由度（ $u, v$ ），共 6 个全局自由度。我们先对自由 DOF 编号，然后对支座（约束）DOF 编号。这将简化后续计算。

自由 DOF（节点 3 处）：
- $q_{1}$ ：节点 1 的垂直位移
- $q_{2}$ ：节点 3 的水平位移
- $q_{3}$ ：节点 3 的垂直位移
支座 DOF（节点 1 和 2 处）：
- $q_{4}$ ：节点 1 的水平位移
- $q_{5}$ ：节点 2 的水平位移
- $q_{6}$ ：节点 2 的垂直位移

因为节点 1 为辊轴支座，节点 2 为固定铰支座，其相应位移必须为零： $q_{4} = q_{5} = q_{6} = 0$

步骤 2：推导各单元在全局坐标系下的刚度矩阵

现在，我们分别分析三个桁架单元，计算其在全局坐标系下的 4×4 刚度矩阵，并用相应的全局自由度编号显式标记行和列。

单元 ①（节点 1 至 2）：

长度： $L_{1} = 2 m$
角度： 该单元为竖直，故 $θ = 90^{\circ}$ 。
余弦： $c = \cos (90^{\circ}) = 0$ ， $s = \sin (90^{\circ}) = 1$ 。
自由度映射： 该单元的局部自由度（ $u_{1}, v_{1}, u_{2}, v_{2}$ ）对应于全局自由度 ( $q_{4}, q_{1}, q_{5}, q_{6}$ )。
刚度矩阵： $\begin{array}{cc} \begin{array}{cccc} q_{4} & q_{1} & q_{5} & q_{6} \end{array} \\ 𝐊_{e}^{(1)} = \frac{E A}{L} & [\begin{array}{cccc} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & - 1 \\ 0 & 0 & 0 & 0 \\ 0 & - 1 & 0 & 1 \end{array}] \end{array} \begin{matrix} q_{4} \\ q_{1} \\ q_{5} \\ q_{6} \end{matrix}$

单元 ②（节点 2 至 3）：

长度： $L_{2} = 2 m$
角度： 该单元为水平，故 $θ = 0^{\circ}$ 。
余弦： $c = \cos (0^{\circ}) = 1$ ， $s = \sin (0^{\circ}) = 0$ 。
自由度映射： 该单元的局部自由度（ $u_{2}, v_{2}, u_{3}, v_{3}$ ）对应于全局自由度 ( $q_{5}, q_{6}, q_{2}, q_{3}$ )。
刚度矩阵： $\begin{array}{cc} \begin{array}{cccc} q_{5} & q_{6} & q_{2} & q_{3} \end{array} \\ 𝐊_{e}^{(2)} = \frac{E A}{L} & [\begin{array}{cccc} 1 & 0 & - 1 & 0 \\ 0 & 0 & 0 & 0 \\ - 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{array}] \end{array} \begin{matrix} q_{5} \\ q_{6} \\ q_{2} \\ q_{3} \end{matrix}$

单元 ③（节点 1 至 3）：

长度： $L_{3} = \sqrt{(2 - 0)^{2} + (2 - 0)^{2}} = 2 \sqrt{2} m$
角度： 该单元为 45°，故 $θ = 45^{\circ}$ 。
余弦： $c = \cos (45^{\circ}) = 1 / \sqrt{2}$ ， $s = \sin (45^{\circ}) = 1 / \sqrt{2}$ 。
自由度映射： 该单元的局部自由度（ $u_{1}, v_{1}, u_{3}, v_{3}$ ）对应于全局自由度 ( $q_{4}, q_{1}, q_{2}, q_{3}$ )。
刚度矩阵： 对于该单元， $c^{2} = 1 / 2$ ， $s^{2} = 1 / 2$ ， $c s = 1 / 2$ 。 $\begin{array}{cc} \begin{array}{cccc} q_{4} & q_{1} & q_{2} & q_{3} \end{array} \\ 𝐊_{e}^{(3)} = \frac{E A}{2 \sqrt{2} L} & [\begin{array}{cccc} 1 & 1 & - 1 & - 1 \\ 1 & 1 & - 1 & - 1 \\ - 1 & - 1 & 1 & 1 \\ - 1 & - 1 & 1 & 1 \end{array}] \end{array} \begin{matrix} q_{4} \\ q_{1} \\ q_{2} \\ q_{3} \end{matrix}$

步骤 3：组装全局刚度矩阵（ $𝐊$ ）

我们现在通过将每个单元矩阵的分量按照全局自由度标签添加到正确位置，来组装 6×6 全局刚度矩阵。最终的矩阵 $𝐊$ 是三个单元刚度矩阵扩展至 6×6 全局空间后的总和。

例如，我们来求全局项 $K_{11}$ ，它对应于全局自由度 $q_{1}$ 和 $q_{1}$ 之间的相互作用。我们在各单元矩阵中查找 (q₁, q₁) 项：

单元 ①： 行和列处的项为 $\frac{E A}{L}$ 。
单元 ②： 没有相应的行或列，其对 $K_{11}$ 的贡献为 0。
单元 ③： 行 $q_{1}$ 与列 $q_{1}$ 处的项为 $\frac{E A}{2 \sqrt{2} L}$ 。
合计： $K_{11} = \frac{E A}{L} + 0 + \frac{E A}{2 \sqrt{2} L} = E A (\frac{1}{2} + \frac{1}{4 \sqrt{2}})$ 。

再举一个例子，求非对角线项 $K_{16}$ （第 1 行，第 6 列）：

单元 ①： 行 $q_{1}$ 与列 $q_{6}$ 处的项为 $\frac{- E A}{L}$ 。
单元 ②： 没有 $q_{6}$ 和 $q_{1}$ 自由度，贡献为 0。
单元 ③： 没有 $q_{6}$ 和 $q_{1}$ 自由度，贡献为 0。
合计： $K_{16} = - \frac{E A}{L} + 0 + 0 = - \frac{E A}{L}$ 。

根据全局行、列索引（直接叠加）添加数值的过程，需对全局刚度矩阵的全部 36 个元素重复进行。

步骤 4：划分与求解系统

完整系统为 $𝐏 = Kq$ 。该组装后的 6×6 $𝐊$ 矩阵是奇异的，因为我们尚未施加支座条件。我们通过将系统划分为自由 DOF 和支座 DOF 来求解。

${\begin{matrix} 𝐏_{f} \\ 𝐏_{s} \end{matrix}} = [\begin{matrix} 𝐊_{f f} & 𝐊_{f s} \\ 𝐊_{s f} & 𝐊_{s s} \end{matrix}] {\begin{matrix} 𝐪_{f} \\ 𝐪_{s} \end{matrix}}$

其中各分块为： * 自由位移： $𝐪_{f} = {\begin{matrix} q_{1} \\ q_{2} \\ q_{3} \end{matrix}}$ （这些是未知量）。 * 支座位移： $𝐪_{s} = {\begin{matrix} q_{4} \\ q_{5} \\ q_{6} \end{matrix}} = 0$ 。 * 施加的力： $𝐏_{f} = {\begin{matrix} P_{1} \\ P_{2} \\ P_{3} \end{matrix}} = {\begin{matrix} 0 \\ 20000 \\ - 30000 \end{matrix}}$ 。 * 支座反力： $𝐏_{s} = {\begin{matrix} P_{4} \\ P_{5} \\ P_{6} \end{matrix}}$ （这些也是未知量）。

这个单个的分块方程展开为两个独立的矩阵方程：

$𝐏_{f} = 𝐊_{f f} 𝐪_{f} + 𝐊_{f s} 𝐪_{s}$
$𝐏_{s} = 𝐊_{s f} 𝐪_{f} + 𝐊_{s s} 𝐪_{s}$

求解未知位移 应用已知边界条件 $𝐪_{s} = 0$ 。这极大地简化了第一个方程： $𝐏_{f} = 𝐊_{f f} 𝐪_{f} + 𝐊_{f s} (0) ⟹ 𝐏_{f} = 𝐊_{f f} 𝐪_{f}$ 矩阵 $𝐊_{f f}$ 是全局 $𝐊$ 矩阵左上角的 3×3 分块，对应于自由自由度 $q_{1}$ 、 $q_{2}$ 和 $q_{3}$ 。现在，我们可以通过对这个小规模非奇异矩阵求逆 ¹ 来求解未知位移 $𝐪_{f}$ ： $𝐪_{f} = 𝐊_{f f}^{- 1} 𝐏_{f}$ 此计算将给出节点 3 水平和垂直位移的数值。

求解支座反力 一旦位移向量 $𝐪_{f}$ 已知，我们可以利用第二个方程（同样利用 $𝐪_{s} = 0$ ）求出未知的支座反力： $𝐏_{s} = 𝐊_{s f} 𝐪_{f} + 𝐊_{s s} (0) ⟹ 𝐏_{s} = 𝐊_{s f} 𝐪_{f}$ 最后一步是简单的矩阵-向量乘法，得出支座（节点 1 和节点 2）的全部四个反力。至此，桁架分析完成。

以下代码展示了所有计算过程。它表明 $𝐏_{s} = {\begin{matrix} 30 \\ - 50 \\ 30 \end{matrix}} k N 。$

import numpy as np
# --------------------------------------
# Step 1: Input node coordinates
# --------------------------------------
# Example geometry: a triangle truss (3 nodes, 2 elements)
nodes = {
    1: [0.0, 0.0],
    2: [0.0, 2.0],
    3: [2.0, 2.0]
}
# --------------------------------------
# Step 2: Define DOF mapping (ID matrix)
# --------------------------------------
# 4 rows (u1, v1, u2, v2) per element
# Each column corresponds to one element
# Example: two elements (1–2) and (2–3)
ID = np.array([
    [4, 5, 4],   # u1 global DOF for elem1, elem2
    [1, 6, 1],   # v1
    [5, 2, 2],   # u2
    [6, 3, 3]    # v2
])
# Number of global DOF
num_dof = int(np.max(ID))
Ks = np.zeros((num_dof, num_dof))  # global stiffness matrix initialized to zero
# --------------------------------------
# Step 3: Define truss element stiffness
# --------------------------------------
def truss_element_stiffness(E, A, node1, node2, nodes):
    """
    Compute 4x4 stiffness matrix for a 2D truss element in global coordinates.
    Parameters
    ----------
    E : float
        Young's modulus
    A : float
        Cross-sectional area
    node1, node2 : int
        Node IDs of the element
    nodes : dict
        Dictionary of node coordinates, e.g. {1: [x1, y1], 2: [x2, y2], ...}
    """
    x1, y1 = nodes[node1]
    x2, y2 = nodes[node2]
    L = np.sqrt((x2 - x1)**2 + (y2 - y1)**2)
    c = (x2 - x1) / L
    s = (y2 - y1) / L
    k = (E * A / L) * np.array([
        [ c*c,  c*s, -c*c, -c*s],
        [ c*s,  s*s, -c*s, -s*s],
        [-c*c, -c*s,  c*c,  c*s],
        [-c*s, -s*s,  c*s,  s*s]
    ])
    return k
# --------------------------------------
# Step 4: Assembly function (Python version of ADDK)
# --------------------------------------
def ADDK(Ks, Ke, ID, elem_id):
    """
    Assemble element stiffness matrix Ke into global Ks using mapping ID.
    - ID: 4 x N array (u1,v1,u2,v2 per element)
    - elem_id: which element (column index)
    """
    n = 4  # local DOF per truss element
    for i in range(n):
        for j in range(n):
            ig = ID[i, elem_id]
            jg = ID[j, elem_id]
            if ig != 0 and jg != 0:
                val = Ke[i, j]
                Ks[ig - 1, jg - 1] += val
    return Ks
    
# --------------------------------------
# Step 5: Example assembly
# --------------------------------------
E = 100e9  # Young's modulus (Pa)
A = 2e-4   # cross-sectional area (m^2)
# Element 1:
Ke1 = truss_element_stiffness(E, A, 1, 2, nodes)
Ks = ADDK(Ks, Ke1, ID, elem_id=0)
# Element 2: Node 2 → Node 3
Ke2 = truss_element_stiffness(E, A, 2, 3, nodes)
Ks = ADDK(Ks, Ke2, ID, elem_id=1)
# Element 3: Node 2 → Node 3
Ke2 = truss_element_stiffness(E, A, 1, 3, nodes)
Ks = ADDK(Ks, Ke2, ID, elem_id=2)
# --------------------------------------
# Step 6: Print the result
# --------------------------------------
np.set_printoptions(precision=3, suppress=True)
print("Global stiffness matrix Ks: (1e7)")
print(Ks/1e7)
Kff=Ks[0:3,0:3]
Kfs=Ks[0:3, 3:6]
Ksf=Ks[3:6,0:3]
Kss=Ks[3:6,3:6]
Pf=np.array([[0],[20e3],[-30e3]])
qf = np.linalg.inv(Kff) @ Pf
print(qf)
Ps =  Ksf @ qf
Ps

array([[ 30000.],
       [-50000.],
       [ 30000.]])

求解 $𝐏_{f}$ 的方法中，存在比求逆 $𝐊_{f f}$ 更优且计算成本更低的方式。↩︎

局部刚度矩阵（Klocal）