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

doi:10.3785/j.issn.1006-754X.2024.03.305

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

陈振中^,^,¹, 黄冬宇¹, 田娇^,¹, 李晓科², 吴子豪³

1.东华大学机械工程学院，上海 201620

2.郑州轻工业大学机电工程学院，河南郑州 450002

3.上海第二工业大学智能制造与控制工程学院，上海 200135

Structural reliability analysis method based on second order parabolic approximation

CHEN Zhenzhong^,^,¹, HUANG Dongyu¹, TIAN Jiao^,¹, LI Xiaoke², WU Zihao³

1.College of Mechanical Engineering, Donghua University, Shanghai 201620, China

2.College of Mechanical and Electrical Engineering, Zhengzhou University of Light Industry, Zhengzhou 450002, China

3.School of Intelligent Manufacturing and Control Engineering, Shanghai Polytechnic University, Shanghai 200135, China

通讯作者: 田　娇（1988—），女，河北保定人，实验师，硕士，从事结构设计与材料可靠性研究，E-mail: tianjiao@dhu.edu.cn

收稿日期: 2023-10-20 修回日期: 2023-11-10

基金资助:

中央高校基本科研业务费专项资金资助项目.  23D110316
国家自然科学基金资助项目.  52375236.  51905492
上海市自然科学基金资助项目.  19ZR1401600

Received: 2023-10-20 Revised: 2023-11-10

作者简介 About authors

陈振中（1986—），男，河南开封人，副研究员，博士，从事可靠性优化设计研究，E-mail:zhenzh.chen@dhu.edu.cn，https://orcid.org/0000-0002-4101-8966 , E-mail：zhenzh.chen@dhu.edu.cn

摘要

针对工程中常用的一阶可靠性分析方法在求解非线性程度较高的极限状态函数的可靠性时精度不足的问题，在二阶可靠性分析方法的基础上提出了一种基于抛物线近似的结构可靠性分析方法。首先，采用一阶可靠性分析方法迭代求解标准正态空间下的最大可能点。然后，以最大可能点与坐标原点构成的向量作为新坐标轴，基于极限状态函数最大可能点各方向上的曲率构建近似抛物线，以提高边界区域的近似精度。最后，根据标准正态分布概率密度对新坐标轴上的近似抛物线进行积分，以求解结构的可靠概率。通过4个算例来比较一阶可靠性方法、二阶可靠性方法与基于二阶抛物线近似的可靠性分析方法，以验证所提出方法的可行性。结果表明，当面对非线性程度较高的可靠性问题时，一阶可靠性方法的求解精度较低，二阶可靠性方法在特殊情况下会发生求解错误，而通过抛物线近似积分的方法可有效提高结构可靠性分析的精度并保证求解的稳定性。研究结果可为复杂结构的可靠性分析提供参考。

关键词： 一阶可靠性 ; 二阶可靠性 ; 抛物线近似 ; 可靠性分析

Abstract

Aiming at the problem that the first order reliability analysis method commonly used in engineering has insufficient accuracy in solving the reliability of limit state function with high nonlinearity, a structural reliability analysis method based on parabolic approximation is proposed on the basis of second order reliability analysis method. Firstly, the first order reliability analysis method was used to iteratively solve the most probable point in standard normal space. Then, the vector composed of the most probable point and the coordinate origin was taken as the new coordinate axis, and the approximate parabola was constructed based on the curvature of the most probable point in each direction of the limit state function to improve the approximate accuracy of the boundary region. Finally, the approximate parabola on the new axis was integrated according to the standard normal distribution probability density to solve the structural reliability probability. The first order reliability method, the second order reliability method and the reliability analysis method based on the second order parabola approximation were compared by four examples to verify the feasibility of the proposed method. The results showed that when facing the reliability problems with high nonlinearity, the accuracy of the first order reliability method was low and the second order reliability method might make solving errors in special cases, while the parabolic approximate integral method could effectively improve the accuracy of structural reliability analysis and ensure the stability of the solution. The research results can provide reference for reliability analysis of complex structures.

Keywords： first order reliability ; second order reliability ; parabolic approximation ; reliability analysis

