能源是支撑当今社会经济发展的命脉，当今中国经济迅速发展, 能源需求旺盛。传统化石能源过量开采和使用导致生态环境日益恶化。核电作为一种清洁能源在优化当今能源结构，改善大气环境方面起到重要作用。防止核泄露是核电厂设计中首要关注的问题。安全壳作为重要核安全相关结构是防止放射性物质泄露的最后一道屏障^[1]，设计安全壳结构时需要考虑诸如地震、龙卷风等外部事件以及内部失水事故等工况^[2]。安全壳在设计基准工况和超设计基准工况下的完整性可以采用极限内压承载力进行量化^[3-5]。目前对于安全壳在事故工况下内压承载力研究要集中在试验研究和数值分析。

关于安全壳试验研究方面，Hessheime等^[6]在桑迪亚实验室展开了1:4预应力混凝土安全壳模型试验。Parmar等^[7]进行了1:4 BARCOM安全壳模型试验，并且研究了安全壳在不同内压下的力学响应。Twidale等^[8]开展了1:10安全壳模型试验，并校核安全壳设计规范。Rizkalla等^[9]测试了加拿大某核电厂1:14安全壳模型。赵树明等^[10]对秦山核电二期安全壳进行了整体性试验，研究了在安全壳在0.35 MPa内压和145 ℃极端工况下的力学性能。

由于安全壳模型试验成本高并且周期长，采用数值模拟方法模拟安全壳非线性行为通常是不可避免的^[11]。Hora等^[12]提出了安全壳结构长期行为的数值模拟方法。Hu等^[13]在考虑长期预应力损失的情况下采用有限元分析了安全壳的极限承载力和破坏模式。Lee等^[14]在考虑混凝土受拉刚化效应的基础上评估了安全壳极限内压承载力。Saudy^[15]采用有限元方法分析了CANDU堆型安全壳的极限承载力。Lundqvist等^[16]评价了安全壳各种预应力损失计算方法，并提出改进预测方法。Shokoohfar等^[17]采用塑性损伤模型进行安全壳在高温和内压共同作用下的力学性能研究。Tong等^[18]采用ANSYS软件详细模拟了在预应力和温度作用下安全壳内压承载力。Huang等^[19]采用有限元软件真实模拟了安全壳在使用荷载以及假定事故工况下安全壳的性能并且考虑了诸如徐变，收缩等时变因素的影响。Zhang等^[20]采用有限元软件模拟5种蒸汽爆发工况下安全壳结构的完整性。

上述研究均属于安全壳性能确定性研究，这些研究可以确定安全壳的极限承载力和相应的破坏模式，但是这些确定性试验和数值模拟研究无法识别安全壳潜在的风险，以及无法了解安全壳在事故工况下的可靠性。由于安全壳结构在设计阶段、材料制备及施工和交付使用过程中不可避免的存在各种不确定性因素，从而导致安全壳结构的承载力具有很大的不确定性。同时传统安全壳设计规范ASME2007^[21]采用容许应力法进行设计，通过限制安全壳中钢筋、混凝土以及预应力筋的应力来保证安全壳的具有一定安全裕度。这种设计方法无法量化安全壳在各种超设计基准工况下的可靠性。同时目前国内外采用概率风险方法对安全壳在事故工况下性能评估的相关报道较少。Ellingwood等^[22-25]提出将可靠度理论用于核电厂结构分析。Pandey^[26]根据安全壳在役监测数据采用可靠度的方法分析了有粘结预应力安全壳的完整性。Kim等^[27]评价了长期预应力损失下安全壳的可靠性，并且分析了安全壳不同位置处的可靠性。上述关于安全壳可靠性研究是以安全壳结构功能函数为简单显示函数为前提。由于在实际情况下安全壳的功能函数十分复杂且无法显示表达。此外关于预应力混凝土安全壳在事故工况下性能评估的相关报道很少。

为此本文采用拉丁超立方抽样技术对安全壳在事故工况下的易损性进行评估，考虑了材料不确定性下安全壳发生功能性失效、结构性失效模式下的易损性，同时分析不同安全壳变截面位置处发生功能性失效时的易损性。最后采用概率评估方法分析了安全壳在严重事故工况下的性能。

1 安全壳有限元模型 1.1 安全壳几何模型及预应力体系

本文分析的某核电厂预应力混凝土安全壳由穹顶、环梁、扶壁、筒体以及底板等主要部分组成。同时为了保证安全壳的密封性在安全壳内表面设置6 mm钢衬里。安全壳筒体部分内径为22.5 m，壁厚为1.2 m。穹顶由曲率半径为30 m和8 m这2段圆弧构成，穹顶壁厚为1.0 m。安全壳筒体上有设备闸门、普通人员闸门和应急人员闸门3个主要洞口。在筒体的竖向和环向布置了2层普通钢筋和预应力筋。

