基于二阶抛物线近似的结构可靠性分析方法

2024-03-09陈振中黄冬宇李晓科吴子豪

工程设计学报 2024年1期

陈振中，黄冬宇，田娇，李晓科，吴子豪

（1. 东华大学机械工程学院，上海 201620；2. 郑州轻工业大学机电工程学院，河南郑州 450002；3. 上海第二工业大学智能制造与控制工程学院，上海 200135）

传统的结构设计通常采用安全系数来保证产品达到预期的可靠度[1]。但是，安全系数过小会使结构可靠性过低，导致产品使用寿命不足；安全系数过大会使结构可靠性过高，导致产品生产成本增加。因此，为保证结构设计的合理性，开展可靠性分析是十分有必要的。

结构可靠性分析是指将结构尺寸、材料性能等各种不确定性因素视作随机变量[2-3]，并根据随机变量的分布来求解结构的可靠概率或失效概率。常用的可靠性分析方法主要分为2类：数值模拟法和近似解析法。其中，蒙特卡洛仿真（Monte Carlo simulation, MCS）法是最常用的数值模拟法[4-6]。该方法通过仿真、实验得到大量样本点对应的结构响应（如位移、载荷等），并统计失效样本点数量与总样本点数量的比值来得到设计点的失效概率，求解精度较高。但是，计算机仿真（如有限元仿真、流体仿真）及实验的计算成本较高。

常用的近似解析法包括一阶可靠性方法（first order reliability method, FORM）[7-8]和二阶可靠性方法（second order reliability method, SORM）[9-10]。FORM 在最大可能点（most probable point, MPP）处对极限状态函数进行一阶泰勒展开。常用的FORM包括Hasofer等[11]提出的验算点法和Rackwitz等[12]提出的当量正态化法。目前，FORM 已在众多工程领域中得到广泛应用。薛自然等[13]运用一次二阶矩（first order second moment, FOSM）法对折臂机构的运动可靠度进行了研究；郑财等[14]运用FOSM 法分析了三轴数控机床加工精度的可靠性；巴振宁等[15]运用改进的一次二阶矩（advanced FOSM, AFOSM）法对埋地管道不同部位的失效概率进行了计算。然而，对于非线性可靠性问题，FORM的求解精度会降低。Chen等[16]在考虑FORM使用时最坏情况的基础上提出了FORM的精度分析方法。由于极限状态曲面在MPP处的曲率会对结构可靠性分析产生较大的影响，SORM通过对极限状态函数进行二阶泰勒展开来提高可靠性分析的精度，常用的SORM 为Breitung 方法[17-18]。相比于FORM，SORM 的求解精度得到了很大程度的改善。但Cai等[19]指出，Breitung方法的失效概率求解公式在特定情况下可能存在巨大的误差，甚至会发生求解错误。

综上，FORM虽然求解效率高，但其在以简单的超平面代替复杂的极限状态曲面时会造成精度损失；SORM的求解精度虽有所提高，但其近似公式可能会出现奇异解。为了解决上述问题，笔者拟提出一种基于二阶抛物线近似的可靠性分析方法。该方法在SORM的基础上以MPP作为抛物线顶点，根据极限状态函数在MPP 处的曲率来构建近似抛物线，并通过求解抛物线极限状态函数的可靠概率来代替原极限状态函数的可靠概率，旨在为结构可靠性分析提供新思路。

1 可靠性分析方法基础理论

1.1 一阶可靠性分析方法

根据结构随机变量的分布情况，结构的失效概率可表示为：

式中：Pf为结构的失效概率；g(x)为极限状态函数，其中x为随机变量，为随机变量的概率密度函数。

在可靠性分析中，须将相关非正态分布的随机变量x转换为独立的标准正态分布的随机变量u。当利用FORM 计算式(1)时，在标准正态空间中将极限状态函数g(u)在MPP 处进行一阶泰勒展开：

其中：

式中：g(u)为标准正态分布下的极限状态函数；uMPP为由MPP与坐标原点构成的向量。

极限状态函数g(u)的可靠度指标β的几何意义如图1所示。由图1可知，β为坐标原点到极限状态函数的最短距离[20]。在二维可靠性分析中，FORM采用直线方程代替极限状态函数，则极限状态函数g(u)的失效概率可表示为：

图1 二维可靠性分析中的可靠度指标与MPPFig.1 Reliability index and MPP in two-dimensional reliability analysis

uMPP通常采用验算点法来迭代求解，当β满足精度要求时，最后一次迭代得到的向量uk即为uMPP。由此可得，基于验算点法的β的迭代过程如下：

式中：βk为第k次迭代后对应的可靠度指标。

