0 引　言

我国部分水面舰船由于服役海区海况条件较好，实际航行时间也低于设计要求，同时加上维护保养得当，舰船平台状态良好，有很大的潜力超期服役^[1]。此类超役舰船一部分会被改装为实体靶船，用以部队训练、试验等科研活动。极限强度是舰船结构强度的重要方面，也是实体靶船生命力的重要保障，关系到靶船航行安全及试训活动组织。目前，国军标等相关水面舰艇规范及标准缺少实体靶船极限强度评估的明确方法^[2]。

船体的极限承载能力一般指总纵弯曲作用下船体梁的极限弯矩，常用的计算方法有解析法、简化逐步破坏法和非线性有限元法等^[3]。其中，非线性有限元法通过建立结构有限元模型，模拟结构的初始缺陷，并考虑材料非线性和几何非线性，计算精度高，被国际船舶与海洋工程结构大会（ISSC）作为极限承载能力预报的衡准方法。如果考虑材料机械性能、结构尺寸和腐蚀等参数的不确定性，极限承载能力为符合一定概率分布的随机变量^[4]。响应面法结合蒙特卡罗法随机抽样是一种求解随机变量统计特征值的较好方法，该方法先通过有限元数值计算并采用最小二乘法等方法拟合一个响应面来替代未知的、真实的极限状态方程，再根据随机变量的概率分布运用蒙特卡罗法抽样获得目标变量的概率分布，二者的结合使结构分析的工作量大大减少，克服了传统蒙特卡罗法庞大的有限元计算工作量和传统响应面因为线性化带来误差的缺点^[5]。

为了合理地评估实体靶船的承载能力，考虑材料机械性能和结构腐蚀等参数的不确定性，本文将采用响应面-蒙特卡罗法结合非线性有限元法研究建立一种实体靶船极限强度的计算方法，分析响应面设计点的选取方案，给出实体靶船极限承载能力的概率分布。

1 基本原理与方法 1.1 可靠性分析的响应面-蒙特卡罗法

工程结构的失效概率可以表示为^[6]：

${P_f} = P\left\{ {G(X) < 0} \right\} = \int_{{D_f}} {f\left( X \right){\rm{d}}X}\text{。} $

(1)

式中： $G\left( X \right)$ 是极限状态方程，即响应面方程，当 $G\left( X \right) < 0$ ，结构失效； ${D_f}$ 是与 $G\left( X \right)$ 相对应的失效区域； $f\left( X \right) = f\left( {{x_1},{x_2}, \cdots ,{x_n}} \right)$ 是基本随机变量 $X$ 的联合密度函数，当 $X$ 为一组相互独立的随机变量时，有 $f\left( X \right) = $ $ f\left( {{x_1},{x_2}, \cdots ,{x_n}} \right) =\displaystyle \prod\limits_{i = 1}^n {f\left( {{x_i}} \right)} $ 。

用蒙特卡罗法计算结构可靠性时，式（1）可表达为：

${\hat P_f} = \frac{1}{N}\sum\limits_{i = 1}^N {I\left[ {G\left( {{{{\hat X}}_i}} \right)} \right]} \text{。}$

(2)

式中： $N$ 为抽样模拟总数；当 $G\left( {{{\hat X}_i}} \right) < 0$ 时， $I\left[ {G\left( {{{{\hat X}}_i}} \right)} \right] = 1$ 。用相对误差来表示蒙特卡罗法的抽样误差有：

$\varepsilon = \frac{{\left| {{{\hat P}_f} - {P_f}} \right|}}{{{P_f}}} < 2\sqrt {\frac{{\left( {1 - {{\hat P}_f}} \right)}}{{N{{\hat P}_f}}}} \text{，}$

(3)

当给定 $\varepsilon {\rm{ = }}0.2$ 时，抽样数目 $N = {{100} /{{{\hat P}_f}}}$ ，由于 ${\hat P_f}$ 通常为一个较小的量，因此抽样数目 $N$ 要取很大才能保证蒙特卡罗法分析结果的准确性。如果将蒙特卡罗法与有限元法直接结合进行结构可靠性分析，将带来庞大的有限元计算工作量。

