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

特定热层温度下入射角对前驱波特性的影响

李琦 程帅 刘文祥 金龙 童念雪 张德志

李琦, 程帅, 刘文祥, 金龙, 童念雪, 张德志. 特定热层温度下入射角对前驱波特性的影响[J]. 爆炸与冲击, 2024, 44(7): 073202. doi: 10.11883/bzycj-2023-0114
引用本文: 李琦, 程帅, 刘文祥, 金龙, 童念雪, 张德志. 特定热层温度下入射角对前驱波特性的影响[J]. 爆炸与冲击, 2024, 44(7): 073202. doi: 10.11883/bzycj-2023-0114
LI Qi, CHENG Shuai, LIU Wenxiang, JIN Long, TONG Nianxue, ZHANG Dezhi. Influence of incident angle on precursor wave characteristics at specific thermal-layer temperature[J]. Explosion And Shock Waves, 2024, 44(7): 073202. doi: 10.11883/bzycj-2023-0114
Citation: LI Qi, CHENG Shuai, LIU Wenxiang, JIN Long, TONG Nianxue, ZHANG Dezhi. Influence of incident angle on precursor wave characteristics at specific thermal-layer temperature[J]. Explosion And Shock Waves, 2024, 44(7): 073202. doi: 10.11883/bzycj-2023-0114

特定热层温度下入射角对前驱波特性的影响

doi: 10.11883/bzycj-2023-0114
详细信息
    作者简介:

    李 琦(1997- ),男,学士,助理工程师,liqi@nint.ac.cn

    通讯作者:

    张德志(1973- ),男,博士,研究员,zhangdezhi@nint.ac.cn

  • 中图分类号: O383.3

