基于四阶矩法的CRDM承压壳体共因失效可靠性分析

陈 鹏，朱翊洲，谢永诚

(上海核工程研究设计院有限公司，上海 200233)

0 引言

目前，在设计核电主设备时，通常采用分析法设计进行分析，并基于ASME BPVC 或RCC-M等规范中的相关设计准则进行结构完整性评定，尽管评定准则中会考虑一定的安全系数，但确定性分析法忽视了结构和运行环境中存在的客观不确定因素，如加工导致的材料力学性能和结构尺寸的随机性，不同工况下温度、压力以及机械载荷波动的不确定性等。为了研究这些不确定因素对设备安全性的影响，并为设备维修和优化设计提供理论指导，有必要在设计过程中充分考虑潜在失效模式和随机因素、开展概率可靠性研究。

近几十年来，国内外学者对概率可靠性理论和应用开展了大量的研究。其中，基于随机采样的Monte-Carlo方法[2-3]被认为是最精确的方法，但当可靠性模型规模较大时，该方法的求解效率难以满足对大样本的随机采样[1-2]；SONG等[4]基于Markov链技术和重要采样法，并通过引入失效域子集和条件失效概率提高了Monte-Carlo法在高可靠性问题中的求解效率，但对于工程中的复杂结构和非正态分布可靠性问题仍然很难适用。

为了满足工程需求，学者们提出了以一次二阶矩法[5]为代表的一阶可靠性方法(FORM)，并基于对验算点求解、非正态分布参数标准化及功能函数线性化等关键问题的研究进行了修正，但由于FORM中对功能函数作了正态分布假设，导致该方法在求解基本随机变量和功能函数偏度较大的可靠性问题时精度较低。在此基础上，学者们将FORM拓展到了高阶矩方法。杨周等[6-7]考虑非正态分布的随机变量，用随机摄动理论计算功能函数的前四阶中心矩，基于Edgeworth级数法给出功能函数的假设分布，并计算可靠度和灵敏度，该方法在功能函数和基本变量随机部分所做的摄动假设限制了其在高非线性功能函数和离散性较大的可靠性问题中的应用；李云贵等[8-9]基于最大熵原理，以功能函数前四阶原点矩作为约束条件，通过Lagrange乘子法计算功能函数的概率分布和可靠度，但对于复杂结构问题，优化求解过程难以收敛；在四阶矩理论中，Pearson分布族也被广泛应用于可靠性问题的求解[10]，根据功能函数的前四阶中心矩，可以直接计算对应的分布函数，然而在实际应用中，有12种类型的概率密度函数，且参数计算需要通过数值积分，难以在工程中应用；此外，ZHAO等[11]基于三次多项式近似，提出以四阶矩为拟合目标的标准正态化技术，在标准正态空间中计算结构可靠度，并证明在较大的偏度范围内均有很高的计算精度。

针对核电设备中失效形式复杂、非正态分布变量较多的结构可靠性问题，本文基于四阶矩标准正态化技术，提出一种适用于复杂核电设备多失效模式的结构可靠性分析方法。基于有限元和序列响应面法建立考虑共因失效的多失效模式可靠性功能函数，用点估计法计算功能函数四阶矩，在此基础上通过四阶矩标准正态化变换得到标准正态空间中的可靠度指标。通过对控制棒驱动机构(CRDM)承压壳体结构多失效模式的可靠度分析，验证本文方法的适用性，并与Monte-Carlo方法结果的对比，证明该方法的计算精度。

1 基于响应面法的多失效模式功能函数

对于工程中的结构可靠性问题，描述结构失效的功能函数通常可以表示为：

g(x)=S-σ(x)

(1)

式中，S为基于相关强度理论和设计规范的设计应力强度值；x为影响应力分量σ的所有随机变量，当S也是随机变量时，功能函数对应的基本随机变量x应包含S。

结构可靠度的定义可以用功能函数描述为：

R=P[g(x)>0]

(2)

对于工程中的大多数结构的可靠性问题，式(1)所定义的目标应力分量通常只能通过有限元计算，无法给出显式的功能函数，因此难以通过FORM或高阶矩方法计算结构可靠度。本文基于响应面法(RSM)，用不含交叉项的二次多项式h(x)对隐式功能函数g(x)进行拟合：

(3)

式中，xi为基本随机变量(i=1,…,n)；a,bi,ci(i=1,…,n)为2n+1个待定系数，定义为：υ=[ab1bn…c1…cn]T。

(4)

