主题005：多自由度系统建模与质量矩阵构建

1. 引言

在前面的主题中，我们详细研究了单自由度（SDOF）系统的动力学特性。然而，实际工程结构往往具有复杂的形态和分布质量，不能简单地用一个质点来代表。例如，高层建筑、桥梁、飞机机翼、汽车底盘等结构都需要采用多自由度（MDOF）模型才能准确描述其动力学行为。

多自由度系统是指具有多个独立运动自由度的系统，每个自由度对应一个独立的广义坐标。描述n自由度系统运动需要n个二阶微分方程，这些方程通常以矩阵形式表示，其中质量矩阵、刚度矩阵和阻尼矩阵是系统的核心特征矩阵。

本主题将重点介绍多自由度系统的建模方法，特别是质量矩阵的构建技术。我们将从离散化建模的基本概念出发，详细讨论集中质量法和一致质量法的原理与应用，并通过丰富的案例展示如何在Python中实现多自由度系统的建模与仿真。

2. 多自由度系统的基本概念

2.1 从单自由度到多自由度

单自由度系统虽然简单，但只能描述结构的整体运动。当结构的不同部分具有相对运动时，就需要引入多自由度模型。考虑一个三层剪切型框架结构，每一层都可以相对于其他层产生水平位移，这就构成了一个三自由度系统。

多自由度系统的运动方程可以表示为矩阵形式：

$$M\ddot{u}+C\dot{u}+Ku=F(t)$$

其中：

• $M$ 是质量矩阵（n×n）
• $C$ 是阻尼矩阵（n×n）
• $K$ 是刚度矩阵（n×n）
• $u$ 是位移向量（n×1）
• $F(t)$ 是外力向量（n×1）

这个方程是单自由度运动方程的自然推广，但矩阵形式使得我们能够同时描述多个自由度之间的耦合关系。

2.2 自由度的定义与选择

自由度是指完全确定系统位形所需的独立坐标数目。在多自由度系统建模中，自由度的选择至关重要，它直接影响模型的复杂度和计算精度。

自由度的类型包括：

平移自由度：描述质点在x、y、z方向的直线运动。对于平面问题，通常有2个平移自由度；对于空间问题，有3个平移自由度。

转动自由度：描述刚体绕x、y、z轴的旋转运动。在梁和板结构中，转动自由度对于准确描述变形至关重要。

广义自由度：在某些情况下，可以使用模态坐标或其他广义坐标来描述系统的运动，这可以显著减少自由度数目。

自由度的选择原则：

1. 足够描述主要振动模态：自由度数目应能捕捉到关心的振动模态
2. 考虑计算资源：过多的自由度会增加计算成本
3. 关注关键位置：在应力集中区域或响应关注区域应适当增加自由度
4. 利用对称性：对于对称结构，可以只建模一半并施加对称边界条件

2.3 多自由度系统的分类

根据质量分布特性，多自由度系统可以分为以下几类：

集中质量系统：质量集中在离散的质点上，通过无质量的结构元件连接。这种模型在高层建筑、多跨桥梁等结构中广泛应用。

分布质量系统：质量沿结构连续分布，如梁、板、壳等连续体结构。这类系统通常需要离散化处理才能进行数值分析。

混合质量系统：部分质量集中，部分质量分布。实际工程结构往往属于这种类型。

3. 质量矩阵的物理意义

3.1 质量矩阵的定义

质量矩阵M是一个对称正定矩阵，其元素$m_{ij}$表示在第j个自由度产生单位加速度时，在第i个自由度上需要施加的力。这一定义来源于达朗贝尔原理，体现了系统的惯性特性。

质量矩阵的一般形式为：

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

对角元素$m_{ii}$称为广义质量或等效质量，表示与第i个自由度相关的惯性。非对角元素$m_{ij}$（$i\ne j$）称为质量耦合项，表示不同自由度之间的惯性耦合。

3.2 质量矩阵的性质

质量矩阵具有以下重要性质：

对称性：$M=M^T$，即$m_{ij}=m_{ji}$。这一性质来源于功的互等定理，是保守系统的基本特征。

正定性：对于任意非零向量$v$，有$v^{T}Mv>0$。这意味着系统的动能总是正的，符合物理实际。

惯性特性：质量矩阵完全描述了系统的惯性分布，决定了系统在给定外力作用下的加速度响应。

坐标依赖性：质量矩阵的具体形式依赖于广义坐标的选择。不同的坐标系会导致不同的质量矩阵形式。

