Structural reliability analysis of icebreaker structure based on subset simulation method
-
摘要:目的
船体结构总纵强度的小失效概率分析是一个复杂、高维和计算耗时长的问题，难以采用传统可靠性分析方法精确计算。为解决复杂结构可靠性分析中小失效概率难以精确计算的问题，建立一套求解结构小失效概率的精确计算流程。
方法该流程通过采用子集模拟法以减少总体样本需求，为覆盖样本空间的全部情况，每个子集内样本通过MCMC（马尔科夫链蒙特卡洛法）抽样生成，结合随机有限元法对样本进行精确计算，得到精确样本响应值，以此提高失效概率的计算精度。同时，引入Kriging动态代理模型，显著减少子集内样本调用有限元计算次数。
结果将该流程运用于一个破冰船船体结构可靠性分析，精确计算全寿期船体结构失效概率为9.58×10−6。
结论提出的基于子集模拟法的结构可靠性分析流程可以精确计算小失效概率。
Abstract:ObjectiveThe analysis of small failure probability in the overall longitudinal strength of a ship's hull structure is a complex, high-dimensional and time-consuming task which is challenging to calculate accurately using traditional reliability analysis methods. To address the precise calculation of small failure probabilities in the reliability analysis of complex structures, this paper establishes a computational process for assessing structural small failure probabilities.
MethodThe Subset Simulation method is used to reduce the overall sample demand. Samples within each subset are generated through MCMC (Markov Chain Monte Carlo) sampling to cover all scenarios within the sample space, and precise computations are performed on these samples using the Stochastic Finite Element Method (SFEM), resulting in accurate sample response values which enhance the accuracy of the failure probability calculations. Furthermore, the introduction of the dynamic Kriging surrogate model significantly reduces the number of FEM computations within each subset.
ResultsThe process is applied to a case study on the reliability analysis of an icebreaker's hull structure, and the analyzed lifecycle structural failure probability is 9.58 × 10-6.
ConclusionThe accuracy and efficiency of the proposed method are validated through comparisons with the computational results under various parameter settings.
-
0. 引 言
随着北极航线的持续开发，极地重型破冰船越来越受关注，其在恶劣环境下的结构可靠性是破冰船设计的一项关键指标。较运输船舶而言，破冰船结构承载能力更强、失效概率更低，为准确评估破冰船在航行和破冰期间的失效概率，需建立一套针对小失效概率精确计算的可靠性分析方法。
结构可靠性分析方法可以分为三类：近似法（如一阶、二阶矩法）、数值模拟法（如蒙特卡洛法）和代理模型法（如响应面法、Kriging）[1]。在分析破冰船结构可靠性时，由于破冰船结构复杂且功能函数为隐式非线性函数，采用近似分析的方法不能精确表达结构可靠性和安全性，而数值模拟法可以在解决隐式功能函数问题的同时得到更精确的结果。其中，蒙特卡洛（Monte Carlo, MC）法因操作方便、计算效率高被研究者大量使用[2]。MC法需要抽取大量样本以满足计算结果的精确性，对于破冰船这类结构复杂的小失效概率问题，在实际工程运用中计算成本将成倍增加。MC法可通过引入重要抽样（IS）密度函数改进，通过增加落入失效区域的概率来达到减小计算量的目的，但同时Au等[3]的研究发现IS不适用于高维随机变量问题，高维联合概率密度与重要抽样概率密度的比值非常小，往往导致结构失效概率被低估。子集模拟法是研究小失效概率的一种高效的计算方法[4-5]，通过引入中间状态将失效概率转化为一系列条件概率的乘积，减少抽样样本数量需求，且适用于高维、复杂随机变量的情况，极大提升了计算效率。
对于计算结构复杂且庞大的物理模型响应值，有限元法（Finite Element Method, FEM）是一种灵活且通用的高效方法。由于工程结构涉及的材料特性和环境载荷具有随机性，为评估这些不确定性因素对结构响应的影响，需在FEM中引入随机变量和过程，量化随机因素对结构产生的影响，这种方法被称为随机有限元法（Stochastic Finite Element Method, SFEM）[6-7]。尽管子集模拟法已极大减少样本需求，但大量有限元计算仍是一个十分耗时的工作。因此，采用代理模型学习实际有限元模型，可快速获取有限元模型的响应值。在结构可靠性分析中，代理模型法是一种常用的快速预测方法[8-11]，对于高度非线性和高维特性的问题都展现出良好的拟合能力，可在每个子集内部结合使用代理模型法减少计算耗时。
本文将结合数值模拟法、有限元法和代理模型法，充分利用其在复杂、高维、高度非线性下计算优势，针对隐式非线性的小失效概率问题，提出一套既能满足计算精度要求，又能在实际计算能力下实现的适用于结构小失效概率问题的可靠性分析流程。首先阐述可靠性计算方法的基本原理，然后选取一艘破冰船进行案例分析。先建立破冰船船体结构可靠性模型，再根据破冰船有限元模型和计算结果建立随子集空间变化而更新的动态代理模型，随后采用子集模拟法从动态代理模型中快速抽样计算，对其在航行和破冰工况下的小失效概率进行分析，准确高效地获得可靠性计算结果。
1. 小失效概率计算基本原理
1.1 子集模拟法
结构的失效概率用Pf表示，可靠性分析中结构失效概率表示为
{P}_{\text{f}}={\displaystyle \iint \cdots }{\displaystyle {\iint }_{G({\boldsymbol{X}}) < 0}f({\boldsymbol{X}})\text{d}{\boldsymbol{X}}} (1) 式中： {\boldsymbol{X}}\left( {{x_1},{x_2}, \ldots ,{x_n}} \right) 为n个随机变量； f({\boldsymbol{X}}) 为X的联合概率密度函数； G({\boldsymbol{X}}) 为状态极限函数。F表示结构的失效域，找到一系列满足条件的子失效域 {F_m} \subset {F_{m - 1}} \cdots \subset {F_{m - i}} \subset {F_1} \subset {F_0} = F ，则结构的失效概率 {P_{\text{f}}}(F) 可表示为
{P_{\text{f}}}(F) = P({F_m}) = \prod\limits_{i = 1}^m {P({F_i}|{F_{i - 1}})} (2) 小失效概率事件可以表示为一系列更大失效事件条件概率的连乘形式，第i个状态的条件失效概率 {P_i}({F_i}|{F_{i - 1}}) 表示为
{P_i}({F_i}|{F_{i - 1}}) = \frac{1}{{{N_{i - 1}}}}\sum {{I_{{F_i}}}[g({x^{i - 1}})]} (3) 式中： {N_{i - 1}} 为第i个状态的样本数量； g({x^{i - 1}}) 为新样本值。
第i个状态的示性函数 {I_{{F_i}}}[ \cdot ] 表示为：
{I_{{F_i}}}[ \cdot ] = \left\{ \begin{aligned} & 0 && g({x^{i - 1}}) \leqslant {b_i} \\& 1 && g({x^{i - 1}}) > {b_i} \end{aligned}\right. (4) 中间状态 {F_i} 所需样本由马尔科夫链蒙特卡洛（MCMC）法生成，运用MCMC法可由已有样本生成符合目标分布的新样本。由于示性函数 {I_{{F_i}}}[ \cdot ] 的每次判断均需调用一次功能函数，为降低计算成本，本文采用Kriging代理模型快速预测部分新样本值 g({x^{i - 1}}) 。
1.2 马尔科夫链蒙特卡洛
MCMC抽样可以高效地从子集中抽取符合目标分布的新样本，对于联合概率密度函数为 \pi ( \cdot ) 的目标分布，马尔科夫链相邻状态间的转移概率密度函数为 p({x_{k + 1}}|{x_k}) ， {x_k} 为随机变量的当前状态， {x_{k + 1}} 为随机变量的下一状态，若马尔科夫链满足细致平稳条件：
\pi ({x_k})p({x_{k + 1}}|{x_k}) = \pi ({x_{k + 1}})p({x_k}|{x_{k + 1}}) (5) 且 p( \cdot ) 符合对称的提议分布，例如高斯分布，则 {x_k} 和 {x_{k + 1}} 均服从目标分布，相邻状态间的转移概率相同，马尔科夫链可收敛为概率密度函数 \pi ( \cdot ) 的稳态分布。
在子集模拟法中进行样本抽样，需选取合适的提议分布 p( \cdot ) ，保证马尔科夫链的稳定性[12]。Metropolis-Hastings（MH）算法是MCMC应用最广泛的算法，但由于MH算法从向量角度出发，抽样效率偏低，故采用改进的MMH算法。该算法基于随机游走思想，对随机变量每一维分量进行抽样，可以弥补高维随机变量抽样效率低、无法遍历整个样本空间的缺点，提高马尔科夫链生成效率，在子集模拟中有很好的适应度。
MMH算法根据当前状态样本中的一维分量 x_k^j(j = 1,2, \ldots ,n) 进行抽样，其中相应一维提议分布的概率密度函数为 q( \cdot |x_k^j) ，生成的备选样本为 \varepsilon _k^j ，则该备选样本的接受率 \alpha (x_k^j,\varepsilon _k^j) 可表示为：
\alpha ({x}_{k}^{j},\text{ }{\epsilon }_{k}^{j})=\mathrm{min}\left[1\text{，}\frac{q({x}_{k}^{j}|{\epsilon }_{k}^{j})\pi ({\epsilon }_{k}^{j})}{q({\epsilon }_{k}^{j}|{x}_{k}^{j})\pi ({x}_{k}^{j})}\right] (6) 随机数 u ～ U(0,1) ，若 u < \alpha (x_k^j,\varepsilon _k^j) ，则接受备选样本 \varepsilon _k^j 作为下一状态样本分量 x_{k + 1}^j ，否则拒绝该样本，继续以 x_k^j 为下一状态样本：
x_{k + 1}^j = \left\{ \begin{aligned} & \varepsilon _k^j,&& \alpha (x_k^j,{\text{ }}\varepsilon _k^j) > u \\& x_k^j,&& \alpha (x_k^j,{\text{ }}\varepsilon _k^j) \leqslant u \end{aligned}\right. (7) 本文将采用改进MMH算法，即对随机变量的每一维度分量进行随机抽样，若备选样本不满足接受率条件，则重新生成备选样本，直到满足条件。重复生成马尔科夫链的下一状态样本，直到遍历每个子集中的中间状态分布情况。
1.3 Kriging代理模型
Kriging代理模型是一种基于变异函数理论构造的最小估计方差的最优线性无偏估计模型[13]，通过对一组已知数据点的空间自相关性进行建模，来预测未知点的值，主要包括线性回归和随机过程两部分，模型表达式如下：
y(x) = \sum\limits_{i = 1}^p {{f_i}(x) \cdot {\beta _i}} + {\textit{z}}(x) (8) 其中， {f_i}(x) 为变量的多项式函数， {\beta _i} 为回归系数， {\textit{z}}(x) 是随机过程，其协方差方程为：
{Cov} [{\textit{z}}({x_i}),{\textit{z}}({x_j})] = \sigma _k^2R(\theta ,{x_i},{x_j}) (9) R(\theta ,{x_i},{x_j}) 为 {x_i} ， {x_j} 之间的相关函数，只与两点间空间距离有关，表示两样本点间的相关性。以高斯函数为例，未知点x处的预测值为：
y(x) = \beta + {[R(x,{x_1}), \ldots ,R(x,{x_n})]^{\text{T}}}{R^{ - 1}}(Y - F\beta ) (10) 在子集模拟当中，示性函数 {I_{{F_i}}}[ \cdot ] 通过Kriging代理模型预测，对预测值 {g_p}(\varepsilon _k^j) 所有大于阈值的样本采取拒绝策略，下一状态样本 x_{k + 1}^j 将表示为：
x_{k + 1}^j = \left\{ \begin{aligned} & \varepsilon _k^j,&&\alpha (x_k^j,\varepsilon _k^j) > u\& {g_p}(\varepsilon _k^j) \leqslant {b_i} \\& x_k^j,&&\alpha (x_k^j,\varepsilon _k^j) \leqslant u \\& x_k^j,&& {g_p}(\varepsilon _k^j) \leqslant {b_i} \end{aligned}\right. (11) 同时，代理模型在子集内部通过若干个计算样本点进行自更新，并预测子集内部其余马尔科夫链生成的新样本响应值。
2. 船体结构可靠性模型
2.1 载荷及结构随机性影响分析
破冰船遭受的外载荷包括静水载荷、波浪载荷和冰载荷。外载荷和破冰船自身结构的不确定性构成了一个不确定性系统，用于建立破冰船可靠性分析模型。分析破冰船结构可靠性时，由于中垂工况失效概率远大于中拱工况，因此本文以中垂工况为研究对象。
2.1.1 静水载荷
静水载荷主要指由船体自重、承载物体和船体浮力等引起的静水弯矩 {M_{{\text{sw}}}} ，静水弯矩的不确定性来源于不同航次、工况下装载量的随机性[14]，最大静水弯矩（ {M_{{\text{sw,max}}}} ）一般出现在船舯区域（0.4～0.6）L，其中L为船长，根据BV船级社规范[15]计算得最大静水弯矩为894.18 kN·m。Hørte等[16]提出静水弯矩的长期分布形式为正态分布：
f(x) = \frac{1}{{\sqrt {2\pi } \sigma }}\exp \left( - \frac{{{{(x - \mu )}^2}}}{{2{\sigma ^2}}}\right) (12) 式中：\mu 为均值；\sigma 为标准差。
中垂工况下，静水弯矩极值的均值与标准差可分别定义为最大静水弯矩的70%和20%，即 0.7{M_{{\text{sw,max}}}} 和 0.2{M_{{\text{sw,max}}}} ，具体数值见表1。
表 1 弯矩分布类型及统计量Table 1. Bending moment distributions and parameters载荷 分布类型 统计量参数 静水弯矩 Normal μ = 625.93 σ = 178.84 — 波浪弯矩 Gumbel xn = 757.06 λ = 53.03 k = 1 垂向冰载荷 Weibull λ = 21.05 k = 0.98 \theta = 1 2.1.2 波浪载荷
随机产生的垂向波浪载荷为船舶结构强度可靠性研究中考虑的主要载荷，对结构强度影响最大。海船的波浪弯矩（ {M_{\text{w}}} ）长期分布近似服从Weibull分布：
{F_{M{\text{w}}}}(x) = P\left[ {{M_{\text{w}}} \leqslant x} \right] = 1 - \exp \left[ { - {{\left( {\frac{x}{\lambda }} \right)}^k}} \right] (13) 式中：k为形状参数；\lambda 为尺度参数。
Chen等[17]、 Soares等[18]、Gaspar等[19]确定了波浪弯矩Weibull分布的参数，以设计寿命20年内遭遇波浪108次计算，即{v_{\text{w}}}{T_0} = {10^8}，Weibull累计分布函数可表示为
{F_{\text{w}}}({M_{\text{w}}}) = 1 - \exp \left[ - \ln ({v_{\text{w}}}{T_0}){\left(\frac{{{M_{\text{w}}}}}{{{M_{{\text{w}},\max }}}}\right)^{{h_{\text{w}}}}}\right] (14) 根据BV船级社规范[15]计算可得，该船的最大波浪弯矩（ {M_{{\text{w}},\max }} ）为976.84 kN·m。可靠性计算中，将载荷幅值的极值作为随机变量，当样本足够多时，波浪弯矩的极值分布趋于极值Ⅰ型，即Gumbel分布：
{F_{\text{w}}}(x) = \exp \left[ { - \exp {{\left( { - \frac{{x - {x_n}}}{\lambda }} \right)}^k}} \right] (15) 其中，参数 {x_n} ， \lambda 表示为
{x_n} = {M_{{\text{w,}}\max }}\frac{{\ln N}}{{\ln ({v_{{\text{w}}0}}{T_0})}} (16) \lambda = \frac{{{M_{{\text{w}},\max }}}}{{\ln ({v_{\text{w}}}{T_0})}} (17) 式中：N为波浪循环次数。计算得出的参数 {x_n} 和 \lambda 见表1。
2.1.3 冰载荷
根据BV船级社规范[15]，在船舯区域[（0.35～0.65）L]的垂向冰致弯矩 {M_{\text{i}}} 与艏部的垂向冰载荷Fi呈线性关系，根据目标船的几何参数可计算得出：Mi = 5.25 Fi。
根据Suyuthi等[20]和Wu等[21]研究的冰载荷长、短期统计分布规律，垂向冰载荷（Fi）的实际极值可用三参数Weibull分布进行描述：
f(x;\lambda ,k,\theta ) = \frac{k}{\lambda }{\left(\frac{{x - \theta }}{\lambda }\right)^{k - 1}}{{e} ^{ - {{\left(\tfrac{{x - \theta }}{\lambda }\right)}^k}}} (18) 其中，\theta 为位置参数。
垂向冰载荷分布的参数值见表1。相应的静水弯矩、波浪弯矩和冰致弯矩概率密度分布见图1。
2.2 船体有限元模型
由于中垂工况下船体的最大应力出现在船舯区域，因此取船舯舱段作为船体结构可靠性分析的研究对象（图2）。在舱段的两个端面施加x，z向线位移约束和y，z向角位移约束，然后采用MPC点施加弯矩。
船体结构的不确定影响因素主要考虑材料不确定性，可采用Lognormal分布进行描述：
f({\log _a}x) = \frac{1}{{\sqrt {2\pi } \sigma }}\exp ( - \frac{{{{({{\log }_a}x - \mu )}^2}}}{{2{\sigma ^2}}}) (19) 船体结构建造主要采用A，AH36和DH36级钢，相关概率分布参数见表2。
表 2 材料属性的不确定性Table 2. Uncertainty in material properties材料 分布类型 统计量参数 A Lognormal μ = 235 σ = 14.1 AH36, DH36 Lognormal μ = 355 σ = 21.3 3. 可靠性计算流程与参数设置
运用Abaqus软件进行有限元计算，为实现自动计算，采用Python编写脚本调用Abaqus后处理自动计算样本，单次样本平均计算时长为5 min。
自动采样计算流程见图3，首先对有限元模型进行参数化处理并编译成inp文件格式，每次MC、MCMC采样过程将替换inp文件中相应参数，并将该inp文件传递给Abaqus后处理进行计算，对完成计算的样本调用后处理接口从odb文件中提取最大Mises应力（ {g_{m,\max }} ）、有效塑性应变（ P{E_{{\text{value}}}} ）数值。
为提高计算效率，在整体计算流程中加入Kriging动态代理模型预测响应值。首先以MC法计算的初始子集内N个样本结果为基础建立一个静态代理模型，N为自定每个子集内样本数。随后以每个子集内部实际计算出新的样本结果对原本Kriging代理模型更新，以更新后的代理模型预测新产生样本，直到满足子集模拟停止条件。这里采用了序列采样方法，逐步添加和更新样本点，在每个子集内部按一定步长采样一个或多个样本点做有限元计算，传递响应值并更新代理模型，然后通过动态代理模型预测当前和其他子集内新样本的响应值。这里的步长指计算样本和预测样本之比，需合理调整步长以达到合适的计算效率和计算精度。
以船体结构某处出现塑性应变作为失效判断条件，采用子集模拟法计算失效概率方法如下：
初始通过MC法生成N个样本，按功能函数的响应值升序排列，第 （1-{P}_{b}）\cdot N 个样本为第1个子集的中间失效事件 {F_1} 界限值 {b_1} ， {P_b} 为中间状态条件失效概率，则此时 P({F_1}) = {P_0} ；然后使用MCMC法生成当前子集内的中间样本，使用马尔科夫链可以生成符合目标分布新的N个样本，计算该子集中第 （1-{P}_{b}）\cdot N 个样本的条件概率 P({F_2}|{F_1}) 。重复以上计算步骤 P({F_m}|{F_{m - 1}}) ，直到满足条件停止计算。
子集模拟的停止条件：对于失效状态 {F_m} ，在第m个子集内，把 {x^{m - 1}} 对应的有效塑性应变值按升序排列，其中至少存在 {P_b} \cdot N 个样本落入失效区域内，即 P{E_{{\text{value}}}} > 0 的样本个数大于 {P_b} \cdot N 个，停止子集模拟。根据此时第m个状态的条件失效概率 {P_m} ，可计算出小失效概率：
{P_f} = {P_m} \cdot {P_b}^{(m - 1)} (20) Zuev等 [22] 的研究表明，条件失效概率 {P_b} 处于0.1～0.3时，子集模拟计算效率较高。通过多次对该模型计算得出： {P_b} 取值越大，计算精度越高，计算效率越低，随 {P_b} 值取值增大，计算时长将会倍增，但计算精度并未显著提高。反之， {P_b} 取值越小，计算精度越低，计算效率越高，尤其当 {P_b} ＜0.05会导致计算误差大幅增高。为最优化本文所研究的小失效概率计算效率，取中间状态的条件失效概率 {P_b} 为0.1。
设定每一阶段样本数为n个，每一阶段阈值 {b_i} 为第i个子集中，按最大Mises应力值 {g_{m,\max }}({x_1},{x_2}, \ldots ,{x_n}) 升序排列后的第900个响应值 {g_{bi}} ，以该值作为下一阶段子集的阈值 {b_{i + 1}} 。
4. 计算结果分析
子集模拟法中对计算精准性影响最大的3个因素为真实值计算次数、Kriging模型预测次数、MCMC方法中MMH算法的提议分布。为研究这三个参数的设置影响，基于航行工况开展算例分析，变参数的计算结果见表3，其中算例1的计算过程见表4，取多次计算结果的中位数为最终失效概率。结果表明提议分布与实际计算样本数量对失效概率计算有一定影响，但波动范围可控，需根据实际情况选取合适的提议分布和计算样本数量。
表 3 变参数的计算结果Table 3. Computational results with variable parameters算例 计算样本 预测样本 提议分布标准差 失效概率 相对误差/% 1 10 990 20 1.06×10−7 31.27 2 10 990 30 2.47×10−7 27.75 3 20 980 20 5.89×10−7 17.20 表 4 算例1的计算过程Table 4. Computational process for case study 1子集 迭代代数 静水弯矩/MPa 波浪弯矩/MPa 许用应力/MPa 最大应力/MPa 塑性应变/mm 1 1 162.199 187.170 278.137 126.285 0 1 2 172.380 189.688 234.217 131.684 0 … … … … … … 2 1 001 154.950 278.609 219.545 163.251 0 2 1 500 160.946 277.533 215.232 165.410 0 … … … … … … … 6 5 567 181.863 421.689 211.107 237.501 4.03×10−7 6 6 000 159.870 378.419 222.995 229.333 9.72×10−8 以算例1为例，共产生子集数6个，由MC产生初始样本1 000个，每个子集内通过MCMC生成1 000个样本，总共7 000个样本。由表4可知，算例1在第5 567个样本出现了第一次失效。
图4为子集4和5中MCMC样本的生成情况。MCMC样本呈现以某值为中心的发散分布状态，并未出现明显的方向趋向，证明样本可靠，马尔科夫链收敛、稳定。
图5展示了MCMC生成的样本值与响应值的相关性。载荷特征与响应值存在一定线性关系，而船体结构所表现出的随机性没有明显地随响应值变化而变化的趋势，在整个运算中始终稳定于初始值附近，符合实际情况。其次，相比服从正态分布的静水弯矩，服从Gumbel分布的波浪弯矩对样本响应值的影响程度更大。
图6（a）中给出了在算例1设置参数下，6个子集的响应值按升序排列后的情况，图6（b）给出存在失效样本的子集内部的失效样本分布情况。具体各子集阈值 {b_i} 及失效样本数量见表5。
表 5 算例1子集阈值Table 5. Subset threshold of Case 1子集 阈值/MPa 失效样本数 MCS 141.022 0 1 156.346 0 2 169.657 0 3 183.236 3 4 201.157 24 5 215.634 107 图7展示了响应值频率图，图7（a）为破冰船船体结构最大应力的整体分布情况，图7（b）采用对数坐标刻度对图7（a）图中尾部分布情况进行局部放大，可更清晰地显示每个子集的阈值及分布情况。从图7可以看出，破冰船结构承载力的响应值并不是可以直接求解的已知分布，很难通过直接计算获得准确的失效概率值。
表6给出了航行、破冰工况下失效概率计算值，从表中可以看出，航行、破冰工况完成计算分别需要计算样本总数为7 861和5 404，相比MC法所需的 100/{P_{\text{f}}} 个样本，计算量有明显减少。
表 6 破冰船可靠性计算结果Table 6. Icebreaker reliability calculations工况 计算样本 预测样本 总样本 总子集数 失效概率Pf 航行 221 7 640 7 861 7 1.06×10−7 破冰 184 5 220 5 404 5 3.17×10−5 对在航行和破冰工况下的破冰船船体结构进行可靠性评估，最终需要求得破冰船总体运行情况下的失效概率。以航行工况下所运行的时间 {S_1} 对应失效概率 {P_{{\text{f}}1}} ，破冰工况下所运行的时间 {S_2} 对应的失效概率 {P_{{\text{f}}2}} 计算，破冰船船体结构总体失效概率 {P_{\text{f}}} 可以表示为：
{P_{\text{f}}} = {P_{{\text{f}}1}} \cdot \frac{{S1}}{{S1 + S2}} + {P_{{\text{f}}2}} \cdot \frac{{S2}}{{S1 + S2}} (21) 破冰船船体结构总体失效概率，按运行时间占比7∶3计算，得出最终破冰船失效概率为9.58×10−6。
5. 结 论
本文提出了一种子集模拟法、有限元法和代理模型相结合的方法，可以用于分析结构可靠性分析中的小失效概率问题。该方法根据子集模拟法计算原理和中间计算过程，对计算流程做出优化和改进，通过加入Kriging代理模型预测新样本响应值，降低了计算成本，极大提高了计算效率。
以某破冰船为研究对象，计算了其在航行和破冰工况下的精确失效概率。计算结果表明，子集模拟法可达到MC法计算的精确程度，也不会出现计算量过大导致无法计算的情况。通过对比发现，以本文研究所涉及的样本量为例，子集模拟法结合代理模型法相比MC法直接模拟最高可减少106倍计算时长。对于小失效概率事件，子集模拟法的计算更为高效。
子集模拟法的中间样本生成采用了改进MMH抽样算法，对于计算高维复杂模型是稳定且适用的。对比提议分布变化导致的失效概率计算误差，证明了在正态分布下提议标准差所导致的计算误差可控。对比实际计算样本数量与预测样本数量比值对计算结果的影响，表明实际计算样本数量与预测样本数量比值越大，结果越精确，所以应在计算成本可接受范围内尽可能增加实际计算样本数量，并采取更优的抽样方法。
在计算过程中马尔科夫链有可能陷入局部收敛无法跳出，产生明显异常值。这些数值会导致子集模拟法的最终计算结果与期望值偏差很大，需提早发现并人为剔除出样本集。在上百次对简化模型的计算中得出，剔除明显异常值后取数据集中位数作为最终结果与MC法直接模拟结果之间的误差不超过20%，属于随机变量正常的波动范围，可以认为子集模拟法的结果是准确可靠的。
该计算流程适用于各类结构可靠性分析计算，本文使用的Kriging代理模型依赖大量实际样本导致有限元计算过程消耗大量的计算时间，后续改进可采用更高效的抽样方法或样本量需求更少的机器学习方法替代Kriging代理模型进行快速预测，提高计算效率。
-
表 1 弯矩分布类型及统计量
Table 1 Bending moment distributions and parameters
载荷 分布类型 统计量参数 静水弯矩 Normal μ = 625.93 σ = 178.84 — 波浪弯矩 Gumbel xn = 757.06 λ = 53.03 k = 1 垂向冰载荷 Weibull λ = 21.05 k = 0.98 \theta = 1 
下载: 导出CSV
表 2 材料属性的不确定性
Table 2 Uncertainty in material properties
材料 分布类型 统计量参数 A Lognormal μ = 235 σ = 14.1 AH36, DH36 Lognormal μ = 355 σ = 21.3
下载: 导出CSV
表 3 变参数的计算结果
Table 3 Computational results with variable parameters
算例 计算样本 预测样本 提议分布标准差 失效概率 相对误差/% 1 10 990 20 1.06×10−7 31.27 2 10 990 30 2.47×10−7 27.75 3 20 980 20 5.89×10−7 17.20
下载: 导出CSV
表 4 算例1的计算过程
Table 4 Computational process for case study 1
子集 迭代代数 静水弯矩/MPa 波浪弯矩/MPa 许用应力/MPa 最大应力/MPa 塑性应变/mm 1 1 162.199 187.170 278.137 126.285 0 1 2 172.380 189.688 234.217 131.684 0 … … … … … … 2 1 001 154.950 278.609 219.545 163.251 0 2 1 500 160.946 277.533 215.232 165.410 0 … … … … … … … 6 5 567 181.863 421.689 211.107 237.501 4.03×10−7 6 6 000 159.870 378.419 222.995 229.333 9.72×10−8
下载: 导出CSV
表 5 算例1子集阈值
Table 5 Subset threshold of Case 1
子集 阈值/MPa 失效样本数 MCS 141.022 0 1 156.346 0 2 169.657 0 3 183.236 3 4 201.157 24 5 215.634 107
下载: 导出CSV
表 6 破冰船可靠性计算结果
Table 6 Icebreaker reliability calculations
工况 计算样本 预测样本 总样本 总子集数 失效概率Pf 航行 221 7 640 7 861 7 1.06×10−7 破冰 184 5 220 5 404 5 3.17×10−5
下载: 导出CSV
-
[1] MELCHERS R E, BECK A T. Structural reliability analysis and prediction[M]. Hoboken: John Wiley & Sons, 2018.
[2] 李雪剑, 秦斌, 肖艺峰, 等. 改进随机森林−蒙特卡罗法在A型液舱支座结构可靠性分析中的应用[J]. 中国舰船研究, 2022, 17(1): 147–153,165. doi: 10.19693/j.issn.1673-3185.02181 LI X J, QIN B, XIAO Y F, et al. An improved random forest-Monte Carlo method and application for structural reliability analysis of A-type independent liquid tank support structure[J]. Chinese Journal of Ship Research, 2022, 17(1): 147–153,165 (in Chinese). doi: 10.19693/j.issn.1673-3185.02181
[3] AU S K, WANG Y. Engineering risk assessment with subset simulation[M]. Singapore: John Wiley & Sons, 2014.
[4] AU S K, BECK J L. Estimation of small failure probabilities in high dimensions by subset simulation[J]. Probabilistic Engineering Mechanics, 2001, 16(4): 263–277. doi: 10.1016/S0266-8920(01)00019-4
[5] CHENG K, LU Z, XIAO S, et al. Estimation of small failure probability using generalized subset simulation[J]. Mechanical Systems and Signal Processing, 2022, 163: 108114. doi: 10.1016/j.ymssp.2021.108114
[6] 郑志宝. 一种新的随机有限元求解方法及其应用[D]. 哈尔滨: 哈尔滨工业大学, 2019. ZHENG Z B. A new method for solving stochastic finite element equations and its applications[D]. Harbin: Harbin Institute of Technology, 2019 (in Chinese).
[7] 雷宽. 基于随机有限元的结构可靠性分析及优化设计研究[D]. 武汉: 武汉理工大学, 2015. LEI K. The reliability analysis and optimal design of metal structure based on stochastic element method[D]. Wuhan: Wuhan University of Technology, 2015 (in Chinese).
[8] KAYMAZ I, MCMAHON C A. A response surface method based on weighted regression for structural reliability analysis[J]. Probabilistic Engineering Mechanics, 2005, 20(1): 11–17. doi: 10.1016/j.probengmech.2004.05.005
[9] CHENG J, LI Q S, XIAO R C. A new artificial neural network-based response surface method for structural reliability analysis[J]. Probabilistic Engineering Mechanics, 2008, 23(1): 51–63. doi: 10.1016/j.probengmech.2007.10.003
[10] LIU H T, CAI J F, ONG Y S. An adaptive sampling approach for Kriging metamodeling by maximizing expected prediction error[J]. Computers & Chemical Engineering, 2017, 106: 171–182.
[11] 史朝印, 吕震宙, 李璐祎, 等. 基于自适应Kriging代理模型的交叉熵重要抽样法[J]. 航空学报, 2020, 41(1): 223123. SHI Z Y, LYU Z Z, LI L Y, et al. Cross-entropy importance sampling method based on adaptive Kriging model[J]. Acta Aeronautica et Astronautica Sinica, 2020, 41(1): 223123 (in Chinese).
[12] 兰成明, 徐震乾, 马君明, 等. 子集模拟中马尔可夫链蒙特卡洛抽样算法比较[J]. 土木工程学报, 2022, 55(10): 33–45,79. LAN C M, XU Z Q, MA J M, et al. Comparison of Markov Chain Monte Carlo sampling algorithms in subset simulation[J]. China Civil Engineering Journal, 2022, 55(10): 33–45,79 (in Chinese).
[13] 刘婧, 王德禹. 基于SMOTE算法和动态代理模型的船舶结构可靠性优化[J]. 中国舰船研究, 2020, 15(5): 114–123. doi: 10.19693/j.issn.1673-3185.01657 LIU J, WANG D Y. Reliability-based design optimization of ship structure using SMOTE algorithm and dynamic surrogate model[J]. Chinese Journal of Ship Research, 2020, 15(5): 114–123 (in Chinese). doi: 10.19693/j.issn.1673-3185.01657
[14] SOARES C G, MOAN T. Statistical analysis of still water load effects in ship structures[J]. Transactions of the Society of Naval Architects and Marine Engineers (SNAME), 1988, 96(4): 129–156.
[15] Bureau Veritas. Rules for the classification of ships operating in polar waters and icebreakers: NR527 DT R04 E[S]. Paris, France: Marine & Offshore, 2021.
[16] HØRTE T, WANG G, WHITE N. Calibration of the hull girder ultimate capacity criterion for double hull tankers[C]//Proceedings of the 10th International Symposium on Practical Design of Ships and Other Floating Structures (PRADS). Houston, Texas, USA: American Bureau of Shipping, 2007: 1−5.
[17] CHEN N Z, SOARES C G. Reliability assessment for ultimate longitudinal strength of ship hulls in composite materials[J]. Probabilistic Engineering Mechanics, 2007, 22(4): 330–342. doi: 10.1016/j.probengmech.2007.05.001
[18] SOARES C G, IVANOV L D. Time-dependent reliability of the primary ship structure[J]. Reliability Engineering & System Safety, 1989, 26(1): 59–71.
[19] GASPAR B, TEIXEIRA A P, SOARES C G. Effect of the nonlinear vertical wave-induced bending moments on the ship hull girder reliability[J]. Ocean Engineering, 2016, 119: 193–207. doi: 10.1016/j.oceaneng.2015.12.005
[20] SUYUTHI A, LEIRA B J, RISKA K. Short term extreme statistics of local ice loads on ship hulls[J]. Cold Regions Science and Technology, 2012, 82: 130–143. doi: 10.1016/j.coldregions.2012.05.017
[21] WU G, KONG S, TANG W Y, et al. Statistical analysis of ice loads on ship hull measured during Arctic navigations[J]. Ocean Engineering, 2021, 223: 108642. doi: 10.1016/j.oceaneng.2021.108642
[22] ZUEV K M, BECK J L, AU S K, et al. Bayesian post-processor and other enhancements of Subset Simulation for estimating failure probabilities in high dimensions[J]. Computers & Structures, 2012, 92-93: 283–296.
计量
- 文章访问数: 71
- HTML全文浏览量: 11
- PDF下载量: 15