安全壳的预应力体系由水平预应力钢束、纯竖向钢束以及由部分竖向预应力沿着穹顶延伸的倒U型钢束组成。安全壳的预应力体系采用后张法施工。水平预应力钢束采用两端张拉，并锚固在环梁上。竖向预应力钢束采用一端张拉，一端锚固在环梁上，另一端锚固在基础底板的预应力筋廊道上。倒U型钢束也采用一端张拉，一端锚固在对边的环梁上，另一端也锚固在预应力筋廊道上。安全壳结构简图如图 1所示。

1.2 预应力损失计算及模拟方法

核电厂安全壳的预应力损失会对安全壳产生不利影响，因此有限元分析需要考虑预应力损失。核电厂预应力混凝土安全壳预应力体系采用后张法施工，预应力损失由长期损失和短期损失2个部分构成。短期预应力损失包含混凝土弹性压缩、摩擦损失、预应力筋锚具变形造成的损失；长期损失包含徐变、收缩、松弛造成的损失。本文采用实体力筋法模拟混凝土预应力。预应力筋向混凝土的传递通过降温法模拟，预应力筋的有效预应力施加在预应力筋的每一个节点上。预应力筋每个节点的降温值为：

$ \Delta T = {\sigma _{{\rm{pe}}}}/(\alpha {E_p}) $

(1)

式中：σ_pe为考虑短期损失和长期损失后预应力筋的有效预应力；α为预应力筋的线膨胀系数，一般取为10^-5/℃；E_p为预应力筋弹性模量。

1.3 材料本构关系模型 1.3.1 混凝土本构关系

混凝土本构关系采用ABAQUS软件中经典的塑性损伤模型。塑性损伤模型假设混凝土2种破坏模式：一种是混凝土受拉开裂破坏，另一种是混凝土压碎破坏。混凝土破坏面以及刚度的折减通过等效塑性拉应变${\widetilde \varepsilon _t}^{pl}$和等效塑性压应变${\widetilde \varepsilon _c}^{pl}$ 2个变量控制。混凝土塑性损伤模型为：

$ {\sigma _{\rm{t}}} = (1 - {d_t}){E_0}({\varepsilon _t} - {\widetilde \varepsilon _t}^{pl}) $

(2)

$ {\sigma _{\rm{c}}} = (1 - {d_c}){E_0}({\varepsilon _c} - {\widetilde \varepsilon _c}^{pl}) $

(3)

式中：σ_t、σ_c分别表示混凝土的拉应力和压应力；d_t和d_c分别表示混凝土受拉和受压损伤因子，具体计算方法见文献[28]；E₀表示混凝的初始弹性模量；ε_t和ε_c分别表述混凝土拉应变和压应变。混凝土单轴受压本构关系参考文献[29]。混凝土受拉本构关系采用文献[30]建议的模型。混凝土单轴受压本构关系如图 2所示。

1.3.2 钢材本构关系

普通钢筋、钢衬里、电气贯穿件均采用理想弹塑 性本构关系模型^[31-32]。预应力筋本构关系采用双线性本构关系。同时普通钢筋、钢衬里、电气贯穿件和预应力筋的泊松比均取为0.3。非预应力钢材本构关系和预应力筋本构关系如图 3所示(其中f_y、ε_y为非预应力钢材屈服强度和屈服应变；E_y为非预应力钢材弹性模量；f_p、f_u为预应力筋的屈服强度和极限强度；E_p为预应力筋的硬化模量)。

1.4 有限元网格划分、接触及网格敏感性分析

安全壳混凝土采用以C3D8R单元为主，洞口少量采用C3D6单元划分网格。预应力筋采用T3D2单元，钢衬里和电气贯穿件采用以S4R单元为主，洞口位置采用S3类型单元划分网格。普通钢筋采用SFM3D4和SFM3D3单元来划分网格。所用钢筋层间距、方位均根据安全壳实际施工图确定。普通钢筋单元和预应力筋单元全部嵌入到混凝土单元中，不考虑它们与混凝土的相对滑移。钢衬里采用ABAQUS软件中蒙皮的方法来模拟，钢衬里与混凝土共用节点且两者没有相对滑移^[33]。安全壳各个部件网格划分如图 4所示。由于混凝土开裂特性，混凝土结构有限元计算结果对网格尺寸敏感，为此本文进行了安全壳网格尺寸敏感性分析。图 5给出了3种网格尺寸下安全壳最大位移响应随内压变化的曲线。从图 5来看，网格尺寸0.6 m与网格尺寸0.8 m计算结果相差很小，网格尺寸1.0 m计算结果与网格尺寸0.6 m和0.8 m在内压小于0.8 MPa时，计算结果相差甚小，但在内压超过0.8 MPa以后差异明显。因此，本文选择0.8 m网格尺寸进行有限元分析。

