空爆载荷作用下固支弹塑性圆板的动力学模型

吴迪; 米国; 郭香华; 张庆明

doi:10.11858/gywlxb.20220525

空爆载荷作用下固支弹塑性圆板的动力学模型

doi: 10.11858/gywlxb.20220525

北京理工大学机电学院, 北京　100081

详细信息

作者简介:
吴　迪（1997－），女，硕士研究生，主要从事结构冲击动力学研究. E-mail：15010059608@163.com

通讯作者:
郭香华（1974－），男，博士，副教授，主要从事爆炸与冲击仿真、材料与结构动力学研究.E-mail：guoxh@bit.edu.cn

中图分类号: O383
计量
- 文章访问数: 263
- HTML全文浏览量: 179
- PDF下载量: 44
出版历程
- 收稿日期: 2022-03-07
- 修回日期: 2022-03-21
- 刊出日期: 2022-10-11

Dynamic Model of Clamped Elastoplastic Circular Plate under Air Blast Loading

School of Mechatronical Engineering, Beijing Institute of Technology, Beijing 100081, China

摘要

摘要: 基于等效单自由度模型理论，考虑变形过程中弯矩和膜力的联合承载，建立了空中近爆作用下固支弹塑性圆板的动力学模型，实现了圆板加载和卸载的全过程描述。基于文献中的试验工况建立有限元模型，利用有限元数值模拟对空爆作用下固支圆板的动态响应进行分析，并将动力学模型计算结果与试验结果及有限元数值模拟结果进行对比，验证了动力学模型计算结果的准确性。结果表明，理论计算结果与试验结果及数值模拟结果吻合较好，建立的动力学模型可以应用于固支圆板在空中近爆载荷作用下的大变形预测，为结构抗爆提供技术支持。
- 单自由度模型 /
- 弹塑性响应 /
- 空爆载荷 /
- 膜力 /
- LS-DYNA
Abstract: Based on equivalent single degree-of-freedom (SDOF) theory, a dynamic model of clamped elastoplastic circular plate under near air blast loading is established by considering the bending moment and the membrane force during the deformation process, and the whole process containing loading and unloading for circular plate is described. The finite element model is established using test parameters in the literature, and the dynamic response of clamped circular plate under air blast loading is analyzed by finite element simulation. After comparing the calculation results of the dynamic model and the finite element simulation results, the accuracy of the calculation results of the dynamic model is verified. The results show that the theoretical calculation results are in good agreement with the test results and the finite element simulation results. The proposed dynamic model can be applied to predict the large deformation of clamped circular plate under near air blast loading, and provides technical support for blast resistant structure.
- single degree-of-freedom model /
- elastoplastic response /
- air blast loading /
- membraneforce /
- LS-DYNA

HTML全文

在结构的动力响应中，通常既有弹性变形，又有塑性变形，并且两种变形场的分界面随时间不断发生变化，因此，在求解弹塑性动力响应时，不仅需要对不同区域采用不同的本构关系，还要处理复杂的动边界问题。现有的关于结构在强动载荷下响应的理论研究^[1-3]大多把材料假定为理想刚塑性，同时忽略材料的应变强化效应和应变率效应，从而简化理论分析过程。Duffey^[4]通过构建板的初始动能等于塑性变形耗散能的能量守恒方程来计算固支圆板在爆炸载荷下的最终挠度，这种方法假定板的位移分布始终满足形函数，并且径向位移的影响可以忽略。然而这一假设仅在小挠度情况下是合理的，当板的中心点挠度与板的厚度为同量级甚至存在几十倍差异时，膜力效应将逐渐显著，甚至超过弯矩效应而在抵抗板的变形中起主导作用。

Nurick等^[5]指出，忽略径向位移对估算平板挠度的影响不大，但会导致计算的膜应变分布与试验获得的应变分布相反，即中心小、外围大，为此提出了一个考虑径向位移的理论模型，得到了与试验结果吻合的膜应变分布。Jones^[6]提出了一种同时考虑弯矩和膜力的理论方法，用以描述在弯矩和膜力都很重要的挠度范围内刚塑性圆板的力学行为，比较准确地预测了Florence^[7]记录的简支圆板在均匀载荷作用下的最终变形。然而这种近似理论分析虽然同时考虑了板的膜效应和弯曲效应，但没有考虑它们的联合效应，即弯曲效应和膜效应是解耦的。Cloete等^[8]考虑径向位移并将弯曲应变和膜应变耦合处理，提出了刚塑性薄圆板在冲击载荷下的解析解，但没有考虑结构的弹性响应。

近年来，研究人员开始将研究重点转移到薄板大挠度变形的弹性效应上。Campbell等^[9]研究了集中力作用下固支梁响应的理论模型，将梁的响应分为3个阶段，即弹性阶段、弹塑性阶段和塑性阶段，并基于力矩平衡构建了梁在不同阶段的载荷-位移关系式，但模型中没有考虑卸载阶段。Karagiozova等^[10]研究了固支圆板在准静态均布载荷下的弹塑性响应，基于能量法得到了圆板在加载和卸载阶段的载荷-位移关系式。

等效单自由度（single degree-of-freedom，SDOF）模型简单实用，输入参数相对较少，在工程中常用于预测结构在爆炸载荷作用下的动力响应。SDOF模型由Biggs^[11]提出，其基本物理量是等效质量、等效刚度和等效抗力等。

战场上装甲车辆防护的根本目标是避免乘员受到致命伤害，因此车辆结构在爆炸载荷作用下的变形不能过大，否则舱内超压及结构碰撞形成的破片都将危及乘员生命。研究钢板在爆炸载荷下的弹塑性响应对于装甲车辆防雷结构设计具有重要意义。本研究针对固支弹塑性圆板在空中近爆载荷作用下的变形，建立其加载和卸载全过程的动力学模型，研究弯矩和膜力联合承载范围内的挠度。钢板承受的载荷并非均匀分布，爆轰产物首先作用于固支圆板的中心，随后扩散并作用于整个圆板表面，因此，圆板中心挠度最大且最先进入塑性阶段，并迅速向板边界扩散。本研究的工况中圆板最终仅发生弹塑性大变形，而未发生破坏。基于文献[12]中的试验工况建立有限元模型，通过有限元数值模拟得到空爆作用下固支圆板的动态响应，对比动力学模型计算结果、试验结果及数值模拟结果，验证动力学模型的有效性，研究结果可以应用于固支圆板在空中近爆载荷作用下的大变形预测。

1. 动力学模型推导

1.1 载荷等效

空中近爆作用下，作用在板上的爆炸冲击载荷具有明显的时空分布。可采用能量等效的方法将不均匀载荷转化为等效均匀载荷^[13]。近爆作用下结构上的爆炸冲击载荷达到峰值的时间差别很小，因此，可以假设载荷只有空间分布，即圆板上各点载荷随时间变化的函数一致，如图1所示。

图 1 载荷等效示意图

Figure 1. Equivalent load diagram

下载: 全尺寸图片幻灯片

将载荷分布函数 $p\left(r,t\right)$ 转化为能量等效均匀载荷

${\int }_{0}^{2\text{π}}{\int }_{0}^{R}p(r,t)\varphi \left(r\right)r\mathrm{d}r\mathrm{d}\theta ={p}_{\mathrm{e}\mathrm{e}}(t){\int }_{0}^{2\text{π}}{\int }_{0}^{R}\varphi (r)r\mathrm{d}r\mathrm{d}\theta$

(1)

${p}_{\mathrm{e}\mathrm{e}}(t)=\frac{{\displaystyle\int }_{0}^{2\text{π}}{\displaystyle\int }_{0}^{R}p(r,t)\varphi (r)r\mathrm{d}rd\theta }{{\displaystyle\int }_{0}^{2\text{π}}{\displaystyle\int }_{0}^{R}\varphi (r)r\mathrm{d}r\mathrm{d}\theta }$

(2)

式中： ${p}_{\mathrm{e}\mathrm{e}}(t)$ 为能量等效均匀载荷， $\varphi \left(r\right)$ 为形函数，R为圆板半径，r为板元距圆心距离，t为时间，θ为角度。

假设加载阶段圆板挠度沿半径的分布满足^{[4, 8, 14]}

$w(r,t)=W(t)\left[1-{\left(\frac{r}{R}\right)}^{2}\right]$

(3)

式中： $W(t)$ 为圆板中心点挠度。

圆板的初始面积 ${A}_{1}=\text{π}{R}^{2}$ 。在加载过程中，由于横向位移增加，假定固支圆板均匀变薄，由形函数计算得到曲面积分