考虑到基本随机变量在每个试验点的分布概率不同，因此文中通过加权最小二乘法求解待定系数向量：

υ=(HTWH)-1HTWg

(5)

其中，W为试验点权重矩阵,为对角阵。

(6)

在核电设备的设计分析中，通常需要对结构不同位置和不同工况下的应力分量进行评定。在结构可靠性分析中，同样需要考虑多失效模式的结构可靠性问题。由于共因失效在功能函数描述上的复杂性，以往在多失效模式的结构可靠性分析中，学者们通常保守地忽略了共因失效的影响，假定各失效模式相互独立，而直接将各点或各失效模式进行串联系统处理，即结构可靠度可以表示为：

RInd=R1R2…Rk

(7)

式中，Ri(i=1,…,k)为第i种失效模式单独计算的可靠度。

张义民等[13]在多失效模式机械零部件可靠性设计的结构分析中，也提出了另一种近似可靠度，即考虑最危险环节的多失效模式可靠度：

RMin=Min[R1,R2,…,Rk]

(8)

这两种方法得到的可靠度实际上是两种假设的极端，即各失效模式完全相互独立和完全相关。

本文在考虑多失效模式结构可靠性问题时，假定各失效模式对应功能函数为gi(x)，i=1,2,…,k，则结构功能函数可以表示为：

g(x)=min[g1(x),g2(x),…,gk(x)]

(9)

即结构中功能函数值最小的失效点都可靠时，整个结构才算可靠。显然，在构建功能函数时，考虑多失效模式比单独计算各失效模式可靠度之后再考虑相关性更加合理。由于式(9)所描述的功能函数是不连续的，因此在传统基于功能函数Taylor展开的可靠性方法中难以适用。本文基于响应面法，在有限元基础上拟合得到各失效模式功能函数gi(x)的近似功能函数hi(x)，并以如下功能函数来描述多失效模式的结构可靠性：

g(x)=min[h1(x),h2(x),…,hk(x)]

(10)

2 可靠性功能函数的四阶矩描述

在结构可靠性问题中，首先认为表征结构失效的功能函数也是一个随机变量，结构失效概率或可靠度由功能函数的分布信息确定，实际工程中，由于功能函数往往比较复杂，且多数情况下基本随机变量不能给出具体分布形式，因此得到功能函数的概率密度函数(PDF)或累积分布函数(CDF)通常是不现实的。在此基础上，矩法被广泛应用于功能函数的随机性描述和可靠性分析中。

当给定基本随机变量的有限阶中心矩时，功能函数的中心矩通常有两种求解方式。其一是基于随机摄动理论和功能函数一阶Taylor展开的摄动法；其二是标准正态空间的点估计法。

2.1 随机摄动法

随机摄动方法假定所有随机变量x和功能函数g(x)都可以表示为确定部分和随机部分，且随机部分远小于确定部分：

x=xd+εxp

(11)

g(x)=gd(x)+εgp(x)

(12)

其中，下标d表示变量或函数的确定部分，下标p表示随机部分，ε表示一个小量。

当随机变量的随机部分比确定部分小得多时：

(13)

则该功能函数的前四阶矩可以写为：

(14)

(15)

(16)

(17)

可以看出，摄动法需要知道功能函数在随机变量均值处的偏导数，此外，功能函数的各阶中心矩完全与基本随机变量对应阶中心矩相关，对于非线性功能函数显然会带来一定的误差。

2.2 点估计法

为了提高适用性和求解精度，本文采用另一种可靠性功能函数四阶矩求解方法，即点估计法，采用一系列标准正态空间的估计点和权函数近似计算功能函数的各阶中心矩。

当极限状态函数g=g(x1)只含有一个服从任意分布的基本随机变量时，g的均值和k阶中心矩为：

(18)

(20)

式中，T-1(u)为随机变量x1到标准正态分布变量u的Rosenblatt变换的逆变换过程[14]；uI，PI为估计点和对应的权值，可以通过对权函数为exp(-x2)的Hermite迭代得到[15]。

对于7点估计法，估计点和权值分别为：

u0=0，u1=-u2=1.154 405 4，u3=-u4=2.366 759 4，u5=-u6=3.750 439 7；P0=16/35,P1=P2=0.240 123 3，P3=P4=3.075 71×10-2，P5=P6=5.482 69×10-4。

当考虑n个基本随机变量时，功能函数g(x)=g(x1,x2,…,xn)可以近似表示为：

(21)

则g(x)的前四阶中心矩可以近似表示为：

(22)

(23)

(24)

(25)