3.3 质量耦合的概念

当质量矩阵存在非零的非对角元素时，我们说系统存在质量耦合。质量耦合意味着一个自由度的加速度会在其他自由度上产生惯性力。

考虑一个两自由度系统，如果质量矩阵为对角阵：

$$M=\left[\begin{array}{cc}{m}_1& 0\\ 0& m_2\end{array}\right]$$

则两个自由度的运动方程是解耦的，可以独立求解。

如果质量矩阵为非对角阵：

$$M=\left[\begin{array}{cc}{m}_1& m_{12}\\ m_{21}& m_2\end{array}\right]$$

则运动方程存在耦合：

$$m_1{\ddot{u}}_1+m_{12}{\ddot{u}}_2+\cdots =F_1$$
$$m_{21}{\ddot{u}}_1+m_2{\ddot{u}}_2+\cdots =F_2$$

这种情况下，两个自由度的加速度相互影响，必须联立求解。

4. 集中质量法

4.1 集中质量法的基本原理

集中质量法（Lumped Mass Method）是最简单的质量离散化方法。其基本思想是将结构的分布质量集中到选定的节点上，形成离散的质点，而连接这些质点的结构元件则视为无质量。

集中质量法的实施步骤：

1. 选择节点：在结构的关键位置布置节点，如楼层中心、质量中心、支座等
2. 分配质量：将结构的质量按一定规则分配到各个节点
3. 确定自由度：为每个节点分配适当的自由度（平移、转动）
4. 构建质量矩阵：根据节点质量和自由度配置形成质量矩阵

4.2 质量分配方法

在集中质量法中，质量的分配需要遵循质量守恒原则。常用的分配方法包括：

几何分配法：根据节点周围的几何区域分配质量。例如，对于梁结构，可以将梁段质量平均分配给两端的节点。

形函数分配法：基于形函数的思想，将分布质量按形函数的值加权分配到节点。这种方法与有限元的一致质量矩阵有相似之处。

质心法：将构件的质量集中到其质心位置，然后将质心位置的质量分配到相邻节点。这种方法适用于复杂形状的结构。

工程经验法：根据工程经验和规范要求，将质量分配到特定节点。例如，在高层建筑中，通常将每层的质量集中到楼层中心。

4.3 集中质量矩阵的形式

在集中质量法中，如果忽略转动惯量，质量矩阵通常是对角矩阵：

$$M=\text{diag}(m_1,m_2,\dots ,m_n)$$

这种对角形式大大简化了计算，因为：

1. 矩阵求逆简单（只需对对角元素求倒数）
2. 矩阵乘法计算量小
3. 存储需求低
4. 便于并行计算

如果考虑转动自由度，质量矩阵可能包含转动惯量项：

$$M=\text{diag}(m_1,m_1,J_1,m_2,m_2,J_2,\dots )$$

其中$m_i$是平移质量，$J_i$是转动惯量。

4.4 集中质量法的优缺点

优点：

• 概念简单，易于理解和实现
• 计算效率高，特别适合大规模问题
• 质量矩阵为对角阵，便于数值计算
• 物理意义明确，便于工程应用

缺点：

• 精度相对较低，特别是对于高阶模态
• 转动惯量的处理较为粗糙
• 对于质量分布不均匀的结构，近似误差较大
• 在某些情况下可能产生虚假的数值振荡

4.5 应用案例：三层剪切框架

考虑一个三层剪切型框架结构，每层质量为$m_i$，层间刚度为$k_i$。采用集中质量法，将每层质量集中到楼层中心，得到三自由度系统。

质量矩阵为：

$$M=\left[\begin{array}{ccc}{m}_1& 0& 0\\ 0& m_2& 0\\ 0& 0& m_3\end{array}\right]$$

假设$m_1=1000$ kg，$m_2=800$ kg，$m_3=600$ kg，则质量矩阵为：

$$M=\left[\begin{array}{ccc}1000& 0& 0\\ 0& 800& 0\\ 0& 0& 600\end{array}\right]\text{ kg}$$

这种对角形式的质量矩阵使得运动方程的求解更加高效。

5. 一致质量法

5.1 一致质量法的理论基础

一致质量法（Consistent Mass Method）基于有限元法的思想，通过形函数将分布质量离散化。这种方法保持了质量矩阵与刚度矩阵的一致性，通常能提供更精确的结果。

一致质量矩阵的定义为：

$$M={\int }_V\rho N^{T}NdV$$

其中：

