疲劳损伤谱时域、频域计算方法 及其等效性验证

高博，张忠，王帅，钟嫄，李海波

（北京强度环境研究所，北京 100076）

运载火箭、导弹、卫星等航天飞行器在发射过程中要经历严酷的力学环境。美国NASA的一项研究结果表明，在首次发射卫星产生的破坏或失效现象中，45%是由于火箭主动段飞行时恶劣的动力学环境引起的。对于低频瞬态环境，目前主要采用冲击响应谱转化为正弦扫频振动试验来进行模拟，仅考虑了峰值破坏等效效应。正弦扫频试验的扫描率对于验收试验取4 oct/min，对于鉴定试验取2 oct/min。即使采用4 oct/min扫描率，持续时间也需1 min，而实际低频瞬态持续时间在秒量级，因此在累计损伤方面存在严重的过试验。本研究拟从累积损伤等效方面考虑，提出基于疲劳损伤谱的环境分析方法，为提升环境试验精细化水平提供支撑。

1 概述

1.1 响应谱

响应谱是一种表征冲击和振动环境严酷度的方法。响应谱描述的是环境对一组单自由度弹簧振子的影响，通过建立一组具有一系列固有频率（fn）和统一品质因子（Q）的单自由度振子模型，计算各振子对振动环境的响应，并绘制响应的某些特征（例如峰值响应）与固有频率的关系图，即得到响应谱。其中应用最广泛的响应谱是冲击响应谱（SRS），冲击响应谱是振子对冲击环境瞬态时间历程的峰值响应。通常情况下，冲击响应谱是通过对振子离散时间模型的时域仿真计算出来的[1]。另一个与之类似的响应谱是由功率谱密度函数描述的高斯随机环境下，振子的绝对加速度峰值谱，称为振动响应谱（VRS）[2]。通过对振子的频率响应进行积分，得到均方根响应，可在频域内计算出振动响应谱。

用冲击响应谱描述振动环境的一个不足在于，一个基础加速度环境与计算得到的响应谱不是一一对应的，即不同的基础加速度环境可能计算出相同的冲击响应谱。如果由冲击响应谱描述冲击试验条件，则可能使用了比预期的更良性的时间历程来测试部件，导致待测部件“欠试验”。针对此问题，冲击环境的冲击响应谱描述通常需要补充关于瞬态环境的时间信息。通常以带限时间矩的形式出现。a时刻的第i阶时间矩mi(a)，是关于时间历程f(t)的函数，由式（1）给出。在加速度环境为高斯随机且加速度功率谱密度函数在半功率带宽中平坦的假设下，振动响应谱具有与其相关联的唯一的由加速度功率谱密度函数描述的环境。在这种情况下，对于具有自然频率fn、品质因子Q以及自然频率处APSD量级为Pa(fn)的振子，基础加速度的APSD和绝对峰值加速度VRS之间的关系由熟悉的Miles公式（2）给出。

1.2 疲劳损伤谱

虽然SRS和VRS可以很好地描述基础加速度环境诱导结构峰值响应的趋势，但它们缺乏环境诱导结构振动循环次数的信息。疲劳损伤谱与冲击响应谱相似，区别在于FDS描述的是疲劳损伤与频率之间的关系。引入该概念是为了检验由加速度功率谱密度（APSD）曲线和持续时间描述的频域环境的严苛程度。该方法要求与APSD关联的过程必须是一个“强耦合”随机过程。文献[3]提出了一种用于评价电动、电液振动台疲劳性能的方法，称为损伤势法，该方法被广泛用于各种振动台疲劳性能的评价[4]。随着理论的发展，有学者提出了基于疲劳损伤谱的随机振动加速试验方法，主要用于非平稳环境的载荷谱制定及加速试验，该方法已广泛应用于各类机械电子产品的加速试验中[5-7]。

2 疲劳损伤谱机理研究

2.1 疲劳损伤谱时域计算方法

在时域中计算FDS的步骤并不复杂，但计算量较大。方法与SRS类似，计算步骤如下：

1）获得描述振动环境的加速度时程a(t)。去除趋势项，如有必要对数据增采样。对于每个具有固有频率fn和质量因子Q的单自由度振子进行后续计算。

