• ISSN 1001-1455  CN 51-1148/O3
  • EI、Scopus、CA、JST、EBSCO、DOAJ收录
  • 力学类中文核心期刊
  • 中国科技核心期刊、CSCD统计源期刊

柱形装药空中爆炸冲击波荷载研究

王明涛 程月华 吴昊

徐世烺, 吴平, 周飞, 李庆华, 曾田, 蒋霄. 活性粉末混凝土抗多次侵彻实验研究及数值预测[J]. 爆炸与冲击, 2021, 41(6): 063301. doi: 10.11883/bzycj-2020-0165
引用本文: 王明涛, 程月华, 吴昊. 柱形装药空中爆炸冲击波荷载研究[J]. 爆炸与冲击, 2024, 44(4): 043201. doi: 10.11883/bzycj-2023-0197
XU Shilang, WU Ping, ZHOU Fei, LI Qinghua, ZENG Tian, JIANG Xiao. Experimental investigation and numerical prediction on resistance of reactive powder concrete to multiple penetration[J]. Explosion And Shock Waves, 2021, 41(6): 063301. doi: 10.11883/bzycj-2020-0165
Citation: WANG Mingtao, CHENG Yuehua, WU Hao. Study on blast loadings of cylindrical charges air explosion[J]. Explosion And Shock Waves, 2024, 44(4): 043201. doi: 10.11883/bzycj-2023-0197

柱形装药空中爆炸冲击波荷载研究

doi: 10.11883/bzycj-2023-0197
基金项目: 国家自然科学基金(52078379)
详细信息
    作者简介:

    王明涛(1996-  ),男,博士研究生,wangmingtao@tongji.edu.cn

    通讯作者:

    程月华(1994-  ),女,博士,博士后,yhcheng@tongji.edu.cn

  • 中图分类号: O383

