序列型地震作用下AP1000核反应堆屏蔽厂房损伤灾变精细化分析

汪大洋， 包嗣海， 陈婉若， 朱 勇， 张永山

(1. 广州大学 土木工程学院， 广州 510006； 2. 广州大学 结构力学分析与测试研究中心， 广州 510006)

近年来核电站事故频发，对人民生命财产造成了严重的威胁。2007年日本新潟县发生里氏6.8级地震，超过柏崎·刈羽核电厂7台机组的抗震设计基准；2011年日本发生9.0级大地震，并继发海啸，导致福岛第一、第二核电站严重受损，如图1所示。鉴于核电结构遭遇地震侵袭的事实可能性及其导致的严重后果，越来越多的学者开始关注地震作用下核电站结构的破坏机理和抗震性能。Kussuga等[1]利用有限元软件ANSYS建立核岛结构有限元模型，针对水平地震动作用下对核岛结构开展非线性动力时程分析；Daogang等[2]对水箱不同水位时结构在地震荷载作用下的加速度反应、最大等效应力分布开展研究，结果表明结构的最大等效应力出现在屏蔽结构的屋顶与重力水箱的相交处。我国核电站结构抗震性能研究虽起步较晚，但也取得的很多具有指导意义的研究成果。李小军等[3]以CAP1400型核岛结构为研究对象，开展非基岩场地下核岛结构振动台试验，研究在非基岩场地条件下核岛结构的适用性和动力特性。钱稼茹等[4-5]针对我国CNP1000核电厂安全壳，完成了一个缩尺比为1∶10的预应力混凝土模型，以此开展单自由度拟动力试验，试验结果表明，在峰值加速度为3g时，地震作用下安全壳的筒体底部出现少量裂缝，结构处于弹性阶段，仍可以按线弹性结构进行地震反应计算。李建波等[6]开展核电厂的楼层反应谱分析在Biggs和Kapur方法的基础上提出一种修正模态法，既考虑了结构-地基动力相互作用的影响,又计算简便易行，也可获得比较满意的精确度。Wang等[7-8]基于核岛有限元数值模型开展了强主余震对结构塑性损伤行为的调查研究，结果表明屏蔽厂房和附属厂房在强烈地震中损伤效应严重，钢安全壳结构基本完好无损，余震的影响强烈余震引起的作用可能导致结构构件的严重损坏，余震的幅度越大加剧效应越明显。

图1 地震海啸引发福岛核事故Fig.1 Earthquake tsunamiinduced a nuclear accident

可以看到，现有研究多从整体表征角度探讨核电结构的抗震性能，实际上核电结构体型庞大且复杂，从整体角度表征的损伤发展并不能精确反映关键控制部位的实际损伤状态，有必要深入探究不同部位、尤其是损伤集中局部区域在强震及其主余震作用下的塑性损伤发展过程。为此，本文以AP1000核反应堆屏蔽厂房结构为研究对象，建立三维弹塑性损伤有限元计算模型，在主震和序列型地震动整体损伤调查的基础上进一步探讨区域局部损伤力学特征，揭示损伤集中区域关键部位的损伤发展，建立有效的精细化损伤评价技术方法，为核反应结构抗震设计和后续研究提供设计指导和参考。

1 屏蔽厂房结构有限元模型建立与验证

1.1 AP1000屏蔽厂房结构概况

AP1000核反应堆屏蔽厂房直径为44.2 m，高度为83.37 m，壁厚为0.914 m，屏蔽厂房上部为冷却水箱，其外半径和内半径分别为11.75 m和5.334 m，水箱壁厚为0.6 m，通风口高度为1.45 m。屏蔽厂房上各有两个圆形设备闸门和两个圆形人员闸门，设备闸门的内径为4.88 m，人员闸门的内径为3.05 m。在标高25.01 m处安装一个设备闸门和一个人员闸门；在标高35.44 m处安装另一个设备闸门和另一个人员闸门。设备闸门和人员闸门的主要贯穿方位设计，如图2所示。结构内部分别由重力水箱结构、斜面结构、等效钢梁和等效设备层组成，如图3所示。

图2 闸门贯穿方位图Fig.2 Gate penetration diagram

图3 AP1000核反应堆屏蔽厂房结构剖面图Fig.3 AP1000 nuclear reactor shielded plant structure section

