基于GPC的环肋耐压圆柱壳结构失稳概率分析

张毅博， 孙志礼， 赵中强， 赵经武

(1. 东北大学 机械工程与自动化学院, 辽宁 沈阳 110819； 2. 中车唐山机车车辆有限公司， 河北 唐山 063000；3. 中国人民解放军第93107部队, 辽宁 沈阳 110141)

随着对隐蔽性要求的不断提高，潜艇下潜深度逐渐增加，环肋耐压圆柱壳结构作为潜艇的基本结构，其稳定性问题随静水压力的增大而愈发突出[1].因此，开展深潜环肋耐压圆柱壳结构失稳概率分析具有重要意义和工程价值.

可靠性分析的是产品在特定条件下完成特定功能的能力.因此，本文采用可靠性方法评估环肋耐压圆柱壳结构的失稳概率.在可靠性领域，结构状态通常被分为两类：安全和失效.近年来，为提高可靠性分析效率，诸多学者提出了各种基于代理模型(如响应面[2]、Kriging模型[3]、支持向量机SVM[4]、人工神经网络ANN等)的分析方法.然而，现有可靠性分析方法主要针对具有连续响应的结构，而对于具有失稳(即屈曲)破坏的拱结构、桁架结构、薄壳结构等的分析研究甚少.尽管可以采用SVM和ANN处理结构状态二分类问题，然而SVM和ANN均需要大量已标注样本才能得到较为精确的分类结果.此外，SVM仍存在核函数、核参数以及损失函数较难选取等问题.

作为高斯过程的一个分支，高斯过程分类(GPC)不仅具有严格的统计理论基础，还能自适应获取超参数[5].然而，GPC执行推理的时间随训练样本数的增加呈立方增长.为此，学者们提出了各种自适应试验设计(DoE)策略，如Kapoor等[6]将最易分类错误的点(MEMP)作为新训练样本；Peng等[7]将最可能失效点(MPP)添加到DoE中.然而，现有DoE策略选择出来的训练样本容易聚集，即引起不必要的功能函数评估.此外，对于具有小失效概率结构的可靠性评估，尽管Yang等[8]提出了一种双循环且不依赖于MPP的基于Krging和重要抽样(IS)的可靠性分析方法，然而目前尚未发现相关基于GPC和IS的可靠性分析方法.

为此，本文提出了一种单循环的基于GPC和IS的自适应分析方法，既避免了双循环中优化算法易陷入局部最优解的缺陷，又保证了在分类面附近区域均匀取样.通过某一分段函数验证了所提方法的高效性与准确性，并将其用于评估某深潜环肋耐压圆柱壳结构的失稳概率.

1 基本理论

1.1 高斯过程分类

高斯过程分类通常使用类标签+1或-1处理二元分类问题，+1代表结构安全，-1代表结构失效.给定N个训练样本X(X=[x1,…,xN])及其类标签Y(Y=[y1,…,yN])，GPC通过建立X和Y间的映射关系(即潜函数f(x))来进行分类.

为对未知样本x*分类，需由式(1)计算x*相应潜变量f*的后验分布.

(1)

式中：f为训练样本X对应的潜向量(f=[f(x1),f(x2),…,f(xN)]T)；p(f*|X,x*,f)为f*的条件先验分布；p(f|X,Y)为f的后验分布.

为求上述积分，GPC假定条件分布p(f|X)服从高斯分布，即

p(f|X)=N(0,K)

(2)

式中，K为协方差矩阵(Kij=k(xi,xj))；k(·)为协方差函数，本文选用应用最广泛的平方指数协方差函数(也称高斯核函数)来计算K.

(3)

式中：σf用来控制局部相关程度；指数项用来表征xi和xj间的距离相关性，即如果xi和xj的距离相对于距离尺度l很小，则它们高相关，否则低相关.

由高斯过程可知，f*和f的联合分布也服从高斯分布，即

(4)

式中：Kx*=[k(x*,x1),…,k(x*,xN)]T；kx*=k(x*,x*).

因此，条件先验分布p(f*|X,x*,f)为

(5)

由贝叶斯规则可得f的后验分布p(f|X,Y)为

