基于土-结构相互作用的核电厂房抗飞机撞击耦合分析与模型验证

2022-03-27王友刚孙运轮李建波李晓红梅润雨牛燕如

振动与冲击 2022年6期

王友刚，孙运轮，李建波，李晓红，陈岩，梅润雨，牛燕如，林皋

(1. 中核能源科技有限公司，北京 100193；2. 大连理工大学海岸与近海工程国家重点实验室，辽宁大连 116024；3. 大连理工大学工程抗震研究所，辽宁大连 116024)

发展核电是国家重要战略需求，确保核电安全是未来内陆核电发展的先决性条件[1]。受“911”事件的重大启示，大型飞机撞击安全壳及其可能的严重后果日益受到人们的关注，尤其近年来各类大型飞机在航空业应用广泛，翼展达数十米，与安全壳直径几乎相当[2]，在以100～200 m/s的速度产生的巨大冲量作用下，可能对安全壳产生极大的冲击荷载。国内外相继出台各项法规[3-4]，正式将抵御大型商用飞机恶意撞击列入核电站建造的考虑范围。

飞机撞击安全壳的研究最早开始于1968年Riera公式[5]，主要是聚焦在不同因素对飞机撞击荷载曲线的影响。随着数值分析技术的快速发展，基于飞射物-靶体相互作用的精细化分析方法，日益在核电结构抗飞机撞击中得到广泛应用[6]，这一领域研究的广度和深度也在逐渐扩展和完善[7]，如撞击位置，安全壳墙体曲率，撞击过程飞机速度衰减特征等。

Sadique等[8]使用4种不同型号飞机撞击安全壳，并对其最大变形和整体性能进行了对比；Jeon等[9]对采用纤维增强混凝土的核安全壳结构抗冲击性能进行了分析；黄涛等[10]建立了精细化的大型商用飞机与核电站安全壳(厂房)有限元模型，并对飞机撞击规律与特点进行了总结；Zou等[11]提出一种将比例边界有限元和八叉树技术相结合的跨尺度分析方法，对飞机撞击核安全壳问题进行了精确且高效的模拟。

安全壳等核电厂房结构在飞机撞击作用下的动力响应，是一种短时、多因素影响的强耦合过程。与采用简单的撞击荷载曲线方法相比，精细化的相互作用分析方法具有显著的优势[12]，尤其是便于从飞机-厂房-地基耦合体系的角度开展机理模拟。内陆非岩基厂址条件下，地基对上部结构的刚度约束作用大幅减弱，大飞机撞击这一冲击载荷的耗散过程及其对安全壳-地基相互作用体系的整体动力影响是值得关注的重要课题，对这一问题的系统性分析，将对核电结构抗飞机撞击领域的发展起到积极的促进作用。

土-结构相互作用框架下，飞机撞击和地震载荷所采用的地基动力模型，具有明显的异同，飞机撞击对于地基外边界来说是典型的内源动力问题，而地震则是外源动力问题[13-15]，如图1、图2所示。从地基中散射波动的单向传播角度看，两者地基外边界的透射条件要求本质上是一致的。

图1 地震荷载

图2 飞机撞击荷载

相互作用体系抗大飞机撞击的数值模拟主要涉及以下几方面问题：上部精细化飞机模型和安全壳模型的建立；相互作用模型的建立；非线性材料本构的选择与参数调整；冲击破坏全过程的数值模拟。上述数值模型较为复杂，片面追求精细化的有限元模型在提高分析精度的同时，也会带来工作量的巨大压力，在保证精度和效率的双重要求下，如何有效平衡两者，是本文的关注点之一。

目前从相互作用体系角度探讨核电厂房结构的抗飞机撞击过程仍缺乏标准化的分析方法。本文提供了一种考虑相互作用条件下飞机撞击分析的方法，主要聚焦在精细模型建立及参数的有效性验证、可模拟无限地基辐射阻尼效应的人工边界设置、以及不同软硬程度土体对上部结构振动响应的影响，并基于此建立了飞机-厂房-地基全耦合分析系统。

1 有限元模型

1.1 飞机模型

1.1.1 飞机模型参数

本文飞机模型使用的是某四引擎大型客机，如图3所示。机身长63.6 m，翼展60.3 m，总质量为234 t，由190 290个单元组成。在建模时使用MASS166单元为其增加配重，使之满足真实的质量分布。

图3 飞机有限元模型

飞机蒙皮(机身、机翼部位)为航空铝材料，引擎的材料是刚，均使用Jackson-Cook本构模型，飞机内部的钢筋仍考虑采用与核岛厂房钢筋相同的本构，材料参数参考文献[16]。

1.1.2 飞机模型验证

为确保用于计算的飞机模型真实有效，通过比较采用数值方法计算的撞击刚性靶板的撞击力和采用Riera方法计算的撞击力对飞机模型的合理性进行比较。分别建立网格尺寸差异明显的两种网格模型，对其网格敏感性进行分析。