1.2 塑性损伤有限元模型建立

1.2.1 材料本构关系模型

AP1000核反应堆屏蔽厂房为钢筋混凝土结构，混凝土压缩行为采用Popovics[9]和Yip[10]提出的本构关系，混凝土受压应力应变关系如下

(1)

混凝土受拉应力应变关系如下

(2)

式中：fct是混凝土抗拉强度；εtr是与抗拉强度fct相对应的混凝土峰值拉应变；σt和εt分别是混凝土拉应力和混凝土拉应变；α是衰变因子。具体关系式如下

α=0.017+0.255(nρ)-0.106(nρ)2+

0.0160(nρ)3

(3)

式中：n=ES/EC，ES是钢筋的弹性模量；ρ是对应的配筋率。考虑混凝土模量折减的塑性损伤模型公式

E=(1-d)Ec

(4)

式中：E是折减的弹性模量；d是损伤因子。通过等效余能理论计算[11]，假设损伤的余能W0等于未损伤的余能Wd，如下式所示

W0=Wd

(5)

(6)

(7)

结合式(5)～(7)，推导出Ed和Ec的关系式如下

Ed=Ec(1-d)2

(8)

由于混凝土的切线模量等于应力应变之比，将其代入式(8)中，可推导出下式

(9)

因此，可以推出混凝土压缩损伤因子dt和拉伸损伤因子dc的计算公式，具体如下

(10)

钢筋本构采用理想弹塑性模型，设计参数参考GB 50010—2010《混凝土结构设计规范》[12]。屏蔽厂房的混凝土抗压强度为27.6 MPa，拉伸强度为2.125 MPa，弹性模量为26 267 MPa，屏蔽厂房中的钢筋采用HRB400E型钢筋，配筋率为0.003 06。混凝土在受拉和受压两种情况下的应力和应变值及损伤因子如图4所示。

(a) 拉伸本构

(b) 压缩本构

(d) 压缩损伤因子图4 混凝土拉压本构与损伤因子Fig.4 Tension and compression constitutive laws and the corresponding damage factors

1.2.2 边界条件与网格划分

通过ABAQUS有限元软件建立AP1000核反应堆屏蔽厂房模型(见图5)，采用四节点曲壳单元(S4R)单元进行建模，S4R是一种通用壳单元，由四个具有弯曲和膜能力的节点组成，可用于薄壳或厚壳结构建模，采用减缩积分方式，包含沙漏模式控制，容许有限薄膜应变。相对于实体单元建模，采用S4R分析效率更高。考虑到屏蔽厂房和安全壳有圆形开洞的人员闸门和设备闸门，如果采用四边形的网格划分可能会产生畸形网格，进而导致模型计算不收敛，因此配合使用三角形网格和四边形网格对屏蔽厂房和安全壳进行分区划分。在圆形开洞处，采用1.0 m的三角形网格加密划分，屏蔽厂房竖直筒壁的其他区域采用边长为1.0 m×1.0 m的四边形网格划分；由于重力水箱内壁直径、外壁直径和基础底面直径均不相同，为了提高网格划分质量，对水箱顶面和斜面结构采用按扇形展开的四边形网格划分，水箱内壁四边形网格尺寸为1.0 m×0.24 m，水箱外壁网格尺寸为1.0 m×0.59 m，水箱顶面网格尺寸由1.0 m×0.24 m至1.0 m×0.59 m逐渐变大，斜坡屋顶网格尺寸由1.0 m×0.24 m至1.0 m×1.0 m逐渐变大，其他区域均采用1.0 m×1.0 m的四边形网格划分，得到AP1000核岛结构ABAQUS有限元模型的网格划分单元数为45 039个，节点总数为38 618个。

1.3 模态分析与验证

采用ABAQUS软件Lanczos法进行模态分析，其内部自动采用稀疏矩阵直接求解器，适合大型结构。表1给出了所建屏蔽厂房结构前八阶模态计算结果，前八阶振型如图6和表1所示。可知，结构第一阶、第二阶自振频率分别为3.609 7 Hz、3.611 6 Hz，相差较小，说明结构两个水平方向的抗侧刚度较为接近，属于对称结构体系。进一步将本文计算结果与既有文献结果进行对比分析(见表2)，可以看出模型第一阶振型均为梁式平动，因参考图纸的细部差异，导致本文基频计算结果较现有文献偏小，一定程度上降到了反应堆结构的抗侧刚度，进而导致基频偏小，为确保计算结果进一步接近实际情况，应考虑闸门洞口效应。总体而言，本文计算模型具有较高的可信度，与既有文献基频最大相差10.73%、最小仅为0.25%。

