核电站高能管道冲击过程数值模拟研究

彭建， 刘浪， 石媛媛

(1.中广核工程有限公司 核电安全监控技术与装备国家重点实验室，广东 深圳 518172； 2.深圳中广核工程设计有限公司，广东 深圳 518172)

国际原子能机构(IAEA)和美国国家标准学会(ANSI)在20世纪80年代对核电站高能管道破裂防护设计开展了大量研究，编制了标准NS-G-1.11[1]和“ANSI/ANS58.2”[2]，明确要求核电站设计必须考虑管道甩动效应，提供适当的预防和缓解措施，以保证核安全不受到损害。如果核电厂高能管道断裂后，从破口释放出来的流体会对断裂管道产生反作用力，从而引起管道甩动，即管道甩动效应。甩击管道撞击附近的其他管道，它会产生管道对管道的冲击，这可能造成更严重的多米诺效应。管道甩击管道可分为2个阶段，即管道自由甩动和管道冲击管道阶段。Reid等[3-5]国外学者对管道自由甩动进行了大量的实验研究，设置了由储气罐、试验管道和爆破片组成的试验装置，创建了大挠度动力梁理论模型，并通过试验论证了模型的有效性。Alzheimer等[6]国外学者对管道冲击管道也开展了实验研究，首先基于美国某核电厂系统布置模型和压力温度情况，识别出100组潜在的管道甩击管道工况；然后采用弹射器装置模拟管道冲击管道试验；最后结论为甩击管道不会引起名义尺寸较大且壁厚也较大的靶物管道发生破裂，可能会引起名义尺寸较小的靶物管道(不论其壁厚如何)发生破裂，会造成名义尺寸相同或较大而壁厚相等或较薄的管道发生贯穿裂纹。上述结论广泛应用于工程设计，但变形机理有待进一步研究。Baum等[7]专家学者采用试验方法模拟了管道自由甩动和管道冲击管道的全过程，研究了其变形机理，结果发现靶物管道存在压缩和弯曲两种机理，支撑方式对靶物管道的变形影响较大。但是，由于试验基于1寸的小直径管道和理想刚性约束，因此对于大管径和柔性管道系统，还有待进一步研究。

国内标准EJ/T 335和NB/T 20516也提出核电站设计必须考虑管道甩击效应。学者对管道甩击也开展了研究，Yang等[8]采用试验和数值模拟的方法研究了冲击位置对靶物管道的变形模式的影响。但是，国内关于管道甩击的研究主要集中在防甩限制件方面，王春霖等[9-12]研究了防甩限制件的分析方法及影响因素。而管道抗冲击性能的研究主要集中在海洋管道[13-20]。杨秀娟等[13]采用非线性动态有限元法模拟了坠物撞击海底管道的过程，分析了物体形状、撞击角度、摩擦、内压、混凝土厚度、埋深及悬空长度与管道凹陷深度的关系。谢丽媛等[14]利用落锤实验和有限元模型得出海底管道在冲击载荷作用下的破坏模式为整体弯曲变形与冲击凹痕部位的局部弯曲耦合形成，管道钢材的屈服强度的增加可有效减小冲击作用下的局部凹陷变形。张荣等[15]通过落锤实验证明圆钢管失效形式均为剪切失效，随冲击能量增加圆钢管变形依次为弹性小变形阶段、塑性大变形阶段、剪切失效阶段等3个阶段。

甩击管道和靶物管道都可能产生较大的变形，并受到很多因素的影响，其中管道的管径比、径厚比、冲击能量、冲击位置、支撑方式已进行了研究。但是，冲击角度的影响暂未见文献报道。工程上，对于甩击管道一般没有要求，主要关注靶物管道的变形。因此,本文采用数值模拟方法，通过模拟管道自由甩动和冲击管道全过程，主要研究冲击角度对靶物管道变形的影响，并分析了约束方式和径厚比的影响。

1 管道冲击数值模型

1.1 数值算法

