基于参数化-超单元法的HFETR乏燃料贮存格架抗震分析

刘金龙，章航洲，王 帅，孙志军，隆 涛

(中国核动力研究设计院，四川 成都 610213)

乏燃料贮存格架(简称格架)是用于贮存乏燃料组件的关键设备，在地震中可能出现变形、碰撞、倾倒等现象，影响乏燃料组件的贮存安全。目前，格架抗震分析一般采用数值模拟，常用的分析模型有梁模型和壳模型。梁模型计算速度快，但精度较低[1-3]；壳模型精度较高，但计算时间长，占用资源多，很难实现多格架同时分析[4-6]。另外，格架结构在满足抗震安全的同时还须符合临界、热工等要求，设计中可能需要多次调整结构尺寸参数、重复建模和抗震计算，而传统的CAD设计、CAE分析建模效率较低，导致格架迭代设计周期大幅增加。

本文在对格架抗震分析模型及建模方式调研和分析的基础上，提出基于超单元法建立高通量工程试验堆(HFETR)格架分析模型，解决模型精度和计算效率不能兼顾的问题[7]；然后基于ANSYS APDL程序实现格架分析模型的参数化，通过VB平台对格架程序进行封装和可视化开发，建立格架分析界面，解决格架建模效率低的问题；最后，利用该方法对整池HFETR格架开展抗震分析。

1 超单元分析模型的建立

HFETR乏燃料组件贮存格架主要部件有贮存单元、栅板、围板、加固梁、底板、底座等，最大外形尺寸为2.4 m×1.8 m×5.2 m(长×宽×高)，格架结构如图1a所示。贮存单元为3×5排列，内部存放乏燃料组件，底部与底板焊接，中部和上部与栅板孔间隙配合。围板与加固梁、底板、栅板之间通过焊接连接。共设计6台乏燃料贮存格架，每台格架单独放置在水池中，相互间不连接，格架与池底间不固接。

1.1 单元选择及网格划分

格架贮存单元、围板、栅板、底板以及底座底板采用壳单元SHELL181模拟；格架外围加固梁、底座与凸台之间的支撑柱以及组件采用梁单元BEAM188模拟；格架与底板接触部位使用实体单元SOLID185模拟；格架围板与底板、栅板与围板、贮存单元与底板、凸台与底板设为绑定接触，接触单元使用CONTA175、TARGE170，为避免接触单元影响超单元生成，设KEYOPT(2)=1(罚函数法)[8-9]。网格划分后的有限元模型(称普通有限元模型)如图1b所示。

1.2 超单元模型建立

根据格架结构和接触边界情况对普通有限元模型进行子结构划分和主节点选取，将底座上半部分及以上的底板、围板、栅板、贮存单元作为超单元，底座部分支撑柱和底垫板部分保留为非超单元，格架外部轮廓节点作为主节点，生成的超单元模型如图1c所示。

a——格架几何模型；b——普通有限元模型；c——超单元模型

1.3 模型验证

分别对超单元模型和普通有限元模型底座底面施加固定约束并进行模态分析，两种模型的模态频率列于表1，振型示于图2、3。由表1可见，超单元模型与普通有限元模型的模态频率相差较小，最大相对偏差仅1.122%。由图2、3可见，两种模型的振型基本一致。表明本文所建超单元模型精度满足要求，可用于格架抗震分析。

图2 普通有限元模型前3阶振型Fig.2 First three mode shapes of ordinary finite element model

表1 两种模型的模态频率比较Table 1 Comparison of modal frequencies between two models

2 基于APDL和VB的格架模型参数化

格架模型参数化主要包括几何模型参数化、边界条件参数化以及可视化分析界面设计。

2.1 几何模型参数化

本文基于APDL程序建立格架几何尺寸变量和各部件间尺寸关系，实现格架从内到外、从上到下的尺寸驱动建模，如图4所示[10-11]。为使各部件之间保持关联，利用APDL关键点命令建立关键点坐标，通过点驱动方式实现格架部件之间尺寸的驱动，部件连接部位共用关键点，在模型修改时可避免部件之间出现交叉或分离等情况，保证了格架整体结构的连续性。组件的阵列主要借助APDL的AGEN、VGEN、LGEN等命令实现，通过输入组件排列方式，自动阵列贮存单元；多格架的排列基于超单元模型，通过SETRN、SE、*DO等命令可实现超单元模型提取、摆放距离调节和摆放组合变化等。最终通过程序设计实现格架模型的建立、尺寸修改和驱动、格架内组件装载和多格架摆放方式的变化等，为格架结构及布置迭代设计提供有力支持。

