主题006：刚度矩阵与阻尼矩阵构建

1. 引言

在结构动力学中，质量矩阵、刚度矩阵和阻尼矩阵是描述系统动力特性的三个核心要素。上一主题我们详细讨论了质量矩阵的构建方法，本主题将重点介绍刚度矩阵和阻尼矩阵的构建技术。

刚度矩阵描述了结构的弹性恢复特性，决定了系统在变形时产生的恢复力。阻尼矩阵则描述了系统的能量耗散机制，反映了振动过程中能量的衰减特性。这两个矩阵与质量矩阵一起，构成了结构动力学的完整数学模型。

本主题将从刚度矩阵的基本概念出发，详细介绍各类结构元件的刚度矩阵推导方法，然后讨论阻尼矩阵的各种模型及其构建技术，最后通过丰富的工程案例展示如何在Python中实现这些矩阵的构建与仿真。

2. 刚度矩阵的基本概念

2.1 刚度矩阵的物理意义

刚度矩阵K是描述结构弹性特性的核心矩阵，其元素$k_{ij}$表示在第j个自由度产生单位位移时，在第i个自由度上需要施加的力。这一定义来源于胡克定律的推广形式。

刚度矩阵的一般形式为：

$$K=\left[\begin{array}{cccc}{k}_{11}& k_{12}& \cdots & k_{1n}\\ k_{21}& k_{22}& \cdots & k_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ k_{n1}& k_{n2}& \cdots & k_{nn}\end{array}\right]$$

对角元素$k_{ii}$称为直接刚度或主刚度，表示该自由度方向的刚度。非对角元素$k_{ij}$（$i\ne j$）称为耦合刚度，表示不同自由度之间的力-位移耦合关系。

2.2 刚度矩阵的性质

刚度矩阵具有以下重要性质：

对称性：$K=K^T$，即$k_{ij}=k_{ji}$。这一性质来源于功的互等定理（Maxwell-Betti互等定理），是线弹性系统的基本特征。

正定性或半正定性：

• 对于约束充分的结构，K是正定矩阵，即对于任意非零向量$v$，有$v^{T}Kv>0$
• 对于存在刚体位移的结构，K是半正定矩阵，即$v^{T}Kv\ge 0$

稀疏性：对于大型结构，刚度矩阵通常是稀疏的，非零元素主要集中在对角线附近（带状矩阵）。

坐标依赖性：刚度矩阵的具体形式依赖于广义坐标的选择。

2.3 刚度矩阵与势能

刚度矩阵与系统的弹性势能密切相关。对于线弹性系统，势能可以表示为：

$$U=\frac{1}{2}{u}^{T}Ku$$

其中$u$是位移向量。这个二次型表达式表明，势能是位移的二次函数，刚度矩阵决定了势能曲面的形状。

根据Castigliano定理，广义力可以通过势能对广义位移的偏导得到：

$$F_i=\frac{\partial U}{\partial u_i}=\sum _j{k}_{ij}{u}_j$$

这正是刚度矩阵的定义。

3. 单元刚度矩阵

3.1 轴向杆单元

考虑一个均匀轴向杆单元，长度为L，截面积为A，弹性模量为E。单元有两个节点，每个节点有一个轴向自由度。

根据材料力学，杆的轴向力-位移关系为：

$$F=\frac{EA}{L}\Delta$$

其中$\Delta$是轴向变形。

对于杆单元，设节点1和节点2的轴向位移分别为$u_1$和$u_2$，则单元刚度矩阵为：

$$K_e=\frac{EA}{L}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]$$

这个矩阵的物理意义：

• $k_{11}=EA/L$：节点1产生单位位移时，在节点1需要施加的力
• $k_{12}=-EA/L$：节点2产生单位位移时，在节点1产生的力（负号表示方向相反）

3.2 扭转杆单元

对于承受扭矩的圆截面杆，扭转刚度矩阵与轴向杆类似：

$$K_e=\frac{GJ}{L}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]$$

其中G是剪切模量，J是极惯性矩。

3.3 Euler-Bernoulli梁单元

梁单元是结构动力学中最重要的单元类型之一。考虑横向弯曲的Euler-Bernoulli梁单元，每个节点有两个自由度：横向位移$v$和转角$\theta$。

梁单元的刚度矩阵为4×4矩阵：

$$K_e=\frac{EI}{L^3}\left[\begin{array}{cccc}12& 6L& -12& 6L\\ 6L& 4L^2& -6L& 2L^2\\ -12& -6L& 12& -6L\\ 6L& 2L^2& -6L& 4L^2\end{array}\right]$$

其中：

• 第1、3行/列对应横向位移
• 第2、4行/列对应转角

这个矩阵的推导基于Euler-Bernoulli梁理论，假设：

1. 变形前垂直于中性轴的平面，变形后仍保持平面且垂直于中性轴
2. 忽略剪切变形
3. 小变形假设

3.4 Timoshenko梁单元

当梁的高跨比较大时，剪切变形不可忽略，需要采用Timoshenko梁理论。Timoshenko梁单元的刚度矩阵考虑了剪切变形的影响：

$$K_e=K_b+K_s$$

其中$K_b$是弯曲刚度矩阵，$K_s$是剪切刚度矩阵。

Timoshenko梁单元更适合描述短粗梁和复合材料梁的动力特性。

