CN109902404B - 不同阻尼形式的结构时程响应积分的统一递推计算方法 - Google Patents

Abstract

本发明公开了不同阻尼形式的结构时程响应积分的统一递推计算方法，包括如下步骤：对实际结构进行有限元离散，建立结构时程响应分析的有限元离散方程模型并进行模态分析，根据模态分析结果进行解耦，获取模态坐标系下的单自由度系统方程；获得单自由度系统方程的位移理论解，根据位移理论解的形式，构造与其形式一致的中间变量，与位移理论解共同作为中间离散递推变量；获取中间离散递推变量在欠阻尼、临界阻尼以及过阻尼形式下计算统一离散递推表达式；计算在不同时程载荷作用下的动力学响应并通过模态组合获得结构实际响应。本发明可以结合并行程序架构很容易实现并行计算，从而提升数值分析的计算规模和可扩展性。

Description

不同阻尼形式的结构时程响应积分的统一递推计算方法

技术领域

本发明属于结构动力学领域，主要解决各类工程结构承受随时间变化外力载荷下的结构动力学响应问题，尤其涉及针对不同阻尼形式的结构在时程载荷作用下的动力学响应计算方法。

背景技术

结构动力分析常用来确定时变载荷对整个结构或部件的影响，同时需要考虑阻尼及惯性效应的作用。结构时程响应分析尤其涉及规模较大的动力时程响应分析一般都采用模态叠加法开展相关分析，通过模态解耦后的各个单自由度时程响应分别计算，再进行模态叠加得到真实结构的的动力学响应，针对模态解耦后的单自由度系统，需要通过数值积分获取其单自由度系统的响应特征，常用的积分算法有Newmark、Wilson-θ等方法，这些积分算法受离散格式等影响，随着计算时间的增加会出现大量的累积误差，计算时间长并且难以获取准确的结构响应特征。

发明内容

本发明的目的就在于为了解决上述问题而提供不同阻尼形式的结构时程响应积分的统一递推计算方法，包括如下步骤：

S1:对实际结构进行有限元离散，建立结构时程响应分析的有限元离散方程模型并进行模态分析，根据模态分析结果进行解耦，获取模态坐标系下的单自由度系统方程；

S2:获得单自由度系统方程的位移理论解，根据位移理论解的形式，构造与其形式一致的中间变量，与位移理论解共同作为中间离散递推变量；

S3:利用数值积分获取中间离散递推变量在欠阻尼、临界阻尼以及过阻尼形式下计算统一离散递推表达式；

S4:计算在不同时程载荷作用下的动力学响应,获得速度和加速度不同时刻的递推解，并通过模态组合获得结构实际响应。

进一步的，所述S1包括如下步骤：

S11:设M：n×n阶质量矩阵；C：n×n阶阻尼矩阵；K：n×n阶刚度矩阵；x(t)：n×1阶结构位移响应向量；f(t)：n×1阶离散载荷向量；

建立结构时程响应分析的有限元离散方程，经过有限元离散后，获取n个自由度的工程结构在外载荷作用下的振动时域运动方程：

S12:进行模态分析，设：ω为系统固有圆频率，φ为n×1阶固有振型向量，对应广义特征值方程为：

Mφ＝ω²Kφ；

求解后得到m阶最小特征对：

(ω₁,φ₁),(ω₂,φ₂),...,(ω_m,φ_m)；

设：Φ为n×m阶振型矩阵；u(t)为m×1阶模态坐标向量；φ_i为第i阶模态振型向量(n×1)；u_i(t)为第i阶模态坐标；

利用模态振型进行坐标变换：

设：第i阶振型阻尼比为ξ_i，M^*为m×m阶模态质量矩阵，M^*＝Φ^TMΦ；K^*为m×m阶模态刚度矩阵，

C^*为m×m阶模态阻尼矩阵，C^*＝Φ^TCΦ＝diag{2ω₁ξ₁,2ω₂ξ₂,...,2ω_mξ_m}；p_i(t)为第i阶模态广义力；