(a) 整体模型

(b) 剖面图图5 设置闸门AP1000核岛结构精细化有限元模型Fig.5 AP1000 NPP structure with gate holes

表1 AP1000模型模态分析

(a) 第1阶

(b) 第2阶

(c) 第3阶

(d) 第4阶

(e) 第5阶

(f) 第6阶

(g) 第7阶

(h) 第8阶图6 AP1000核反应堆屏蔽厂房结构前8阶振型图Fig.6 AP1000 nuclear reactor shield building of the first 8 modes

表2 核反应堆模型第一振型对比Tab.2 Comparison the first mode of the nuclear reactor model

2 地震动输入和工况设置

基于美国太平洋地震工程研究中心地震波数据库PEER (Pacific Earthquake Engineering ReseaMSh Center)选取地震动，考虑三向地震输入，其中水平方向与竖向加速度峰值比设为H1∶H2∶V=1∶1∶0.67(H1、H2、V分别代表X、Z水平方向、Y竖向)。核岛结构一般建在土体剪切波速大于700 m/s的岩土场地上，选取的3条地震波剪切波速分别为804 mm/s、829 mm/s和891 mm/s，满足核电站场地土剪切波速要求。各地震波的信息见表3，图7为所选地震波反应谱与RG1.60[15]设计反应谱的对比，可见所选地震波与RG1.60谱具有吻合性较好。

表3 地震波信息

(a) EA水平向

(b) EB波水平向

(c) EC波水平向

(d) EA波竖向

(f) EC波竖向

根据我国GB 50267—1997《核电厂抗震设计规范》[14]中规定极限安全地震动不大于0.5g的地区，核电站设计基准地震(极限安全地震)的加速度峰值不得低于0.15g，法国将安全设计地震的加速度峰值为0.2g，以0.15g～0.25g作为设计地震动的加速度峰值，而从美国西屋公司引进的AP1000第三代核电机组的抗震设计标准为0.3g，在我国台湾地区7台核电机组的设计地震动峰值加速度均超过0.3g，为此将0.3g作为本文的设计地震动输入。此外，考虑到近年来频发的罕遇大地震，以及多次余震作用，如：汶川5·12地震里氏震级8.0级、地震烈度11度，发生27次5级以上地震；日本3·11里氏震级9.0级，地震发生后10个小时内共发生13次5级以上地震。在设计地震动调查的基础上进一步研究超设计地震动作用，探究结构损伤发展过程。将所选的3条地震动的分析工况按PGA分别为0.3g-2.1g进行设置，PGA间隔为0.2g，单一主震和主余震各30个工况，共60个工况，见表4。采用重复法构造主余震，主余震之间间隔60 s，主余震峰值比为1∶0.852 6[16]，代表性主余震工况AA8如图8所示。

表4 单一主震和主余震下分析工况表

图8 AA8主余震序列Fig.8 AA8 main-aftershock sequence

3 屏蔽厂房整体动力灾变分析

钢筋混凝土屏蔽厂房结构在地震作用下拉伸塑性损伤效应明显，为探究其灾变过程，提取拉伸损伤因子指标，根据1.2.1节的损伤因子计算原理进行计算，反映混凝土材料在逐渐增大的地震动输入下的拉伸损伤累积变化过程。图9给出了工况EA1-EA10、EB1-EB10、EAA1-EAA10、EBB1-EBB10拉伸损伤因子云图。由图可知，在设计地震动作用下(0.3g)结构损伤效应不明显，但随着地震动强度的增加，灾变损伤效应逐渐加剧，损伤区域不断增大；损伤首先出现在闸门洞口处，因刚度突变导致应力集中而形成损伤，随后逐步加剧并向四周扩散，集中分布在高度为22～39 m区域范围内；上部洞口和斜面区域亦随着地震动PGA的增大，动力损伤灾变效应逐渐加剧，为结构第二损伤发展区域，同样因刚度突变和应力集中导致；其余部位，如底部基础区域、中部无洞口区域以及顶部水箱区域，随着输入PGA的增大虽呈现不断增大的损伤发展，但其损伤程度要远远小于闸门洞口区域和上部洞口与斜面区域。