响应面法通过拟合一个响应面 $\hat Z = \hat G\left( X \right) = \hat G( {X_1}, $ $ {X_2}, \cdots ,{X_n} )$ 来替代未知的、真实的极限状态方程 $Z = G\left( X \right)$ ，进而利用JC等方法进行可靠度分析。用三次多项式表达包含2个随机变量的响应面函数的形式为：

$\hat Z = \hat G(x,y) \approx \sum\limits_{{\begin{array}{l} {i + j = 0}\\ {i + j \leqslant 3} \end{array}}}^3 {{p_{ij}}{x^i}{y^j}} \text{。}$

(4)

式中： ${p_{ij}}$ 是待定系数； $x$ 和 $y$ 分别是基本随机变量。

为使响应面函数更好地逼近真实状态，设计点选取非常重要。一方面尽量要少，以降低计算工作量，另一方面还要尽量包含较多极限状态信息，以保证响应面的拟合精度。

响应面-蒙特卡罗法^[7]首先根据随机变量的概率分布选取确定功能函数值的设计点。设计点一般以随机变量的均值点为中心，根据 $3\sigma $ 原则在 ${\ \mu _i} \pm k{\sigma _i}$ 范围内选取。进而，根据设计点方案运用有限元数值计算得到设计点对应的功能函数值。根据有限元计算获得的一系列功能函数值，采用最小二乘法等方法拟合响应面，来替代真实的极限状态方程。在此基础上，运用蒙特卡罗法对随机变量进行抽样，采用响应面进行 $N$ 次数值计算，得到 $N$ 个功能函数值 ${Z_j}\left( {j = 1,2, \cdots ,N} \right)$ ，进而确定功能函数的概率分布，或根据统计 ${Z_j} < 0$ 的个数，计算失效概率或可靠度指标。因此，响应面与蒙特卡罗法结合使结构分析的工作量大大减少，提高效率。

1.2 极限承载能力分析的非线性有限元法

基于非线性有限元法的船体极限强度计算精度主要取决于模型构建的准确性和非线性求解器的选择与参数设置^[8]。

为了避免边界条件的影响，模型范围一般选择舱段模型。板和主要纵向构件应采用4节点单元，且尽量保证网格形状接近正方形。为了准确模拟极限状态结构的塑性大变形和失效模式演化，骨材腹板一般要求至少布置3～4个单元，并以此为基准划分整个舱段模型的网格往往可以取得较好的效果。同时，极限强度计算需要考虑结构的初始挠度^[9]，包括：

板架的整体变形

${w_{S0}} = {B_0}\sin \frac{{\text{π} x}}{a}\sin \frac{{\text{π} y}}{S}\text{；}$

板的局部变形

${w_{P0}} = {A_0}\sin \frac{{m\text{π} x}}{a}\sin \frac{{\text{π} y}}{b}\text{；}$

骨材的侧倾变形

${w_{T0}} = {C_0}\frac{y}{{{h_w}}}\sin \frac{{\text{π} x}}{a}\text{。}$

式中： $a$ 和 $b$ 分别为横向和纵向构件的间距； $S$ 为纵桁间距； ${h_w}$ 为骨材腹板高度； $m$ 为半波数，取满足 ${a /b} \leqslant $ $ \sqrt {m\left( {m + 1} \right)} $ 的最小正整数； ${A_0}$ ， ${B_0}$ 和 ${C_0}$ 为初始挠度的幅值，分别取 ${b / {200}}$ ， ${a /{1\;000}}$ 和 ${a/ {1\;000}}$ ，mm。

由于极限强度计算需要考虑材料的非线性，在缺少船用钢材料拉伸曲线的情况下，一般采用理想弹塑性或理想塑性强化模型以考虑材料的非线性效应。