3.5 平面桁架单元

平面桁架单元是轴向杆单元在平面内的推广。单元与x轴的夹角为$\theta$，局部坐标系下的刚度矩阵需要转换到整体坐标系：

$$K_e^{(global)}=T^T{K}_e^{(local)}T$$

其中变换矩阵T为：

$$T=\left[\begin{array}{cccc}\cos \theta & \sin \theta & 0& 0\\ -\sin \theta & \cos \theta & 0& 0\\ 0& 0& \cos \theta & \sin \theta \\ 0& 0& -\sin \theta & \cos \theta \end{array}\right]$$

3.6 平面框架单元

平面框架单元同时考虑轴向变形和弯曲变形，是轴向杆单元和梁单元的组合。每个节点有三个自由度：两个平移和一个转动。

局部坐标系下的刚度矩阵为6×6矩阵：

$$K_e^{(local)}=\left[\begin{array}{cccccc}\frac{EA}{L}& 0& 0& -\frac{EA}{L}& 0& 0\\ 0& \frac{12EI}{L^3}& \frac{6EI}{L^2}& 0& -\frac{12EI}{L^3}& \frac{6EI}{L^2}\\ 0& \frac{6EI}{L^2}& \frac{4EI}{L}& 0& -\frac{6EI}{L^2}& \frac{2EI}{L}\\ -\frac{EA}{L}& 0& 0& \frac{EA}{L}& 0& 0\\ 0& -\frac{12EI}{L^3}& -\frac{6EI}{L^2}& 0& \frac{12EI}{L^3}& -\frac{6EI}{L^2}\\ 0& \frac{6EI}{L^2}& \frac{2EI}{L}& 0& -\frac{6EI}{L^2}& \frac{4EI}{L}\end{array}\right]$$

4. 整体刚度矩阵的组装

4.1 直接刚度法

整体刚度矩阵的组装采用直接刚度法（Direct Stiffness Method），基本步骤如下：

1. 单元分析：计算各单元在局部坐标系下的刚度矩阵
2. 坐标变换：将单元刚度矩阵转换到整体坐标系
3. 对号入座：根据单元节点编号，将单元刚度矩阵的元素叠加到整体刚度矩阵的对应位置
4. 边界条件处理：施加约束条件，消除刚体位移

4.2 组装算法

整体刚度矩阵的组装算法：

初始化整体刚度矩阵 K_global = 0
对于每个单元 e:
    计算单元刚度矩阵 K_e（局部坐标系）
    如果单元方向与整体坐标系不一致:
        计算坐标变换矩阵 T
        K_e = T^T * K_e * T
    获取单元节点编号 [i, j, ...]
    对于 K_e 中的每个元素 K_e[p][q]:
        确定在整体矩阵中的位置 [I, J]
        K_global[I][J] += K_e[p][q]

4.3 边界条件的处理

在组装完成后，需要施加边界条件来约束结构的刚体位移。常用方法：

直接法：直接删除约束自由度对应的行和列

• 优点：简单直接，减小矩阵规模
• 缺点：需要重新编号自由度

置大数法：在约束自由度对应的对角元素上置一个很大的数（如$10^{20}$）

• 优点：实现简单，保持矩阵结构
• 缺点：可能引入数值误差

罚函数法：通过罚函数将约束条件引入能量泛函

• 优点：数学基础严谨
• 缺点：需要选择适当的罚参数

拉格朗日乘子法：引入拉格朗日乘子处理约束

• 优点：精确满足约束条件
• 缺点：增加自由度数目

5. 弹簧-质量系统的刚度矩阵

5.1 串联弹簧系统

考虑n个质量通过弹簧串联的系统。设第i个弹簧的刚度为$k_i$，连接质量$m_{i-1}$和$m_i$。

刚度矩阵为三对角矩阵：

$$K=\left[\begin{array}{ccccc}{k}_1+k_2& -k_2& 0& \cdots & 0\\ -k_2& k_2+k_3& -k_3& \cdots & 0\\ 0& -k_3& k_3+k_4& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & k_n\end{array}\right]$$

5.2 并联弹簧系统

当多个弹簧并联连接两个质量时，等效刚度为各弹簧刚度之和：

$$k_{eq}=\sum _i{k}_i$$

5.3 复杂弹簧网络

对于复杂的弹簧网络，可以通过以下步骤构建刚度矩阵：

1. 识别所有自由度（各质量块的位移）
2. 对每个自由度，分析与其相连的弹簧
3. 对角元素$k_{ii}$为连接到节点i的所有弹簧刚度之和
4. 非对角元素$k_{ij}$为连接节点i和j的弹簧刚度之和的负值

6. 阻尼矩阵的基本概念

6.1 阻尼的物理机制

阻尼是结构振动过程中能量耗散的机制，主要来源包括：

材料阻尼：材料内部摩擦引起的能量耗散，与材料的粘弹性特性有关

结构阻尼：连接部位的摩擦、微滑移等引起的能量耗散

流体阻尼：结构在空气或液体中运动时受到的阻力

辐射阻尼：振动能量以波的形式向外传播而耗散

6.2 阻尼矩阵的定义

阻尼矩阵C描述了阻尼力与速度之间的关系：

$$F_d=C\dot{u}$$

其中$F_d$是阻尼力向量，$\dot{u}$是速度向量。

6.3 阻尼矩阵的性质

