乏燃料容器坠落事故工况下核燃料厂房的安全性分析

2019-04-03郭全全蔡利建

振动与冲击 2019年6期

程刚，郭全全，周耀，马英，蔡利建

(1.北京航空航天大学交通科学与工程学院，北京 100191； 2.中国核电工程有限公司，北京 100840)

由于核电工程的特殊性，核电厂在遭受外部撞击[1-2]和内部事故工况下的安全性越来越受到重视，其中乏燃料容器可能发生的坠落事故已成为一个基准工况，在设计时需要采取专门的工程措施并进行安全性分析。

本文应用ANSYS/LS-DYNA对乏燃料容器的坠落进行了数值模拟。ANSYS/LS-DYNA强大的非线性分析功能已在爆炸、冲击等问题的分析中的得到广泛的应用和工程验证。陆新征等[3]对世界贸易中心受撞后的倒塌进行了力学和仿真分析，并根据结构进行了参数讨论。张涛等通过建立精细化的空客A320模型，模拟了100 m/s冲击带内衬钢板的核电站安全壳。朱秀云等基于荷载时程分析法，进行了波音707-320型号商用飞机撞击某钢板混凝土结构安全壳的响应分析。宋晓滨等[4]对重物高空坠落高速撞击钢筋混凝土楼盖时的结构损伤进行了分析。张斌等[5]模拟了爆炸冲击对带泡沫混凝土回填层隧道衬砌结构的动力响应。何思明等[6]以实际滚石防护结构棚洞为原型，对滚石冲击过程进行了模拟，研究了不同冲击角度下滚石对棚洞结构冲击的动力响应。

根据核电站的工艺布置，满载总重量约为113 t的乏燃料容器在运输过程中可能在装卸口处坠落，经过约17 m的自由落体运动后冲击-0.05 m标高处0.8 m厚的钢筋混凝土楼板。在此冲击过程中，乏燃料容器可能击穿该处楼板而继续以一定初速度作自由落体运动，最终冲击-12.5 m标高处2.2 m厚的核燃料厂房筏板。

为保证厂房整体结构的安全、减小乏燃料水池的变形和振动输入，设计时在装卸口筏板上铺设5 m厚的蒸压加气混凝土(Autoclaved Aerated Concrete,AAC)砌块形成减振层，通过AAC砌块的压溃和开裂耗散冲击能量，减小燃料厂房的冲击响应，如图1所示。

图1 乏燃料容器坠落事故示意图Fig.1 Schematic view of spent fuel container dropping accident

本文应用ANSYS/LS-DYNA对AAC砌块的落锤冲击试验进行了数值模拟，标定了AAC的本构模型参数，进而建立厂房的实体-分层壳有限元模型，模拟了

乏燃料容器坠落对厂房整体结构的冲击动力响应，并从混凝土主应变、水池跨中挠度以及楼层反应谱三个方面进行了安全性分析。

1 AAC材料参数的有限元标定

AAC是一种多孔的脆性材料，其力学特性和混凝土相似但又存在差别，在ANSYS/LS-DYNA众多材料模型中，尚无特定适用于AAC的模型。

本文首先参考美国混凝土学会ACI544.2R标准中关于纤维混凝土(FRC)的落锤试验方法，采用XJL-300B型双管式导向落锤冲击试验系统，测得不同落距下经过减振层后的冲击力时程曲线。试件采用B04级AAC块体，试件为直径d=0.15 m、厚度h=0.1 m的圆柱体。试验时，试件置于厚度为0.02 m板上，在钢板下方布置力传感器。质量5.31 kg的落锤分别从落距H=1.0 m，H=1.5 m，H=2.0 m处自由下落冲击AAC试块，采用量程20 kN，采样频率为100 kHz的传感器采集冲击力时程曲线，具体安装形式如图2(a)所示。

采用与试验试件完全一致的几何尺寸，建立有限元分析模型(见图2(b))，其中，落锤、钢板、传感器均采用*MAT_ELASTIC弹性材料模型，AAC参考Li等对爆炸冲击作用下AAC砌体结构的研究中所采用的MAT96材料模型。各接触面之间均采用自动面面接触，约束下钢板底面节点的所有自由度作为边界条件。通过逐步调整MAT96的参数，最终取模拟冲击力时程曲线与试验曲线最接近的参数作为厂房整体响应分析中AAC的材料参数，见表1。

