Newmark-beta法结构动力响应计算程序

资 源 简 介

本项目开发一套基于Newmark-beta数值积分算法的结构动力方程求解器，旨在精确计算工程结构在地震、风载或机械振动等动态荷载作用下的时程响应。该程序的核心功能详细包括：1. 动力方程求解：针对单自由度（SDOF）及多自由度（MDOF）线性弹性体系，建立通用的二阶常微分动力平衡方程数值解法。2. 算法灵活配置：完整实现Newmark-beta逐步积分公式，允许用户通过调整参数（gamma和beta）来选择不同的积分方案，如无条件稳定的“平均加速度法”（gamma=0.5, beta=0.25）或条件稳定的“线性加速度法”（gamma=0.5, beta=1/6），以适应不同刚度特性的结构分析需求。3. 自动化计算流程：程序初始化阶段自动计算积分常数并组装有效刚度矩阵，在时间步进循环中递归计算每一时刻的有效荷载向量，进而求解当前步的位移，并根据运动学关系更新速度和加速度。4. 阻尼模型处理：内置瑞利阻尼（Rayleigh Damping）计算模块，可根据用户提供的模态阻尼比和参考频率自动构造系统的比例阻尼矩阵。5. 数据分析与可视化：能够完整记录结构关键节点在全时域内的动力响应数据，并绘制位移、速度、加速度的时程曲线，自动捕捉并输出最大响应值，为结构抗震性能评估和设计优化提供直接依据。

详 情 说 明

结构动力反应分析 Newmark-beta 法计算程序

项目简介

本项目是一个基于 MATLAB 开发的结构动力学数值仿真工具，专门用于求解多自由度（MDOF）线性弹性体系在动态荷载作用下的时程响应。核心算法采用 Newmark-beta 逐步积分法，能够精确计算结构在地震、机械振动或简谐激励下的位移、速度和加速度响应。

程序内置了一个典型的三层剪切型框架结构模型，演示了从参数定义、阻尼矩阵构造、动力方程求解到结果可视化的完整流程。

功能特性

• 多自由度体系求解：支持任意自由度数量的线性系统动力分析（示例代码配置为 3 自由度剪切模型）。
• Newmark-beta 数值积分：实现了通用的 Newmark 数值积分算法，默认配置为“平均加速度法”（$gamma=0.5, beta=0.25$），具有无条件稳定性，也可调整参数实现线性加速度法。
• Rayleigh 阻尼自动计算：内置瑞利阻尼（Rayleigh Damping）生成模块，根据用户设定的目标阻尼比（如 5%）和系统的前两阶固有频率，自动求解质量比例系数 $alpha$ 和刚度比例系数 $beta$，构造正交阻尼矩阵。
• 高效矩阵求解：在积分过程中采用了有效刚度矩阵（Effective Stiffness Matrix）的 LU 分解技术，避免了在每个时间步重复求逆，显著提高了计算效率。
• 特定的动态激励模拟：预设了作用于结构顶层的正弦动态荷载，频率设定为避开共振区的特定值（基频的 0.8 倍），用于观察系统的受迫振动响应。
• 全面的后处理与可视化：

* 自动统计并输出每一层的最大位移、最大速度和最大加速度。 * 绘制所有楼层的位移、速度、加速度时程曲线。 * 绘制顶层节点的相平面图（位移-速度关系）。

系统要求

• MATLAB R2016a 或更高版本（需包含基础数学计算功能）。
• 无需额外的工具箱（Toolbox）。

使用方法

1. 将 main.m 文件保存在 MATLAB 的工作路径中。
2. 直接运行 main 函数。
3. 程序将在命令行窗口输出结构的动力特性（频率、阻尼系数）和响应极值统计。
4. 程序将自动弹出两个图形窗口，分别显示时程响应曲线和相平面图。

代码实现逻辑与算法详解

