侧爆作用下土中竖井结构弹塑性变形特征分析及计算

摘要

针对土中钢板‒钢筋混凝土竖井结构在常规武器爆炸作用下塑性铰线形成前的动力响应问题，综合采用模型试验、有限元数值仿真方法，分析了竖井结构中远距离爆炸作用下的变形过程及特征，确定了以量纲一环向相对位移${\alpha }_{\mathrm{1}}$作为结构变形程度的判据。推导了土中竖井结构在常规武器爆炸作用下的动力响应理论解，分析了按照弹性理论计算的量纲一环向相对位移${\alpha }_{\mathrm{1}\mathrm{e}}$和有限元计算得到的弹塑性量纲一环向相对位移${\alpha }_{\mathrm{1}\mathrm{p}}$的关系及其比值${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$的影响因素，构建了基于弹性解的弹塑性量纲一环向相对位移${\alpha }_{\mathrm{1}\mathrm{p}}$工程算法。提出了以量纲一环向相对位移${\alpha }_{\mathrm{1}}$作为结构变形程度的判据，并确定了结构闭环塑性铰线出现、破坏模式改变的阈值为${\alpha }_{\mathrm{1}}=\mathrm{3}\mathrm{\%}$；对比了理论与数值模拟计算结果，验证了提出的理论计算方法可以较好地计算钢板‒混凝土竖井组合结构在常规武器爆炸作用下弹性阶段的动力响应；结构长径比、结构缩尺比例、装药量大小对${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$随${\alpha }_{\mathrm{1}\mathrm{e}}$的变化规律的影响可以忽略；等效厚径比、钢板与混凝土厚度比较大及结构的相对高度较小时，计算得到的${\alpha }_{\mathrm{1}\mathrm{e}}$偏小，${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$较大。混凝土强度较大时，${\alpha }_{\mathrm{1}\mathrm{e}}$相同情况下结构塑性变形较小，${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$较小。

关键词

竖井结构; 爆炸作用; 弹塑性变形; 壳体振动理论

地下防护工程中常常出现竖井结构形式，如地下储藏设施、导弹发射井、防护工程应急出入口等。分析竖井结构的动力响应特征并对结构的变形进行定量计算，对竖井结构抗爆设计计算及武器打击下竖井结构的毁伤评估具有重要意义。

大量研究表明，结构特征变形量准确、快速的定量计算是防护工程快速毁伤评估的关键^［1-5］。针对地下结构受爆炸作用的动力响应计算，通常采用土与结构相互作用的分析方法^［6-8］。李国豪等^［9］、林润德等^［10-11］对非饱和土中应力波传播及土与结构的相互作用、结构的弹塑性动力响应进行了推导，钱七虎等^［12-13］考虑结构及地表卸荷，推导了核爆作用下箱型结构与饱和土、非饱和土的相互作用计算方法及结构变形计算方法，并对比了梁板模型、刚体模型、柔体模型的计算误差。对于地下壳体结构，冯淑芳^［14］、冷星^［15］忽略结构自身柔性，考虑岩土介质的弹簧效应，计算了拱结构在各种爆炸冲击波作用下的动力响应。、张宇等^［16］、陈海龙等^［17］考虑结构卸载及土体耗能，推导了拱结构在爆炸作用下反射系数及动力响应。

由此可见，对常规武器作用下土中钢板‒钢筋混凝土组合竖井结构的动力响应问题，特别是竖井结构弹塑性变形的定量计算方法研究较少。本文以某典型竖井结构为研究对象，对结构弹塑性变形特征进行了试验及有限元计算分析，并通过理论推导、参数分析的方法给出了竖井结构闭环塑性铰线出现前弹塑性响应阶段变形表征量的计算方法。这对常规武器作用下土中钢板‒钢筋混凝土组合竖井结构破坏模式和破坏程度的判别具有重要意义。

1 土中竖井结构动力响应特征分析

1.1　有限元计算及试验验证

1.1.1　试验设置

以某典型防护工程竖井结构为对象，开展了两种缩比尺寸靶体的外场爆炸试验。竖井结构可以简化为底部固端、顶部自由的土中圆柱壳结构。试验弹体、靶体及测点布置如图1所示。

图1  试验设置及测点布置

Fig.1  Experimental setting and arrangement of measuring point

图1中R为装药点距结构的法向距离，h_s为竖井结构钢板厚度，h_c为竖井结构混凝土厚度，r₁为井壁外径，r₂为井壁内径，L为竖井结构高度，a₀为装药点距结构固端截面的距离，φ₀为装药点到圆柱壳结构切线交点（C0点）与迎爆点（C1点）的角度，h_e为装药点埋深。图1中C1~C5为试验中的5个测点，其中C2、C3距C1的距离均为0.5 m。每个测点均布设土压传感器、加速度传感器及速度传感器。C0点为装药中心深度平面内装药到竖井结构切线交点。下文中竖井结构均以柱坐标系（x， φ， r）进行描述，以迎爆面固端支座位置为原点（x=0， φ=0），径向以指向圆心为正。两种缩尺靶体尺寸如表1所示。

表1  竖井结构主要尺寸

Tab.1  Main dimension parameters of shaft structure

靶体	L/m	hc /m	hs /m	r1 /m	（a0+he）/m
靶体1	3.17	0.070	0.001 2	0.366	4.91
靶体2	5.28	0.116	0.002 0	0.610	8.18

靶体混凝土采用C50混凝土，轴向、环向分别布置两层HRB335受力钢筋，环向、轴向钢筋配筋率均为1.5%。钢板屈服强度为235 MPa。试验工况设置如表2所示。试验中采用明挖开坑后支模浇筑混凝土，待混凝土养护结束后安置传感器，而后回填土的施工工序，使用PVC管预留出装药通道。试验过程如图2所示。

表2  试验工况设置

Tab.2  Setting of experimental conditions

炮次	1	2	3	4	5	6	7	8	9
靶体类型	1	1	1	1	1	1	1	2	2
装药量/kg	0.93	0.93	0.93	0.93	2.20	2.20	4.70	10.10	22.00
装药埋深he/m	3.2	3.2	3.2	3.2	3.2	3.2	3.2	5.3	5.3
比例爆距$\overline{R}$/（m·kg-1/3）	1.2	1.0	0.8	0.6	1.2	1.0	1.2	0.8	0.6

图2  试验过程图片

Fig.2  Process of the experiment

1.1.2　有限元计算模型及材料参数

采用Truegrid软件进行建模，LS-DYNA软件进行有限元计算。计算模型如图3所示。

图3  有限元模型

Fig.3  Finite element model

由于土中竖井结构受侧向爆炸问题具有对称性，为提高计算效率，建立1/2模型进行计算。土‒空气‒结构相互作用采用流固耦合算法。土体及空气模型采用欧拉网格，竖井井壁混凝土模型采用拉格朗日网格。钢板采用Shell单元，钢筋采用Beam单元，钢板与混凝土使用*CONSTRAINED_LAGRANGE_IN_SOLID命令耦合连接。为了避免边界应力波的反射，在主要计算模型外设置了带阻尼的土体，削弱传递至边界的应力波。装药使用无网格体积分数填充的方式设置。

有限元模型中混凝土*MAT_PSEUDO_ TENSOR模型为伪张量模型，该模型常用来描述混凝土在爆炸冲击荷载作用下的动力响应^［18］。该模型参数明确，可以通过输入单位转换系数及混凝土强度得到一个通用的混凝土强度模型，以最大强度面σ_max及残余强度面σ_failed的线性组合描述其破坏面，并自动根据泊松比和单轴抗压强度生成一个压力‒体积应变的状态方程。表3给出混凝土材料模型的参数设置。

表3  混凝土材料模型参数

Tab.3  Parameters of concrete material model

密度/（ kg·m-3）	泊松比	单轴抗压强度fc/MPa	单位转化参数a0
2 500	0.2	50	-0.000 145

钢板材料使用*MAT_SIMPLIFIED_ JOHNSON_COOK模型，该模型在J-C模型的基础上忽略热效应及热损伤，适用于钢板‒混凝土结构的动力响应计算。钢筋材料采用*MAT_ PLASTIC_KINEMATIC模型。其材料模型如表4、表5所示。

表4  钢筋材料模型参数^［19］

Tab.4  Parameters of steel material model

密度/（ kg·m-3）	泊松比	强度σy/108Pa	弹性模量E /1011Pa	塑性模量Et/GPa	失效参数
7 850	0.29	3.35	2.05	2	0.14

表5  钢板材料模型参数^［20］

Tab.5  Parameters of steel plate material model

密度/（ kg·m-3）	泊松比	应力函数系数A/MPa	应力函数系数B/MPa	应力函数系数N	应力函数系数C
7 800	0.3	345	336	0.42	0.026

试验中土体为含水率8%~12%的非饱和黏土，有限元中土体采用* MAT_SOIL_AND_FOAM 模型，模型参数可以参照文献［28］取值，如表6所示。

表6  土体材料模型参数^［21］

Tab.6  Parameters of soil material model

密度/（ kg·m-3）	剪切模量G /107Pa	体积卸载模量K/1010Pa	屈服函数系数a0/1012	屈服函数系数a1/106	屈服函数系数a2
1 970	6.38	3	1.39	3.64	0.445

试验采用TNT装药，有限元中装药采用*MAT_HIGH_ EXPLOSIVE_BURN模型，*EOS_JWL状态方程，空气建模采用*MAT_NULL模型，*EOS_ GRUNEISEN状态方程，参数按照文献［22-23］取值。

1.1.3　试验与有限元结果对比验证

图4为4号靶体和9号靶体试验与有限元计算结果的对比。可以发现，有限元计算结果得到的结构毁伤情况与试验结果较为吻合，4号靶体有限元计算结果与试验结果中井壁混凝土均出现环状贯穿裂缝，试验中出现裂缝的环向范围为±49.5°，有限元计算结果中裂缝环向范围为±53.1°，误差为6.7%。9号靶体有限元计算结果与试验结果中井壁混凝土均出现轴向较长的贯穿裂缝，试验中出现裂缝的环向范围为±51.5°，有限元计算结果中裂缝环向范围为±57.2°，误差为9.9%。

图4  试验与有限元计算毁伤效果对比

Fig.4  Comparison of experiment and finite element calculation

图5给出了各炮次靶体迎爆点（C1点）最大位移试验与有限元计算的对比。图中空缺的点为试验未能测到有效数据。

图5  各炮次C1点最大位移试验与有限元计算对比

Fig.5  Comparison of maximum displacement of C1 point

由图5可知，各个炮次得到的C1点最大位移试验及有限元计算结果的误差大部分均小于20%。说明有限元计算模型可以较好地反应钢板‒钢筋混凝土竖井结构在爆炸荷载作用下动态响应过程。

1.2　竖井结构动力响应特征分析

以6号炮次工况为对象研究竖井结构的动力响应及变形特征。图6给出了竖井结构在装药埋深平面内外径圆截面及轴向截面在各个时刻的变形情况，为了更好地表现变形特征，将变形值扩大为原始数值的10倍。

图6  6号炮次工况下不同时刻结构变形情况

Fig.6  Structural deformation of No.6 shot at different times

由图6可以发现，在8.0 ms及9.6 ms时，C1点出现的径向位移，明显大于C5点，且结构自由端还未出现明显变形，结构处于以迎爆点（C1）为中心的壳壁变形。而在20.0 ms时，C1点与C5点均出现较为明显的变形，圆截面变形程度减小，轴向变形中迎爆面与背爆面变形程度相近，结构变形以整体的梁式变形为主。

在6号炮次基础上将比例爆距改为0.6 m·kg^-1/3，研究竖井结构动力响应及变形特征。图7为$\overline{R}$ =0.6 m·kg^-1/3工况下截面的变形情况。

图7  $\overline{R}$ =0.6 m·kg^-1/3时不同时刻结构变形情况

Fig.7  Structural deformation at different times at $\overline{R}$ =0.6 m·kg^-1/3

由图7可以发现，$\overline{R}$ =0.6 m·kg^-1/3时，结构迎爆点（C1）的位移明显变大，结合图4a及图7可知此时结构已经出现了连续的截面贯穿裂缝，形成了塑性铰线。

图8给出了$\overline{R}$ =0.6 m·kg^-1/3工况下结构破坏过程。

图8  $\overline{R}$ =0.6 m·kg^-1/3时结构破坏过程

Fig.8  Failure process of structure at $\overline{R}$ =0.6 m·kg^-1/3

根据有限元分析，可知薄壁竖井结构变形过程中，结构上距C1点一定距离处出现沿轴向开展的贯穿裂缝，而后不断延伸为闭环的塑性铰线。因此，在结构弹塑性动力响应分析中更关注于结构的圆截面变形程度。通过图6、7发现，可以以C1点、C5点径向位移之和的最大值与结构外径的比值作为指标描述竖井结构横截面的变形情况。这个值可以较好地表征结构相对于初始圆形截面的变形，即结构环向变形。将这个量纲一环向相对位移记为${\alpha }_{\mathrm{1}}$，计算方法为

$${\alpha }_{\mathrm{1}}=\frac{W_{\mathrm{C}\mathrm{1}}+W_{\mathrm{C}\mathrm{5}}}{\mathrm{2}{r}_{\mathrm{1}}}$$

（1）

式中：$W_{\mathrm{C}\mathrm{1}}$、$W_{\mathrm{C}\mathrm{5}}$分别为C1点和C5点的径向位移。

在竖井结构表面出现闭环贯穿裂缝，形成塑性绞线之前，结构的破坏模式均处于带裂缝弹塑性动力响应阶段。根据有限元计算结果，发现当${\alpha }_{\mathrm{1}}\ge \mathrm{3}\mathrm{\%}$时结构会形成以C1为中心的闭环塑性铰线，结构的动力响应模式发生改变。

2 竖井结构弹性响应计算方法

2.1　结构参数等效及自振频率简化计算

竖井结构的自振频率是计算竖井结构在爆炸作用下的动力响应的必要参数^［24］，在计算钢板‒钢筋混凝土结构的自振频率时，通常将其等效为单一均质材料结构。

对于壳体结构，基于抗弯刚度和薄膜刚度得到钢板混凝土组合壳的等效参数。弹性协调变形阶段，钢板组合壳体结构的等效弹性模量、等效壳体厚度可以按照式（2）、（3）计算^［25］。钢‒混凝土组合结构的等效泊松比近似取0.2^［26］。

$$E_{\mathrm{e}}=\frac{\sqrt[]{\mathrm{3}}(\mathrm{1}-{{\nu }_{\mathrm{e}}}^{\mathrm{2}})K}{\mathrm{6}\sqrt[]{D/K}}$$

（2）

$$h_{\mathrm{e}}=\mathrm{2}\sqrt[]{\mathrm{3}}\sqrt[]{D/K}$$

（3）

$$r_{\mathrm{e}}=r_{\mathrm{2}}+h_{\mathrm{e}}/\mathrm{2}$$

（4）

式（2）~（4）中：$E_{\mathrm{e}}$为等效弹性模量；${\nu }_{\mathrm{e}}$为等效泊松比；$h_{\mathrm{e}}$为等效壳体厚度；$r_{\mathrm{e}}$为结构等效中面半径；$r_{\mathrm{2}}$为井壁内径；$D$为截面抗弯刚度；$K$为壳体薄膜刚度。

$$D=E{I}_{\mathrm{u}}$$

（5）

$$K=(\mathrm{1}-{\rho }_{\mathrm{p}})\frac{E_{\mathrm{c}}{h}_{\mathrm{c}}}{\mathrm{1}-{\nu }_{\mathrm{c}}^{\mathrm{2}}}+\frac{E_{\mathrm{s}}{h}_{\mathrm{s}}+h_{\mathrm{c}}{\rho }_{\mathrm{p}}}{\mathrm{1}-{\nu }_{\mathrm{s}}^{\mathrm{2}}}$$

（6）

式（5）、（6）中：$I_{\mathrm{u}}$为混凝土未开裂时的截面惯性矩，可按照钢板混凝土组合板的截面惯性矩进行计算^［27］ ；${\rho }_{\mathrm{p}}$为截面配筋率；$E_{\mathrm{c}}$为混凝土弹性模量；$h_{\mathrm{c}}$为混凝土厚度；${\nu }_{\mathrm{c}}$为混凝土泊松比；$E_{\mathrm{s}}$为钢弹性模量；$h_{\mathrm{s}}$为钢板厚度；${\nu }_{\mathrm{s}}$为钢泊松比。

在得到结构的等效参数后，可以结合LOVE壳体理论，按照变量分离法求解竖井结构的自振频率。李学斌^［28］对其进行了较为详细的推导，但其推导过程过于复杂，针对地下竖井结构为长细薄壁圆柱壳结构的特点，在其推导的基础上带入两个简化假设：①对于薄壳应该有厚径比远小于1，即$h_{\mathrm{e}}/r_{\mathrm{e}}$≪1；②对于较为细长的圆柱壳有周向波长远小于轴向波长，即$n^{\mathrm{2}}$≫${\lambda }^{\mathrm{2}}$，则有$n^{\mathrm{2}}\pm {\lambda }^{\mathrm{2}}\approx n^{\mathrm{2}}$。可以得到关于自振频率的方程如式（7）所示。

$$\begin{array}{l}{\lambda }_{mn}^{\mathrm{4}}=(\mathrm{2}{\Omega }_{mn}^{\mathrm{6}}-(\mathrm{2}+(\mathrm{3}-\nu )n^{\mathrm{2}}+\mathrm{2}k{n}^{\mathrm{4}}){\Omega }_{mn}^{\mathrm{4}}+\\ ((\mathrm{1}-\nu )(\mathrm{1}+n^{\mathrm{2}})n^{\mathrm{2}}+(\mathrm{3}-\nu )k{n}^{\mathrm{6}}){\Omega }_{mn}^{\mathrm{2}}-\\ (\mathrm{1}-\nu )k{n}^{\mathrm{8}})/(\mathrm{1}-\nu ){\nu }^{\mathrm{2}}\end{array}$$

（7）

式中：$k=h_{\mathrm{e}}^{\mathrm{2}}/\mathrm{12}{r}_{\mathrm{e}}^{\mathrm{2}}$为厚径比系数；${\Omega }_{mn}^{\mathrm{2}}$为频率系数；m， n分别为竖井结构振动的轴向和环向模态数；$\nu$为泊松比。

带入结构边界条件即可得到式（7）的解即为第（m， n）阶自振频率系数${\Omega }_{mn}$，分析式（7）可以发现${\Omega }^{\mathrm{2}}=\pm {{\Omega }_{\mathrm{1}}}^{\mathrm{2}}\mathrm{、}\pm {{\Omega }_{\mathrm{2}}}^{\mathrm{2}}\mathrm{、}\pm {{\Omega }_{\mathrm{3}}}^{\mathrm{2}}$，且有${\Omega }_{\mathrm{1}}^{\mathrm{2}}<{\Omega }_{\mathrm{2}}^{\mathrm{2}}<{\Omega }_{\mathrm{3}}^{\mathrm{2}}$，其中${\Omega }_{\mathrm{1}}$为弯曲振型的频率系数。

2.2　土中竖井结构弹性响应求解

根据Costantino提出的土与壳体结构相互作用的模型，结合化学爆炸中土体位移对结构上的荷载影响很小，可以将爆炸产生的壳体结构上径向荷载写为 ^［29］

$$P_{\mathrm{r}}=\mathrm{2}{p}_i^{}-{\rho }_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}c\dot{W}$$

（8）

式中：$P_{\mathrm{r}}$为作用在竖井结构上的荷载；$p_i$为结构界面处的自由场地冲击应力径向分量；${\rho }_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}$为土体介质的密度；$c$为土体压缩波速；$W$为竖井结构径向位移。

图9为土中竖井结构受侧向爆炸的示意图。图中r₁为竖井结构外径；$\kappa (\phi )$为结构某点至装药连线与结构迎爆点至装药点连线的夹角；R为装药点至结构迎爆点的距离；$L^{\mathrm{\text{'}}}(\phi ,x)$为装药点至结构上某点的距离。

图9  竖井受爆炸作用模型

Fig.9  Model of shaft subjected to explosion

爆炸产生的自由场地冲击应力如式（9）所示 ^［30］ 。

$$p_i=p_{i-\mathrm{0}}f(t)f(x,\phi )$$

（9）

式中：$p_{i-\mathrm{0}}$为自由场入射压力峰值。

忽略土中应力波竖井结构上的绕射作用，爆炸作用下土与结构界面上的径向应力分布如式（10）~（13）所示。式（10）中$\vartheta$为装药至结构某点连线方向与该点指向结构圆心方向的夹角，式（11）给出了$\mathrm{c}\mathrm{o}\mathrm{s}\vartheta$的计算方法。

$$f(x,\phi )={(\frac{L^{\mathrm{\text{'}}}(x,\phi )}{R})}^{-\mu }(\xi +(\mathrm{1}-\xi )\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\vartheta )$$

（10）

$$\mathrm{c}\mathrm{o}\mathrm{s}\vartheta =\mathrm{c}\mathrm{o}\mathrm{s}(\phi +\kappa (\phi ))\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}(\frac{\left|x-a_{\mathrm{0}}\right|}{R}))$$

（11）

$$\kappa (\phi )=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}(\frac{r_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}\phi }{R+r_{\mathrm{1}}-r_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}\phi })$$

（12）

$$L^{\mathrm{\text{'}}}(x,\phi )=\left\{\begin{array}{l}\sqrt[]{{(\frac{r_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}\phi }{\mathrm{s}\mathrm{i}\mathrm{n}\kappa })}^{\mathrm{2}}+(a_{\mathrm{0}}{-x)}^{\mathrm{2}}}\phi \ne \mathrm{0}\\ \sqrt[]{R^{\mathrm{2}}+(a_{\mathrm{0}}{-x)}^{\mathrm{2}}}\phi =\mathrm{0}\end{array}\right.$$

（13）

在工程计算中，认为以迎爆点（C1）为中心的某一范围内荷载是同时到达的，且忽略这一区域外的荷载。根据荷载分布特点，该区域环向取值为$\phi \in (-{\phi }_{\mathrm{0}},{\phi }_{\mathrm{0}})$，$x\in (\mathrm{m}\mathrm{i}\mathrm{n}((a_{\mathrm{0}}-R),\mathrm{0}),\mathrm{m}\mathrm{a}\mathrm{x}((a_{\mathrm{0}}+R),L))$。

根据振型分解法，及圆柱壳结构的振型正交条件^［31］可以写出结构振动微分方程的解为

$${\ddot{T}}_{mn}(t)+\mathrm{2}{\xi }_{mn}{\omega }_{mn}{\dot{T}}_{mn}(t)+{\omega }_{mn}^{\mathrm{2}}{T}_{mn}(t)=\frac{P_{mn}{p}_i(t)}{M_{mn}}$$

（14）

$${\xi }_{mn}=\frac{{\rho }_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}{c}_{\mathrm{s}}}{\mathrm{2}(\mathrm{1}+{\alpha }_{mn}^{\mathrm{2}}+{\beta }_{mn}^{\mathrm{2}})({\rho }_{\mathrm{s}}{h}_{\mathrm{s}}+{\rho }_{\mathrm{c}}{h}_{\mathrm{c}}){\omega }_{mn}}$$

（15）

$$P_{mn}={\int }_{\mathrm{0}}^L{\int }_{\mathrm{0}}^{\mathrm{2}\pi }(U_{mn}{q}_x+V_{mn}{q}_{\phi }+W_{mn}{q}_r)\mathrm{d}\theta \mathrm{d}x$$

（16）

$$M_{mn}={\int }_{\mathrm{0}}^L{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}\rho h(U_{mn}^{\mathrm{2}}+V_{mn}^{\mathrm{2}}+W_{mn}^{\mathrm{2}})\mathrm{d}\theta \mathrm{d}x$$

（17）

式（14）~（17）中：$P_{nm}$为广义力；$M_{nm}$为广义质量；${\omega }_{mn}$为第（m， n）阶自振频率；$T_{mn}(t)$为第（m， n）阶振型的振动函数；${\xi }_{mn}$为第（m， n）阶振型的土‒结构辐射阻尼系数；${\rho }_{\mathrm{s}},{\rho }_{\mathrm{c}}$分别为钢板、混凝土的密度；${\alpha }_{mn}^{},{\beta }_{mn}^{}$分别为结构环向振型、轴向振型与径向振型的振幅比；$W_{mn},V_{mn},U_{mn}$分别为第（m， n）阶结构径向、环向、轴向振型函数；$q_x,q_{\phi },q_r$分别为结构上轴向、环向、径向荷载分布函数。

式（14）的解可分段写出。

在荷载升压时间t_r前，即0<t<t_r时有

$$\begin{array}{l}{T}_{mn}(t)=\frac{P_{mn}{\int }_{\mathrm{0}}^t(t/t_{\mathrm{r}}){\mathrm{e}}^{-\xi {\omega }_{mn}^{}(t-\tau )}\mathrm{s}\mathrm{i}\mathrm{n}({\omega }_{\mathrm{d},mn}(t-\tau ))\mathrm{d}\tau }{{\omega }_{\mathrm{d},mn}{M}_{mn}},\\ {\omega }_{\mathrm{d},mn}=\sqrt[]{\mathrm{1}-{\xi }^{\mathrm{2}}}{\omega }_{mn}^{}\end{array}$$

（18）

t_r <t时有

$$T_{mn}(t)=\frac{P_{mn}{\int }_{t_{\mathrm{r}}}^t\begin{array}{l}{e}^{-\frac{(\tau -t_{\mathrm{r}})}{t_{\mathrm{a}}}}{e}^{-\xi {\omega }_{mn}^{}(t-\tau )}\\ \mathrm{s}\mathrm{i}\mathrm{n}({\omega }_{\mathrm{d},mn}\tau (t-\tau ))\end{array}\mathrm{d}\tau +{\varphi }_{mn}(t)}{{\omega }_{\mathrm{d}}{M}_{mn}}$$

（19）

$${\varphi }_{mn}(t)=e^{-\xi {\omega }_{mn}^{}(t-t_{\mathrm{r}})}\left(\begin{array}{l}{T}_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}({\omega }_{\mathrm{d},mn}(t-tr)+\\ \frac{{\dot{T}}_{\mathrm{0}}+\xi {\omega }_{mn}^{}{T}_{\mathrm{0}}}{{\omega }_{\mathrm{d},mn}}\mathrm{s}\mathrm{i}\mathrm{n}({\omega }_{\mathrm{d},mn}(t-t_{\mathrm{r}}^{})\end{array}\right)$$

（20）

式（20）中： $t_{\mathrm{a}}$为土中应力波达到结构表面的时间；${\omega }_{\mathrm{d},mn}$为阻尼修正的结构自振频率；$T_{\mathrm{0}}$为升压阶段结束时第（m， n）阶振型的结构位移。

2.3　计算方法的验证

以1号靶体为例，将比例爆距改为$\overline{R}$ =2.5 m·kg^‒1/3进行计算。有限元计算和理论计算的对比如图10所示。

图10  有限元计算与理论计算结果对比

Fig.10  Comparison of finite element calculation and theoretical calculation

由图10可以发现，理论与有限元计算结果的吻合程度较高，误差均在20%以内，说明上述理论计算方法具有一定的可行性。

3 基于弹性理论计算的弹塑性变形

3.1　弹性计算变形与实际弹塑性变形

在竖井结构表面出现闭环贯穿裂缝，形成塑性铰线之前，结构的破坏模式均处于带裂缝弹塑性动力响应阶段。该阶段结构产生的裂缝引起了结构的刚度下降和结构残余位移，实际的结构弹塑性变形比按照弹性理论计算的大很多。为了分析有限元计算得到的实际弹塑性变形值${\alpha }_{\mathrm{1}\mathrm{p}}$与弹性理论计算得到的${\alpha }_{\mathrm{1}\mathrm{e}}$的相关关系，以5号炮次工况为例，计算了比例爆距由2.5 m·kg^-1/3至0.7 m·kg^-1/3之间各个工况的变形情况，图11给出了${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$与${\alpha }_{\mathrm{1}\mathrm{e}}$的关系。

图11  ${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$与${\alpha }_{\mathrm{1}\mathrm{e}}$的关系

Fig.11  ${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$ versus ${\alpha }_{\mathrm{1}\mathrm{e}}$

由图11可以发现，${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$随着${\alpha }_{\mathrm{1}\mathrm{e}}$的增大而增大，这是由于${\alpha }_{\mathrm{1}\mathrm{e}}$虽然是按照弹性解得到的量纲一变形值，但也能说明结构的变形程度，${\alpha }_{\mathrm{1}\mathrm{e}}$越大，结构损伤程度越高，刚度下降程度越大，实际的塑性变形值与按照弹性计算的变形值相差越大。

根据图11及其他有限元计算结果，当结构处于弹性阶段时，实际弹塑性变形值与弹性理论计算值之比${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$应为1，对于弹塑性变形而言，弹性阶段占比很小，可以近似的认为当${\alpha }_{\mathrm{1}\mathrm{e}}$等于1时，${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$等于1，则可以按照式（21）对其进行拟合。

$${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}=\kappa {\alpha }_{\mathrm{1}\mathrm{e}}^{\mathrm{2}}+\mathrm{1}$$

（21）

式中：$\kappa$为拟合系数。

3.2　拟合系数的参数分析及计算方法

土中竖井结构受侧向爆炸作用受诸多因素影响，主要可以分为以下4类：①结构尺寸参数；②结构材料参数；③装药参数；④土体介质参数。

认为钢板‒钢筋混凝土竖井结构配筋率不变，主要研究钢板与混凝土对结构响应的影响，不考虑钢筋参数的影响。

钢板‒钢筋混凝土竖井结构的尺寸主要有结构长度L，内径r₂，钢板及混凝土厚度h_s、h_c。对于结构材料参数，在弹塑性响应阶段，钢板、钢筋一般不会屈服，且钢材的弹性模量变化范围较小。因此，材料参数中仅考虑混凝土强度f_c的影响。对于装药参数，主要考虑装药TNT当量Q，以及装药埋深与竖井底部埋深的差值a₀。同时，现阶段主要研究单一土体介质，不考虑土体介质参数影响。

3.2.1　装药量的影响

图12给出了h_s/h_c=0.017 1，h_e/r_e=0.194，f_c为50 MPa，L/r_e=7.95，a₀/L=0.536，不同装药量情况下拟合系数$\kappa$的变化。0.93、2.20、4.70 kg装药对应了表1中靶体1，4.30、10.10、22.00 kg装药对应了靶体2。

图12  装药量对拟合系数$\kappa$的影响

Fig.12  Effect of charge on fitting coefficient $\kappa$

由图12可以发现，装药量大小及结构的缩尺比例对拟合系数$\kappa$的影响较小，误差小于±5%，因此可以忽略装药量大小的影响。

3.2.2　装药与结构特征尺寸的影响

对h_e/r_e、L/r_e、h_s/h_c、a₀/L 4个特征尺寸参数进行分析，计算工况中混凝土强度为50 MPa，装药量为2.20 kg。分别研究4个参数变化对拟合参数的影响。

首先观察L/r_e对参数$\kappa$的影响，如图13所示。

图13  L/r_e对拟合系数$\kappa$的影响

Fig.13  Effect of L/r_e on fitting coefficient $\kappa$

由图 13可以发现，L/r_e大小对拟合系数$\kappa$的影响较小，误差小于±10%，因此可以忽略L/r_e的影响。

由式（2）~（6）可以发现，r₂、h_s、h_c均会影响等效厚径比h_e/r_e。因此，可以对不同h_s/h_c情况下h_e/r_e的变化对拟合系数$\kappa$的影响进行研究，如图14a所示。a₀/L对拟合系数$\kappa$的影响如图14b所示。

图14  $\kappa$的影响因素及规律

Fig.14  Influencing factors and rules of$\kappa$

可以发现，随着h_e/r_e、h_s/h_c的增大，拟合系数$\kappa$均增大，即${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$增大。这是由于随着h_e/r_e、h_s/h_c的增大，结构的抗力增大，变形能力增强，理论计算时未考虑结构的塑性变形，计算得到的${\alpha }_{\mathrm{1}\mathrm{e}}$较小，而结构塑性变形${\alpha }_{\mathrm{1}\mathrm{p}}$变大，${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$也增大。在a₀/L较小时，系数$\kappa$较大，即${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$较大。这是因为当a₀/L较小时，由于固定端的影响，理论计算得到的${\alpha }_{\mathrm{1}\mathrm{e}}$较小，但实际变形时，由于固定端截面被破坏，实际变形${\alpha }_{\mathrm{1}\mathrm{p}}$较大，${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$也较大。这种端部的影响在a₀/L较大时不再存在，因此当a₀/L大于0.6后，$\kappa$值几乎不再变化。

3.2.3　混凝土强度的影响

由式（2）~（6）可以发现，混凝土强度会通过影响混凝土的弹性模量影响等效厚径比h_e/r_e。因此，应当对不同混凝土强度下h_e/r_e的变化对拟合系数$\kappa$的影响进行研究。计算工况中L/r_e=7.95，h_s/h_c=0.03，a₀/L=0.536。图15给出了不同混凝土强度下h_e/r_e的变化对拟合系数$\kappa$影响。

图15  混凝土强度对拟合系数$\kappa$的影响

Fig.15  Effect of concrete strength on fitting coefficient $\kappa$

由图15可以发现，混凝土强度的增大使$\kappa$减小，即${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$减小。这是由于混凝土强度对弹性理论计算得到的${\alpha }_{\mathrm{1}\mathrm{e}}$影响较小，但混凝土强度较大时实际弹塑性变形会明显减小，从而${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$减小。

3.2.4　拟合系数$\kappa$的确定方法

由上述分析可以发现，拟合系数$\kappa$主要受等效厚径比h_e/r_e的影响。此外，装药埋深与结构的相对高度a₀/L也会影响拟合系数$\kappa$。

根据图11、12、15中的数据拟合，可以将${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$的计算工程算法写为

$${\alpha }_{\mathrm{1}}^{\mathrm{\text{'}}}={\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}=\kappa {\alpha }_{\mathrm{1}\mathrm{e}}^{\mathrm{2}}+\mathrm{1}$$

（22）

$$\kappa =\mathrm{6.14}\times {\mathrm{10}}^{\mathrm{8}}{\beta }_{a_{\mathrm{0}}/L}{(h_{\mathrm{e}}/r_{\mathrm{e}})}^{A_{\mathrm{2}}}$$

（23）

$$A_{\mathrm{2}}={\beta }_{f_{\mathrm{c}}}(\mathrm{5.49}-\mathrm{1.68}(h_{\mathrm{s}}/h_{\mathrm{c}}))$$

（24）

式（23）、（24）中：${\beta }_{a_{\mathrm{0}}/L}$为端部影响系数；${\beta }_{f_{\mathrm{c}}}$为混凝土强度影响系数。${\beta }_{a_{\mathrm{0}}/L}$取值方法如下所示。${\beta }_{f_{\mathrm{c}}}$取值如表7所示。

$${\beta }_{a_{\mathrm{0}}/L}=\mathrm{0.99}+\mathrm{201}{\mathrm{e}}^{-\mathrm{12.45}{a}_{\mathrm{0}}/L}$$

（25）

表7  不同混凝土强度下${\beta }_{f_c}$取值表

Tab.7  Value of ${\beta }_{f_c}$ for different concrete strengths

混凝土强度/MPa	30	40	50	60	70	80
${\beta }_{f_{\mathrm{c}}}$取值	0.978	0.989	1.000	1.010	1.017	1.023

4 结论

本文开展了土中钢板‒钢筋混凝土竖井结构受侧向爆炸的野外模型试验，通过试验与有限元计算结果的验证，得到了可行的有限元计算模型。根据有限元计算结果，讨论了竖井结构的弹塑性变形特征，提出了结构变形的表征量。推导了常规武器爆炸作用下竖井结构的弹性动力响应计算方法，并给出了弹塑性变形值的计算方法。得到了以下结论：

（1）开展了9组外场模型试验，对试验及有限元计算结果进行了对比验证。有限元计算结果得到的中远距离爆炸作用下结构上荷载、结构位移、裂缝分布与试验吻合度较高。验证了有限元计算模型的有效性。

（2）常规武器爆炸作用下竖井结构首先出现壳壁变形，结构圆截面出现明显变形，而后出现悬臂梁式的整体变形。当装药距离较近，埋深距结构固定端截面较远时，结构首先出现由于圆截面变形引起的沿轴向伸展的裂缝，而后扩展为闭环塑性铰线。根据有限元计算竖井结构破坏现象，提出以量纲一环向相对位移${\alpha }_{\mathrm{1}}$作为结构变形程度的判据，确定了${\alpha }_{\mathrm{1}}=\mathrm{3}\mathrm{\%}$为闭环塑性铰线形成的阈值。

（3）基于LOVE壳体理论推导了一端固支、一端自由的竖井结构的自振频率，并针对地下防护工程中常见的长细薄壁竖井进行了自振频率计算方法的简化。引入Costantino模型，推导了考虑土与结构相互作用的爆炸荷载下竖井结构的弹性动力响应计算方法，并通过有限元计算结果的对比验证了理论计算方法的可行性。

（4）分析了按照弹性理论计算的量纲一环向相对位移${\alpha }_{\mathrm{1}\mathrm{e}}$和有限元计算得到的弹塑性量纲一环向相对位移${\alpha }_{\mathrm{1}\mathrm{p}}$的关系，并进行了参数分析。发现结构长径比L/r_e、结构缩尺比例、装药量大小对${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$随${\alpha }_{\mathrm{1}\mathrm{e}}$的变化规律几乎没有影响。等效厚径比h_e/r_e、钢板与混凝土厚度比h_s/h_c较大及结构的相对高度a₀/L较小时理论计算会高估结构抗力，计算得到的${\alpha }_{\mathrm{1}\mathrm{e}}$偏小，${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$较大；混凝土强度较大时，相同${\alpha }_{\mathrm{1}\mathrm{e}}$情况下结构塑性变形较小，${\alpha }_{\mathrm{1}\mathrm{p}}/{\alpha }_{\mathrm{1}\mathrm{e}}$较小。

（5）考虑了等效厚径比h_e/r_e、钢板与混凝土厚度比h_s/h_c、装药埋深与结构的相对高度a₀/L、混凝土强度f_c的影响，建立了基于弹性计算${\alpha }_{\mathrm{1}\mathrm{e}}$的弹塑性量纲一环向相对位移${\alpha }_{\mathrm{1}\mathrm{p}}$的工程算法。

作者贡献声明

孙善政：研究构思，理论推导，数值模拟计算。

卢 浩：稿件撰写与审核。

熊自明：数值模拟计算结果对比验证。

马云飞：算法完善，稿件修订。