1.5 随机变量统计特性

本文选取最为常见的材料不确定性作为不确定性研究对象，研究在考虑材料不确定性下安全壳的内压易损性。由于混凝土弹性模量、抗拉强度与抗压强度相互关联，因此本文在进行易损性分析时考虑了它们的关联关系。混凝土弹性模量E_c、抗拉强度f_t与抗压强度f_c之间关系可以参考欧洲规范计算。为了简化分析，在本文有限元分析中材料密度和泊松比不作为随机变量。所有随机变量概率统计特性列在表 1中。

1.6 安全壳破坏准则

安全壳的破坏模式可以分为功能性失效和结构性失效2种，安全壳的功能性失效可以定义为安全壳失去防止泄露的功能。结构性失效是从结构极限强度方面定义安全壳的破坏状态，结构性失效可以认为是结构的整体变形、应变超过相应的限值^[34]。本文采用文献[35]建议的预应力混凝土安全壳破坏准则。具体破坏准则可以归纳为：

1) 钢衬里应变达到0.4%(安全壳功能性失效准则)；

2) 普通钢筋应变达到0.4%(安全壳结构性失效准则)；

3) 预应力筋总应变达到0.8%(安全壳结构性失效准则)。

2 有限元计算结果分析 2.1 自重和预应力工况下安全壳位移响应

图 6给出了安全壳在施加重力荷载和预应力工况下安全壳整体位移分布云图。在施加重力荷载和预应力工况后，安全壳筒体部位以向内收缩变形为主，穹顶以向外膨胀变形为主。此时安全壳穹顶的位移最大，穹顶产生的向外变形约为18.6 mm，安全壳向内收缩变形大小约为10.8 mm。

2.2 安全壳内压作用下非线性位移响应

安全壳对应混凝土刚发生开裂时，最大位移发生在设备闸门两侧小范围区域，此时最大位移约为20.24 mm，安全壳其他部位的位移相对较小(如图 7(a)所示)。钢衬里达到破坏准则时，安全壳位移发生重分布，同时安全壳整体位移出现大幅度增长，此时安全壳的最大位移发生在穹顶顶部，安全壳最大位移约为120.8 mm(如图 7(b)所示)。在普通钢筋和预应力筋达到相应破坏准则时，安全壳位移分布规律与钢衬里达到相应破坏准则时相同。但在普通钢筋达到破坏准则时，安全壳的最大位移较钢衬里达到相应破坏准则时有较大幅度的增长(如图 7(c)所示)。预应力筋达到相应破坏准则时，安全壳整体最大位移与普通钢筋达到相应破坏准则时对应的最大位移增长幅度较小(如图 7(d)所示)。

3 安全壳易损性分析 3.1 拉丁超立方抽样技术

Monte-Carlo作为一种经典的不确定性分析方法，广泛运用在工程领域。由于传统Monte-Carlo方法需要计算工作量大，计算效率低。尤其对于安全壳这种复杂结构采用传统Monte-Carlo方法进行分析计算工作量巨大。为此本文采用高效率Monte-Carlo抽样方法进行易损性分析。拉丁超立方抽样技术是一种分层抽样方法，可以有效避免重复抽样，能以较小的样本量反映总体的变异性。因此，采用拉丁超立方抽样技术可以大大减少样本计算数量，节约计算时间。分层抽样法与重要抽样方法相似，它们都是使贡献大的抽样更多的出现。但是分层抽样方法不改变原来的概率密度函数，而是将抽样区间分成一些子区间，并使各个子区间的抽样点数不同，对贡献大的子区间抽取更多的样本^[36]。拉丁超立方抽样方法主要分为以下2个步骤：

1) 将每个随机变量x_i等分成K个子区间，即每个子区间的概率为x_i^k，可以取每个子区间的中点值作为样本的代表值为(如图 8(a)所示)：

$ F({x_i}^k) = \left( {k - 1/2} \right)/K, \;\;k = 1, 2, \ldots , K $

(4)

2) 从每个随机变量x_i抽取一个样本代表x_i^k按照随机编号进行排列，对所有随机变量样本都进行随机排列，会形成随机排列的拉丁方(如图 8(b)所示以)，将上述K×N矩阵的每一列为一个随机组合的样本序列。

3.2 安全壳易损性计算结果

根据上述破坏准则可以得到安全壳相应的内压承载力。图 9给出了安全壳发生功能性失效、结构性失效时以及不同变截面位置发生功能性失效对应的内压承载力分布。