2）计算伪速度响应pv(t)。

3）计数获得伪速度的雨流循环谱n(PV)。

4）根据Minor法则和S-N曲线（设疲劳指数b）计算振子累积损伤D(fn,Q,b)。

5）绘制以振子累积损伤D(fn,Q,b)为纵坐标，以振子固有频率fn为横坐标的图。这里的累积损伤谱D(fn,Q,b)就被称为疲劳损伤谱。

2.1.1 加速度时程的获取——去趋势项

在设计结构荷载时，稳态加速度和低频瞬态加速度载荷是需要单独考虑的，通常是在这些载荷的基础上增加随机振动来推导结构荷载。因此，在疲劳损伤谱的计算中，须从时程中去除稳态和极低频瞬态信号。这可以通过名为“detrending”的通用程序来完成。从数据中去除极低频趋势项的方法包括：高通数字滤波、傅里叶或小波滤波、拟合滤波多项式法等。本项目采用分段样条曲线对数据进行趋势拟合。通过从数据中减去拟合曲线，从而消除趋势项。该方法是一种适用于不同类型数据和相关趋势的稳健方法。在MATLAB中，该方法是通过计算20个左右的移动平均点，并使用“spline”选项的interpl命令增采样到原始数据长度来实现的。

2.2.2 振子的伪速度响应计算

研究表明，与绝对加速度等其他响应相比，伪速度与应力的关系更为密切[8-11]。伪速度与应力的关系大致成正比，而疲劳损伤与应力周期相关。因此，必须根据每个单自由度振子的伪速度响应来计算周期，从而得到有意义的疲劳损伤谱。式（3）给出了加速度和伪速度之间的传递函数，其中ζ是阻尼比（ζ= 1(2Q)），ωn是固有圆频率（ 2πfn）。

快速计算伪速度响应的一种方法是将其转换为离散时间滤波模型。伪速度模型由式（4）给出。伪速度响应采样的滤波器模型pv(k)与基础加速度a(k)的关系如式（5）所示。采样标号可以通过t(k)=kTS转换为时间。利用在MATLAB中的filter命令可实现该递推公式。

式中：

其中TS是加速度时间历史的采样间隔。

结果表明，由于斜坡不变法的低通效应，如果输入加速度激励的固有频率大于采样频率的17%（fs=1/TS），则峰值加速度可能有10%的误差[12]。因此，建议采样率大于加速度时程中最大频率分量的10倍（fs≥10fmax）。如果加速度历程没有以足够高的采样速率，则可以通过增采样（插入零法）、FIR滤波、减采样对其重新采样。这可以通过在MATLAB中使用upfirdn或resample函数来实现。

2.2.3 雨流循环计数

如前所述，结构的疲劳损伤需要应力循环计数。由于许多结构的伪速度与应力大致成正比，因此可以通过计算伪速度响应的循环计数得到应力循环计数。周期计数的定义方法有许多，包括峰谷、平交、距离对、雨流等[13]。雨流循环计数定义与结构损伤关系最为密切，它相当于材料应力应变平面上的闭合滞回线。在文献[14-17]中，使用了对窄带信号有效的更简化的周期定义。然而，当振子的固有频率偏离加速度环境中的谱峰时，就违反了窄带假设。这导致窄带周期计数与对宽带信号有效的雨流周期计数器相比，给出的结果较差。

雨流循环有许多等价的定义[13]。对于文中的结果，使用四点算法[18]计算雨流周期。循环计数算法从时间历史的一系列极值（峰和谷）开始。四点算法一次考虑4个连续的时间历程极值（S1,S2,S3,S4）。三个连续范围由以下项构成：。如果 ΔS2小于或等于其相邻范围，则将其计算为一个循环，并从极值序列中移除其极值。不符合此定义的极值留在序列残差中，并作为半周期计算。有关雨流计数的更多信息，可参考文献[13]。

2.2.4 计算累计损伤

由循环谱计算疲劳损伤D需要将Palgrem-Minor规则（式（6））与S-N关系（式（7））相结合。

式中：ni为具有应力幅度Si的循环次数；Ni为具有应力幅度Si的循环引起失效时的次数；c为比例常数，b为疲劳指数，一般为4～25。常数c可以设置为单位值，因为比例不影响FDS用于环境比较。

