核能不仅单位能量大，而且资源丰富，是公认的唯一现实的可大规模替代常规能源的清洁、经济能源。目前我国在建核电机组共29台，在建规模世界第一^[1]。预应力筋混凝土安全壳(简称安全壳)是核反应堆的最后一道屏障，其主要用途是屏蔽发生基准失水事故时产生的辐射物质，防止外物的袭击等，应该具有承受内压而不出现混凝土裂缝的能力^[2]。日本福岛核电厂事故表明，核电厂安全壳承受事故内压的能力十分重要，因此，对核电厂安全壳的研究就显得尤为迫切。

夏祖讽等^[3]采用杆单元模拟预应力筋对先进核电厂安全壳进行了有限元分析，得到了合理的安全壳预应力张拉顺序以及安全壳在设计事故内压、严重事故内压状态下的工作性能和极限承载能力，并与1:10的大比尺结构模型试验结果进行了比较。陈勤等^[4]则采用空间曲面膜单元来模拟预应力筋。Yonezawa等^[5]也采用杆单元模拟预应力筋，建立了2个有限元模型：360°壳单元模型和轴对称模型，通过调整预应力筋和混凝土之间的摩擦影响，分析研究了安全壳在内压作用下的极限承载力。在这些分析中，对安全壳中预应力筋的模拟都采用了简化的方式，在安全壳局部应力方面分析的不够。本文将以中国改进型压水堆核电厂(CPR1000)^[6]实际安全壳为研究对象，严格按照安全壳中预应力钢筋的布置方式，采用更合理的方法模拟预应力钢筋与混凝土之间的相互作用；通过详细的三维有限元计算，给出不同内压水平下预应力筋的内力分布、安全壳的薄弱部位、安全壳的开裂发展过程等。期望本文的研究成果对我国CPR1000核电厂安全壳的设计提供参考。

CPR1000安全壳结构由底板、筒体和穹顶组成，如图 1所示。安全壳内径37m，筒体部分高50.11m，底板底面至穹顶的总高66.68m，筒体的正常壁厚0.9m，穹顶的正常厚度0.8m，在标高22.9m方位角0°位置设有一直径7.4m的设备孔。筒体外侧设置了4个互成90°的扶壁柱。筒体部分设有2层环向及一层竖向预应力钢束，竖向钢束的上端锚固在环梁的顶面，下端则锚固在安全壳底板上。穹顶部分则布置了3层预应力钢束，钢束锚固在环梁上。预应力钢束采用1860级预应力钢绞线，采用后张法施加预应力。竖向预应力钢束横截面积为5400 mm²，数量为144根；环向预应力钢束横截面积为2850mm²，数量为223根；穹顶预应力钢束与筒体的环向预应力钢束相同，数量为174根^[7]。安全壳混凝土内还配置有HRB400普通钢筋，筒体普通钢筋的竖向体积配筋率为1%，环向配筋率为2%，穹顶互成90°角的配筋率各为2%。在安全壳内部还有一层6mm厚的钢衬里，钢衬里采用铆钉与混凝土壳紧密连接。

图 1 安全壳结构组成 Fig. 1 Composition of containment structure

图选项

本文采用ANSYS软件对CPR1000安全壳进行非线性有限元分析。安全壳有限元网格如图 2所示。安全壳的底板、筒体和穹顶采用3D加筋混凝土实体单元(SOLID65)模拟，预应力筋采用只受拉的杆单元(LINK180)模拟，钢衬里采用壳单元(SHELL181)模拟。

图 2 安全壳计算模型 Fig. 2 Numerical model

图选项

混凝土与预应力筋之间相互作用的模拟方法大致有2种：一是等效荷载法，二是实体力筋法^[8]。等效荷载法虽然网格划分简单，不必考虑预应力筋的位置直接建模，但张拉过程难以模拟，无法模拟应力损失引起力筋各处应力不同的情况。当预应力筋布置复杂时，很难用合适的外荷载代替预应力的作用，尤其是无法真正模拟预应力筋与混凝土结构之间的相互作用。实体力筋法是用杆单元在预应力筋的实际位置模拟预应力筋，预应力的模拟可采用降温法和初应变法。实体力筋法可消除等效荷载法的缺点，对预应力混凝土结构的应力分析能够比较精确地模拟。该方法在有限元建模时有三种处理方法，即实体切分法、节点耦合法和约束方程法^[9]。