p(f|X,Y)=p(Y|f)p(f|X)/p(Y|X)

(6)

将式(5)和式(6)代入式(1)中，则样本x*的类标签y*为+1的概率为

(7)

然而，式(1)和式(7)没有解析解，需采用数值近似法求解上述积分，本文选用Laplace近似法.Laplace采用高斯分布近似p(f*|X,Y,x*)，即

(8)

式中：

因此，式(7)中的分类概率近似为

(9)

式中，σ(·)为响应函数，本文选择标准正态分布的累积分布函数作为响应函数.

显然，分类概率0.5为结构安全与否的界限.

1.2 基于重要抽样的可靠性分析理论

在可靠性分析中，失效概率定义为

(10)

式中：G(x)为结构的功能函数；fX(x)为x的联合概率密度函数；x=[x1,…,xM]，M为影响结构性能的不确定因素(即随机变量)的个数.

尽管可以采用各种基于Monte Carlo采样的方法计算失效概率，但当失效概率很小时(如Pf<10-4)，需大量的Monte Carlo样本(>107)才能得到较为精确的结果.为此，需采用重要抽样等方差缩减技术来减少可靠性评估中随机抽样样本的数量.根据重要抽样，式(10)可改写为

(11)

(12)

(13)

2 基于GPC和IS的自适应分析方法

2.1 所提自适应DoE策略

现有基于GPC的自适应DoE策略主要包含两种：基于最可能失效点(GPC+MPP)和基于最易分类错误点(GPC+MEMP).

在可靠性分析中，MPP是对失效概率影响最大的点，且MPP附近区域高斯过程分类器的拟合精度直接影响着失效概率估计的准确性.因此，文献[7]将MPP作为新训练样本，其可由式(14)近似获取.

(14)

由式(9)可知，分类概率越接近0.5，其类标签越容易分类错误.因此，文献[6]将满足式(15)的点，即最易分类错误的点MEMP作为新训练样本.

(15)

然而，上述选点策略存在一些缺陷：1)他们的候补样本集是基于Monte Carlo采样得到的，难以应用于小失效概率评估；2)选择出来的训练样本容易聚集，从而增加不必要的功能函数评估次数.为此，本文提出采用马尔科夫链蒙特卡洛法(MCMC)和欧式距离来避免上述缺陷.

由上述可知，高斯过程分类器在分类概率为0.5处的拟合精度对失效概率估计的准确性起决定性影响.因此，应将分类概率为0.5的样本作为候补样本.为此，本文采用MCMC法产生NCS个满足式(16)的样本并将其作为候补样本集.

|p(y=+1|X,Y,x)-0.5|≤[ε]

(16)

式中，[ε]取为0.01.

为避免训练样本聚集，本文提出采用欧式距离来保证取样均匀性，其中欧式距离为

d=|xc-xt|

(17)

式中：xc为候补样本；xt为现有训练样本.

因此，本文将满足式(18)的点作为新训练样本，

(18)

(19)

式中，dmin为xc距DoE中所有训练样本的最短距离.

2.2 基于核密度估计的重要抽样

重要抽样的核心是重要抽样密度函数h(x)的构造，理论上最优h(x)为

h(x)=If(x)fX(x)/Pf

(20)

(21)

式中：ω为带宽；K(·)为核概率密度函数，本文选取正态分布的概率密度函数作为K(·)，其表达式为

(22)

带宽ω可由最小化平均积分平方误差获取，

(23)

2.3 停止准则

现有停止准则主要有两种，然而它们由于其固有缺陷而难以被工程实际所采纳.

1) 功能函数调用次数达到某一规定的阈值[9].显然，若规定的阈值过小，则易错估失效概率；反之，引起不必要的功能函数评估.

2) 前后两次迭代失效概率估计的相对误差小于某一阈值[7]，即

(24)

显然，若选择的新训练样本对高斯过程分类器的精度改善甚微，则提前终止迭代且得到错误的失效概率估计.

为此，本文基于失效概率估计的稳定性，提出了一种新的、更加精确的停止准则，即

(25)

(26)

2.4 所提自适应分析方法

本文所提自适应分析方法的基本步骤如下：

