Processing math: 0%

基于子集模拟法的破冰船结构可靠性分析

陶楷文, 刘斌

陶楷文, 刘斌. 基于子集模拟法的破冰船结构可靠性分析[J]. 中国舰船研究, 2024, 19(X): 1–9. DOI: 10.19693/j.issn.1673-3185.03688
引用本文: 陶楷文, 刘斌. 基于子集模拟法的破冰船结构可靠性分析[J]. 中国舰船研究, 2024, 19(X): 1–9. DOI: 10.19693/j.issn.1673-3185.03688
TAO K W, LIU B. Structural reliability analysis of icebreaker structure based on subset simulation method[J]. Chinese Journal of Ship Research, 2024, 19(X): 1–9 (in Chinese). DOI: 10.19693/j.issn.1673-3185.03688
Citation: TAO K W, LIU B. Structural reliability analysis of icebreaker structure based on subset simulation method[J]. Chinese Journal of Ship Research, 2024, 19(X): 1–9 (in Chinese). DOI: 10.19693/j.issn.1673-3185.03688
陶楷文, 刘斌. 基于子集模拟法的破冰船结构可靠性分析[J]. 中国舰船研究, 2024, 19(X): 1–9. CSTR: 32390.14.j.issn.1673-3185.03688
引用本文: 陶楷文, 刘斌. 基于子集模拟法的破冰船结构可靠性分析[J]. 中国舰船研究, 2024, 19(X): 1–9. CSTR: 32390.14.j.issn.1673-3185.03688
TAO K W, LIU B. Structural reliability analysis of icebreaker structure based on subset simulation method[J]. Chinese Journal of Ship Research, 2024, 19(X): 1–9 (in Chinese). CSTR: 32390.14.j.issn.1673-3185.03688
Citation: TAO K W, LIU B. Structural reliability analysis of icebreaker structure based on subset simulation method[J]. Chinese Journal of Ship Research, 2024, 19(X): 1–9 (in Chinese). CSTR: 32390.14.j.issn.1673-3185.03688

基于子集模拟法的破冰船结构可靠性分析

基金项目: 国家自然科学基金资助项目(52271326)
详细信息
    作者简介:

    陶楷文,男,1997年生,硕士生。研究方向:船舶与海洋工程结构安全与可靠性。E-mail:1303132929@163.com

    刘斌,男,1985年生,博士,教授。研究方向:船舶与海洋工程结构安全与可靠性。E-mail:liubin8502@whut.edu.cn

    通讯作者:

    刘斌

  • 中图分类号: U663.6

Structural reliability analysis of icebreaker structure based on subset simulation method