ABAQUS是一套功能强大的工程模拟有限元软件，包括ABAQUS/Standard和ABAQUS/Explicit分析模块。ABAQUS/Explicit在处理冲击问题方面非常成熟和方便，用户仅需要提供合理的结构参数、材料参数、边界条件和载荷工况这些工程数据就可以获得精确解。ABAQUS/Explicit采用中心差分法进行显式时间积分：

(1)

式中:M为质量矩阵；C为阻尼矩阵；K为刚度矩阵；u(t)为位移矩阵；F(t)为载荷矩阵；Δt为时间增量。

显式时间积分的求解精度取决于时间积分步长，需要非常小的时间步以保持稳定，而时间步取决于单元尺寸和波速；有稳定的时间增量的限制，完成指定分析一般需要较多的时间增量，不需要求解线性方程组，每个增量步的计算成本较低，适用于冲击等高频动力学问题，能分析各种复杂的接触问题。本文采用ABAQUS/Explicit模拟管道甩击管道，通过体积粘性考虑模型阻尼，体积粘性采用软件推荐值，线性参数取0.06，二次参数取1.2，最小稳定时间步长和接触算法通过验证分析确定。

1.2 材料模型

为描述金属材料的高速冲击效应，利用商业有限元软件ABAQUS建立管道甩击管道的数值模型。管道采用壳单元S4四边形单元进行网格划分，选择Simpson积分方法，厚度方向取5个积分点，单元应力位置在下表面；材料选用Mises弹塑性模型，真应力应变曲线根据ASME BPVC VIII的附录，并基于屈服强度、抗拉强度、弹性模型计算获得，本文采用核电站中的典型碳钢材料P265，其真应力应变曲线见图1。应变率效应采用经典的Cowper-Symonds模型：

图1 真应力应变曲线Fig.1 True Stress-strain curve

(2)

1.3 结构模型

结构模型由甩击管道和靶物管道组成(见图2)，甩击管道和靶物管道参数选择核电站中的典型参数。甩击管道参数见表1；靶物管道长度L3为4 m，外径为88.9 mm，壁厚、约束方式、角度见表2。甩击载荷F根据标准ANSI/ANS-58.2[2]计算获得：

图2 管道甩击管道模型示意Fig.2 Schematic diagram of pipe whipping pipe model

表1 甩击管道几何参数Table 1 Geometric parameters of whipping pipe

表2 靶物管道参数Table 2 Parameters of target pipe

F=CTPA

(3)

式中:CT为推力系数；P为管道正常运行工况压力；A为管道内截面面积。

P取核电站中典型压力15.5 MPa，CT取1.0，则F为92 822 N。甩击管道在锚固端采用全约束，在自由端采用壳载荷类型施加载荷F，该载荷类型的方向伴随单元方向转动。靶物管道的约束方式见表2。

2 模型验证

由于试验数据局限性，采用不同的试验分开验证管道自由甩动和冲击阶段的数值模型有效性。

2.1 管道自由甩动阶段

1) 实验模型

文献[3]开展了管道自由甩动试验，试验装置见图3。本实验主要由储气罐、试验管道和爆破片3个部分组成。储气罐模拟设备提供一个稳定的压力源，该储气罐足够大，整个实验过程中压力下降得非常小。爆破片模拟管道的双端剪切断裂，压力达到指定的限值，在1 ms内完成爆破，触发试验管道自由甩动，试验参数见表3。管道的运动使用高速胶片摄像机(Hadlands HYSPEED)进行记录，它提供了一系列管道甩动的瞬时剖面。仿真算法、材料模型和结构模型，基于表3的试验参数进行数值模拟。

图3 管道自由甩动装置Fig.3 Schematic diagram of pipe free whip experimental facility

表3 管道自由甩动参数Table 3 Parameters of pipe free whipping

2) 结果分析