1) 实体切分法：基本思路是先以混凝土结构的几何尺寸创建实体模型，然后用工作平面和预应力筋线拖拉形成的面切分实体，用切分后体上与预应力筋线相同位置的线定义预应力筋。这种方法是基于几何模型的处理，力筋位置精确，求解结果精确，但当预应力筋线型复杂时，建模麻烦。

2) 节点耦合法：基本思路是分别建立实体和预应力筋的几何模型，各自独立划分单元，之后采用耦合节点自由度将预应力筋单元和实体单元联系起来。采用此种方法需要将混凝土单元划分的很密，这样一来，计算代价将会很高，计算效率很低。由于安全壳中预应力筋较多且密，建立准确的自由度耦合非常不易。

3) 约束方程法：类似于节点耦合法，只是在联系混凝土单元和预应力筋单元的时候是选择预应力筋的节点，通过CEINTF命令自动选择混凝土单元的数个节点与预应力筋的一个节点建立约束方程，通过建立多组约束方程，将预应力筋单元和混凝土单元联系起来。由于安全壳中预应力筋很多(541根)，更由于安全壳形状都为曲面，模拟预应力筋的杆单元的数量非常大，要保证每个预应力筋的节点与混凝土相应的节点之间建立准确的约束方程几乎是不可能的。

为了能较好模拟预应力筋与混凝土之间的相互作用，保证所建立的模型能准确反映工程实际，在忽略混凝土和预应力筋之间的粘结滑移作用情况下，本文提出了一种建立安全壳有限元模型的方法。首先采用杆单元建立预应力筋的有限元模型，然后沿安全壳筒壁和穹顶的厚度方向复制预应力筋节点，在复制的层与层节点之间建立混凝土单元。此种方法保证了力筋节点和混凝土节点之间的共同作用，避免了混凝土单元细化的要求，很好地模拟了预应力筋和混凝土之间的相互作用。

普通钢筋的模拟采用整体式模型，也称分布式模型或弥散式模型，即将钢筋连续均匀分布于实体单元中，整体式模型综合了钢筋对刚度的贡献。对SOLID65单元来说，通过设定2个正交方向的体积配筋率来考虑钢筋的贡献。

安全壳底板底部节点的自由度全部约束，在模拟钢衬里的SHELL181单元上施加均布的内压，内压共分十级，每级0.1MPa。

由于安全壳结构复杂，预应力筋很多，且布置复杂，本文采用较多的单元类型和数量来建立一个比较符合实际的安全壳模型，具体计算规模见表 1。

混凝土强度等级为C50，计算中混凝土弹性模量E_c=34.5GPa，泊松比υ_c=0.2，密度ρ_c=2500kg/m³。混凝土采用多线性随动强化模型，应力应变关系如图 3所示。混凝土的单轴应力应变关系上升段采用GB50010-2010^[10]规定的公式，下降段则采用Hongnestad^[11]的处理办法，即：

图 3 混凝土应力应变关系 Fig. 3 Stress-strain relationship of concrete

图选项

其中f_c为单轴抗压强度，取为40MPa，ε₀和ε_cu为峰值应变和极限压应变，分别取为0.002和0.0033，n为模型指数，取为2。 混凝土破坏准则采用Willam-Warnke五参数模型，单轴抗拉强度f_t=2.85MPa，张开裂缝的剪力传递系数β_t=0.5，闭合裂缝的剪力传递系数β_c=1.0。SOLID65单元是通过修正应力应变曲线关系来模拟混凝土开裂，引入垂直于裂缝表面上的一个缺陷平面来表示在某个积分点上出现了裂缝，如图 4所示，其中ε^ck为f_t对应的主应力方向应变，T_c为拉应力释放因子，取缺省值0.6，R^t为图示直线的斜率，将随着求解的收敛而自适应下降^[12]。

图 4 Rt示意图 Fig. 4 Schematic diagram of Rt

图选项

预应力筋、普通钢筋和钢衬里的应力应变关系均采用理想弹塑性模型，其材料参数见表 2。

整个计算模拟了安全壳的施工过程，分为3步：1)计算安全壳在自重作用下的响应；2)施加预应力；3)施加内压。