各单变量函数的前四阶中心矩可根据式(18)～(20)得到。可以看出，用点估计法计算功能函数的前四阶中心矩不需要功能函数对基本变量的偏导数，因此更适用于复杂功能函数形式。此外，ZHAO等[14]指出，点估计法在估计功能函数高阶矩时需要更多的估计点才能保证精度，因此对于非正态性很强的基本随机变量，可能需要增加估计点个数。

本文基于摄动法四阶矩计算的可靠度结果,以后缀“_SD”表示，基于点估计法四阶矩计算的可靠度结果，以后缀“_PE”表示。

3 基于功能函数四阶矩的可靠度

当无法给出结构可靠度功能函数的概率密度函数，而只能得到前四阶中心矩时，要计算结构可靠度和可靠度指标，在统计学中需要通过有限阶中心矩对功能函数及其标准化正态变换形式做适当的简化。

首先，定义功能函数的标准化形式为：

(26)

一种较为常用的简化形式是采用Edgeworth级数法，将功能函数展开为标准正态分布函数的多项式形式：

(27)

式中，Hi(·)为i阶Hermite多项式；y为标准化功能函数zu；α3g，α4g为功能函数的偏度和峰度，α3g=M3g/(M2g)3/2,α4g=M4g/(M2g)2；Φ(·)和φ(·)为标准正态分布的累积分布函数和概率密度函数。

用Edgeworth级数法可以将结构可靠度近似表示为：

REW=P(g(x)>0)=F(βsm)

(28)

可以看出，Edgeworth级数法在拟合功能函数标准化变量的累积分布时作了尾部截断，试图用有限阶中心矩(如四阶矩)直接拟合变量的分布函数，这也导致了式(27)在很多情况下并不符合一个随机变量的累积分布特征。在计算结构可靠度时，可能会出现R>1的情况，对于这种情况，可以采用以下公式进行可靠度修正：

(29)

在关于四阶矩结构可靠性的研究中，ZHAO等[16]提出了另一种功能函数的标准正态化方法，即高阶矩标准化方法(HOMST)，对于标准化后功能函数分布变量zu，假定存在如下标准正态变换：

(30)

式中，aj为待定系数，通过计算Sx(zu)的有限阶中心矩，使其分别等于标准正态分布u的对应阶中心矩，可以得到各项系数。

当考虑u的三阶矩时，将功能函数的标准化变量zu通过如下关系式表示为标准正态变量：

(31)

在标准正态空间中，四阶矩可靠度指标可以表示为：

(32)

基于四阶矩可靠度指标的结构可靠度为：

RFM=Φ(βfm)

(33)

在本文的研究中，首先基于有限元和响应面法，根据结构的潜在失效模式分别建立各失效点的功能函数，并得到结构的多失效模式功能函数(见式(10))，采用点估计法计算功能函数的前四阶中心矩，最后用式(32)(33)得到可靠度指标和可靠度。

4 CRDM承压壳体结构可靠性分析

控制棒驱动机构(CRDM)是核反应堆中的关键设备，某磁力提升型CRDM结构如图1所示，主要由棒行程壳体、棒位探测器、钩爪部件、磁轭线圈、驱动杆部件和钩爪壳体结构组成。其中，棒行程壳体和钩爪壳体所组成的承压壳体结构不仅为CRDM提供支撑，还是整个反应堆一回路的压力边界，其可靠性对反应堆的安全运行至关重要。

图1 CRDM承压壳体结构示意

本文研究承压壳体在设计工况下的结构可靠性。在Ansys中建立轴对称模型，其中局部结构如图2所示。分别考虑承压壳体与顶盖贯穿件异种金属焊缝处的强度失效(评定截面ASN1)、钩爪壳体中部强度失效(评定截面ASN2)、棒行程壳体和钩爪壳体Canopy焊缝处的强度失效(评定截面ASN3)。失效准则基于ASME BPVC第三卷NB分卷，本文考虑设计工况下的一次薄膜应力的规范限制。

图2 CRDM承压壳体有限元模型和评定截面

为了研究3个失效模式在设计工况下的一次薄膜应力强度的概率可靠性，考虑设计内压P、步跃冲击载荷F、材料设计应力强度Sm三个变量为基本随机变量，则3个失效模式一次薄膜应力可靠度的功能函数可以表示为：

gi(x)=Sm-σi(P,F) (i=1,2,3)

(34)

其中，x=[SmPF]T为基本随机变量，并且假定分别服从双参数Weibull分布、对数正态分布、正态分布，分布参数如表1所示。