模态坐标系下m个解耦的单自由度系统方程：

进一步的，设第i阶的模态参与因子γ_i＝φ_i ^Td，载荷f(t)位置处单位载荷作用对应的载荷向量为d，t时刻的载荷值为f(t)，其特征在于，模态广义力具体表达式为：

p_i(t)＝φ_i ^Tdf(t)＝γ_if(t)。

进一步的，所述S2包括如下步骤：

对解耦后的单自由度模态运动方程，利用Duhamel理论公式，得到位移响应表达式，进行一阶和二阶求导，得到第i阶模态对应的速度响应、加速度响应：

速度响应：

加速度响应：

构造与解耦后的单自由度模态运动方程解的形式一致的中间变量v_i(t)为：

进一步的，S3包括如下步骤：

采用Simpson积分在[a，b]区间进行积分数值求解：

对位移理论解和中间变量中的积分项进行积分，所述离散递推表达式为：

(1)欠阻尼情况下离散递推表达式：

其中：

(2)临界阻尼情况下离散递推表达式：

其中：

(3)过阻尼情况下离散递推表达式：

其中：

进一步的，所述动力学响应为：

本发明的有益效果在于：

(1)将解耦后单自由度系统的Duhamel理论解与Simpson数值积分相结合，通过二次函数的数值积分格式提升了数值求解精度；

(2)通过巧妙构建一个与位移理论解结构形式一致的中间变量，从而更好地建立了数值积分递推表达式。最终形成的两个数值积分递推表达式具有一定的对称性，便于各类阻尼情形下的数值积分的推导，具有较强创新性；

(3)首次提出一种针对不同阻尼形式的结构时程响应积分统一递推算法。针对各种阻尼情形，推导形成了具有统一格式的、结构形式对称、简洁明了的离散递推公式，并针对时程响应分析构建了具体的求解算法和技术实施路线；

(4)提出的计算方法易于程序实现，可以很方便地根据递推公式和算法构建相应的计算程序，实用性强。而且，该方法可以结合并行程序架构很容易实现并行计算，从而提升数值分析的计算规模和可扩展性；

(5)本发明采用并行计算平台予以具体实现，从结构的有限元模型离散、矩阵组装到模态分析，再到结构的时程响应分析全过程均在并行计算环境之中进行，对称统一的积分求解算法结合并行计算平台，大大提升了解决复杂工程问题时程响应分析的能力和计算效率，计算规模可以扩展到数亿自由度以上，可扩展的并行CPU核数可达到上万个，计算能力超目前商业有限元软件时程响应分析能力约两个量级；

(6)本发明提出的技术路线能够广泛应用于典型武器系统发射和再入大气层过程中承受随时间变化的空气脉动压力时程响应分析，显著提升了数值分析能力，对于刻画武器系统飞行过程中关注部位的时间历程响应特征和结构优化发挥了重要作用；

(7)此外，利用该方法可以有效分析航空航天结构承受空气脉动压力时变载荷的响应特征、路面基础激励对车辆等交通工具产生的时变载荷响应特征，以及各类民用设施在承受随时间变化的风载或地震载荷下的结构响应特征等，可为各类军民领域中工程结构的设计改进、优化以及动力学评估等提供了数值仿真工具和模拟依据；根据数值分析结果，可以避免不利的结构设计，防止结构重要关注部位响应过大而导致的产品结构破坏或失效，可有效降低试验成本并加快产品研发速度。

附图说明

图1是不同阻尼形式的结构时程响应积分的统一递归计算方法的流程图。

具体实施方式

在采用模态叠加进行时程响应分析时，模态分析流程负责对结构的模态分析，通过与单元、材料、载荷、约束等有限元分析模块协助，形成并求解特征值方程组，最终得到结构的固有频率和模态振型。在此基础上，再采用不同的时程积分方法进行模态叠加获得结构的响应特征。

