爆炸冲击波的提取和显示方法

李伟; 宁建国; 王成

doi:10.11883/1001-1455(2003)02-0178-5

摘要: 在差分格式中引入人工粘性捕捉激波,处理压力、密度等物理量的间断,使间断得到平滑,在计算结果中,冲击波就表现为物理量具有很大梯度的一个狭窄过渡区。利用边缘检测算子从计算数据中检测出位于边缘上的点,通过处理压力数据可以在图像中清晰地显示冲击波的轮廓。

恶意行径导致的土木工程灾害，越来越频繁。据统计资料显示，2005年至今，印尼、泰国、巴基斯坦、俄罗斯等国家接连发生的恐怖爆炸事件，造成了大量的人员伤亡和财产损失，其爆炸方式主要有直接命中爆炸、近空爆炸和间接爆炸等。而核电结构由爆炸可能引起的核泄漏将产生更严重的后果。目前，AP1000核岛厂房是世界上最先进的第三代核电技术^[1]，相较第二代厂房新增的屏蔽结构成为核厂房的第一层外部防线。研究它在爆炸荷载作用下的动力响应、破坏模式及抗爆性能，能为核厂房的安全评估和提高其抗冲击能力措施的制定提供理论依据，具有重要的军事价值和社会政治意义。

爆炸的物理传播过程极其复杂，由爆炸荷载引起的结构动力响应更复杂。很多学者对爆炸冲击波的传播演化规律和不同结构遭受爆炸冲击荷载时的动力响应进行了研究。Rogers等^[2]、Baker^[3]、Karpp等^[4]研究了爆炸冲击波在框架结构、空气、薄壁壳结构中的传播规律。Hashemi等^[5]、Tang等^[6]对桥梁结构在不同比例距离和不同炸药量爆炸荷载作用下的动力行为和最不利损伤部位进行了探讨；Kelliher等^[7]、Feng等^[8]对高层建筑在爆炸荷载作用下的动力响应和鲁棒性进行了研究；Chen等^[9]、Zhang等^[10]对混凝土重力坝在不同部位、不同介质中爆炸荷载作用下的损伤演化机理进行了研究，探寻最不利爆炸部位。而核岛厂房的兴起发展时间较其他结构晚，因此对新型核岛厂房在爆炸荷载作用下的研究不足。由于核岛厂房模型试验的花费巨大、数据采集困难等，使核岛厂房模型试验的发展有很大的局限性。而随着计算机的发展，对核岛厂房爆炸试验的研究开始从模型试验向数值模拟发展。因此，建立精细的新型核岛厂房数值模型，进行爆炸荷载作用下的动力响应分析，十分必要。

目前，在核电领域，王天运等^[11]根据核电站安全壳的结构形式，采用流固耦合算法，对它在爆炸荷载作用下的动力响应进行了数值模拟，得出了造成安全壳破损的最小安全距离和爆炸当量；Pandey等^[12]研究了安全壳在遭受外部不同距离爆炸时的非线性响应；Bao等^[13]用LS-DYNA非线性有限元软件，模拟了钢筋混凝土在近距离冲击荷载下的残余强度；Cao等^[14]对曲面板在冲击荷载下的动力响应进行了讨论，并对核电结构在冲击荷载下的损伤程度和损伤区域进行了预测；赵春风等^[15]采用流固耦合算法，对钢筋混凝土安全壳在内爆作用下的动力响应进行了数值模拟，证明在比例距离为3.258 m/kg^1/3、爆源距基底48 m的内爆作用下，数值模拟能保持较好的精度。综上所示，流固耦合技术日益成熟，可用于爆炸领域的研究。

但是，以前研究多集中于探讨安全壳在内部爆炸荷载作用下的动力响应，在研究外爆冲击时，均忽略了设备开口引起的结构非对称性的影响，对于新型核岛机型AP1000屏蔽厂房的抗爆能力的研究还处于空白。本文中，采用非线性动力有限元软件AUTODYN，通过炸板实验验证数值算法的可靠性，对屏蔽厂房在不同起爆位置接触爆炸荷载下的动力响应进行研究，探究最不利起爆部位，对厂房的动力响应、失效行为和损伤特性进行分析，并提出一种新的配筋方案来增强结构的抗爆能力。

1. 状态方程和材料参数

钢筋和混凝土材料在爆炸作用下，会经历高达10～10⁶ s⁻¹的应变率^[16]。在这种高应变率情况下，钢筋的强度能够提高50%，混凝土的抗压强度能够提高100%。同时，Batra等^[17]、Kalthoff等^[18-19]、Needleman等^[20]、Ravi-Chandar^[21]发现很多材料在高速荷载作用下形成了延展性剪切破坏。因此，本文中，同时考虑钢筋和混凝土在冲击荷载作用下的高应变率效应。

1.1 炸药