当极限状态函数的非线性程度较低时，上述迭代过程可以稳定收敛；但当极限状态函数为高度非线性时，迭代过程会因发生振荡而难以收敛。Yang 等[21]通过引入混沌控制法（chaos control, CC）法来保证迭代过程收敛，即通过严格控制迭代步长来实现稳定收敛。引入CC 法后uk的迭代过程可表示为：

式中：λ为步长控制因子，C为单位矩阵。

当极限状态函数为线性或远离坐标原点时，利用FORM可以求解得到准确的结果。然而，当极限状态函数的曲率较大时，FORM的近似值会与真实值产生较大偏差，此时须采用更精确的近似方法进行求解。

1.2 二阶可靠性分析方法

SORM对极限状态函数在MPP处进行二阶泰勒展开时考虑了极限状态函数的非线性程度。二阶泰勒展开后的极限状态函数可表示为：

式中：∇2g(uMPP)为MPP处的Hessian矩阵。

令单位向量αu和Q分别为：

利用单位向量αu构造正交矩阵A，使ATA=I。正交矩阵A的表达式如下：

SORM 通过计算Hessian 矩阵来求解极限状态函数在MPP处的主曲率，从而计算失效概率，具体公式如下：

其中：

式中：κi为第i个方向上的主曲率，eig(·)为矩阵的特征向量，(·)n-1为删去第n行和第n列的n-1 阶矩阵。

通常情况下，SORM在利用式(9)计算失效概率时须满足βκi＜1，然而实际应用时并非总是满足该条件。另外，当1-βκi的值过小时，会产生巨大的误差。

2 基于二阶抛物线近似的可靠性分析方法

2.1 抛物线方程的建立

通过上文分析可知，在二维情况下，当极限状态函数呈高度非线性时，简单地以直线方程代替极限状态函数并不能满足可靠性分析的精度要求，因此本文选择在MPP处构建抛物线来近似表示极限状态函数，如图2所示。

图2 二维可靠性分析中极限状态函数的抛物线拟合与旋转Fig.2 Parabola fitting and rotation of limit state functions in two-dimensional reliability analysis

根据定义，向量uMPP与FORM近似的极限状态函数相互垂直，故以极限状态函数的MPP为顶点，以向量uMPP为对称轴建立抛物线。对于二维可靠性问题，只需1个主曲率即可确定一条抛物线；对于n维问题，则需n-1 个主曲率来确定所需的抛物线方程。如图2(b)所示，将拟合得到的抛物线旋转至u1轴上，由于标准正态分布函数具有轴对称性质，旋转后的抛物线与原抛物线具有相同的失效概率。由此可得，基于极限状态函数MPP处曲率近似拟合的抛物线方程可表示为：

式中：ai为抛物线系数，由极限状态函数在MPP处的主曲率κi确定，两者的关系为-2ai=κi。

由此可得，最终结构的失效概率可表示为：

式中：φ(ui)为标准正态分布下的概率密度函数。

2.2 失效概率计算

如图3 所示，在二维情况下，具有相同值的概率密度函数对应的几何图形为圆。令圆的半径为ρ，则可将圆划分为两部分：一部分为位于抛物线外的可行域；另一部分为位于抛物线内的失效域。由角度比例可知，整个圆中失效域所占的比例为。则通过极坐标变化，式(11)可转换为：

图3 二维可靠性分析中可行域与失效域划分示意Fig.3 Schematic of feasible and failure domain division in two-dimensional reliability analysis

其中：

式中：θ为抛物线与圆的相交点与坐标原点构成的向量与坐标轴之间的夹角。

在n维情况下，具有相同值的概率密度函数对应的几何图形为超球体。该超球体同样可分为两部分：一部分是抛物面外的可行域，另一部分是抛物面内的失效域。失效域的面积S可通过球面积分来求解，其表达式为：

式中：fui为f对ui的一阶导数，其中f为超球面方程，f2+u21+u22+…+u2n-1=ρ2；V为投影区域，其为n-1维的超椭球体。

将式(14)转换为以下形式：

基于坐标变换可将V从n-1维椭球体变为n-1维超球面，则式(15)可进一步转换为：

在极坐标中，ui与ζi的转换关系如下：

其中：

式中：l1、l2、…、ln-1为对应椭球体的短半径。

根据计算得到的超球体与抛物面相交区域的失效域面积S，失效概率可视作失效区域面积与超球体面积的比值，即：

其中：

式中：Ss为超球体的面积。

3 实际算例

为了验证本文所提出的基于二阶抛物线近似的可靠性分析方法（下文简称为本文方法）的可行性以及准确性，结合3 个数值算例和1 个工程算例进行对比分析。采用本文方法对各算例进行可靠性分析，并与其他可靠性分析方法进行比较，包括FORM（选用CC法）、SORM（选用Breitung方法）和MCS法。

除了算例2 外，本文以MCS 法（样本数为1×106个）直接仿真模拟计算得到的值作为精确值，再分别采用FORM、SORM和本文方法计算可靠概率Pr（Pr=1-Pf）并进行对比。

