液-液斜界面R-T失稳特征的实验研究

江荣宝 黄熙龙 邹立勇 施红辉 吴鋆

江荣宝, 黄熙龙, 邹立勇, 施红辉, 吴鋆. 液-液斜界面R-T失稳特征的实验研究[J]. 高压物理学报, 2018, 32(5): 054201. doi: 10.11858/gywlxb.20180512
引用本文: 江荣宝, 黄熙龙, 邹立勇, 施红辉, 吴鋆. 液-液斜界面R-T失稳特征的实验研究[J]. 高压物理学报, 2018, 32(5): 054201. doi: 10.11858/gywlxb.20180512
JIANG Rongbao, HUANG Xilong, ZOU Liyong, SHI Honghui, WU Jun. Experimental Investigation on the Characteristics of Unstability at Liquid-Liquid Tilted Interface Induced by Rayleigh-Taylor Instability[J]. Chinese Journal of High Pressure Physics, 2018, 32(5): 054201. doi: 10.11858/gywlxb.20180512
Citation: JIANG Rongbao, HUANG Xilong, ZOU Liyong, SHI Honghui, WU Jun. Experimental Investigation on the Characteristics of Unstability at Liquid-Liquid Tilted Interface Induced by Rayleigh-Taylor Instability[J]. Chinese Journal of High Pressure Physics, 2018, 32(5): 054201. doi: 10.11858/gywlxb.20180512

液-液斜界面R-T失稳特征的实验研究

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

国家自然科学基金 11472253

国家自然科学基金 11672277

国家自然科学基金 11602247

科学挑战专题 TZ2016001

详细信息
    作者简介:

    江荣宝(1993—), 男, 硕士研究生, 主要从事界面不稳定性的实验研究. E-mail: 18367053723@163.com

    通讯作者:

    黄熙龙(1988—), 男, 硕士, 助理研究员, 主演从事界面不稳定性研究. E-mail: xlhuang@caep.cn

  • 中图分类号: O357.5