图3 超单元模型前3阶振型Fig.3 First three mode shapes of superelement model

图4 格架各部件驱动示意图Fig.4 Diagram of rack components drive

2.2 边界条件参数化

格架抗震分析时所需考虑的边界条件主要包括水环境、地震、阻尼、摩擦等。其中水环境的仿真最为复杂和困难。在地震中，格架周围的水对格架运动的影响十分显著，水对格架的影响不能忽略。目前，模拟格架与水、组件与水之间流固耦合作用常用的方法是附加质量法，该方法适用于小位移附加质量的计算，水动力质量矩阵形式[12-13]如下：

(1)

式中：MH为流固耦合水动力质量，本格架x、y方向为水平方向，分别为MHx和MHy；M1为格架所处位置置换水的质量；M2为边界内无格架情况下水的质量。

格架水间隙示意图如图5所示。对于图5所示的两个水下矩形柱体，中间矩形柱体代表格架，外围矩形柱体代表周围池壁或格架，中间间隙充满水。两柱体在水平方向上相对运动而产生的水动力质量在数值上表示为：

图5 格架水间隙示意图Fig.5 Schematic diagram of water gap of rack

MHx=2ρhC2(C/3g1+C/3g3+B/g2+B/g4)

(2)

MHy=2ρhB2(B/3g2+B/3g4+C/g1+C/g3)

(3)

(4)

(5)

式中：ρ为水质量密度；h为格架高度；C、B为考虑水间隙的标称尺寸；g1、g2、g3、g4为格架周围间隙尺寸。各格架外部的附加质量通过APDL编程后实现自动计算，并自动赋值给流固耦合单元FLUID38实现附加质量的模拟。乏燃料组件内含水以及与贮存单元之间的水间隙附加质量，使用MASS21质量点分别施加在组件和贮存单元的水平方向上[14]。

地震作用采用加速度时程模拟[15-16]，考虑运行基准地震(OBE)和安全停堆地震(SSE)两种地震工况，地震加速度时程曲线如图6所示，分别在x、y、z3个方向上同时施加地震载荷，持续时间20 s。地震数据通过*CREATE、*VREAD、ACEL等命令实现自动载入和施加。

图6 地震加速度时程曲线Fig.6 Seismic acceleration time history curve

格架自由放置在乏燃料水池中，为重点关注格架滑移情况，格架底座与水池底面的摩擦系数保守设为0.2[17]。结构阻尼采用瑞利阻尼，OBE时阻尼比取2%，SSE时阻尼比取4%，根据模态分析结果选择感兴趣的2.2～30.1 Hz进行计算[5]。OBE时，α=0.515，β=1.97×10-4；SSE时，α=1.03，β=3.94×10-4。单格架模型及边界条件加载情况如图7所示，格架底座与水池底面采用面面接触(CONTA175、TARGE170)，抗震计算采用完全法瞬态动力分析(ANTYPE、TRANS)。

图7 单格架模型环境设置示意图Fig.7 Schematic diagram of single frame model environment setting

2.3 可视化分析界面设计

格架模型参数化程序编写完成后，利用VB对APDL程序封装并建立参数关联，通过开发可视化界面，使格架参数修改更加便捷，建模和分析时通过调用ANSYS软件(Batch功能)在后台运行，将命令流提交给ANSYS进行分析，结束后自动关闭ANSYS，最后使用VB命令提取计算数据，实现格架分析结果的查看，其流程及原理如图8所示。

图8 VB调用ANSYS原理图Fig.8 VB calling ANSYS schematic diagram

基于VB可视化平台开发的界面可实现格架模型建立和修改、边界条件加载、分析计算以及结果查看等整个流程的操作，主界面如图9所示。

图9 格架分析主界面Fig.9 Main interface of rack analysis

3 乏燃料贮存格架抗震分析

本文分别选取满载、半载、空载3种典型装载形式对整池6台格架开展抗震分析，格架装载情况如图10所示。格架之间按设计间距摆放，6台格架超单元模型如图11所示。根据计算结果对格架应力、碰撞、倾倒情况进行分析。

图10 单格架装载示意图Fig.10 Schematic diagram of single rack loading

图11 6台格架超单元模型示意图Fig.11 Schematic diagram of superelement model of six racks

