多种随机载荷下的风电设备动态可靠度分析

孙 逸，张西龙

（中国电建集团西北勘测设计研究院有限公司，陕西省西安市 710065）

0 引言

可靠度理论，一般是基于应力—强度干涉理论而建立的静态分析理论[1]。然而在实际工况下，设备受到的载荷往往是多次且随机的，并且设备在疲劳、老化等因素的影响下，其强度也必然会变化，因此，基于干涉理论的可靠性模型一般会造成不可忽视的误差[2-3]。国内学者王正[4]运用统计学原理，建立随机载荷多次作用下的等效分布函数，将随机载荷的作用次数视为随时间变化的量，由此提出了动态可靠度计算理论。

对于风电设备，由于其输入载荷、材料状态等均具有强烈的随机性，可靠性模型往往处于非线性受力状态。同济大学李杰教授[5]及其团队提出了一种适用于一般随机系统，可以分析系统随机信息的概率密度演化理论，在复杂多自由度结构系统分析方面取得了较为不错的反响。本文结合王正提出的动态可靠度模型与李杰的概率密度演化法，提出一种新的动态可靠度计算方法，旨在获得真实、有效且计算简便的结果。

1 动态可靠度的概述

一般情况下，风电设备可靠性的设计与核验多采用有限元建模方法，即对风电设备各个零件进行强度分析，继而扩展至整个机组系统。该方法基于商业软件的计算，是现在的主流分析手段[5]。但是，有限元建模方法往往只能对模型进行静力、较短时间范畴的动力可靠度计算。这与实际情况下有着显著的差异，风电机组往往承受着强风、地震等在空间、时间、强度维度上变化的随机载荷。因此，单一地考虑风电设备在某时间节点、受到某一定量的载荷而求得的静态可靠度，无法准确地计算长时间工作的风电设备寿命与可靠度，不具备普适性。

可靠度，作为工程领域中一个非常重要的指标，其对于产品设备的结构设计、维护修理、寿命预测有着不可代替的指导作用。因此，精确计算设备可靠性，从而保障人员生命安全、减少财产损失、节约生产成本的目的，是现代机械领域十分重要的研究方向。

可靠性，是定性指标，一般指设备在一定时间下，正常运转的能力。而可靠度，则是数量指标，指设备在一定时间下，正常运转的概率[6]。

假定某一设备零件所受载荷的应力S和自身的强度γ分别为服从相互独立的某种随机分布的概率函数，其概率密度函数分别为fS（S）、fγ（γ）。则可靠度可定义为：

如图1所示为应力—强度干涉理论的示意图。在同一坐标系下，当载荷S与零件强度γ分别服从某种分布时，二者的概率密度未重叠部分属于绝对安全区，因为该区域内γ＞S；载荷S与零件强度γ重叠部分，即图中阴影为失效区，γ＜S，称之为应力S和强度γ的干涉区。显然，当干涉区越小，意味着可靠度越大，零件越安全。

图1 应力S和强度γ的干涉Figure 1 Interference between S and γ

由积分学可得，事件（γ＞S）发生的概率，即可靠度为：

实际工作中，机械设备往往是在服役期内多次往复运作，若忽视应力作用次数及其带来的疲劳磨损，将会导致模型计算与实际工程严重脱节。

由应力—强度干涉理论可得，载荷多次作用下的零件可靠度可通过如下积分式获得：

设在任意时刻t，载荷出现n次的概率服从参数为λ的分布：

其中，N（t）指在时刻t时载荷出现n次；P（n）为载荷出现n次的概率。

结合可靠度公式（3），零件在t时刻载荷作用n次时的可靠度为：

通过极限求和，可以最终得到零件在任意时刻下的动态可靠度：

进一步对公式（6）中的指数部分进行泰勒展开，可以将其化简为：

以上公式都是建立在零件自身的强度不变的基础上进行。但在实际工程中，机械零件自身由于疲劳损伤、腐蚀老化等原因，其强度会发生退化。本文将零件强度的退化简化为与时间有关的函数，对原强度分布进行加权处理。零件的强度取决于初始强度γ0与加权因子的共同影响。其公式如下：

2 概率密度演化理论

上述动态可靠度模型，将静态可靠度改良为载荷多次作用下的时变动态可靠度。但是，该模型仅仅将输入载荷简单地归纳为一个或多个同种的载荷，再进行简单地叠加、相乘。事实上，风电机组的工况往往是复杂、多变的：无论是输入载荷、零件自身状况，其在时间、空间上的激励均具有随机性，导致结果是系统将以多个维度且非线性的形式输出。若运用干涉模型进行积分，将会带来巨大的计算量。而李杰教授[7]提出的概率密度演化法则很好地解决了这个问题，该理论的诞生就是为了解决土木工程中，强风、地震、海啸等强随机性输入下的非线性动力系统问题。

记Z（t）为需要研究解决的最终结果，系统所有的随机参数用同一函数表征，则有等式如下：

表达为积分形式，则为：

在微小时间增量dt之后的t+dt时刻，Ωt+dt是Ωt及其边界运动叠加的结果，即：

由此可见，无论在任意时刻t≠t0，Ωt均依赖于Ωθ。这也是必须考虑增广系统而非系统Z(t)的原因。

将式（12）代入式（11），先考虑等式的右侧，有：

将式（13）代入式（11）右侧，化简可得：

应用散度定理，并略去高阶项，式（14）可记作：

考虑到Ωt×Ωθ的任意性，并消去等式两边的dt，可最终得到概率密度演化方程：