• $\rho$ 是材料密度
• $N$ 是形函数矩阵
• $V$ 是单元体积

这个积分表示将分布质量通过形函数加权分配到各个节点。

5.2 形函数的概念

形函数（Shape Function）是有限元方法的核心概念，它描述了单元内部位移场与节点位移之间的关系。对于一维杆单元，常用的形函数为线性形函数：

$$N_1(x)=1-\frac{x}{L},N_2(x)=\frac{x}{L}$$

其中$L$是单元长度。单元内任意点的位移可以表示为：

$$u(x)=N_1(x)u_1+N_2(x)u_2$$

对于梁单元，需要考虑横向位移和转角，通常采用三次Hermite形函数：

$$N_1=1-3{\xi }^2+2{\xi }^3$$
$$N_2=L(\xi -2{\xi }^2+{\xi }^3)$$
$$N_3=3{\xi }^2-2{\xi }^3$$
$$N_4=L(-{\xi }^2+{\xi }^3)$$

其中$\xi =x/L$是无量纲坐标。

5.3 一致质量矩阵的推导

以均匀杆单元为例，推导一致质量矩阵。设杆长为$L$，截面积为$A$，密度为$\rho$。

形函数矩阵为：

$$N=[N_1{N}_2]=\left[1-\frac{x}{L}\frac{x}{L}\right]$$

一致质量矩阵为：

$$M={\int }_0^L\rho A{N}^{T}Ndx$$

计算各元素：

$$m_{11}={\int }_0^L\rho A{\left(1-\frac{x}{L}\right)}^{2}dx=\frac{\rho AL}{3}$$

$$m_{12}=m_{21}={\int }_0^L\rho A\left(1-\frac{x}{L}\right)\frac{x}{L}dx=\frac{\rho AL}{6}$$

$$m_{22}={\int }_0^L\rho A{\left(\frac{x}{L}\right)}^{2}dx=\frac{\rho AL}{3}$$

因此，杆单元的一致质量矩阵为：

$$M=\frac{\rho AL}{6}\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right]$$

注意到总质量$m=\rho AL$，所以：

$$M=\frac{m}{6}\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right]$$

这与集中质量法（$M=\text{diag}(m/2,m/2)$）有明显不同，一致质量矩阵存在非对角耦合项。

5.4 梁单元的一致质量矩阵

对于Euler-Bernoulli梁单元，考虑横向振动，每个节点有横向位移和转角两个自由度。一致质量矩阵为4×4矩阵：

$$M=\frac{\rho AL}{420}\left[\begin{array}{cccc}156& 22L& 54& -13L\\ 22L& 4L^2& 13L& -3L^2\\ 54& 13L& 156& -22L\\ -13L& -3L^2& -22L& 4L^2\end{array}\right]$$

这个矩阵体现了平移质量和转动惯量的耦合。与集中质量法相比，一致质量法更准确地描述了梁的惯性特性。

5.5 一致质量法的优缺点

优点：

• 精度高，特别是高阶模态的计算
• 与有限元刚度矩阵一致，保持能量协调
• 能更好地描述转动惯量效应
• 收敛性好，随着网格加密快速收敛到精确解

缺点：

• 质量矩阵为满阵或带状矩阵，计算量大
• 需要数值积分，实现复杂
• 存储需求高
• 对于大规模问题，计算效率较低

5.6 质量矩阵的比较

特性	集中质量法	一致质量法
矩阵形式	对角阵	满阵或带状阵
计算效率	高	较低
精度	较低	较高
存储需求	低	较高
实现难度	简单	复杂
收敛性	一般	好

在实际应用中，集中质量法常用于初步分析和大型问题，而一致质量法用于需要高精度的场合。

6. 多自由度系统建模实例

6.1 两自由度弹簧-质量系统

考虑一个简单的两自由度弹簧-质量系统，如图5.1所示。两个质量$m_1$和$m_2$通过弹簧$k_1$、$k_2$、$k_3$连接。

系统参数：

• $m_1=2$ kg，$m_2=1$ kg
• $k_1=100$ N/m，$k_2=150$ N/m，$k_3=100$ N/m

质量矩阵：

采用集中质量法，质量矩阵为对角阵：

$$M=\left[\begin{array}{cc}2& 0\\ 0& 1\end{array}\right]\text{ kg}$$

刚度矩阵：

根据力的平衡，刚度矩阵为：

$$K=\left[\begin{array}{cc}{k}_1+k_2& -k_2\\ -k_2& k_2+k_3\end{array}\right]=\left[\begin{array}{cc}250& -150\\ -150& 250\end{array}\right]\text{ N/m}$$