${A}_{2}=\frac{\text{π}}{6{R}^{2}{W}^{2}}\left[{\left({R}^{4}+4{R}^{2}{W}^{2}\right)}^{\tfrac{3}{2}}-{R}^{6}\mathrm{c}\mathrm{s}\mathrm{g}\mathrm{n}({R}^{2})\right]$

(4)

由于积分结果较复杂，可将其替换为

${A}'_{2}=\text{π} \left({R}^{2}+{W}^{2}\right)$

(5)

图2给出了式(4)和式(5)的计算结果对比。可见，当中心点挠度为0.15 m时，通过式(4)和式(5)计算出的曲面积分值的相对误差仅为0.44%，可以进行替换。因此，变形板的面积与初始板面积之比约为 $({R}^{2}+{W}^{2})/{R}^{2}$ ，板的厚度h可以表示为

图 2 式

$\left(4\right)$ 和式

$\left(5\right)$ 的计算结果对比

Figure 2. Comparison of calculation results between Eq. (4) and Eq. (5)

下载: 全尺寸图片幻灯片

$h=\frac{H{R}^{2}}{{R}^{2}+{W}^{2}}$

(6)

式中：H为圆板的初始厚度。

1.2 屈服条件

当板挠度的量级超过板厚度时，膜力和弯矩的作用同样重要，板元的受力分析如图3所示，其中： ${N}_{r}$ 为径向膜力， ${N}_{\theta }$ 为周向膜力， ${M}_{r}$ 为径向弯矩， ${M}_{\theta }$ 为周向弯矩。

图 3 板元受力分析

Figure 3. Force analysis of plate microelement

下载: 全尺寸图片幻灯片

因此，如图4所示，采用弯矩和膜力联合承载的屈服条件^[15]

图 4 弯矩和膜力联合承载的屈服条件

Figure 4. Yield condition of combined bending moment and membrane force

下载: 全尺寸图片幻灯片

$\left|\frac{{M}_{\theta }}{{m}_{0}}\right|+{\left(\frac{{N}_{\theta }}{{n}_{0}}\right)}^{2}=1$

(7)

式中： ${m}_{0}={\sigma }_{0}{h}^{2}/4$ 和 ${n}_{0}={\sigma }_{0}h$ 分别为单位长度的塑性极限弯矩和极限膜力，其中 ${\sigma }_{0}$ 为材料的屈服应力。

在每个铰圆上有

$\left|\frac{{M}_{r}}{{m}_{0}}\right|+{\left(\frac{{N}_{r}}{{n}_{0}}\right)}^{2}=1$

(8)

由屈服条件可知，当膜力等于 ${n}_{0}$ 时，弯矩等于零，此时该区域开始进入塑性膜阶段。

1.3 加载阶段的应力场和变形能场

固支圆板发生大变形过程中，除了形成较大的横向挠度外，还会发生径向和周向位移，其应变分布为

${\varepsilon }_{r}=\frac{1}{2}{\left(\frac{\mathrm{d}W}{\mathrm{d}r}\right)}^{2}+\frac{\mathrm{d}u}{\mathrm{d}r},\quad{\varepsilon }_{\theta }=\frac{u}{r}$

(9)

式中： $u$ 为径向位移， ${\varepsilon }_{r}$ 和 ${\varepsilon }_{\theta }$ 分别为径向和周向膜应变。

Taylor^[14]证明了均匀静载下理想塑性膜近似满足 ${{\varepsilon }_{r}=\varepsilon }_{\theta }$ ，即径向和周向膜应力相等。Karagiozova等^[10]在研究均匀准静态压力作用下理想弹塑性圆板的响应时也假定了径向和周向应变相等，并通过分析证明了基于该假设可以获得变形的合理预测，简化了计算过程。本研究也采用此方法对变形几何方程进行简化。

将式 $\left(3\right)$ 代入式 $\left(9\right)$ ，令径向和周向应变相等，并结合边界条件，得到

$u=\frac{{W}^{2}}{R}\left[\frac{r}{R}-{\left(\frac{r}{R}\right)}^{3}\right]$

(10)

于是径向和周向应变分布为

${\varepsilon }_{r}={\varepsilon }_{\theta }={\left(\frac{W}{R}\right)}^{2}\left(1-\frac{{r}^{2}}{{R}^{2}}\right)$

(11)

径向和周向曲率为

${\kappa }_{r}=-\frac{{\mathrm{d}}^{2}w}{\mathrm{d}{r}^{2}}=\frac{2W}{{R}^{2}},\quad{\kappa }_{\theta }=-\frac{1}{r}\frac{\mathrm{d}w}{\mathrm{d}r}=\frac{2W}{{R}^{2}}$

(12)

因此，对于膜应变和弯曲应变都不能忽略的情况，径向和周向应变分布为

${\varepsilon }_{r}={\varepsilon }_{\theta }={\left(\frac{W}{R}\right)}^{2}\left(1-\frac{{r}^{2}}{{R}^{2}}\right)+2\frac{W\textit{z}}{{R}^{2}}$

(13)

通过板厚度对应变进行积分，在弹性加载阶段，膜力和弯矩的表达式如下

${N}_{r}={N}_{\theta }=\frac{E}{1-\mu}{\int }_{-\tfrac{h}{2}}^{\tfrac{h}{2}}{\varepsilon }_{r}\mathrm{d}\textit{z}=\frac{Eh}{1-\mu}{\left(\frac{W}{R}\right)}^{2}\left(1-\frac{{r}^{2}}{{R}^{2}}\right)$

(14)

${M}_{r}={M}_{\theta }=\frac{E}{1-\mu}{\int }_{-\tfrac{h}{2}}^{\tfrac{h}{2}}{\varepsilon }_{r}\textit{z}\mathrm{d}\textit{z}=\frac{EW{h}^{3}}{6{R}^{2}\left(1-\mu\right)}$

(15)

式中：E为材料的弹性模量，μ为泊松比，z为圆板的横向坐标。

在加载过程中，板的变形能变化为

$\delta {U}_{\mathrm{E}}={\int }_{0}^{2\text{π}}{\int }_{0}^{R}\left({M}_{r}\mathrm{\delta }{\kappa }_{r}+{M}_{\theta }\mathrm{\delta }{\kappa }_{\theta }+{N}_{r}\mathrm{\delta }{\varepsilon }_{r}+{N}_{\theta }\mathrm{\delta }{\varepsilon }_{\theta }\right)r\mathrm{d}r\mathrm{d}\theta =2\text{π}{\int }_{0}^{R}\left({2M}_{r}\mathrm{\delta }{\kappa }_{r}+{2N}_{r}\mathrm{\delta }{\varepsilon }_{r}\right)r\mathrm{d}r$

(16)

当固支圆板发生足够大的横向位移时，根据弯矩与膜力之间的不同关系，沿板的半径可以定义3个区域：弹性区、塑性膜区、弯矩与膜力联合承载的塑性区，如图5所示。最靠近边界的区域，即 ${R}_{\mathrm{P}}$ <r≤R区域，为弹性区。弹性区以外的区域都是塑性区，满足屈服条件，即式 $\left(7\right)$ 和式 $\left(8\right)$ ，其中：0≤r≤ ${R}_{\mathrm{P}\mathrm{N}}$ 区域满足N = n₀，M = 0，结构进入塑性流动，此时仅膜力参与承载，也就是塑性膜区域；而 ${R}_{\mathrm{P}\mathrm{N}}$ < $r$ ≤ ${R}_{\mathrm{P}}$ 区域则是弯矩与膜力联合参与承载的塑性区。

图 5 固支圆板发生弹塑性大变形后的应力场分布

Figure 5. Stress field of clamped circular plate after large elastic-plastic deformation

下载: 全尺寸图片幻灯片

半径 ${R}_{\mathrm{P}}$ 沿圆板半径r分隔弹性和塑性区域，将式 $\left(14\right)$ 和式 $\left(15\right)$ 代入式 $\left(7\right)$ 可以求得

${R}_{\mathrm{P}}=R\sqrt{1-\sqrt{\frac{{\left(1-\mu\right)}^{2}{\sigma }_{0}^{2}}{{E}^{2}}{\left(\frac{R}{W}\right)}^{4}\left(1-\frac{2E}{3\left(1-\mu\right){\sigma }_{0}}\frac{Wh}{{R}^{2}}\right)}}$

(17)

将 ${R}_{\mathrm{P}}=0$ 代入式 $\left(17\right)$ ，可得中心点屈服时的挠度 ${W}_{\mathrm{P}\mathrm{M}}$ ；将 ${R}_{\mathrm{P}}=R$ 代入式 $\left(17\right)$ ，可得边界屈服时的中心点挠度 ${W}_{\mathrm{P}\mathrm{B}}$ 。

