动态分析理论和计算方法

时间：2022-11-06 百科知识版权反馈

【摘要】：因而结构设计时应避免使系统固有频率与可能的载荷频率过分接近。在ABAQUS中，“＊FREQUENCY”过程用来求解系统的振型和固有频率。在求解过程中只要给出所需振型的数目和所关心的最高频率。阻尼是针对模态动力学过程定义的，阻尼是ABAQUS分析步中定义的一部分，每阶模态可以定义不同数量的阻尼。其中ω1为计算体系的基本固有频率，λ1为该频率时的阻尼比。“＊MODAL DAMPING”选项中的rayleigh参数，表示采用瑞利阻尼。

7.1.2　动态分析理论和计算方法

1）线性动态分析

如果仅考虑系统承受载荷后的长期响应，静力分析（Static analysis）就可以解决问题。但是如果加载时间很短（如桥头跳车）或者载荷本身是动态载荷（如移动载荷），就必须采用动态分析（Dynamic analysis）。

由牛顿第二运动定律可得动力学平衡方程：

M＋I－P＝0　　　　　　　　　　　　　　　（7－39）

式中　M——系统的质量；

　　　——系统的加速度；

　　　I——系统的内力；

　　　P——施加的外力。

静态分析与动态分析有两个不同点：其一，动态分析中存在惯性力M，静态分析不存在惯性力；其二，内力I的定义不同，动态分析中的内力主要源于运动（如阻尼）和结构的变形，而静态分析中内力的产生仅源于结构的变形。

图7.7　单自由度的质量－弹簧系统

（1）固有频率

单自由度振动系统，是振动分析中最简单的一种，是指在振动过程中，振系的任一瞬间形态由一个独立坐标即可确定的系统。图7.7为单自由度无阻尼的质量－弹簧系统。

单自由度的质量－弹簧系统的动力学方程为

m＋kx－P＝0　　　　　　　　　　　　　　　（7－40）

其固有频率为

如果移动质量块然后释放，该系统将以这个频率振动。假如以此频率施加外力，位移的幅度将剧烈增加，即发生了共振。

实际系统具有很多固有频率，对于多自由度系统，其固有频率和主振型可以根据系统的无阻尼自由振动方程得到，即

M＋Kx＝0　　　　　　　　　　　　　　　（7－42）

式中　Kx——弹簧的内力，其余参数同式（7－39）。

假设方程的解为，其中A为系统自由振动时的振幅向量（列阵），A＝［A₁A₂…A_n］^T，将其代入方程式（7－42），并消去x，得到主振型方程：

（K－M）A＝0　　　　　　　　　　　　　　　（7－43）

令H＝K－M，主振型方程存在非零解的条件是A的系数矩阵H的行列式必须等于0，即

上式称为特征方程。将行列式展开，得到关于（）的n阶多项式。由于质量矩阵M为正定矩阵，刚度矩阵K为正定或半正定矩阵，因此，一般可从特征方程得到）的n个大于零的正实根，称为系统的特征值。将特征值开方后得到n个ω_nr（r＝1，2，…，n），称为系统的n个固有频率（Natural frequency）。

在多数情况下，n个固有频率互不相等，特殊情况下存在零根或重根。如果n个特征值互不相等，可以将它们按照从小到大的次序排列为0＜，并分别称为第一阶（基本）固有频率、第二阶固有频率、……、第n阶固有频率。

系统的固有频率只与系统的固有物理特性（惯性和弹性）有关，而与其他条件无关。因而结构设计时应避免使系统固有频率与可能的载荷频率过分接近。

在ABAQUS中，“＊FREQUENCY”过程用来求解系统的振型和固有频率。在求解过程中只要给出所需振型的数目和所关心的最高频率。

（2）振型叠加

