基于贝叶斯估计的核电厂安全壳内压概率安全评估

2024-04-24田澳楠潘晓兰苏春阳

原子能科学技术 2024年4期

田澳楠,郑志,潘晓兰,苏春阳,王勇

(太原理工大学土木工程学院,山西太原 030024)

近些年来,频发的地震灾害给核电厂的安全性能带来了严重威胁,这些超出设计基准的地震事故可能会对核电厂的主要系统和部件带来不可预测的损坏。为此,许多研究人员对核电厂的安全性能进行了概率评估,以抵御核电厂在严重事故下的潜在影响。Prinja等[1]使用一阶可靠度分析法预测了安全壳在超压荷载下的失效概率,并分析了带钢衬里的安全壳和不带钢衬里的安全壳的结构可靠性。Kim等[2]对考虑长期钢筋束退化的安全壳结构进行了内压荷载下的可靠性分析,结果表明钢筋束数量与安全壳的可靠性呈正相关关系。Hoseyni等[3]采用中值压力法对安全壳在超压作用下的失效概率进行了综合评估,并给出了安全壳在特定压力下的失效概率和极限承压能力。Granger等[4]采用置信度的方法对安全壳进行了在内压荷载下的概率安全评估,并对安全壳的抗泄漏能力进行了分析。金松等[5]考虑材料不确定性、统计不确定性以及建模不确定性等因素并对安全壳开展了概率安全评估,最后使用泰勒展开方法得到了安全壳的总失效概率。Zhao等[6]在考虑结构不确定性的情况下提出了一种基于PE-IDA的地震易损性分析方法,并通过单自由度体系对该方法进行了验证。结果表明,结构参数的不确定性对易损性的影响与结构损伤程度呈正比关系。Kwag等[7]提出了一种基于完全采样的地震概率评估量化方法,该方法能够准确地表示安全壳结构地震易损性信息的相关性。王晓磊等[8-9]对我国核电厂安全壳进行了地震易损性研究,并基于地震风险解析函数和风险卷积函数评估了我国某核电厂安全壳地震风险,提出了适用于我国核电厂地震易损性及风险评估的方法。上述学者在对安全壳进行事故荷载或地震作用下的易损性分析时均假定易损性函数服从对数正态分布模型,然而该模型未能综合考虑认知、经验、统计等多种不确定性因素的影响,这可能会极大影响安全壳概率安全评估的准确性。当前,Gardoni等[10-11]提出了完整的易损性评估理论,其通过采用贝叶斯估计建立结构的概率能力模型和需求模型,并综合考虑模型所有相关的不确定性,最大限度地提高了易损性评估和概率安全评估的准确性,然而该理论还未应用于核电厂安全壳的易损性和概率安全评估。应用贝叶斯估计理论进行安全壳的易损性和概率安全评估面临两大难点:1) 建立安全壳内压概率需求模型时首先需要确定能够反映结构整体反应变化的参数;2) 准确的概率需求模型还需要兼顾模型的简洁性,这需要对模型进行修正。

为此,本文提出一种新的混凝土损伤参数来评估安全壳混凝土的整体损伤性能,同时利用该参数建立确定性的内压需求模型。然后采用贝叶斯估计方法建立安全壳的概率需求模型,并对模型参数进行优化。在此基础上对不同损伤状态下的安全壳开展易损性评估及概率安全评估。

1 安全壳结构有限元模型

1.1 安全壳几何模型介绍

预应力混凝土安全壳由穹顶、筒体、钢衬里、扶壁柱、普通钢筋、预应力筋等组成。安全壳结构总高度为69.0 m,其中筒体高度为48.0 m,穹顶高度为21.0 m。穹顶可以减小安全壳在极高内压下的不利影响,穹顶的跨度和曲率半径分别为40.0 m和20.0 m;筒体的内外径分别为40.0 m和42.2 m。双层普通钢筋分别为竖直和水平分布。为了保证安全壳结构的完整性,在安全壳筒壁上布置了预应力筋系统。预应力筋系统由330根预应力筋组成,其中190根水平预应力筋锚固在扶壁柱上,140根倒“U”型预应力筋锚固在底板上。为了保证安全壳的密封性,在安全壳的内表面设置了6 mm厚的钢衬里。此外,为了满足安全壳的功能要求,在筒壁上留有一个直径为7.0 m的设备洞口。安全壳的几何简图如图1所示。

