机械结构可靠性分析的主要目的是在考虑结构各输入参数随机性情况下计算其失效概率.近年来国内外学者将一阶矩、二阶矩、Monte Carlo(MC)^[1]、代理模型等方法应用于机械结构可靠性分析.一阶矩及二阶矩方法难以适应工程结构非线性隐式的特点, 精度有限.MC方法鲁棒性最好, 然而耗时过长.代理模型方法是目前应用最广泛的机械结构可靠性分析方法.

众多学者通过调整代理模型函数形式^[2-6]及选取样本点策略^[7-10]提高代理模型精度及效率.Blatman^[2], Roussouly^[3]等将稀疏多项式应用于结构可靠性分析, 通过不同方法筛选“重要项”.此外, 为提高模型的非线性程度, Kriging^[4-5]、神经网络^[6]、SVM^[7]也被广泛采用.Kriging模型是一种高效的差值模型.因此, 本文采用Kriging作为基本代理模型.

样本点的选取策略对代理模型精度和效率有很大影响.为提高样本点的质量, 很多自适应抽样方法被提出^{[3, 5, 7-10]}.现有抽样方法中鲜有考虑在迭代抽样过程中可能出现由于个别样本点间距离过近甚至矩阵奇异问题, 尤其是当算法接近收敛时, 设计点附近很小范围内往往集中过多样本点.本文提出一种新的自适应抽样方法，避免对“不重要区域”没必要抽样; 保证抽取的样本点保持与现有样本点的距离.

假设某机械结构的输入变量为X(X为N维随机向量), X的联合概率密度函数为f(x), 该结构的功能函数为G(x), 其安全域和失效域分别为S_s={x|G(x)>0, x∈R^N}和S_f={x|G(x)≤0, x∈R^N}, 则该结构失效概率为

(1)

若已知样本集Ω={(x_i, y_i), i=1, 2, …, M}, 且功能函数G(x)可表示为

(2)

其中：g(x)为多项式基函数；β为g(x)的系数向量, 本文中取g(x)的次数为1；z(x)为零均值同方差高斯过程, z(x_i), z(x_j)的协方差为

(3)

其中：σ_z²是高斯过程的方差；R(x_i, x_j; θ)表示z(x_i)和z(x_j)的相关系数, θ为参数.应用最广泛的是高斯相关函数, 其形式为

(4)

其中, x_{i, n}表示向量x_i的第n个元素.

在给定样本集Ω情况下, 式(2)~式(4) 中的未知参数β, σ_z², θ可通过极大似然法估计得到,

(5)

式中:为广义最小二乘法估计的残差；R为M×M相关系数矩阵；分别表示θ, β, σ_z²的极大似然估计值.点x处结构的功能函数可表示为Y的线性组合形式,

(6)

式中，Y=[y₁, y₂, …, y_M]^T.

要实现的最小方差无偏性^[11], 则c(x)需满足

(7)

求解式(7), 可得的Kriging表达式为

(8)

则该结构失效概率P_f的估计值为

(9)

当X的维数较高时, 直接求式(9) 的积分值很困难.本文采用MC法通过随机抽样近似计算, 根据大数定律,

(10)

其中：N_MC为随机抽样次数；x_{MC, i}(i=1, 2, …, N_MC)是来自f(x)的独立同分布随机序列.由中心极限定理知, 当N_MC足够大时, 的均值依分布收敛于正态分布, 可得

其中:Ф(·)为标准正态分布函数; λ₀为正数.

若要保证在95%的置信水平下与间的相对误差不超过5%, N_MC需满足

(11)

通常情况下工程机械结构的失效概率数量级很小, (1－)应接近1, 因此忽略分母中的(1－), 式(11) 可写成

(12)

N_MC表示当结构的极限状态函数为=0时N_MC次随机抽样中失效样本个数的期望.

本文提出在=0上选取具有“代表性的”点, 计算这些点的结构功能函数值, 更新=0, 在迭代过程中使=0逐渐接近G(x)=0，直至满足收敛条件.

本文提出的选取样本点方法的主要步骤为

步骤1  t=0, 应用拉丁超立方抽样方法^[12]随机产生初始样本点并计算各点结构功能函数值, 计算初始样本点个数为M₀, 令

步骤2  t=t+1, 随机抽取K个=0上的点.给定, 生成N维服从f(x)的随机向量, 若该随机向量满足式(13), 则认为其在=0上.当满足式(13) 随机向量数达到K个时, 停止生成随机向量过程, 并记

(13)

步骤3  对作聚类分析, 将其分成k类并得到k个类别的中心, 并将k个中心点映射到=0上.令{s_{t－1, 1}, s_{t－1, 2}, …, s_{t－1, k}}表示k个聚类中心.通常情况下=0是非线性曲面, 因此不能保证k个聚类中心点在=0上, 因此要将其映射到=0上.映射的方法是在=0上寻找满足式(14) 的点, 得到, 其中为的设计点.

(14)

其中, i=1, 2, …, k.

步骤4  调整集合S_t-1中各点位置.如式(15) 定义D₀, 本文中认为两个样本点间距离小于D₀时是不可接受的, 需要调整个别点位置.

(15)

式中，e₀为给定常数.

步骤5  计算S_t-1中各样本点的结构功能函数值.令