表1 各基本随机变量分布参数

需要说明的是，设计内压P和步跃冲击载荷F的均值取自CRDM设计规范书。设计应力强度Sm的均值根据承压壳体材料牌号，取ASME规范中SA-182 F304LN在设计温度下的设计应力强度值，各参数变异系数均取0.1。考虑到商密要求，表1中的各分布参数与真实值之间均有一定的偏差。

为了得到显式的功能函数，根据一次应力有限元模型和响应面法可以分别拟合得到3个失效模式的近似功能函数hi(x),i=1,2,3。在实际运行过程中，三种失效模式都与基本随机变量x相关，即存在共因失效关系，则承压壳体共因失效模式的可靠度功能函数可以表示为：

g(x)=min[h1(x),h2(x),h3(x)]

(35)

表2列出本文点估计法所计算的功能函数前四阶中心矩，对应的功能函数偏度为-0.384，峰度为3.254，可以看出，当忽略三阶及以上中心矩时，用传统二阶矩方法得到的可靠度将会偏大。

表2 功能函数前四阶中心矩

为了研究几种常用可靠度计算方法的精度，并验证串联系统假设在共因失效问题中的保守性，表3列出了5种可靠度计算结果。

表3 CRDM承压壳体多失效模式可靠度

表3中，MC法为基于功能函数式(35)的106次Monte-Carlo仿真结果，假定为精确值并作为本文近似方法的验证对象；EW法为采用Edgeworth级数法，根据式(27)(28)得到的四阶矩可靠度；SM法为根据可靠度指标βsm得到的二阶矩可靠度；SS法为基于传统串联系统假设，分别以四阶矩可靠度指标法单独计算三种失效模式可靠度，并根据公式(7)得到的可靠度；FM法为本文采用的四阶矩可靠度指标法，根据公式(32)(33)得到的可靠度。

从表3可以看出，相对于Edgeworth级数法和二阶矩法，本文所采用的四阶矩可靠度指标法具有更高的计算精度。在考虑3个基本随机变量的情况下，本文方法只需要调用21次功能函数就能得到与Monte-Carlo法调用106次精度相接近的结果。此外，串联系统假设得到的多失效模式可靠度明显误差更大，且比本文考虑共因失效的可靠度保守。

图3示出了本文考虑共因失效的多失效模式可靠性分析方法(FM)和串联系统假设(SS)的可靠度的对比，可以看出，串联系统假设只有在结构整体可靠度较高时才能适用，结构可靠度越小，用串联系统假设得到的可靠度误差越大。而本文方法在不同分布参数下均与Monte-Carlo法的仿真结果相吻合。

为了研究非正态分布参数对二阶矩和四阶矩方法精度的影响，图4示出了只考虑失效模式2(ASN2强度失效)的可靠度与应力强度Sm形状参数间的关系曲线，图中虚线为Sm的偏度，实线为不同方法得到的可靠度(其中_SM表示二阶矩可靠度、_FM表示四阶矩可靠度、_SD表示摄动法求解功能函数中心矩、_PE表示点估计法求解功能函数中心矩)。

图3 不同Sm均值下FM和SS方法可靠度结果对比

图4 失效模式2可靠度与Sm形状参数间的关系曲线

由图4可以看出，对于服从Weibull分布的应力强度Sm，当形状参数逐渐增大时，其偏度也逐渐增大(负值表示分布向左偏)；当Sm的偏度较小时，分别用摄动法和点估计法得到的二阶矩可靠度及四阶矩可靠度都与Monte-Carlo结果有相近的精度；当偏度达到一定值时，由于忽略了功能函数的偏度，二阶矩法逐渐偏离Monte-Carlo结果，此外，基于摄动法的四阶矩可靠度与Monte-Carlo法的可靠度误差也明显变大，甚至大于二阶矩可靠度的误差。

5 结语

本文考虑CRDM承压壳体共因失效条件下的多失效模式可靠性问题，基于有限元和响应面法建立了非连续的多失效模式功能函数，在考虑任意分布参数的情况下，用点估计法计算功能函数前四阶中心矩，并依此计算四阶矩可靠度指标和可靠度。

分析结果表明，该方法能够有效求解复杂结构多失效模式的可靠性问题，相比于Edgeworth级数四阶矩法和传统二阶矩法都有更高的精度，在相同精度的条件下也比Monte-Carlo方法有更小的求解规模；在求解多失效模式结构可靠性问题时，本文方法也比传统基于独立失效假设的分析方法有更高的精度，适用于核电主设备的结构可靠性分析。