通过改变网格尺寸开展最小稳定时间步长的敏感性分析，网格长度与管道周长的比值为1/30、1/15、1/10、1/6，分别对应最小稳定时间步长为5.0×10-7、1.1×10-6、1.7×10-6、2.2×10-6。图4为不同网格尺寸下管道自由端的变形时程图，从图可知，不同的时间步长，变形最大差异为2%，因此，后续自由甩动分析的最小稳定时间步长达到10-7的量级，变形结果的精度是足够的。图5为文献[3]中试验和ABAQUS数值模型的自由甩动变形结果对比图，从图可知，变形结果基本吻合。

图4 管道自由端的变形时程曲线Fig.4 Deformation-time curves at pipe free end

图5 试验和ABAQUS的甩动变形结果对比Fig.5 Comparison of whipping deformation between experiment and ABAQUS calculation

2.2 管道冲击管道阶段

1) 实验模型

文献[8]开展了管道冲击管道试验，试验装置见图6。本实验主要由弹簧弹射器、甩动管道和靶物管道3部分组成。弹簧弹射器提供甩动管道的能量，在弹簧作用下，甩动管道对靶物管道产生冲击作用，试验参数见表4。管道甩动和冲击过程采用高速摄像机(Kodak Vidicon)记录，记录速度为2 000帧/s。试验后按轴向和周向方向进行切片，测量了冲击段管道中心线的直径减小和弯曲角。采用第1章的仿真算法、材料模型和结构模型，基于表4的试验参数进行数值模拟。

图6 管道冲击管道装置Fig.6 Schematic diagram of pipe-to-pipe impact setup

表4 管道冲击管道参数Table 4 Parameters of pipe-to-pipe impact

2) 结果分析

通过改变网格尺寸开展最小稳定时间步长的敏感性分析，网格长度与管道周长的比值为1/33、1/22、1/13，分别对应最小稳定时间步长为6.5×10-7、8.5×10-7、1.4×10-6，接触约束选择拉格朗日乘子法，接触摩擦采用罚函数，摩擦系数取0.1。如图7所示，最小稳定时间步长达10-7的量级时，冲击截面的变形结果基本无差异。为了分析接触约束算法的影响，保持网格长度为1/33，分析了接触约束为罚函数的变形情况，如图7所示，与拉格朗日乘子法的变形结果基本一致。另外，还分析了接触摩擦的影响，从图7可知，摩擦因素对该数值模型也基本无影响。因此，对于该数值模型，接触约束算法和摩擦的影响可忽略，但最小稳定时间步长需达到10-7的量级。将文献[8]中试验数据和ABAQUS计算结果进行对比，分别从靶物管道轴向、径向变形进行分析，从图8可知，试验结果与数值模拟基本吻合。进一步定量分析发现，冲击截面管道直径变化量Δd的计算值和实验值分别为22.88 mm和23.7 mm，差异仅为3.5%(见图9)。

图7 冲击截面直径变化量的敏感性分析Fig.7 Sensitivity analysis of diameter reduction at impact section

图8 试验和ABAQUS的靶物管道变形结果对比Fig.8 Comparison of target pipe deformation between experiment and ABAQUS calculation

图9 靶物管道冲击截面直径变化图Fig.9 Diameter reduction of target pipe at the impact section

3 工程问题分析

3.1 甩击角度效应

采用以上经验证的数值方法和模型，研究不同甩击角度对靶物管道变形的影响。以靶物管道外径88.9 mm、壁厚11.13 mm、两端固定为例，工况1～4、11～12对甩击角度从30°～90°的6种工况进行了数值模拟。从图10和图11可知，靶物管道的变形可分为2种模式，压缩(即局部变形)和弯曲(即整体变形)，最后的变形状态是压缩和弯曲的耦合作用。压缩主要发生在冲击截面，可以通过冲击截面的直径变化量Δd表征。弯曲可以通过弯曲角θ表征(见图10)。图12展示了冲击截面Δd随时间的变化，冲击瞬间，冲击截面发生了较大的变形，然后随着θ的增加而增加，从而进一步说明了Δd和θ之间存在一定的相互作用关系。Alzheimer等[6]为了综合考虑压缩和弯曲效应，定义了一个变形参数Td来表征管道的变形行为：