针对具体工程结构和承受的时程载荷条件，采用有限元建模方式建立其有限元模型，结合材料参数进行有限元离散，根据模型离散构建n×n阶质量矩阵M和刚度矩阵K，基于质量矩阵M和刚度矩阵K构建模态分析广义特征值方程，进行模态分析计算，获取前m阶模态频率和模态振型(ω₁,φ₁),(ω₂,φ₂),...,(ω_m,φ_m)，该步骤可以借助各类有限元分析软件很容易实现。

根据时程载荷作用区域进行载荷有限元离散，构造各载荷分布的方向向量d，结合第一步求解得到的各阶模态振型，计算时程载荷对应的各阶模态参与因子γ_i＝φ_i ^Td,结合时程载荷具体大小和给定的各阶模态阻尼比，通过解耦构建模态坐标下的m个单自由度系统方程。

对各个单自由度系统方程，根据构造的中间变量v_i(t)，选取时间计算离散步长，按照结构t＝0时刻的初始条件和不同阻尼类型，根据位移和中间变量的离散递推表达式，计算模态坐标系下各单自由度系统在不同离散时刻的位移、速度、加速度等响应特征。

根据模态坐标下计算获取的各单自由度系统位移、速度和加速度等响应特征，通过模态叠加获取结构关注点在物理坐标系下不同时刻的实际响应特征，输出关注点的时域响应曲线。

下面结合附图对本发明作进一步说明：

如图1所示，本发明不同阻尼形式的结构时程响应积分的统一递推计算方法，包括如下步骤：

S1:对实际结构进行有限元离散，建立结构时程响应分析的有限元离散方程模型并进行模态分析，根据模态分析结果进行解耦，获取模态坐标系下的单自由度系统方程：

动力学有限元分析的基本思路是通过对结构物理模型进行有限元离散，并对动力学离散方程进行求解，最终获取所需的物理量并对结构特性进行评估。

设M：n×n阶质量矩阵；C：n×n阶阻尼矩阵；K：n×n阶刚度矩阵；x(t)：n×1阶结构位移响应向量；f(t)：n×1阶离散载荷向量；经过有限元离散后，n个自由度的工程结构在外载荷作用下的振动时域运动方程可用下式进行描述：

对上式进行模态解耦，需先进行模态分析，其对应的广义特征值方程为：

Mφ＝ω²Kφ；

式中：ω为系统固有圆频率；φ为n×1阶固有振型向量；

通过对式(2)所对应的广义特征值方程求解，可以得到m阶最小特征对：

(ω₁,φ₁),(ω₂,φ₂),...,(ω_m,φ_m)；

设：Φ为n×m阶振型矩阵；u(t)为m×1阶模态坐标向量；φ_i为第i阶模态振型向量(n×1)；u_i(t)为第i阶模态坐标。

利用振型进行坐标变换：

当m＝n时等式严格成立，但大规模计算时需要计算的模态阶数一般远小于矩阵维数(自由度数)，因此取近似。模态振型取为关于质量矩阵M的正交归一化振型，即：

经过上述模态坐标变换后，将物理坐标系中的多自由度运动方程解耦转换成模态坐标系下的多个单自由度运动方程。

设：第i阶振型阻尼比为ξ_i，M^*为m×m阶模态质量矩阵，M^*＝Φ^TMΦ；K^*为m×m阶模态刚度矩阵，

C^*为m×m阶模态阻尼矩阵，C^*＝Φ^TCΦ＝diag{2ω₁ξ₁,2ω₂ξ₂,...,2ω_mξ_m}；p_i(t)为第i阶模态广义力；设第i阶的模态参与因子γ_i＝φ_i ^Td，载荷f(t)位置处单位载荷作用对应的载荷向量为d，t时刻的载荷值为f(t)，其特征在于，模态广义力具体表达式为：p_i(t)＝φ_i ^Tdf(t)＝γ_if(t)。

模态坐标系下m个解耦的单自由度系统方程：

S2:获得单自由度系统方程的位移理论解，根据位移理论解的形式，构造与其形式一致的中间变量，与位移理论解共同作为中间离散递推变量：

对于解耦后的单自由度模态运动方程，初始条件下t＝0，u_i(t)＝u_i0，