Study on blast loadings of cylindrical charges air explosion

  • 摘要: 受比例距离、装药长径比、起爆方式、方位角、冲击波入射角以及反射面相对位置等多种因素的影响,球形装药空中爆炸冲击波荷载的计算方法不适用于柱形装药。为探究柱形装药空中自由场爆炸冲击波入射和反射荷载,首先,开展柱形TNT装药单端起爆的空中爆炸试验,并基于显式动力学分析软件AUTODYN进行数值模拟,通过与试验和规范进行对比,验证了采用的有限元分析方法的适用性。进一步开展考虑比例距离、长径比、起爆方式、方位角和刚性反射等因素的1000余组柱形装药空中爆炸工况的数值模拟。基于模拟结果,揭示了柱形装药空中爆炸入射冲击波峰值超压和最大冲量及其形状因子的分布特征,提出峰值超压和最大冲量临界比例距离的判定准则和确定方法,阐明了刚性反射冲击波峰值超压和反射系数的变化规律。最后,提出柱形装药空中爆炸入射和反射冲击波荷载的计算方法,并得到360余组试验数据的验证。该方法可快速计算作用于建筑结构上的爆炸荷载,并为弹药毁伤效能评估、结构动态响应和破坏分析及其抗爆设计提供参考。
  • 近年来,随着各种破坏性武器的不断发展以及恐怖主义袭击事件不断增多[1-2],提高现有工程结构的综合防护能力已经成为防护工作的当务之急。迄今为止,大部分工程结构中使用最广泛的材料仍然是普通混凝土。然而,普通混凝土在高速冲击荷载作用下会发生非常严重的脆性破坏,在结构的背面出现大面积的震塌、崩落以及碎片飞溅现象[3-6],对附近人员和设备的安全构成严重威胁。因此,找到合适的工程防护材料成为解决上述问题的关键。活性粉末混凝土(reactive powder concrete,RPC)[7]是根据最紧密堆积原理制备成的具有高强度和良好韧性的水泥基材料。由于其本身的高强度以及钢纤维和基体之间的桥连作用,使得RPC相较于传统混凝土能够有效地减小迎弹面的开坑尺寸、侵彻深度,降低结构发生整体性破坏的机率[8-9]。截至目前,已有很多学者针对RPC的抗侵彻性能开展了实验研究[10-11],但是关于RPC的实验数据仍然有限,无法系统地指导现有工程结构设计。进行大型抗侵彻实验存在实验成本高昂、周期长、动态数据采集困难等问题。

    近些年,随着有限元技术的不断发展,通过数值模拟方法预测混凝土类材料在高速冲击荷载作用下的破坏形态日益受到学者们的青睐。充分发挥有限元软件的优势,可以有效地指导实验、工程顺利进行。然而,目前对RPC材料进行动态数值模拟所选用的本构模型仍然是传统的混凝土本构模型,例如:HJC模型[12]、RHT模型[13]、CSCM模型[14]、K&C模型[15],这些模型都无法准确地反映RPC的拉伸延性特征以及应变率效应。因此,要想准确地描述RPC材料的基本力学特性,需要对上述模型的参数作出系统地调整。在以往的RPC高速冲击数值模拟研究中,运用较多的是K&C本构模型,该模型综合考虑了混凝土类材料的应变率效应、应变硬化软化现象、剪胀效应、围压效应,能够准确地描述混凝土类材料在大变形、高应变率和高静水压力下的力学特性。

    本文中,首先对RPC靶体进行多次侵彻实验,得到RPC靶体破坏数据,并计算出RPC的Forrestal[16-17]公式中与靶体材料相关的系数;然后基于K&C模型和现有RPC的基本力学性能实验数据,对K&C本构模型的强度面参数、损伤参数、损伤演化模型、应变率效应、状态方程参数进行修正,采用修正的K&C模型参数模拟上述多次侵彻实验并与实验结果进行对比,验证模拟结果的正确性;最后在质量为10 kg、直径为80 mm弹体正侵彻情况下,对长2 200 mm、宽2 200 mm、高1 260 mm的RPC靶体进行数值预测,得到弹体速度与侵彻深度的关系,并利用Forrestal公式验证数值预测结果的合理性,得到实验过程中弹体的最大加速度和弹体贯穿靶体时的极限速度,以期为实验的顺利进行提供指导。

    实验所用原料为硅酸盐水泥(PC),二氧化硅质量分数大于95%、比表面积为18 200 m2/kg的硅灰(SF),减水率大于40%的聚羧酸高效减水剂,比表面积大于10 000 m2/kg的矿渣,细度模数为2.5的普通砂,长度为13 mm、直径为0.175 mm的钢纤维,普通自来水。将以上原料按表1的配合比配制钢纤维掺量为2%的RPC。RPC靶体的浇筑采用强制式卧轴搅拌机搅拌,先将混合好的水泥、粉煤灰、硅灰、精细砂与减水剂干拌2 min,随后加入水搅拌2 min,然后把钢纤维均匀加入继续搅拌5 min,并浇筑到预先制作好的钢桶当中。同时也浇筑了基本力学性能测试试件,常温养护28 d后的测试结果如表2所示,其中fcft分别为材料的无侧限抗压强度和无侧限抗拉强度,E为弹性模量,μ为泊松比,ρ为密度。

    表  1  RPC材料配合比
    Table  1.  Mixture proportions of RPC
    kg/m3
    材料胶凝材料减水剂钢纤维
    RPC1 23892817.7234.3159
    下载: 导出CSV 
    | 显示表格
    表  2  RPC基本力学性能参数
    Table  2.  Basic mechanical performance parameters of RPC
    材料fc/MPaft/MPaE/GPaμρ/(g·cm−3
    RPC1209.2746.20.222.44
    下载: 导出CSV 
    | 显示表格

    实验所用弹体材质为30GrMnSiNi2A低合金超高强度钢,淬火后抗拉屈服强度可以达到1 650 MPa,由25 mm口径滑膛炮发射,在580 m/s弹速下弹体为刚性侵彻,弹体外形与尺寸如图1所示。弹体直径为25 mm,头部曲径比rCRH=3,长径比为6,壁厚与弹径比为0.14,底托采用尼龙材质。弹体内部填入惰性材料以调整质量,平均发射质量为352.9 g。RPC靶体为圆柱体,尺寸为600 mm×600 mm,侧面用4 mm厚钢板箍紧,直径为20倍弹径,可以忽略边界效应。

    图  1  弹体尺寸
    Figure  1.  Projectile sizes

    实验布置如图2所示。发射装置为25 mm口径滑膛炮,靶体安放在支架上,调整木楔让靶体表面与炮口轴线垂直,炮口距离靶体3 m。由滑膛炮分别发射3发弹丸对RPC靶体进行多次侵彻实验。着靶速度由测速板与电子计时仪测出。高速摄影机记录弹体着靶姿态,拍摄弹体飞行和弹靶初始作用过程。

    图  2  侵彻实验设备布置
    Figure  2.  Arrangement of penetration experiment equipments

    对RPC靶体分别进行3次打击,打击位置均为靶体正中心,弹体着靶时的速度分别为511.5、552.5、560.0 m/s,图3给出了高速摄影机记录的511.5 m/s速度下典型的侵彻过程。

    图  3  高速摄影机记录的弹体撞击靶体过程
    Figure  3.  Process of the projectile impacting the targetrecorded by a high-speed camera

    侵彻的实验结果如表3所示,其中v0为弹体的初速度,h为弹体侵入靶体沿轴线方向的深度,SH分别为迎弹面弹坑的面积和深度,N为迎弹面裂纹条数,Wmax为最大裂纹宽度。

    表  3  靶体多次侵彻实验结果
    Table  3.  Experimental results of targets subjected to multiple penetrations
    侵彻次数v0/(m·s−1h/mmS/cm2H/mmNWmax/mm
    1511.5129.1329.7 59.5 0 0
    2552.5257.4344.1 79.812 2
    3560.0290.3354.8114.91313
    下载: 导出CSV 
    | 显示表格

    表3可以看出:随着侵彻次数的增加,弹体的侵彻深度、迎弹面的弹坑深度、裂纹条数、最大裂纹宽度增加明显,但迎弹面弹坑面积基本没有变化。此外,为了准确预测混凝土类材料受弹体侵彻的侵彻深度,目前已有几十种经验公式可供参考,其中Forrestal公式较符合本次实验条件。选用第1次侵彻的实验结果来确定Forrestal公式中与靶体材料相关的系数,本实验侵彻深度大于2倍弹径,符合Forrestal公式的适用条件,侵彻深度L的计算公式为:

    L=2Mπd2ρBln(1+ρBv21Sfc)+2d
    (1)

    式中:M为弹体质量(kg),d为弹体直径(m),S为实验确定的常数,参数Bv21由下式给出:

    B=8rCRH124r2CRH,v21=2Mv20πd3Sfc2M+πd3Bρ
    (2)

    联立式(1)~(2)可得:

    S=2ρBMv20fc(2M+πd3Bρ)exp[(L2d)πd2ρB/(2M)+ρBπd3/(2M+πd3Bρ)1]
    (3)

    为确定RPC Forrestal侵彻深度计算公式中的抗侵彻系数,将弹体第1次侵彻靶体的实验结果与Wu等[11]、崔亚男[18]、Feng等[19]的实验结果代入式(3),分别得到不同强度的RPC Forrestal公式中的抗侵彻系数,如图4所示。可以看出,根据本次实验计算得到的Forrestal公式中的抗侵彻系数与Feng等、Wu等、崔亚男的实验结果较接近,可以认为本次实验结果有效,因此,本文中将120 MPa RPC的Forrestal公式中的抗侵彻系数定为S*=7.454 3。

    图  4  2%钢纤维掺量下不同强度RPC的Forrestal公式中的S*
    Figure  4.  S* in the Forrestal formula of RPC with differentcompressive strengths at 2% steel fiber content

    图5(a)中可以看到,弹体第1次侵彻RPC靶体后,靶体迎弹面未出现任何裂缝,只是形成了一个较小的弹坑,并且实验过程中发现,第1发弹体侵入靶体一定深度后弹出,这主要是由于钢纤维与RPC基体之间的桥联作用以及RPC本身强度较高,加上弹速较低,导致弹体未能直接侵入到靶体当中。弹体第2次侵入靶体后,迎弹面出现了许多微小的裂缝,但是弹坑面积未出现明显的增加,弹体直接侵入到靶体中,未出现弹出现象,如图5(b)所示。弹体第3次侵入后,迎弹面可见的裂缝数量相较于第2次侵彻没有明显增多,但是裂缝的宽度都有一定程度的增加,出现了一条贯穿的径向裂缝,弹坑面积基本没有发生变化,如图5(c)所示。

    图  5  RPC多次侵彻实验结果
    Figure  5.  Experimental results of RPC multiple penetrations

    为了观察弹体在靶体中的运动轨迹并回收弹体,对实验后的靶体进行了切割解剖,剖面如图6所示。由于第1发弹体在侵入靶体一定深度后直接弹出靶体,图6仅标出了第2发弹体和第3发弹体的运动轨迹,图6中的红色虚线为第2发和第3发弹体侵彻的弹道。可见,由于靶体中存在初始缺陷以及受到第1次侵彻的影响,第2次和第3次侵彻的弹道都发生了明显的偏转。

    图  6  多次侵彻实验后的靶体剖面
    Figure  6.  Profile of targets after multiple penetration experiments

    根据实验情况建立如图7所示的有限元模型,弹体、靶体、钢箍均采用SOLID164实体单元,并选用全尺寸模型进行计算。弹体与靶体之间设为面面侵蚀接触,靶体与钢箍设为面面自动接触。为了避免计算过程中单元发生畸变,引入MAT_ADD EROSION侵蚀准则,选用最大主应变控制单元失效[20],对于本次模拟侵蚀应变的取值,选取第1次侵彻模拟结果和实验数据吻合最佳时的侵蚀应变作为该组实验的侵蚀应变取值,即当单元的最大主应变超过0.15时,单元立即被删除。

    图  7  有限元模型
    Figure  7.  Finite element model

    材料模型选择和模型参数确定是决定模拟结果可靠的关键因素。弹体采用刚体材料模型,钢箍和钢筋均选用LS-DYNA中的P-K随动模型,弹体、钢箍以及钢筋的模型参数如表4所示,其中σy为屈服强度。

    表  4  弹体、钢箍以及钢筋材料模型参数
    Table  4.  Model parameters of projectile, steel culvert and steel bar material
    材料ρ/(kg·m−3E/GPaμσy/MPa
    弹体7 8502100.31 650
    钢箍/钢筋7 8002100.3 300
    下载: 导出CSV 
    | 显示表格

    活性粉末混凝土靶体的材料模型选用K&C模型,但K&C模型的参数是基于45.4 MPa混凝土实验数据由程序自动生成的,无法准确描述RPC的拉伸韧性、应变率关系和压力相关性,因此本文中将基于现有的RPC实验数据,通过调整K&C模型的参数,使之可以更好地描述RPC的基本力学特性。

    3.2.1   强度面参数

    K&C[15]模型引入了3个失效面:初始屈服面、极限强度面、残余强度面,分别描述混凝土的初始强度、极限强度及残余强度的变化规律,其表达式为:

    Δσy={a0y+p/(a1y+a2yp)p0.15fc1.35ft+3p(13ft/fc)0p0.15fc1.35(p+ft)p0
    (4)
    Δσm={a0+p/(a1+a2p)pfc/3(1.5/φ)(p+ft)0pfc/33(p/η+ft)p0
    (5)
    Δσr=a0f+p/(a1f+a2fp)
    (6)

    式中:aiaiyaifi=0,1,2)为由材料三轴围压实验数据确定的常数,p为静水压力,η(λ)为损伤演化关系,φ为拉压子午比。φ具体定义为:

    φ(p)={0.5p00.5+1.5ft/fcp=fc/31.15fc/[a0+2.3/(3a1+2.3a2fc)]p=2.3fc/30.753p=3fc1p8.45fc
    (7)

    当前失效面定义为:

    Δσ=rfΔσ(p/rf)=3J2={[η(ΔσmΔσy)+Δσy]rHardening[η(ΔσmΔσr)+Δσr]rSofting
    (8)

    式中:rf为应变率效应系数(ψ),r为当前子午与压缩子午的比值。

    在K&C模型中,强度面的9个参数均参照45.4 MPa混凝土的三轴围压实验数据由LS-DYNA程序自动生成,它们与材料抗压强度的关系表示为:

    {a0=0.2956fc,a1=0.4463,a2=0.0808/fca0y=0.2232fc,a1y=0.6250,a2y=0.2575/fca0f=0,a1f=0.4417,a2f=0.1183/fc
    (9)

    搜集了近些年与RPC和混凝土相关的三轴围压实验数据[21-31],如图8所示,发现在低归一化静水压力下,混凝土与RPC归一化等效应力值基本一致。由于RPC高归一化静水压力的实验数据很难获得,因此本文中假定RPC在高归一化静水压力下归一化等效应力值与混凝土一致,以此获得RPC的极限强度面参数。从图8中还可以看出,采用K&C自动生成的强度面参数明显低于实验数据点,而修正后的RPC强度面参数与实验数据非常吻合。而RPC弹性强度面参数的确定依然按照K&C模型中的方法,即假定Δσy=0.45Δσm[15]。Kong等[32]在修改的K&C模型中认为混凝土类材料的残余强度面平行于极限强度面,并应用于混凝土靶体的贯穿数值模拟当中,取得了不错的模拟效果。本文中也采取相同的方法,这样就可以全面地确定RPC的强度面的参数,具体表示为:

    图  8  混凝土和RPC三轴围压实验数据
    Figure  8.  Triaxial confining pressure experimental data of concrete and RPC
    {a0=0.4426fc,a1=0.5698,a2=0.02516/fca0y=0.2797fc,a1y=0.8989,a2y=0.06850/fca0f=0,a1f=0.5698,a2f=0.02516/fc
    (10)
    3.2.2   损伤演化模型

    K&C模型的当前失效面将随着损伤因子λ的增大先强化后软化,λ为等效塑性应变的函数,定义为:

    λ={t0d¯εp/[rf(1+p/ft)b1]p0t0d¯εp/[rf(1+p/ft)b2]p0
    (11)

    式中:d¯εp=23dεpijdεpij为有效塑性应变增量,dεpij为有效塑性应变增量张量,b1b2分别为压缩和拉伸损伤参数。

    当混凝土类材料三向受拉时,只产生体积应变而不产生形状应变,此时偏应变为零,因此损伤变量也为零。为考虑应力路径接近三向受拉时体积损伤应变的影响,K&C模型中引入了体积损伤变量Δλ

    Δλ=b3fdkd(εVεV,y)
    (12)

    式中:b3为控制混凝土三向受拉时软化段的损伤参数;kd为内变量因子;εV为体积应变,εV, y为屈服点附近的体积应变;fd则表征应力路径与三向受拉状态靠近的程度。fd由下式给出:

    fd={1|3J2/p|/0.10|3J2/p|0.10|3J2/p|0.1
    (13)

    在K&C模型中,损伤演化模型由13对(λ, η)数据点来定义,这些相邻的数据点由插值函数来连接。其定义是针对单一强度(45.4 MPa)混凝土实验得到的数据,与之吻合较好,但无法应用到其他强度的RPC中,更不能描述RPC材料具有的拉伸延性特征,因此需要提出一个新的、可靠的损伤演化关系来模拟RPC损伤关系。通过不断地对单个单元进行拉伸、压缩模拟,最终确定损伤参数b1b2分别为2.2和–7.25,将反映ηλ关系的13对点定义如表5所示。

    表  5  损伤演化函数η(λ)
    Table  5.  Damage evolution function η(λ)
    λη λη
    00 4.0×10−60.51
    2.7×10−50.626.7×10−40.37
    6.8×10−50.921.2×10−30.27
    8.0×10−50.992.0×10−30.20
    1.0×10−41.005.5×10−30.10
    1.4×10−40.961.6×10−20
    2.6×10−40.66
    下载: 导出CSV 
    | 显示表格

    图9展示了采用自动生成参数的K&C模型、新参数的K&C模型以及RPC单轴拉伸压缩实验数据,结果表明,新的K&C参数模拟的结果与实验数据更吻合。

    图  9  拉伸和压缩应力应变曲线
    Figure  9.  Tensile and compressive stress-strain curves

    图10给出了体积损伤参数b3对三轴拉伸应力应变曲线的影响,可以看出,三轴拉伸应力应变曲线的软化梯度随b3的减小而减小。对三轴拉伸体积损伤参数b3的取值,Weerhijm等[33]假定:单轴拉伸的断裂应变和三轴拉伸的断裂应变相等。因此本文中b3取值为0.018。

    图  10  不同b3值对应的三轴拉伸应力应变曲线
    Figure  10.  Triaxial tension strain stress curvescorresponding to different b3 values
    3.2.3   应变率参数

    材料的动态力学特性会随着加载速率的变化而变化。学者们针对普通混凝土的这种特性已经开展了深入的研究,发现在高加载速率下,混凝土的压缩强度会增加100%,拉伸强度会增加600%[34],并且提出了可靠的应变率公式加以描述,如K&C模型中的应变率公式被广泛应用于混凝土受到冲击、爆炸荷载作用下的数值模拟中。RPC材料的拉伸实验表明[35],RPC材料拥有较普通混凝土更明显的率依赖性。因此,在RPC材料动态数值模拟中考虑应变率效应就显得非常重要。此前,许多学者[36-38]对RPC受冲击和爆炸荷载的数值模拟都采用基于普通混凝土提出的应变率公式,这显然无法较准确地模拟RPC在高应变率荷载下的动态力学响应。Park等[38]对基体强度分别为56、81和180 MPa的RPC进行了拉伸应变率敏感性试验研究,并提出了拉伸应变率公式(式(14))。Lin[39]基于现有RPC的动态压缩实验数据,提出了动态压缩应变率公式(式(15)),并将其应用于RPC板抗爆数值模拟中,取得了不错的模拟效果。本文中也采用相同的计算公式(式(14)~(15)),并据此计算出RPC的应变率效应数据点,如表6所示。此外,最新版本的K&C模型中将混凝土动态拉伸增强因子和动态压缩增强因子上限分别设置为9.7和2.94,以避免高估混凝土材料在高应变率下的应变率效应。本文中也将这一设置引入到RPC的动态数值模拟:

    表  6  活性粉末混凝土的K&C模型应变率效应特征点取值
    Table  6.  K&C model strain rate characteristic points of reactive powder concrete
    ˙ε/s1ψ˙ε/s1ψ˙ε/s1ψ
    –30 0009.97–101.27–1×10–41.07
    –4 7829.97–31.24–1×10–51.03
    –1 0005.41–11.2201.00
    –3003.45–0.11.18301.00
    –1002.29–0.011.142652.94
    –251.28–1×10–31.1130 0002.94
    下载: 导出CSV 
    | 显示表格
    ψt={(˙ε/˙ε0)0.01465˙ε25s10.002352(˙ε/˙ε0)0.3735˙ε25s1
    (14)
    ψc={1˙ε30s10.728+0.008337(˙ε/˙ε0)˙ε30s1
    (15)

    式中:ψtψc分别为动态拉伸增强因子和动态压缩增强因子,˙ε0=106s1为静态应变率。

    3.2.4   状态方程参数

    K&C模型采用8号状态方程来定义混凝土压力与体积应变之间的关系:

    p=C(εV)+γ0T(εV)E0
    (16)

    式中:E0为初始体积内能,γ0为温度特征系数,C(εV)和T(εV)分别为压力和温度与体积应变之间的关系。在加载阶段,压力由式(4)确定,卸载沿卸载体积模量至压力截断点(图11(a)),重新加载首先沿卸载路径至卸载开始点,然后沿加载路径继续上升。

    图  11  状态方程
    Figure  11.  State equation

    尽管K&C模型中的8号状态方程能够很好地描述混凝土压力与体积应变之间的关系,但是由于RPC的孔隙率远小于普通混凝土,导致原来的8号状态方程无法准确地描述RPC的压力与体积应变的关系。材料的状态方程所需的实验数据通常需要进行飞片撞击试验。因此本文中首先通过基本力学性能实验测试求得RPC的弹性模量和泊松比,并按照式(17)求得RPC材料的体积模量KV,然后基于Marsh[40]、高乐[41]和严少华等[42]的Hugoniot试验数据校核K&C模型中的状态方程,以获得适合于RPC的8号状态方程参数:

    KV=E3(12μ)
    (17)

    图11(b)展示了校核结果,可以看出,自动生成的8号状态方程明显低于试验数据,而校核后的8号状态方程与试验数据基本一致。校核后的状态方程参数如表7所示,其中εV为体积应变,σV为相应体积应变所对应的体积应力,Kav代表卸载模量。其他参数对计算结果基本没有影响,因此选用K&C模型自动生成的参数,至此关于RPC的K&C模型参数完全确定。

    表  7  RPC的K&C模型8号状态方程参数
    Table  7.  Parameters of No. 8 equation of state in the K&C model of RPC
    εV
    εV1εV2εV3εV4εV5εV6εV7εV8εV9εV10
    00.00150.00430.01010.03050.05130.07260.09430.1740.208
    σV/GPa
    σV1σV2σV3σV4σV5σV6σV7σV8σV9σV10
    00.0410.0940.2920.8811.6222.5113.5738.71411.579
    Kav/GPa
    Kav1Kav2Kav3Kav4Kav5Kav6Kav7Kav8Kav9Kav10
    27.527.527.88529.28834.84340.42545.98050.186112.915137.5
    下载: 导出CSV 
    | 显示表格

    为进一步观察靶体内部损伤状况并验证所选K&C模型参数的准确性,利用LS-DYNA中的重启动方法分别对3次侵彻进行数值模拟,该方法可以记录第1次弹体冲击RPC靶体后的应力、应变和损伤状态并输出重启动文件。在分析弹体第2次冲击RPC靶体时,LS-DYNA会读取弹体第1次冲击RPC靶体后的应力、应变和损伤状态等数据并继续进行计算,这样便可以实现对弹体多次冲击受损后RPC靶体的数值模拟,模拟结果如图1214所示。利用数值模拟得到弹体第1次的侵彻深度和开坑面积分别为13.2 cm、315.2 cm2,而实验结果分别为12.9 cm、329.7 cm2。第2次侵彻深度和开坑面积的模拟结果分别为25.9 cm、336.4 cm2,而实验结果分别为25.7 cm、344.1 cm2。可知前两次模拟得到的侵彻深度与实验结果基本一致,误差不超过5%。图14给出了第3次侵彻的模拟结果,第3次模拟得到的侵彻深度为34.0 cm,与实验得到的侵彻深度(29.0 cm)存在较大误差。并且模拟得到的弹体偏转角度小于实验结果,这样就使得模拟得到的侵彻深度大于实验测量的结果。这主要是由于在进行有限元分析时将RPC靶体视为完全均质材料,因此在第1次和第2次模拟结果中可以看到,RPC靶体内部的损伤也是均匀分布的,弹体侵彻RPC靶体的弹道没有明显偏转。然而实验中RPC靶体并非完全均质,一部分区域存在初始缺陷,造成弹体第1次冲击后靶体内部的损伤更加不均匀。当弹体第2次冲击靶体时,弹体四周受到的阻力不一致,导致弹体沿损伤更严重的区域偏转。因此才会看到实验中第2发弹体和第3发弹体都发生了明显的偏转,而模拟结果无法再现实验中弹体发生的明显偏转。

    图  12  第1次侵彻RPC靶体的模拟结果
    Figure  12.  Simulation results of the first penetration of the RPC target
    图  13  第2次侵彻RPC靶体的模拟结果
    Figure  13.  Simulation results of the second penetration of the RPC target
    图  14  第3次侵彻RPC靶体的模拟结果
    Figure  14.  Simulation results of the third penetration of the RPC target

    此外,K&C模型中将损伤定义为D=2λ/(λ+λm),当D从0到1时,表示材料进入了强化阶段,当D从1到2时,表示材料进入了软化阶段,而当D=2时认为单元完全失效。从损伤云图1214中可以看到,随着侵彻次数增多,靶体内部的损伤也不断扩大,同时迎弹面的裂缝数量也在增多。在迎弹面,D>1.8的损伤面积没有发生太大变化,这与多次侵彻RPC靶体后迎弹面开坑面积基本相同的实验结果一致。因此可以认为RPC的K&C模型参数和所用的模拟方法具有可靠性。

    基于上述RPC的K&C模型参数和模拟方法,采用LS-DYNA软件对RPC侵彻实验进行数值预测,为实验测试系统的选用和弹速设计提供指导。该实验的靶体尺寸为2 200 mm×2 200 mm ×1 260 mm,弹体直径为80 mm,长为361.35 mm,弹体质量为10 kg并为靶体配置钢筋,具体参数如图15所示。

    图  15  弹体与钢筋分布示意图(单位:mm)
    Figure  15.  Schematic diagram of the distribution of projectiles and reinforcement (unit: mm)

    根据实验情况的对称性,为了提高计算效率和节省计算资源,建立如图16所示的1/4有限元模型,弹体与RPC靶体均采用SOLID164实体单元,钢筋采用BEAM161梁单元,靶体与钢筋之间选用CONSTRAINED_LAGRANGE_IN_SOLID耦合方式进行相互作用,弹体与靶体之间采用CONTACT_ERODING_SURFACE_TO_SURFACE面面侵蚀接触关键字控制,模型单位制m-s-kg,共包含节点2 395 628个、六面体单元1 592 608个、梁单元7 267个。弹体和钢筋均选用P-K随动强化模型,靶体采用修正的K&C模型参数。

    图  16  有限元模型
    Figure  16.  Finite element model

    图17(a)(c)为卵形弹以850 m/s的速度正侵彻RPC靶的数值模拟结果。可以看到,弹体高速侵彻对RPC产生了明显的局部破坏,形成了直径为0.356 m、深度为0.192 m的弹坑,同时形成了直径略大于弹径的圆柱形孔道,迎弹面产生了数十条主裂缝,弹体侵彻靶体所造成的损伤主要集中在柱形孔道。图17(d)(e)分别显示了弹体的加速度、速度、位移时程曲线,可以看到弹体的最大加速度为63 000g,弹体的速度在侵入靶体后一直减小,直至零。

    图  17  850 m/s弹速下的侵彻结果
    Figure  17.  Penetration results at the projectile velocity of 850 m/s

    图18显示了弹体以1 150 m/s的速度贯穿RPC靶体的模拟结果,同850 m/s的侵彻结果一样,靶体形成了直径为0.340 m、深度为0.213 m的漏斗坑。不同的是迎弹面产生了更多裂缝,且靶体内部损伤更严重,此外背弹面出现了弹体贯穿形成了直径为0.505 m、深度为0.264 m的贯穿坑和宏观主裂纹。从图18(d)(e)中可以看出,弹体加速度在0.2 ms附近达到峰值,最大值为73 000g,相较于850 m/s的侵彻结果,加速度的峰值提高了16%,同时弹体在贯穿靶体后的剩余速度为100 m/s。

    图  18  1 150 m/s弹速下的贯穿结果
    Figure  18.  Penetration results at the projectile velocity of 1 150 m/s

    通过上述数值模拟方法,可以得到不同速度弹体侵入RPC靶体后的侵彻深度,不同弹速下弹体侵入RPC靶体的深度与RPC-Forrestal公式的预测结果如图19所示,可以发现,采用本文中的数值模拟方法预测的侵彻深度与RPC-Forrestal公式预测的结果较吻合。因此可以认为本次数值模拟的结果有效。另外,根据数值模拟和RPC-Forrestal公式可以发现:当弹体速度超过1 100 m/s时,弹体会贯穿RPC靶体。

    图  19  弹速与侵彻深度之间的关系
    Figure  19.  Relationship between projectile velocity and penetration depth

    开展了直径25 mm弹体多次侵彻活性粉末混凝土靶实验,得到了靶体的破坏数据;确定了活性粉末混凝土的K&C模型参数;通过多次侵彻实验验证了改进后K&C模型参数的正确性,并对直径为80 mm弹体冲击2 200 mm×2 200 mm×1 260 mm活性粉末混凝土靶体的实验结果进行数值预测,得到以下结论:

    (1)在相同弹速下沿同一点逐次侵彻RPC靶体,随着侵彻次数的增加,弹体的侵彻深度、迎弹面的弹坑深度、裂纹条数、最大裂纹宽度明显增大,但是迎弹面弹坑面积基本没有变化;

    (2)通过侵彻实验确定了抗压强度为120 MPa RPC的Forrestal侵彻深度计算公式中抗侵彻系数S*=7.454 3;

    (3)基于RPC的准静态单轴压缩、拉伸实验、三轴围压实验、拉伸压缩应变率实验和飞片撞击实验数据,提出了一组适用于RPC的K&C本构模型参数;

    (4)采用非线性有限元方法对侵彻实验进行了数值模拟,模拟结果与实验结果吻合较好,并采用同样的方法对实验进行了数值预测,得到了两种预设弹速下弹体的加速度、速度、位移时程曲线以及不同弹速下弹体侵入靶体的深度,并利用Forrestal经验公式验证了数值预测结果的准确性,最终确定了弹体贯穿靶体的速度至少要达到1 100 m/s。

  • 图  1  柱形装药爆炸试验传感器布置

    Figure  1.  Cylindrical charge explosion test setup and sensor layout

    图  2  柱形装药空中爆炸试验有限元模型

    Figure  2.  Finite element model of cylindrical charges air explosion test

    图  3  不同网格尺寸下的压力时程曲线

    Figure  3.  Pressure-time histories corresponding to different grid sizes

    图  4  入射和反射超压时程曲线的试验和数值模拟结果对比

    Figure  4.  Comparisons of the test and simulated incident and reflected overpressure-time histories

    图  5  Shi等[13]爆炸试验设置和有限元模型

    Figure  5.  Explosion test setup, finite element model of Shi et al.[13]

    图  6  典型测点处超压时程曲线

    Figure  6.  Overpressure-time histories at typical measurement points

    图  7  UFC 3-340-02规范[4]球形装药中心起爆有限元模型

    Figure  7.  Finite element model of spherical charges ignited at the central point in UFC 3-340-02 specification[4]

    图  8  入射峰值超压的规范[4]和模拟结果对比

    Figure  8.  Comparisons of simulated and specified[4] peak incident overpressure

    图  9  柱形装药长径比、起爆方式和方位角

    Figure  9.  Length-to-diameter ratio, initiation method, and azimuth angle of cylindrical charges

    图  10  柱形装药空中爆炸自由场有限元模型

    Figure  10.  Finite element model of cylindrical charges air explosion

    图  11  柱形装药轴向和径向冲击波反射示意图

    Figure  11.  Schematic diagram of axial and radial blast wave reflection of cylindrical charge

    图  12  典型柱形装药空中爆炸的反射场有限元模型

    Figure  12.  Typical finite element model for reflections of cylindrical charges air explosion

    图  13  柱形装药3种起爆方式下的压力云图(L/D=1)

    Figure  13.  Pressure contours of cylindrical charge under three ignition methods (L/D=1)

    图  14  典型长径比柱形装药中心起爆压力云图

    Figure  14.  Pressure contours of central ignited cylindrical charge with typical L/D

    图  15  柱形装药中心起爆不同时刻的压力云图(L/D=1)

    Figure  15.  Instantaneous pressure contours of central ignited cylindrical charge (L/D=1)

    图  16  典型长径比柱形装药3种起爆方式的峰值超压分布

    Figure  16.  Peak overpressure distributions of cylindrical charge with typical L/D and three ignition methods

    图  17  典型长径比柱形装药3种起爆方式的峰值超压形状因子分布

    Figure  17.  Peak overpressure shape factor distributions of cylindrical charge with typical L/D and three ignition methods

    图  18  典型峰值超压形状因子与最大峰值超压形状因子

    Figure  18.  Typical shape factor and maximum shape factor of peak overpressure

    图  19  3种起爆方式下不同长径比柱形装药的峰值超压临界比例距离

    Figure  19.  Critical scaled distance of peak overpressure for cylindrical charge with different L/D under three ignition methods

    图  20  典型长径比柱形装药3种起爆方式的最大冲量分布

    Figure  20.  Maximal impulse distributions of cylindrical charge with typical L/D and three ignition methods

    图  21  典型长径比柱形装药3种起爆方式的最大冲量形状因子分布

    Figure  21.  Maximal impulse shape factor distributions of cylindrical charge with typical L/D and three ignition methods

    图  22  典型最大冲量形状因子与最大的最大冲量形状因子

    Figure  22.  Typical shape factor and maximum shape factor of maximal impulse

    图  23  3种起爆方式下不同长径比柱形装药的最大冲量临界比例距离

    Figure  23.  Critical scaled distance of maximal impulse for cylindrical charge with different L/D under three ignition methods

    图  24  理想爆炸波超压时程曲线

    Figure  24.  Ideal blast wave overpressure-time histories

    图  25  3种起爆方式下柱形装药轴向反射峰值超压与投影距离的关系

    Figure  25.  Dependence of axial reflected peak overpressure on projection distance of cylindrical charge with three ignition methods

    图  26  典型长径比下柱形装药轴向反射峰值超压与投影距离的关系

    Figure  26.  Dependence of axial reflected peak overpressure on projection distance of cylindrical charge with typical L/Ds

    图  27  典型长径比柱形装药中心起爆的轴向和径向超压反射系数

    Figure  27.  Axial and radial overpressure reflection coefficients of central ignited cylindrical charge with typical L/Ds

    图  28  L/D=0.25柱形装药单/双端起爆的轴向和径向峰值超压反射系数

    Figure  28.  Axial and radial peak overpressure reflection coefficients of single/double-end ignited cylindrical charge with L/D=0.25

    图  29  爆炸荷载试验与计算结果对比

    Figure  29.  Comparisons of test data and calculated results of blast loadings

    表  1  L/D1柱形装药入射峰值超压拟合公式系数

    Table  1.   Fitted formula coefficients for peak incident overpressure of central ignited cylindrical charges with L/D1

    起爆方式 α/(°) 0.3 m/kg1/3Z≤1 m/kg1/3 1 m/kg1/3<Z≤15 m/kg1/3
    A B C R2 A B C R2
    中心起爆 0 173 9 399 0.943 7 −211 390 31 0.612 6
    30 −33 503 −166 0.904 5 −200 424 40 0.902 7
    60 −103 958 −491 0.988 3 −202 644 −11 0.969 0
    90 −242 1 352 −141 0.926 4 281 501 −11 0.986 3
    单端起爆 0 −56 1 013 −900 0.884 5 −491 619 22 0.616 4
    30 −39 365 −96 0.800 2 −284 546 1 0.880 1
    60 60 −131 693 0.940 7 −114 509 9 0.961 8
    90 −225 1 438 −244 0.941 8 306 512 −19 0.976 8
    120 −205 1 732 −1 360 0.977 9 −238 721 −16 0.969 6
    150 −27 580 −261 0.943 1 −171 431 43 0.915 8
    180 77 1 246 −1 057 0.909 9 −531 574 34 0.775 3
    双端起爆 0 −25 865 −796 0.875 3 −887 1 096 −56 0.733 9
    30 −29 304 −24 0.793 5 −313 593 4 0.914 4
    60 29 22 602 0.954 3 −69 505 18 0.968 2
    90 −255 2 054 −1 068 0.923 7 187 694 −74 0.932 8
    下载: 导出CSV

    表  2  最大冲量临界比例距离拟合公式系数

    Table  2.   Coefficients of fitted formula for the maximal impulse critical scaled distance

    起爆方式 拟合公式系数
    k a b c R2
    中心起爆 2.13 0.25 1.14 0.32 0.945 0
    单端起爆 2.48 0.24 1.09 0.61 0.981 2
    双端起爆 2.60 0.23 0.97 0.46 0.927 8
    下载: 导出CSV

    表  3  最大冲量拟合公式系数

    Table  3.   Coefficients of the fitted formula for the maximal impulse

    起爆方式 公式拟合系数
    A1,1 A1,2 A1,1,1 A1,1,2 A1,2,1 A1,2,2 A1,3,1 A1,3,2 A1,4,1 A1,4,2
    中心 2 088 −2 044 1 222 −1 043 −4 422 4 503 2 583 −2 715 −971 1 101
    单端 341 −110 276 −257 186 −163 247 −236 −27 105
    双端 640 −384 1 304 −783 −2077 1 352 1 026 −620 85 −190
    起爆方式 A2,1 A2,2 A2,1,1 A2,1,2 A2,2,1 A2,2,2 A2,3,1 A2,3,2 A2,4,1 A2,4,2
    中心 −543 568 −550 459 1 321 −1 436 −677 824 173 −293
    单端 −135 21 −127 101 −43 22 −106 95 −40 −23
    双端 −88 −123 −516 225 528 16 −154 −126 −174 250
    起爆方式 A3,1 A3,2 A3,1,1 A3,1,2 A3,2,1 A3,2,2 A3,3,1 A3,3,2 A3,4,1 A3,4,2
    中心 45 −50 60 −49 −110 134 42 −73 3 21
    单端 18 −1 17 −12 7 −2 13 −11 11 2
    双端 15 28 57 −16 −63 −33 12 34 25 −38
    下载: 导出CSV
  • [1] US Department of the Army. Fundamentals of protective design for conventional weapons: TM 5-855-1 [S]. Washington, USA: US Department of the Army, 1986.
    [2] American Society of Civil Engineers. Blast protection of buildings: ASCE 59-11 [S]. Reston, Virginia, USA: American Society of Civil Engineers, 2011.
    [3] Canadian Standards Association. Design and assessment of buildings subjected to blast loads: CSA/S850-23 [S]. Toronto, Canda: Canadian Standards Association, 2023.
    [4] US Department of Defense. Structures to resist the effects of accidental explosions, with change 2: UFC 3-340-02 [S]. Washington, USA: US Department of Defense, 2014.
    [5] 中华人民共和国国家质量监督检验检疫总局, 中国国家标准化管理委员会. 爆破安全规程: GB 6722-2014 [S]. 北京: 中国标准出版社, 2014.
    [6] STONER R G, BLEAKNEY W. The attenuation of spherical shock waves in air [J]. Journal of Applied Physics, 1948, 19(7): 670–678. DOI: 10.1063/1.1698189.
    [7] BRODE H L. Numerical solutions of spherical blast waves [J]. Journal of Applied Physics, 1955, 26(6): 766–775. DOI: 10.1063/1.1722085.
    [8] BAKER W E. Explosions in air [M]. Austin, USA: University of Texas Press, 1974: 6–10.
    [9] HENRYCH J, ABRAHAMSON G R. The dynamics of explosion and its use [M]. Amsterdam New York, USA: Elsevier Science Publishing Company, 1979: 218.
    [10] MILLS C A. The design of concrete structures to resist explosions and weapon effects [C]// The 1st International Conference on Concrete for Hazard Protections. Edinburgh, UK: European Cement Association, 1987: 11–15.
    [11] WU C Q, HAO H. Modeling of simultaneous ground shock and airblast pressure on nearby structures from surface explosions [J]. International Journal of Impact Engineering, 2005, 31(6): 699–717. DOI: 10.1016/j.ijimpeng.2004.03.002.
    [12] KINNEY G F, GRAHAM K J, KENNETH J. Explosive shocks in air [M]. Berlin, Germany: Springer Verlag, 1985: 1–17.
    [13] SHI Y C, WANG N, CUI J, et al. Experimental and numerical investigation of charge shape effect on blast load induced by near-field explosions [J]. Process Safety and Environmental Protection, 2022, 165: 266–277. DOI: 10.1016/j.psep.2022.07.018.
    [14] ISMAIL M M, MURRAY S G. Study of the blast waves from the explosion of nonspherical charges [J]. Propellants, Explosives, Pyrotechnics, 1993, 18: 132–138. DOI: 10.1002/prep.19930180304.
    [15] SIMOENS B, LEFEBVRE M H, MINAMI F. Influence of different parameters on the TNT-equivalent of an explosion [J]. Central European Journal of Energetic Materials, 2011, 8(1): 53–67.
    [16] ANASTACIO A C, KNOCK C. Radial blast prediction for high explosive cylinders initiated at both ends [J]. Propellants, Explosives, Pyrotechnics, 2016, 41(4): 682–687. DOI: 10.1002/prep.201500302.
    [17] KNOCK C, DAVIES N. Predicting the peak pressure from the curved surface of detonating cylindrical charges [J]. Propellants, Explosives, Pyrotechnics, 2011, 36(3): 203–209. DOI: 10.1002/prep.201000001.
    [18] KNOCK C, DAVIES N. Predicting the impulse from the curved surface of detonating cylindrical charges [J]. Propellants, Explosives, Pyrotechnics, 2011, 36(2): 105–109. DOI: 10.1002/prep.201000002.
    [19] KNOCK C, DAVIES N, REEVES T. Predicting blast waves from the axial direction of a cylindrical charge [J]. Propellants, Explosives, Pyrotechnics, 2015, 40(2): 169–179. DOI: 10.1002/prep.201300188.
    [20] GAO C, KONG X Z, FANG Q, et al. Numerical investigation on free air blast loads generated from center-initiated cylindrical charges with varied aspect ratio in arbitrary orientation [J]. Defence Technology, 2022, 18(9): 1662–1678. DOI: 10.1016/j.dt.2021.07.013.
    [21] WU C Q, FATTORI G, WHITTAKER A, et al. Investigation of air-blast effects from spherical-and cylindrical-shaped charges [J]. International Journal of Protective Structures, 2010, 1(3): 345–362. DOI: 10.1260/2041-4196.1.3.345.
    [22] HU Y, CHEN L, FANG Q, et al. Blast loading model of the RC column under close-in explosion induced by the double-end-initiation explosive cylinder [J]. Engineering Structures, 2018, 175: 304–321. DOI: 10.1016/j.engstruct.2018.08.013.
    [23] SHERKAR P, SHIN J, WHITTAKER A, et al. Influence of charge shape and point of detonation on blast-resistant design [J]. Journal of Structural Engineering, 2016, 142(2): 1–11. DOI: 10.1061/(asce)st.1943-541x.0001371.
    [24] XIAO W F, ANDRAE M, GEBBEKEN N. Effect of charge shape and initiation configuration of explosive cylinders detonating in free air on blast-resistant design [J]. Journal of Structural Engineering, 2020, 146(8): 1–13. DOI: 10.1061/(asce)st.1943-541x.0002694.
    [25] THAM C Y. Numerical simulation on the interaction of blast waves with a series of aluminum cylinders at near-field [J]. International Journal of Impact Engineering, 2009, 36(1): 122–131. DOI: 10.1016/j.ijimpeng.2007.12.011.
    [26] SIMOENS B, LEFEBVRE M. Influence of the shape of an explosive charge: quantification of the modification of the pressure field [J]. Central European Journal of Energetic Materials, 2015, 12(2): 195–213.
    [27] PAPE R, MNISZEWSKI K R, LONGINOW A, et al. Explosion phenomena and effects of explosions on structures Ⅲ: methods of analysis (explosion damage to structures) and example cases [J]. Practice Periodical on Structural Design and Construction, 2010, 15(2): 153–169. DOI: 10.1061/(ASCE)SC.1943-5576.0000040.
    [28] KNOCK C, DAVIES N. Blast waves from cylindrical charges [J]. Shock Waves, 2013, 23(4): 337–343. DOI: 10.1007/s00193-013-0438-7.
    [29] YANG T C, LUO Y Z, HU G Q, et al. Probability distribution and determination of blast loading during structural blast resistant study [J]. Shock and Vibration, 2022. DOI: Artn 736728810.1155/2022/7367288.
    [30] WISOTSKI J, SNYER. W H. Characteristics of blast waves obtained from cylindrical high explosive charges: 80210, DRI-2286 [R]. Denver, USA: University of Denver, Denver Research Institute, 1965.
    [31] PLOOSTER M N. Blast effects from cylindrical explosive charges: experimental measurements: NWC TP 6382 [R]. China Lake, USA: Naval Report Centre, 1982.
    [32] SU Q, WU H, SUN H S, et al. Experimental and numerical studies on dynamic behavior of reinforced UHPC panel under medium-range explosions [J]. International Journal of Impact Engineering, 2021, 148: 1–23. DOI: 10.1016/j.ijimpeng.2020.103761.
    [33] TIAN S Z, YAN Q S, DU X L, et al. Experimental and numerical studies on the dynamic response of precast concrete slabs under blast load [J]. Journal of Building Engineering, 2023, 70: 1–18. DOI: 10.1016/j.jobe.2023.106425.
    [34] WHARTON R K, FORMBY S A, MERRIFIELD R. Airblast TNT equivalence for a range of commercial blasting explosives [J]. Journal of Hazardous Materials, 2000, 79(1): 31–39. DOI: 10.1016/S0304-3894(00)00168-0.
  • 期刊类型引用(4)

    1. 康耕新,颜海春,张亚栋,刘明君,郝礼楷. 接触爆炸下混凝土墩破坏效应试验与数值模拟. 兵工学报. 2024(01): 144-155 . 百度学术
    2. 常慧珠,张冬梅,麻旭东,李世中,贾保峰,赵紫良,何隆. 弹丸侵彻混凝土目标的热固耦合数值分析. 中北大学学报(自然科学版). 2024(02): 185-193 . 百度学术
    3. 崔冰,王景全,刘加平. UHPC桥梁研究进展与规模化应用技术路径分析. 中国公路学报. 2023(09): 1-19 . 百度学术
    4. 钟锐,张峰领. 超高性能混凝土、纤维增强高强及高延性混凝土抗侵彻性能比较. 硅酸盐学报. 2021(11): 2423-2434 . 百度学术

    其他类型引用(9)

  • 加载中
图(29) / 表(3)
计量
  • 文章访问数:  561
  • HTML全文浏览量:  200
  • PDF下载量:  288
  • 被引次数: 13
出版历程
  • 收稿日期:  2023-05-25
  • 修回日期:  2023-12-27
  • 网络出版日期:  2024-01-23
  • 刊出日期:  2024-04-07

目录

/

返回文章
返回