预应力筋与混凝土之间相互作用的模拟方法是否合理可以从预应力筋单元的内力上反映出来，如果每根预应力筋杆单元的内力基本一致，表明模拟方法是准确的。各类型所有预应力筋单元在仅有预应力、预应力+0.4MPa内压(设计内压值)、预应力+0.6MPA内压、预应力+0.8MPa内压4种工况下的内力沿环向的分布情况如图 5所示。由图 5可以看出，各预应力筋单元的内力分布基本一致。各类型预应力筋单元内力的最大值、最小值、均值以及相对偏差(最大值与最小值之差)如表 3所示，最大相对偏差基本都在6%以内。最大相对偏差11.4%发生在竖向预应力筋单元上，主要发生在设备孔周围(如图 5所示)，这是由于这里的预应力筋弯曲比较厉害造成的。计算结果表明本文采用的预应力筋与混凝土相互作用的模拟方法是可行且合理的。随着内压的施加，预应力筋单元的内力逐步增加，表明预应力筋对安全壳承受内压起到了很大的作用。

图 5 预应力筋单元内力分布 Fig. 5 The distribution of axial forces of the tendon elements

图选项

图 6给出了内压作用下安全壳最大主应力沿360°展开的分布规律。在0.4MPa内压作用下，拉应力较大的区域集中在设备孔周围、扶壁柱、环梁附近和筒体底部的位置，最大拉应力出现在设备孔上下两侧，但小于混凝土的抗拉强度。穹顶最大主应力分布比较均匀。当内压增大到0.5MPa时，设备孔上下两侧的拉应力进一步增大，达到了抗拉强度，此时会有裂缝的产生，见图 7(a)。当内压达到0.8MPa时，设备孔几乎被高拉应力区所环绕。穹顶上拉应力区集中程度更加明显了，这些高应力区域是穹顶预应力筋比较少的位置，如图 2(d)所示区域。

图 6 内压作用下筒体最大主应力云图 Fig. 6 Maximum stress contours of cylinder

图选项

在各级内压作用下安全壳混凝土的开裂情况如图 7 所示。在内压为0.4MPa(设计内压值)时，安全壳混凝土没有开裂发生；当内压达到0.5MPa时，设备孔上下侧出现局部开裂，范围很小；随着内压增大，设备孔上下侧的裂缝进一步扩展；当内压达到0.8MPa时，设备孔左右两侧也出现了裂缝，上下侧裂缝进一步扩展，开裂严重，安全壳混凝土失去了它的功能。

图 7 安全壳混凝土开裂 Fig. 7 Concrete crack in containment

图选项

安全壳内部设有6mm厚的钢衬里，其作用是在安全壳发生事故时防止辐射物质的扩散，即保证安全壳的密闭性。图 8给出了各级内压作用下钢衬里第一主应力沿360度展开的分布情况。在仅有预应力作用下，钢衬里第一主应力分布均匀(如图 8(a)所示)，处于受压状态。最大应力发生在安全壳底板与柱体斜坡交界处，这是由于此处的突然转折引起的。在环梁与筒体交界处，由于穹顶预应力筋和竖向预应力筋的作用，使得此处所受压应力较大。；当内压达到0.8MPa时，孔口部位的最大拉应力达到120MPa (图 8(b))；而当内压达到1.0MPa时，孔口部位的最大拉应力已达到294MPa (图 8(c))，但还低于钢衬里材料的屈服极限，还能够保持它的设计功能。钢衬里应力的快速增长，是与混凝土的损伤开裂密切相关的。由于在0.8MPa下，孔口附近的混凝土出现了加大范围的损伤开裂 (图 6(b))，因此在1.0MPa下，钢衬里的应力较0.8MPa时的应力增加了近1.5倍。

图 8 钢衬里应力云图 Fig. 8 Stress contours of steel lining

图选项

本文给出了能较好模拟预应力混凝土安全壳的有限元分析方法，针对实际的CPR1000安全壳进行了详细的内压分析，从计算结果分析可得出以下结论：

1) 本文所采取的安全壳建模方法能很好地模拟预应力筋与混凝土之间的相互作用，能够较好的反映安全壳在内压作用下的应力分布及开裂位置。

2) 在内压作用下应力较大的区域集中在扶壁柱、设备孔周围、环梁附近和筒体底部的位置，这些部位应是核电厂安全壳的薄弱部位，尤其是设备孔周围和基底部位应予以重视。

3) 在0.4MPa设计事故内压作用下，安全壳没有开裂，满足核电厂安全壳的设计使用要求。随着内压的增大，安全壳首先在设备孔上下侧开裂，并逐步扩展，然后设备孔左右两侧也出现了裂缝，在0.8MPa时，孔口附近的混凝土开裂严重，但钢衬里的最大应力仍然小于材料的屈服极限，安全壳还可以起到它防辐射的功能。