在线性问题中，可以应用系统的固有频率和振型来确定结构在动载荷作用下的动态响应。采用振型叠加（Modal superposition）技术，通过振型组合可以计算结构的变形，每阶振型都乘以一个标量因子。模型的位移矢量x定义为

式中　α_i——振型φ_i的标量因子。振型叠加技术在模拟小变形、线弹性材料和无接触条件的情况下有效，即必须是线性问题。

在结构动力学问题中，结构的响应往往被相对较少的几阶振型控制，在计算这类系统的响应时，应用振型叠加为较优的方法。一般在最初计算振型和频率时需花费一些时间，但在计算响应时将节省大量时间。

如果模拟非线性问题，由于分析时固有频率可能发生明显的变化，因此振型叠加将不适用。这时须对平衡方程直接积分，在振型分析时将花费很多时间。

以下问题适合于采用线性瞬态动力学分析：

①系统是线性的：线性材料，无接触条件，无非线性的几何效应。

②响应只受较少频率支配。当响应中各频率成分增加时，例如冲击，振型叠加方法的效果将大大降低。

③载荷的主要频率在可得到的固有频率范围内，以确保对载荷描述足够精确。

④由任何突加载荷所产生的初始加速度能用特征模型精确描述。

⑤系统阻尼不大。

（3）阻尼

当给予一个无阻尼系统自由振动时，则其振幅保持恒定。实际上任何系统在运动时，结构结合处的摩擦和材料的迟滞效应等均会耗散能量，振幅将逐渐减小直至振动停止，这种能量耗散即为阻尼。通常假定阻尼为黏性阻尼或阻尼正比于速度。式（7－39）转换为

式中　C——系统阻尼矩阵；

　　　——系统的速度。

在ABAQUS中，针对非阻尼系统计算其振型。对于每一振型，含阻尼固有频率和无阻尼固有频率的关系是：

式中　ω_d——阻尼特征值；

　　　ξ——临界阻尼比，ξ＝；

　　　c——该振型阻尼；

　　　c₀——临界阻尼。

对ξ较小的情形（ξ＜0.1），含阻尼系统的特征频率非常接近无阻尼系统的对应值。当ξ增大时，采用无阻尼系统的特征频率就不太准确，当ξ接近1时，采用无阻尼系统的特征频率失效。

当结构处于临界阻尼（ξ＝1）时，施加一个扰动后，结构很快就会恢复到静止的起始位置，而且不会摆动，如图7.8所示。

①ABAQUS中阻尼的定义

图7.8　结构系统中的阻尼

在ABAQUS的瞬时模态分析中，阻尼可分为：直接模态阻尼（Direct modal damping）、瑞利阻尼（Rayleigh damping）和复合模态阻尼（Composite modal damping）。

阻尼是针对模态动力学过程定义的，阻尼是ABAQUS分析步（Step）中定义的一部分，每阶模态可以定义不同数量的阻尼。

（a）直接模态阻尼

采用直接模态阻尼可以定义对应于每阶振型的阻尼比ξ。典型的取值范围是1％～10％的临界阻尼。直接模态阻尼允许精确定义每阶振型的阻尼。

“＊MODAL DAMPING”选项中的“MODAL＝DIRECT”参数，表示采用直接模态阻尼。例如，对于前十阶振型的阻尼定义为4％的临界模态阻尼，11～20阶振型的阻尼为5％的临界阻尼，则在分析步（Step）中定义为

＊MODAL DAMPING，MODAL＝DIRECT

1，10，0.04

11，20，0.05

（b）瑞利阻尼

瑞利阻尼根据瑞利（Rayleigh）提出的线性阻尼假设来确定，即［C］＝α［M］＋β［K］。其中α和β是与结构固有频率和阻尼比有关的比例常数，实际中常采用简化的取值形式α＝λ₁ω₁，β＝λ₁/ω₁。其中ω₁为计算体系的基本固有频率，λ₁为该频率时的阻尼比。例如散粒体材料（黏土）ω₁为8.2rad/s，路面材料（包括沥青混凝土、水泥稳定碎石与二灰土）ω₁为18.6rad/s。下表列出了部分材料阻尼比与瑞利阻尼系数的关系。

