基于电流密度的分段计算电爆炸模型研究

于红新 冉汉政 杜涛 丰伟

王殿玺, 郭香华, 张庆明. 聚脲涂覆钢板在爆炸载荷作用下的动态响应[J]. 高压物理学报, 2019, 33(2): 024103. doi: 10.11858/gywlxb.20180650
引用本文: 于红新, 冉汉政, 杜涛, 丰伟. 基于电流密度的分段计算电爆炸模型研究[J]. 高压物理学报, 2018, 32(6): 062401. doi: 10.11858/gywlxb.20180537
WANG Dianxi, GUO Xianghua, ZHANG Qingming. Dynamic Response of Polyurea Coated Steel Plate under Blast Loading[J]. Chinese Journal of High Pressure Physics, 2019, 33(2): 024103. doi: 10.11858/gywlxb.20180650
Citation: YU Hongxin, RAN Hanzheng, DU Tao, FENG Wei. Electrical Explosion Model Based on Current Density in Piecewise Calculation[J]. Chinese Journal of High Pressure Physics, 2018, 32(6): 062401. doi: 10.11858/gywlxb.20180537

基于电流密度的分段计算电爆炸模型研究

doi: 10.11858/gywlxb.20180537
基金项目: 

中国工程物理研究院电子工程研究所创新发展基金 S20160801

详细信息
    作者简介:

    于红新(1988-), 男, 硕士, 研究实习员, 主要从事脉冲功率技术研究. E-mail:yuhongxin1221@126.com

  • 中图分类号: TJ450.2;TM133