图9 单一地震作用下屏蔽厂房累积拉伸损伤因子云图Fig.9 Cumulative tension damage factor cloud of shielded factory building under single earthquake

在整体结构未损伤或损伤不明显的工况下，余震对经历主震后的结构损伤加剧效应较小，如图10所示(以损伤最为严重的22～39 m高度范围为例)，可见工况EB3和EBB3(PGA为0.7g)以下，二者损伤程度差别不大，余震后损伤区域扩展有限。然而，当输入工况PGA达到0.9g以上时，余震对经历主震后的结构损伤加剧效应，不论是损伤因子大小还是损伤区域发展都大大增加，如图11工况EB6和EBB6。究其原因，结构在地震作用下本质上是能量的传递、耗散和吸收的过程，以损伤耗能为例，核岛结构的损伤耗能可表示为[18-19]

图10 单一地震和序列型地震作用下屏蔽厂房拉伸损伤因子云图Fig.10 Tension damage factor cloud of shield building under single earthquake and sequence earthquake

(a) EA1&EAA1

(b) EA2&EAA2

(c) EA3&EAA3

(d) EA4&EAA4

(e) EA5&EAA5

(f) EA6&EAA6

(g) EA7&EAA7

(h) EA8&EAA8图11 屏蔽厂房结构地震损伤耗能Fig.11 Earthquake damage energy consumption of shield building

虽然随着地震动强度增加单一地震和主余震下屏蔽厂房损伤耗能都随之增加，但增量有所不同，如工况EAA1、EAA4、EAA7、EAA8相对于对应的工况EA1、EA4、EA7、EA8，其损伤耗能增量依次为0.05 MJ、1.09 MJ、2.61 MJ、3.21 MJ，耗能增量呈现逐渐趋势增大，直接导致结构损伤灾变效应不断加剧。

此外，结构灾变损伤在闸门洞口、通风洞口与斜面处集中发展，呈现非常不规则的损伤分布，如图9工况EA6上部设备洞口附近的损伤因子达到0.86，而距离该洞口上部5 m位置的损伤因子仅为0.81。同时需要看到，屏蔽厂房结构随着地震PGA的增长，灾变过程由闸门洞口开始，首先向环带扩展，然后逐步向下部延伸，直至基础平台顶板。可见，闸门洞口损伤发展直接决定了屏蔽厂房地震灾变过程，因而有必要探究损伤集中区域的灾变过程，有助于掌握屏蔽厂房在不同等级地震激励下结构性态。

4 屏蔽厂房区域动力灾变分析

4.1 核岛结构区域损伤划分

为探究结构关键区域灾变损伤，基于分区思想对结构区域进行划分。沿高度分为基础区域(0～10 m，简称BS区)和主体区域(10～65.6 m，简称MS区)，同时考虑沿高度方向损伤差异较大，进一步细分为BS-1、BS-2、MS-1～MS-7区，如图12(a)所示；考虑沿环向的损伤亦有较大差异，环向每隔15°划分为一个区域，如图12(b)所示。由前述分析，屏蔽厂房人员和设备闸门处存在应力集中现象，设备闸门对应编号HB1、HB2，人员闸门分别对应HS1、HS2，则闸门HB1、HS1分别在MS-2区域的20和17区块中；闸门HB2、HS2分别在MS-4区域的16和17区块中。

(a) 竖向分区

(b) 环向分区图12 屏蔽厂房区域划分示意图Fig.12 Division of the shield building

4.2 区域灾变损伤分析