表7.1　瑞利阻尼系数α和β与阻尼比λ₁的关系

“＊MODAL DAMPING”选项中的rayleigh参数，表示采用瑞利阻尼。例如，对于前十阶振型规定α＝0.2525和β＝2.9×10^－3，对于11～20阶振型规定α＝0.2727和β＝3.03×10^－3，则在分析步（Step）中定义为

＊MODAL DAMPING，RAYLEIGH

1，10，0.2525，2.9e－3

11，20，0.2727，3.03e－3

提示：质量矩阵［M］和阻尼矩阵［C］的求法如下：

应用弹性动力学的Hamilton变分原理，可得到弹性系统动力问题的有限元基本方程如下：

式中　［M］、［C］、［K］——分别为系统的质量矩阵、阻尼矩阵和刚度矩阵；

　　　｛δ｝、｝、｝——分别为位移、速度和加速度向量；

　　　｛F（t）｝——动载荷。

（1）质量矩阵

从运动的结构中取出一微小部分，根据达朗贝尔原理，在它的单位体积上作用的惯性力为

式中　ρ——材料的密度。

在对结构进行离散化以后，取出一个单元，并采用如下形式的位移函数：

｛r｝＝［N］｛δ｝^e　　　　　　　　　　　　　　　（7－50）

将上式代入式（7－49），得到

再利用载荷移置的一般公式｛P［N］^T｛p_i｝dxdydz，求得作用于单元结点上的惯性力如下：

可见单元质量矩阵为

［m］＝∫［N］^Tρ［N］dV　　　　　　　　　　　　　　　（7－53）

（2）阻尼矩阵

如果假定阻尼力正比于质点的运动速度，从运动的结构中取出一微小部分，在它单位体积上作用的阻尼力是：

式中　a——比例常数；

　　　ρ——材料密度；

　　　N——形函数。

利用载荷移置的一般公式｛P［N］^T｛p_i｝dxdydz，可求得单元结点上的阻尼力如下：

那么有：

［c］＝a∫［N］^Tρ［N］dV＝a［m］　　　　　　　　　　　　　　　（7－56）

由上式可见，单元阻尼矩阵［c］正比于单元质量矩阵［m］。

如果假定阻尼力正比于应变速度，则阻尼力可表示为

作用于单元结点上的阻尼力为

其中，

［c］＝β∫［B］^T［D］［B］dV＝β［k］　　　　　　　　　　　　　　　（7－59）

从上式可以看出，单元阻尼矩阵［c］正比于单元刚度矩阵［k］。

通常根据实测资料，由振动过程中结构整体的能量消耗来决定阻尼的近似值，因此不是计算单元的阻尼矩阵，而直接计算结构的整体阻尼矩阵［C］。一般采用如下的线性关系式，称为瑞利阻尼，即

［C］＝α［M］＋β［K］　　　　　　　　　　　　　　　（7－60）

（c）复合模态阻尼

在复合模态阻尼中，对应于每种材料的阻尼，定义其与临界阻尼的比例，并且复合模态阻尼是对应于整体结构的阻尼。

②选择阻尼值

在大多数线性动力学问题中，为了获得精确的结果，阻尼的确定非常重要。但是阻尼只是在结构吸收能量特性意义上的近似，而不是去模拟造成这种效果的物理机制。所以，确定模型中需要的阻尼数据非常困难，除了偶尔从动态试验中获得这些数据，大多数需要通过经验或参考资料获得数据。这时应通过参数分析来评价阻尼系数对于模拟过程的敏感性。

提示：ABAQUS中所有单元均可用于动态分析。选取单元的一般原则与静态模拟相同。但是，在模拟冲击和爆炸载荷时，应选用一次单元。因为它们具有集中质量格式，在模拟应力波效果方面优于采用二次单元的一致质量格式。