Electrical Explosion Model Based on Current Density in Piecewise Calculation

  • 摘要: 瞬态电阻是爆炸丝在电爆炸过程中最重要的参数之一,准确描述爆炸丝电阻的非线性时变特性是电爆炸仿真的关键,决定了电爆炸过程的仿真精度。为解决传统仿真模型精度差的问题,提出了根据电流密度采用分段计算的方法校正传统电阻率-比作用量关系,建立基于电流密度的分段计算电爆炸仿真模型。结果表明,在大脉冲电流条件下,基于电流密度的分段计算电爆炸模型精度较传统的仿真模型精度明显提高,可用于高精度的爆炸丝电爆炸过程研究。

     

  • 在当今世界的大环境下,无论是民用领域还是军事领域,爆炸冲击的防护都变得越来越重要。传统防护爆炸冲击的方法是提高单层金属防护介质的厚度,但是单纯增加防护介质的厚度不仅会增加成本,还会增加整个防护结构的重量,降低装备的机动性。因此,研究由多层材料所组成的防护结构具有非常重要的意义。聚脲弹性体作为一种新兴的聚合物材料,具有良好的抗爆炸冲击性能,被广泛应用在建筑和舰艇等大型防护结构中。美国海军陆战队研究部门通过研究后发现,在军车底部喷涂聚脲可以在很大程度上减少爆炸冲击对人员和车辆的破坏[1]。宋彬等[2]对比了聚脲弹性体夹层、无夹层和橡胶夹层3种防爆罐在相同质量炸药爆炸载荷作用下的抗爆性能,发现聚脲弹性体夹层防爆罐的抗爆性能优于其他两种防爆罐。甘云丹等[3]对比分析了单层钢板与聚脲涂覆钢板在水中爆炸下的抗爆性能,发现聚脲弹性体可以明显提高钢板的抗爆能力。Bahei-El-Din等[4]发现,在用纤维复合材料为夹层、泡沫为芯层的传统三明治结构中加入聚脲层可以提高其防护性能。Amini等[5]通过分析单位厚度的破坏能量后得出:在钢板的背面涂覆聚脲时,可以较明显地提高钢板的抗爆性能。赵鹏铎等[6]对不同涂覆方式下单钢板和箱体结构的抗爆性能进行了较详细的试验研究,结果表明:在等面密度条件下,钢板迎爆面涂覆聚脲并不能提高其抗爆性能;而在等钢板厚度条件下,在单层钢板和箱体结构上涂覆聚脲能够提高其抗爆性能,且聚脲涂覆在钢板背面时其抗爆性能更好。Tekalur等[7]对由乙烯酯玻璃纤维、聚脲等材料组成的复合结构进行了抗爆性能研究,结果表明相较于单一的乙烯酯玻璃纤维板,聚脲涂层可以明显提高其抗爆性能,且采用聚脲作为夹心层的三明治结构具有最好的抗爆性能。LeBlanc等[8]研究了聚脲层的厚度和位置对聚脲、E-Glass/Epoxy复合结构防护性能的影响,发现相对于聚脲涂覆在受冲击测,聚脲涂覆在背侧更能大幅度提高防护结构的抗爆性能。Ackland等[9]采用局部爆炸加载方法研究了相同面密度聚脲涂覆钢板复合结构的抗爆性能,结果表明单一钢板结构的抗爆性能优于聚脲涂覆结构的抗爆性能。大多数学者的研究显示,在等厚度钢板上涂覆聚脲层可以在一定程度上提高钢板的抗冲击性能,但是定量地研究爆心距、炸药质量以及涂覆聚脲厚度对聚脲涂覆钢板结构抗爆性能的响应规律还比较少。

    本研究利用LS-DYNA有限元软件分析聚脲涂覆位置对钢板抗爆性能的影响,并结合量纲分析方法讨论爆心距、炸药质量、涂覆聚脲厚度对钢板抗爆性能的响应规律,以期为聚脲弹性体涂覆结构抗爆性能的评估与防护结构设计提供参考。

    聚脲涂覆钢板在爆炸载荷作用下的响应涉及多个物理量,可通过量纲分析方法得到各物理量之间的关系。图1为球形TNT炸药与聚脲涂覆钢板的相互作用示意。

    图  1  球形TNT与聚脲涂覆钢板
    Figure  1.  Schematic of spherical TNT and polyurea coated steel plate

    在没有发生断裂的情况下,钢板在爆炸载荷作用下的最大位移可作为其抗爆性能的重要指标之一,而决定最大位移ω的物理量包含以下几个方面。

    (1)球形炸药的物理量参数,包括:炸药量Q,装药密度ρe、单位质量炸药所释放的化学能Ee、爆轰产物的膨胀指数γe

    (2)钢板的物理量参数,包括:密度ρs、弹性模量Es、泊松比vs、屈服极限σs、厚度D、钢板尺寸l

    (3)聚脲的物理量参数,包括:密度ρp、弹性模量Ep、泊松比vp、拉伸强度σp、厚度d、尺寸lp

    (4)钢板与聚脲之间的粘结力,包括:拉伸强度σ、剪切强度τ

    (5)空间几何参数:球形装药质心到靶板的垂直距离h

    基于以上各物理量参数,钢板最大位移ω与各物理量之间的关系函数可写为

    ω=f(Q,ρe,Ee,γe;ρs,Es,vs,σs,D,l;ρp,Ep,vp,σp,d,lp;σ,τ;h)
    (1)

    取装药密度ρe、单位质量炸药所释放的化学能Ee、钢板厚度D作为基本量,(1)式可以化为

    ωD=f(3Q/ρeD,γe;ρsρe,EsρeEe,vs,σsρeEe,lD;ρpρe,EpρeEe,vp,σpρeEe,dD,lpD;σρeEe,τρeEe;hD)
    (2)

    如果保持钢板、聚脲和炸药的材料性质以及钢板的几何尺寸不变,即

    (ρe,Ee,γe;ρs,Es,vs,σs,D,l;ρp,Ep,vp,σp)=const
    (3)

    (1)式就能变为更简单的无量纲关系式

    ωD=f(3Q/ρeD;dD,lpD;σρeEe,τρeEe;hD)
    (4)

    当保持聚脲的几何尺寸及厚度、钢板与聚脲之间的粘结力、炸高不变时,(4)式化为简单的无量纲关系

    ωD=f1(3Q/ρeD)
    (5)

    当保持聚脲的尺寸及厚度、钢板与聚脲之间的粘结力、炸药量不变时,(4)式可变为简单的无量纲关系,即

    ωD=f2(hD)
    (6)

    当保持聚脲涂层的尺寸lp、钢板与聚脲之间的粘结力、炸药量、炸高不变时,(4)式可化为

    ωD=f3(dD)
    (7)

    上面几个式子表明,当保持钢板与聚脲之间的粘结力不变时,钢板的最大位移只和炸高、炸药质量、聚脲涂层厚度有关,具体的无量纲关系式可通过试验或数值模拟确定。下面将用一系列数值仿真来确定相关函数关系。

    为了验证所用模拟方法的有效性,基于LS-DYNA有限元软件对文献[6]中爆炸载荷作用下聚脲涂覆钢板的试验结果进行仿真分析。试验中聚脲涂覆底材为边长400 mm的方形Q235A钢板,板厚1.2 mm,钢板迎爆面涂覆6 mm厚聚脲层。聚脲涂覆钢板结构固定在两个夹板之间,40 g柱形炸药在50 mm爆心距下对聚脲涂覆钢板的加载区域为250 mm×250 mm,更详细的试验方案见文献[6]。

    在数值模拟建模中,由于爆炸载荷和聚脲涂覆钢板结构具有对称性,所以只建立1/4模型,并且在模型的边界处施加固支约束和对称约束。试验中聚脲涂覆钢板的变形区域为250 mm×250 mm,为了节省计算时间,分析时只考虑聚脲涂覆钢板的实际变形区域。钢板、聚脲层用Lagrange网格离散,聚脲层采8节点实体单元SOLID164,而钢板采用壳单元的建模方式。聚脲涂覆钢板与炸药的有限元模型及尺寸如图2所示。

    图  2  聚脲涂覆钢板结构与炸药模型
    Figure  2.  Schematic of polyurea-coated structure and explosive

    Q235A钢采用动力学塑性本构模型(*MAT_PLASTIC_KINEMATIC)。由于聚脲是一种超黏弹性体,具有很明显的应变率效应,故采用6参数的Mooney-Rivlin模型[10],具体的材料参数列于表1

    表  1  聚脲弹性体材料参数[6]
    Table  1.  Parameters of polyurea elastomer material[6]
    Density/(kg·m–3)Tensile strength/MPaElongation at break/%Tear strength/(kN·m–1)Shore hardness (HA)
    1020244008585–95
    下载: 导出CSV 
    | 显示表格

    用关键字*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE_TIEBREAK[11]定义钢板与聚脲层之间的黏结力,设置两种材料之间的黏结强度。计算中当钢板与聚脲界面上的应力大于它们之间的黏结强度时,钢板与聚脲发生分层破坏。

    图3给出文献[6]中工况3的靶板变形情况,图4为相应的数值模拟位移云图。可见,中心处最大位移的试验结果与仿真结果吻合得非常好,验证了本研究所采用模型的有效性。

    图  3  试验中聚脲涂覆钢板的变形[6]
    Figure  3.  Deformation of polyurea coated steel plate in test[6]
    图  4  数值模拟中聚脲涂覆钢板的变形
    Figure  4.  Deformation of polyurea coated steel plate in simulation

    首先对聚脲与钢板组成的具有相同面密度的单一钢板结构、正面涂层、背面涂层和正背面均有涂层4种结构在爆炸载荷作用下的动态响应进行分析,并从钢板的最大位移、塑性应变两个角度分析不同结构的抗爆性能,具体结构设计方案见表2。复合结构的尺寸均为800 mm×800 mm(方形板),钢板厚度为3 mm,聚脲涂层的厚度为4 mm,球形炸药位于方形板中心位置的正上方,炸高为200 mm。

    表  2  结构设计方案
    Table  2.  Structural design schemes
    Condition
    number
    Structural
    diagram
    Steel plate
    thickness/mm
    Polyurea
    thickness/mm
    TNT
    dose/kg
    Detonation
    distance/mm
    A1
    Contrast condition
    301200
    A2
    Single steel plate
    3.5201200
    A3
    Front coating
    341200
    A4
    Back coating
    341200
    A5
    Double sided coating
    32+21200
     Note: Q235A steel plate; Polyurea.
    下载: 导出CSV 
    | 显示表格

    钢板为Q235A钢,采用动力学塑性本构模型(*MAT_PLASTIC_KINEMATIC),强化方式为随动强化,并且考虑材料的应变率效应,其屈服应力放大系数的计算公式为

    σdσY=1+(˙εC)1P
    (8)

    式中:σY为屈服应力,σd为动态屈服应力,˙ε为应变率,PC为材料应变率效应常数。Q235A钢的具体参数见表3

    表  3  Q235A钢材料参数[12]
    Table  3.  Parameters of Q235A steel material[12]
    Density/(kg·m–3)Young’s modulus/GPaPoisson’s ratioσY/MPaTangent modulus/GPaC/s–1P
    78502100.324021055
    下载: 导出CSV 
    | 显示表格

    上述5种防护结构在相同当量爆炸载荷作用下的变形情况如图5所示。球形炸药位于防护结构中心位置的正上方,所以钢板中心区域的位移最大。为了更好地对比钢板的最大位移和塑性变形情况,提取钢板中心区域单元的位移、塑性应变随时间的变化,如图6图7所示。

    图  5  不同工况中钢板的位移云图
    Figure  5.  Displacement cloud chart of steel plate in different working conditions
    图  6  钢板最大变形处的位移时程曲线
    Figure  6.  Displacement-time curves at maximum deflections of steel plates
    图  7  钢板最大变形处的塑性应变时程曲线
    Figure  7.  Plastic strain-time curves at maximum deflections of steel plates

    图6图7可知,纯钢板、无涂覆层的A1工况中的钢板变形大于其余工况,且A3>A5>A4>A2, A3、A5、A4、A2工况中钢板的最大变形分别比A1减小了7.0、7.7、7.9、14.0 mm。比较钢板中心区域的最大塑性应变得到,A1>A3>A5>A4>A2,且A3、A5、A4、A2工况中钢板的最大塑性应变分别比A1减小了11.2%、12.4%、15.5%、16.9%。这说明在等厚度钢板上涂覆一定厚度的聚脲材料能够提高钢板的抗爆性能,并且聚脲涂覆在钢板背爆面的抗爆效果优于其他两种方式的涂覆。背面涂层结构受到爆炸冲击时,聚脲层的变形、破裂等弥散了较多冲击能量,从而减小钢板的变形;而当聚脲涂覆在迎爆面时,聚脲层首先受到爆炸冲击,聚脲层与钢板发生大面积的分层现象,且聚脲层发生撕裂脱落,从而在一定程度上减弱了爆炸冲击对钢板的作用:这两种涂覆方式都对爆炸冲击能量有一定的耗散作用,只是两者相比,背面涂覆方式更好。在等面密度条件下(A2~A5),纯钢板结构比聚脲涂覆钢板复合结构的抗爆性能更好,说明聚脲更应该涂覆在已有结构表面而不是用于设计等面密度的复合结构。

    由2.3节可知,在等厚度钢板上涂覆一定厚度的聚脲层能够提高钢板的抗爆性能,且聚脲涂覆在钢板背面的方式更好。本节以A4工况为例,用量纲方法分析爆心距、炸药质量、聚脲层厚度对聚脲涂覆钢板变形的影响。

    保持炸药质量、聚脲层厚度不变,改变球形炸药质心距钢板中心位置的距离(0.2~0.8 m)。首先对得到的钢板中心区域最大位移数据进行无量纲化,再用指数函数对无量纲数据进行拟合,得到如图8所示的无量纲化关系曲线。

    图  8  爆心距对聚脲涂覆钢板最大位移的影响
    Figure  8.  Relationship between the detonation distance and the maximum displacement of polyurea coated steel plate

    图8中可以看出,在保持其他变量不变的条件下,随着炸药质心与钢板之间距离的增加,钢板的最大位移近似呈指数形式递减,拟合关系式为ω/D=2.55+40.3e0.005h/D。由球形炸药在无限空气介质中爆炸的超压、比冲量经验公式可知,当保持炸药质量不变时,随着爆心距的增加,结构受到的超压峰值和冲量下降,爆炸载荷传递到结构上的能量也随之减少,因此聚脲涂覆钢板的最大变形量随爆心距的增加而近似呈指数形式递减。

    同样以A4工况为例,保持涂层聚脲厚度和炸高不变,改变炸药的质量,分析炸药质量对聚脲涂覆钢板最大位移的影响规律。和3.1节一样,首先对得到的聚脲涂覆钢板最大位移数据进行无量纲化,再用线性函数拟合无量纲数据,得到如图9所示的关系曲线。

    图  9  炸药质量对聚脲涂覆钢板最大位移的影响
    Figure  9.  Relationship between the mass of explosive and the maximum displacement of polyurea coated steel plate

    图9中可以看出,在保持其他变量不变的情况下,随着球形炸药质量的增加,钢板中心区域的最大位移近似呈线性增加趋势,拟合线性关系式为ω/D=2.333QρeD332.75。同样从球形炸药在无限空气介质中爆炸的超压、比冲量经验公式可知,当保持爆心距不变时,随着炸药质量的增加,爆炸载荷作用在结构上的超压和比冲量增大,传递到复合结构上的能量也随之增加,因此聚脲涂覆钢板的最大位移随炸药质量的增加近似呈线性增长趋势。

    保持炸药质量、炸高不变,改变涂覆在钢板背面的聚脲层的厚度。同样对该条件下得到的所有钢板的最大位移数据进行无量纲化,用线性函数对无量纲数据进行拟合,得到如图10所示的关系曲线。

    图  10  聚脲层厚度对聚脲涂覆钢板最大位移的影响
    Figure  10.  Relationship between the thickness of polyurea layer and the maximum displacement of polyurea coated steel plate

    图10中可以看出,在保持其他变量不变的情况下,随着涂覆在钢板背面聚脲层厚度的增加,钢板的最大位移近似呈线性减小趋势,拟合关系式为ω/D=33.981.372(d/D)。钢板的最大变形随着涂覆聚脲层厚度的增加而减小的主要原因是:增加聚脲层厚度在一定程度上增加了结构的刚度,从而整体上提高了结构的抗爆性能;将聚脲涂覆在背面时,由于聚脲材料的黏弹性特性,复合结构在受到爆炸冲击后,一部分以破裂碎片动能的形式带走大量冲击能量,一部分能量被聚脲层弥散和吸收,且聚脲层越厚,其破裂的程度越严重,破片的数量越多,以动能形式耗散的能量越大,所以聚脲层越厚,其吸收和弥散的能量越多。

    通过数值仿真分析了聚脲涂覆方式对钢板抗爆性的影响,结合量纲分析方法讨论了爆心距、炸药质量、涂覆聚脲厚度对钢板中心最大位移的影响,得到以下结论:

    (1)在等厚度钢板上涂覆聚脲能够提高钢板整体的抗爆性能,且聚脲涂覆在钢板背爆面的效果优于其他两种方式,而在等面密度条件下,纯钢板结构比聚脲涂覆钢板复合结构的抗爆性能更好;

    (2)在保持其他变量不变的情况下,随着爆心距的增大,聚脲涂覆钢板的最大位移近似呈指数形式递减;

    (3)改变炸药质量,聚脲涂覆钢板的最大位移近似随着炸药质量的增加而呈线性增加;

    (4)改变涂覆在钢板背爆面的聚脲厚度,其他条件保持不变,则聚脲涂覆钢板的最大位移随着涂层厚度的增加近似呈线性递减趋势。

  • 图  RLC放电电路

    Figure  1.  Discharge circuit containing R, L and C parameters

    图  典型的RLC放电电流波形

    Figure  2.  Typical RLC discharge current waveform

    图  典型电阻率-比作用量曲线

    Figure  3.  Typical resistivity-specific action curves

    图  不同电流密度下电阻率-比作用量曲线

    Figure  4.  Resistivity-specific action curves at different current densities

    图  电阻率曲线(1 TA/m2)

    Figure  5.  Resistivity curve(1 TA/m2)

    图  分段的电阻率-比作用量曲线

    Figure  6.  Resistivity-specific action curves at different current densities

    图  Saber仿真模型

    Figure  7.  Simulation model in Saber

    图  实验示意图

    Figure  8.  Schematic diagram of experimental setup

    图  短路实验电流波形

    Figure  9.  Current waveform of short-circuit discharge experiment

    图  10  短路仿真模型

    Figure  10.  Simulation model of short-circuit

    图  11  爆炸丝结构示意图

    Figure  11.  Schematic diagram of exploding wire structure

    图  12  电流波形

    Figure  12.  Current waveform

    图  13  控制电流曲线

    Figure  13.  Curve of control current

    图  14  仿真与实验电流波形

    Figure  14.  Simulation and experimental current waveforms

    表  1  爆炸丝单点起爆实验与仿真结果

    Table  1.   Experiment and simulation results of explosive wire with single point initiation

    State Peak current/kA Error/% Time of peak current/μs Error/% Burst current/kA Error/% Burst time/μs Error/%
    Experiment 2.088 0.638 1.315 0.251
    Simulation 1 2.066 0.95 0.623 2.35 1.309 0.76 0.253 0.80
    Simulation 2 2.138 2.39 0.594 6.90 1.439 9.42 0.243 3.19
    Simulation 3 1.936 7.18 0.712 11.60 1.298 1.29 0.225 10.36
    Simulation 4 1.600 23.37 0.532 16.61 1.263 3.95 0.206 17.93
    下载: 导出CSV
  • [1] TUCKER T J.Behavior of exploding glod wires[J].Journal of Applied Physics, 1961, 32(10):1894-1900. doi: 10.1063/1.1728259
    [2] TUCKER T J, STANTON P L.Electrical gurney energy: a new concept in modeling of energy transfer from electrically exploded conductors: SAND-75-0244[R].Albuquerque, NM: Sandia National Laboratories, 1975.
    [3] 杨汉武, 钟辉煌.PSpice模型用于电爆炸丝的数值模拟[J].国防科技大学学报, 2000, 22(增刊):38-42. http://d.old.wanfangdata.com.cn/Periodical/gfkjdxxb2000Z1011

    YANG H W, ZHONG H H.Numerical simulation of electrical exploding wires via PSpice models[J].Journal of National University of Defense Technology, 2000, 22(Suppl):38-42. http://d.old.wanfangdata.com.cn/Periodical/gfkjdxxb2000Z1011
    [4] 杨集, 杜涛, 赵翔, 等.基于电容放电单元的爆炸丝起爆特性数值分析[J].太赫兹科学与电子信息学报, 2014, 12(3):461-465. http://d.old.wanfangdata.com.cn/Periodical/xxydzgc201403029

    YANG J, DU T, ZHAO X, et al.Numerical analysis of exploding wire initiation characteristic based on capacitor discharge unit[J].Journal of Terahertz Science and Electronic Information Technology, 2014, 12(3):461-465. http://d.old.wanfangdata.com.cn/Periodical/xxydzgc201403029
    [5] LEE R S.Fireset: UCID-21322[R].Livermore, CA: Lawrence Livermore National Laboratory, 1988.
    [6] FURNBERG C M.Computer modeling of detonators[C]//Proceedings of the 37th Midwest Symposium on Circuits and Systems.Piscataway, NJ: Institute of Electrical and Electronics Engineers, 1994, 1: 646-649.
    [7] IVANENKOV G V, STEPNIEWSKI W.Three-temperature model for the dynamics of a plasma produced by exploding metal wires[J].Plasma Physics Reports, 2000, 26(1):21-32. doi: 10.1134/1.952818
    [8] 王刚华, 胡熙静.单丝Z箍缩的一维单温辐射磁流体力学计算[J].高压物理学报, 2002, 16(1):57-60. doi: 10.3969/j.issn.1000-5773.2002.01.009

    WANG G H, HU X J.One-dimensional simulations of single temperature radiative magnetic field dynamics for liner implosion experiments[J].Chinese Journal of High Pressure Physics, 2002, 16(1):57-60. doi: 10.3969/j.issn.1000-5773.2002.01.009
    [9] 张振涛, 汤铁钢, 王健, 等.爆炸丝起爆系统研制[J].高能量密度物理, 2013, 3(1):35-40. http://www.cnki.com.cn/Article/CJFDTOTAL-QJGY201403053.htm

    ZHANG Z T, TANG T G, WANG J, et al.Development of exploding wire initiation system[J].High Energy Density Physics, 2013, 3(1):35-40. http://www.cnki.com.cn/Article/CJFDTOTAL-QJGY201403053.htm
    [10] TUCKER T J, TOTH R P.EBW1: a computer code for the prediction of the behavior of electrical circuit containing exploding wire elements: SAND-75-0041[R].Albuquerque, NM: Sandia National Laboratories, 1975.
    [11] YING G F, LI X W, JIA S L.Electrical characteristics of microsecond electrical explosion of Cu wires in air under various parameters[J].Transactions on Plasma Science, 2018, 46(4):972-981. doi: 10.1109/TPS.2018.2805351
  • 加载中
图(14) / 表(1)
计量
  • 文章访问数:  6170
  • HTML全文浏览量:  2797
  • PDF下载量:  17
出版历程
  • 收稿日期:  2018-04-09
  • 修回日期:  2018-05-31

目录

/

返回文章
返回