机械产品的动态可靠性是指其在给定的时间区间上能够正常运行的概率^[1]。输入的不确定性必然导致实际输出具有不确定性, 而输出的不确定性在一定程度上便会导致机械产品的破坏。可靠性分析能够得到输出响应量的统计规律, 从而能够衡量机械产品安全运行的概率, 为保证产品的安全提供参考。

目前, 已经形成的比较有效的求解机械产品可靠性方法包括一次可靠性方法(FORM)^[2-3]、改进的一次可靠性方法(FOSM)^[4-5]和蒙特卡罗方法(MCS)等^[6]。蒙特卡罗法虽然对功能函数的形式及变量的概率密度函数和维数没有限制且一般精度较高, 但其计算量巨大。改进的一次方法能够很好地权衡计算精度和计算效率, 本文在该方法基础上研究时空动态可靠性及其全局灵敏度分析的高效方法。对于工程中常见的输入变量变异性较小的问题, 该类方法具有较高的精度。

对于动态可靠性问题, 目前使用较多的方法可以分为2类:基于极值响应的方法^[7-9]和基于首次超越的方法^[10-11]。基于极值响应的方法针对极值分布可求得的情况, 其能够将动态可靠性问题转变为静态可靠性问题进行求解。目前基于首次超越的思想方法中使用较多的方法是使用Rice公式^[12], 但获得跨越率依然很难。

包络函数法常应用于电子工程领域^[13], 有学者^[14-15]将包络函数法应用于结构优化问题。Du^[16]提出了动态可靠性分析问题的包络函数法, 能够很好地克服以上动态可靠性方法的缺点。该方法通过引入失效边界上失效概率最大点处的包络面来近似失效域, 从而将动态可靠性问题转化为扩展点处的静态可靠性问题。由于该方法只针对包含一维参数的动态可靠性问题, 因此并不能用来分析多参数动态可靠性问题。

可靠性全局灵敏度分析可以得到影响产品可靠性的输入变量的重要性, 如静态问题的可靠性全局灵敏度分析^[17]。为了能够分析得到影响机械产品动态可靠性的基本输入变量的重要性, 十分有必要研究动态可靠性的全局灵敏度分析方法。

本文在文献[16]基础上提出一种计算含有时间和空间多维参数的动态可靠性问题的包络函数法, 该方法能够很好地克服以上方法的缺陷, 得到较为准确的结果。对于多输出的问题, 该方法亦具有一定的适用性。同时, 在该可靠性方法的基础上, 本文结合三点估计提出了一种计算可靠性全局灵敏度指标的新方法。

1 时空动态可靠性分析定义

时间t和空间z的多参数的动态输出响应为f(x, z, t)。其中, x=(x₁, …, x_n)是n维随机输入变量, z∈[z₀, z_f]是空间位置参数, t∈[t₀, t_f]是时间参数。为了保证机械产品能够正常工作, 要求输出响应必须不大于要求的阀值ε。因此, 多参数动态可靠性定义为在给定的时间和空间内输出响应满足要求阀值的概率, 定义式如下

(1)

而失效概率则定义为至少存在一个时刻或者空间位置使输出响应大于阀值的概率, 即

(2)

式中, 并集(z∈[z₀, z_f])∪(t∈[t₀, t_f])表示时间与空间组成的参数取值域, 若以z为横坐标, t为纵坐标, 则其表达为图 1中的矩形域。

2 时空动态可靠性的包络函数法 2.1 包络函数的定义

分别使用F₁(x)=0和F₂(x)=0作为随时间t和空间z连续变化的动态失效边界f(x, z, t)=ε和f(x, z, t)=-ε上的包络函数, 首先考虑正的失效边界f(x, z, t)=ε的包络函数F₁(x)=0。f(x, z, t)=ε可以看作是当参数z和t在其区间变化时的一系列随机输入变量空间上的超曲面。

考虑空间和时间的2个邻近超曲面分别为f(x, z, t)=ε与f(x, z+δz, t)=ε及f(x, z, t)=ε与f(x, z, t+δt)=ε, 其中δz>0, δt>0, 则有

(3)

(4)

当δz和δt无限趋近于零时, 可以得到如下的偏微分方程

(5)

(6)

上述两式可以简单记为如下形式

(7)

(8)

理论上, 通过求解以上2个偏微分方程可以得到空间参数z和时间参数t关于输入变量x的函数, 记(7) 式和(8) 式的解为z=h₁(x)和t=h₂(x), 并将z和t的解代入f(x, z, t)-ε=0便可以得到正失效边界包络函数F₁(x)的表达式如下

(9)

同理可以求得负失效边界的包络函数表达式F₂(x), 此时动态可靠性问题转化成如下的形式

(10)

由上面的公式可以看出, 动态可靠性的求解已经被转化成一个静态可靠性的问题。以下将结合一次可靠性的近似方法来估计包络函数。