2）显式非线性动态分析

ABAQUS/Standard中的非线性动态程序采用隐式的时间积分。这适用于非线性结构的动态问题。例如，某一突发事件激发了结构的动态响应，像冲击，以及在结构响应中，包括由于塑性或黏性阻尼造成大量的能量耗散。在这些研究时，初始的高频响应非常重要，振动由模型的耗散机制而迅速衰减。

非线性动态分析的另一种办法是ABAQUS/Explicit提供的显式动态程序。显示算法以应力波的传播方式在模型中传播结果，一次一个单元地传播。因而，它最适合于求解应力波影响为主的问题以及模拟时间很短（不超过1s）的问题。

显式算法的另一优点是它能模拟不连续的非线性，例如接触和失效问题，并且比采用ABAQUS/Standard更容易。而对于大型、高度不连续的问题，即使响应是准静态的（Quasi－static）的，采用ABAQUS/Explicit来模拟也更容易些。

从用户的观点来看，显示算法与隐式算法的区别是：

①显式算法需要很少的时间增量步，它仅依赖于模型的最高固有频率，而与载荷的类型和持续的时间无关。通常的模拟需要取10 000～1 000 000个增量步，每个增量步的计算成本相对较低。

②隐式算法对时间增量步的大小没有内在的限制，增量的大小通常取决于精度和收敛情况。典型的隐式模拟所采用的增量步数目要比显式模拟小几个数量级。然而，在每个增量步中必须求解一套完整的方程组，所以对于每一增量步的成本，隐式算法远高于显式算法。

（1）动力学显式有限元方法

ABAQUS/Explicit应用中心差分法对运动方程进行显式的时间积分，应用一个增量步的动力学条件计算下一个增量步的动力学条件。在增量步开始时，程序求解动力学平衡方程，表示为用结点质量矩阵M乘以结点加速度等于结点的合力（所施加的外力P与单元内力I之间的差值）：

＝P－I　　　　　　　　　　　　　　　（7－61）

在当前增量步开始时（t时刻），结点的计算加速度为

采用中心差分方法对结点加速度在时间上进行积分，在计算速度的变化时假定加速度为常数。应用速度的变化量加上前一个增量步中点的速度来确定当前增量步中点的速度：

速度对时间的积分加上在增量步开始时的位移以确定增量步结束时的位移：

可以看出，在增量步开始时提供了满足动力学平衡条件的加速度。根据加速度，可以显式地前推速度和位移。显式即指在增量步结束时的状态仅依赖于该增量步开始的位移、速度和加速度。为了使该方法产生精确的结果，时间增量必须相当小，这样增量步的加速度几乎为常数。

（2）自动时间增量和稳定性

稳定性限制了ABAQUS/Explicit求解器所采用的最大时间步长，这是应用ABAQUS/Explicit进行计算的一个重要因素。

应用显式方法，基于在增量步开始时t的模型状态，通过时间增量Δt前推到当前时刻的模型状态。如果时间增量大于这个最大的时间步长，则此时间增量已经超出了稳定性限制（Stability limit）。超过稳定性限制的一个可能后果就是数值不稳定，它可能导致解答不收敛。由于一般不可能精确地确定稳定性限制，因而采用保守的估计值。为了提高计算效率，ABAQUS/Explicit选择时间增量，使其尽可能地接近而且不超过稳定性限制。

以在系统中的最高频率（ω_max）的形式定义稳定性限制。无阻尼的稳定性限制由下式定义：

而有阻尼的稳定性限制由下式定义：

式中　ξ——最高频率振型的临界阻尼值。

临界阻尼定义了在自由的和有阻尼的振动关系中有振荡运动与无振荡运动之间的限制。为了控制高频振荡，ABAQUS/Explicit总是以体积黏性的形式引入一个小量阻尼。