阻尼矩阵通常具有以下性质：

对称性：$C=C^T$，对于大多数实际结构成立

正定性：${\dot{u}}^{T}C\dot{u}>0$，保证能量耗散

稀疏性：与刚度矩阵类似，通常是稀疏矩阵

7. 阻尼模型

7.1 粘性阻尼模型

粘性阻尼模型假设阻尼力与速度成正比：

$$F_d=c\dot{u}$$

其中c是粘性阻尼系数。这是最简单的阻尼模型，也是结构动力学中最常用的模型。

粘性阻尼的优点：

• 数学处理简单
• 与线性微分方程理论兼容
• 便于实验测定

7.2 瑞利阻尼（Rayleigh Damping）

瑞利阻尼假设阻尼矩阵是质量矩阵和刚度矩阵的线性组合：

$$C=\alpha M+\beta K$$

其中$\alpha$和$\beta$是瑞利阻尼系数。

瑞利阻尼的优点：

• 自动满足正交性条件（与模态振型正交）
• 只需确定两个参数
• 便于模态分析

瑞利阻尼系数与模态阻尼比的关系：

$${\zeta }_i=\frac{\alpha }{2{\omega }_i}+\frac{\beta {\omega }_i}{2}$$

通过指定两个模态的阻尼比，可以求解$\alpha$和$\beta$。

7.3 模态阻尼

模态阻尼假设各模态之间没有阻尼耦合，阻尼矩阵在模态坐标下是对角阵：

$${\Phi }^{T}C\Phi =\text{diag}(2{\zeta }_1{\omega }_1,2{\zeta }_2{\omega }_2,\dots ,2{\zeta }_n{\omega }_n)$$

其中$\Phi$是模态矩阵。

模态阻尼的优点：

• 物理意义明确
• 便于实验测定
• 与模态叠加法兼容

7.4 结构阻尼（滞回阻尼）

结构阻尼假设阻尼力与位移成正比，但相位滞后90度：

$$F_d=i\eta Ku$$

其中$\eta$是结构阻尼系数（损耗因子），$i$是虚数单位。

结构阻尼更适合描述材料内部摩擦引起的能量耗散。

7.5 库仑阻尼

库仑阻尼（干摩擦阻尼）假设阻尼力的大小恒定，方向与速度相反：

$$F_d=-\mu N\cdot \text{sgn}(\dot{u})$$

其中$\mu$是摩擦系数，N是正压力。

库仑阻尼是非线性阻尼，分析较为复杂。

7.6 非线性阻尼

实际工程中，阻尼往往具有非线性特性。常见的非线性阻尼模型包括：

幂律阻尼：$F_d=c|\dot{u}{|}^n\cdot \text{sgn}(\dot{u})$

分段线性阻尼：在不同速度范围采用不同的阻尼系数

流体动力阻尼：考虑流体动力学效应的复杂阻尼模型

8. 阻尼矩阵的构建方法

8.1 基于实验数据的阻尼矩阵

通过模态试验可以测定各阶模态的阻尼比，然后构建阻尼矩阵：

$$C=M\Phi \text{diag}(2{\zeta }_i{\omega }_i){\Phi }^{T}M$$

其中$\Phi$是质量归一化的模态矩阵。

8.2 单元阻尼矩阵

类似于刚度矩阵，可以定义单元阻尼矩阵，然后组装整体阻尼矩阵。

对于粘性阻尼，单元阻尼矩阵通常与单元刚度矩阵成正比：

$$C_e=\beta K_e$$

8.3 比例阻尼矩阵的确定

确定瑞利阻尼系数的步骤：

1. 选择两个控制频率${\omega }_i$和${\omega }_j$（通常对应关心的频率范围）
2. 指定这两个频率处的阻尼比${\zeta }_i$和${\zeta }_j$
3. 求解方程组：

$$\left[\begin{array}{cc}\frac{1}{2{\omega }_i}& \frac{{\omega }_i}{2}\\ \frac{1}{2{\omega }_j}& \frac{{\omega }_j}{2}\end{array}\right]\left[\begin{array}{c}\alpha \\ \beta \end{array}\right]=\left[\begin{array}{c}{\zeta }_i\\ {\zeta }_j\end{array}\right]$$

解得：

$$\alpha =\frac{2{\omega }_i{\omega }_j({\zeta }_i{\omega }_j-{\zeta }_j{\omega }_i)}{{\omega }_j^2-{\omega }_i^2}$$

$$\beta =\frac{2({\zeta }_j{\omega }_j-{\zeta }_i{\omega }_i)}{{\omega }_j^2-{\omega }_i^2}$$

9. 数值实现与Python编程

9.1 单元刚度矩阵函数

import numpy as np

def bar_stiffness_matrix(E, A, L):
    """
    计算轴向杆单元刚度矩阵
    
    参数:
        E: 弹性模量 (Pa)
        A: 截面积 (m²)
        L: 杆长 (m)
    
    返回:
        K: 2×2刚度矩阵
    """
    k = E * A / L
    return k * np.array([[1, -1], [-1, 1]])