图1 安全壳几何简图

1.2 安全壳材料属性及其本构关系

混凝土的强度等级为C50。一般来说,混凝土被认为是各向异性的材料,其本构关系是影响裂缝演化的关键因素。为了尽可能准确地模拟安全壳的非线性行为,采用了经典的混凝土塑性损伤(CDP)模型[12]。该模型采用各向同性损伤和各向同性拉压弹性的概念来模拟混凝土的非线性行为,可以很好地表示混凝土在受力情况下的损伤演变。CDP模型通过式(1)定义拉伸破坏和压缩破坏;非弹性应变以式(2)表示;等效塑性应变可用式(3)来实现。损伤参数的确定采用了Sidoroff破坏模型[13]。

(1)

(2)

(3)

普通钢筋、钢衬里及预应力筋采用理想弹塑性模型。泊松比均取0.3,钢衬里和预应力筋的弹性模量为2×105MPa,钢筋的弹性模量为1.95×105MPa。钢衬里和钢筋的屈服强度分别为320 MPa、400 MPa,预应力筋的极限抗拉强度为1 860 MPa。

1.3 安全壳有限元分析方法

本研究通过3个加载步骤来获得安全壳在内压作用下的易损性评估结果。第1步施加重力荷载;第2步通过降温法[14]给预应力筋施加预应力,如式(4)所示;第3步施加线性增加的内压。

(4)

式中:ΔT为预应力筋的降温温度;Δσ为控制应力;λ为线膨胀系数;E为预应力筋的弹性模量。

1.4 有限元网格划分及边界条件设定

安全壳有限元模型计算时混凝土容易产生局部应力集中现象,网格单元过密则会产生频繁的不收敛问题。因此在划分网格时,在保持精度的前提下选取适当的网格尺寸。表1列出不同网格尺寸下完成安全壳预应力筋张拉时的相关统计信息。

表1 不同网格尺寸的统计信息

由表1可看出,采用0.8 m的网格尺寸可以在保证精度的前提下节约模型运行时间。为精确评估安全壳洞口的破坏表现,在洞口区域对网格进行加密处理。在划分网格单元时,混凝土采用实体单元(C3D8R)、钢衬里采用壳单元(S4R)模拟,两者通过共节点连接,且不考虑它们之间的相对滑移[15]。为减小钢筋建模的复杂性,普通钢筋采用表面单元(SFM3D4)模拟,预应力筋通过桁架单元(T3D2)模拟。安全壳的结构网格划分如图2所示。钢筋、预应力筋与混凝土采用分离式建模,其中钢筋和预应力筋完全嵌入混凝土中,不考虑它们与混凝土之间的相对滑移[16]。

a——混凝土;b——普通钢筋;c——钢衬里;d——预应力筋

1.5 安全壳损伤参数定义

安全壳结构的完整性和密闭性一直是众多学者关注的重点。当前,工程界常采用应力、应变或位移来定义安全壳的功能或结构失效,这类定义能够有效反映结构的局部开裂或损伤情况,却不能直接反映结构整体的开裂或损伤情况。在持续增长的内压作用下,安全壳各部件的强度、刚度会发生退化和损伤,结构损伤不断累积,最终产生破坏。综合来说,安全壳混凝土的损伤行为是一个不断累积的过程,在承压过程中会产生贯穿裂缝并不断扩展至整个筒壁,已经产生贯穿裂缝的区域会发生强度和刚度的退化,这会对安全壳整体损伤和失效产生影响。基于此,本研究提出了“损伤面积比”的参数来量化安全壳混凝土的损伤行为,以完整地反映安全壳在损伤演化过程中的失效情况。

(5)

式中:DRc为混凝土损伤面积比;AreaD为混凝土产生贯穿裂缝的区域面积;AreaT为混凝土的总面积。DRc的取值范围为[0,1],当DRc=0时代表安全壳混凝土尚未产生贯穿裂缝,当DRc=1时代表安全壳混凝土均已产生贯穿裂缝。