目前，极限强度分析的非线性求解方法主要有Newton-Raphson迭代法、Riks法和显式动态算法等。其中，Newton-Raphson法以载荷增量为控制量，根据节点力平衡判断收敛，但结构大变形引起的节点力不平衡容易造成不收敛。Riks法采用弧长作为控制变量，改进了Newton-Raphson法的迭代策略，计算效率和收敛性主要取决于弧长相关参数的设置。显式动态算法采用中心差分法，用上一步和当前步的结果计算下一步的计算结果，无需迭代运算，收敛性好，计算效率较高，目前被广泛地应用于结构静态和准静态的非线性响应分析。

2 实体靶船极限强度可靠性分析 2.1 计算模型

以某型超役水面舰艇为例，该船已服役近30年，将改装成为实体靶船。考虑材料屈服限和结构腐蚀的随机变化，采用响应面-蒙特卡罗法结合非线性有限元法，开展了极限强度可靠性分析。

选取船中舱段为分析对象，采用Abaqus通用有限元软件建立有限元模型，网格尺寸整体上取50 mm×50 mm，同时纵骨腹板划分3个单元，单元类型选用4节点薄壳单元S4R，如图1所示。该舰材料屈服限为235 MPa，杨氏模量为206 GPa，泊松比0.3

2.2 随机变量的概率分布

由于是实体靶船，因此结构腐蚀是重点考虑的一个随机变量，同时由于材料屈服限对于舰船极限强度影响较大，因此主要考虑腐蚀量和材料屈服限的随机变化。材料屈服限服从正态分布，变异系数取0.05，即 ${\sigma _y} \sim N\left( {235,\;138.06} \right)$ ^[4]。

国内外学者通过对结构腐蚀统计分析和试验研究，提出了一些腐蚀模型^[10]。秦圣平等^[11]在现有腐蚀模型比较分析的基础上，提出了一种适用于船体结构时变可靠性分析的非线性腐蚀模型，该模型可以较好模拟钢结构在腐蚀环境下的腐蚀损伤过程。表达式为：