板的塑性区从板中心向边界传播。半径 ${R}_{\mathrm{P}\mathrm{N}}$ 表示膜力达到极限值（ $N({R}_{\mathrm{P}\mathrm{N}})={n}_{0}$ ）的板半径，通过极限膜力可以求得

${R}_{\mathrm{P}\mathrm{N}}=R\sqrt{1-\frac{{R}^{2}}{{W}^{2}}\frac{\left(1-\mu\right){\sigma }_{0}}{E}}$

(18)

将 ${R}_{\mathrm{P}\mathrm{N}}=0$ 代入式 $\left(18\right)$ ，可得中心点进入塑性膜状态时的挠度 ${W}_{\mathrm{P}\mathrm{N}\mathrm{M}}$ 。

取圆板中面圆心为基点建立质量-弹簧-阻尼模型，如图6所示，其中： $p\left(t\right)$ 为等效载荷， $m$ 为等效质量， ${K}_{\mathrm{E}{\mathrm{M}}_{1}}$ 和 ${K}_{\mathrm{E}{\mathrm{N}}_{1}}$ 为等效刚度， ${f}_{\mathrm{P}{\mathrm{M}}_{1}}$ 、 ${f}_{\mathrm{P}{\mathrm{N}}_{1}}$ 、 ${f}_{\mathrm{P}{\mathrm{N}}_{2}}$ 和 ${f}_{\mathrm{B}}$ 为等效抗力。根据质量-弹簧-阻尼模型的变形能与圆板应变能相等，可以求解动力学模型的等效刚度和等效抗力。

图 6 加载阶段的等效质量-弹簧-阻尼模型

Figure 6. Equivalent mass-spring-damper model in loading phase

下载: 全尺寸图片幻灯片

1.3.1 弹性区加载

弹性区 ${R}_{\mathrm{P}}\left(W\right)$ <r≤ $R\left(W\right)$ 内，满足

$\delta {\varepsilon }_{r}=\delta {\varepsilon }_{\theta }=2\delta W\frac{W}{{R}^{2}}\left(1-\frac{{r}^{2}}{{R}^{2}}\right)$

(19)

${\mathrm{\delta }\kappa }_{r}=\delta {\kappa }_{\mathrm{\theta }}=\frac{2\mathrm{\delta }W}{{R}^{2}}$

(20)

${N}_{r}={N}_{\theta }=\frac{Eh}{1-\mu}{\left(\frac{W}{R}\right)}^{2}\left(1-\frac{{r}^{2}}{{R}^{2}}\right)$

(21)

${M}_{r}={M}_{\theta }=\frac{EW{h}^{3}}{6{R}^{2}\left(1-\mu\right)}$

(22)

板的变形能 $\mathrm{\delta }{U}_{{\mathrm{E}}_{1}}$ 由弯矩引起的弯曲应变能 $\delta {U}_{\mathrm{E}{\mathrm{M}}_{1}}$ 和膜力引起的伸长应变能 $\delta {U}_{\mathrm{E}{\mathrm{N}}_{1}}$ 组成

$\delta {U}_{{\mathrm{E}}_{1}}=\delta {U}_{\mathrm{E}{\mathrm{M}}_{1}}+\delta {U}_{\mathrm{E}{\mathrm{N}}_{1}}$

(23)

$\delta {U}_{\mathrm{E}{\mathrm{M}}_{1}}=2\text{π} {\int }_{{R}_{\mathrm{P}}}^{R}\left({2M}_{r}\mathrm{\delta }{\kappa }_{r}\right)r{\rm d}r,\quad\delta {U}_{\mathrm{E}{\mathrm{N}}_{1}}=2\text{π} {\int }_{{R}_{\mathrm{P}}}^{R}\left({2N}_{r}\mathrm{\delta }{\varepsilon }_{r}\right)r{\rm d}r$

(24)

分别计算其等效刚度

${K}_{\mathrm{E}{\mathrm{M}}_{1}}=\frac{{\partial }^{2}{U}_{\mathrm{E}{\mathrm{M}}_{1}}}{\partial W{\left(t\right)}^{2}},\quad {K}_{\mathrm{E}{\mathrm{N}}_{1}}=\frac{{\partial }^{2}{U}_{\mathrm{E}{\mathrm{N}}_{1}}}{\partial W{\left(t\right)}^{2}}$

(25)

当圆板完全处于弹性状态时， ${R}_{\mathrm{P}}=0$ ；当圆板完全处于塑性状态时， ${R}_{\mathrm{P}}=R$ 。

1.3.2 弯矩和膜力联合承载的塑性区加载

弯矩和膜力联合参与承载的塑性区为 ${R}_{\mathrm{P}\mathrm{N}}\left(W\right)$ <r≤ ${R}_{\mathrm{P}}\left(W\right)$ ，此时 ${N}_{r}\left({N}_{\theta }\right) < {n}_{0}$ ， ${M}_{r}\left({M}_{\theta }\right) < {m}_{0}$ 。区域内满足

${N}_{r}={N}_{\theta }=\frac{Eh}{1-\mu}{\left(\frac{W}{R}\right)}^{2}\left(1-\frac{{r}^{2}}{{R}^{2}}\right)$

(26)

${M}_{r}={M}_{\theta }={m}_{0}\left[1-{\left(\frac{{N}_{r}}{{n}_{0}}\right)}^{2}\right]=\frac{{h}^{2}{\sigma }_{0}}{4}\left[1-{\left(\frac{E}{1-\mu}\right)}^{2}{\left(\frac{W}{R}\right)}^{4}{\left(1-\frac{{r}^{2}}{{R}^{2}}\right)}^{2}\right]$

(27)

板的变形能 ${\mathrm{\delta }U}_{{\mathrm{P}}_{1}}$ 由弯矩引起的弯曲应变能 $\mathrm{\delta }{U}_{\mathrm{P}{\mathrm{M}}_{1}}$ 和膜力引起的伸长应变能 $\mathrm{\delta }{U}_{\mathrm{P}{\mathrm{N}}_{1}}$ 组成

${\mathrm{\delta }U}_{{\mathrm{P}}_{1}}=\delta {U}_{\mathrm{P}{\mathrm{M}}_{1}}+\delta {U}_{\mathrm{P}{\mathrm{N}}_{1}}$

(28)

$\delta {U}_{\mathrm{P}{\mathrm{M}}_{1}}=2\text{π} {\int }_{{R}_{\mathrm{P}\mathrm{N}}}^{{R}_{\mathrm{P}}}\left(2{M}_{r}\mathrm{\delta }{\kappa }_{r}\right)r{\rm d}r,\quad\delta {U}_{\mathrm{P}{\mathrm{N}}_{1}}=2\text{π} {\int }_{{R}_{\mathrm{P}\mathrm{N}}}^{{R}_{\mathrm{P}}}\left({2N}_{r}\delta {\varepsilon }_{r}\right)r{\rm d}r$

(29)

分别计算其等效抗力

${f}_{\mathrm{P}{\mathrm{M}}_{1}}=\frac{\partial {U}_{\mathrm{P}{\mathrm{M}}_{1}}}{\partial W(t)},\quad{f}_{\mathrm{P}{\mathrm{N}}_{1}}=\frac{\partial {U}_{\mathrm{P}{\mathrm{N}}_{1}}}{\partial W(t)}$

(30)

当圆板完全处于塑性状态时， ${R}_{\mathrm{P}}=R$ 。

1.3.3 塑性膜区加载

塑性膜区为0≤r≤ ${R}_{\mathrm{P}\mathrm{N}}(W)$ ，此时 $N={n}_{0}$ ， $M=0$ ，区域内满足

$\delta {U}_{{\mathrm{P}}_{2}}=\delta {U}_{\mathrm{P}{\mathrm{N}}_{2}}=2\text{π} {\int }_{0}^{{R}_{\mathrm{P}\mathrm{N}}}\left(2{n}_{0}\mathrm{\delta }{\varepsilon}_{r}\right)r{\rm d}r$

(31)

其等效抗力为

${f}_{\mathrm{P}{\mathrm{N}}_{2}}={\partial {U}_{\mathrm{P}{\mathrm{N}}_{2}}}/{\partial W(t)}$

(32)

1.3.4 边　界

板的固支边界存在结构能量耗散。如果边界达到其塑性状态，即 ${M}_{\mathrm{B}}={m}_{0}$ 时，边界成为塑性铰线，将发生有限旋转，引起塑性能量耗散。因此，通过在 $r=R$ 处引入弯矩为 ${M}_{\mathrm{B}}$ 的理想弹塑性弹簧^{[10, 16]}来模拟固支边界