本文分别使用不同大小的网格对模型进行剖分，模型1主要采用的壳单元网格尺寸约为(280×155)mm2，梁单元网格尺寸约为155 mm，该模型的节点数为110 601，单元数为190 290，如图4所示。模型2主要采用的壳单元网格尺寸约为(140×75)mm2，梁单元网格尺寸约为75 mm，该模型的节点数为310 663，单元数为681 586，如图5所示。

图4 模型1网格

图5 模型2网格

采用数值方法进行飞机撞击刚性靶板的撞击力计算，是通过对飞机进行精细化的有限元建模，以获得真实的飞机撞击刚性靶板过程。在进行分析时，假定飞机以200 m/s的速度垂直撞击刚性靶板，如图6所示。

图6 飞机垂直撞击刚性靶板

通过以上两种方法，分别得到飞机撞击刚性靶板的撞击力时程曲线和冲量时程曲线对比图，如图7、图8所示。

图7 撞击力时程对比

图8 撞击力冲量时程对比

通过以上结果可以看出，模型1和模型2与Riera方法的结果吻合得很好，以此验证了两种有限元模型的合理性。同时，通过结果对比也可以看出，两种模型的计算结果相差不大，然而模型1的计算时间约为4 h，模型2的计算时间约为20 h。因此，采用模型1的飞机模型进行后续分析。

1.2 核岛厂房模型

以某核电站示范项目为背景建立精细化的厂房模型。主要由反应堆厂房、乏燃料厂房、核辅助厂房和电气厂房构成，如图9(a)所示。反应堆厂房细致模拟了外部安全壳结构与内部舱室结构，如图9(b)、图9(c)所示。

核岛厂房整体由混凝土浇筑，在安全壳穹顶内侧与舱室结构外侧有钢板，安全壳壳体内部布设有环向钢筋与纵向拉筋，如图9(d)所示。利用ANSYS/LS-DYNA软件建立了如图9(a)～图9(d)所示有限元模型，其中混凝土使用SOLID164单元，钢板使用SHELL163单元，钢筋使用LINK160单元，钢筋与混凝土采用分离式建模方法，并通过*CONSTRAINED_LAGRANGE_IN_SOLID关键字耦合在一起，单元尺寸为0.3 m，单元总数为4 232 520个。

图9 有限元模型

1.2.1 核岛厂房本构模型参数

不同载荷作用下混凝土材料的率敏感性差异较大，如图10所示。冲击荷载下应变率约100～101s-1。在动力学问题分析中，混凝土结构的损伤及演化是十分复杂的过程，若还使用混凝土的静态本构则计算结果的可信度无法保证，可有效模拟材料的率相关动力特性是DYNA软件的突出特征。基于LS-DYNA中丰富的材料库，考虑使用混凝土连续面盖帽模型进行飞机撞击问题的分析。该模型考虑了材料的硬化、损伤以及率相关性，目前在钢筋混凝土结构低速冲击领域应用广泛。

图10 不同荷载下混凝土应变率

核岛厂房使用C40混凝土材料，根据GB 50010—2010《混凝土结构设计规范》，抗压强度为27 MPa，密度为2 400 kg/m3，这两个参数为混凝土本构所需输入参数。核岛厂房中钢板厚度为20 mm，环向钢筋直径为40 mm，纵向拉筋直径为20 mm，密度均为7 800 kg/m3,泊松比为0.3，使用可以考虑应变率与失效的Plastic-Kinematic(塑性随动硬化)本构模型，材料参数如表1所示。

表1 核岛厂房钢筋、钢板材料参数

1.2.2 核岛厂房网格敏感性分析

在数值模拟之前，需要对有限元模型进行网格敏感性分析，以保证计算的效率和准确性。为了选出最合适的网格尺寸，以网格尺寸和反应堆厂房沿径向划分网格个数为标准，构建5种不同类型的模型。各个模型的节点、单元个数等信息，如表2所示。

表2 不同类型厂房模型网格信息

在进行分析时，使用1.1.2节分析中选择的飞机模型1，飞机以100 m/s速度垂直撞击反应堆厂房墙体，撞击点坐标为：(52 015.0,61 116.9,16 434.0)，如图11所示。通过对比这5种模型的最大变形位移，进而选择合适的网格尺寸。

图11 核岛厂房网格敏感性分析飞机撞击位置图

将计算得到的不同模型反应堆厂房墙体在飞机撞击过程中，最大变形位移结果如表3所示。从表3数据可知，300_R006模型和400_R006模型最大变形值相差不大，而500_R006比这两个模型增加了10 mm多。反应堆厂房模型最大变形值沿着径向6层单元和5层差不多，4层则比这两个模型增加了50 mm多。从计算效率来看，300_R006的网格数量超过400_R006约20万，如果采用300 mm网格尺寸会大大增加计算耗时，而6层比5层则超过5万多，对计算时间影响不是很大，最终选择了400_R006模型，整体网格尺寸400 mm，反应堆厂房外墙沿着径向分布有6层单元。