步骤1 通过Nataf变换将非独立随机变量转换为独立标准正态变量，采用拉丁超立方抽样随机生成N个初始训练样本SDoE=[x1,x2,…,xN]并调用真实功能函数获取类标签YDoE=[y1,y2,…,yN].

步骤3 基于当前高斯过程分类器，应用MCMC法生成NCS个满足式(16)的候补样本(本文令NCS=3 000).根据式(18)选择新训练样本，调用功能函数获取其对应的真实类标签，并将新训练样本添加到DoE中，返回步骤2.

3 算例验证

以一分段函数为例，验证所提自适应分析方法的准确性及高效性，其表达式为[10-11]

(27)

式中，x1，x2相互独立且均服从标准正态分布.

取12个拉丁超立方样本作为初始训练样本，且每种方法各运行10次以消除由不同初始训练样本及候补样本的随机性所产生的影响，其中GPC+MPP和GPC+MEMP以基于核密度估计的重要抽样(KDE-IS)生成的重要抽样样本作为候补样本.不同方法运行10次的平均结果如表1所示.此外，图1a，1b还分别显示了由不同方法得到的失效概率估计和最大相对稳定性ε随迭代次数Nit的变化趋势.显然，由表1和图1b可知，所提方法在满足失效概率分析精度要求的同时，还需要更少的功能函数评估次数.图2比较了不同方法在调用真实功能函数28次时DoE中训练样本的分布及预测的高斯过程分类边界.显然，所提自适应分析方法选择的训练样本分布更均匀且能更好地收敛到真实分类边界，即需要调用真实功能函数次数更少.

表1 不同方法运行10次的平均结果Table 1 Average results of 10 runs of different methods

4 环肋耐压圆柱壳失稳概率分析

某深潜环肋耐压圆柱壳结构采用矩形肋骨，其截面如图3所示，其中，r为耐压壳体内径，t为耐压壳体厚度，l为肋骨间距，a为肋骨高度，b为肋骨宽度.各参数及环肋耐压圆柱壳所受外部压力p均服正态分布，相应分布参数如表2所示[12].

表2 各随机变量的分布特征Table 2 Distribution of random variables

为评估该型环肋耐压圆柱壳结构的失稳概率，需采用ANSYS中的更符合工程实际和更具应用价值的非线性屈曲分析来分析其稳定性.由于环肋耐压圆柱壳结构的对称性，仅建立其1/4有限元模型，如图4所示.模型材料为某高强度钢，其弹性模量为2.0×105MPa，泊松比为0.3，屈服强度为785 MPa.圆柱壳采用壳单元shell181建模，矩形肋骨采用beam188单元建模.模型一端固定约束，另一端允许轴向位移，1/4截面处施加对称约束，壳结构外表面上施加压力载荷.

由ANSYS中非线性屈曲分析得到的该型环肋耐压圆柱壳结构在稳定和失稳时的节点位移云图分别如图5、图6所示.其中，图6中该型环肋耐压圆柱壳结构的失稳类型为整体失稳.

通过建立该型环肋耐压圆柱壳结构参数化有限元模型及MATLAB和ANSYS的联合仿真，基于所提分析方法，评估其失稳概率.为此，在标准正态空间内通过拉丁超立方抽样随机选取30个初始训练样本并调用ANSYS仿真得到其相应最大节点位移；构造初始高斯过程分类器，根据2.1节所提自适应DoE策略选择新训练样本，调用ANSYS获取最大节点位移并加入到DoE中.经过238次循环仿真后，满足收敛条件，即迭代停止，所得该型深潜环肋耐压圆柱壳结构的失稳概率约为8.242×10-5.

5 结 论

1) 通过某一分段函数表明，所提方法在满足失效概率分析精度要求的同时还极大地提高了分析效率，避免了不必要的功能函数评估.

2) 将所提自适应分析方法应用于评估某型深潜环肋耐压圆柱壳结构的失稳概率，结果表明：所提方法能够很好地适用于工程实际，这为预防潜艇失稳、缩短新型潜艇的研制周期及解决潜艇下潜的稳定性评估提供了一种切实可行的分析方法.