TNT炸药采用标准JWL状态方程描述：

$p = A\left( {1 - \frac{\omega }{{{R_1}V}}} \right){{\rm e} ^{ - {{{R}}_{{1}}}{{V}}}} + B\left( {1 - \frac{\omega }{{{R_2}V}}} \right){{\rm e} ^{ - {{{R}}_{\rm{2}}}{{V}}}} + \frac{{\omega E}}{V}$

(1)

式中：p为压力，V为相对体积，E为体积内能，A、B、R₁、R₂、ω为常数。

炸药的材料参数分别为^[22]：ρ=1 630 kg/m³，A=371.2 MPa，B=32.3 MPa，R₁=4.15，R₂=0.95，ω=0.3，E=8.9 GJ/m³。

1.2 空气

空气采用理想气体状态方程描述：

$p = {C_0} + {C_1}\mu + {C_2}{\mu ^2} + {C_3}{\mu ^3} + ({C_4} + {C_5}\mu + {C_6}{\mu ^2})E$

(2)

式中：C₀、C₁、C₂、C₃、C₄、C₅、C₆均为材料参数，E为初始体积内能。

空气的材料参数分别为^[22]：C₀=C₁=C₂=C₃=0，C₄=C₅=0.4，C ₆=0，E=2.5 MJ/m³，ρ=1.225 kg/m³。

1.3 混凝土

在侵蚀和爆炸荷载作用下，混凝土具有大变形、高应变率和高压行为。为此，Riedel等^[23]提出了RHT混凝土本构模型，包括强度模型和损伤模型，通常与p-α状态方程一起使用。

RHT的强度模型分为三个阶段，即弹性阶段（等效应力 $\bar \sigma {\text{≤}} {\sigma _{\rm{e}}}$ ）、线性强化阶段（ ${\sigma _{\rm{e}}} {\text{≤}} \bar \sigma {\text{≤}}{\sigma _{\rm{f}}}$ ）和损伤软化阶段，如图1所示。在材料的损伤软化过程中，屈服应力随着材料损伤的增加而降低，最终降至残余应力σ_r。失效应力σ_f、弹性极限应力σ_e和残余应力σ_r由当前的应力、偏应力和应变率决定。RHT本构通过引入三个极限面，即屈服面、失效面和残余强度面，分别描述了混凝土的初始屈服强度、失效强度及残余强度的变化规律。

图 1 RHT本构应力应变曲线

Figure 1. Stress-strain curve of RHT constitutive model

下载: 全尺寸图片幻灯片

失效面由标准化静水压力 ${p^*}$ 、Lode角 $\theta$ 和应变速率 $\dot \varepsilon$ 共同表达：

${Y_{{\rm{fail}}}}({p^*},\theta ,\dot \varepsilon ) = {Y_{\rm{c}} }\left( {{p^*}} \right){R_3}(\theta )f_{\rm{di}}(\dot \varepsilon )$

(3)