${M}_{\mathrm{B}}=\left\{ \begin{array}{l}\dfrac{EW{h}^{3}}{6{R}^{2}\left(1-\mu\right)}\quad W < {W}_{\mathrm{E}\mathrm{M}\mathrm{A}\mathrm{X}}\\ \dfrac{{\sigma }_{0}{h}^{2}}{4}\quad\quad\;\;\quad W{\text{≥}} {W}_{\mathrm{E}\mathrm{M}\mathrm{A}\mathrm{X}}\end{array}\right.$

(33)

式中： ${W}_{\mathrm{E}\mathrm{M}\mathrm{A}\mathrm{X}}$ 为边界成为塑性铰线时的中心点挠度。

${W}_{\mathrm{E}\mathrm{M}\mathrm{A}\mathrm{X}}={3{R}^{2}\left(1-\mu\right)}/({2EWh})$

(34)

因此，边界处耗散能量的变化为

$\delta {U}_{\mathrm{B}}=2\text{π} R{M}_{\mathrm{B}}\delta {\theta }_{\mathrm{B}}$

(35)

${\theta }_{\mathrm{B}}={2W}/{R}$

(36)

式中： ${\theta }_{\mathrm{B}}$ 为边界处的旋转角度。

其等效抗力为

${f}_{\mathrm{B}}={\partial {U}_{\mathrm{B}}}/{\partial W(t)}$

(37)

1.4 卸载阶段的应力场和变形能场

在爆炸加载过程中，弹性变形和塑性变形均随外力增大而增大，外力卸载后，弹性变形可以恢复，塑性变形不能恢复，且塑性变形不积累变形能。与弹塑性有限元原理类似，本研究假设弹性变形与塑性变形不耦合，弹性变形积累的弹性势能与塑性变形耗散能互不影响，具有各自独立的分布规律，变形也有各自独立的分布规律，即

$\varepsilon ={\varepsilon }_{\mathrm{E}}+{\varepsilon }_{\mathrm{P}}$

(38)

由于假设圆板横向变形时刻满足式 $\left(3\right)$ ，因此满足

$W={W}_{\mathrm{E}\mathrm{E}}+{W}_{\mathrm{E}\mathrm{P}}$

(39)

式中： ${W}_{\mathrm{E}\mathrm{E}}$ 为中心点等效弹性挠度， ${W}_{\mathrm{E}\mathrm{P}}$ 为中心点等效塑性挠度。 ${w}_{\mathrm{E}\mathrm{E}}$ 和 ${w}_{\mathrm{E}\mathrm{P}}$ 分别满足

${w}_{\mathrm{E}\mathrm{E}}={W}_{\mathrm{E}\mathrm{E}}\left[1-{\left(\frac{r}{R}\right)}^{2}\right],\quad{w}_{\mathrm{E}\mathrm{P}}={W}_{\mathrm{E}\mathrm{P}}\left[1-{\left(\frac{r}{R}\right)}^{2}\right]$

(40)

卸载过程为纯弹性变形过程，因此在考虑结构塑性变形后的几何形状基础上，仅考虑弹性变形引起的应力场大小与分布。由于卸载阶段弹性能的释放是通过膜力和弯矩的减小实现的，且卸载阶段应力-应变关系仍符合材料原有弹性阶段的规律，因此卸载过程可以看作对圆板施加假想的等大反向的膜力和弯矩加载过程^[17-18]，其等效质量-弹簧-阻尼模型如图7所示，其中： ${K}'_{\mathrm{E}{\mathrm{M}}_{1}}$ 、 ${K}'_{\mathrm{E}{\mathrm{N}}_{1}}$ 、 ${K}'_{\mathrm{E}{\mathrm{M}}_{2}}$ 、 ${K}'_{\mathrm{E}{\mathrm{N}}_{2}}$ 和 ${K}'_{\mathrm{E}{\mathrm{N}}_{3}}$ 为卸载阶段的等效刚度， ${f}'_{\mathrm{c}}$ 和 ${f}'_{\mathrm{B}}$ 为卸载阶段的等效抗力，上标“′”表示卸载过程。

图 7 卸载阶段的等效质量-弹簧-阻尼模型

Figure 7. Equivalent mass-spring-damper model in unloading phase

下载: 全尺寸图片幻灯片

1.4.1 弹性区卸载

弹性区卸载的应力场满足

${\varepsilon '_{r}}={\varepsilon '_{\theta }}=\frac{{W}_{\mathrm{E}\mathrm{E}\mathrm{M}\mathrm{A}\mathrm{X}}^{2}-{\left(W-{W}_{\mathrm{E}\mathrm{P}\mathrm{M}\mathrm{A}\mathrm{X}}\right)}^{2}}{{R}^{2}}\left(1-\frac{{r}^{2}}{{R}^{2}}\right)$

(41)

${N}'_{r}={N}'_{\theta }=\frac{EH\left[{W}_{\mathrm{E}\mathrm{E}\mathrm{M}\mathrm{A}\mathrm{X}}^{2}-{\left(W-{W}_{\mathrm{E}\mathrm{P}\mathrm{M}\mathrm{A}\mathrm{X}}\right)}^{2}\right]}{\left(1-\mu\right)\left({R}^{2}+{W}^{2}\right)}\left(1-\frac{{r}^{2}}{{R}^{2}}\right)$

(42)

${\kappa }'_{r}={\kappa }'_{\theta }=\frac{2\left({W}_{\mathrm{M}\mathrm{A}\mathrm{X}}-W\right)}{{R}^{2}}$

(43)

${M}'_{r}={M}'_{\theta }=\frac{E{H}^{3}{R}^{4}\left({W}_{\mathrm{M}\mathrm{A}\mathrm{X}}-W\right)}{6\left(1-\mu\right){\left({R}^{2}+{W}^{2}\right)}^{3}}$

(44)

式中： ${\varepsilon }'_{r}$ 和 ${\varepsilon }'_{\theta }$ 为卸载阶段的应变， ${N}'_{r}$ 和 ${N}'_{\theta }$ 为卸载阶段的膜力， ${\kappa }'_{r}$ 和 ${\kappa }'_{\theta }$ 为卸载阶段的曲率， ${M}'_{r}$ 和 ${M}'_{\theta }$ 为卸载阶段的弯矩， ${W}_{\mathrm{E}\mathrm{E}\mathrm{M}\mathrm{A}\mathrm{X}}$ 为中心点等效弹性挠度的最大值， ${W}_{\mathrm{E}\mathrm{P}\mathrm{M}\mathrm{A}\mathrm{X}}$ 为中心点等效塑性挠度的最大值， ${W}_{\mathrm{M}\mathrm{A}\mathrm{X}}$ 为中心点挠度的最大值。

板弹性变形能的释放分别通过弯曲应变能的释放和伸长应变能的释放进行

$\delta {U}'_{{\mathrm{E}}_{1}}={\mathrm{\delta }U}'_{\mathrm{E}{\mathrm{N}}_{1}}+\delta {U}'_{\mathrm{E}{\mathrm{M}}_{1}}$

(45)

${\mathrm{\delta }U}'_{\mathrm{E}{\mathrm{N}}_{1}}=2\text{π} {\int }_{{R}_{\mathrm{P}\mathrm{M}\mathrm{A}\mathrm{X}}}^{R}\left(2{N}'_{r}\mathrm{\delta }{\varepsilon }'_{r}\right)r\mathrm{d}r,\quad\delta {U}'_{\mathrm{E}{\mathrm{M}}_{1}}=2\text{π} {\int }_{{R}_{\mathrm{P}\mathrm{M}\mathrm{A}\mathrm{X}}}^{R}\left(2{M}'_{r}\mathrm{\delta }{\kappa }'_{r}\right)r\mathrm{d}r$

(46)

式中： ${R}_{\mathrm{P}\mathrm{M}\mathrm{A}\mathrm{X}}$ 为卸载前 ${R}_{\mathrm{P}}$ 的最终值。分别计算其等效刚度

${K}'_{\mathrm{E}{\mathrm{N}}_{1}}=\frac{{\partial }^{2}{U}'_{\mathrm{E}{\mathrm{N}}_{1}}}{\partial W{\left(t\right)}^{2}},\quad{K}'_{\mathrm{E}{\mathrm{M}}_{1}}=\frac{{\partial }^{2}{U}'_{\mathrm{E}{\mathrm{M}}_{1}}}{\partial W{\left(t\right)}^{2}}$

(47)

1.4.2 膜力和弯矩联合承载的塑性区卸载

板弹性变形能的释放分别通过弯曲应变能的释放和伸长应变能的释放进行

$\delta {U}'_{{\mathrm{E}}_{2}}={\mathrm{\delta }U}'_{\mathrm{E}{\mathrm{N}}_{2}}+\delta {U}'_{\mathrm{E}{\mathrm{M}}_{2}}$