运动方程：

$$\left[\begin{array}{cc}2& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}{\ddot{x}}_1\\ {\ddot{x}}_2\end{array}\right]+\left[\begin{array}{cc}250& -150\\ -150& 250\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}{F}_1(t)\\ F_2(t)\end{array}\right]$$

6.2 简支梁的离散化

考虑一个长度为$L$、抗弯刚度为$EI$、单位长度质量为$\rho A$的简支梁。将其离散为n个单元。

集中质量法离散化：

将梁等分为n段，每段质量$m=\rho AL/n$集中在单元节点上。对于简支梁，内部节点质量为相邻两段质量之和，端部节点质量为半段质量。

对于n=4的离散化：

• 节点1（左端）：$m_1=\rho AL/8$
• 节点2：$m_2=\rho AL/4$
• 节点3：$m_3=\rho AL/4$
• 节点4：$m_4=\rho AL/4$
• 节点5（右端）：$m_5=\rho AL/8$

一致质量法离散化：

使用梁单元的一致质量矩阵组装整体质量矩阵。对于n个单元，有n+1个节点，每个节点2个自由度（横向位移和转角），总自由度为2(n+1)。

6.3 三层框架结构

考虑一个三层钢筋混凝土框架结构，如图5.2所示。每层高度$h=3$ m，柱截面尺寸$0.4\times 0.4$ m，混凝土弹性模量$E=30$ GPa，密度$\rho =2500$ kg/m³。

参数计算：

柱截面惯性矩：
$$I=\frac{0.4^4}{12}=2.133\times 10^{-3}{\text{ m}}^4$$

柱的线刚度：
$$i_c=\frac{EI}{h}=\frac{30\times 10^9\times 2.133\times 10^{-3}}{3}=2.133\times 10^7\text{ N⋅m}$$

假设每层有4根柱，层间侧移刚度近似为：
$$k=\frac{12\times 4\times i_c}{h^2}=\frac{48\times 2.133\times 10^7}{9}=1.138\times 10^8\text{ N/m}$$

质量计算：

假设每层面积$200$ m²，楼板厚度$0.15$ m，则每层质量：
$$m=\rho \times A\times t=2500\times 200\times 0.15=75000\text{ kg}$$

考虑活荷载和隔墙重量，取$m=100000$ kg/层。

质量矩阵：

$$M=10^5\times \left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]\text{ kg}$$

7. 质量矩阵的组装方法

7.1 有限元组装原理

对于复杂结构，需要将各个单元的质量矩阵组装成整体质量矩阵。组装过程遵循以下原则：

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

7.2 坐标变换

当单元方向与整体坐标系不一致时，需要进行坐标变换。变换公式为：

$$M^{(g)}=T^T{M}^{(l)}T$$

其中：

• $M^{(l)}$ 是局部坐标系下的单元质量矩阵
• $M^{(g)}$ 是整体坐标系下的单元质量矩阵
• $T$ 是坐标变换矩阵

对于平面杆单元，与x轴夹角为$\theta$时，变换矩阵为：

$$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]$$

7.3 组装算法

整体质量矩阵的组装算法如下：

初始化整体质量矩阵 M_global = 0
对于每个单元 e:
    计算单元质量矩阵 M_e
    获取单元节点编号 [i, j, ...]
    对于 M_e 中的每个元素 M_e[p][q]:
        确定在整体矩阵中的位置 [I, J]
        M_global[I][J] += M_e[p][q]

7.4 边界条件的处理

在组装完成后，需要施加边界条件。常用的方法有：

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

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

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

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

8. 数值实现与Python编程

8.1 质量矩阵的数据结构

在Python中，质量矩阵可以用多种数据结构表示：

密集矩阵：使用NumPy的ndarray，适合小规模问题

稀疏矩阵：使用scipy.sparse，适合大规模问题，节省内存

对角矩阵：对于集中质量法，可以只存储对角元素

8.2 集中质量矩阵的实现

import numpy as np

def lumped_mass_matrix(masses):
    """
    构建集中质量矩阵
    
    参数:
        masses: 各节点质量列表或数组
    
    返回:
        M: 对角质量矩阵
    """
    n = len(masses)
    M = np.diag(masses)
    return M

# 示例：三自由度系统
masses = [1000, 800, 600]  # kg
M = lumped_mass_matrix(masses)
print("集中质量矩阵:")
print(M)

8.3 一致质量矩阵的实现