步骤6  根据Ω_t计算, .若满足式(16) 的收敛条件, 则停止迭代过程, 即为P_f估计值; 否则返回步骤2, 直至满足收敛条件.

当迭代过程收敛时, 可以认为

由中心极限定理, 若能保证抽样次数满足

则

根据式(11) 和式(12),

收敛条件为

(16)

选用文献[10]中某二维输入变量结构作为分析对象, 该算例的主要目的是进一步用于说明本文所提出的方法.该结构的功能函数为

其中, x₁, x₂相互独立且均服从标准正态分布.

(17)

其中,

通过MATLAB软件, 采用数值积分的方法计算式(17), 得到失效概率为

采用拉丁超立方抽样方法在[-5, 5]×[-5, 5]内随机产生6个初始样本点(M₀=6), 计算其功能函数值, 得到Ω₀, 建立初始Kriging代理模型, 如图 1a所示.取K=200, ε=0.01, 随机产生K个=0上的点.取e₀=0.15, k=2, 依据2.1节步骤3的方法, 对K个随机点进行聚类分析, 并将k个聚类中心映射到=0上.图 1为算法迭代过程中G(x)=0与=0的对比及每步迭代产生的新样本点.算法在第4次迭代时收敛, 计算得到=3.59×10^－3, 图 2为收敛后G(x)=0与=0的对比.

图 1(Fig. 1)

图 1 迭代过程中G(x)=0与=0对比 Fig.1 Comparison between G(x)=0 and =0 during iteration (a)—第1次迭代；(b)—第2次迭代; (c)—第3次迭代；(d)—第4次迭代.

图 2(Fig. 2)

图 2 收敛后G(x)=0与=0对比 Fig.2 Comparison between G(x)=0 and =0 after convergence

本节选用文献[2-3]中的桁架结构作为分析对象, 如图 3所示.该结构有10个输入随机变量, 分别为A₁, A₂, E₁, E₂, P₁, …, P₆.其中A₁, A₂分别为水平杆及倾斜杆的横截面积；E₁, E₂为水平杆及倾斜杆的杨氏模量；P₁, …, P₆为随机载荷.涉及到的10个随机变量的分布类型及参数如表 1所示.

图 3(Fig. 3)

图 3 桁架结构示意图 Fig.3 Figure of the truss structure

表 1(Table 1)

表 1 实例2随机变量分布类型及参数 Table 1 Distribution patterns and parameters for 2nd example 变量 分布类型 均值 标准差 P1~P6/N Gumbel 5×104 7.5×103 A1/m2 对数正态 2×10－3 2×10－4 A2/m2 对数正态 1×10－3 1×10－4 E1/Pa 对数正态 2.1×1011 2.1×1010 E2/Pa 对数正态 2.1×1011 2.1×1010

表 1 实例2随机变量分布类型及参数 Table 1 Distribution patterns and parameters for 2nd example

该结构在P₁, …, P₆随机载荷作用下点E在竖直方向位移最大, 因此以E点位移(s(x))作为结构响应, 为与文献[2]保持一致, 位移阀值设为0.14 m, 因此该结构功能函数为

此桁架结构的失效概率约为3.45×10^-5, 该结果来自500 000次重要抽样^[2].

令K=3 000, ε=0.001, k=8, e=0.15, 首先将10个随机变量变换为标准正态随机变量, 再应用第2节所提出方法计算该结构的失效概率.拉丁超立方抽样在[-5, 5]¹⁰立方体内进行(M₀=15).图 4为应用本文所提方法运行5次所得的迭代次数与失效概率预测值折线图, 由于失效概率预测值数量级变化比较大, 纵坐标采用对数形式.表 2为本文方法与其他方法精度对比.

图 4(Fig. 4)

图 4 Pf收敛折线图 Fig.4 Broken lines of Pf during convergence

表 (Table )

表 实例2结果对比 Table Comparison of results for 2nd example 来源 ×105 相对误差/% Mcall 参考值[2] 3.45 — — 文献[2]方法 2.35 -31.9 207 文献[3]方法 3.30 -4.3 142 本文方法 3.36~3.64 -2.6~5.5 81~89

表 实例2结果对比 Table Comparison of results for 2nd example

实例1主要用来具体说明本文所提出的结构可靠性分析方法.图 1表明, 本文所提出方法能够保证选取的样本点落在“重要区域”, 避免了对非重要区域的不必要抽样.同时, 各样本点间能够保持一定距离, 这也在很大程度上提高了算法的效率.实例2选用的桁架结构在其他文献[2-3]中已有提到, 通过对比可知, 本文所提出方法能够在较少调用结构功能函数情况下得到较高精度失效概率估计值.通过图 4可知, 本文方法收敛速度快且稳定.

1) 提出了一种自适应的选点方法, 通过实例表明该方法能够保证所选样本点落在对失效概率贡献较大的重要区域, 有效减少了对非重要区域的抽样次数, 提高各样本点的质量.

2) 提出的自适应方法能够保证各样本点间距离不会太小, 结合Kriging模型, 能够更充分利用各点统计信息, 提高实效概率计算精度.

3) 提出的机械结构可靠性分析方法能够在较少调用结构功能函数情况下实现对单设计点机械结构失效概率的较准确估计.