Experimental Investigation on the Characteristics of Unstability at Liquid-Liquid Tilted Interface Induced by Rayleigh-Taylor Instability

  • 摘要: 实验研究了不同加速度下不相溶流体斜界面Rayleigh-Taylor(R-T)不稳定性诱导湍流混合区的发展特征。利用高速阴影测试技术研究了硅油/碘化钾溶液界面R-T不稳定性诱导界面的演化规律,定量分析了混合区发展宽度、界面倾角与时间的关系。通过实验结果分析得出:水平界面与斜界面混合区宽度演化的差别主要体现在后期,前、中期两者演化规律基本一致;界面倾角的演化呈现出两种不同的趋势,随时间先呈抛物线增长,随后线性增长,演化后期流体界面发生翻转;由R-T不稳定性诱导的界面混合与K-H不稳定性诱导的界面混合呈相互竞争的关系。

     

  • 当轻流体加速重流体或者在重力场中支撑重流体时,流体界面将产生Rayleigh-Taylor(R-T)不稳定性[1]。R-T不稳定性的演化及其诱导的湍流混合是一种复杂的非线性多尺度问题,在惯性约束核聚变(Inertial Confinement Fusion, ICF)中靶丸的设计、内燃机燃油喷雾等技术领域有重要的应用价值。此外,R-T不稳定性对超新星爆发、地球物理学中地下盐丘及火山岛的形成起着重要作用。因此,R-T不稳定性问题得到众多学者广泛的关注[2-5]

    Rayleigh[6]首次研究了R-T不稳定性的发展,分析了重力场中密度不均匀流体的稳定性,获得了增长率与波数的特征关系。Read[7]通过实验研究竖直加速箱体内两种流体R-T不稳定的演化,发现了扰动振幅hi与演化时间t的关系

    hi=αiAgt2
    (1)

    式中:αi是尖钉或气泡的增长系数;A是无量纲参数Atwood数,A=(ρ2-ρ1)/(ρ2+ρ1),ρ1ρ2为轻、重流体的密度;g是恒定加速度。而对于非恒定加速度,用箱体位移Z代替gt2

    Z=g(t)dtdt
    (2)

    式中:g(t)为非恒定加速度。

    Dimonte等[8]利用直线电机(Liner Electric Motor, LEM)实验研究了复杂加速度下不同密度比流体R-T不稳定性的发展过程,得出对于非恒定加速度且满足Ag>0,扰动振幅

    hi=2αiAS
    (3)

    其中

    S=12[gdt]2
    (4)

    随后Dimonte等[9]利用直线电机研究了流体在连续不同阶段的加速度(加速-减速-加速)下的界面失稳发展特征,结果表明:流体在第一次R-T不稳定性阶段(加速阶段),扰动振幅随时间二次增长;在R-T稳定阶段(减速阶段),由于流体发生翻转,扰动振幅降低。Waddell等[10]利用重物和滑轮组系统作为驱动手段实验研究了低Atwood数下具有单模初始扰动的界面R-T不稳定性的演化,发现初期扰动振幅的增长规律与线性稳定理论推导结果较为吻合。Holford等[11]利用挡板实验研究了重力场中流体斜界面R-T不稳定性的发展,但他们主要研究倾斜角度对流体混合效率的影响。McFarland等[12-13]实验研究了激波作用下Air/He-SF6斜界面Richtmyer-Meshkov(R-M)不稳定性受初始扰动的影响,定量研究了混合区宽度增长规律,发现在不同初始界面倾角(30°~60°)下,混合区宽度都随着无量纲时间线性增长。Dimonte等[14]采用不同的数值方法二维或三维数值模拟了小扰动波长下流体界面R-T不稳定性诱导的混合区发展,发现只有气泡直径及增长宽度与实验结果较为吻合。国内黄文斌等[15]利用果冻实验技术研究了气-液界面R-T不稳定性,分析了不同波长扰动对混合区内部结构特征的影响规律。施红辉等[16]利用激波管实验研究了气-液界面R-T不稳定性演化过程。刘金宏等[17]利用竖直激波管高压气体驱动箱体实验研究了一定倾角下混合区的不对称性特征。综上所述,国内外对R-T不稳定性进行了大量详尽的实验及数值模拟研究,但主要集中在水平界面,实际工程应用中界面更加复杂,加速度与密度梯度不总是共线的,如ICF中,由于加工工艺的因素靶丸并不严格呈球状,造成密度的不均匀性,此时烧蚀面与加速度处于不共线的状态,而加速度与密度梯度不共线会伴随K-H不稳定性的产生,使R-T失稳演化更加复杂。因此研究此类问题不仅在靶丸的设计等物理应用上有重要价值,更有助于提升对界面不稳定物理机制的认识。

    本研究基于改造后的激波管,采用高压气体膨胀驱动箱体,箱体内盛有不相溶的两种流体,在重力作用下形成具有自然随机扰动的初始界面。实验中,将实验管道倾斜一定角度形成斜界面,研究界面加速后的R-T不稳定性发展过程,利用阴影测试技术获得不同加速度下初始倾角分别为0°和5°时斜界面的演化规律。

    图 1为实验装置示意图,其中:A为高压驱动段,B为活塞运动段,C为隔离A与B的膜片,D为在实验管道运动的活塞,E为活塞与箱体的连接杆,F为存放实验液体的箱体,G为缓冲垫。实验时,往A中充入一定压力的气体,鼓起的膜片C与预先布置的电热丝贴在一起,电热丝通电后,膜片沿着电热丝的方向破开,气体膨胀驱动活塞及箱体向下运动。箱体上端安装加速度传感器,将加速度信号传送至示波器,示波器将延时信号传至高速相机实现同步拍摄。当加速度方向与流体界面上的密度梯度方向不共线时,称界面为斜界面。实验中,将装置倾斜至一定角度,加速度方向与流体界面具有一定夹角,形成初始斜界面。重力加速度g0与箱体加速度a的夹角为流体界面初始倾角γ0表 1为实验的工况,包括初始界面厚度、初始界面倾角及高压段初始压力。箱体空腔尺寸为60 mm(长)×30 mm(宽)×120 mm(高),即箱体容积V=216 mL,轻、重流体各占一半。采用的轻流体为硅油(密度0.917 g/mL,折射率1.389,黏度11.1 mPa·s),重流体为碘化钾溶液(密度1.289 g/mL,折射率1.389,黏度1.27 mPa·s),轻/重流体界面的Atwood数为0.169。实验在常温(25 ℃)、标准大气压下完成。

    图  1  实验装置示意
    Figure  1.  Schematic of experimental apparatus
    表  1  实验参数
    Table  1.  Experimental parameters
    Case Atwood number Initial interfacial thickness/mm γ0/(°) Initial pressure at high pressure section/MPa
    1 0.169 0.721 0 0.25
    2 0.169 0.752 5 0.25
    3 0.169 0.734 0 0.40
    4 0.169 0.761 5 0.40
    下载: 导出CSV 
    | 显示表格

    流体的折射率与其密度相关,可以用阴影法测试流体折射率的变化。实验流体在加速过程中,由R-T不稳定性诱导的湍流混合区有较大的密度变化,阴影法能对混合区宽度及特征结构进行清晰的显示。图 2为反射式阴影法示意图,点光源光线经球面反射镜M1后得到一束均匀的平行光,平行光经过实验测试窗口,由球面反射镜M2产生汇聚光路,高速相机在焦点后适当位置记录成像信息。本实验高速相机幅频为3 000幅/秒,曝光时间为1 μs。相比于直接拍照,阴影法采用的平行光避免了因俯仰角产生的测量误差,保证了实验数据的准确性。

    图  2  阴影法示意
    Figure  2.  Schematic of shadowgraph

    当高压气体推动箱体运动时,轻流体向重流体加速。由于惯性力的作用,轻流体“下沉”渗透重流体中形成“气泡”,重流体“上浮”混入轻流体中形成“尖钉”。这些结构具有较大的密度梯度,在测试图像中均显示为阴影区域。在该区域中,“尖钉”与“气泡”之间的流体混合区称为湍流混合区。混合区宽度反映了由R-T不稳定性诱导的湍流混合的范围和强度。测得的加速度信号数据经平滑处理后如图 3所示,可知界面加速度在破膜后瞬间达到最大加速度,由于气体驱动力降低及存在管道摩擦力,加速度衰减至零后又反向增大。图 3中加速度正向为重力加速度的方向。在加速度从最大值衰减到零时,加速度在界面法向的分量与密度梯度反向,此时界面是R-T不稳定的;而当箱体加速度进入负增长阶段(加速度方向与重力加速度方向相反)时,加速度在界面法向的分量与密度梯度同向,此时界面是R-T稳定的。本实验结果的时间范围处于从零时刻到箱体第一次撞击缓冲垫之间。

    图  3  加速度曲线
    Figure  3.  Acceleration curves

    图 4是高压段初始压力为0.25 MPa时,水平界面及斜界面(5°)的混合区发展阴影图。在图 4(a)中,随着时间的推移,流体界面处产生“气泡”和“尖钉”,且在前、中期“气泡”和“尖钉”结构的体积逐渐增大,数量逐渐增多,混合区宽度增大。后期界面处“尖钉”和“气泡”的体积达到最大值,但两种结构的数量减少,其原因是“尖钉”和“气泡”在混合区发展后期进行合并形成体积更大的结构,与文献[10]中提到的“气泡合并理论”相吻合。由于箱体壁面的约束,流体在壁面处生成壁面涡旋,且壁面涡旋随着时间的增长逐渐增大。

    图  4  高压段初始压力为0.25 MPa时界面湍流混合区发展阴影图
    Figure  4.  Evolution of turbulent mixing zone when the initial pressure of the high pressure section is 0.25 MPa

    当界面倾斜时,如图 4(b)所示,流体界面倾角随着时间的发展逐渐增大,界面呈翻转趋势。对比图 4(a)发现,同一时刻,壁面涡更加明显,且斜界面处形成的“气泡”和“尖钉”的体积及数量较小,界面更平滑。这是由于加速度方向与流体界面不共线,根据力的合成原理,可将加速度a分解为垂直于流体界面的分量aR-T和平行于流体界面的分量aK-H(见图 5),其中:aR-T诱导界面R-T不稳定性的产生;aK-H导致流体切向流动,使界面更平滑。同时,aK-H导致流体界面存在速度差从而诱导Kelvin-Helmholtz(K-H)不稳定性。R-T不稳定性与K-H不稳定性的耦合作用使流体界面混合区内部结构更复杂。

    图  5  界面加速度分解
    Figure  5.  Acceleration decomposed at the interface

    图 6为高压段初始压力为0.40 MPa时,水平界面及斜界面(5°)的混合区发展阴影图。对比图 4(a)图 6(a)可知:当高压段初始压力增大时,水平界面上“气泡”和“尖钉”的数量增加,混合区结构发展更剧烈;同一时刻,混合区宽度随着高压段初始压力的增大而增大,壁面涡尺寸随着加速度的增大而增大。对比图 4(b)图 6(b)可知,斜界面、高压段初始压力为0.40 MPa时的混合区结构比水平界面、初始压力为0.25 MPa时的更加明显。这是因为加速度垂直于界面的分量aR-T变大,由aR-T引发R-T不稳定性诱导的混合区结构尺寸更大;同时aK-H增大会导致流体界面切向流动增大,产生涡,流体结构发生改变。发展后期,图 4(b)图 6(b)混合区内尖钉数量减少。这是因为初始加速度降低,界面处aR-T降低,界面倾角随着演化时间逐渐增大,界面加速度在垂直界面上的分量aR-T变得更小,aK-H衰减较为缓慢,此时K-H不稳定性在诱导界面混合中占主导地位。

    图  6  高压段初始压力为0.40 MPa时界面湍流混合发展阴影图
    Figure  6.  Evolution of turbulent mixing zone when the initial pressure of the high pressure section is 0.40 MPa

    为了定量地对比高压段初始压力与初始倾斜角对混合区宽度演化的影响,将最大“尖钉”与最大“气泡”切线的垂直距离定义为混合区的宽度W,水平线与斜界面的夹角定义为斜界面倾角γ[18],如图 7所示。图 8为混合区宽度W随时间的演化图,可以看出,在0~20 ms之间,4种工况下混合区宽度增长缓慢,增长曲线基本重合。当高压段初始压力为0.25 MPa时,在20~80 ms之间混合区宽度进入快速增长阶段;但当高压段初始压力为0.40 MPa时,混合区宽度在20~70 ms内进入快速增长阶段,混合区宽度进入后期缓慢增长的时刻提前约10 ms。混合区增长速率随着加速度的增大而增大,同一时刻混合区宽度随着加速度的增大而增大。后期湍流混合区宽度增长速率降低,是由于界面所受到的加速度的衰减导致的。同时因为加速度呈负增长,此时流体加速度在界面法向的分量aR-T与密度梯度同向,界面是R-T稳定的,混合区内部“尖钉”和“气泡”结构不再发展,混合区宽度增长速率降低。

    图  7  斜界面混合区参数定义
    Figure  7.  Parameter definiton of mixing zonein obilique interface
    图  8  混合区宽度W随时间的演化
    Figure  8.  Evolution of width of turbulentmixing zone over time

    对比可知,初始界面倾角对混合区宽度发展的影响主要体现在后期。水平界面混合区宽度最大值大于斜界面混合区宽度最大值。这一方面是因为界面倾斜时,随着时间的发展,界面倾角逐渐增大,加速度在界面切向上的分量aK-H增大,降低了在密度梯度方向混合的能量,导致混合区增长速率变小。另一方面可能因为由R-T不稳定性诱导的界面混合与K-H不稳定性诱导的界面混合呈相互竞争关系,二者的耦合作用限制了最大“气泡”及“尖钉”的发展。

    图 9所示为斜界面倾角γ随时间的演化及其拟合曲线。从图 9中可发现,斜界面倾角为5°时,两种高压段初始压力下界面倾角的变化趋势相同,且初始压力较大时,同一时刻斜界面倾角的增长率较高。高压段初始压力为0.25 MPa时:在0~83 ms之间,界面倾角随时间近似呈抛物线增长; 随后界面倾角随时间线性增长,拟合方程为γ=2.87+0.188t。高压段初始压力为0.40 MPa时:在0~76 ms之间,界面倾角也随时间近似呈抛物线增长;随后呈线性增长,拟合方程为γ=-5.525+0.365t。这两个时间点对应着加速度从最大值衰减到零过程中的时刻,即界面处于R-T不稳定阶段。当高压段初始压力由0.25 MPa提高到0.40 MPa时,界面倾角在线性增长阶段的斜率增加了约1倍。结合图 8可知:混合区宽度随时间线性增长时,混合区界面倾角随时间近似呈抛物线增长; 但当后期混合区宽度不再增长时,混合区界面倾角却随着时间不断增长,且增长趋势近似线性。

    图  9  界面倾角γ随时间的演化
    Figure  9.  Evolution of obliquity of turbulentmixing zone over time

    采用高压气体驱动箱体的方式实验研究了两种不同初始倾角界面在不同加速度下的R-T不稳定性。利用阴影法获得了液-液斜界面R-T不稳定性诱导的界面失稳和混合的演化特征,分析了湍流混合区宽度及界面倾角随时间的定量关系。实验结果发现:加速度越大,混合区宽度增长率及界面倾角增长率越高;倾角效应在后期影响着混合区宽度的发展;混合区宽度不再增长时,其界面倾角一直增加,界面呈翻转趋势。阴影法可以清晰地显示混合区的外部特征,但对混合区内部的涡结构演化无能为力,后续工作将采用PIV-PLIF同步测量技术对混合区速度分布及密度分布进行测试,并将对由R-T不稳定性和K-H不稳定性诱导的湍流混合的耦合效应进行更深入的分析。

  • 图  实验装置示意

    Figure  1.  Schematic of experimental apparatus

    图  阴影法示意

    Figure  2.  Schematic of shadowgraph

    图  加速度曲线

    Figure  3.  Acceleration curves

    图  高压段初始压力为0.25 MPa时界面湍流混合区发展阴影图

    Figure  4.  Evolution of turbulent mixing zone when the initial pressure of the high pressure section is 0.25 MPa

    图  界面加速度分解

    Figure  5.  Acceleration decomposed at the interface

    图  高压段初始压力为0.40 MPa时界面湍流混合发展阴影图

    Figure  6.  Evolution of turbulent mixing zone when the initial pressure of the high pressure section is 0.40 MPa

    图  斜界面混合区参数定义

    Figure  7.  Parameter definiton of mixing zonein obilique interface

    图  混合区宽度W随时间的演化

    Figure  8.  Evolution of width of turbulentmixing zone over time

    图  界面倾角γ随时间的演化

    Figure  9.  Evolution of obliquity of turbulentmixing zone over time

    表  1  实验参数

    Table  1.   Experimental parameters

    Case Atwood number Initial interfacial thickness/mm γ0/(°) Initial pressure at high pressure section/MPa
    1 0.169 0.721 0 0.25
    2 0.169 0.752 5 0.25
    3 0.169 0.734 0 0.40
    4 0.169 0.761 5 0.40
    下载: 导出CSV
  • [1] SHARP D H.An overview of Rayleigh-Taylor instability[J].Physica D:Nonlinear Phenomena, 1984, 12(1):3-18. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=JJ0235350510
    [2] ANDRONOV V A, BAKHRAKH S M, MOKOHOV V N, et al.Effect of turbulent mixing on the compression of laser targets[J].JETP Letters, 1979, 29(1):56-59. http://adsabs.harvard.edu/abs/1979ZhPmR..29...62A
    [3] ALON U, HECHT J, OFER D, et al.Power laws and similarity of Rayleigh-Taylor and Richtmyer-Meshkov mixing fronts at all density ratios[J].Physical Review Letters, 1995, 74(4):534-536. doi: 10.1103/PhysRevLett.74.534
    [4] DIMONTE G.Nonlinear evolution of the Rayleigh-Taylor and Richtmyer-Meshkov instabilities[J].Physics of Plasmas, 1999, 6(5):2009-2015. doi: 10.1063/1.873491
    [5] ATZENI S, MEYER-TER-VEHN J. 惯性聚变物理[M]. 沈百飞, 译. 北京: 科学出版社, 2008: 198-253.

    ATZENI S, MEYER-TER-VEHN J. The physics of inertial fusion[M]. Translated by SHEN B F. Beijing: Science Press, 2008: 198-253.
    [6] RAYLEIGH L.Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density[J].Proceeding of the London Mathematical Society, 2009, 14(1):200-207. doi: 10.1112-plms-s1-14.1.170/
    [7] READ K I.Experiment investigation of turbulent mixing by Rayleigh-Taylor instability[J].Physica D:Nonlinear Phenomena, 1984, 12(1):45-58. http://adsabs.harvard.edu/abs/1984PhyD...12...45R
    [8] DIMONTE G, SCHNEIDER M.Density ratio dependence of Rayleigh-Taylor mixing for sustained and impulsive acceleration histories[J].Physics of Fluids, 2000, 12(2):304-321. doi: 10.1063/1.870309
    [9] DIMONTE G, RAMAPRABHU P, ANDREWS M.Rayleigh-Taylor instability with complex acceleration history[J].Physical Review E, 2007, 76(4):1-6. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=JJ0231719232
    [10] WADDELL J T, NIEDERHAUS C E, JACOBS J W.Experimental study of the Rayleigh-Taylor instability:low atwood number liquid system with single-mode initial perturbations[J].Physics of Fluids, 2001, 13(5):1263-1273. doi: 10.1063/1.1359762
    [11] HOLFORD J M, DALZIEL S B, YOUNGS D.Rayleigh-Taylor instability at a tilted interface in a laboratory experiments and numerical simulations[J].Laser and Particle Beams, 2003, 21(4):419-423. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=JJ025589876
    [12] MCFARLAND J A, GREENOUGH J A, RANJAN D.Investigation of the initial perturbation amplitude for the inclined interface Richtmyer-Meshkov instability[J].Physica Scripta, 2013, 155(21):14-20. http://iopscience.iop.org/1402-4896/2013/T155/014014
    [13] MCFARLAND J A, REILLY D, CREEL S, et al.Experimental investigation of the inclined interface Richtmyer-Meshkov instability before and after reshock[J].Experiments in Fluids, 2014, 55(3):16-40. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=JJ0232341654
    [14] DIMONTE G, YOUNGS D L, DIMITS A.A comparative study of the turbulent Rayleigh-Taylor instability using high-resolution three-dimensional numerical simulations:the alpha-group collaboration[J].Physics of Fluids, 2004, 16(5):1668-1693. doi: 10.1063/1.1688328
    [15] 黄文斌, 邹立勇, 刘金宏, 等.初始扰动对于气液界面Rayleigh-Taylor不稳定性发展的影响[J].实验流体力学, 2010, 24(3):39-41. doi: 10.3969/j.issn.1672-9897.2010.03.008

    HUANG W B, ZOU L Y, LIU J H, et al.Effects of initial perturbations on Rayleigh-Taylor instability growth at gas liquid interface[J].Journal of Experiments in Fluid Mechanics, 2010, 24(3):39-41. doi: 10.3969/j.issn.1672-9897.2010.03.008
    [16] 施红辉, 卓启威.Richtmyer-Meshkov不稳定性流体混合区发展的实验研究[J].力学学报, 2007, 39(3):417-421. doi: 10.3321/j.issn:0459-1879.2007.03.016

    SHI H H, ZHUO Q W.Evolution of the fluid mixing zone in Richtmyer-Meshkov instability at a gas/liquid interface[J].Chinese Journal of the Theoretical and Applied Mechanics, 2007, 39(3):417-421. doi: 10.3321/j.issn:0459-1879.2007.03.016
    [17] 刘金宏, 谭多望, 张旭, 等.斜界面Rayleigh-Taylor不稳定性混合实验研究[J].高压物理学报, 2012, 26(6):688-690. http://www.gywlxb.cn/CN/abstract/abstract1523.shtml

    LIU J H, TAN D W, ZHANG X, et al.Experimental investigation of mixing at tilted interface induced by Rayleigh-Taylor instability[J].Chinese Journal of High Pressure Physics, 2012, 26(6):688-690. http://www.gywlxb.cn/CN/abstract/abstract1523.shtml
    [18] YOUNGS D L.Experimental investigation of turbulent mixing by Rayleigh-Taylor instability[J].Physica D:Nonlinear Phenomena, 1991, 3(12):586-589. http://linkinghub.elsevier.com/retrieve/pii/016727898490513X
  • 加载中
图(9) / 表(1)
计量
  • 文章访问数:  7334
  • HTML全文浏览量:  3145
  • PDF下载量:  178
出版历程
  • 收稿日期:  2018-01-31
  • 修回日期:  2018-03-28

目录

/

返回文章
返回