def beam_stiffness_matrix(E, I, L):
    """
    计算Euler-Bernoulli梁单元刚度矩阵
    
    参数:
        E: 弹性模量 (Pa)
        I: 截面惯性矩 (m⁴)
        L: 梁长 (m)
    
    返回:
        K: 4×4刚度矩阵
    """
    k = E * I / L**3
    return k * np.array([
        [12, 6*L, -12, 6*L],
        [6*L, 4*L**2, -6*L, 2*L**2],
        [-12, -6*L, 12, -6*L],
        [6*L, 2*L**2, -6*L, 4*L**2]
    ])

def truss_stiffness_matrix(E, A, L, theta):
    """
    计算平面桁架单元刚度矩阵（整体坐标系）
    
    参数:
        E: 弹性模量 (Pa)
        A: 截面积 (m²)
        L: 杆长 (m)
        theta: 与x轴夹角 (rad)
    
    返回:
        K: 4×4刚度矩阵
    """
    k = E * A / L
    c = np.cos(theta)
    s = np.sin(theta)
    
    return k * np.array([
        [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]
    ])

9.2 整体刚度矩阵组装

def assemble_stiffness_matrix(n_nodes, elements, element_stiffness_func):
    """
    组装整体刚度矩阵
    
    参数:
        n_nodes: 节点总数
        elements: 单元列表，每个元素为(node_i, node_j, *params)
        element_stiffness_func: 计算单元刚度矩阵的函数
    
    返回:
        K: 整体刚度矩阵
    """
    ndof = n_nodes  # 假设每个节点1个自由度
    K = np.zeros((ndof, ndof))
    
    for elem in elements:
        i, j = elem[0], elem[1]  # 节点编号
        params = elem[2:]  # 单元参数
        
        # 计算单元刚度矩阵
        K_e = element_stiffness_func(*params)
        
        # 组装到整体矩阵
        K[i, i] += K_e[0, 0]
        K[i, j] += K_e[0, 1]
        K[j, i] += K_e[1, 0]
        K[j, j] += K_e[1, 1]
    
    return K

9.3 阻尼矩阵构建

def rayleigh_damping_matrix(M, K, alpha, beta):
    """
    构建瑞利阻尼矩阵
    
    参数:
        M: 质量矩阵
        K: 刚度矩阵
        alpha: 质量比例系数
        beta: 刚度比例系数
    
    返回:
        C: 阻尼矩阵
    """
    return alpha * M + beta * K

def calculate_rayleigh_coefficients(omega_i, omega_j, zeta_i, zeta_j):
    """
    计算瑞利阻尼系数
    
    参数:
        omega_i, omega_j: 两个控制频率 (rad/s)
        zeta_i, zeta_j: 对应的阻尼比
    
    返回:
        alpha, beta: 瑞利阻尼系数
    """
    A = np.array([[1/(2*omega_i), omega_i/2],
                  [1/(2*omega_j), omega_j/2]])
    b = np.array([zeta_i, zeta_j])
    
    coeffs = np.linalg.solve(A, b)
    return coeffs[0], coeffs[1]

10. 案例研究

10.1 三自由度弹簧-质量系统

考虑一个三自由度弹簧-质量系统，参数如下：

• 质量：$m_1=2$ kg, $m_2=1.5$ kg, $m_3=1$ kg
• 弹簧刚度：$k_1=100$ N/m, $k_2=150$ N/m, $k_3=100$ N/m, $k_4=80$ N/m

质量矩阵：

$$M=\left[\begin{array}{ccc}2& 0& 0\\ 0& 1.5& 0\\ 0& 0& 1\end{array}\right]$$

刚度矩阵：

$$K=\left[\begin{array}{ccc}{k}_1+k_2& -k_2& 0\\ -k_2& k_2+k_3& -k_3\\ 0& -k_3& k_3+k_4\end{array}\right]=\left[\begin{array}{ccc}250& -150& 0\\ -150& 250& -100\\ 0& -100& 180\end{array}\right]$$

阻尼矩阵（瑞利阻尼，假设$\alpha =0.1$, $\beta =0.01$）：

$$C=0.1M+0.01K$$

10.2 简支梁的刚度矩阵

考虑一个简支梁，离散为3个单元。每个单元长度L=1m，EI=1000 N·m²。

单元刚度矩阵：

$$K_e=\frac{1000}{1^3}\left[\begin{array}{cccc}12& 6& -12& 6\\ 6& 4& -6& 2\\ -12& -6& 12& -6\\ 6& 2& -6& 4\end{array}\right]$$

组装后的整体刚度矩阵（考虑简支边界条件）为8×8矩阵。

10.3 平面桁架结构

考虑一个简单的平面桁架，由3根杆件组成。通过坐标变换和组装，得到整体刚度矩阵。

11. 高级主题

11.1 几何非线性刚度矩阵

对于大变形问题，需要考虑几何非线性效应。切线刚度矩阵包括：

$$K_T=K_E+K_G$$

其中$K_E$是弹性刚度矩阵，$K_G$是几何刚度矩阵（初应力刚度矩阵）。

几何刚度矩阵考虑了当前应力状态对结构刚度的影响，在屈曲分析和后屈曲分析中非常重要。

11.2 刚度矩阵的缩聚

对于大型结构，可以通过静力缩聚或动力缩聚减少自由度：

Guyan缩聚：

$$K_r=K_{mm}-K_{ms}{K}_{ss}^{-1}{K}_{sm}$$

其中下标m表示主自由度，s表示从自由度。

11.3 复刚度矩阵