设ω_Di为第i阶模态阻尼固有频率，

利用Duhamel理论公式，可以得到其理论解为：

进行一阶求导，可以得到第i阶模态对应的速度响应为：

进行二阶求导，可以得到第i阶模态对应的加速度响应为：

上述得到模态解耦后单自由度系统的精确理论解，但在实际有限元分析计算中，由于模态广义载荷的复杂性和多样性，通常对Duhamel积分项无法获取积分理论表达式，仍需要借助数值积分进行离散，获取积分项的数值解。

采用的数值求解方法一般有矩形法、梯形法和Simpson数值积分法三种。矩形法采用一系列水平直线段(0次函数)近似代替原来的各曲线段，梯形法采用斜/直线段(1次函数)代替曲线段，而Simpson积分则是采用一系列二次函数来近似代替被积函数曲线段，故又叫抛物线法。Simpson积分和其他积分方法相比具有更高的求解精度。

Simpson积分数值计算表达式如下所示：

本发明通过巧妙引入一个形式与位移理论解结构形式相同的中间变量，实现数值积分求解并形成针对时间历程的递推表达式。

设中间变量v_i(t)：

则：

S3:利用数值积分获取中间离散递推变量在欠阻尼、临界阻尼以及过阻尼形式下计算统一离散递推表达式：

(1)欠阻尼情况下离散递推表达式：

其中：

(2)临界阻尼情况下离散递推表达式：

其中：

(3)过阻尼情况下离散递推表达式：

其中：

S4:计算在不同时程载荷作用下的动力学响应,获得速度和加速度不同时刻的递推解，并通过模态组合获得结构实际响应：

获得中间变量v_i(t+Δt)的数值解后，便可根据中间变量v_i(t)得到速度的数值解：

根据加速度响应获得递推时刻t+Δt的加速度表达式：

设φ_ik为第i阶模态对应第k个自由度的元素值，对于第k个自由度，可以得到模态叠加后的响应为：

有限元模型离散和模态分析是本发明进行时域积分计算的前提条件，可以借助普通的有限元软件或自主程序轻松实现。本发明的核心是振动时域积分求解，可以计算结构关注点的时程响应曲线。若在有限元模型离散和模态分析过程中采用区域分解的并行计算模式，整个流程可以很容易扩展到并行求解中。

本发明具有以下优点：

(1)将解耦后单自由度系统的Duhamel理论解与Simpson数值积分相结合，通过二次函数的数值积分格式提升了数值求解精度；

(2)通过巧妙构建一个与位移理论解结构形式一致的中间变量，从而更好地建立了数值积分递推表达式。最终形成的两个数值积分递推表达式具有一定的对称性，便于各类阻尼情形下的数值积分的推导，具有较强创新性；

(3)首次提出一种针对不同阻尼形式的结构时程响应积分统一递推算法。针对各种阻尼情形，推导形成了具有统一格式的、结构形式对称、简洁明了的离散递推公式，并针对时程响应分析构建了具体的求解算法和技术实施路线；

(4)提出的计算方法易于程序实现，可以很方便地根据递推公式和算法构建相应的计算程序，实用性强。而且，该方法可以结合并行程序架构很容易实现并行计算，从而提升数值分析的计算规模和可扩展性；

(5)对于以上技术路线，本发明采用并行计算平台予以具体实现，从结构的有限元模型离散、矩阵组装到模态分析，再到结构的时程响应分析全过程均在并行计算环境之中进行，对称统一的积分求解算法结合并行计算平台，大大提升了解决复杂工程问题时程响应分析的能力和计算效率，计算规模可以扩展到数亿自由度以上，可扩展的并行CPU核数可达到上万个，计算能力超目前商业有限元软件时程响应分析能力约两个量级；

(6)本发明提出的技术路线能够广泛应用于典型武器系统发射和再入大气层过程中承受随时间变化的空气脉动压力时程响应分析，显著提升了数值分析能力，对于刻画武器系统飞行过程中关注部位的时间历程响应特征和结构优化发挥了重要作用；