3.1 算例1

算例1中的非线性极限状态函数g1(u)为：

式中：c、b为未知系数、参数。

极限状态函数g1(u)中的变量u1、u2均服从标准正态分布且相互独立，即u1～N(0,12),u2～N(0, 12)。极限状态函数g1(u)的图像如图4所示。当g1(u)=0时，极限状态函数g1(u)在标准坐标系下的图形为抛物线，且所有抛物线均经过相同的顶点(b, 0)，抛物线顶点是与坐标原点距离最近的点，故β=b。

图4 算例1的极限状态函数Fig.4 Limit state function of example 1

需要注意的是，随着极限状态函数中系数c和参数b的改变，FORM 对极限状态函数的近似结果始终为过抛物线顶点且与以b为半径的圆相切的一条直线，不随系数c的变化而变化；而SORM 对极限状态函数的近似结果为过顶点(b, 0)的曲线，会随着系数c的变化而变化。基于不同方法的算例1的可靠性分析结果如表1所示。由表1可知，在参数b相同、系数c不同的情况下，FORM 求解的可靠概率始终保持不变：当b=2 时，Pr=0.977 250；当b=3时，Pr=0.998 650。显然，该结果对于不同的极限状态函数来说是不合理的。当抛物线的非线性程度较高时，FORM与MCS法的可靠性求解结果的误差较大。对于不同的系数，本文方法求解得到的可靠概率比SORM的精确，且与精确解几乎一致。当取c=-0.25，b=2，3时，FORM 与MCS法的求解结果存在较大误差。当取c=- 0.25，b= 3 时，由于曲率过大，使得βκ1≥1，因此采用SORM 的近似公式求解可靠概率时，产生了虚数解，导致求解发生错误；当取c=- 0.25，b= 2 时，由于βκi=1，SORM求解得到的可靠概率为无穷大，同样发生求解错误。而本文方法在每种情况下均未发生求解错误，且求解得到的可靠概率与MCS法的求解结果相近，误差较小。

表1 算例1的可靠性分析结果Table 1 Reliability analysis results of example 1

3.2 算例2

算例2为圆轴在随机力矩下的可靠性分析。如图5所示，圆轴的一端处于固定状态，另一端受到外部力矩My、Mz（My、Mz的合力矩为M，M与y轴的夹角为ω）和扭矩T的作用。假设上述3个随机变量服从相互独立的标准正态分布，即：My～N(μ1,δ21)，Mz～N(μ2,δ22)，T～N(μ3,δ23)，其中μi为随机变量的均值，δi为随机变量的标准方差。假设拉应力和剪应力（σx为x方向的拉应力，τxy′为xy′面的剪应力）在圆轴固定端A处达到最大。根据最大剪应力准则，可得极限状态方程：

图5 随机力矩下的圆轴示意Fig.5 Schematic of circular axis under random torque

式中：Meq为作用在圆轴上的等效弯矩，[M]为许用弯矩，[σ]为许用应力。

假设上述3 个随机变量的标准方差满足δ1=δ2=δ3=δ，则可将外部力矩My、Mz和扭矩T标准正态化为变量u1、u2和u3，则式(23)可以转化为：

由式(24)可知，g2(u)服从三自由度的非中心卡方分布，其非中心参数，故算例2 无需采用MCS 法进行仿真。根据非中心卡方分布的性质，可直接得到g2(u)的可靠概率Pr的真实解析解：

其中：

利用解析公式(25)、FORM、SORM 和本文方法对随机力矩下的圆轴进行可靠性分析，结果如图6所示。由图6可知，随着η的逐渐减小，SORM求解得到的可靠概率Pr与真实解析解之间的误差逐渐增大，其求解精度甚至低于FORM，还出现了奇异解，因此不具有普适性。当采用本文方法求解可靠概率时，在η=2.5～3的情况下，其求解精度优于FORM和SORM，且本文方法在η≥3的情况下也能保证理想的精度。由此可知，本文方法在进行可靠性分析时更为稳定。

3.3 算例3

在算例3中，非线性极限状态函数g3(x)为：

极限状态函数g3(x)中各随机变量的分布如表2所示。

表2 算例3中随机变量的分布情况Table 2 Random variables distribution of example 3

算例3中的随机变量服从正态、对数正态和极值Ⅰ型分布，其可靠性分析结果如表3所示。由表3可知，与MCS 法相比，FORM 求解得到的可靠概率的误差较大；SORM 求解时未发生错误，求解精度较高，SORM 求解的可靠概率与MCS 法的精确解之间的相对误差为0.073 2%；当面向非正态分布的随机变量时，本文方法的求解精度仍较高，与SORM 相比，本文方法求解得到的可靠概率的精度更高，其与精确解的相对误差仅为0.001 0%。