对于频域分析，复刚度矩阵考虑了阻尼的频变特性：

$$K^{\ast }(\omega )=K+i\omega C-{\omega }^{2}M$$

复刚度矩阵的实部代表储能特性，虚部代表耗能特性。

11.4 随机结构的刚度矩阵

当结构参数具有随机性时，刚度矩阵也是随机的。可以采用：

• 蒙特卡洛模拟：大量随机抽样计算
• 随机有限元法：基于摄动理论的分析方法
• 谱分析方法：基于随机场的展开方法

12. 工程应用实例

12.1 高层建筑结构的刚度建模

某30层框架-剪力墙结构：

框架部分：

• 梁单元模拟框架梁和柱
• 考虑轴向、弯曲和剪切变形

剪力墙部分：

• 平面应力单元或壳单元
• 考虑面内和面外刚度

连接部分：

• 刚性楼板假设
• 考虑节点刚域效应

12.2 桥梁结构的阻尼建模

大跨度斜拉桥的阻尼特性：

结构阻尼：

• 钢箱梁：$\zeta =0.5\%-1\%$
• 混凝土桥塔：$\zeta =1\%-2\%$

附加阻尼：

• 粘滞阻尼器
• 调谐质量阻尼器(TMD)
• 调谐液体阻尼器(TLD)

12.3 旋转机械的刚度与阻尼

汽轮机转子系统的建模：

轴段刚度：

• 考虑剪切变形和转动惯量
• 陀螺效应的影响

轴承刚度与阻尼：

• 油膜轴承的动力特性系数
• 刚度系数$k_{xx},k_{xy},k_{yx},k_{yy}$
• 阻尼系数$c_{xx},c_{xy},c_{yx},c_{yy}$

13. 误差分析与验证

13.1 刚度矩阵的精度评估

离散化误差：

• 单元尺寸的影响
• 形函数阶次的影响
• 收敛性分析

边界条件误差：

• 约束条件的合理施加
• 弹性支承的模拟

13.2 阻尼模型的验证

模态阻尼比的测定：

• 半功率带宽法
• 对数衰减法
• 随机减量法

阻尼矩阵的修正：

• 基于试验数据的修正
• 模型更新技术

14. 总结与展望

14.1 主要内容回顾

本主题系统介绍了刚度矩阵和阻尼矩阵的构建方法：

1. 刚度矩阵基本概念：定义、性质和物理意义
2. 单元刚度矩阵：杆单元、梁单元、桁架单元等
3. 整体刚度矩阵组装：直接刚度法和边界条件处理
4. 阻尼矩阵基本概念：阻尼机制和阻尼模型
5. 阻尼模型：粘性阻尼、瑞利阻尼、模态阻尼等
6. 数值实现：Python编程实现
7. 工程应用：高层建筑、桥梁、旋转机械等

14.2 关键知识点

• 刚度矩阵是对称矩阵，描述结构的弹性特性
• 阻尼矩阵描述能量耗散机制
• 瑞利阻尼是最常用的阻尼模型
• 单元刚度矩阵通过坐标变换和组装形成整体刚度矩阵
• 边界条件的正确处理对分析结果至关重要

14.3 后续学习方向

掌握了刚度矩阵和阻尼矩阵构建方法后，下一步将学习：

1. 模态分析：求解固有频率和振型
2. 时程响应分析：计算动力响应
3. 频域分析：频率响应函数
4. 非线性分析：几何非线性和材料非线性
5. 随机振动分析：随机荷载作用下的响应

15. 练习题

15.1 理论题

1. 解释刚度矩阵的物理意义，为什么刚度矩阵是对称的？

2. 推导Euler-Bernoulli梁单元的刚度矩阵。

3. 比较瑞利阻尼和模态阻尼的优缺点。

4. 说明刚度矩阵缩聚的基本原理。

5. 如何处理结构中的刚体位移模式？

15.2 计算题

1. 一个两自由度系统，弹簧刚度$k_1=100$ N/m, $k_2=150$ N/m, $k_3=100$ N/m。构建其刚度矩阵。

2. 一个悬臂梁长L=2m，EI=1000 N·m²，离散为2个单元。计算整体刚度矩阵。

3. 已知前两阶固有频率为${\omega }_1=10$ rad/s, ${\omega }_2=30$ rad/s，对应的阻尼比为${\zeta }_1=0.02$, ${\zeta }_2=0.03$。计算瑞利阻尼系数。

15.3 编程题

1. 编写Python函数，实现平面框架单元的刚度矩阵计算。

2. 编写程序，实现整体刚度矩阵的自动组装。

3. 编写程序，根据模态试验数据构建模态阻尼矩阵。

16. 参考文献

1. Chopra, A.K. (2019). Dynamics of Structures: Theory and Applications to Earthquake Engineering, 5th Edition. Pearson.

2. Clough, R.W. and Penzien, J. (2003). Dynamics of Structures, 3rd Edition. Computers & Structures, Inc.

3. Bathe, K.J. (2014). Finite Element Procedures, 2nd Edition. Klaus-Jürgen Bathe.

4. Craig, R.R. and Kurdila, A.J. (2006). Fundamentals of Structural Dynamics, 2nd Edition. Wiley.

5. 刘晶波, 杜修力. (2005). 结构动力学. 机械工业出版社.