3.1 格架应力分析

6台格架满载、半载、空载时在OBE、SSE中底座对水池底的最大应力如图12所示，由于格架重心向-x、-y方向倾斜，格架左上角的底座对水池的应力强度较大，格架1在半载、SSE情况下与水池发生最大相互作用，发生时刻为9.79 s。对格架1超单元模型进行扩展，得到格架1最大应力发生在底座。底座结构如图13所示，格架小幅度倾斜时，底垫板不会传递竖直方向的扭矩，螺柱底部受到的作用力由螺柱传递给螺套，竖直方向的力由螺套肩部传递给底板，水平方向的力由螺套颈部传递给底板。格架底座最大压应力发生在螺柱与螺套连接处，如图14所示，最大值为230 MPa，满足应力要求。

图12 水池底应力云图Fig.12 Stress intensity cloud image at bottom of pool

图13 底座示意图Fig.13 Schematic diagram of base

图14 底座应力强度最大部位云图Fig.14 Maximum stress intensity cloud image of support feet

3.2 格架滑移分析

格架在OBE、SSE工况下x、y方向的最大滑移位移曲线如图15所示。在x方向，格架空载情况下在SSE工况时产生的滑移距离最大，为-15.4 mm；同理，y方向最大滑移距离发生在格架空载且处于SSE工况时，为-27.4 mm，均小于格架之间、格架与水池之间的间隙，不会发生碰撞。在OBE、SSE两种地震工况下均显示格架空载时对地震的响应最剧烈，在同一时刻等加速度条件下，滑移幅度从小到大依次为满载、半载、空载，即随着组件装载量的增加，格架整体质量增加，格架滑移幅度有减小的趋势，因此在格架滑移分析时空载情况应作为重点工况进行验证。

图15 地震工况下格架在x、y方向的最大滑移位移曲线Fig.15 Maximum slip displacement curves of rack in x and y directions under seismic condition

另外，格架在半载情况下累积的滑移位移与满载、半载情况下有方向相反的趋势，这是由于受格架结构及装载情况影响，质量大小及分布不同，格架受地震影响产生的倾斜姿态、底座摩擦力分布等情况各异，使得格架滑移位移累积效果各不相同。

3.3 格架倾倒分析

为保守估计格架倾倒情况，将抗震计算得到的格架底部最大起跳高度与格架倾倒限值比较。格架倾斜示意图如图16所示。由图16可知，格架倾斜时重心抬高，外侧底座为主要支撑点，当重心位于支撑点正上方时，重心抬升至最高值，超过该点即会发生倾倒。格架倾斜时重心抬高距离Δh按下式计算：

图16 格架倾斜示意图Fig.16 Diagram of rack tilt

(6)

式中：h为格架重心高度；l为支撑点与重心的垂直距离。格架重心抬高Δh时，支撑点另一侧底座将抬高Δz(该高度视为格架倾倒限值)，l越小倾倒限值越小。根据以上分析，分别计算得到格架满载、半载和空载时最小的倾倒限值分别为531、472、557 mm。

地震工况下格架在z方向的最大起跳高度曲线示于图17。可见，在OBE、SSE工况中，格架均在半载情况下获得最大起跳高度，格架半载时OBE工况下最大起跳高度为0.907 mm，SSE工况下最大起跳高度为1.21 mm，两者均远小于格架倾倒限值，因此格架在OBE和SSE工况下均不会发生倾倒。分析表明，格架在半载时倾倒限值最小，而在地震中格架半载时倾斜幅度最大，因而在地震中格架半载时发生倾倒的可能性最大，在格架倾倒分析中半载情况应作为重点工况进行验证。

图17 地震工况下格架在z方向的最大起跳高度曲线Fig.17 Maximum jump height curve of rack in z direction under seismic condition

4 结论

本研究基于超单元法建立了HFETR乏燃料贮存格架抗震分析模型，并基于APDL程序和VB平台对格架分析模型进行了参数化处理，建立了格架专用分析界面，解决了多格架抗震分析困难以及结构迭代设计流程复杂的问题，为格架结构的迭代设计及抗震分析提供了有力的支撑。利用参数化-超单元法对整池HFETR格架进行了抗震分析，对比分析了格架在空载、半载和满载3种情况下的表现。结果表明：格架在3种情况下应力均满足要求，不会发生碰撞和倾倒；在滑移分析和倾倒分析中，空载和半载工况需重点关注。