当混凝土拉应变超过峰值拉应变时,混凝土损伤超过损伤限值,即代表安全壳混凝土层发生受拉破坏,关于计算混凝土产生贯穿裂缝的区域面积在ABAQUS程序的后处理模块完成,详细计算步骤如下:1) 依据《混凝土设计规范》(GB50010—2010),C50混凝土的标准抗拉强度为2.64 MPa,弹性模量为34 500 MPa,由此可计算得到C50混凝土的峰值拉应变和相应的损伤因子;2) 在ABAQUS中的Contour Plot模块中定义该损伤因子为损伤限值,当混凝土的损伤超过损伤限值时即代表安全壳混凝土层发生受拉破坏;3) 在ABAQUS中选定混凝土损伤超过损伤限值的区域,定义为AreaD。

图3示出安全壳混凝土产生贯穿裂缝的判别方式,左侧区域混凝土损伤穿过整个混凝土层,其为混凝土贯穿裂纹区域;右侧区域混凝土损伤未穿过整个混凝土层,因此不能判别为混凝土贯穿裂纹区域。

图3 安全壳混凝土裂缝贯穿

1.6 安全壳概率有限元分析

开展安全壳概率安全评估时,有必要综合考量多种误差因素对结果造成的影响。这些误差因素主要包括安全壳几何不确定性、模型不确定性、材料不确定性以及荷载不确定性等。文献[17]认为安全壳是施工质量受到严格控制的结构,可以忽略由几何不确定性带来的影响。因此,本文在进行安全壳概率安全评估中重点考察了材料不确定性、统计不确定性、荷载不确定性以及模型不确定性的影响。安全壳各种材料的统计特性[15,18]列于表2。表2中:fc为混凝土单轴抗压强度;fy400为钢筋屈服强度;Ey400为钢筋弹性模量;fl为钢衬里屈服强度;El为钢衬里弹性模量;fp为预应力筋屈服强度;Ep为预应力筋弹性模量。利用表2给出的7个材料参数的分布形式和范围,本文采用拉丁抽样方法在Matlab中生成了100组样本,依次导入ABAQUS中运行得到计算结果。

表2 随机变量的统计特性

2 易损性分析理论

传统的易损性分析通常假定易损性函数服从对数正态分布的形式,如式(6)所示。这种假定使得进行易损性评估时仅需要估计关键的易损性参数即可,大大简化了易损性评估的计算量。然而,这种简化使得易损性评估缺乏严格的理论分析,并且不能全面考虑概率分析中存在的认知、经验、统计等相关不确定性因素,这可能会对易损性评估结果产生较大影响。

(6)

式中:F为失效概率;Φ为标准正态累积分布函数;pm为安全壳结构承压能力中值;βs为对数标准差;p为安全壳事故压力。

贝叶斯估计是一种广泛应用于结构易损性评估和可靠性分析的统计学方法。它通过建立结构的概率能力和概率需求模型综合考虑概率分析中的相关不确定性因素,从而得到更可靠的易损性分析结果。

不论安全壳材料、承受的荷载如何变化,安全壳混凝土的损伤面积比情况反映了混凝土整体的损伤情况以及可能的泄漏水平,因此本研究中假定安全壳概率能力模型是确定的损伤面积比。而安全壳概率需求模型可通过下式[11]确定:

(7)

(8)

由于安全壳在内压作用下的力学行为较为复杂,很难建立损伤与内压的确切物理力学关系。因此,确定性模型通过拟合100组样本的损伤中值与内压的数学关系得到,如图4所示。可以发现,安全壳混凝土的损伤中值与内压近似服从均值为1.351 8、方差为0.112 8的累积正态分布,如式(9)所示。

(9)

图4 混凝土损伤面积比与内压的数学关系

式中,P为内压。

在构建安全壳概率需求模型时,混凝土抗压强度(h1(x)=fc)、钢筋强度(h2(x)=fy400)、钢筋弹性模量(h3(x)=Ey400)、预应力筋强度(h4(x)=fp)、预应力筋弹性模量(h5(x)=Ep)、钢衬里强度(h6(x)=fl)、钢衬里弹性模量(h7(x)=El)及内压(h8(x)=P)被选择作为确定性需求模型的修正项。然而,修正因子过多可能会导致模型过拟合,模型变得过于复杂也会导致其无法很好地推广到新样本,因此有必要对修正项进行参数优化和校正。