2.2 包络函数的估计

假设输入变量相互独立且均服从正态分布, 对于输入变量非正态或相关的情况可以使用R-F方法与Rosenblatt方法[6]将其转化成正态变量。在输入变量变异性较小的基础上使用一阶泰勒公式将输出响应函数在输入变量均值处进行线性展开^[6]

(11)

式中, μ=(μ₁, …, μ_n)^T表示输入变量x的均值向量, a₀(z, t)和a(z, t)分别由下式决定

(12)

(13)

为了简化运算, 将随机输入变量x转化为标准正态变量u

(14)

式中，σ=(σ₁, …, σ_n)^T为x的标准差列阵。则(11) 式可以改写成如下形式

(15)

式中, b₀(z, t)和b(z, t)分别为

(16)

(17)

到此, 就把原求解包络函数的问题转化为分别求解线性函数族L(u, z, t)=ε和L(u, z, t)=-ε的包络函数F₁(u)=0和F₂(u)=0的问题。

首先考虑正失效边界的包络函数F₁(u)=0, 结合(7) 式和(8) 式, 其可以由下式决定

(18)

式中, 分别表示b₀(z, t)和b(z, t)对参数z求偏导数, 分别表示b₀(z, t)和b(z, t)对参数t求偏导数。

由包络函数的意义可知, 包络函数应该在参数区间[z₀, z_f]和[t₀, t_f]上包含所有的失效域。则包络函数F₁(u)=0上的扩展点应该在线性函数L(u, z, t)=ε上, 并且在失效边界上具有最大的概率密度。假定u^*为包络函数上的扩展点, 则u^*为线性函数L(u, z, t)=ε上离原点最近的点。因此, u^*向量与线性函数L(u, z, t)=ε相互垂直, 即u^*与线性函数L(u, z, t)=ε在u^*处的梯度共线。如图 2所示在标准正态空间下, 小方点即为扩展点。

据该性质可得以下等式

(19)

式中, c为常数, 表示线性函数L(u, y, t)=ε梯度方向上的单位向量。(19) 式满足(18) 式, 将(19) 式代入(18) 式中的第一个等式有

(20)

通过以上等式可以解得常数c的值如下

(21)

将常数c的值代入(19) 式可得扩展点的值为

(22)

将求得的扩展点代入(18) 式的第二和第三个等式得到关于参数z和t的方程组如下

(23)

通过求解以上方程组, 便可以得到参数z和t的一组解。由于此时是对应正失效边界包络函数进行求解, 因此其输出响应必为正值。将求得的解代入(12) 式计算均值输出, 排除掉对应均值输出为负值的解, 并且将最后得到的参数值记为(z_i⁽¹⁾, t_i⁽¹⁾)(i=1, 2, …, m₁), 则扩展点为

(24)

同理, 对应负失效边界的包络函数F₂(u)=0, 可求得其对应的参数值解为(z_j⁽²⁾, t_j⁽²⁾)(j=1, 2, …, m₂), 其扩展点为

(25)

此时, 正失效边界包络函数F₁(u)=0就可以使用超平面L(u, z_i⁽¹⁾, t_i⁽¹⁾)=ε来估计, 其中i=1, 2, …, m₁; 负失效边界包络函数F₂(u)=0可以使用超平面L(u, z_j⁽²⁾, t_j⁽²⁾)=-ε来估计, 其中, j=1, 2, …, m₂。结合(10) 式可知, 求解可靠性此时转化为如下形式

(26)

上述包络面考虑了空间和时间参数分别在给定区间[z₀, z_f]与[t₀, t_f]内变化的情况, 其空间与时间边界以及端点处的情况也需要考虑。考虑响应在时间和空间参数区间4个边界上包络函数的情况, 即在边界(z∈(z₀, z_f), t₀), (z∈(z₀, z_f), t_f), (z₀, t∈(t₀, t_f))和(z_f, t∈(t₀, t_f))上求解其正负失效边界包络函数对应的参数解。此时, 问题相当于求解单参数动态可靠性的问题, 详细过程可参考文[16]。使用文[16]中单参数包络函数法求得4个边界上的包络函数F₁(u)=0的m′₁个参数解(z_i⁽¹⁾, t_i⁽¹⁾)(i=m₁+1, m₁+2, …, m₁+m′₁); 以及包络函数F₂(u)=0的m′₂个参数解为(z_j⁽²⁾, t_j⁽²⁾)(j=m₂+1, m₂+2, …, m₂+m′₂)。最后还需要引入参数区间端点的包络函数, 考虑响应在参数区间端点处的情况, 则需要分别加上(z₀, t₀), (z₀, t_f), (z_f, t₀)和(z_f, t_f)处的失效边界L(u, z^*, t^*)(其中z^*为z₀或z_f, t^*为t₀或t_f)。