6. 克拉夫, R.W., 彭津, J. (2017). 结构动力学 (第二版修订版). 高等教育出版社.

---

附录：Python完整代码

以下是本主题所有案例的完整Python实现代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eig, eigh
import matplotlib.animation as animation
from matplotlib.patches import Rectangle, Circle, FancyBboxPatch, FancyArrowPatch
import matplotlib.patches as mpatches

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

def bar_stiffness_matrix(E, A, L):
    """计算轴向杆单元刚度矩阵"""
    k = E * A / L
    return k * np.array([[1, -1], [-1, 1]])

def beam_stiffness_matrix(E, I, L):
    """计算Euler-Bernoulli梁单元刚度矩阵"""
    k = E * I / L**3
    return k * np.array([
        [12, 6*L, -12, 6*L],
        [6*L, 4*L**2, -6*L, 2*L**2],
        [-12, -6*L, 12, -6*L],
        [6*L, 2*L**2, -6*L, 4*L**2]
    ])

def truss_stiffness_matrix(E, A, L, theta):
    """计算平面桁架单元刚度矩阵（整体坐标系）"""
    k = E * A / L
    c = np.cos(theta)
    s = np.sin(theta)
    return k * np.array([
        [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]
    ])

def spring_system_stiffness(k_list, n):
    """构建串联弹簧-质量系统刚度矩阵"""
    K = np.zeros((n, n))
    for i in range(n):
        if i == 0:
            K[i, i] = k_list[0] + k_list[1] if len(k_list) > 1 else k_list[0]
            if len(k_list) > 1:
                K[i, i+1] = -k_list[1]
                K[i+1, i] = -k_list[1]
        elif i < n - 1:
            K[i, i] = k_list[i] + k_list[i+1]
            K[i, i-1] = -k_list[i]
            K[i-1, i] = -k_list[i]
            K[i, i+1] = -k_list[i+1]
            K[i+1, i] = -k_list[i+1]
        else:
            K[i, i] = k_list[i]
            K[i, i-1] = -k_list[i]
            K[i-1, i] = -k_list[i]
    return K

def rayleigh_damping_matrix(M, K, alpha, beta):
    """构建瑞利阻尼矩阵"""
    return alpha * M + beta * K

def calculate_rayleigh_coefficients(omega_i, omega_j, zeta_i, zeta_j):
    """计算瑞利阻尼系数"""
    A = np.array([[1/(2*omega_i), omega_i/2],
                  [1/(2*omega_j), omega_j/2]])
    b = np.array([zeta_i, zeta_j])
    coeffs = np.linalg.solve(A, b)
    return coeffs[0], coeffs[1]

# 案例1：三自由度弹簧-质量系统
print("=" * 60)
print("案例1：三自由度弹簧-质量系统")
print("=" * 60)

m = np.array([2.0, 1.5, 1.0])
k = np.array([100, 150, 100, 80])

M = np.diag(m)
K = spring_system_stiffness(k, 3)

print("\n质量矩阵 M:")
print(M)
print("\n刚度矩阵 K:")
print(K)

# 瑞利阻尼
alpha, beta = 0.1, 0.01
C = rayleigh_damping_matrix(M, K, alpha, beta)
print(f"\n瑞利阻尼矩阵 (alpha={alpha}, beta={beta}):")
print(C)

# 求解特征值问题
eigenvalues, eigenvectors = eig(K, M)
idx = np.argsort(eigenvalues)
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

omega = np.sqrt(np.real(eigenvalues))
f = omega / (2 * np.pi)

print("\n固有频率:")
for i in range(3):
    print(f"  ω{i+1} = {omega[i]:.4f} rad/s, f{i+1} = {f[i]:.4f} Hz")

# 案例2：简支梁分析
print("\n" + "=" * 60)
print("案例2：简支梁刚度矩阵分析")
print("=" * 60)

E = 2e11  # 钢弹性模量
I = 8.33e-6  # 矩形截面0.1x0.2m的惯性矩
L_total = 3.0  # 梁总长
n_elem = 3
L = L_total / n_elem

print(f"\n梁参数: E={E/1e9:.1f} GPa, I={I*1e6:.2f}e-6 m^4, L={L_total}m")
print(f"离散为 {n_elem} 个单元，每个单元长 {L:.2f}m")

# 单元刚度矩阵
K_e = beam_stiffness_matrix(E, I, L)
print("\n单元刚度矩阵 (kN/m, kN, kN·m):")
print(K_e / 1000)

# 组装整体刚度矩阵（简化，只考虑横向位移）
# 简支梁边界条件：两端位移为0，转角自由
n_nodes = n_elem + 1
ndof = 2 * n_nodes  # 每个节点2个自由度

K_global = np.zeros((ndof, ndof))

for i in range(n_elem):
    # 单元自由度映射到整体自由度
    dof_i = [2*i, 2*i+1, 2*(i+1), 2*(i+1)+1]
    for p in range(4):
        for q in range(4):
            K_global[dof_i[p], dof_i[q]] += K_e[p, q]

print("\n整体刚度矩阵 (部分):")
print("矩阵维度:", K_global.shape)

# 应用边界条件（简支：两端位移为0）
# 删除第0行/列和第2*(n_nodes-1)行/列（两端位移约束）
constrained_dofs = [0, 2*(n_nodes-1)]  # 两端位移约束
free_dofs = [i for i in range(ndof) if i not in constrained_dofs]
K_reduced = K_global[np.ix_(free_dofs, free_dofs)]