模型优化和校正过程分为以下3步:1) 采用贝叶斯估计方法计算所有模型参数和标准差的后验分布;2) 对比解释函数的变异系数,去除变异系数最大的一项,该项被认为是解释函数中包含信息量最少的一项;3) 去除某项后如果后验均值未发生明显变化,则认为去除该项是合理的,重复此流程,直到后验均值明显增大,认为去除该项影响了模型的精度,则该项不能去除,因此,返回上一步,模型优化停止。

按照上述的参数优化方法对参数进行取舍,如表3所列。第1次去除了参数θ3(1.867),第2次去除了参数θ7(21.419),随后依次去除了θ5(2.044)、θ2(5.864)、θ6(1.848),当去除参数θ1(0.646)后,模型标准差显著增大,这表明去掉该项会使模型精度降低,因此,参数优化过程停止,取上一次优化的结果作为最优模型。

表3 模型参数的统计信息

经贝叶斯估计方法优化后,参数θ1、θ4、θ8最终保留为修正项,因此,安全壳混凝土的损伤-内压概率需求模型可简化为式(10)。计算模型每个参数的标准差和变异系数,以及后验分布的均值,如表4所列。

(10)

表4 模型参数的后验统计

极限状态方程定义了安全壳的安全状态与失效状态的边界,这在进行易损性评估中起着关键作用。通常,极限状态方程的形式表达为能力与需求之差,形式如下:

Zk=Ck(x)-Dk(x)

(11)

式中:Zk为安全壳混凝土的第k个失效状态的安全余量;Ck为第k个状态下的能力函数;Dk为第k个状态下的需求函数。

需求模型是包含随机变量的函数,因此具有概率分布,而能力模型被认为是安全壳特定状态下的恒定损伤值,它反映了安全壳的整体损伤与泄漏水平,因此极限状态方程也具有概率分布。当Zk<0时表明该状态下安全壳的损伤能力值小于损伤需求值,即安全壳达到该状态下的失效水平。安全壳的易损性评估可采用极限状态方程的积分形式来表达:

(12)

式中:Fk为第k个状态下的安全壳的失效概率;P为第k个状态下安全壳损伤能力值小于损伤需求值的概率;fzk(Z)为概率密度函数;Z为功能函数。

概率需求模型由多个具有不确定性因素的随机变量组成,这导致易损性评估存在不确定性。因此,有必要给出具有置信边界的易损性曲线,该曲线可以表征结构在特定损伤状态下的失效概率范围,获得更加精确的易损性评估结果。易损性曲线的置信边界可通过1阶分析法获得[10],已知失效概率和可靠性指标存在如下关系:

β(c,d,Θ)=Φ-1[1-Fk(c,d,Θ)]

(13)

式中:β为可靠性指标;c为能力模型中的变量;d为需求模型中的变量;Θ为模型参数。

可靠性指标β的方差可近似由下式获得:

(14)

置信边界可被定为{Φ[-β-σβ]、Φ[-β+σβ]},该置信边界对应于15%和85%的置信水平。

3 有限元计算结果

3.1 安全壳样本差异性分析

图5示出考虑材料不确定性情况下混凝土损伤演化的差异。根据混凝土的抽样分布情况将结果样本分为Mean-1Std、Mean和Mean+1Std 3个区域,Mean表示混凝土抗拉强度均值,Std表示混凝土抗拉强度标准差。从每个区域内随机挑选出样本,比较其损伤演化和样本间的差异性。在0.8 MPa内压下,混凝土裂缝集中在设备洞口上下两侧,此时,由材料不确定性引起的变异性并不明显。当内压增长至1.0 MPa时,选自Mean-1Std区域的安全壳样本的混凝土裂缝已经扩展至穹顶顶部,而选自Mean及Mean+1Std区域的安全壳样本尚未发生明显变化。随着内压继续增长至1.2 MPa,选自Mean-1Std区域的安全壳混凝土裂缝广泛开展至穹顶、洞口周围及底部,样本混凝土损伤已经达到18.61%,而选自Mean及Mean+1Std区域的安全壳样本损伤仅达到4.63%和1.99%。内压增长至1.4 MPa时,选自Mean-1Std区域的安全壳混凝土裂缝扩展至整个筒体,损伤程度远超过选自Mean及Mean+1Std区域的安全壳样本。当内压增长至1.6 MPa时,安全壳混凝土绝大部分开裂,样本差异性则明显减小。结果表明,由材料不确定性引起的安全壳样本损伤的差异性不可忽略,尤其当内压处于[1.2 MPa,1.6 MPa]时,这种差异性更为明显。采用混凝土损伤面积比的形式可以反映出安全壳混凝土在内压作用下的损伤演化行为,并为研究其在不同损伤水平下的失效概率奠定了基础。