def consistent_mass_matrix_bar(rho, A, L):
    """
    计算杆单元的一致质量矩阵
    
    参数:
        rho: 材料密度 (kg/m³)
        A: 截面积 (m²)
        L: 杆长 (m)
    
    返回:
        M: 2×2一致质量矩阵
    """
    m = rho * A * L  # 单元总质量
    M = (m / 6) * np.array([[2, 1],
                            [1, 2]])
    return M

def consistent_mass_matrix_beam(rho, A, L):
    """
    计算Euler-Bernoulli梁单元的一致质量矩阵
    
    参数:
        rho: 材料密度 (kg/m³)
        A: 截面积 (m²)
        L: 梁长 (m)
    
    返回:
        M: 4×4一致质量矩阵
    """
    m = rho * A * L  # 单元总质量
    M = (m / 420) * np.array([[156, 22*L, 54, -13*L],
                               [22*L, 4*L**2, 13*L, -3*L**2],
                               [54, 13*L, 156, -22*L],
                               [-13*L, -3*L**2, -22*L, 4*L**2]])
    return M

8.4 质量矩阵组装函数

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

9. 案例研究：多自由度系统自由振动分析

9.1 问题描述

考虑一个四自由度弹簧-质量系统，如图5.3所示。各质量块通过弹簧连接，分析其自由振动特性。

系统参数：

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

9.2 质量矩阵构建

采用集中质量法：

$$M=\left[\begin{array}{cccc}2& 0& 0& 0\\ 0& 1.5& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 0.5\end{array}\right]\text{ kg}$$

9.3 刚度矩阵构建

根据力的平衡：

$$K=\left[\begin{array}{cccc}{k}_1+k_2& -k_2& 0& 0\\ -k_2& k_2+k_3& -k_3& 0\\ 0& -k_3& k_3+k_4& -k_4\\ 0& 0& -k_4& k_4+k_5\end{array}\right]$$

$$K=\left[\begin{array}{cccc}350& -150& 0& 0\\ -150& 250& -100& 0\\ 0& -100& 180& -80\\ 0& 0& -80& 130\end{array}\right]\text{ N/m}$$

9.4 特征值问题求解

自由振动问题转化为广义特征值问题：

$$(K-{\omega }^{2}M)\varphi =0$$

通过求解特征值问题，得到系统的固有频率和振型。

9.5 Python实现

import numpy as np
from scipy.linalg import eig

# 系统参数
m = np.array([2.0, 1.5, 1.0, 0.5])  # kg
k = np.array([200, 150, 100, 80, 50])  # N/m

# 构建质量矩阵（集中质量法）
M = np.diag(m)

# 构建刚度矩阵
K = np.array([[k[0]+k[1], -k[1], 0, 0],
              [-k[1], k[1]+k[2], -k[2], 0],
              [0, -k[2], k[2]+k[3], -k[3]],
              [0, 0, -k[3], k[3]+k[4]]])

# 求解特征值问题
eigenvalues, eigenvectors = eig(K, M)

# 排序
idx = np.argsort(eigenvalues)
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

# 计算固有频率
omega = np.sqrt(eigenvalues.real)
f = omega / (2 * np.pi)

print("固有频率 (Hz):")
for i in range(len(f)):
    print(f"  模态 {i+1}: {f[i]:.4f} Hz")

print("\n振型:")
for i in range(len(f)):
    print(f"  模态 {i+1}: {eigenvectors[:, i]}")

10. 高级主题

10.1 质量矩阵的缩聚技术

对于大型结构，自由度数目可能非常庞大。质量矩阵缩聚技术通过静力或动力缩聚，减少自由度数目，提高计算效率。

静力缩聚（Guyan缩聚）：

将自由度分为主自由度（保留）和从自由度（缩聚）。通过静力关系，从自由度可以用主自由度表示：

$$u_s=-K_{ss}^{-1}{K}_{sm}{u}_m$$

缩聚后的质量矩阵：

$$M_r=T^{T}MT$$

其中$T$是变换矩阵。

动力缩聚：

考虑惯性效应，通过迭代方法获得更精确的缩聚矩阵。

10.2 转动惯量的处理

对于需要考虑转动效应的结构，质量矩阵应包含转动惯量项。

对于刚体，转动惯量矩阵为：

$$J=\left[\begin{array}{ccc}{J}_{xx}& J_{xy}& J_{xz}\\ J_{yx}& J_{yy}& J_{yz}\\ J_{zx}& J_{zy}& J_{zz}\end{array}\right]$$