如前所述，伪速度与应力的关系可表示为：σ(t)=k∙pv(t)。由此得出以下损伤方程：

式中：ni是伪速度振幅PVi的循环数。计算每个振子固有频率的伪速度响应、周期计数和损伤，最终形成FDS、D(fn,Q,b)。

2.2 疲劳损伤谱频域计算方法

本节将从加速环境的APSD计算疲劳损伤谱D(fn,Q,b)。频域FDS的计算量远小于时域FDS，为了使计算有效，与APSD相关联的响应历程必须是“强耦合”的平稳随机过程。此外，在每个弹簧振子的半功率带宽上，APSD必须近似均匀，且阻尼很小（ς≤0.1）。对于弱阻尼，半功率带宽近似为Br=2ςfn。在这种情况下，SDOF振子的响应将是窄带的。

将平稳信号的损伤方程式（8）中的伪速度替换为应力，并写成连续形式，得到式（9）。

式中：p(S)为应力最大值的概率密度函数（PDF）；T为暴露在应力环境中的总时间；vm+为应力历史中每单位时间正最大值的个数。对于窄带弱阻尼振子响应，最大值每 1fn秒出现一次，所以vm+可以替换为fn。

此外，在假设条件下，无论加速度环境的PDF如何，振子的响应都将接近高斯分布，峰值的PDF将接近瑞利形式[19-22]。方程（10）显示了峰值的瑞利分布。

式中：S为峰值的应力值；σS为应力时程的均方根。

将式（10）代入式（9）得到式（11）。

式中：Γ为伽马函数。回顾伪速度和应力之间的关系，公式（11）可以重写为：

利用Parseval定理可以计算出均方根伪速度振子响应[23]。首先，通过将伪速度FRF平方量级乘以APSD来计算频率响应量级的平方。然后在所有频率上对结果进行积分，以根据APSD量级Pa计算振子伪速度RMS、σpv。假设APSD环境在每个振子的半功率带宽内相对平坦，则允许使用闭合形式近似（类似于Mile′s方程），见式（13）。

均方根伪速度响应通过式（13）计算，每个振子损伤通过式（12）计算，最终形成疲劳损伤谱D(fn,Q,b)。计算过程中可能需要插值APSD。由于APSD通常由对数坐标给出，因此APSD谱中两点之间的插值将位于对数曲线图中的直线上。在频率fi（位于f1和f2之间，对应APSD量级A1和A2）处，插值APSD可以通过式（14）来完成，其中N为对数斜率。

3 时域、频域算法等效性验证

上文给出了疲劳损伤谱在时域和频域中的计算方法，为验证两种方法计算结果的等效性，构建加速度功率谱，如图1所示。APSD谱的具体参数见表1。由此通过逆傅里叶变换并指定随机相角可生成遍历平稳、高斯的随机时间历程，如图2所示，持续时间为10 s。此时间历程的功率谱与给定的加速度功率谱对比如图3所示，可以看出，二者在低频部分吻合较好，高频部分存在一定误差，但在±1.5 dB之内。

图2 等效的高斯平稳随机时间历程 Fig.2 Equivalent gaussian, random time-histories

图3 APSD谱对比 Fig.3 APSD Spectrum Comparison

表1 APSD谱参数 Tab.1 APSD Spectrum

分别由图1所示的APSD谱及时间历程计算疲劳损伤谱。在计算中取参数Q=10，b=8，持续时间为10 s。用两种方法计算的FDS如图4所示。

图1 APSD谱 Fig.1 APSD Spectrum

图4 疲劳损伤谱对比 Fig.4 Fatigue Damage Spectrum Comparison

可以看出，时域法在100 Hz以前的低频部分，计算出的损伤略大于频域法。在100 Hz以后，时域法与频域法结果几乎相同。

4 结语

上述研究表明，高斯随机振动在时域和频域计算的疲劳损伤谱是具有等效性的。因此，可以通过比较不同加速度时间历程的FDS来比较环境严酷程度。在满足随机性和平坦性假设的前提下，也可通过不同环境的APSD谱（包含持续时间）的FDS来比较环境严酷程度。