图2 AAC材料模型参数的有限元标定Fig.2 FE calibrating of constitutive model parameters for AAC

模型名称密度/(kg·m-3)弹性模量/Pa泊松比抗剪强度/Pa抗拉强度/Pa断裂能/(N·m-1)剪力保留系数黏度/(Pa·s-1)抗压强度/Pa失效应变Li et al[7]6255.3×1080.21.0×1060.7×106800.037.17×1053.07×1060.01本文4001.5×1090.153.5×1051.5×105900.035×1052.44×1060.25

图3对比了落锤冲击试验和有限元计算所得到的冲击力时程曲线。选取冲击力时程曲线中的峰值力、冲击持续时间、冲量作为关键参数，其中试件编号的含义为：以B04-0.1-5.31-1.0为例，三组数据分别代表质量为5.31 kg的落锤从1.0 m高处冲击直径0.15 m、厚度0.1 m的B04等级AAC试块。对其试验值和模拟值作进一步对比，详见表2。

图3 试验-模拟冲击力时程曲线对比Fig.3 Time-history curve of impact forcebetween test and FE analysis

试件编号冲击力峰值/kN持续时间/ms冲量/(N·s)B04-0.1-5.31-1.0试验值5.8411.3226.72模拟值5.6510.1226.15相对误差/%3.2010.602.12B04-0.1-5.31-1.5试验值6.8412.0033.05模拟值6.8813.0030.47相对误差/%-0.58-8.307.81B04-1.0-5.31-2.0试验值8.5612.2140.19模拟值8.0011.7032.98相对误差/%6.464.1817.94注:相对误差=(试验值-模拟值)/试验值×100%

有限元分析结果与试验结果对比可知，0.1 m厚的AAC试件在5.31 kg落锤1.0 m，1.5 m，2.0 m落距冲击下，冲击力峰值、冲击持续时间、冲量均表现出良好的一致性，验证了参数取值的有效性，能够较好地模拟落锤冲击作用下AAC材料的动力响应。

2 厂房有限元模型的建立

由于燃料厂房结构庞大、复杂，本项目研究的是冲击作用下的厂房总体响应，因此，在不影响分析精度的基础上对模型作适当简化。

2.1 单元及材料选择

乏燃料容器选用Solid164单元，将其简化成均质弹性材料，按照质量和体积不变的原则作等效折算，密度为2 234 kg/m3，弹模和泊松比按钢材设置为2.0×1011Pa和0.3。

厂房主体选用Shell163单元。根据分层壳原理(见图4)，将壳单元划分为5层：保护层-钢筋层-混凝土层-钢筋层-保护层，保护层取0.05 m，以实际配筋率折算其中钢筋层的厚度，混凝土层为总厚度减去其他层的厚度。C40混凝土选用*MAT_PLASTIC_KINEMATIC(MAT3)材料模型，密度2 450 kg/m3，弹性模量3.25×1010Pa，屈服强度1.91×107Pa，断裂失效塑性变形设定为混凝土最大受压变形,取0.003。HRB500钢筋也选用MAT3材料模型，密度7 800 kg/m3，弹性模量2.0×1011Pa，屈服强度4.35×108Pa，断裂失效塑性变形设定为0.075。

图4 分层壳原理示意图Fig.4 Schematic view of the layered shell

冲击处楼板选用Solid164单元，材料选用*MAT WINFRITH_CONCRETE(MAT84)模型[8]，其断裂能取值见文献[9]。通过关键字MAT_ADD_EROSION控制混凝土材料失效。以主压变作为混凝土材料失效准则，阈值设为0.144。按照混凝土结构设计规范GB 50010—2010选择C40混凝土的对应参数选取，具体取值见表3。

表3 混凝土底板参考和使用的MAT84材料参数Tab.3 Parameters in model MAT84 for concrete plate

Solid164楼板单元和Shell163墙体单元通过局部细化墙体单元网格，与楼板单元网格划分完全重合，通过共节点实现楼板与墙体的协调变形(见图5)。