令(z₁^*, t₁^*)=(z₀, t₀), (z₂^*, t₂^*)=(z₀, t_f), (z₃^*, t₃^*)=(z_f, t₀), (z₄^*, t₄^*)=(z_f, t_f), 则考虑所有端点及边界后原时空动态可靠性问题可简化为

(27)

由于(27) 式的每一个极限状态方程均是线性正态的, 因此可以先求解所有等价输出的联合概率密度函数, 然后再由数值积分法求得可靠度。此处使用多元正态分布函数Φ_m(ε, μ, Σ)来求解(27) 式中的可靠性问题, 其中μ和Σ分别为均值向量和协方差向量。引入如下指示函数

(28)

式中, i=1, 2, …, m₁+m′₁, j=1, 2, …, m₂+m。考虑正失效边界包络面上发生负向失效(即其值小于-ε)与负失效边界包络面上发生正向失效(即其值大于ε)的可能性较小, 可将(27) 式写成如下形式。

(29)

式中, m=m₁+m′₁+m₂+m′₂+4。则可通过如下积分求解该动态可靠性问题。

(30)

式中, 。

2.3 包络函数求解时空动态可靠性的基本步骤

综合上述包络函数求解时空动态可靠性的过程, 其主要步骤可以给出如下:

步骤1 通过(23) 式求得正失效边界包络函数对应的参数值解, 并代入(12) 式计算均值输出, 剔除掉对应均值输出为负值的解, 最后得到参数值解(z_i⁽¹⁾, t_i⁽¹⁾)(i=1, 2, …, m₁); 同理计算出负失效边界包络函数对应的参数值解(z_j⁽²⁾, t_j⁽²⁾)(j=1, 2, …, m₂)。

步骤2 计算时间和空间参数边界以及端点对应的包络函数。

步骤3 采用数值积分或是数字模拟等方法求解时空动态可靠性问题的近似解, 此处采用(30) 式来求解时空多参数动态可靠性。

2.4 计算多输出时空动态可靠性

考虑如下所示的多参数多输出串联响应问题

(31)

式中, l为输出响应的个数。分别对上式左边输出响应函数f_k(x, z, t)(k=1, 2, …, l)在输入变量均值处进行一阶泰勒展开并对输入变量标准正态化

(32)

式中, b_0k(z, t)和b_k(z, t)(k=1, 2, …, l)的定义如(16) 式和(17) 式所示, u的定义见(14) 式。

使用包络函数法求解多输出问题的步骤与以上求解单输出问题的步骤类似。使用(18)~(25) 式分别计算输出响应f_k(x, z, t)(k=1, 2, …, l)在失效边界上对应的参数解。以m₁^(k)+m′₁^(k)和m₂^(k)+m′₂^(k)分别表示第k(k=1, 2, …, l)个输出L_k(u, z, t)≤ε和L_k(u, z, t)≥-ε所对应的参数解(z_i^(k1), t_i^(k1))(i=1, 2, …, m₁^(k)+m′₁^(k))和(z_j^(k2), t_j^(k2)) (j=1, 2, …, m₂^(k)+m′₂^(k))的个数, 且令端点的参数解为(z₁^(k)*, t₁^(k)*)=(z₀, t₀), (z₂^(k)*, t₂^(k)*)=(z₀, t_f), (z₃^(k)*, t₃^(k)*)=(z_f, t₀)和(z₄^(k)*, t₄^(k)*)=(z_f, t_f) (k=1, 2, …, l)。由于是串联系统, 则其总的可靠性为每个输出响应均安全时对应的可靠性, 即多输出的可靠性可以表示成如下形式

(33)

由(33) 式可以看出, 多输出时空动态可靠性问题此时已经被转化为静态的可靠性问题, 其值可采用显式极限状态函数的可靠性求解的数值积分或数字模拟来求得。

3 动态可靠性全局灵敏度分析的包络函数法 3.1 方法思路

对于一个动态可靠性的问题, 其失效域定义为

(34)

由文献[17]可知, 输入变量对失效概率P_f影响的可靠性全局灵敏度因子可以定义为

(35)

式中

(36)

(37)

式中, P_f表示机械产品的失效概率, P_f|Xi表示X_i取其实现值时的失效概率, E(·)和V(·)分别为期望算子和方差算子, I_ΩF为失效域指示函数, 当x∈{Ω_F}时, I_ΩF(x)=1, 否则I_ΩF(x)=0。且对(36) 式经过进一步证明可得^[18]

(38)

由上式可以看出, δ_i^P的求解包含两部分, 即内层条件失效概率和外层方差的求解。

本文将在统计矩方面具有较强优越性的三点估计法^[19]和包络函数法结合起来计算可靠性全局灵敏度指标。根据三点估计法, 可以得到动态失效概率P_f和δ_i^P如下^{[6, 20]}：