(7)此外，利用该方法可以有效分析航空航天结构承受空气脉动压力时变载荷的响应特征、路面基础激励对车辆等交通工具产生的时变载荷响应特征，以及各类民用设施在承受随时间变化的风载或地震载荷下的结构响应特征等，可为各类军民领域中工程结构的设计改进、优化以及动力学评估等提供了数值仿真工具和模拟依据；根据数值分析结果，可以避免不利的结构设计，防止结构重要关注部位响应过大而导致的产品结构破坏或失效，可有效降低试验成本并加快产品研发速度。

本发明的技术方案不限于上述具体实施例的限制，凡是根据本发明的技术方案做出的技术变形，均落入本发明的保护范围之内。

Claims (3)

1.不同阻尼形式的结构时程响应积分的统一递推计算方法，其特征在于，包括如下步骤：

S1:对实际结构进行有限元离散，建立结构时程响应分析的有限元离散方程模型并进行模态分析，根据模态分析结果进行解耦，获取模态坐标系下的单自由度系统方程；

S2:获得单自由度系统方程的位移理论解，根据位移理论解的形式，构造与其形式一致的中间变量，与位移理论解共同作为中间离散递推变量；

S3:利用数值积分获取中间离散递推变量在欠阻尼、临界阻尼以及过阻尼形式下计算统一离散递推表达式；

S4:计算在不同时程载荷作用下的动力学响应，获得速度和加速度不同时刻的递推解，并通过模态组合获得结构实际响应；

所述S1包括如下步骤：

S11:设M：n×n阶质量矩阵；C：n×n阶阻尼矩阵；K：n×n阶刚度矩阵；x(t)：n×1阶结构位移响应向量；f(t)：n×1阶离散载荷向量；

建立结构时程响应分析的有限元离散方程，经过有限元离散后，获取n个自由度的工程结构在外载荷作用下的振动时域运动方程：

S12:进行模态分析，设：ω为系统固有圆频率，φ为n×1阶固有振型向量，对应广义特征值方程为：

Mφ＝ω²Kφ；

求解后得到m阶最小特征对：

(ω₁,φ₁),(ω₂,φ₂),...,(ω_m,φ_m)；

设：Φ为n×m阶振型矩阵；u(t)为m×1阶模态坐标向量；φ_i为第i阶模态振型向量(n×1)；u_i(t)为第i阶模态坐标；

利用模态振型进行坐标变换：

设：第i阶振型阻尼比为ξ_i，M^*为m×m阶模态质量矩阵，M^*＝Φ^TMΦ；K^*为m×m阶模态刚度矩阵，

C^*为m×m阶模态阻尼矩阵，C^*＝Φ^TCΦ＝diag{2ω₁ξ₁,2ω₂ξ₂,...,2ω_mξ_m}；p_i(t)为第i阶模态广义力；

模态坐标系下m个解耦的单自由度系统方程：

所述S2包括如下步骤：

对解耦后的单自由度模态运动方程，利用Duhamel理论公式，得到位移响应表达式，进行一阶和二阶求导，得到第i阶模态对应的速度响应、加速度响应：

速度响应：

加速度响应：

构造与解耦后的单自由度模态运动方程解的形式一致的中间变量vi(t)为：

所述S3包括如下步骤：

采用Simpson积分在[a，b]区间进行积分数值求解：

对位移理论解和中间变量中的积分项进行积分，所述离散递推表达式为：

(1)欠阻尼情况下离散递推表达式：

其中：

(2)临界阻尼情况下离散递推表达式：

其中：

(3)过阻尼情况下离散递推表达式：

其中：

2.根据权利要求1所述不同阻尼形式的结构时程响应积分的统一递推计算方法，设第i阶的模态参与因子γ_i＝φ_i ^Td，载荷f(t)位置处单位载荷作用对应的载荷向量为d，t时刻的载荷值为f(t)，其特征在于，模态广义力具体表达式为：

3.根据权利要求1所述不同阻尼形式的结构时程响应积分的统一递推计算方法，其特征在于，所述动力学响应为：

Images