图5 安全壳混凝土损伤演化

3.2 易损性曲线分析结果

图6示出安全壳在不同损伤状态下的概率密度曲线。随着混凝土损伤面积比增大,概率密度曲线的峰值逐渐向右移动,这表明样本损伤中值逐渐增加。曲线的峰值逐渐减小且曲线的轮廓逐渐变宽,这表明样本分布的标准差逐渐增大,样本的可变性逐渐增加,由非弹性行为引起的不确定性增加。

图6 不同损伤状态下的概率密度曲线

表5列出相同损伤下采用传统方式和贝叶斯估计获得的易损性参数,并绘制了相应的易损性曲线,如图7所示。可以看出,在相同损伤下,采用贝叶斯估计和传统方式获得的易损性参数中值相差不大,采用贝叶斯估计获得的易损性参数变异系数大于传统方式,使得贝叶斯易损性曲线相对平缓。这是由于贝叶斯估计在构建易损性曲线时经过了严格的理论分析,并且综合考虑了各项材料参数对需求模型的影响,最终获得的易损性曲线具有较高的准确度。

表5 传统方式和贝叶斯估计易损性参数对比

DRc:a——0.1;b——0.2;c——0.3;d——0.4;e——0.5;f——0.6;g——0.7;h——0.8;i——0.9

图8示出不同损伤状态下的贝叶斯易损性曲线及相应的置信边界。随着混凝土损伤面积比增大,易损性曲线的能力中值逐渐增大,这表明在内压一定时,安全壳越难以发生损伤程度严重的失效。如在内压增长至1.2 MPa时,安全壳混凝土有49.7%的可能发生损伤程度为0.1的失效,而仅有7.4%的可能发生损伤程度为0.5的失效。通过量化不同损伤状态下安全壳的失效概率有利于掌握安全壳在内压增长过程中的损伤及失效情况,这有利于监管者针对不同的损伤情况做出相应的防护措施。

DRc:a——0.1;b——0.2;c——0.3;d——0.4;e——0.5;f——0.6;g——0.7;h——0.8;i——0.9

4 安全壳概率安全评价

概率安全评估被认为是核工程界评估易损性的重要标准。通常情况下,安全壳结构的事故压力遵循对数正态分布,中值为0.663 MPa,对数标准差为0.3[19]。安全壳结构的总失效概率(CCFP)可采用下式进行计算:

(15)

式中:F(fail|p)为安全壳结构在特定损伤状态下的易损性曲线;fp(p)为事故压力的概率密度函数。根据式(15)可计算得到安全壳在任一损伤状态下的总失效概率,表6列出采用传统易损性评估方法和贝叶斯估计方法得到的安全壳在不同损伤状态下的总失效概率的相关信息。表6中:μCCFP为总失效概率的均值;σCCFP为总失效概率的标准差;δCCFP为总失效概率的变异系数。

表6 安全壳结构的总失效概率的统计信息

由表6可见,总体而言安全壳的总失效概率随着混凝土损伤面积比的增大而减小,且安全壳在任一损伤状态下的总失效概率均小于0.1[20],这表明本文分析的安全壳结构在损伤演化的全过程中均满足严重事故下的性能要求。不同之处在于,采用贝叶斯估计方法获得的总失效概率的均值大于传统易损性评估方法,而变异系数小于传统易损性评估方法。这表明传统易损性评估可能低估了安全壳在损伤状态下的总失效概率,采用贝叶斯估计方法可以取得更加保守的概率安全评估结果。