（3）动态振荡的阻尼

在模型中加入阻尼有两个原因：限制数值振荡或为系统增加物理的阻尼。ABAQUS/Explicit提供了在分析中加入阻尼的方法。

①体黏性

体黏性引入了与体积应变相关的阻尼，其目的是改进对高速动力学事件的模拟。ABAQUS/Explicit包括体黏性的线性与二次的形式。可以在定义分析步时修改默认的体黏性参数。

（a）线性体黏性

默认情况下，线性体黏性来自阻尼在单元最高阶频率中的振荡。由下式可得到与体积应变率成线性关系的体黏性压力：

p₁＝b₁ρc_dL^e　　　　　　　　　　　　　　　（7－67）

式中　b₁——阻尼系数，默认值为0.06；

　　　ρ——材料密度；

　　　c_d——膨胀波速；

　　　L^e——单元的特征长度；

　　　——体积应变速率。

（b）二次体黏性

仅在实体单元中（除平面应力单元CPS4R外）包括二次体黏性，并且只有当体积应变速率可压缩时采用。由下式得到二次体黏性压力：

式中　b₂——阻尼系数，默认值为1.2。

（c）基于体黏性的临界阻尼比

体黏性压力只是基于每个单元的膨胀模式。在最高阶单元模式中的临界阻尼比由下式给出：

式中　ξ——临界阻尼比。线性项单独代表了6％的临界阻尼，而二次项一般是更小的量。

②黏性压力

黏性压力载荷一般应用在结构问题或准静态问题，以阻止低阶频率的动态影响，从而以最少数目的增量步达到静态平衡。由下式的分布载荷形式施加：

式中　p——施加的压力；

　　　c_v——黏度，在数据行中作为载荷的量值给出；

　　　——施加在黏性压力面上点的速度矢量；

　　　——该点处表面上的单位外法线矢量。

在典型的情况下，使当前动力影响最小化的有效方法是将c_v设置为ρc_d的很小的百分数（1％或2％）。

③材料阻尼

材料模型可能以塑性耗散或黏弹性的形式提供阻尼，也可以使用瑞利阻尼。与瑞利阻尼相关的阻尼系数有两个：质量比例阻尼α_R和刚度比例阻尼β_R。

（a）质量比例阻尼α_R

α_R因子定义了与单元质量矩阵成比例的阻尼。引入的阻尼力源于在模型中结点的绝对速度。合理的质量比例阻尼不会明显地降低稳定极限。

（b）刚度比例阻尼β_R

β_R因子定义了一个与弹性材料刚度成比例的阻尼。对于任何非线性分析，都可以引入阻尼；而对于线性分析，提供了标准的瑞利阻尼。对于线性分析，刚度比例阻尼等于β_R乘以刚度矩阵。必须慎重地使用刚度比例阻尼，因为它可能使稳定极限明显地降低。为了避免大幅度地降低稳定时间增量，刚度比例阻尼应该小于或等于未考虑阻尼时的初始时间增量的量级。

（c）离散减振器

阻尼定义的另一种选择是定义单独的减振器单元。每个减振器单元提供了一个与它的两个结点之间的相对速度成正比的阻尼力。其优点是能够把阻尼施加在有必要施加的结点上。减振器应当与其他单元并行使用，例如弹簧，其不会引起稳定极限的明显下降。

（4）能量平衡

①能量平衡的表述

整体模型的能量平衡可写为

E_total＝E_t＋E_V＋E_FD＋E_KE－E_W　　　　　　　　　　　　　　　（7－71）

式中　E_t——内能；

　　　E_V——黏性耗散能；

　　　E_FD——摩擦耗散能；

　　　E_KE——动能；

　　　E_W——外加载荷所做的功。

这些能量分量的总和为E_total，是个常数。在数值模拟中，E_total只是近似的常数，一般有小于1％的误差。