本项目完全基于 main.m 文件中的逻辑实现，主要包含以下几个核心模块：

1. 模型参数定义与组装

程序首先初始化一个三层剪切型框架结构模型。

• 质量矩阵 (M)：构建对角矩阵，假设每层质量均为 4000 kg。
• 刚度矩阵 (K)：基于剪切变形假设，通过层间刚度（$k_{story} = 2 times k_{col}$）组装三对角刚度矩阵。刚度矩阵反映了层间耦合关系，自对角项为相邻层刚度之和，非对角项为负的层间刚度。
• 积分参数：设定 Newmark 方法的积分参数 $gamma=0.5$ 和 $beta=0.25$，对应无条件稳定的平均加速度法；时间步长设定为 0.02 秒。

2. 瑞利阻尼 (Rayleigh Damping) 构造

为了模拟真实的能量耗散，程序实现了比例阻尼矩阵的自动化计算：

• 特征值分析：调用 eig(K, M) 求解广义特征值问题，获取系统的固有频率（$omega$）。
• 系数求解：选取前两阶模态频率（$omega_1, omega_2$）和目标阻尼比（$zeta=0.05$），求解方程组 $0.5(alpha/omega_i + betaomega_i) = zeta$，得到 Rayleight 阻尼系数 $alpha_R$ 和 $beta_R$。
• 矩阵生成：利用公式 $C = alpha_R M + beta_R K$ 组装最终的阻尼矩阵。

3. 荷载定义

程序模拟了一个单点激励工况：

• 在结构的最高层（第 3 层）施加正弦荷载 $P(t) = 10000 sin(0.8 omega_1 t)$。
• 其他楼层不受外力作用。

4. Newmark-beta 求解器 (核心算法)

核心求解逻辑封装在 NewmarkSolver 子函数中，执行步骤如下：

• 初始化常数：预先计算积分所需的常数 $a_0$ 至 $a_7$，这些常数仅依赖于时间步长 $dt$ 和积分参数 $beta, gamma$。
• 初始状态计算：根据初始位移 $u_0$、初始速度 $v_0$ 和初始荷载 $P_0$，通过动力平衡方程求解初始加速度 $a_0 = M^{-1}(P_0 - C v_0 - K u_0)$。
• 有效刚度矩阵：构建有效刚度矩阵 $K_{eff} = K + a_0 M + a_1 C$。
• 矩阵分解：利用 lu 函数对 $K_{eff}$ 进行 LU 分解，得到 $L_{eff}$ 和 $U_{eff}$，以便在后续时间步中快速回代求解。
• 时间步进循环：

1. 根据当前时刻的 $u, v, a$ 计算有效荷载向量 $P_{eff}$。 2. 利用分解后的刚度矩阵求解下一时刻的位移 $u_{next}$。 3. 根据 Newmark 运动学公式更新下一时刻的加速度 $a_{next}$ 和速度 $v_{next}$。 4. 保存结果并进入下一时间步。

5. 数据分析与可视化

计算完成后，程序通过两个子函数处理结果：

• visualizeResults：创建一个包含三个子图的窗口，分别绘制各层的位移、速度和加速度随时间的变化曲线，并在同一坐标系中通过不同颜色区分楼层。额外绘制一个窗口显示顶层的相平面轨迹（位移 vs 速度）。
• calculateAndPrintStats：遍历计算结果矩阵，提取每一自由度的绝对最大值，并格式化输出到控制台，便于用户快速评估结构的最大响应。

输出示例说明

运行程序后，您将看到如下形式的文本输出：

• 结构动力特性：显示基频（f1）、二阶频率（f2）以及计算得到的瑞利阻尼系数 alpha 和 beta。
• 计算耗时：显示求解过程所消耗的 CPU 时间。
• 响应极值统计：表格形式列出每一层（Floor 1 至 Floor 3）在全时程范围内的最大位移（m）、最大速度（m/s）和最大加速度（m/s^2）。