Influence of incident angle on precursor wave characteristics at specific thermal-layer temperature

  • 摘要: 空中强爆炸会释放热辐射使地表形成热层,冲击波进入热层后传播速度加快,形成前驱波。冲击波入射角是影响前驱波特性的重要因素,但目前相关工作大多依赖理论推导,实验研究较少。利用爆炸波模拟激波管平台,开展了热层温度为300 ℃时入射角对前驱波形成影响的研究实验,结合实验构型建立了数值仿真模型,并将实验、数值仿真结果与已有理论结果进行了对比。结果表明:实验获得的前驱波形成临界角范围与理论计算结果一致;入射角越大,前驱波超过马赫杆的距离越大,到时提前越多;前驱波会导致超压峰值减小,且随着入射角的增大,超压峰值减小程度先增大后减小;整体上看,前驱波动压峰值随入射角的增大而增大,当入射角达到一定阈值后,动压峰值增大程度开始在一定范围内波动,这是由气流密度和粒子速度的峰值到时不同所导致的;动压冲量的增大程度随入射角的增大逐渐增大。
  • 空中强爆炸会释放热辐射,在冲击波到达前,地面附近会形成热空气层,即热层。热层温度在400 K以上,最高可达2 000 K[1]。由于热层气体的温度、声速均高于大气环境中气体的温度、声速,强爆炸冲击波经过热层时可能形成热层前驱波[2],后文简称前驱波。前驱波可能导致地面附近的动压峰值、动压冲量大幅增大[3-6],使地面目标遭受更严重的毁伤。

    了解前驱波的形成机理、载荷历程特性是开展前驱波毁伤效应研究的前提。目前,关于前驱波的研究,特别是关于冲击波以一定角度入射热层这类问题的相关研究较少,已有的研究多以理论分析和数值模拟为主。Miller等[7]针对SMOKY实验进行了系列仿真计算,发现冲击波与实验场地的斜坡地形作用形成的马赫杆在热层的影响下发生明显前倾,该前倾特征可作为前驱波形成的证明。Zaslavskii等[8]通过理论计算研究了斜入射平面冲击波与水平热层的相互作用,发现当冲击波与水平方向的夹角小于某一临界值时不会出现前驱波,并进一步指出形成前驱波的必要条件是入射波波速的水平分量小于热层中冲击波的波速。乔登江[9]依据冲击波在不同性质气体界面的折射规律,分析得到形成前驱波的临界入射角的正弦等于热层内外声速的比。贾雷明等[10]通过仿真计算研究了前驱波特性,发现入射角为60°时,热层中靠近壁面处的压力曲线出现双峰结构,其中第1个峰值远小于无热层时的压力峰值,第2个峰值与无热层时的压力峰值相近,认为第1个峰是前驱波到达导致,第2个峰与波后的漩涡结构相关。关于前驱波特性的实验研究,大多针对入射角为90°时冲击波与水平热层相互作用的特殊工况。例如:Griffith[11]研究了激波管中入射波马赫数小于1.14时前驱波结构的演化过程,发现冲击波运动经过10倍热层厚度后,波阵面的形状趋于稳定。Gion[12]采用与文献[11]中相同的装置观测了更高热层温度的前驱波,发现前驱波与主激波之间存在一个过渡区,并利用纹影图像计算得到了过渡区内的温度分布。

    综上,目前对于冲击波以一定角度入射热层这类问题的研究较少,且研究方法以理论计算和数值仿真为主,实验研究仅涉及入射角为90°时冲击波与水平热层相互作用的特殊工况,主要关注前驱波的结构特征及产生条件,对于前驱波的超压、动压历程特性及其影响因素的研究还未见公开报道。本文中,利用可同时模拟冲击波超压、动压特征的激波管平台,开展入射角对热层前驱波影响的研究实验,并建立数值仿真模型,分析入射角对热层前驱波形成和超压、动压历程特性的影响机理,以期研究结果可为强爆炸热层前驱波毁伤效应研究提供支撑。

    实验中使用的爆炸波模拟激波管如图1所示。该激波管驱动段为直径100 mm的圆柱体,实验段横截面为边长234 mm的正四边形。在实验段上设有测试光窗,光窗中心距实验段起点1.75 m。在光窗前15 cm处设有压力测点,用于监测入射冲击波的强度。光窗两侧设有200 mm的纹影仪。实验时将楔形模型固定在光窗位置,实验现场如图2所示。加热楔形模型产生厚约1 cm的热层,利用纹影系统记录冲击波与热层的作用过程,冲击波波阵面与热层之间的夹角β即为入射角,如图3所示。

    图  1  实验设计
    Figure  1.  Experimental design
    图  2  实验现场
    Figure  2.  Experimental site
    图  3  纹影图像
    Figure  3.  Schlieren image

    在入射波强度和热层状态一致的前提下,本文中通过改变楔形模型冲击波来流方向一侧的角度改变冲击波入射角,进行了4组8次实验,实验条件如表1所示。通过对比不同入射角时有热层和无热层的纹影图像,判断是否出现前驱波,得到热层温度300 ℃对应的临界入射角(简称临界角)的范围,并进一步分析入射角对前驱波的影响。

    表  1  实验条件
    Table  1.  Experimental conditions
    实验
    编号
    入射冲击波
    压力/kPa
    入射角
    β/(°)
    温度/
    有无
    热层
    1-1507515
    1-2300
    2-1506015
    2-2300
    3-1504515
    3-2300
    4-1503015
    4-2300
    下载: 导出CSV 
    | 显示表格

    图3可以看出:冲击波进入热层时,楔形模型表面的热层已基本稳定,高温楔形模型对周围空气的影响较小,仅在图像的右上角产生了轻微的扰动,而冲击波与热层相互作用的区域并未受到影响。因此,分析过程中对流以及高温楔形模型对周围空气的影响忽略不计。

    截取4组实验所得纹影图像中波阵面周围50 mm×50 mm范围进行局部放大,如图4所示,依据Miller等[7]得到的前驱波出现的证明判断是否出现前驱波。当入射角为75°和60°时,冲击波与楔形模型作用发生马赫反射,在有热层时形成了清晰的前驱波,如图4(a)~(b)所示。根据冲击波反射理论[13],超压为50 kPa的冲击波与楔形模型作用发生马赫反射的临界角为46°。当入射角为45°时,由于接近马赫反射临界角,因此无热层时,只观察到楔形模型表面的反射波波阵面较厚,三波点位置不清晰;但当楔形模型表面有热层时,可以看到有较弱的前驱波形成,说明入射角45°接近前驱波形成的临界角,如图4(c)所示。当入射角为30°时,冲击波与楔形模型作用发生规则反射,有热层时无前驱波出现,说明形成前驱波的临界角大于30°。按照形成前驱波的临界角的正弦等于热层内外声速的比[9],计算得到热层温度为300 ℃时形成前驱波的临界角为43°。通过实验发现,形成前驱波的临界角为30°~45°,这与理论结果可以互相验证。

    图  4  不同入射角度、有无热层实验纹影图像对比
    Figure  4.  Comparison of experimental schlieren images with thermal layer with ones without thermal layer at different incident angles

    借鉴Ethridge等的方法[3],以前驱波与马赫杆的相对位置作为前驱波特性的表征量之一。对比有热层和无热层时纹影图像的差异,以无热层时马赫杆所在位置为基准,将前驱波超前马赫杆的距离(leading distance,LD)定义为Dl,如图5所示。以马赫杆到达斜面中点的时刻为准,入射角为45°且无热层时,三波点的位置不清晰,Dl无法确定。入射角为60°和75°时,Dl分别为3.7和7.9 mm。

    图  5  前驱波超前马赫杆的距离示意图
    Figure  5.  Diagram of the distance of the precursor wave leading the Mach stem

    入射角为60°和75°时,冲击波与楔形模型作用发生马赫反射。记马赫杆到达斜面中点的时刻为t1,此时马赫杆与楔形模型斜面的交点位置为X1,记t1t时刻马赫杆与楔形模型斜面的交点位置为X2,则2个时刻的位置差即为马赫波在Δt内的位移D。位移D和时间差Δt已知,且时间差Δt足够小,即可得到t1时刻马赫波的波速vM。同理,可得相同时刻的前驱波波速vp,进一步得到前驱波波速与马赫波波速之差Δv=vpvM

    在本实验中,纹影视场直径为200 mm,所用高速相机的拍摄帧频为39 000 s−1。入射角为75°和60°时前驱波和马赫波的波速及波速差如表2所示。在有热层的情况下,入射角为75°和60°时冲击波前驱波波速相近,但在无热层的情况下入射角为75°的冲击波马赫波的波速远低于入射角为60°的冲击波马赫波的波速。这表明:随入射角的增大,前驱波波速变化较小,马赫波波速大幅降低,前驱波与马赫波的波速差增大。这是导致入射角为75°时的Dl大于入射角为60°时的Dl的原因。

    表  2  不同入射角时前驱波、马赫波波速及波速差
    Table  2.  Precursor and Mach wave velocities as well as their differences at different incident angles
    β/(°)vp/(m∙s−1)vM/(m∙s−1)Δv/(m∙s−1)
    75493.36447.5345.83
    60500.77478.9821.79
    下载: 导出CSV 
    | 显示表格

    为深入研究入射角对热层前驱波压力特征的影响,本文中结合实验构型建立仿真模型,采用分步计算法开展研究。第1步,如图6所示,建立与实验完全相同的完整二维轴对称激波管模型,并在模型中与实验时压力测点对应的位置设数据监测点,得到该位置的总压、静压和总温数据。第2步,以数据监测点为起点,从完整模型中截取部分建立二维平面模型,部分模型长度L为2 m。部分模型左边界为压力入口边界,第一步计算得到的总压、静压和总温数据为该边界的输入;部分模型右边界为无反射压力出口边界;上下边界为绝热壁面,温度与实验时激波管壁温保持一致。在部分模型下壁面距压力入口15 cm处添加倾斜固壁边界模拟楔形模型,斜壁上方厚度1 cm的空气域为热层,如图7所示。依据Griffith[11]得到的热层温度分布公式,设置热层区域内的空气温度沿垂直斜面方向呈指数分布,靠近固壁的空气温度最高为300 ℃。在计算过程中,对气体采用理想气体模型[14],对湍流选用标准k-ε湍流模型[14],忽略对流对冲击波的影响。

    图  6  在完整二维轴对称模型中截取部分模型
    Figure  6.  Interception of the partial model from a complete two-dimensional axisymmetrical model
    图  7  部分模型中添加斜面及高温热层
    Figure  7.  Addition of inclined plane and high-temperature thermal layer in the partial model

    数值仿真方法的验证结果如图89所示。从图8可以看出,在相同位置,数值仿真所得入射波超压和动压曲线与实验所得压力测点处超压和动压曲线均吻合较好。图9为入射角75°时数值仿真所得压力云图,可以看到有热层时出现清晰的前驱波,与实验现象一致。对比实验和数值仿真的LD发现:数值仿真结果略大于相同入射角时的实验结果。入射角为60°时,LD的实验值为3.7 mm,LD的数值仿真值为4.8 mm,两者的相对误差约为22%;入射角为75°时,LD的实验值为7.9 mm,LD的数值仿真值为8.3 mm,两者的相对误差约为5%。入射角60°和75°对应的LD的大小关系数值仿真结果与实验结果一致。综上说明,利用分步计算法研究入射角对前驱波的影响是可行的。

    图  8  入射波超压和动压随时间的演化
    Figure  8.  Evolution of incident overpressure and dynamic pressure with time
    图  9  有无热层时压力云图数值仿真结果
    Figure  9.  Numerically-simulated pressure clouds with and without thermal layer

    本文中分别计算了入射角为90°、75°、60°、45°和30°的5组工况。与文献[15]中分析不同距离处前驱波参数的方法相似,选取斜面中心处距壁面0.1 mm的位置为压力测点,以冲击波到时、超压峰值、动压峰值和动压冲量为参数分析入射角对冲击波超压和动压特性的影响。

    入射角对冲击波到时及其提前程度的影响如图10所示。从图10(a)可以看出:随着入射角的增大,冲击波到时逐渐延迟,延迟增速则随入射角的增大而减小;与无热层时冲击波到时相比,有热层时前驱波到时提前且提前量随入射角的增大而增大。从图10(b)可以看出:入射角为30°时,无前驱波出现,有热层和无热层时冲击波到时相同;入射角为45°时,冲击波到时提前0.015 ms;入射角为90°时,冲击波到时提前0.065 ms。

    图  10  入射角对冲击波到时及其提前程度的影响
    Figure  10.  Influences of incident angle of shock wave on arrival time and its advance degree

    此外,从图10(a)可以看出,有热层时入射角60°和75°的冲击波到时相差仅0.005 ms,但无热层时二者相差0.035 ms。这表明入射角60°和75°的前驱波波速相近,而入射角为75°时的马赫波波速小于入射角为60°时的马赫波波速。相较于入射角为60°时,入射角为75°时前驱波波速与马赫波波速的差值更大,所以冲击波到时提前量更大,这与在第2节得出的结论一致。

    入射角对超压峰值及其减小程度的影响如图11所示。从图11(a)可以发现,随入射角的增大,超压峰值减小;有热层时测点处的超压峰值小于无热层时该点的超压峰值。出现前驱波后,超压峰值减小的程度(超压峰值差)随入射角的增大而减小,如图11(b)所示。入射角为45°时,超压峰值差最大,约为8 kPa;入射角为90°时,超压峰值差最小,约为2 kPa;当入射角为30°、无前驱波出现时,测点处超压峰值仍减小,超压峰值差约为4 kPa。

    图  11  入射角对超压峰值及其减小程度的影响
    Figure  11.  Influence of incident angle of shock wave on peak overpressure and its decreasing degree

    入射角对动压峰值的影响如图12所示。从图12(a)可以看出,有热层存在时,测点处动压峰值增大。从图12(b)可以看出,动压峰值增大的程度(动压峰值差)整体上逐渐增大,当入射角达到一定阈值后开始在一定范围内波动。Ekler等[16]研究发现,有热层时,动压峰值增大的主要原因是热层中的波后粒子速度大幅度提高。图13为入射角为60°时粒子速度和气流密度随时间的演化曲线。可以看到,测点处的粒子速度在冲击波到达后迅速上升到最大值然后逐渐降低,如图13(a)所示;气流密度与之不同,是先上升到一个较小的平台,然后继续增大,如图13(b)所示。这是由于,有热层时前驱波到达以后,经过一段时间入射波后的高密度气体才到达[7]。这导致了当气流密度达到峰值时,粒子速度已经在下降过程中了。粒子速度与气流密度上升过程不同步,使有热层时测点处动压的峰值时刻和峰值大小产生了偏差。因此,动压峰值差在入射角达到一定阈值后开始在一定范围内波动。

    图  12  入射角对动压峰值及其增大程度的影响
    Figure  12.  Influence of incident angle of shock wave on peak dynamic pressure and its increase degree
    图  13  入射角为60°时测点的粒子速度和气流密度随时间的演化
    Figure  13.  Evolutions of particle velocity and airflow density with time at the measured point when the incident angle of shock wave is 60°

    图14为入射角对动压冲量的影响。有热层存在时,测点处动压冲量增大。动压冲量增大的程度(动压冲量差)随入射角的增大而增大。入射角为30°时,无前驱波出现,动压冲量差仅为0.57 Pa·s;入射角为90°时,动压冲量差达到了3.92 Pa·s,约为无热层时动压冲量的17.7%。

    图  14  入射角对动压冲量及其增大程度的影响
    Figure  14.  Influences of incident angle of shock wave on dynamic pressure impulse and its increase degree

    利用爆炸波模拟激波管平台开展实验,研究了热层温度为300 ℃时前驱波出现的临界角,并与理论值进行了比较,进一步采用分布计算法,开展数值仿真研究了热层温度为300 ℃时入射角对热层前驱波压力特征的影响,得到了以下主要结论。

    (1)热层温度为300 ℃时,冲击波与热层相互作用产生前驱波的临界角介于30°~45°,与理论结果一致。入射角越大,前驱波超过马赫杆的距离越大。

    (2)有热层时冲击波到时明显提前,入射角越大,到时提前量越大,入射角为90°时,冲击波到时提前量最大。

    (3)热层会导致超压峰值减小,随入射角的增大,超压峰值差先增大后减小。无前驱波出现时,超压峰值减小的现象仍然存在。

    (4)动压峰值差整体上逐渐增大,当入射角达到一定阈值后,动压峰值差开始在一定范围内波动。这是由于,粒子速度与气流密度上升过程不同步,导致动压峰值时刻和峰值大小不确定,而气流密度的特殊变化则是前驱波与入射波先后到达所致。

    (5)前驱波会导致动压冲量增大,动压冲量差随入射角的增大而增大,入射角为90°时动压冲量差达到最大值。

  • 图  1  实验设计

    Figure  1.  Experimental design

    图  2  实验现场

    Figure  2.  Experimental site

    图  3  纹影图像

    Figure  3.  Schlieren image

    图  4  不同入射角度、有无热层实验纹影图像对比

    Figure  4.  Comparison of experimental schlieren images with thermal layer with ones without thermal layer at different incident angles

    图  5  前驱波超前马赫杆的距离示意图

    Figure  5.  Diagram of the distance of the precursor wave leading the Mach stem

    图  6  在完整二维轴对称模型中截取部分模型

    Figure  6.  Interception of the partial model from a complete two-dimensional axisymmetrical model

    图  7  部分模型中添加斜面及高温热层

    Figure  7.  Addition of inclined plane and high-temperature thermal layer in the partial model

    图  8  入射波超压和动压随时间的演化

    Figure  8.  Evolution of incident overpressure and dynamic pressure with time

    图  9  有无热层时压力云图数值仿真结果

    Figure  9.  Numerically-simulated pressure clouds with and without thermal layer

    图  10  入射角对冲击波到时及其提前程度的影响

    Figure  10.  Influences of incident angle of shock wave on arrival time and its advance degree

    图  11  入射角对超压峰值及其减小程度的影响

    Figure  11.  Influence of incident angle of shock wave on peak overpressure and its decreasing degree

    图  12  入射角对动压峰值及其增大程度的影响

    Figure  12.  Influence of incident angle of shock wave on peak dynamic pressure and its increase degree

    图  13  入射角为60°时测点的粒子速度和气流密度随时间的演化

    Figure  13.  Evolutions of particle velocity and airflow density with time at the measured point when the incident angle of shock wave is 60°

    图  14  入射角对动压冲量及其增大程度的影响

    Figure  14.  Influences of incident angle of shock wave on dynamic pressure impulse and its increase degree

    表  1  实验条件

    Table  1.   Experimental conditions

    实验
    编号
    入射冲击波
    压力/kPa
    入射角
    β/(°)
    温度/
    有无
    热层
    1-1507515
    1-2300
    2-1506015
    2-2300
    3-1504515
    3-2300
    4-1503015
    4-2300
    下载: 导出CSV

    表  2  不同入射角时前驱波、马赫波波速及波速差

    Table  2.   Precursor and Mach wave velocities as well as their differences at different incident angles

    β/(°)vp/(m∙s−1)vM/(m∙s−1)Δv/(m∙s−1)
    75493.36447.5345.83
    60500.77478.9821.79
    下载: 导出CSV
  • [1] GLASSTONE S, DOLAN P J. The effects of nuclear weapons: TID-28061 [R]. 3rd ed. Washington: Department of Defense, 1977.
    [2] ARTEM’EV V I, BERGEL’SON V I, MEDVEDYUK S A, et al. Amplification of intense blast pressure in a medium with this rarefied channels [J]. Fluid Dynamics, 1996, 31(1): 121–126. DOI: 10.1007/BF02230756.
    [3] ETHRIDGE N H, KEEFER J H, CREPEAU J E, et al. Real surface (non-ideal) effects on nuclear explosion airblast from PRISCILLA-type events, part 1: comparison and evaluation of ideal and non-ideal airblast from priscilla computations; and part 2: sharc hydrocode calculations of the priscilla event (phase 1). Final report, 12 March 1992-30 May 1995: AD-A-302079/9/XAB [R]. Aberdeen: Applied Research Associates, Inc. , 1995.
    [4] FRY M A, KAMATH P, BOOK D L. Priscilla calculations and comparison with data: AD-A155 666 [R]. Washington: Naval Research Laborarory, 1985.
    [5] SWIFT L M, SACHS D C, KRIEBEL A R. Air-blast phenomena in the high-pressure region: WT-1403 [R]. Menlo Park, California: Stanford Research Institute, 1960.
    [6] BRYANT E J , ALLEN F. Dynamic pressure impulse for near ideal and non-ideal blast waves: DNA-82-C-0287 [R]. Alexandria, VA: Defense Nuclear Agency, 1983.
    [7] NEEDHAM C E. Blast waves [M]. Berlin, Heidelberg: Springer, 2010.
    [8] ZASLAVSKII B I, MOROZKIN S Y, PROKOF’EV A A, et al. Flow of a planar shock wave around a thermal layer at a rigid wall [J]. Journal of Applied Mechanics and Technical Physics, 1990, 31(3): 354–361. DOI: 10.1007/BF00864563.
    [9] 乔登江. 核爆炸物理概论 [M]. 北京: 国防工业出版社, 2003.
    [10] 贾雷明, 田宙. 平面冲击波与热层作用数值研究 [J]. 现代应用物理, 2015, 6(4): 305–311. DOI: 10.3969/j.issn.2095-6223.2015.04.013.

    JIA L M, TIAN Z. Numerical simulation of the interactions between planar shocks and thermal layer [J]. Modern Applied Physics, 2015, 6(4): 305–311. DOI: 10.3969/j.issn.2095-6223.2015.04.013.
    [11] GRIFFITH W C. Interaction of a shock wave with a thermal boundary layer [J]. Journal of the Aeronautical Sciences, 1956, 23(1): 16–22. DOI: 10.2514/8.3495.
    [12] GION E J. Plane shock interacting with thermal layer [J]. Physics of Fluids, 1977, 20(4): 700–702. DOI: 10.1063/1.861928.
    [13] BEN-DOR G. Shock wave reflection phenomena [M]. 2nd ed. Berlin: Springer, 2007. DOI: 10.1007/978-3-540-71382-1.
    [14] 王福军. 计算流体动力学分析: CFD软件原理与应用[M]. 北京: 清华大学出版社, 2004.
    [15] NEEDHAM C E, EKLER R G, KENNEDY L W. Extended desert calculation results with comparisons to PRISCILLA experimental data and a near-ideal calculation: ARL-CR-235 [R]. Albuquerque: S-Cubed, 1995.
    [16] EKLER R G, NEEDHAM C E, KENNEDY L W. Extended grassland calculation results with comparisons to PRISCILLA experimental data and a near-ideal calculation: AD-A-298081/1/XAB [R]. Albuquerque: S-Cubed, 1995.
  • 加载中
图(14) / 表(2)
计量
  • 文章访问数:  164
  • HTML全文浏览量:  43
  • PDF下载量:  44
  • 被引次数: 0
出版历程
  • 收稿日期:  2023-04-03
  • 修回日期:  2024-05-13
  • 网络出版日期:  2024-05-14
  • 刊出日期:  2024-07-15

目录

/

返回文章
返回