${Y_{\rm{c}}}\left( {{p^*}} \right) = {f'_{\rm{c}}}{A_{\rm{f}}}{({p^*} - {p^*_{{\rm{spall}}}}f_{\rm{di}}(\dot \varepsilon ))^N}$

(4)

式中：Y_fail(p^*,θ, $\dot \varepsilon$ )为失效面函数，Y_c(p^*)、R₃(θ)和f_di( $\dot \varepsilon$ )分别为压力相关项、Lode角相关项和应变率函数。 ${p^*} = p/f'_{\rm{c}}$ 、 $p_{{\rm{spall}}}^* = {f_{\rm{t}}}/{f_{\rm{c}}}$ 分别为标准静水压力和用单轴压缩强度f_c归一化的标准层裂强度，p为实际静水压力，θ、 ${\rm{\dot \varepsilon }}$ 分别为Lode角和应变率， $f'_{\rm c}$ 为归一化的标准压力，A_f和N分别为失效面常数和失效面指数。R₃(θ)为Lode角的相关函数，可以由下式表示：

${R_3}(\theta ) = \frac{{2(1 - Q_2^2)\cos \theta + (2{Q_2} - 1)\sqrt {4(1 - Q_2^2){{\cos }^2}\theta + 5Q_2^2 - 4{Q_2}} }}{{4(1 - Q_2^2){{\cos }^2}\theta + {{(1 - 2{Q_2})}^2}}}$

(5)

${Q_2} = {Q_0} + {B_{\rm{Q}}}{p^*}{\rm{ }}\begin{array}{*{20}{c}} {}&{0.51 {\text{≤}} {Q_2} {\text{≤}} 1.0} \end{array}$

(6)

$\theta = \frac{1}{3}\arccos \left(\frac{{3\sqrt 3 {J_3}}}{{2J_2^{1/2}}}\right){\rm{ }}\begin{array}{*{20}{c}} {}&{}&{0 {\text{≤}} \theta {\text{≤}} \dfrac{{\text{π}} }{3}} \end{array}$

(7)

式中：Q₀为材料的初试拉压强度比，B_Q为过度因子，J₂、J₃为偏应力张量的第二、第三不变量。对于混凝土的动态抗拉强度，Malvar等^[16]研究了大量的试验数据，并提出了 $f_{\rm{dit}}(\dot \varepsilon )$ 为：

$f_{\rm{dit}}(\dot \varepsilon ) = \frac{{{f_{\rm{td} }}}}{{{f_{\rm{t} 0}}}} = \left\{ {\begin{array}{*{20}{l}} {{{\left( {\dot \varepsilon /{{\dot \varepsilon }_{{\rm t}0}}} \right)}^\gamma }}&\quad\quad{{\rm{ }}\dot \varepsilon {\text{≤}} 1\;{{\rm s}^{ - 1}}} \\ {\beta {{\left( {\dot \varepsilon /{{\dot \varepsilon }_{\rm{t} 0}}} \right)}^{1/3}}}&\quad\quad{{\rm{ }}\dot \varepsilon {\text{＞}} 1\;{{\rm s}^{ - 1}}} \end{array}} \right.$

(8)

式中：f_td为应变率为 $\dot \varepsilon$ 时的动力抗拉强度，f_t0为应变率为 ${\dot \varepsilon _{{\rm{t}}0}}$ 时的静力抗拉强度， ${\dot \varepsilon _{{\rm{t}}0}}$ =3.0×10⁻⁶ s⁻¹。其中，logβ=6γ−2，γ=1/(1+8f_c/f_c0)，f_c是混凝土的静态单轴抗压强度，f_c0=10 MPa。

对于混凝土的动态抗压强度，有：

$f_{\rm{dic}}(\dot \varepsilon ) = \frac{{{f_{{\rm{cd}}}}}}{{{f_{{\rm{c0}}}}}} = \left\{ {\begin{array}{*{20}{l}} {{{\left( {\dot \varepsilon /{{\dot \varepsilon }_{{\rm{c0}}}}} \right)}^{1.026\zeta }}}&\quad\quad{{\rm{ }}\dot \varepsilon {\text{≤}} 30\;{{\rm s}^{ - 1}}} \\ {\delta {{\left( {\dot \varepsilon /{{\dot \varepsilon }_{{\rm{c0}}}}} \right)}^{1/3}}}&\quad\quad{{\rm{ }}\dot \varepsilon {\text{＞}} 30\;{{\rm s}^{ - 1}}} \end{array}} \right.$

(9)

式中：f_cd为应变率为 $\dot \varepsilon$ 时的动力抗压强度，f_c0为应变率为 ${\dot \varepsilon _{{\rm{c}}0}}$ 时的静力抗压强度， ${\dot \varepsilon _{{\rm{c}}0}}$ =30×10⁻⁶ s⁻¹。其中，logδ=6.156ζ−0.49，ζ=(5+3f_cu)⁻¹，f_cu为混凝土的静力抗压强度。

RHT本构通过引入损伤变量D_a描述残余失效面：

${D_{\rm{a}}} = \sum {\frac{{\Delta {\varepsilon _{\rm{p}}}}}{{\varepsilon _{\rm{p}}^{{\rm{fail}}}}}}$

(10)

$\varepsilon _{\rm{p}}^{{\rm{fail}}} = {D_1}{({p^*} - p_{{\rm{spall}}}^{\rm{*}})^{{D_2}}} {\text{≥}} \varepsilon _{\rm{f}}^{{\rm{min}}}$

(11)

式中：∆ε_p为塑性应变， $\varepsilon _{\rm{p}}^{{\rm{fail}}}$ 为失效应变，D₁和D₂为损伤常数， $\varepsilon _{\rm{f}}^{{\rm{min}}}$ 为最小失效应变。Hulatt等^[24]指出，在辐射防护领域，使用重混凝土可以有效地屏蔽各种射线对环境的影响。因此，在核岛研究中，混凝土密度一般大于2 500 kg/m³。

混凝土RHT本构模型的材料参数分别为：密度2 750 kg/m³，泊松比0.167，过度因子1.220，参考温度300 K，比热654 J /(kg·K)，损伤因子D₁=0.040，损伤因子D₂=1.000，最小失效应变0.001，弹性强度/抗拉强度=0.700，弹性强度/抗压强度=0.530，侵蚀应变0.002，剪切模量16.700 GPa，抗压强度40.000 MPa，抗拉强度/抗压强度=0.100，抗剪强度/抗压强度=0.180，完整失效面常数A_f=1.920，完整失效面指数0.760，子午线比率0.690，断裂强度常数0.700，断裂强度指数0.800，压缩应变率0.032，拉伸应变率0.036。

1.4 钢筋

钢筋采用Johnson-Cook材料模型，该模型适合描述材料承受大应变、高应变率和高温下的强度性能，这些性能可能出现在由于高速碰撞或者炸药爆炸引起的强烈冲击荷载问题中。在爆炸载荷作用下，热传导时间远大于爆轰持续时间，且由于变形热导致的温升不大，材料出现强化继而接近失效时，应力减小的幅度较小，所以温度软化效应不明显。材料的屈服和硬化模型为：

$\sigma = (A' + B'\varepsilon _{\rm{p}}^n)[1 + C'\ln({\dot \varepsilon _{\rm{p}}}/{\dot \varepsilon _0})]$

(12)

式中：ε_p为等效塑性应变， $\dot \varepsilon _{\rm p}$ 为等效塑性应变率， $\dot \varepsilon _0$ 为参考塑性应变率， $A'$ 为材料的屈服强度， $B'\varepsilon _{\rm{p}}^n$ 为对材料强化段的描述， $C'$ 为应变率敏感因数。

钢筋材料参数分别为^[25]：密度7 830 kg/m³， $E'$ =200 MPa，μ=0.28， $A'$ =300 MPa， $B'$ =384 MPa，n=0.26， $C'$ =0.014， ${{\rm{\dot \varepsilon }}_0}$ =1 s⁻¹，失效应变0.1。

2. CEL算法和模型验证

2.1 CEL算法

爆炸问题的计算，一般采用拉格朗日（Lagrange）算法和欧拉（Euler）算法。在Lagrange算法中，材料附在网格上，可以追踪每个质量点的运动，所以网格和材料一起变形，因而这种算法适用于追踪小变形、大位移的材料运动，即容易确定时间历程、材料和结构界面及其内部的应力应变状态，且与其他算法相比计算速度较快，可以较好模拟固体材料行为。在Euler算法中，自由边界面和材料的交界面通过固定的欧拉网格表达，材料可通过固定的网格边界流进流出，所以大变形或有流动的情形不会导致网格畸变，适合描述液体和气体的行为，但由于必须采用复杂的算法追踪材料的运动，计算效率较低下。在数值模拟流固耦合及大变形的问题时，可将欧拉单元和拉格朗日单元耦合以发挥两种算法的优势，即耦合的欧拉-拉格朗日（coupled Euler-Lagrange，CEL）算法^[26-27]。本文中，采用Lagrange算法模拟钢筋混凝土板，采用Euler算法模拟炸药和空气。

2.2 模型验证

由于核电屏蔽厂房结构的复杂性、大型化、数据信息采集的难度和误差大，在现有条件下很难进行物理实验，并且实验具有一定的局限性，所以目前的研究多通过数值计算对结构进行爆炸动力响应分析。为了验证分析屏蔽厂房在爆炸作用下动力响应及抗爆措施的数值模型能否有效重演破坏现象，本文中对Li等^[28]的钢筋混凝土接触爆炸试验进行数值计算。

运用AUTODYN软件，采用CEL算法，严格参照实验条件，通过分离式配筋建立三维钢筋混凝土板模型。图2为混凝土板爆炸冲击试验的现场安装和布置图，比对冲击试验，取炸药量为0.1 kg，钢筋混凝土板尺寸为2 000 mm×800 mm×120 mm，板用9 $\varnothing$ 12 mm的纵向钢筋和11 $\varnothing$ 8 mm@200 mm的箍筋加固。将0.1 kg的TNT炸药放置在板的表面上。先用起重机将钢筋混凝土板放在平台上，再用角钢夹板将板的两端用螺栓固定。图3为全尺寸的三维数值模型，采用分离式配筋方法建立钢筋混凝土板。采用CEL算法，混凝土、钢筋采用Lagrange网格，空气、炸药采用Euler网格。其中，炸药直径d=40 mm，高度h=46 mm。混凝土板和钢筋网格大小均为5 mm，单元数量分别为3 938 221个和5 200个。与炸药接触的表面为迎爆面，相对的表面为背爆面。

图 2 钢筋混凝土板

Figure 2. Reinforced concrete slab

下载: 全尺寸图片幻灯片

图 3 三维数值模型

Figure 3. Three-dimensional numerical model

下载: 全尺寸图片幻灯片

图4为钢筋混凝土板迎爆面和背爆面的实验现象和数值模型的计算结果，损伤值由0（蓝：从未损坏）到1（红：完全损坏）逐渐增大。试验结果表明，在迎爆面有许多径向裂纹和环向裂纹，在背爆面有混凝土脱落。这是混凝土剥落和穿孔破坏的典型结果，因为在接触爆炸荷载作用下，迎爆面主要产生剪切破坏，背爆面主要产生弯曲破坏。通过测量，混凝土爆坑直径为20 cm，剥落直径为33 cm。由于混凝土材料具有较大的抗压能力和较弱的抗拉能力，因此背爆面的损伤大于迎爆面的损伤。数值结果与试验结果保持高度的一致性。混凝土爆坑和剥落的直径分别为19和31.5 cm，数值计算结果的混凝土爆坑和剥落直径误差分别为5%和4.5%。比较结果说明，钢筋混凝土板在接触爆炸荷载作用下的动力响应用数值模拟的方式研究是可行的。文献[29-31]中验证了这种数值方法也适用于曲壳，但在以往的研究中，还没有清楚地解释结构曲率对冲击载荷下结构损伤程度的影响。因此，本文采用的数值模型和CEL算法能够有效地再现AP1000核岛厂房在爆炸冲击荷载作用下的失效过程。

图 4 炸药起爆冲击作用下钢筋混凝土板的试验和数值结果

Figure 4. Experimental and numerical results of reinforce concrete slab under explosion

下载: 全尺寸图片幻灯片

3. 非线性动力响应分析

3.1 数值模型建立

核岛厂房由屏蔽厂房、安全壳结构和辅助厂房构成。屏蔽厂房作为核厂房的第一道外部防线，是本文主要的研究对象。图5(a)为AP1000型核岛厂房的构成。屏蔽厂房半径22.1 m，壁厚0.92 m，整体高度82.95 m。采用分离式配筋建立精细化的三维数值模型。空气模型的尺寸为50 m×50 m×83 m，可以完全包裹住核电厂房。在爆炸荷载作用下，屏蔽厂房的破坏往往在很短的时间发生，因此暂不考虑地基与厂房之间的相互作用及地基能量反射的影响，假设它坐落在刚性基础上，因此核岛厂房的底部施加全约束。此外，为了防止爆炸冲击波产生的能量在有限元模型空气的边界处产生折射、反射等现象，在空气模型的四周施加了流出边界条件模拟辐射条件。采用CEL算法，混凝土、钢筋采用Lagrange网格，空气、炸药采用Euler网格，数值计算模型如图5(b)所示。

图 5 AP1000核岛厂房有限元模型

Figure 5. Finite element model of AP1000 nuclear power plant

下载: 全尺寸图片幻灯片

钢筋混凝土承受动力荷载时，他们之间的黏结滑移现象明显且不能忽略，因此运用AUTODYN的滑动接触模型对他们之间的黏结滑移进行模拟。也就是说，在一定的黏结力限制内，钢筋随着混凝土的变形而变形，且黏结力与变形成正比关系。当黏结力超过这个限制时，钢筋与混凝土即会分离。钢筋与混凝土的最大黏结力取为18 MPa。

3.2 网格敏感性分析

在爆炸冲击波传播及其与结构发生相互作用时，数值模拟的精度在很大程度上取决于结构网格尺寸。因此，在数值模拟前，有必要进行网格尺寸的收敛性分析。对结构网格采取300、500、800 mm三种方式进行剖分，用模态分析对结果进行对比。由表1可以看到，500、800 mm网格相较于300 mm网格的一阶模态相差1.1%和1.6%，差别微乎其微。为了进一步说明网格尺寸的合理性，验证爆炸荷载起爆后到达目标点的超压峰值对网格的敏感度。图6为三种网格尺寸下，在位置Ⅰ-1起爆时，目标点（测点5）的超压曲线。由图6可以看出，500、800 mm网格较300 mm网格的超压峰值最大值相差0.9%和6.3%，超压最大值到达时间分别晚0.125和0.625 s。因此，在设备孔、空气孔和爆源附近等容易产生应力集中的区域划分300 mm网格，其他区域采用500 mm网格，在提高计算效率的同时，可以达到足够的计算精度要求。

表 1 不同网格尺寸下核岛厂房的自振频率

Table 1. Frequency of nuclear island under different element sizes

网格尺寸/mm	频率/Hz
网格尺寸/mm	一阶	二阶	三阶
300	3.316 9	3.326 5	5.141 0
500	3.366 6	3.377 1	5.191 8
800	3.372 0	3.382 4	5.184 8

下载: 导出CSV

| 显示表格

图 6 不同网格尺寸下测点5的超压峰值曲线

Figure 6. Peak overpressure curves at point 5under different element sizes

下载: 全尺寸图片幻灯片

3.3 爆源位置和测点

考虑屏蔽厂房受爆炸冲击荷载作用位置的随机性和厂房结构本身的非对称性，对爆源位置处于不同轴向高度和环向角度，进行了多组数值计算。

为便于描述工况和解释现象，对屏蔽厂房各区域按轴向高度、环向初始角度进行定义，如图7所示。在轴向，选取厂房的5个典型部位作为爆源高度：（Ⅰ）设备孔轴向中间部位（20.75 m）；（Ⅱ）圆柱高2/3 处（42.5 m）；（Ⅲ）空气孔下方5 m处（55 m）；（Ⅳ）锥形屋面中心1/2 处（68 m）；（Ⅴ）水箱中心处（76.75 m）。考虑屏蔽厂房设备孔的存在而导致的厂房非对称性，在4个环向部位用1（0°）、2（90°）、3（180°）、4（270°）表示。TNT炸药均建立在屏蔽厂房对应部位的外侧。对爆源附近区域的网格进行局部加密。表2为20种工况的组合方式。在空气模型边界采用透射边界条件，厂房底部采用固定约束保证底部无刚体位移。每种工况下，在爆源轴向布置12个测点，测点布设如图8所示。

表 2 爆炸位置和工况

Table 2. Locations of explosion and working conditions

爆炸位置	工况	爆炸位置	工况	爆炸位置	工况	爆炸位置	工况	爆炸位置	工况
Ⅰ-1	1	Ⅱ-1	5	Ⅲ-1	9	Ⅳ-1	13	Ⅴ-1	17
Ⅰ-2	2	Ⅱ-2	6	Ⅲ-2	10	Ⅳ-2	14	Ⅴ-2	18
Ⅰ-3	3	Ⅱ-3	7	Ⅲ-3	11	Ⅳ-3	15	Ⅴ-3	19
Ⅰ-4	4	Ⅱ-4	8	Ⅲ-4	12	Ⅳ-4	16	Ⅴ-4	20

下载: 导出CSV

| 显示表格

图 7 轴向、环向的炸点布设

Figure 7. Layout of axial and circumferential explosion points

下载: 全尺寸图片幻灯片

图 8 测点布置

Figure 8. Arrangement of observation points

下载: 全尺寸图片幻灯片

选取TNT等效爆炸当量为815 kg TNT，相当于文献[32]中某导弹装药质量的2～3倍，对每种工况进行接触爆炸数值模拟。为了简化计算，假设：（1）整个屏蔽厂房的材料是各向同性的均匀材料，底部为固定；（2）不计厂房和其内部设备的相互影响。

3.4 冲击波压力演化分析

以工况1为例，进行压力演化分析，各测点的超压曲线如图9所示。由图9可见，到达正超压峰值时间的先后顺序依次为测点3、4、2、5、1、6、7、8、12、9、10、11，压力幅值大小依次为测点3、4、2、5、1、6、7、8、9、10、11、12。测点3的压力在发生爆炸后很短的时间内上升到最大值，其他各测点到达正超压的时刻随着偏离爆源位置距离的增加而有所延迟，压力峰值随着偏离爆源处距离的增加而减小。测点12到达压力峰值时间在测点9前，这是由于特殊的结构位置造成的超压峰值提前。同时也发现，各测点正超压持续时间随着偏离爆源处距离的增加而有所延长。

图 9 工况1各测点的压力曲线

Figure 9. Pressure curves of observation points of working condition 1

下载: 全尺寸图片幻灯片

图10为屏蔽厂房在工况1下的应力波传播过程。当t=1 ms时，炸药在极短时间释放了大量能量，爆源处的压力迅速上升，产生强烈的应力波，正超压峰值已超过混凝土的动态抗压强度；当t=3 ms时，冲击波所携带的能量急剧衰弱，作用在厂房的单位面积能量衰弱，压力幅值沿着设备孔的方向偏移；当t=7 ms时，爆源中心区域空气大量排出，压力降低至负值，该区域开始受拉力作用；当t=12 ms时，结构的压力幅值减小，压力幅值继续向孔洞部位偏移，但未达到致其损伤的强度；当t=30 ms时，压力幅值继续降低，压力幅值除在爆源附近集中，开始在空气孔和设备孔的区域集中。

图 10 工况1的压力分布

Figure 10. Pressure distributions of working condition 1

下载: 全尺寸图片幻灯片

3.5 损伤演化分析

图11为工况1的损伤演化云图（损伤值由0（蓝）到1（红）逐渐增大）。为了讨论损伤的程度，以结构贯穿破坏的损伤质量作为评价指标。贯穿破坏表示单元的损伤值达到1，而结构贯穿破坏的损伤质量，可以通过有限元模型中所有损伤值达到1的单元的质量和计算。当t=5 ms时，爆炸产物急剧膨胀，瞬间形成的空腔充满了高温高压气体，形成强烈的应力波，对爆源区域的钢筋混凝土挤压，形成了椭圆破碎区，即所谓爆坑，此时钢筋混凝土失去了承载能力，结构主要破损区域的形状基本形成。爆坑形成的椭圆区域以环向为长轴，因为屏蔽结构的环向惯性矩大于轴向惯性矩，相同能量产生的轴向应力大于环向应力，所以爆源区域破坏以轴向拉压为主。当t=20 ms时，应力波继续向周围钢筋混凝土介质传播，此时爆炸冲击波的能量大部分消耗于混凝土的压缩和粉碎，作用在混凝土的单位面积能量迅速降低，但仍继续扩大损伤区域，至此，结构完全损伤区域基本形成。在t=20 ms到t=50 ms间，厂房的主要破碎区域不再继续扩展。这说明，冲击波在传播过程中急剧衰减，应力波变成弱的压缩应力波，其强度已低于混凝土的动抗压强度，虽不能直接压碎混凝土，但其负应力可以引起外围钢筋混凝土层产生径向扩张和切向拉伸应变，使厂房形成浅层破坏区。爆坑在竖直向的损伤最大长度为8.295 m，环向损伤最大长度12.32 m。

图 11 工况1的损伤分布

Figure 11. Damage distributions of working condition 1

下载: 全尺寸图片幻灯片

图12为工况1的损伤质量曲线，结构在爆炸发生后3和20 ms时，损伤质量达到稳定时损伤质量的88%和98%。图13为工况1时，结构底部测点1、爆源附近测点4和5、圆筒折坡处测点8、水箱折坡处测点9、结构顶部测点12的位移、速度、加速度曲线。位移变形只在爆源附近测点4、5处呈线性增长，发生塑性变形，在其他各测点均产生1 mm左右的可恢复变形；距离爆源最近的测点4、5的速度、加速度幅值远大于其他各测点，测点4的速度和加速度幅值分别为0.364 m/s和100 m/s²。

图 12 工况1的损伤质量曲线

Figure 12. Damage mass curve of working condition 1

下载: 全尺寸图片幻灯片

图 13 工况1各测点的位移、速度、加速度曲线

Figure 13. Displacement, velocity and acceleration curves of different points of working condition 1

下载: 全尺寸图片幻灯片

通过以上的直观定性描述和定量分析，得知整个爆炸过程在较短的时间内结构的贯穿破坏损伤质量趋于稳定，说明接触爆炸产生了局部的冲击破坏。通过对损伤演化过程和破碎区形状的分析，可以发现，当采用增加配筋的方式提高屏蔽厂房的抗爆性能时，纵向配筋的效果强于环向配筋的效果。

3.6 不同爆源位置的破坏程度和模式

3.6.1 相同轴向高度、不同环向角度

对爆源位置在相同轴向高度、不同环向角度的工况进行分析。图14为在爆炸发生50 ms后20种工况的结构贯穿破坏的损伤质量曲线。当轴向工况为Ⅱ、Ⅲ、Ⅳ、Ⅴ，而环向工况不同时，损伤质量变化很小，最大损伤质量变化不足2%；当轴向工况为Ⅰ时，不同环向角度的损伤质量最大变化为21.5%。环向工况2的损伤质量最小，在该工况下，炸药爆炸后迅速膨胀的应力波在初始的传播过程中很快到达设备孔区域，其耗散的能量远大于其他环向工况。

图 14 不同轴向高度、不同环向角度爆源位置的损伤质量曲线

Figure 14. Damage mass curves with different axial and circumferential explosion points

下载: 全尺寸图片幻灯片

限于篇幅，只给出在轴向高度Ⅱ、不同环向角度起爆后，结构在能量稳定后的损伤云图，如图15所示。可以直观地对比发现，在该轴向高度、不同环向角度下，最终的贯穿破坏区域、浅层破坏区域的形状和损伤质量均有很大程度的相似性。其他轴向高度同样符合这个规律。

图 15 不同环向角度爆源位置的损伤云图

Figure 15. Damage distributions with different circumferential explosion points

下载: 全尺寸图片幻灯片

综上所述，屏蔽厂房在接触爆炸荷载作用下，除了爆源位置在设备孔高度区的不同环向角度产生的损伤程度差异较大，在其他轴向高度、不同环向角度产生的损伤程度差异很小，当应力波传播到设备孔的位置时，能量发生耗散，相应的破碎区损伤质量减小。

3.6.2 不同轴向高度、相同环向角度

现在讨论不同轴向高度、相同环向角度的损伤模式差异（以环向角度为0°的工况1、5、9、13、17为例）。由于上节得出，相同轴向高度、不同环向角度的损伤差异较小，所以其他环向角度工况规律可类比。图16为各工况在1、5、10 ms三个时刻的压力和损伤演化云图。

图 16 不同轴向高度爆源位置的压力和损伤分布

Figure 16. Pressure and damage distributions with different axial explosion points

下载: 全尺寸图片幻灯片

当t=1 ms时，应力波集中于爆源区域，且压力以爆源为水平轴对称，损伤区域均呈以环向为长轴的椭圆形状，不同轴向工况的应力波特点无明显差异，结构破坏以轴向拉压破坏为主。

当t=5 ms时：工况1的应力波以水平方向为对称轴对称，工况5的应力波不再对称，在爆源偏上部有明显的受拉区域；工况9的正超压幅值出现在离爆源区域最近的设备孔区域，爆源中心承受负压作用；工况13的正超压峰值出现在爆源区域，在爆源处和水平成45°方向呈现对称的负压，负超压值出现在水箱轴向中心处；工况17的应力波幅值出现在爆源区域，且断裂面与垂直方向成45°，这是由于剪力在结构破坏中的影响逐渐增大导致的，表现出了明显的剪切破坏模式。此时工况1、5、9、13的损伤区域质量基本稳定。

当t=10 ms时，工况17的损伤质量继续增加，且裂纹方向沿着轴向扩展，究其原因是当爆炸物位于水箱中心区域时，应力波传播到水箱上下边界时发生反射，轴向方向应力波震荡剧烈，扩散受阻，在轴向形成了二次损伤。

由不同轴向高度起爆的损伤演化特性可以得出：当采用增加配筋率的方式提高结构的抗爆性能时，在圆柱区域和楔形体区域应重点提高纵向配筋率；在水箱区域起爆，由于水箱上下边界阻碍应力波的耗散，会产生二次损伤，应同时提高纵向配筋率和环向配筋率。

由上节，接触爆炸位置在相同轴向高度、不同环向角度时，损伤程度差异不大，所以该部分只给出不同轴向高度、相同环向角度1爆源的von Mises应力曲线，如图17所示。可以发现，在接触爆炸荷载作用下，von Mises应力均呈现双峰值现象，这是由结构在受应力波的正压和负压状态下呈现不同的响应导致的。爆源处von Mises应力峰值由大到小依次为工况17、9、13、5、1。工况17的von Mises应力远大于其他轴向工况，是因为应力波在该区域的传播过程中，传播到水箱上下边界时发生反射，应力波在爆源区震荡剧烈，难以耗散导致的。工况9的von Mises幅值＞工况5的von Mises幅值＞工况1的von Mises幅值，是由于结构底部施加全约束，结构呈以底部固支的悬臂梁特性。而工况13的爆炸部位相对特殊的几何形状，在探寻规律时与其他工况不具有可比性，单从幅值角度看，工况13的von Mises应力幅值介于工况5和工况9的von Mises应力幅值之间。

图 17 不同轴向高度爆源位置的von Mises应力曲线

Figure 17. Von Mises stress curves with different axial explosion points

下载: 全尺寸图片幻灯片

图18(a)为工况1、5、9、13、17在爆炸冲击波稳定后的损伤质量曲线。可以发现，接触爆炸发生在圆柱区域的损伤程度大于发生在楔形体和水箱区域的损伤程度。在各轴向工况中，接触爆炸发生在轴向工况Ⅱ的损伤程度最严重。此外，在爆炸发生20 ms后，轴向工况为Ⅰ、Ⅱ、Ⅲ、Ⅳ的结构破碎区损伤质量均趋于稳定。由图18(b)可以发现，轴向工况Ⅴ的结构破碎区损伤质量还在持续增长，原因如前。

图 18 不同轴向高度爆源位置的损伤质量曲线

Figure 18. Damage mass curves with different axial explosion points

下载: 全尺寸图片幻灯片

综上所述，接触爆炸发生在屏蔽厂房不同轴向高度造成的损伤程度有明显的差异。在轴向高度四种工况中，接触爆炸发生在轴向工况Ⅱ即42.5 m高度处造成的结构损伤程度最大。由于应力波在传播过程中遇到孔洞部位时，能量发生耗散，可以推断相同质量TNT的接触爆炸发生在轴向高度位于设备孔和空气孔之间的某个位置，屏蔽厂房破碎区的损伤质量最大。

4. 结　论

基于显示非线性动力分析程序AUTODYN，对AP1000核岛厂房在接触爆炸荷载作用下的动力响应及损伤行为进行分析，得到以下结论。

（1）AP1000核岛屏蔽厂房在接触爆炸的高超压、短持时的荷载作用下，结构发生局部破坏，结构的损伤质量在爆炸发生后较短的时间内趋于稳定。屏蔽厂房任意部位在遭受815 kg TNT接触爆炸荷载作用下，爆源区域失去承载能力，第一层防线失效。

（2）爆源位置在相同轴向高度、不同环向角度产生的损伤程度差异很小，在制定相应的抗爆措施时可以不考虑结构的非对称性。

（3）接触爆炸发生在圆柱区域的损伤程度大于发生在楔形体和水箱区域的损伤程度，且发生在轴向高度趋于设备孔和空气孔之间的位置，屏蔽厂房破碎区的损伤程度最严重。

（4）从结构的损伤演化和失效模式可以得出，当采用增加配筋率的方式提高屏蔽厂房的抗爆性能时，在圆柱区域、楔形体区域应重点增强纵向配筋率，在水箱区域，应同时增强纵向配筋率和环向配筋率。