其中对角元素是绕各轴的转动惯量，非对角元素是惯性积。

10.3 非结构质量的处理

实际工程中，结构往往承受附加质量，如：

• 设备质量
• 液体质量（考虑晃动效应）
• 覆冰质量
• 人群荷载

这些非结构质量需要以适当的方式添加到质量矩阵中。

10.4 变质量系统

某些结构的质会随时间变化，如：

• 火箭发射（燃料消耗）
• 建筑拆除
• 积雪融化

变质量系统的运动方程需要考虑质量变化率的影响。

11. 工程应用实例

11.1 高层建筑抗震分析

某30层钢筋混凝土框架-剪力墙结构，层高3m，总高90m。采用集中质量法建模，每层质量约1500吨。

建模要点：

1. 每层简化为一个质点，具有水平平移自由度
2. 考虑偶然偏心，每层增加扭转自由度
3. 考虑P-Δ效应的几何刚度
4. 采用Rayleigh阻尼模型

质量矩阵：

$$M=\text{diag}(m_1,m_1,J_1,m_2,m_2,J_2,\dots ,m_{30},m_{30},J_{30})$$

其中$m_i$是第i层质量，$J_i$是第i层转动惯量。

11.2 桥梁结构动力分析

某三跨连续梁桥，跨径布置为60m+100m+60m。采用一致质量法建模，每个节点考虑竖向位移和转角。

建模要点：

1. 主梁离散为梁单元
2. 考虑桥墩的柔性
3. 考虑支座刚度
4. 考虑车辆质量（移动荷载）

11.3 旋转机械转子动力学

汽轮机转子系统，具有多个叶轮和轴承支撑。采用集中质量法，将叶轮简化为刚性圆盘。

建模要点：

1. 轴段采用梁单元
2. 叶轮简化为刚性圆盘，具有质量和转动惯量
3. 轴承采用弹簧-阻尼元件
4. 考虑陀螺效应

12. 误差分析与精度评估

12.1 离散化误差

离散化误差来源于将连续系统近似为离散系统。主要影响因素：

单元尺寸：单元越小，误差越小，但计算成本增加

形函数阶次：高阶形函数可以提高精度

质量分布近似：集中质量法的近似误差

12.2 收敛性分析

通过网格加密研究，评估解的收敛性。对于一致质量法，固有频率的收敛速度通常为$h^{2p}$，其中$h$是单元尺寸，$p$是形函数多项式阶次。

12.3 精度改进方法

1. 网格加密：增加单元数目
2. 高阶单元：采用二次或三次形函数
3. 混合质量法：结合集中质量和一致质量的优点
4. 自适应网格：根据误差估计自动调整网格

13. 总结与展望

13.1 主要内容回顾

本主题系统介绍了多自由度系统的建模方法和质量矩阵构建技术：

1. 多自由度系统基本概念：自由度的定义、分类和选择原则
2. 质量矩阵的物理意义：定义、性质和质量耦合概念
3. 集中质量法：原理、实现和优缺点
4. 一致质量法：基于有限元的质量离散化方法
5. 质量矩阵组装：坐标变换和组装算法
6. 数值实现：Python编程实现质量矩阵构建
7. 工程应用：高层建筑、桥梁、旋转机械等实例

13.2 关键知识点

• 质量矩阵是对称正定矩阵，描述系统的惯性特性
• 集中质量法简单高效，质量矩阵为对角阵
• 一致质量法精度高，但计算成本较高
• 质量矩阵的形式依赖于坐标系的选择
• 质量矩阵与刚度矩阵协调使用，保证能量守恒

13.3 后续学习方向

掌握了质量矩阵构建方法后，下一步将学习：

1. 刚度矩阵构建：建立系统的弹性恢复力模型
2. 阻尼矩阵构建：描述能量耗散机制
3. 模态分析：求解固有频率和振型
4. 时程响应分析：计算系统在外荷载作用下的动力响应
5. 频域分析：频率响应函数和传递函数

13.4 工程实践建议

1. 模型简化：在保证精度的前提下，合理简化模型
2. 参数识别：通过试验数据识别质量参数
3. 模型验证：与理论解或试验结果对比验证
4. 灵敏度分析：研究参数变化对结果的影响
5. 不确定性分析：考虑参数不确定性的影响

14. 练习题

14.1 理论题

1. 解释质量矩阵的物理意义，为什么质量矩阵是对称正定的？