文献[37]认为安全壳的内压承载力服从对数正态分布，因此安全壳在事故工况下的内压易损性可以表达为：

$ \begin{array}{*{20}{c}} {{p_f}(p) = \mathit{\Phi }\left( {\frac{{\ln \left( {p/{p_m}} \right)}}{\beta }} \right) = }\\ {\int\limits_{ - \infty }^p {\frac{1}{{\sqrt {2\pi } \beta p}}} \exp \left\{ { - \frac{{\ln \left( {p/{p_m}} \right)}}{{2{\beta ^2}}}} \right\}{\rm{d}}p} \end{array} $

(5)

式中：p为安全壳的内压荷载；p_m为安全壳内压承载力均值；β为对数标准差。根据上述安全壳内压承载力服从对数正态分布假定，本文采用最大似然法得到安全壳功能性失效、结构性失效以及不同变截面位置功能性失效准则下对应易损性曲线统计参数。表 2、表 3分别给出了安全壳发生功能性失效、结构性失效以及安全壳不同变截面位置发生功能性时对应易损性曲线统计参数。

同时，不同置信水平的安全壳内压承载力计算公式为：

$ {p_Q} = {p_m}\cdot{{\rm{e}}^{{\mathit{\Phi }^{ - 1}}}}^{\left( Q \right)\beta } $

(6)

式中：p_Q为不同分位值的内压承载力；Q为分位值水平；Φ^-1(·)表示正态分布的反函数。

图 10、表 4和表 5分别给出安全壳功能性失效、结构性失效以及安全壳不同变截面位置处功能性失效对应的易损性曲线以及相应5%分位、95%分位对应的内压承载力。从图 10(a)和表 4综合来看，安全壳的破坏失效由功能性失效控制。安全壳发生功能失效对应的内压承载力与发生结构失效对应的内压承载力差异较大。基于普通钢筋失效准则与预应力筋失效准则对应的内压承载力相差较小。从图 10(b)和表 5结合来看，不同变截面位置处发生功能性失效的次序为环梁底部、环梁顶部、穹顶变截面位置、截椎体部位。变截面部位发生功能性失效比安全壳整体发生功能性失效时的内压承载力大，安全壳环梁底部、穹顶变截面位置、截锥体部位发生功能失效对应的内压承载力大于安全壳整体发生功能性失效对应的内压承载力。(同时本文偏于安全，取安全壳整体发生功能性失效时对应5%分位和95%分位的内压承载力分别作为安全壳的上限和下限内压承载力)。因此可以得到安全壳的下限和上限内压承载力分别为1.234 9 MPa和1.362 6 MPa。

3.3 安全壳严重事故工况下性能评估

目前关于安全壳严重事故工况下的性能指标评价的相关研究报道较少。文献[39]指出安全壳确定性的性能指标为在堆芯熔化24 h后安全壳仍然能保持其完整性和安全壳在严重事故工况下的条件失效概率(CCFP)不超过10%。安全壳在事故工况下的条件失效概率可以表达为：

$ {\rm{CCFP}} = \int\limits_0^\infty {{F_{{\rm{fail}}}}} \left( p \right)f\left( p \right){\rm{d}}p $

(9)

式中：F_fail(p)为安全壳的内压易损性曲线；f(p)为安全壳事故压力分布曲线。

由于本文研究安全壳的设计压力与文献[40]中设计压力不同，因此本文采用修正文献[40]中给出的事故内压分布曲线分布的参数得到适合本文安全壳的事故内压分布易损性曲线。修正后适合本文使用的安全壳事故压力服从对数正态分布，事故压力均值为0.663 7 MPa，对数标准差为0.3。本文偏于安全的采用3.2节安全壳整体发生功能性失效对应的易损性曲线作为安全壳严重事故工况下的内压易损性曲线，对卷积积分格式的式(9)进行数值计算得到安全壳的条件失效概率为0.01，说明本文分析的安全壳可以很好地满足事故工况下的性能指标。

4 结论

1) 安全壳破坏由功能性失效控制，且对应安全壳发生功能性失效时对应的内压承载力与发生结构性失效时对应的内压承载力相差较大，就结构性失效而言，普通钢筋失效准则、预应力筋失效准则对应的内压承载力相差较小。

2) 安全壳不同变截面位置发生功能失效的次序为环梁底部、环梁顶部、穹顶变截面位置、截椎体部位，环梁底部、环梁顶部、穹顶变截面位置、截椎体部位发生功能性失效时对应的内压承载力比安全壳整体发生功能性失效对应的内压承载力高，说明这些部位相对安全壳整体性失效是偏于安全的。

3) 安全壳在严重事故工况下的条件失效概率为0.01远小于性能指标0.1，说明该安全壳可以很好的满足严重事故工况下要求的性能要求。