地基基础选用Solid164单元，在厂房结构的长、宽、埋深方向分别扩展三倍，选用*MAT_ELASTIC弹性材料模型。密度2 150 kg/m3，弹性模量为2.15×109Pa，泊松比为0.39。

AAC减振层选用Solid164单元，选用冲击试验标定的MAT96材料模型，具体参数见表1。

水池水体选用Solid164单元建模，通过更改K文件关键字*SECTION_SOLID将水体设定为流体，*MAT_NULL和*EOS_GRUNEISEN定义水体材料和状态方程。

图5 Solid164单元与Shell163单元共用节点示意图Fig.5 Co-node of Solid164 and Shell163 elements

2.2 约束条件和荷载施加

接触：Shell单元与Solid单元、Solid单元之间采用自动面面接触。混凝土底板、减振材料、混凝土墙体接触面静摩擦因素为0.5，动摩擦因素0.3。

约束：在地基底面施加竖直方向的位移约束，在地基侧面和底面节点设置无反射边界，自动吸收膨胀波和剪切波，模拟半无限大地基边界[11]。

流固耦合：添加算法控制关键字*CONTROL_ALE和流固耦合关键字*CONSTRAINED_LAGRANGE_IN_SOLID实现流固耦合[12-13]。

荷载施加：向乏燃料容器施加竖直方向-18.25 m/s的初速度来模拟乏燃料容器从17 m高度坠落的冲击。

2.3 厂房结构阻尼参数确定

采用Rayleigh阻尼，可表达为质量矩阵M和刚度矩阵K的线性组合，即：C=αM+βK，α和β分别为质量比例系数和刚度比例系数。

传统方法是通过指定两阶“参考”频率(ωi和ωj)及其阻尼比(ξi和ξj)进行计算。此法缺点在于两阶参考频率之间的模态阻尼比偏小而在参考频率以外的模态阻尼比偏大。为了得到更加贴合结构实际的阻尼系数α和β，刘红石[14]通过考虑结构振型与阻尼矩阵的正交性，利用最小二乘法原理，提出了一种综合考虑前n阶频率的Rayleigh阻尼参数确定方法，即

采用次空间法求解厂房的50阶模态见表4。

表4 核燃料厂房前50阶模态表Tab.4 The first 50 order model of nuclear fuel plant

代入前50阶频率，求得

α=2.725β=0.000 75

2.4 厂房结构整体模型

厂房结构整体模型见图6，整个模型划分了22 832个实体单元， 80 382个壳单元。

图6 核燃料厂房整体模型三视图Fig.6 Three view drawing of nuclear fuel plant model

2.5 减振层厚度的选择

图7为AAC立方体(边长100 mm)的静力压缩应力应变关系曲线，可以近似认为，在峰值应变前为弹性阶段，峰值应变后迅速进入压溃平台阶段，即有效的减振吸能阶段。本项目参照秦山核电站二期工程将减振层厚度设为5 m，考虑到工程的经济性，分析了厚度为3 m和5 m时减振层的主压应变云图(见图8)，可以看出减振层厚度为3 m时，整个厚度均已进入塑性吸能阶段(主压应变超过2×10-2)，而5 m时，底部还有一定厚度减振层未进入塑性吸能阶段，因此偏于安全地仍选减振层厚度为5 m。

图7 AAC立方体压缩应力应变关系曲线Fig.7 The static compression stress-strain relationship curve of AAC cubic

图8 乏燃料容器冲击下减振层压应变云图Fig.8 Compression strain of the AAC layer under the impact of the spent fuel container

3 厂房结构冲击响应及分析

3.1 混凝土应变

乏燃料容器坠落过程(见图9)：首先冲击-0.05 m处0.8 m厚的钢筋混凝土楼板，此时混凝土楼板达到设定的失效应变(0.144)而删除(见图10)；周边侧墙混凝土最大拉应变、最大压应变分别为1.45×10-3和1.65×10-3，墙体局部拉裂(见图11)。乏燃料容器在穿透楼板后继续坠落冲击-7.5 m处5 m厚减振层，此时筏板混凝土最大拉应变和压应变分别为1.37×10-4和1.03×10-4，筏板轻微拉裂(见图12)；减振层周边侧墙最大拉应变和压应变分别为6.04×10-4和1.02×10-3，出现局部拉裂(见图13)。