2. 比较集中质量法和一致质量法的优缺点，说明各自的适用场合。

3. 什么是质量耦合？在什么情况下会出现质量耦合？

4. 推导杆单元的一致质量矩阵，并与集中质量矩阵比较。

5. 说明质量矩阵缩聚的基本原理和应用场合。

14.2 计算题

1. 一个三自由度系统，质量分别为$m_1=3$ kg, $m_2=2$ kg, $m_3=1$ kg。构建其集中质量矩阵和一致质量矩阵（假设为均匀杆离散化）。

2. 一个悬臂梁长$L=2$ m，截面为矩形$0.1\times 0.2$ m，材料密度$\rho =7850$ kg/m³。将其离散为2个单元，计算一致质量矩阵。

3. 对于第2题的悬臂梁，分别用集中质量法和一致质量法计算前两阶固有频率，并比较结果。

14.3 编程题

1. 编写Python函数，实现任意自由度集中质量系统的质量矩阵构建。

2. 编写Python程序，实现梁单元一致质量矩阵的自动组装。

3. 编写程序比较集中质量法和一致质量法在计算悬臂梁固有频率时的精度差异。

15. 参考文献

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
import matplotlib.patches as mpatches

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

def lumped_mass_matrix(masses):
    """构建集中质量矩阵"""
    return np.diag(masses)

def consistent_mass_matrix_bar(rho, A, L):
    """杆单元一致质量矩阵"""
    m = rho * A * L
    return (m / 6) * np.array([[2, 1], [1, 2]])

def consistent_mass_matrix_beam(rho, A, L):
    """梁单元一致质量矩阵"""
    m = rho * A * L
    return (m / 420) * np.array([
        [156, 22*L, 54, -13*L],
        [22*L, 4*L**2, 13*L, -3*L**2],
        [54, 13*L, 156, -22*L],
        [-13*L, -3*L**2, -22*L, 4*L**2]
    ])

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

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

m1, m2 = 2.0, 1.0  # kg
k1, k2, k3 = 100, 150, 100  # N/m

# 集中质量矩阵
M_lumped = lumped_mass_matrix([m1, m2])
print("\n集中质量矩阵 M:")
print(M_lumped)

# 刚度矩阵
K = np.array([[k1+k2, -k2], [-k2, k2+k3]])
print("\n刚度矩阵 K:")
print(K)

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

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

print("\n固有频率:")
print(f"  ω₁ = {omega[0]:.4f} rad/s, f₁ = {f[0]:.4f} Hz")
print(f"  ω₂ = {omega[1]:.4f} rad/s, f₂ = {f[1]:.4f} Hz")

print("\n振型:")
for i in range(2):
    mode = eigenvectors[:, i]
    mode = mode / np.max(np.abs(mode))  # 归一化
    print(f"  模态 {i+1}: [{mode[0]:.4f}, {mode[1]:.4f}]ᵀ")

# 案例2：四自由度系统
print("\n" + "=" * 60)
print("案例2：四自由度弹簧-质量系统")
print("=" * 60)

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

M = lumped_mass_matrix(m)
K = build_stiffness_matrix_spring(k, 4)

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

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(4):
    print(f"  ω{i+1} = {omega[i]:.4f} rad/s, f{i+1} = {f[i]:.4f} Hz")

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

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

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

# 图3：固有频率
ax3 = axes[1, 0]
mode_numbers = np.arange(1, 5)
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]
x_pos = np.arange(4)
width = 0.2
colors = ['red', 'green', 'blue', 'orange']
for i in range(4):
    mode = eigenvectors[:, i]
    mode = mode / np.max(np.abs(mode))
    ax4.bar(x_pos + i*width, mode, width, label=f'模态{i+1}', color=colors[i], alpha=0.7)
ax4.set_xlabel('质量块编号')
ax4.set_ylabel('归一化位移')
ax4.set_title('系统振型', fontsize=12)
ax4.set_xticks(x_pos + width * 1.5)
ax4.set_xticklabels(['m₁', 'm₂', 'm₃', 'm₄'])
ax4.legend()
ax4.axhline(y=0, color='k', linestyle='-', linewidth=0.5)
ax4.grid(True, alpha=0.3)

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

# 案例3：集中质量法 vs 一致质量法
print("\n" + "=" * 60)
print("案例3：集中质量法 vs 一致质量法")
print("=" * 60)

# 均匀杆参数
rho = 7850  # kg/m³ (钢)
A = 0.01    # m²
L_total = 10  # m
E = 2e11    # Pa