知识共享许可协议
基于子集模拟法的破冰船结构可靠性分析陶楷文,采用知识共享署名4.0国际许可协议进行许可。
  • 摘要:
    目的 

    船体结构总纵强度的小失效概率分析是一个复杂、高维和计算耗时长的问题,难以采用传统可靠性分析方法精确计算。为解决复杂结构可靠性分析中小失效概率难以精确计算的问题,建立一套求解结构小失效概率的精确计算流程。

    方法 

    该流程通过采用子集模拟法以减少总体样本需求,为覆盖样本空间的全部情况,每个子集内样本通过MCMC(马尔科夫链蒙特卡洛法)抽样生成,结合随机有限元法对样本进行精确计算,得到精确样本响应值,以此提高失效概率的计算精度。同时,引入Kriging动态代理模型,显著减少子集内样本调用有限元计算次数。

    结果 

    将该流程运用于一个破冰船船体结构可靠性分析,精确计算全寿期船体结构失效概率为9.58×10−6

    结论 

    提出的基于子集模拟法的结构可靠性分析流程可以精确计算小失效概率。

    Abstract:
    Objective 

    The 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.

    Method 

    The 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.

    Results 

    The 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.

    Conclusion 

    The accuracy and efficiency of the proposed method are validated through comparisons with the computational results under various parameter settings.

  • 随着北极航线的持续开发,极地重型破冰船越来越受关注,其在恶劣环境下的结构可靠性是破冰船设计的一项关键指标。较运输船舶而言,破冰船结构承载能力更强、失效概率更低,为准确评估破冰船在航行和破冰期间的失效概率,需建立一套针对小失效概率精确计算的可靠性分析方法。

    结构可靠性分析方法可以分为三类:近似法(如一阶、二阶矩法)、数值模拟法(如蒙特卡洛法)和代理模型法(如响应面法、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],对于高度非线性和高维特性的问题都展现出良好的拟合能力,可在每个子集内部结合使用代理模型法减少计算耗时。

    本文将结合数值模拟法、有限元法和代理模型法,充分利用其在复杂、高维、高度非线性下计算优势,针对隐式非线性的小失效概率问题,提出一套既能满足计算精度要求,又能在实际计算能力下实现的适用于结构小失效概率问题的可靠性分析流程。首先阐述可靠性计算方法的基本原理,然后选取一艘破冰船进行案例分析。先建立破冰船船体结构可靠性模型,再根据破冰船有限元模型和计算结果建立随子集空间变化而更新的动态代理模型,随后采用子集模拟法从动态代理模型中快速抽样计算,对其在航行和破冰工况下的小失效概率进行分析,准确高效地获得可靠性计算结果。

    结构的失效概率用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}})

    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算法,即对随机变量的每一维度分量进行随机抽样,若备选样本不满足接受率条件,则重新生成备选样本,直到满足条件。重复生成马尔科夫链的下一状态样本,直到遍历每个子集中的中间状态分布情况。

    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)

    同时,代理模型在子集内部通过若干个计算样本点进行自更新,并预测子集内部其余马尔科夫链生成的新样本响应值。

    破冰船遭受的外载荷包括静水载荷、波浪载荷和冰载荷。外载荷和破冰船自身结构的不确定性构成了一个不确定性系统,用于建立破冰船可靠性分析模型。分析破冰船结构可靠性时,由于中垂工况失效概率远大于中拱工况,因此本文以中垂工况为研究对象。

    静水载荷主要指由船体自重、承载物体和船体浮力等引起的静水弯矩 {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
    波浪弯矩Gumbelxn = 757.06λ = 53.03k = 1
    垂向冰载荷Weibullλ = 21.05k = 0.98\theta = 1
    下载: 导出CSV 
    | 显示表格

    随机产生的垂向波浪载荷为船舶结构强度可靠性研究中考虑的主要载荷,对结构强度影响最大。海船的波浪弯矩( {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

    根据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

    图  1  弯矩值的概率密度
    Figure  1.  The probability density of bending moment value

    由于中垂工况下船体的最大应力出现在船舯区域,因此取船舯舱段作为船体结构可靠性分析的研究对象(图2)。在舱段的两个端面施加xz向线位移约束和yz向角位移约束,然后采用MPC点施加弯矩。

    图  2  有限元计算模型
    Figure  2.  Finite element calculation model

    船体结构的不确定影响因素主要考虑材料不确定性,可采用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
    材料分布类型统计量参数
    ALognormalμ = 235σ = 14.1
    AH36, DH36Lognormalμ = 355σ = 21.3
    下载: 导出CSV 
    | 显示表格

    运用Abaqus软件进行有限元计算,为实现自动计算,采用Python编写脚本调用Abaqus后处理自动计算样本,单次样本平均计算时长为5 min。

    自动采样计算流程见图3,首先对有限元模型进行参数化处理并编译成inp文件格式,每次MC、MCMC采样过程将替换inp文件中相应参数,并将该inp文件传递给Abaqus后处理进行计算,对完成计算的样本调用后处理接口从odb文件中提取最大Mises应力( {g_{m,\max }} )、有效塑性应变( P{E_{{\text{value}}}} )数值。

    图  3  计算流程图
    Figure  3.  Calculation flowchart

    为提高计算效率,在整体计算流程中加入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}}

    子集模拟法中对计算精准性影响最大的3个因素为真实值计算次数、Kriging模型预测次数、MCMC方法中MMH算法的提议分布。为研究这三个参数的设置影响,基于航行工况开展算例分析,变参数的计算结果见表3,其中算例1的计算过程见表4,取多次计算结果的中位数为最终失效概率。结果表明提议分布与实际计算样本数量对失效概率计算有一定影响,但波动范围可控,需根据实际情况选取合适的提议分布和计算样本数量。

    表  3  变参数的计算结果
    Table  3.  Computational results with variable parameters
    算例计算样本预测样本提议分布标准差失效概率相对误差/%
    110990201.06×10−731.27
    210990302.47×10−727.75
    320980205.89×10−717.20
    下载: 导出CSV 
    | 显示表格
    表  4  算例1的计算过程
    Table  4.  Computational process for case study 1
    子集迭代代数静水弯矩/MPa波浪弯矩/MPa许用应力/MPa最大应力/MPa塑性应变/mm
    11162.199187.170278.137126.2850
    12172.380189.688234.217131.6840
    21 001154.950278.609219.545163.2510
    21 500160.946277.533215.232165.4100
    65 567181.863421.689211.107237.5014.03×10−7
    66 000159.870378.419222.995229.3339.72×10−8
    下载: 导出CSV 
    | 显示表格

    以算例1为例,共产生子集数6个,由MC产生初始样本1 000个,每个子集内通过MCMC生成1 000个样本,总共7 000个样本。由表4可知,算例1在第5 567个样本出现了第一次失效。

    图4为子集4和5中MCMC样本的生成情况。MCMC样本呈现以某值为中心的发散分布状态,并未出现明显的方向趋向,证明样本可靠,马尔科夫链收敛、稳定。

    图  4  样本分布情况图
    Figure  4.  Sample distribution chart

    图5展示了MCMC生成的样本值与响应值的相关性。载荷特征与响应值存在一定线性关系,而船体结构所表现出的随机性没有明显地随响应值变化而变化的趋势,在整个运算中始终稳定于初始值附近,符合实际情况。其次,相比服从正态分布的静水弯矩,服从Gumbel分布的波浪弯矩对样本响应值的影响程度更大。

    图  5  样本数据散点图
    Figure  5.  Scatter plot of sample data

    图6(a)中给出了在算例1设置参数下,6个子集的响应值按升序排列后的情况,图6(b)给出存在失效样本的子集内部的失效样本分布情况。具体各子集阈值 {b_i} 及失效样本数量见表5

    图  6  算例1计算结果
    Figure  6.  Calculation result of Case 1
    表  5  算例1子集阈值
    Table  5.  Subset threshold of Case 1
    子集阈值/MPa失效样本数
    MCS141.0220
    1156.3460
    2169.6570
    3183.2363
    4201.15724
    5215.634107
    下载: 导出CSV 
    | 显示表格

    图7展示了响应值频率图,图7(a)为破冰船船体结构最大应力的整体分布情况,图7(b)采用对数坐标刻度对图7(a)图中尾部分布情况进行局部放大,可更清晰地显示每个子集的阈值及分布情况。从图7可以看出,破冰船结构承载力的响应值并不是可以直接求解的已知分布,很难通过直接计算获得准确的失效概率值。

    图  7  响应值堆积直方图
    Figure  7.  Histogram of response value stacking

    表6给出了航行、破冰工况下失效概率计算值,从表中可以看出,航行、破冰工况完成计算分别需要计算样本总数为7 861和5 404,相比MC法所需的 100/{P_{\text{f}}} 个样本,计算量有明显减少。

    表  6  破冰船可靠性计算结果
    Table  6.  Icebreaker reliability calculations
    工况计算样本预测样本总样本总子集数失效概率Pf
    航行2217 6407 86171.06×10−7
    破冰1845 2205 40453.17×10−5
    下载: 导出CSV 
    | 显示表格

    对在航行和破冰工况下的破冰船船体结构进行可靠性评估,最终需要求得破冰船总体运行情况下的失效概率。以航行工况下所运行的时间 {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

    本文提出了一种子集模拟法、有限元法和代理模型相结合的方法,可以用于分析结构可靠性分析中的小失效概率问题。该方法根据子集模拟法计算原理和中间计算过程,对计算流程做出优化和改进,通过加入Kriging代理模型预测新样本响应值,降低了计算成本,极大提高了计算效率。

    以某破冰船为研究对象,计算了其在航行和破冰工况下的精确失效概率。计算结果表明,子集模拟法可达到MC法计算的精确程度,也不会出现计算量过大导致无法计算的情况。通过对比发现,以本文研究所涉及的样本量为例,子集模拟法结合代理模型法相比MC法直接模拟最高可减少106倍计算时长。对于小失效概率事件,子集模拟法的计算更为高效。

    子集模拟法的中间样本生成采用了改进MMH抽样算法,对于计算高维复杂模型是稳定且适用的。对比提议分布变化导致的失效概率计算误差,证明了在正态分布下提议标准差所导致的计算误差可控。对比实际计算样本数量与预测样本数量比值对计算结果的影响,表明实际计算样本数量与预测样本数量比值越大,结果越精确,所以应在计算成本可接受范围内尽可能增加实际计算样本数量,并采取更优的抽样方法。

    在计算过程中马尔科夫链有可能陷入局部收敛无法跳出,产生明显异常值。这些数值会导致子集模拟法的最终计算结果与期望值偏差很大,需提早发现并人为剔除出样本集。在上百次对简化模型的计算中得出,剔除明显异常值后取数据集中位数作为最终结果与MC法直接模拟结果之间的误差不超过20%,属于随机变量正常的波动范围,可以认为子集模拟法的结果是准确可靠的。

    该计算流程适用于各类结构可靠性分析计算,本文使用的Kriging代理模型依赖大量实际样本导致有限元计算过程消耗大量的计算时间,后续改进可采用更高效的抽样方法或样本量需求更少的机器学习方法替代Kriging代理模型进行快速预测,提高计算效率。

  • 图  1   弯矩值的概率密度

    Figure  1.   The probability density of bending moment value

    图  2   有限元计算模型

    Figure  2.   Finite element calculation model

    图  3   计算流程图

    Figure  3.   Calculation flowchart

    图  4   样本分布情况图

    Figure  4.   Sample distribution chart

    图  5   样本数据散点图

    Figure  5.   Scatter plot of sample data

    图  6   算例1计算结果

    Figure  6.   Calculation result of Case 1

    图  7   响应值堆积直方图

    Figure  7.   Histogram of response value stacking

    表  1   弯矩分布类型及统计量

    Table  1   Bending moment distributions and parameters

    载荷分布类型统计量参数
    静水弯矩Normalμ = 625.93σ = 178.84
    波浪弯矩Gumbelxn = 757.06λ = 53.03k = 1
    垂向冰载荷Weibullλ = 21.05k = 0.98\theta = 1
    下载: 导出CSV

    表  2   材料属性的不确定性

    Table  2   Uncertainty in material properties

    材料分布类型统计量参数
    ALognormalμ = 235σ = 14.1
    AH36, DH36Lognormalμ = 355σ = 21.3
    下载: 导出CSV

    表  3   变参数的计算结果

    Table  3   Computational results with variable parameters

    算例计算样本预测样本提议分布标准差失效概率相对误差/%
    110990201.06×10−731.27
    210990302.47×10−727.75
    320980205.89×10−717.20
    下载: 导出CSV

    表  4   算例1的计算过程

    Table  4   Computational process for case study 1

    子集迭代代数静水弯矩/MPa波浪弯矩/MPa许用应力/MPa最大应力/MPa塑性应变/mm
    11162.199187.170278.137126.2850
    12172.380189.688234.217131.6840
    21 001154.950278.609219.545163.2510
    21 500160.946277.533215.232165.4100
    65 567181.863421.689211.107237.5014.03×10−7
    66 000159.870378.419222.995229.3339.72×10−8
    下载: 导出CSV

    表  5   算例1子集阈值

    Table  5   Subset threshold of Case 1

    子集阈值/MPa失效样本数
    MCS141.0220
    1156.3460
    2169.6570
    3183.2363
    4201.15724
    5215.634107
    下载: 导出CSV

    表  6   破冰船可靠性计算结果

    Table  6   Icebreaker reliability calculations

    工况计算样本预测样本总样本总子集数失效概率Pf
    航行2217 6407 86171.06×10−7
    破冰1845 2205 40453.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.

图(7)  /  表(6)
计量
  • 文章访问数:  71
  • HTML全文浏览量:  11
  • PDF下载量:  15
  • 被引次数: 0
出版历程
  • 收稿日期:  2023-12-15
  • 修回日期:  2024-04-16
  • 网络出版日期:  2024-04-22

目录

/

返回文章
返回