图10 靶物管道弯曲应变Fig.10 Bending strain of target pipe

图11 靶物管道冲击截面压缩应变Fig.11 Crushing strain of target pipe at the impact section

图12 直径变化量时程曲线Fig.12 Diameter reduction-time curves

(4)

式中:d为靶物管道直径；dm为甩击管道直径；θ为弯曲角。

从图13可知，靶物管道冲击截面的直径变化量Δd和应变随着甩击角的变大而变大，这主要是因为随着甩击角度的变大，甩击管道和靶物管道的接触面积变小，从而导致冲击更加局域化。图14说明弯曲角θ仅仅在甩击角较小的30°时较小，其他甩击角情况下的弯曲角基本一致；变形参数Td随甩击角的变大而增大，说明其考虑了压缩效应的变化趋势。因此，甩击角度主要影响压缩变形，对弯曲变形不敏感，变形参数Td综合表征了靶物管道的变形特性，且在甩击角为90°时变形最大。

图13 直径变化量和应变随甩击角的变化Fig.13 Diameter reduction and strain vs. whipping angle

图14 弯曲角和变形参数随甩击角的变化Fig.14 Bending angle and deformation parameter vs. whipping angle

3.2 约束方式影响

为研究不同约束方式对靶物管道变形的影响，以靶物管道外径88.9 mm、壁厚11.13 mm、甩击角为90°，工况4～7对不同的约束方式进行了数值模拟。从表5中可知，约束方式对压缩基本没有影响，主要影响弯曲。工况5和6的弯曲角基本一致，而大于工况4，小于工况7，这是因为靶物管道轴向约束不同，工况4两端轴向都被约束，工况5和6只有一端轴向被约束，而工况7两端轴向都没有约束。变形参数与弯曲角的变化趋势一致。因此，约束方式主要通过管道轴向是否约束影响管道的变形，轴向约束对靶物管道的变形影响较大，且两端简支时变形最大。

表5 靶物管道的变形结果Table 5 Deformation results of target pipe

3.3 径厚比的影响

为了进一步研究靶物管道的变形特性，选择了更小壁厚的靶物管道进行研究。图15给出了工况8的靶物管道冲击截面的变形过程，随着时间的推移，靶物管道逐渐被压扁。塑性应变主要发生在冲击面，冲击背面的塑性应变非常小，冲击点两侧凸起，形成类似“心形”的截面形状，这可能主要是由于惯性作用导致。工况9和10的冲击截面变形与工况8基本一致，只是变形更加大，最后压成了饼状。

图16展示了冲击点的应变变化情况：冲击瞬间，应变发生突变，随后增长趋势变缓，在30～40 ms时段，应变基本无变化，这主要因为靶物管道两端约束处形成了塑性铰，变形模式由冲击截面压缩变成了整体弯曲模式；40～60 ms时段，由于管道两端轴向约束的作用，整体弯曲刚度大于冲击截面的压缩刚度，变形模式又变成了压缩为主，从而使冲击点的应变随时间的增加而增加。应变变化的趋势与冲击截面直径变化的趋势基本一致。这说明冲击截面直径的变化可以表征冲击截面的应变特性。

图16 工况8的直径变化量和应变时程曲线Fig.16 Diameter reduction and strain-time curves at test 8

4 结论

1) 靶物管道的变形包括压缩和弯曲2种变形模式，且压缩和弯曲是相互作用的，综合考虑了压缩和弯曲的变形参数可以很好的表征靶物管道的变形行为。

2) 靶物管道的变形随甩击角度的增加而增加，正交甩击时变形最大；甩击角度主要影响靶物管道的压缩变形，对弯曲变形影响较小。

3) 靶物管道的变形与轴向是否约束相关性大，轴向无约束时变形较大；轴向约束主要影响靶物管道的弯曲变形，对压缩变形不敏感。

4) 随着径厚比的增加，靶物管道变形越大；冲击截面直径变化量与应变随时间的变化趋势一致。