# 离散化为不同数量的单元
n_elements_list = [1, 2, 4, 8, 16]
print("\n悬臂梁固有频率比较 (Hz)")
print("-" * 60)
print(f"{'单元数':<10} {'集中质量法':<20} {'一致质量法':<20} {'理论解':<15}")
print("-" * 60)

# 理论解（悬臂梁第一阶固有频率）
omega_theory = 3.5156 * np.sqrt(E * 0.0001 / (rho * A * L_total**4))
f_theory = omega_theory / (2 * np.pi)

for n_elem in n_elements_list:
    L = L_total / n_elem
    m_elem = rho * A * L
    
    # 集中质量法
    M_lumped_diag = np.full(n_elem, m_elem)
    M_lumped_diag[0] = m_elem / 2  # 端部质量减半
    M_lumped_full = np.diag(M_lumped_diag)
    
    # 一致质量法
    M_consistent = np.zeros((n_elem, n_elem))
    M_e = consistent_mass_matrix_bar(rho, A, L)
    for i in range(n_elem):
        M_consistent[i, i] += M_e[0, 0] if i == 0 else M_e[0, 0] + M_e[1, 1]
        if i < n_elem - 1:
            M_consistent[i, i+1] = M_e[0, 1]
            M_consistent[i+1, i] = M_e[1, 0]
    
    # 简化刚度矩阵（假设）
    k_elem = E * A / L
    K_simple = np.zeros((n_elem, n_elem))
    for i in range(n_elem):
        K_simple[i, i] += k_elem if i == 0 else 2 * k_elem
        if i < n_elem - 1:
            K_simple[i, i+1] = -k_elem
            K_simple[i+1, i] = -k_elem
    K_simple[n_elem-1, n_elem-1] = k_elem
    
    # 求解
    try:
        _, eigvec_lumped = eig(K_simple, M_lumped_full)
        f_lumped = np.sqrt(np.sort(np.real(eigvec_lumped))[0]) / (2 * np.pi)
    except:
        f_lumped = 0
    
    try:
        _, eigvec_cons = eig(K_simple, M_consistent)
        f_consistent = np.sqrt(np.sort(np.real(eigvec_cons))[0]) / (2 * np.pi)
    except:
        f_consistent = 0
    
    print(f"{n_elem:<10} {f_lumped:<20.6f} {f_consistent:<20.6f} {f_theory:<15.6f}")

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

# 初始条件
u0 = np.array([1.0, 0.5, 0.0, 0.0])  # 初始位移
v0 = np.zeros(4)  # 初始速度

# 模态叠加法求解响应
dt = 0.01
t_max = 10
t = np.arange(0, t_max, dt)
n_steps = len(t)

# 计算模态坐标初始条件
n_modes = 4
q0 = np.zeros(n_modes)
dq0 = np.zeros(n_modes)

for i in range(n_modes):
    phi_i = eigenvectors[:, i]
    # 质量归一化
    phi_i = phi_i / np.sqrt(np.dot(phi_i, np.dot(M, phi_i)))
    q0[i] = np.dot(phi_i, np.dot(M, u0))
    dq0[i] = np.dot(phi_i, np.dot(M, v0))

# 计算响应
u_history = np.zeros((n_steps, 4))
for i in range(n_steps):
    u = np.zeros(4)
    for j in range(n_modes):
        phi_j = eigenvectors[:, j]
        phi_j = phi_j / np.sqrt(np.dot(phi_j, np.dot(M, phi_j)))
        omega_j = omega[j]
        # 自由振动响应
        q_j = q0[j] * np.cos(omega_j * t[i]) + (dq0[j] / omega_j) * np.sin(omega_j * t[i])
        u += q_j * phi_j
    u_history[i, :] = u

# 创建动画
fig_anim, ax_anim = plt.subplots(figsize=(12, 6))
ax_anim.set_xlim(-0.5, 4.5)
ax_anim.set_ylim(-1.5, 1.5)
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(5):
    line, = ax_anim.plot([], [], 'k-', linewidth=2)
    spring_lines.append(line)

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

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

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_history[frame, :]
    x_positions = np.arange(4) + 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(3):
        spring_lines[i+1].set_data([x_positions[i], x_positions[i+1]], [0, 0])
    # 右端弹簧（假设连接到固定点）
    spring_lines[4].set_data([x_positions[3], 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('mdof_vibration_animation.gif', writer='pillow', fps=20)
print("动画已保存: mdof_vibration_animation.gif")

plt.close('all')

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