如图2所示为不同时刻下典型概率密度演化函数。可以清楚地看到，随着时间的演化，可靠度关于载荷强度的概率密度有着明显的变化。随着零件工作时长的增加，可靠度的概率密度会逐渐“扁平化”，说明该零件在更广泛的载荷范围内都有着失效的风险，即随着零件工作时间的增加，其失效的可能性越大，这与我们的常识相符合。

图2 不同时刻下典型概率密度演化函数Figure 2 Evolution function of typical probability density at different time

3 概率密度演化理论下的动态可靠度

3.1 可靠度的数学模型

基于应力强度干涉理论的可靠度概率密度函数表达式为：

由于结构在服役过程中，不再有其他随机载荷介入，此时整个系统处于概率守恒状态，故pR(t)应满足概率密度演方程：

其中，a(t)为Dirac函数，它在除零以外的点上都等于零，且其在整个定义域上的积分等于1。

联立式（18）与上述初始条件，可以求得动态可靠度的密度函数p(r，t)。而通过对该密度函数关于载荷定义域上的不定积分，可得最终的动态可靠度：

式（19）即为基于概率密度演化理论的零件动态可靠度的理论解析表达式。由方程（18）可知，该联立方程的求解可以视为于解偏微分方程的初值问题。

若考虑零件强度随时间退化，可将零件强度进行加权处理。一般零件强度退化模式有：线性退化、对数退化、指数退化。结合式（17）与强度退化函数，可得在退化条件下，多次随机载荷下的可靠度概率密度函数表达式：

同式（18），其概率密度演方程为：

最终的动态可靠度Rdec(t)的表达式为：

3.2 演化方程的数值计算

简单的偏微分方程，一般的商业数学软件就能满足计算需要，但在面对存在巨大计算量的实际工程需要时，数值求解是现实而又必须的[8]。本文拟用数值分析中最常用的单边差分格式有限差分法，来解决该问题。概率密度演化理论的数值计算可采用如图3所示的流程。

图3 计算流程图Figure 3 Calculation flow chart

如图4所示，对求解域进行网络划分，即将r-t空间进行离散：rn= nΔr，ti=iΔt。u在点（n，i）的值为，根据微分的数学定义，可以近似得到的值为。同理，。因此，偏微分方程在点（nΔr，iΔt）处可近似为：

图4 网络划分Figure 4 Network division

3.3 算例

某一零件强度服从均值为μr=300MPa、标准差为σr=30MPa的正态分布；载荷S服从均值为μs=200MPa、标准差为σs=40MPa的正态分布。在载荷多次作用时，取λ=1，va(t)=0.5，误差不超过10-3。

根据式（22），求得的动态可靠度密度函数如图5所示。

图5 随时间演化的概率密度函数Figure 5 Probability density function evolving with time

可以清晰看到，经由概率密度演化法计算的动态可靠度其概率密度会随着时间的演化而扁平化，其函数图峰值之间变低，在载荷强度定义域上逐渐分散。纵观整个演化过程，其形如缓缓流动的河流；在三维视角上，则如连绵不断又高低起伏的山脉。该现象为概率守恒的系统在状态空间中流动演化的结果。回归到零件可靠度的问题是，可以解读为随着零件工作的运作，可靠度的密度函数在载荷域上的分布逐渐变宽，致使零件失效的载荷域渐渐变大，零件失效的可能性在不断提高。

如图6、表1所示，为利用概率密度演化法计算得到的动态可靠度，与前文的动态可靠度和蒙特卡罗法模拟10000次的对比。可见，用本文方法求得的可靠度，仍保留了失效规律的“盆浴曲线”，并且比应力—强度干涉理论的结果更好地拟合了蒙特卡洛模拟结果[9]。由图可知，构件磨合期结束后，零件的主要失效形式由磨损失效变为疲劳失效，构件的失效率会显著降低，可靠度下降趋势则逐渐变缓。当构件工作时长超过其额定工作时长后，零件的主要失效形式又由疲劳失效变为磨损失效，失效率显著上升，可靠度下降趋势随之增大。利用概率密度演化法求得的动态可靠度，在构件进入疲劳失效期后，存在“磨损—疲劳”状态改变的阶变时间点，可靠度下降趋势变化更加明显，说明本文方法在描述构件动态可靠度的状态改变方面更为准确。分别利用PDEM和干涉理论计算出的结果，与蒙特卡洛法做对比散点图，如图7所示。显然，相比较干涉理论，PDEM计算的动态可靠度可以更好地拟合蒙特卡洛模拟结果，这足以说明密度演化模型的准确性。

图6 零件的动态可靠度Figure 6 Dynamic reliability of parts

图7 散点图Figure 7 Scatter chart

表1 不同时刻下单一零件的动态可靠度Table 1 Dynamic reliability of a single part at different times

若考虑零件为强度随时间退化，其退化因子为μ(t)=1-0.00004t的线性函数，则图8为退化条件下单一零件的动态可靠度，其中红色为强度退化，黑色为强度不退化的对照组。表2为两种条件下可靠度的数值记录。

图8 退化条件下单一零件的动态可靠度Figure 8 Dynamic reliability of single part under degradation condition

表2 不同时刻强度退化零件的动态可靠度Table 2 Dynamic reliability of parts with strength degradation at different times

4 结论

风电设备的可靠度计算，多采用有限元建模的方式。该方法较为准确、直观，但由于有限元方法只能计算短时间内的结果，无法预测其长时间的状态或可靠度。这显然对于价值高、运转时间长的风电设备是不利的。本文结合应力—强度干涉理论与概率密度演化方法，提出了一种可以考虑多种随机输入载荷，预测设备实时可靠度的数学模型。通过与仿真数据与Monte Carlo模拟数据对比，该方法可以有效预测设备的寿命及可靠性，可以为有限元建模与工程实际提供必要的辅助作用。