print(f"\n应用边界条件后刚度矩阵维度: {K_reduced.shape}")

# 案例3：瑞利阻尼系数计算
print("\n" + "=" * 60)
print("案例3：瑞利阻尼系数计算")
print("=" * 60)

# 假设前两阶固有频率和阻尼比
omega_1, omega_2 = 10.0, 30.0  # rad/s
zeta_1, zeta_2 = 0.02, 0.03

alpha, beta = calculate_rayleigh_coefficients(omega_1, omega_2, zeta_1, zeta_2)

print(f"\n控制频率: ω1={omega_1} rad/s, ω2={omega_2} rad/s")
print(f"目标阻尼比: ζ1={zeta_1*100}%, ζ2={zeta_2*100}%")
print(f"\n计算得到的瑞利阻尼系数:")
print(f"  α = {alpha:.6f}")
print(f"  β = {beta:.6f}")

# 验证阻尼比
omega_range = np.linspace(5, 50, 100)
zeta_range = alpha / (2 * omega_range) + beta * omega_range / 2

print("\n各阶模态阻尼比验证:")
for i in range(3):
    zeta_i = alpha / (2 * omega[i]) + beta * omega[i] / 2
    print(f"  模态 {i+1} (ω={omega[i]:.2f} rad/s): ζ = {zeta_i*100:.2f}%")

# 可视化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1：刚度矩阵结构
ax1 = axes[0, 0]
im1 = ax1.imshow(K, cmap='Reds', aspect='auto')
ax1.set_title('刚度矩阵结构', fontsize=12)
ax1.set_xlabel('自由度')
ax1.set_ylabel('自由度')
for i in range(3):
    for j in range(3):
        text = ax1.text(j, i, f'{K[i, j]:.0f}', ha="center", va="center", color="black")
plt.colorbar(im1, ax=ax1)

# 图2：阻尼矩阵结构
ax2 = axes[0, 1]
im2 = ax2.imshow(C, cmap='Blues', aspect='auto')
ax2.set_title('瑞利阻尼矩阵结构', fontsize=12)
ax2.set_xlabel('自由度')
ax2.set_ylabel('自由度')
for i in range(3):
    for j in range(3):
        text = ax2.text(j, i, f'{C[i, j]:.2f}', ha="center", va="center", color="black", fontsize=8)
plt.colorbar(im2, ax=ax2)

# 图3：固有频率
ax3 = axes[1, 0]
mode_numbers = np.arange(1, 4)
bars = ax3.bar(mode_numbers, f, color='steelblue', edgecolor='navy')
ax3.set_xlabel('模态阶数')
ax3.set_ylabel('固有频率 (Hz)')
ax3.set_title('系统固有频率', fontsize=12)
ax3.set_xticks(mode_numbers)
for i, bar in enumerate(bars):
    height = bar.get_height()
    ax3.text(bar.get_x() + bar.get_width()/2., height,
             f'{f[i]:.3f}', ha='center', va='bottom')

# 图4：瑞利阻尼曲线
ax4 = axes[1, 1]
ax4.plot(omega_range, zeta_range * 100, 'b-', linewidth=2, label='瑞利阻尼')
ax4.scatter([omega_1, omega_2], [zeta_1*100, zeta_2*100], color='red', s=100, zorder=5, label='控制点')
ax4.axvline(x=omega_1, color='r', linestyle='--', alpha=0.5)
ax4.axvline(x=omega_2, color='r', linestyle='--', alpha=0.5)
ax4.set_xlabel('频率 (rad/s)')
ax4.set_ylabel('阻尼比 (%)')
ax4.set_title('瑞利阻尼特性曲线', fontsize=12)
ax4.legend()
ax4.grid(True, alpha=0.3)
ax4.set_xlim([0, 50])

plt.tight_layout()
plt.savefig('stiffness_damping_analysis.png', dpi=150, bbox_inches='tight')
print("\n图片已保存: stiffness_damping_analysis.png")

# 案例4：阻尼对响应的影响
print("\n" + "=" * 60)
print("案例4：阻尼对系统响应的影响")
print("=" * 60)

# 使用Newmark-beta方法计算时程响应
def newmark_beta(M, C, K, F, dt, u0, v0, beta=0.25, gamma=0.5):
    """Newmark-beta方法求解动力响应"""
    n_dof = M.shape[0]
    n_steps = F.shape[0]
    
    u = np.zeros((n_steps, n_dof))
    v = np.zeros((n_steps, n_dof))
    a = np.zeros((n_steps, n_dof))
    
    u[0] = u0
    v[0] = v0
    
    # 初始加速度
    a[0] = np.linalg.solve(M, F[0] - C @ v0 - K @ u0)
    
    # 等效刚度矩阵
    K_eff = K + (gamma / (beta * dt)) * C + (1 / (beta * dt**2)) * M
    
    for i in range(n_steps - 1):
        # 等效荷载
        F_eff = (F[i+1] + 
                 M @ ((1/(beta*dt**2))*u[i] + (1/(beta*dt))*v[i] + (1/(2*beta)-1)*a[i]) +
                 C @ ((gamma/(beta*dt))*u[i] + (gamma/beta-1)*v[i] + dt*(gamma/(2*beta)-1)*a[i]))
        
        # 求解位移
        u[i+1] = np.linalg.solve(K_eff, F_eff)
        
        # 计算速度和加速度
        a[i+1] = (1/(beta*dt**2)) * (u[i+1] - u[i]) - (1/(beta*dt)) * v[i] - (1/(2*beta)-1) * a[i]
        v[i+1] = v[i] + (1-gamma)*dt*a[i] + gamma*dt*a[i+1]
    
    return u, v, a