(48)

${\mathrm{\delta }U}'_{\mathrm{E}{\mathrm{N}}_{2}}=2\text{π} {\int }_{{R}_{\mathrm{P}\mathrm{N}\mathrm{M}\mathrm{A}\mathrm{X}}}^{{R}_{\mathrm{P}\mathrm{M}\mathrm{A}\mathrm{X}}}\left(2{N}'_{r}\mathrm{\delta }{\varepsilon }'_{r}\right)r\mathrm{d}r,\quad\delta {U}'_{\mathrm{E}{\mathrm{M}}_{2}}=2\text{π} {\int }_{{R}_{\mathrm{P}\mathrm{N}\mathrm{M}\mathrm{A}\mathrm{X}}}^{{R}_{\mathrm{P}\mathrm{M}\mathrm{A}\mathrm{X}}}\left(2{M}'_{r}\mathrm{\delta }{\kappa }'_{r}\right)r\mathrm{d}r$

(49)

式中： ${R}_{\mathrm{P}\mathrm{N}\mathrm{M}\mathrm{A}\mathrm{X}}$ 为卸载前 ${R}_{\mathrm{P}\mathrm{N}}$ 的最终值。分别计算其等效刚度

${K}'_{\mathrm{E}{\mathrm{N}}_{2}}=\frac{{\partial }^{2}{U}'_{\mathrm{E}{\mathrm{N}}_{2}}}{\partial W{\left(t\right)}^{2}},\quad{K}'_{\mathrm{E}{\mathrm{M}}_{2}}=\frac{{\partial }^{2}{U}'_{\mathrm{E}{\mathrm{M}}_{2}}}{\partial W{\left(t\right)}^{2}}$

(50)

1.4.3 塑性膜区卸载

塑性膜区仅考虑伸长应变能的释放

$\delta {U}'_{{\mathrm{E}}_{3}}={\mathrm{\delta }U}'_{\mathrm{E}{\mathrm{N}}_{3}}=2\text{π} {\int }_{0}^{{R}_{\mathrm{P}\mathrm{N}\mathrm{M}\mathrm{A}\mathrm{X}}}\left(2{N}'_{r}\mathrm{\delta }{\varepsilon }'_{r}\right)r\mathrm{d}r$

(51)

计算其等效刚度

${K}'_{\mathrm{E}{\mathrm{N}}_{3}}={{\partial }^{2}{U}'_{\mathrm{E}{\mathrm{N}}_{3}}}/{\partial W{\left(t\right)}^{2}}$

(52)

1.4.4 边　界

边界处弯矩

${M}'_{\mathrm{B}}=\left\{\begin{array}{l}\dfrac{EW{h}^{3}}{6{R}^{2}\left(1-\mu\right)}\;\;\;\;{W}_{\mathrm{M}\mathrm{A}\mathrm{X}} < {W}_{\mathrm{E}\mathrm{M}\mathrm{A}\mathrm{X}}\\ \dfrac{{\sigma }_{0}{h}^{2}}{4}\;\;\;\;\;\;\;\;\;\;\;\;\;\;{W}_{\mathrm{M}\mathrm{A}\mathrm{X}}\text{≥} {W}_{\mathrm{E}\mathrm{M}\mathrm{A}\mathrm{X}}\end{array}\right.$

(53)

边界处旋转角度

${\mathrm{\delta }\theta }'_{\mathrm{B}}=\frac{2\mathrm{\delta }\left({W}_{\mathrm{M}\mathrm{A}\mathrm{X}}-W\right)}{R}=-\frac{2\mathrm{\delta }W}{R}$

(54)

边界处耗散能量的变化

$\delta {U}'_{\mathrm{B}}=2\text{π} R{M}'_{\mathrm{B}}\delta {\theta }'_{\mathrm{B}}$

(55)

其等效抗力为

${f}'_{\mathrm{B}}={\partial {U}'_{\mathrm{B}}}/{\partial W(t)}$

(56)

1.4.5 黏性阻尼

采用黏性阻尼模型描述圆板振动中受到的介质阻尼、结构阻尼和材料阻尼。黏性阻尼力与相对速度成正比，即

${f}'_{\mathrm{c}}=c\dot{W}(t)$

(57)

式中：c为阻尼系数。

