基于有限单元的动力响应分析法在MATLAB 中的应用研究

1 引言

采用有限单元对结构进行振动分析的方法可分为动力响应分析与振动特性分析2 种。动力响应分析主要求解的是结构在外力作用下发生振动时，结构的位移、速度与加速度随时间的变化关系；振动特性分析主要求解的是结构的固有属性，其中包括固有频率、模态振型等。

基于有限单元的动力响应分析法首先通过d’Alembert 原理建立单元体的运动方程，然后按照有限元的集合方法，将单元体的运动微分方程集合成结构整体的运动方程，最后通过逐步积分的方法求解出结构的位移、速度与加速度随时间的变化关系，从而能够研究出结构在动力荷载作用下的内力变化以及破坏过程。

其中，单元体的运动方程中包含有单元体的刚度矩阵、质量矩阵和阻尼矩阵。在将单元体运动方程集合成结构整体运动方程的过程中需要将单元体刚度矩阵、质量矩阵和阻尼矩阵集合成相应的整体矩阵。当结构需要得到更高的计算精度时，可以通过增加结构的单元划分数量来实现，但是，势必也会增加矩阵之间的运算量。

而MATLAB 作为一款优异的数值计算软件，其基础数据单位为矩阵，能够将工程问题与数值计算之间实现更好的连接。结构的刚度矩阵与质量矩阵中往往含有大量的零元素，而MATLAB 中内置的稀疏矩阵正好可以解决上述矩阵的运算量问题，在节约计算机储存空间的同时加快矩阵的运算速度，对于处理此类结构振动分析问题具有优异的性能[1,2]。

2 基于有限单元的动力响应分析法

2.1 简介

有限单元法通过将结构离散成一组指定数量的单元体，并在每个离散的单元体的指定点处设置结点使得相邻单元体的有关参数拥有一定的连续性，然后，用每个单元上假设的函数来近似表示整个结构上待求的未知场函数[3]。

在动力响应分析所处理的结构动力学问题中，未知场函数一般为结构自由振动或者受迫振动产生的位移、速度与加速度。结构整体的位移、速度与加速度时程数据可以从每个单元体结点处的响应数据近似求得。理论上结构划分单元体的数量越多，最后的结果便越精确。此时，结构的动力响应分析便从一个连续体无限自由度问题转换为多自由度问题，自由度的数量即为单元体数量。

2.2 Newmark 法

Newmark 法是一种常用的动力响应分析法。Newmark 法通过假设时间步内的加速度分布，依据积分的方法来获得时间步内结构的速度与位移表达式，从而得到时间步结束点处的速度与位移。

多自由度系统在外力荷载作用下的整体运动方程如下[4]：

式中，M、C 和K 分别为多自由度系统的整体质量矩阵、阻尼矩阵和刚度矩阵；V（··t）、V（·t）、V（t）分别为多自由度系统的单元结点加速度、速度和位移；F（t）为多自由度系统的外力矩阵。

假设该系统具有初始位移以及初始速度如下：

将公式（2）代入公式（1）中可立刻求得多自由度系统初始加速度：

式中，M-1为整体质量矩阵的逆矩阵。

Newmark 法假设的加速度分布模式有如下公式[5]：

当式（5）中β=1/4、γ=1/2 时，Newmark 法假设的加速度分布模式为常加速度分布；当式（5）中的β=1/6、γ=1、2 时，Newmark 法假设的加速度分布模式为线加速度分布。以下有关Newmark 内容的讨论均按照常加速度分布进行。

当β=1/4、γ=1/2 时，式（4）、（5）的解有如下形式：

在这里定义7 个常数：

代入系统的初始位移、初始速度和初始加速度后，（6）式中的等效荷载F~与等效刚度K~可以由下式（8）、（9）分别求得：

在依据式（6）、（7）、（8）求解出多自由度系统任意时间步的位移V（t+Δt）后，系统任意时刻的速度即可由式（10）求得：

上式中的多自由度系统任意时刻的加速度可由下式求得：

基于有限单元的Newmark 法的具体过程为：（1）连续体结构离散化。划分单元体的分割方案需要综合考虑结构受力情况、约束情况以及求解结果的精度要求，划分单元体的形状应根据具体工程问题选择为杆单元、梁单元（包含平面梁单元以及空间梁单元）、三角单元以及矩形单元等中的一种或多种。（2）根据具体选择的单元形状确定对应的单元刚度矩阵和单元质量矩阵。（3）将单元刚度矩阵和单元质量矩阵拼接成结构整体刚度矩阵和整体质量矩阵。对于单元所处坐标系不同的工程问题，还需要进行坐标转换。（4）外力荷载的移置。针对外力的作用情况，依据静力等效原理将外力向单元结点处移置。（5）获取结构的初始位移、速度和加速度向量以及选取时间步长Δt。（6）根据式（6）、（8）、（9）、（10）、（11）求解出结构任意步长时刻的位移、速度和加速度向量。

3 Newmark 法在MATLAB 中的应用及算例

以一座移动荷载作用下的简支梁作为算例对象（见图1），利用Newmark 法编写MATLAB 程序来得到该结构的动力响应数据。

图1 移动荷载作用下的简支梁

采用移动力模型[6]来描述移动荷载作用下的简支梁。假设梁的长度为16m，截面积为7.2m2，弹性模量为3.25×1010N/m2，惯性矩为0.864m4，密度为摘要kg/m3。假设移动荷载为6.38×105N，从0 时刻出发以60km/h 的速度匀速移动，将该梁划分为160 个梁单元，取Newmark 法时间步长为0.摘要s。

使用MATLAB 作为编程工具求解上述问题的流程如图2所示[7]。

图2 Newmark 法MATLAB 程序流程图

依据上述流程图过程设计出的MATLAB 程序得到了简支梁在移动荷载作用下的动力响应数据。如图3a 所示，约0.5s后，梁的跨中截面处出现最大挠度变形，此时荷载正好移动到跨中截面处。

图3 简支梁跨中截面不同自由度的挠度时程曲线

设计MATLAB 程序时，对每一个单元结点处划分成3 个自由度方向，分别是X方向（梁的水平方向）、Y方向（梁的挠度方向）以及梁的弯矩方向。图3b 是梁跨中截面处单元1 自由度（X方向）上的位移时程曲线，由于梁没有发生水平方向上的位移，故所有时刻的位移均为0。由此可从一定程度上验证该方法的有效性。

4 讨论与展望

基于有限单元的动力响应分析法的求解精度与

结构划分的单元数量密切相关，单元数量增加的同时也会增加未知场函数矩阵的计算量，此时选择MATLAB 作为处理此类问题的工具将极大地缩短计算时长，并可以提升求解精度。在考虑移动荷载作用下结构的时程分析时，移动荷载的等效也将成为影响时程数据求解精度的影响因素，本文中MATLAB 程序设计时考虑在单元数量足够的情况下将移动荷载视作为相应时刻各对应单元的中心荷载，然后移置成结点等效荷载，从而形成整个时间域上结点力矩阵，以此对移动荷载的移动过程进行模拟。

在实际工程应用中，该程序的精度要求仍存在疑问，在未来工作中，可考虑依据相关精度要求对单元划分数量进行调整，或者在移动荷载移置成为单元结点等效荷载的过程中添加除中心荷载外的其余荷载位置来逐步提高该方法的精度。