# 初始条件
u0 = np.array([0.1, 0.0, 0.0])
v0 = np.zeros(3)

# 时间参数
dt = 0.01
t_max = 10
t = np.arange(0, t_max, dt)
n_steps = len(t)

# 外力（脉冲荷载）
F = np.zeros((n_steps, 3))
F[:50, 0] = 100  # 前0.5秒在第一个质量上施加100N的力

# 不同阻尼比的响应
zeta_cases = [0.0, 0.02, 0.05, 0.1]
colors = ['red', 'blue', 'green', 'orange']

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

for idx_zeta, zeta_target in enumerate(zeta_cases):
    if zeta_target == 0:
        C_case = np.zeros_like(M)
        label = '无阻尼'
    else:
        # 计算瑞利阻尼系数（基于前两阶频率）
        alpha_c = 2 * zeta_target * omega[0] * omega[1] / (omega[0] + omega[1])
        beta_c = 2 * zeta_target / (omega[0] + omega[1])
        C_case = rayleigh_damping_matrix(M, K, alpha_c, beta_c)
        label = f'ζ={zeta_target*100:.0f}%'
    
    # 计算响应
    u_resp, v_resp, a_resp = newmark_beta(M, C_case, K, F, dt, u0, v0)
    
    # 绘制第一个质量的位移响应
    ax = axes[idx_zeta // 2, idx_zeta % 2]
    ax.plot(t, u_resp[:, 0], color=colors[idx_zeta], linewidth=1.5)
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel('位移 (m)')
    ax.set_title(f'{label} - 质量1位移响应', fontsize=12)
    ax.grid(True, alpha=0.3)
    ax.set_xlim([0, t_max])

plt.tight_layout()
plt.savefig('damping_response_comparison.png', dpi=150, bbox_inches='tight')
print("\n图片已保存: damping_response_comparison.png")

# 动画演示：三自由度系统自由振动
print("\n正在生成动画...")

# 使用有阻尼情况生成动画
alpha_anim = 0.2
beta_anim = 0.005
C_anim = rayleigh_damping_matrix(M, K, alpha_anim, beta_anim)

# 初始条件
u0_anim = np.array([0.5, 0.3, 0.1])
v0_anim = np.zeros(3)

# 自由振动（无外力）
F_zero = np.zeros((n_steps, 3))
u_anim, v_anim, a_anim = newmark_beta(M, C_anim, K, F_zero, dt, u0_anim, v0_anim)

# 创建动画
fig_anim, ax_anim = plt.subplots(figsize=(12, 6))
ax_anim.set_xlim(-0.5, 3.5)
ax_anim.set_ylim(-1.0, 1.0)
ax_anim.set_xlabel('质量块位置')
ax_anim.set_ylabel('位移')
ax_anim.set_title('三自由度系统有阻尼自由振动动画')
ax_anim.grid(True, alpha=0.3)
ax_anim.axhline(y=0, color='k', linestyle='-', linewidth=0.5)

# 绘制弹簧
spring_lines = []
for i in range(4):
    line, = ax_anim.plot([], [], 'k-', linewidth=2)
    spring_lines.append(line)

# 绘制质量块
mass_patches = []
mass_colors = ['red', 'green', 'blue']
for i in range(3):
    circle = Circle((i+1, 0), 0.15, color=mass_colors[i], ec='black', linewidth=2)
    ax_anim.add_patch(circle)
    mass_patches.append(circle)

# 添加标签
for i in range(3):
    ax_anim.text(i+1, -0.4, f'm{i+1}', ha='center', fontsize=10)

# 固定端
ax_anim.plot([0, 0], [-0.5, 0.5], 'k-', linewidth=4)
ax_anim.plot([4, 4], [-0.5, 0.5], 'k-', linewidth=4)

def init():
    for line in spring_lines:
        line.set_data([], [])
    for patch in mass_patches:
        patch.center = (patch.center[0], 0)
    return spring_lines + mass_patches

def animate(frame):
    # 更新质量块位置
    positions = u_anim[frame, :]
    x_positions = np.arange(3) + 1 + positions
    
    for i, patch in enumerate(mass_patches):
        patch.center = (x_positions[i], 0)
    
    # 更新弹簧
    spring_lines[0].set_data([0, x_positions[0]], [0, 0])
    for i in range(2):
        spring_lines[i+1].set_data([x_positions[i], x_positions[i+1]], [0, 0])
    spring_lines[3].set_data([x_positions[2], 4], [0, 0])
    
    return spring_lines + mass_patches

anim = animation.FuncAnimation(fig_anim, animate, init_func=init,
                               frames=n_steps, interval=50, blit=True)

anim.save('damped_vibration_animation.gif', writer='pillow', fps=20)
print("动画已保存: damped_vibration_animation.gif")

plt.close('all')

print("\n" + "=" * 60)
print("程序运行完成！")
print("=" * 60)