$c=2\zeta \sqrt{Km}=2\zeta \sqrt{\left({K}'_{\mathrm{E}{\mathrm{N}}_{1}}+{K}'_{\mathrm{E}{\mathrm{M}}_{1}}+{K}'_{\mathrm{E}{\mathrm{N}}_{2}}+{K}'_{\mathrm{E}{\mathrm{M}}_{2}}+{K}'_{\mathrm{E}{\mathrm{N}}_{3}}\right)m}$

(58)

式中： $\zeta$ 为相对阻尼系数。综合考虑材料阻尼、结构阻尼以及介质阻尼，将 $\zeta$ 取为 $0.16$ 。

1.4.6 能量守恒

${W}_{\mathrm{E}\mathrm{E}\mathrm{M}\mathrm{A}\mathrm{X}}$ 和 ${W}_{\mathrm{E}\mathrm{P}\mathrm{M}\mathrm{A}\mathrm{X}}$ 为未知量，需要通过能量守恒计算得到。假设加载过程中积累的弹性能在卸载过程中全部转化为释放的弹性变形能、边界耗散和黏性耗散，黏性阻尼力通过卸载过程的平均速度近似计算得到，由此求解出 ${W}_{\mathrm{E}\mathrm{E}\mathrm{M}\mathrm{A}\mathrm{X}}$ 和 ${W}_{\mathrm{E}\mathrm{P}\mathrm{M}\mathrm{A}\mathrm{X}}$ ，进而得到卸载阶段各等效刚度。

1.5 外力功等效

根据SDOF模型，等效载荷所做的虚功与圆板表面等效均匀载荷所做虚功相等，求解加载的广义力

$Q(t)\delta W={\int }_{0}^{2\text{π}}{\int }_{0}^{R}p(t)\mathrm{\delta }W\left(1+\frac{{W}^{2}}{{R}^{2}}\right)\left(1-\frac{{r}^{2}}{{R}^{2}}\right)r\mathrm{d}r\mathrm{d}\theta =\frac{\text{π}\mathrm{\delta }W{R}^{2}p(t)}{2}\left(1+\frac{{W}^{2}}{{R}^{2}}\right)$

(59)

$Q(t)=\frac{\text{π}{R}^{2}}{2}\left(1+\frac{{W}^{2}}{{R}^{2}}\right)p(t)$

(60)

式中：p $\left(t\right)$ 为载荷经验公式描述的爆炸冲击波的压力-时间关系。炸药爆炸后，冲击波由爆炸源向外传播，压力先迅速增加到 ${p}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ ，随后呈指数衰减。常用的Friedlander方程可以较好地描述爆炸冲击波的压力-时间关系

$p(t)=\left\{\begin{array}{l}{p}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(1-\dfrac{t-{t}_{\mathrm{a}}}{{t}_{\mathrm{d}}}\right){\rm {exp}}\left(-b\dfrac{t-{t}_{\mathrm{a}}}{{t}_{\mathrm{d}}}\right)\quad0\text{≤} t\text{≤} {t}_{\mathrm{d}}\\ 0\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\; t > {t}_{\mathrm{d}}\end{array}\right.$

(61)

式中： ${p}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ 为压力峰值， ${t}_{\mathrm{a}}$ 为冲击波到达时间， ${t}_{\mathrm{d}}$ 为冲击波正相持续时间， $b$ 为波形参数。

1.6 质量等效

根据SDOF模型，动能与圆板动能相等，求解等效质量

$\frac{1}{2}m{\dot{W}}^{2}(t)=\frac{1}{2}\rho H{\int }_{0}^{2\text{π}}{\int }_{0}^{R}{\left[\dot{W}(t)\left(1-\frac{{r}^{2}}{{R}^{2}}\right)\right]}^{2}r\mathrm{d}r\mathrm{d}\theta =\frac{1}{6}\text{π} \rho H{R}^{2}{\dot{W}}^{2}(t)$

(62)

$m=\frac{\text{π}\rho H{R}^{2}}{3}$

(63)

针对变形的不同阶段，分析其应力状态，求出SDOF模型的等效载荷、等效质量、等效刚度、等效抗力，代入SDOF模型的运动方程进行求解，通过多个阶段的叠加，可以得到固支弹塑性圆板在空爆作用下的完整动态响应。

2. 有限元数值模拟

基于高强度圆形钢板空中近爆试验工况^[12]，建立有限元模型，进行数值模拟计算，并与试验结果^[12]对比，验证有限元模拟结果的有效性。

试验中使用球形TNT裸装药，从炸药中心点起爆。靶板由上下压板夹紧固支，上压板的厚度为50 mm，下压板的厚度为80 mm，靶板与上下压板的尺寸均为 $1.5\;\mathrm{m}\mathrm{m}\times 1.5\;\mathrm{m}$ ，靶板受空爆载荷作用的范围为半径0.5 m的圆。试验布置见图8。

图 8 试验布置

Figure 8. Test setup

下载: 全尺寸图片幻灯片

试验工况如表1所示。

表 1 3种高强度钢的试验参数^[12]

Table 1. Test parameters of three kinds of high strength steel^[12]

Case	Material	Thickness/mm	Yield strength/MPa	Mass of TNT/kg	Detonation distance/mm
1	Weldox700E steel	8	800	6	250
2	B1 steel	8	1070	6	250
3	B2 steel	8	1600	6	250
4	Weldox700E steel	12	800	10	250
5	B1 steel	12	1070	10	250
6	B2 steel	12	1600	10	250

下载: 导出CSV

| 显示表格

按照试验工况，通过显式有限元动力分析软件LS-DYNA建立有限元模型，对TNT爆炸冲击波在空气中的传播及其与固支圆板的相互作用进行模拟。空气和上下压板采用六面体网格划分，压板采用刚体材料模型描述，靶板采用二维网格划分，基于任意拉格朗日-欧拉（arbitrary Lagrangian-Eulerian，ALE）算法进行数值模拟分析。图9为有限元模型示意图。各工况下模拟得到的中心点位移-时间曲线如图10所示。

图 9 有限元模型

Figure 9. Finite element model

下载: 全尺寸图片幻灯片

图 10 各工况下有限元数值模拟得到的中心点位移-时间曲线

Figure 10. Displacement-time curves of center point by finite element simulation under various calculation cases

下载: 全尺寸图片幻灯片

数值模拟计算得到的靶板中心点最大位移和最终位移与试验结果的对比如表2所示。所有工况下有限元数值模拟结果与试验结果的相对误差均在10%以内，说明有限元模拟的空中近爆作用下固支圆板的变形结果与试验结果较吻合。

表 2 各工况下试验与有限元模拟得到的中心点位移对比

Table 2. Comparison of center point displacements between test and finite element simulation under various calculation cases

Case	Maximum displacement of center point			Final displacement of center point
Case	Test/mm	Sim./mm	Error/%	Test/mm	Sim./mm	Error/%
1	145	148.84	2.65	129	137.42	6.53
2	128	134.67	5.21	102	111.88	9.69
3	103	99.52	−3.38	85	79.07	−6.98
4	166	152.44	−8.17	153	141.95	−7.22
5	148	133.87	−9.55	118	120.12	1.80
6	115	106.74	−7.18	83	79.28	−4.48

下载: 导出CSV

| 显示表格

3. 动力学模型验证

根据有限元计算结果，选取6和10 kg TNT空爆作用于圆板时距圆心0～0.5 mm范围的峰值压力载荷，拟合得到的峰值压力分布函数如图11所示。

图 11 峰值压力点及其拟合函数

Figure 11. Peak pressure points and fitting functions

下载: 全尺寸图片幻灯片

将拟合出的函数代入式 $\left(2\right)$ 得到等效均匀载荷峰值，并代入式 $\left(61\right)$ 可得等效均匀载荷 $p\left(t\right)$ 。经过拟合，得到 ${t}_{\mathrm{a}}=0$ ， $b=1$ ， ${t}_{\mathrm{d}}=5\times {10}^{-4}$ s，6和10 kg TNT 空爆下的等效均匀峰值压力分别为48.3和83.7 MPa。

以工况1为例，通过计算 ${W}_{\mathrm{P}\mathrm{M}}$ 、 ${W}_{\mathrm{P}\mathrm{B}}$ 、 ${W}_{\mathrm{P}\mathrm{N}\mathrm{M}}$ ，可将变形过程分为5个阶段。

(1) 阶段1（0～W_PM）

在此阶段，整个圆板均未屈服，完全处于弹性状态，将等效载荷、等效质量、等效刚度、等效抗力代入运动方程

$m\ddot{W}(t)+\left({K}_{\mathrm{E}{\mathrm{M}}_{1}}+{K}_{\mathrm{E}{\mathrm{N}}_{1}}\right)W(t)+{f}_{\mathrm{B}}=\frac{\text{π}{R}^{2}}{2}\left(1+\frac{{W}^{2}}{{R}^{2}}\right)p(t)$

(64)

初始条件为 $W\left(t=0\right)=0$ ， $\dot{W}\left(t=0\right)=0$ 。

(2) 阶段2（W_PM～W_PNM）

在此阶段，膜力弯矩联合承载的塑性区逐渐由圆板中心向边界传播，运动方程为

$m\ddot{W}(t)+\left({K}_{\mathrm{E}{\mathrm{M}}_{1}}+{K}_{\mathrm{E}{\mathrm{N}}_{1}}\right)W(t)+{f}_{\mathrm{P}{\mathrm{M}}_{1}}+{f}_{\mathrm{P}{\mathrm{N}}_{1}}+{f}_{\mathrm{B}}=\frac{\text{π}{R}^{2}}{2}\left(1+\frac{{W}^{2}}{{R}^{2}}\right)p(t)$

(65)

初始条件为 $W({t}_{\mathrm{P}\mathrm{M}})={W}_{\mathrm{P}\mathrm{M}},\dot{W}({t}_{\mathrm{P}\mathrm{M}})={v}_{\mathrm{P}\mathrm{M}}$ ，其中 ${t}_{\mathrm{P}\mathrm{M}}$ 为阶段1的结束时刻， ${v}_{\mathrm{P}\mathrm{M}}$ 为阶段1结束时圆板中心点速度。

(3) 阶段3（W_PNM～W_PB）

在此阶段，塑性膜区逐渐由圆板中心向边界传播，但弹性区仍未消失，运动方程为

$m\ddot{W}(t)+\left({K}_{\mathrm{E}{\mathrm{M}}_{1}}+{K}_{\mathrm{E}{\mathrm{N}}_{1}}\right)W(t)+{f}_{\mathrm{P}{\mathrm{M}}_{1}}+{f}_{\mathrm{P}{\mathrm{N}}_{1}}+{f}_{\mathrm{P}{\mathrm{N}}_{2}}+{f}_{\mathrm{B}}=\frac{\text{π}{R}^{2}}{2}\left(1+\frac{{W}^{2}}{{R}^{2}}\right)p(t)$

(66)

初始条件为 $W({t}_{\mathrm{P}\mathrm{N}\mathrm{M}})={W}_{\mathrm{P}\mathrm{N}\mathrm{M}},\dot{W}({t}_{\mathrm{P}\mathrm{N}\mathrm{M}})={v}_{\mathrm{P}\mathrm{N}\mathrm{M}}$ ，其中 ${t}_{\mathrm{P}\mathrm{N}\mathrm{M}}$ 为阶段2的结束时刻， ${v}_{\mathrm{P}\mathrm{N}\mathrm{M}}$ 为阶段2结束时圆板中心点速度。

(4) 阶段4（W_PB～W_MAX）

在此阶段，塑性区完全替代弹性区，运动方程为

$m\ddot{W}(t)+{f}_{\mathrm{P}{\mathrm{M}}_{1}}+{f}_{\mathrm{P}{\mathrm{N}}_{1}}+{f}_{\mathrm{P}{\mathrm{N}}_{2}}+{f}_{\mathrm{B}}=\frac{\text{π}{R}^{2}}{2}\left(1+\frac{{W}^{2}}{{R}^{2}}\right)p(t)$

(67)

初始条件为 $W({t}_{\mathrm{P}\mathrm{B}})={W}_{\mathrm{P}\mathrm{B}},\dot{W}({t}_{\mathrm{P}\mathrm{B}})={v}_{\mathrm{P}\mathrm{B}}$ ，其中 ${t}_{\mathrm{P}\mathrm{B}}$ 为阶段3的结束时刻， ${v}_{\mathrm{P}\mathrm{B}}$ 为阶段3结束时圆板中心点速度。

(5) 阶段5（卸载阶段）

在此阶段开始时刻，圆板中心点位移已经运动到最大，随后在弹簧和阻尼的作用下发生欠阻尼自由振动，运动方程为

$m\ddot{W}(t)+\left({K}'_{\mathrm{E}{\mathrm{N}}_{2}}+{K}'_{\mathrm{E}{\mathrm{M}}_{2}}+{K}'_{\mathrm{E}{\mathrm{N}}_{3}}\right)W(t)+2\zeta \sqrt{m\left({K}'_{\mathrm{E}{\mathrm{N}}_{2}}+{K}'_{\mathrm{E}{\mathrm{M}}_{2}}+{K}'_{\mathrm{E}{\mathrm{N}}_{3}}\right)}\dot{W}(t)+{f}'_{\mathrm{B}}=0$

(68)

初始条件为 $W({t}_{\mathrm{M}\mathrm{A}\mathrm{X}})={W}_{\mathrm{M}\mathrm{A}\mathrm{X}},\dot{W}({t}_{\mathrm{M}\mathrm{A}\mathrm{X}})=0$ ，其中 ${t}_{\mathrm{M}\mathrm{A}\mathrm{X}}$ 为阶段4的结束时刻。

由此，通过多个阶段的叠加可以得到工况1时固支弹塑性圆板在空爆作用下的完整动态响应。分别计算6种工况下弹塑性圆板中心点位移随时间的变化曲线，并与数值模拟得到的圆板中心点挠度进行对比，如图12所示。从图12可以看出，在相同当量的爆炸载荷作用下，弹塑性圆板材料的屈服强度越高，其中心点的最大位移和最终位移越小。这是由于材料的屈服强度越高，圆板加载后期的等效刚度和等效抗力越大，中心点速度衰减得越快，从而导致中心点的最大位移和最终位移越小。近爆条件下固支弹塑性圆板的变形过程受多种因素的影响，因此，有限元数值模拟得到的弹塑性圆板中心点位移随时间变化曲线不是振幅逐渐衰减的理想的简谐振动曲线。

图 12 各工况下动力学模型计算结果与试验及有限元数值模拟结果的对比

Figure 12. Comparison of dynamic model results with test and finite element simulation results under various cases

下载: 全尺寸图片幻灯片

通过动力学模型计算各个工况的中心点位移，并与试验结果进行对比，结果如表3所示。从表3可以看出，6种工况下SDOF模型计算结果与试验结果之间的误差表现出相似的趋势，即计算得到的圆板中心点最大位移相对于试验值偏小，而最终位移相对于试验值偏大，导致这些误差的原因可能有以下3方面：(1) SDOF模型采用固支边界条件，然而在试验中很难实现对靶板的理想刚性约束，往往会在边界处产生一定的滑移，另外，多次爆炸冲击作用后，试验台架的约束作用不可避免地减弱，导致试验测量的中心点最大位移偏大；(2) 试验中空爆载荷并非均匀载荷，爆轰产物首先作用于圆板中心，随后逐渐向周围扩散，且圆板的实际变形过程并非严格遵循形函数，导致动力学模型计算得到的中心点最大位移小于实际最大位移；(3) 由于实际的金属材料存在不同程度的应变强化效应和应变率效应，而动力学模型将材料简化为理想弹塑性本构模型，计算出的等效抗力偏小，导致求解的中心点最终位移相对于试验位移偏大。

表 3 各工况下动力学模型计算结果与试验结果的对比

Table 3. Comparison between dynamic model results and test results under various cases

Case	Maximum displacement			Final displacement
Case	Test/mm	SDOF model/mm	Error/%	Test	SDOF model/mm	Error/%
1	145	141.48	−2.43	129	129.22	0.17
2	128	126.03	−1.54	102	110.06	7.90
3	103	101.49	−1.47	85	86.21	1.42
4	166	160.24	−3.47	153	148.57	−2.90
5	148	141.27	−4.55	118	127.66	8.19
6	115	114.19	−0.70	83	95.03	14.49

下载: 导出CSV

| 显示表格

将动力学模型计算的中心点最终位移代入式 $\left(3\right)$ ，得到靶板变形后的轮廓线，并与试验轮廓线进行对比，结果如图13所示。从图13可以看出，各工况下动力学模型计算得到的靶板轮廓线与试验结果拟合较好，因此，式 $\left(3\right)$ 可以较好地描述空中近爆作用下固支圆板的变形轮廓。

图 13 各工况下动力学模型与试验得到的靶板轮廓线的对比

Figure 13. Comparison of circular plate deformations between dynamic model and test under various cases

下载: 全尺寸图片幻灯片

总的来说，动力学模型计算的中心点位移与试验结果的最大相对误差未超过15%，符合工程模型的误差要求，可为结构抗爆设计提供快速算法。

4. 结　论

(1) 通过6种工况的对比分析可知，等效SDOF模型计算结果与试验及数值模拟结果的偏差在合理范围内，采用等效SDOF方法可以较准确地预测空爆载荷作用下固支圆板的变形情况。考虑到SDOF模型的计算效率比有限元数值模拟的计算效率高，因此，本研究发展的单自由度方法可为装甲车辆防爆结构设计提供技术支撑。

(2) 质量-弹簧-阻尼模型的计算结果是振幅逐渐衰减的简谐振动曲线。然而，由于固支圆板的实际运动过程受多种因素影响，如试验中很难做到刚性约束，台架多次爆炸冲击后难免有变形，边界约束能力减弱等，实际钢板变形并非严格按照形函数发生，并且实际材料并非理想的弹塑性，存在不同程度的应变强化效应和应变率效应，因此，SDOF模型对系统固有频率和振幅的估计存在一定误差。

(3) 如果将固支圆板看作一个系统，空爆载荷视为系统的输入，固支圆板上各点发生的振动视为系统的输出，那么，系统动力响应分析的准确与否不仅取决于系统描述的准确性，同时还取决于输入载荷描述的准确性。由于精确描述近场爆炸载荷十分困难，因此，通过LS-DYNA分析结果拟合得到作用在圆板上的爆炸载荷是一种可行的方法。当然，改善爆炸载荷模型精度的方法都有益于提升本研究所建立的SDOF模型的预测结果。

图 1 载荷等效示意图

Figure 1. Equivalent load diagram

下载: 全尺寸图片幻灯片

图 2 式 $\left(4\right)$ 和式 $\left(5\right)$ 的计算结果对比

Figure 2. Comparison of calculation results between Eq. (4) and Eq. (5)

下载: 全尺寸图片幻灯片

图 3 板元受力分析

Figure 3. Force analysis of plate microelement

下载: 全尺寸图片幻灯片

图 4 弯矩和膜力联合承载的屈服条件

Figure 4. Yield condition of combined bending moment and membrane force

下载: 全尺寸图片幻灯片

图 5 固支圆板发生弹塑性大变形后的应力场分布

Figure 5. Stress field of clamped circular plate after large elastic-plastic deformation

下载: 全尺寸图片幻灯片

图 6 加载阶段的等效质量-弹簧-阻尼模型

Figure 6. Equivalent mass-spring-damper model in loading phase

下载: 全尺寸图片幻灯片

图 7 卸载阶段的等效质量-弹簧-阻尼模型

Figure 7. Equivalent mass-spring-damper model in unloading phase

下载: 全尺寸图片幻灯片

图 8 试验布置

Figure 8. Test setup

下载: 全尺寸图片幻灯片

图 9 有限元模型

Figure 9. Finite element model

下载: 全尺寸图片幻灯片

图 10 各工况下有限元数值模拟得到的中心点位移-时间曲线

Figure 10. Displacement-time curves of center point by finite element simulation under various calculation cases

下载: 全尺寸图片幻灯片

图 11 峰值压力点及其拟合函数

Figure 11. Peak pressure points and fitting functions

下载: 全尺寸图片幻灯片

图 12 各工况下动力学模型计算结果与试验及有限元数值模拟结果的对比

Figure 12. Comparison of dynamic model results with test and finite element simulation results under various cases

下载: 全尺寸图片幻灯片

图 13 各工况下动力学模型与试验得到的靶板轮廓线的对比

Figure 13. Comparison of circular plate deformations between dynamic model and test under various cases

下载: 全尺寸图片幻灯片

表 1 3种高强度钢的试验参数^[12]

Table 1. Test parameters of three kinds of high strength steel^[12]

Case	Material	Thickness/mm	Yield strength/MPa	Mass of TNT/kg	Detonation distance/mm
1	Weldox700E steel	8	800	6	250
2	B1 steel	8	1070	6	250
3	B2 steel	8	1600	6	250
4	Weldox700E steel	12	800	10	250
5	B1 steel	12	1070	10	250
6	B2 steel	12	1600	10	250

下载: 导出CSV

表 2 各工况下试验与有限元模拟得到的中心点位移对比

Table 2. Comparison of center point displacements between test and finite element simulation under various calculation cases

Case	Maximum displacement of center point			Final displacement of center point
Case	Test/mm	Sim./mm	Error/%	Test/mm	Sim./mm	Error/%
1	145	148.84	2.65	129	137.42	6.53
2	128	134.67	5.21	102	111.88	9.69
3	103	99.52	−3.38	85	79.07	−6.98
4	166	152.44	−8.17	153	141.95	−7.22
5	148	133.87	−9.55	118	120.12	1.80
6	115	106.74	−7.18	83	79.28	−4.48

下载: 导出CSV

表 3 各工况下动力学模型计算结果与试验结果的对比

Table 3. Comparison between dynamic model results and test results under various cases

Case	Maximum displacement			Final displacement
Case	Test/mm	SDOF model/mm	Error/%	Test	SDOF model/mm	Error/%
1	145	141.48	−2.43	129	129.22	0.17
2	128	126.03	−1.54	102	110.06	7.90
3	103	101.49	−1.47	85	86.21	1.42
4	166	160.24	−3.47	153	148.57	−2.90
5	148	141.27	−4.55	118	127.66	8.19
6	115	114.19	−0.70	83	95.03	14.49

下载: 导出CSV

参考文献(18)

[1]	王芳, 冯顺山, 俞为民. 爆炸冲击波作用下靶板的塑性大变形响应研究 [J]. 中国安全科学学报, 2003, 13(3): 58–61. doi: 10.16265/j.cnki.issn1003-3033.2003.03.016 WANG F, FENG S S, YU W M. Study on large plastic deformation response of target plate under explosive blast wave [J]. China Safety Science Journal, 2003, 13(3): 58–61. doi: 10.16265/j.cnki.issn1003-3033.2003.03.016
[2]	吴成, 金俨, 李华新. 固支方板对水中爆炸作用的动态响应研究 [J]. 高压物理学报, 2003, 17(4): 275–282. doi: 10.3969/j.issn.1000-5773.2003.04.006 WU C, JIN Y, LI H X. A study on square plate dynamic response under underwater explosion [J]. Chinese Journal of High Pressure Physics, 2003, 17(4): 275–282. doi: 10.3969/j.issn.1000-5773.2003.04.006
[3]	何建, 肖玉凤, 陈振勇, 等. 空爆载荷作用下固支矩形钢板的塑性极限变形 [J]. 哈尔滨工业大学学报, 2007, 39(2): 310–313, 329. doi: 10.3321/j.issn:0367-6234.2007.02.036 HE J, XIAO Y F, CHEN Z Y, et al. Plastic limited deformation analysis of the clamped rectangular steel plate subjected to air non-contact explosions [J]. Journal of Harbin Institute of Technology, 2007, 39(2): 310–313, 329. doi: 10.3321/j.issn:0367-6234.2007.02.036
[4]	DUFFEY T A. Large deflection dynamic response of clamped circular plates subjected to explosive loading [R]. Albuquerque: Sandia National Laboratories, 1967.
[5]	NURICK G N, PEARCE H T, MARTIN J B. Predictions of transverse deflections and in-plane strains in impulsively loaded thin plates [J]. International Journal of Mechanical Sciences, 1987, 29(6): 435–442. doi: 10.1016/0020-7403(87)90004-X
[6]	JONES N. Impulsive loading of a simply supported circular rigid plastic plate [J]. Journal of Applied Mechanics, 1968, 35(1): 59–65. doi: 10.1115/1.3601174
[7]	FLORENCE A L. Circular plate under a uniformly distributed impulse [J]. International Journal of Solids and Structures, 1966, 2(1): 37–47. doi: 10.1016/0020-7683(66)90005-9
[8]	CLOETE T J, NURICK G N. On the influence of radial displacements and bending strains on the large deflections of impulsively loaded circular plates [J]. International Journal of Mechanical Sciences, 2014, 82: 140–148. doi: 10.1016/j.ijmecsci.2014.02.026
[9]	CAMPBELL T I, CHARLTON T M. Finite deformation of a fully fixed beam comprised of a non-linear material [J]. International Journal of Mechanical Sciences, 1973, 15(5): 415–422. doi: 10.1016/0020-7403(73)90040-4
[10]	KARAGIOZOVA D, YU T X, SHI S Y, et al. On the influence of elasticity on the large deflections response of circular plates to uniform quasi-static pressure [J]. International Journal of Mechanical Sciences, 2017, 131/132: 894–907. doi: 10.1016/j.ijmecsci.2017.07.032
[11]	BIGGS J M. Introduction to structural dynamics [M]. New York: McGraw-Hill, 1964.
[12]	钟贤哲. 高强度钢Weldox700E的动态本构关系及其抗爆特性 [D]. 北京: 北京理工大学, 2019. ZHONG X Z. Dynamic constitutive relations and explosion resistance properties of high strength steel Weldox700E [D]. Beijing: Beijing Institute of Technology, 2019.
[13]	RIGBY S E, TYAS A, BENNETT T. Single-degree-of-freedom response of finite targets subjected to blast loading: the influence of clearing [J]. Engineering Structures, 2012, 45: 396–404. doi: 10.1016/j.engstruct.2012.06.034
[14]	TAYLOR G I. The distortion under pressure of diaphragm which is clamped along its edge and stressed beyond the elastic limit [J]. Underwater Explosion Research, 1950.
[15]	KYOHEI K, PIAN T H H. Large deformations of rigid-plastic circular plates [J]. International Journal of Solids and Structures, 1981, 17(11): 1043–1055. doi: 10.1016/0020-7683(81)90012-3
[16]	WANG X D, YU T X. Parkes revisited: effect of elastic deformation at the root of a cantilever beam [J]. International Journal of Impact Engineering, 1991, 11(2): 197–209. doi: 10.1016/0734-743X(91)90006-2
[17]	汤安民, 许明达, 赵蕾. 自由变形结构弹塑性分析的新方法 [J]. 广西大学学报(自然科学版), 2009, 34(1): 36–39. doi: 10.13624/j.cnki.issn.1001-7445.2009.01.008 TANG A M, XU M D, ZHAO L. New method for elasto-plastic analysis of free-form deformation structure [J]. Journal of Guangxi University (Natural Science Edition), 2009, 34(1): 36–39. doi: 10.13624/j.cnki.issn.1001-7445.2009.01.008
[18]	许明达. 自由变形结构弹塑性分析的新方法及应用 [D]. 西安: 西安理工大学, 2009. XU M D. The structure of free-form deformation elasto-plastic analysis and application of new method [D]. Xi’an: Xi’an University of Technology, 2009.

施引文献

资源附件(0)

访问统计

点击查看大图

图(13) / 表(3)

计量

文章访问数: 263
HTML全文浏览量: 179
PDF下载量: 44

1. 动力学模型推导
1.1 载荷等效
1.2 屈服条件
1.3 加载阶段的应力场和变形能场
1.4 卸载阶段的应力场和变形能场
1.5 外力功等效
1.6 质量等效
2. 有限元数值模拟
3. 动力学模型验证
4. 结　论

空爆载荷作用下固支弹塑性圆板的动力学模型

doi: 10.11858/gywlxb.20220525

作者简介: 吴 迪（1997－），女，硕士研究生，主要从事结构冲击动力学研究. E-mail：15010059608@163.com

通讯作者: 郭香华（1974－），男，博士，副教授，主要从事爆炸与冲击仿真、材料与结构动力学研究.E-mail：guoxh@bit.edu.cn

计量

出版历程

Dynamic Model of Clamped Elastoplastic Circular Plate under Air Blast Loading

1. 动力学模型推导

1.1 载荷等效

1.2 屈服条件

1.3 加载阶段的应力场和变形能场

1.3.1 弹性区加载

1.3.2 弯矩和膜力联合承载的塑性区加载

1.3.3 塑性膜区加载

1.3.4 边 界

1.4 卸载阶段的应力场和变形能场

1.4.1 弹性区卸载

1.4.2 膜力和弯矩联合承载的塑性区卸载

1.4.3 塑性膜区卸载

1.4.4 边 界

1.4.5 黏性阻尼

1.4.6 能量守恒

1.5 外力功等效

1.6 质量等效

2. 有限元数值模拟

3. 动力学模型验证

4. 结 论

计量

出版历程

目录

1. 动力学模型推导

1.1 载荷等效

1.2 屈服条件

1.3 加载阶段的应力场和变形能场

1.4 卸载阶段的应力场和变形能场

1.5 外力功等效

1.6 质量等效

2. 有限元数值模拟

3. 动力学模型验证

4. 结 论

作者简介:
吴　迪（1997－），女，硕士研究生，主要从事结构冲击动力学研究. E-mail：15010059608@163.com

通讯作者:
郭香华（1974－），男，博士，副教授，主要从事爆炸与冲击仿真、材料与结构动力学研究.E-mail：guoxh@bit.edu.cn

1.3.4 边　界

1.4.4 边　界

4. 结　论

4. 结　论