PDF (1981KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

陈振中, 黄冬宇, 田娇, 李晓科, 吴子豪. 基于二阶抛物线近似的结构可靠性分析方法[J]. 工程设计学报, 2024, 31(1): 50-58 doi:10.3785/j.issn.1006-754X.2024.03.305

CHEN Zhenzhong, HUANG Dongyu, TIAN Jiao, LI Xiaoke, WU Zihao. Structural reliability analysis method based on second order parabolic approximation[J]. Chinese Journal of Engineering Design, 2024, 31(1): 50-58 doi:10.3785/j.issn.1006-754X.2024.03.305

传统的结构设计通常采用安全系数来保证产品达到预期的可靠度^[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　一阶可靠性分析方法

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

\begin{array}{l} P_{f} = P (g (x) \leq 0) = \\ \underset{g (x) \leq 0}{\underset{︸}{\int \dots \int}} f (x_{1}) f (x_{2}) \dots f (x_{n}) d x_{1} d x_{2} \dots d x_{n} \end{array}

(1)

式中：P_f为结构的失效概率； $g (x)$ 为极限状态函数，其中 x 为随机变量， $x = {[x_{1} x_{2} \dots x_{n}]}^{T}$ ； $f (x_{i})$ 为随机变量的概率密度函数。

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

g (u) = (u - u_{M P P}) \nabla g (u_{M P P})

(2)

其中：

u = {[u_{1} u_{2} \dots u_{n}]}^{T}

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

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

P_{f} = Φ (- β)

(3)

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 二维可靠性分析中的可靠度指标与MPP

Fig.1 Reliability index and MPP in two-dimensional reliability analysis

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

β_{k} = \frac{g (u_{k - 1}) - \nabla g {(u_{k - 1})}^{T} u_{k - 1}}{{‖g (u_{k - 1})‖}^{2}}

(4)

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

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

\{\begin{array}{l} u_{k} = u_{k - 1} + λ C (f (u_{k - 1}) - u_{k - 1}) \\ f (u_{k - 1}) = - \frac{g (u_{k - 1}) - \nabla g {(u_{k - 1})}^{T} u_{k - 1}}{{‖\nabla g (u_{k - 1})‖}^{2}} \nabla g (u_{k - 1}) \end{array}

(5)

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

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

1.2　二阶可靠性分析方法

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

\begin{array}{l} g (u) = {(u - u_{M P P})}^{T} \nabla g (u_{M P P}) + \\ \frac{1}{2} {(u - u_{M P P})}^{T} \nabla^{2} g (u_{M P P}) (u - u_{M P P}) \end{array}

(6)

式中： $\nabla^{2} g (u_{M P P})$ 为MPP处的Hessian矩阵。

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

\{\begin{array}{l} α_{u} = - \frac{\nabla g (u_{M P P})}{‖\nabla g (u_{M P P})‖} \\ Q = - \frac{\nabla^{2} g (u_{M P P})}{‖\nabla^{2} g (u_{M P P})‖} \end{array}

(7)

利用单位向量 $α_{u}$ 构造正交矩阵 A，使 $A^{T} A = I$ 。正交矩阵 A 的表达式如下：

A = [A_{1} A_{2} \dots A_{n - 1} α_{u}]

(8)

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

P_{f} = \frac{Φ (- β)}{\sqrt[]{\prod_{i = 1}^{n - 1} (1 - β κ_{i})}}

(9)

其中：

κ = [κ_{1} κ_{2} \dots κ_{n - 1}] = e i g ({(A^{T} Q A)}_{n - 1})

式中： $κ_{i}$ 为第i个方向上的主曲率， $e i g (\cdot)$ 为矩阵的特征向量， ${(\cdot)}_{n - 1}$ 为删去第n行和第n列的n $-$ 1阶矩阵。

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

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

2.1　抛物线方程的建立

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

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 二维可靠性分析中极限状态函数的抛物线拟合与旋转

Fig.2 Parabola fitting and rotation of limit state functions in two-dimensional reliability analysis

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

g (u) = - u_{1} + \sum_{i = 1}^{n - 1} a_{i} u_{i}^{2} + β

(10)

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

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

\begin{array}{l} P_{f} = \int_{- \infty}^{\infty} \dots \int_{- \infty}^{\infty} φ (u_{2}) φ (u_{3}) \dots φ (u_{n}) \\ d u_{2} d u_{3} \dots d u_{n} \int_{β + \sum_{i = 2}^{n} a_{i} u_{i}^{2}}^{\infty} φ (u_{1}) d u_{1} \end{array}

(11)

式中： $φ (u_{i})$ 为标准正态分布下的概率密度函数。

2.2　失效概率计算

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

P_{f} = \int_{β}^{\infty} ρ e x p (- \frac{ρ^{2}}{2}) \frac{θ}{π} d ρ

(12)

其中：

θ = a r c c o s [(\sqrt[]{1 + 4 a_{i} β + 4 a_{i}^{2} ρ^{2}} - 1) / 2 a_{i} ρ]

(13)

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

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 二维可靠性分析中可行域与失效域划分示意

Fig.3 Schematic of feasible and failure domain division in two-dimensional reliability analysis

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

S = \underset{V : \frac{u_{1}^{2}}{a_{1}^{2}} + \dots + \frac{u_{n - 1}^{2}}{a_{n - 1}^{2}} \leq 1}{\underset{︸}{\int \dots \int}} \sqrt[]{1 + f_{u_{1}}^{2} + \dots + f_{u_{n - 1}}^{2}} d u_{1} \dots d u_{n - 1}

(14)

式中： $f_{u_{i}}$ 为f对 $u_{i}$ 的一阶导数，其中f为超球面方程， $f^{2} + u_{1}^{2} + u_{2}^{2} + \dots + u_{n - 1}^{2} = ρ^{2}$ ；V为投影区域，其为n $-$ 1维的超椭球体。

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

S = \underset{V : \frac{u_{1}^{2}}{a_{1}^{2}} + \dots + \frac{u_{n - 1}^{2}}{a_{n - 1}^{2}} \leq 1}{\underset{︸}{\int \dots \int}} \frac{ρ}{\sqrt[]{ρ^{2} - (u_{1}^{2} + u_{2}^{2} + \dots + u_{n - 1}^{2})}} d u_{1} \dots d u_{n - 1}

(15)

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

\begin{array}{l} S = \int_{r = 0}^{1} \int_{ζ_{1} = 0}^{π} \dots \int_{ζ_{n - 3} = 0}^{π} \int_{ζ_{n - 2} = 0}^{2 π} \frac{ρ}{\sqrt[]{ρ^{2} - (u_{1}^{2} + u_{2}^{2} + \dots + u_{n - 1}^{2})}} l_{1} \dots l_{n - 1} \\ r^{n - 1} s i n^{n - 3} ζ_{1} s i n^{n - 4} ζ_{1} \dots s i n ζ_{n - 2} d r d ζ_{1} \dots d ζ_{n - 2} \end{array}

(16)

在极坐标中，u_i 与 $ζ$ _i 的转换关系如下：

\{\begin{array}{l} u_{1} = l_{n - 1} ρ c o s ζ_{1} \\ u_{2} = l_{n - 2} ρ s i n ζ_{1} c o s ζ_{2} \\ u_{3} = l_{n - 3} ρ s i n ζ_{1} s i n ζ_{2} c o s ζ_{3} \\ ⋮ \\ u_{n - 2} = l_{2} ρ s i n ζ_{1} \dots s i n ζ_{n - 2} c o s ζ_{n - 2} \\ u_{n - 1} = l_{1} ρ s i n ζ_{1} \dots s i n ζ_{n - 2} s i n ζ_{n - 2} \end{array}

(17)

其中：

l_{i} = ρ s i n θ

(18)

式中：l₁、l₂、…、l_n $-$ ₁为对应椭球体的短半径。

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

P_{f} = {(2 π)}^{- n / 2} \int_{β}^{+ \infty} e x p (- \frac{ρ^{2}}{2}) S_{s} \frac{S}{S_{s}} d ρ

(19)

其中：

S_{s} = A_{n - 2} A_{n - 1} \dots A_{0} ρ^{n - 1}

(20)

A_{Z} = \{\begin{array}{l} π \frac{(2 m - 1)!!}{2 m!!}, Z = 2 m, m \geq 1 \\ 2 \frac{2 m!!}{(2 m + 1)!!}, Z = 2 m + 1, m \geq 1 \end{array}

(21)

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

3 实际算例

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

除了算例2外，本文以MCS法（样本数为 $1 \times 10^{6}$ 个）直接仿真模拟计算得到的值作为精确值，再分别采用FORM、SORM和本文方法计算可靠概率P_r（P_r=1-P_f）并进行对比。

3.1　算例1

算例1中的非线性极限状态函数 $g_{1} (u)$ 为：

g_{1} (u) = - u_{1} + c u_{2}^{2} + b

(22)

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

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

图4

新窗口打开| 下载原图ZIP| 生成PPT

图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时，P_r=0.977 250；当b=3时，P_r=0.998 650。显然，该结果对于不同的极限状态函数来说是不合理的。当抛物线的非线性程度较高时，FORM与MCS法的可靠性求解结果的误差较大。对于不同的系数，本文方法求解得到的可靠概率比SORM的精确，且与精确解几乎一致。当取 $c = - 0.25$ ， $b =$ 2，3时，FORM与MCS法的求解结果存在较大误差。当取 $c = - 0.25$ ， $b = 3$ 时，由于曲率过大，使得 $β κ_{1} \geq 1$ ，因此采用SORM的近似公式求解可靠概率时，产生了虚数解，导致求解发生错误；当取 $c = - 0.25$ ， $b = 2$ 时，由于 $β κ_{i} = 1$ ，SORM求解得到的可靠概率为无穷大，同样发生求解错误。而本文方法在每种情况下均未发生求解错误，且求解得到的可靠概率与MCS法的求解结果相近，误差较小。

表1 算例1的可靠性分析结果

Table 1 Reliability analysis results of example 1

参数b	系数c	FORM	SORM	本文方法	MCS法
2	-0.25	0.977 250	—	0.948 829	0.948 827
	0.01		0.977 691	0.977 774	0.977 783
	0.1		0.980 772	0.981 377	0.981 457
	0.3		0.984 661	0.985 658	0.985 667
	1		0.989 825	0.990 830	0.990 851
	5		0.995 035	0.995 625	0.995 632
3	-0.25	0.998 650	—	0.994 863	0.993 330
	0.01		0.998 688	0.998 694	0.998 695
	0.1		0.998 932	0.998 956	0.998 957
	0.3		0.999 193	0.999 227	0.999 226
	1		0.999 489	0.999 519	0.999 516
	5		0.999 757	0.999 774	0.999 773

新窗口打开| 下载CSV

3.2　算例2

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

M_{e q}^{2} = M_{y}^{2} + M_{z}^{2} + T^{2} \leq {[M]}^{2} = {(\frac{π r^{3} [σ]}{4})}^{2}

(23)

式中： $M_{e q}$ 为作用在圆轴上的等效弯矩， $[M]$ 为许用弯矩， $[σ]$ 为许用应力。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 随机力矩下的圆轴示意

Fig.5 Schematic of circular axis under random torque

假设上述3个随机变量的标准方差满足 $δ_{1} = δ_{2} = δ_{3} = δ$ ，则可将外部力矩 $M_{y}$ 、 $M_{z}$ 和扭矩T标准正态化为变量 $u_{1}$ 、 $u_{2}$ 和 $u_{3}$ ，则式(23)可以转化为：

g_{2} (u) = \sum_{i = 1}^{3} {(u_{i} + \frac{μ_{i}}{δ})}^{2} \leq η^{2} = {(\frac{[M]}{δ})}^{2}

(24)

由式(24)可知， $g_{2} (u)$ 服从三自由度的非中心卡方分布，其非中心参数 $γ^{2} = \sum_{i = 1}^{3} μ_{i}^{2} / δ^{2}$ ，故算例2无需采用MCS法进行仿真。根据非中心卡方分布的性质，可直接得到 $g_{2} (u)$ 的可靠概率P_r的真实解析解：

\begin{array}{l} P_{r} = P (g_{2} (u) \leq η^{2}) = Φ (γ + η) - Φ (γ - η) + \\ \frac{e x p [- {(γ + η)}^{2} / 2] - e x p [- {(γ - η)}^{2} / 2]}{2 \sqrt[]{π} γ} \end{array}

(25)

其中：

\{\begin{array}{l} η = β + γ \\ γ = (1 - β κ_{1}) η = (1 - β κ_{2}) η \end{array}

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

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 算例2的可靠性分析结果

Fig.6 Reliability analysis results of example 2

3.3　算例3

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

g_{3} (x) = x_{2} - 8 100 \frac{(x_{1} + x_{4})}{x_{3}^{2}}

(26)

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

表2 算例3中随机变量的分布情况

Table 2 Random variables distribution of example 3

变量	均值	标准差	分布类型
$x_{1}$	60	6.0	正态
$x_{2}$	2 000	74.0	对数正态
$x_{3}$	24	1.2	对数正态
$x_{4}$	50	10.0	极值Ⅰ型

新窗口打开| 下载CSV

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

表3 算例3的可靠性分析结果

Table 3 Reliability analysis results of example 3

方法	可靠概率	可靠度指标 $β$	相对误差/%
MCS法	0.956 4	1.709 8
FORM	0.961 0	1.761 8	0.481 0
SORM	0.957 1	1.717 7	0.073 2
本文方法	0.956 5	1.711 7	0.001 0

新窗口打开| 下载CSV

3.4　算例4

算例4为一个工程算例。如图7所示，该算例所分析的对象为一个对称的屋顶桁架结构，其上杆为AD、DC、CF、FB，压杆为DE、FG，底部杆为AE、EG、GB，张力杆为CE、CG。假设屋顶桁架的框架上承受了1个均匀分布的载荷q，将载荷q转换为节点载荷P（ $P = q L / 4$ ）。根据结构力学，节点C处桁架的垂直挠度W可表示为：

W = \frac{q L^{2}}{2} (\frac{3.81}{A_{c} E_{c}} + \frac{1.13}{A_{s} E_{s}})

(27)

式中：L为底部杆AE、EG、GB的总长度，A_c和A_s分别为钢筋混凝土和钢的横截面积，E_c和E_s分别为钢筋混凝土和钢的弹性模量。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7 屋顶桁架结构示意

Fig.7 Schematic of roof truss structure

鉴于W应满足 $W \leq 0.03 m$ ，则该屋顶桁架结构挠度的极限状态函数可表示为：

g_{4} (x) = 0.03 - W = 0.03 - \frac{q L^{2}}{2} (\frac{3.81}{A_{c} E_{c}} + \frac{1.13}{A_{s} E_{s}})

(28)

式中： $x = [q L A_{s} A_{c} E_{s} E_{c}]$ ，所有随机变量相互独立，其分布如表4所示

表4 算例4中随机变量的分布情况

Table 4 Random variables distribution of example 4

变量	单位	均值	标准偏差	分布类型
$q$	N/m	20 000	1 400	正态
L	m	12	0.12	正态
$A_{s}$	m²	$9.82 \times 10^{- 4}$	$5.985 2 \times 10^{- 5}$	正态
$A_{c}$	m²	0.04	0.004 8	正态
$E_{s}$	Pa	$1.0 \times 10^{11}$	$6.0 \times 10^{9}$	正态
$E_{c}$	Pa	$2.0 \times 10^{10}$	$1.2 \times 10^{9}$	正态

新窗口打开| 下载CSV

基于不同方法的算例4的可靠性分析结果如表5所示。由表5可知，FORM的求解误差较大；本文方法的求解精度较高，优于SORM，由此验证了该方法的可行性。

表5 算例4的可靠性分析结果

Table 5 Reliability analysis results of example 4

方法	可靠概率	可靠度指标 $β$	相对误差/%
MCS法	0.990 5	2.344 5
FORM	0.992 3	2.421 2	0.181 0
SORM	0.991 7	2.397 3	0.128 1
本文方法	0.991 6	2.392 0	0.111 1

新窗口打开| 下载CSV

4 结论

1）本文利用CC法迭代求解了极限状态函数的MPP，并根据极限状态函数在MPP处各方向上的曲率来构建抛物线方程，以近似表示原极限状态函数。结果表明，所构建的抛物线能够很好地反映极限状态函数在MPP处的边界区域。

2）FORM在求解非线性程度较高的极限状态函数的可靠性时精度较低，SORM在特殊情况下会发生求解错误。通过3个数值算例和1个工程实例来比较FORM、SORM与本文基于二阶抛物线近似的可靠性分析方法，对比结果验证了所提出方法的有效性和可行性。

3）本文基于二阶抛物线近似的可靠性分析方法要求解极限状态函数的二阶导数，计算量较大，具有一定的局限性。在下一阶段，将进一步深入研究如何在减少计算量的同时提高可靠性的求解精度。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

陈振中

基于可靠性的设计优化中精确解耦与高效抽样技术研究

［D］.武汉：华中科技大学，2013：1-2.