图13和图14分别显示了单一和序列型地震作用下屏蔽厂房区域拉伸损伤因子平均值(Mean)和最大值(Max)的对比变化梯度图，图15给出了单一地震作用下BS和MS区域平均和最大拉伸损伤因子沿环向变化图，区域平均值为所在区域所有提取点拉伸损伤因子的平均值，区域最大值为所在区域所有提取点拉伸损伤因子的峰值。由图可见，屏蔽厂房随着地震动强度的增大，损伤部位主要集中在闸门与上部洞口附近及其所处环带上。当地震动PGA较小时，拉伸损伤因子沿环向呈现明显的梯度变化(如图13(a)、15(a)所示)，说明材料损伤首先发生在刚度突变处引起的应力集中效应处，如图15(a)中EA2工况拉伸损伤因子为0.340(C13区)，而该区域环带最小值仅为0.074(C18区)，二者相差4.6倍；然而，随着PGA不断增大，损伤不断沿环向扩展，最终贯通整个环带，平均拉伸损伤因子在同一高度的环带上逐渐趋于均匀(如图13(d)、15(b)所示)，如图15(b)中EA8工况拉伸损伤因子最大值为0.803(C10区)，而该区域环带最小值仅为0.738(C18区)，二者相当。

(a) EA1(Mean)

(b) EA3(Mean)

(c) EA5(Mean)

(d) EA8(Mean)

(e) EA1(Max)

(f) EA3(Max)

(g) EA5(Max)

(h) EA8(Max)图13 单一地震作用下屏蔽厂房区域拉伸损伤因子平均值(Mean)和最大值(Max)对比Fig.13 Comparison of mean (Mean) and maximum (Max) of tension damage factor in shield building area under single earthquake

(a) EAA1(Mean)

(b) EAA3(Mean)

(c) EAA5(Mean)

(d) EAA8(Mean)

(e) EAA1(Max)

(f) EAA3(Max)

(g) EAA5(Max)

(h) EAA8(Max)图14 序列型地震作用下屏蔽厂房区域拉伸损伤因子平均值(Mean)和最大值(Max)对比Fig.14 Comparison of mean and maximum of tension damage factor in shield building area under sequence earthquake

图13(e)～(h)和图14(e)～(h)为序列型地震作用下屏蔽厂房区域拉伸损伤因子的平均值和最大值，可以看出其与同一PGA输入下的损伤梯度图相比，主要差别在于洞口附近的损伤进一步加剧，如图13(c)、14(c)所示，说明余震对于主震后结构损伤的加剧效应仍集中在已损伤区。为了探究余震对结构的损伤加剧影响，图16分别给出了MS-2、MS-4中HB1和HB2两个洞口的损伤平均值和损伤最大值的损伤因子比曲线(主余震下损伤因子/单一主震下损伤因子)。由图16可知，HB1、HB2的损伤因子比均大于1，表明余震加剧了洞口处结构的损伤破坏程度，随着地震动PGA增大，损伤因子之比整体呈现先上升后下降的趋势，当地震动PGA超过1.5g时损伤因子之比趋于1，说明结构在单一主震下已发生破坏，基本无加剧效应。对比损伤平均值和损伤最大值的损伤因子比，可以看出在地震动较小时，损伤最大值的损伤因子比高于损伤平均值的损伤因子比，而随着地震动强度的增加，损伤因子比趋向于1。

(a) BS-1(Mean)

(b) BS-2(Mean)

(d) MS-4(Mean)

(e) MS-2(Max)

(f) MS-4(Max)图15 单一地震作用下BS和MS区域平均和最大拉伸损伤因子沿环向变化图Fig.15 Changes of the average and maximum tension damage factor in the BS and MS regions along hoop direction under single earthquake

(a) HB1

(b) HB2图16 损伤因子比Fig.16 Value of damage factor ratio

对比区域拉伸损伤因子平均值和最大值，可以看出在洞口高度所处的环带上，最大值与平均值差距明显，如图14(c)、(g)所示，在22～39 m、60～65.6 m高度范围内，损伤因子云图存在明显的梯度区别，且梯度变化主要集中在洞口附近；而在其他高度范围内，拉伸损伤因子的平均值和最大值相差不大。因此，屏蔽厂房结构在超设计地震激励下除各类洞口及其所在环带存在明显的应力集中与塑性损伤外，其余区域的动力损伤灾变过程在同一高度环带上趋于均匀，说明屏蔽厂房结构在设计、超设计、甚至巨震作用下的损伤变形特征均以剪切型破坏为主，动力灾变过程为洞口刚度突变处首先由于应力集中效应导致损伤失效，随后向洞口所在高度的环带方向拓展、贯穿，最后沿高度向上下延伸，损伤程度以闸门洞口环带最为严重、其次为上部洞口环带、再次为基础和上下洞口之间区域、最后为顶部水箱区域。