$d(t) = \left\{ \begin{array}{l} 0,0 \leqslant t \leqslant {T_{st}}, \\ {d_\infty }\left\{ {\left. {1 - \exp \left[ - {{\left(\frac{{t - {T_{st}}}}{\eta }\right)}^\beta }\right]} \right\}} \right. ,{\rm{ }}{T_{st}} \leqslant t \leqslant {T_L} \text{。} \\ \end{array} \right.$

(5)

式中： ${d_\infty }$ ， $\ \beta $ ， $\eta $ ， ${T_{st}}$ 为待定参数，均服从正态分布，其统计特征值见表1。

将表中相关参数代入腐蚀模型，采用蒙特卡罗抽样，抽样结果如图2所示，获得了服役30年的腐蚀量服从正态分布， $d\left( {30} \right) \sim N\left( {1.67,\;4.36{\rm{E}} - 3} \right)$ 。

2.3 极限弯矩计算与响应面构建

采用显式动态算法分析该舰中垂状态的极限强度，图3为舱段的弯矩-转角曲线，图4为极限状态结构应力云图。可以发现，该舰的中垂极限弯矩为36.21 MN·m，极限状态该舰强横梁间的甲板板架发生了屈曲破坏。

根据材料屈服限和结构腐蚀量的概率分布及其统计特征值，采用 $3\sigma $ 原则确定49个响应面设计点，采用非线性有限元法计算获得了不同腐蚀量和材料屈服限组合下的极限强度，如表2所示。

从舱段极限强度的计算结果可以发现：随着材料屈服限的提高，极限强度逐步增加；随着结构腐蚀量的增加，极限强度不断降低。总体而言，材料屈服限对极限强度的影响较腐蚀量的影响要大。

在此基础上，采用三次多项式表达式对响应面进行拟合，设计了10种不同设计点方案，分别包括49，25，21，21，17，17，17，13，13和13个设计点，如图5所示。表3为相关设计点方案的响应面拟合结果，图6为部分设计点方案拟合得到的响应面。

为了评价响应面拟合效果，引入了局部拟合优度 $R_{Loc}^2$ 和整体拟合优度 $R_G^2$ 两个评价指标。其中，局部拟合优度反映当前设计点方案 $n$ 个设计点与该方案下拟合得到的响应面之间的拟合优度；总体拟合优度指全部49个设计点与 $n$ 个设计点方案拟合得到的响应面之间的拟合优度，表达式为：

$R_{Loc}^2 = 1 - \frac{{\displaystyle\sum\limits_{i = 1}^k {{{({z_i} - {{\hat z}_i})}^2}} }}{{\displaystyle\sum\limits_{i = 1}^k {{{({z_i} - \bar z)}^2}} }}\text{，}$

$R_G^2 = 1 - \frac{{\displaystyle\sum\limits_{i = 1}^{49} {{{({z_i} - {{\hat z}_i})}^2}} }}{{\displaystyle\sum\limits_{i = 1}^{49} {{{({z_i} - \bar z)}^2}} }}\text{。}$

式中： ${z_i}$ 为设计点的有限元计算结果，即原始数据； $\hat z$ 为根据响应面函数得到设计点的极限强度； $\bar z$ 表示原始数据的均值，即49个有限元计算结果的均值。

从响应面的拟合结果可以发现：前9个方案的局部拟合效果均大于0.99；同时，前8个方案的整体拟合优度均大于0.99，且方案1、方案4、方案6和方案8达到了0.9998，对应的设计点数量分别为49，21，17和13，说明响应面的整体拟合优度与设计点的选取息息相关，并不是设计点数量越多越好。方案10的设计点选取没有考虑材料屈服限和腐蚀量的相关性，拟合得到的响应面已失真，如图6（d）所示，局部拟合优度仅0.482。方案8虽然只选取了13个设计点，但其整体拟合优度与49个设计点一致，达到了0.9998，局部拟合优度甚至达到了1.0，拟合优度比方案2（25个设计点）和方案3（21个设计点）都要好。

因此，响应面设计点的合理选取，不仅可以有效降低计算工作量，同时可以保证响应面拟合精度。

2.4 极限强度的可靠性分析

获得极限强度的响应面后，结合相关参数的概率分布，采用蒙特卡罗法随机抽样，进一步分析实体靶船极限强度的概率分布。

首先，开展了蒙特卡罗法的收敛性分析，抽样次数分别取1000，10000，100000和1000000次，样本的统计特征值见表5。可以发现，当样本数量达到10000次，极限强度的均值与标准差已趋于稳定；当样本数量增加到100000次或更高，极限强度的统计特征值变化不大。

在正态分布的基础上，进一步采用对数正态分布和Weibull分布拟合极限强度的概率分布，拟合结果如图7所示。其中，正态分布和对数正态分布的分布函数分别为 $\hat Z \sim N\left( {36.33,\;1.75} \right)$ 和 $\ln \left( {\hat Z} \right) \sim N\left( {3.59,\;0.048} \right)$ ；Weibull分布的形状参数为37.16，尺度参数为21.47，分布函数为 $\hat Z \sim W\left( {37.16,\;21.47} \right)$ 。

可以发现：极限强度的随机分布与Weibull分布偏差较大。因此，实体靶船极限强度的概率分布更趋近于正态分布或对数正态分布。

3 结　语

本文采用蒙特卡罗法与响应面法相结合的可靠性分析方法，开展考虑材料屈服限和腐蚀影响的舱段极限强度研究，研究表明：在保证拟合优度的前提下，响应面设计点的合理选取，可以有效降低计算工作量；考虑材料屈服限和腐蚀联合影响的实体靶船极限强度概率分布更趋近于正态分布或对数正态分布。本文研究建立的基于响应面-蒙特卡罗法的可靠性分析方法，对于实体靶船极限强度评估具有指导意义。