表3 不同模型最大变形位移

2 考虑土-结构相互作用效应的飞机撞击分析

在工程应用中，通常将地基假定为刚性，在核岛厂房底板底部施加固端约束。这种假定在柔性地基场地条件下，与实际情况存在较大差异。因此，本文考虑使用直接法，对完整的地基场地进行建模，并在地基外边界设置虚拟人工边界，以模拟无限地基辐射阻尼效应。

2.1 结构地基相互作用模型

该相互作用模型是在核岛厂房模型的基础上，增加了地基模型而构建的。地基模型长300 m，宽250 m，深45 m，如图12所示。

图12 地基模型(m)

由于核岛厂房外轮廓并非规则几何图形，因此在网格剖分时，靠近核岛厂房局部范围内，根据核岛厂房外轮廓进行合理剖分，其余整体网格尺寸为2 m。整个地基模型共包含约433万个节点和493万个单元。在地基外侧5个立面上，设置虚拟人工边界——完美匹配层(perfect matching layer，PML)，以考虑地基无限半空间的影响。

考虑到场地类型的多样性，选取了AP1000 3种标准设计的场地类型[17]以及一种满足剪切波速包络性的软土场地，进行算例对比分析。地基土采用理想弹塑性Drucker-Prager本构模型，该模型考虑了土的非线性。不同场地土体参数参考GB 50218—1994《工程岩体分级标准》进行设置，如表4所示。

表4 不同场地土体参数

同时，还考虑了3种重要撞击位置对结构响应的影响，分别是飞机水平撞击核岛厂房外壳(撞击位置1,)，飞机水平撞击空冷塔(撞击位置2)和飞机斜向垂直撞击穹顶(撞击位置3)，如图13、图14、图15所示。

图13 水平撞击核岛厂房

图14 水平撞击空冷塔

图15 斜向垂直撞击穹顶

2.2 虚拟人工边界

PML是在计算区域外人为添加的一具有一定厚度、含任意吸收介质的吸收层。为验证PML的精度和可靠性，对图16所示具有解析解的半无限介质波动问题，分别通过解析和数值方法进行计算。

图16 内源问题

地基模型在构建时，由于基础底板形状不规则，在网格剖分时，难免会出现三棱柱网格。为分析这样的网格是否对计算结果有影响，构建了如图17所示两种模型。

图17 不同网格PML模型

模型1包含三棱柱单元，模型2全是六面体单元。该内源问题的对比结果，如图18所示。从图18可知，模型1和模型2结果非常吻合，可见三棱柱单元的影响可忽略。虽然PML和黏弹边界的结果同解析解有些差别，但这两种方法的计算结果相差不大，PML的计算结果精度同黏弹边界的计算精度相当。

图18 结果对比

2.3 结果分析

2.3.1 不同场地条件的影响

模型中设置了两个监测点对计算结果进行分析，点A设置在反应堆厂房底板中心处，点B设置在核岛厂房某个舱室的内壁上，如图19所示。

图19 监测点

在分析不同场地对上部结构振动响应的影响时，飞机以100 m/s速度撞击。飞机撞击位置1时，点A的三向加速度反应谱结果如图20所示，点B的三向加速度反应谱结果如图21所示。

图20 撞击位置1点A的加速度反应谱

图21 撞击位置1点B的加速度反应谱

从不同场地的对比结果可以看出，在约4 Hz频率处，加速度反应谱存在第一个峰值，该频率与结构基频相当，在此频率下结构共振响应较大。

加速度峰值作为衡量振动响应程度的重要指标，其值可通过加速度反应谱在极高频处的谱值进行估计。由于Y向是飞机撞击方向，点A在不同工况下的Y向加速度反应谱峰值及1 000 Hz处的加速度峰值，如表5所示。

表5 撞击位置1点A处撞击正向Y动力响应

从监测点的对比结果可以看出，场地对反应堆厂房振动响应有明显影响。工况4场地的加速度反应谱普遍高于其他，此时场地条件最为柔软，因此我们可以判断场地越软，加速度反应谱将会越大，此时，土结构相互作用效应较为显著。同时，加速度反应谱在极高频处的谱值等于加速度响应的峰值，从工况1～工况4逐渐增大反映出，场地条件越柔，结构的振动响应越剧烈。因此，在核电厂房抗飞机撞击引起的核电结构整体动力响应中，需要考虑土-结构相互作用。

2.3.2 撞击位置影响

飞机撞击不同位置时，点A的三向加速度反应谱结果对比如图22所示，点B的三向加速度反应谱结果对比如图23所示。