5 闸门洞口损伤发展分析

由上述分析可知，屏蔽厂房损伤发展在闸门和上部洞口区域最为集中，并以闸门洞口最显著，本节将进一步探讨设备闸门HB1、HB2和人员闸门HS1、HS2的损伤影响范围及其随PGA输入大小的变化规律。表5给出了闸门洞口0.3 m、0.5 m、1 m半径区域范围内(如图12(a)所示)的拉伸损伤因子衰减率δ的对比情况(δ=(Max-Mean)/Max×100%)。

表5 单一地震作用下闸门洞口1 m范围内拉伸损伤因子

由表5可以，在相同PGA输入下，随着洞口半径方向区域的逐渐增大，拉伸损伤因子的衰减率大幅提升，说明拉伸损伤因子随半径增大呈快速下降趋势，进一步证明闸门洞口为损伤最严重区域，如PGA为0.5g时，设备闸门HB1在0.3 m、0.5 m、1 m半径范围内拉伸损伤因子平均值依次为0.742、0.664、0.544。在相同半径区域范围内，随着PGA输入量级增大，拉伸损伤因子的衰减率不断降低，说明拉伸损伤因子随PGA增大呈快速上升趋势，亦直接证明结构损伤灾变随着PGA输入的增大而不断加剧，如设备闸门HB1在0.3 m半径范围内拉伸损伤因子平均值从0.3g～1.1g依次为0.608、0.742、0.845、0.889、0.933。可以看出，当输入PGA低于0.5g时，闸门洞口对0.5 m范围内灾变损伤的影响最为严重；当输入PGA超过0.9g时，闸门洞口对周围材料损伤严重影响的区域范围达到2.0 m。因此，当采用平均值描述闸门所在高度环带的损伤时，宜将该区域范围内的损伤情况单独考虑，同时建议在设计中有针对性地对闸门洞口采取一定的加强处理措施，以提升屏蔽厂房在超设计地震作用下抵抗损伤灾变的能力。环向0.3 m平均值、环向0.5 m平均值、环向1 m平均值分别为0.804、0.742、0.664、0.544，对应的差率分别为8%、21%、48%。

6 结 论

本文根据不同地震动下拉压损伤分布情况对核岛结构进行区域划分并对各区域进行损伤分析，得到以下结论：

(1) 屏蔽厂房动力灾变效应随PGA不断加剧，损伤带首先出现在闸门洞口处，随后逐步加剧并向环带扩散，集中分布在高度为22～39 m区域范围内；上部洞口和斜面区域动力损伤灾变效应亦逐渐加剧，为第二损伤带；底部基础区域、中部无洞口区域以及顶部水箱区域虽呈现不断增大的损伤趋势，但损伤程度远小于洞口区域和上部斜面区域。

(2) 余震对主震后结构损伤的加剧效应仍集中在已损伤区(闸门和通风洞口、上部斜面区)，随地震动PGA增大加剧效应呈现先增长后下降的趋势，当PGA超过1.5g时，主震下结构损伤已贯穿整个环带，此时余震进一步作用对结构的加剧效应不明显。

(3) 屏蔽厂房结构在设计基准地震、超设计基准地震、甚至巨震作用下的损伤变形特征均以剪切型破坏为主，动力灾变过程为洞口刚度突变处首先由于应力集中效应导致损伤失效，随后向洞口所在高度的环带方向拓展、贯穿，最后沿高度向上下延伸。

(4) 区域损伤因子平均值和最大值在洞口所处环带上差距明显，梯度变化主要集中在洞口附近，尤以22～39 m、60～65.6 m高度范围最强，而在其他高度区域，损伤因子的平均值和最大值相差不大。

(5) 应对闸门洞口一定半径范围内的区域动力灾变进行单独损伤评估，并宜采用最大值描述其损伤行为，其余部位可采用损伤平均值描述。当PGA低于0.5g时，闸门洞口对0.5 m半径范围内灾变损伤的影响最为严重；当输入PGA超过0.9g时，闸门洞口对周围材料损伤严重影响的半径区域可达2.0 m。在设计中应有针对性地对闸门洞口采取加强处理措施。