(39)

(40)

式中, l_Xi^(k)和P_Xi^(k)(k=1, 2, 3) 分别为三点离散分布的取值点和相应的概率值^[20], (k=1, 2, 3) 表示把变量X_i固定在名义值(k=1, 2, 3) 时的动态条件失效概率。

3.2 基本步骤

步骤1 通过文献[20]中公式分别求得三点离散分布的取值点和相应的概率值, 即l_Xi^(k)和P_Xi^(k)(k=1, 2, 3)。

步骤2 分别将变量X_i固定在名义值处, 即令X_i^(k)=l_Xi^(k), 然后使用包络函数法计算其条件失效概率值(k=1, 2, 3)。

步骤3 将步骤2中得到的条件失效概率值代入(39)~(40) 式, 求得P_f和δ_i^P的值。

步骤4 将步骤3中求得的P_f和δ_i^P的值代入(35) 式, 便可以得到可靠性全局灵敏度指标S_i的值。

4 算例 4.1 悬臂梁

如图 3所示为一长方体悬臂梁结构, 其横截面宽为b, 高为h。受力形式如图所示, 其中, 不规则分布载荷, 悬臂梁末端受到随时间t变化的力F=F₀cost, 固定弯矩M=12 167 N·m。悬臂梁截面上下表面为其正应力最大处, 且在该载荷作用下上下表面应力对称分布。构造上表面应力σ_max为输出响应函数, 则有:σ_max=, 该结构失效为最大应力大于给定的临界阀值σ^*。在本算例中变量b、h、q₀和F₀为随机输入变量, 其分布形式如表 1所示; z和t分别为位置和时间参数, 且z∈[0,9], t∈[0,9]。

图 4给出了随机输入变量取均值时输出响应σ_max随时间和空间位置变化的曲面图, 从图中可以看出, 点(5.22, 0)、(5.22, 6.28) 和(8.18, 3.14) 为输入变量取均值时输出响应的峰值。表 2列出了不同阀值下的失效概率计算结果, 通过与蒙特卡罗方法的结果进行比较, 可以看出该方法计算的准确性。表 3还列出了包络函数法和蒙特卡罗法的计算量, 通过对比可以看出包络函数方法能够很大程度地提高计算效率, 减少计算成本。表 4中的输入变量全局灵敏度计算结果表明, 通过包络函数法得到的输入变量全局灵敏度指标的大小排序与蒙特卡罗法结果相同, 即:h>b>F₀>q₀, 因此, 通过减少输入变量h的不确定性能够最大程度改善结构的安全性能。

4.2 多输出问题

对于如下所示的串联多输出系统:, 该系统的输出响应为f(x, z, t)=max{|f₁(x, z, t)|, |f₂(x, z, t)|}, 系统失效定义为输出响应超过临界阀值a*。其中, x₁, x₂和x₃是相互独立的正态输入变量, 其分布为x_i~N(10, 0.1²)(i=1, 2, 3);z和t是位置和时间参数, 且z∈[1,9], t∈[1,9]。在此选取临界阀值a*=130, 使用包络函数法计算的系统失效概率和输入变量全局灵敏度指标如下。

由以上的计算结果可以看出, 所提方法对于多输出问题仍然适用。从表 5可以看出包络函数法计算出的多输出动态可靠性问题的失效概率与蒙特卡罗方法基本相同, 但是却能大大地减少计算量。从表 6可以看出, 使用包络函数法得出的输入变量全局灵敏度指标排序与蒙特卡罗法得出的排序结果相同, 即:x₁>x>x₃, 因此, 与其它输入变量相比, x₁的随机性对失效概率的影响最大, 可以通过减小x₁的不确定性来提高系统的安全性。

5 结论

本文提出了一种计算时间和位置多参数的动态可靠性分析的包络函数方法, 该方法通过使用多个展开点处的分段超平面来近似原随时间和空间连续变化的失效超曲面, 计算精度高而计算量较小。同时, 该方法亦适用于多输出系统。算例结果显示, 所提方法能够将计算量减小4个量级左右, 因此所提方法具有高效性。

鉴于基于失效概率的可靠性全局灵敏度指标在结构安全及可靠性评估领域中的重要作用, 本文提出了基于包络函数法和三点估计法的时空动态可靠性全局灵敏度指标的高效算法。所提方法能够在保证计算精度的同时, 很大程度上减小了计算量, 从而高效地找出影响结构动态安全的重要性变量并为后续的结构优化设计提供指导。

本文方法计算精度与输入变量的变异性及功能函数的非线性程度有关。对于输入变量变异性较大或功能函数的非线性程度较大的情况, 使用该方法可能会得到较大的误差。不过针对该问题, 可以采取一定的手段来提高它的精度, 比如使用更多的展开点等, 这将在后续的工作中进行研究。