• ISSN 0258-2724
  • CN 51-1277/U
  • EI Compendex
  • Scopus 收录
  • 全国中文核心期刊
  • 中国科技论文统计源期刊
  • 中国科学引文数据库来源期刊

基于离散单元法的铁路板结道床动力特性分析

崔旭浩 肖宏

高睿, 石知政, 刘洋泽鹏, 陈静, 张荣隆. 土工格栅对受污道砟直剪特性影响的试验研究[J]. 西南交通大学学报, 2021, 56(6): 1185-1191. doi: 10.3969/j.issn.0258-2724.20191142
引用本文: 崔旭浩, 肖宏. 基于离散单元法的铁路板结道床动力特性分析[J]. 西南交通大学学报, 2021, 56(6): 1197-1204. doi: 10.3969/j.issn.0258-2724.20200113
GAO Rui, SHI Zhizheng, LIU Yangzepeng, CHEN Jing, ZHANG Ronglong. Experimental Study on Effect of Geogrid on Direct Shear Behavior of Contaminated Ballast[J]. Journal of Southwest Jiaotong University, 2021, 56(6): 1185-1191. doi: 10.3969/j.issn.0258-2724.20191142
Citation: CUI Xuhao, XIAO Hong. Dynamic Characteristics Analysis of Hardening Railway Ballast Bed Based on Discrete Element Method[J]. Journal of Southwest Jiaotong University, 2021, 56(6): 1197-1204. doi: 10.3969/j.issn.0258-2724.20200113

基于离散单元法的铁路板结道床动力特性分析

doi: 10.3969/j.issn.0258-2724.20200113
基金项目: 中央高校基本科研业务费专项资金(2019YJS112)
详细信息
    作者简介:

    崔旭浩(1994—),男,博士研究生,研究方向为轨道工程,E-mail:xhcui01@126.com

    通讯作者:

    肖宏(1978—),男,教授,研究方向为轨道工程与工务管理, E-mail:xiaoh@bjtu.edu.cn

  • 中图分类号: U213.7

Dynamic Characteristics Analysis of Hardening Railway Ballast Bed Based on Discrete Element Method

  • 摘要:

    为研究铁路有砟道床板结对其自身动力特性的影响,利用离散单元法并考虑道砟颗粒的真实外形建立了板结道床的仿真分析模型,分析了脏污板结及板结程度对有砟道床动态响应的影响,并探讨了板结引起道床弹性损失、刚度变大的细观机理. 研究结果表明:道床板结会增大道床中道砟颗粒的振动水平,并且板结越严重影响越明显,板结可将道砟颗粒的振动加速度增大20%~30%;板结会加强列车荷载对道床的冲击作用,增大道砟颗粒之间的接触力容易引起道砟颗粒破碎劣化;板结会放大不同位置处道床工作性能的差异,增大道床刚度的不均匀性;道砟颗粒之间的脏污板结材质会抑制道砟颗粒的相对移动,可将道砟颗粒的滑动分数降低至正常的50%左右,减小列车荷载作用下道床整体的宏观变形,从而呈现出刚度增大的特性.

     

  • 道砟层是铁路轨道的重要组成部分,它通常由20~65 mm的中粗颗粒组成[1],具有减小枕木传递的荷载、提供横纵向阻力以及优化排水条件等功能[2]. 在风力、雨水等自然条件的作用下,铁路沿线的扬尘砂土等粉尘会散落在道床上,列车运行时的荷载会引起道砟的磨损和下层路基中小土石颗粒的上涌. 污染物的沉积会弱化道砟层的力学性能,影响铁路运营时基础的稳定性和安全性,因此必须花费大量资源对受污染的道砟进行清理和替换[3],所以研究污染物对道砟层的影响具有十分重要的工程和经济意义.

    当前,国内外学者已经对受污染道砟做了大量的研究. Akbar等[4]研究道砟颗粒级配对含砂污染物道砟剪切力学性能的影响,发现砂脏污会减小道砟的剪切强度;Kashani等[5]采用排水三轴试验研究了不同含水率条件下受污染道砟的变形特征;徐旸等[6]对多种道床脏污评估指标进行综合分析,认为道床抗剪性能的下降幅度与脏污程度有关;高亮等[7]通过直剪试验研究黄砂和煤灰脏污对道砟力学特性的影响效果,发现煤灰材质所引起的道床脏污对散体道床剪切力学性能的削弱作用大于黄砂脏污.

    为了控制道砟层的沉降和横向位移,许多学者开始研究采用土工格栅来加强道砟层的整体性. Hussaini等[8-9]通过对格栅加固道砟进行动力循环荷载试验,发现土工格栅可以限制道砟的纵横向变形;Chen等[10-11]研究了土工格栅在不同形状、位置等条件下的加固效果,确定土工格栅相对于道砟平均粒径的最优孔径比为1.4;Sweta等[12]通过改变道砟直剪试验的剪切速度,发现剪切速度可以在一定程度上影响土工格栅的加固效果.

    现有研究大多仅关注土工格栅的加固作用或污染物的弱化作用,但是污染物进入道砟颗粒骨架空隙中,不仅对道砟本身的承载机制产生影响,更会影响土工格栅对道砟散体的加固效果. 土工格栅和污染物的综合作用对道砟力学特性的影响规律尚不明晰,有待进一步研究. 因此,本研究利用自行研制的大型直剪仪,在考虑黏土污染的情况下,对土工格栅加固的道砟进行试验,分析不同污染程度下道砟混合体剪切强度、变形特性等变化规律,得出了不同污染情况下非线性强度及参数的拟合曲线,对受污染道砟维护和加固具有重要参考.

    本次试验采用的道砟材料是武汉地区附近未服役的一级铁路道砟,母岩为花岗岩. 道砟样本级配如图1所示,符合铁道行业标准TB/T 2140—2008中4.3.1章节关于道砟级配的要求;道砟试样的平均粒径d50为42 mm,不均匀系数Cu为1.59. 试样中加入的污染物是武汉地区Q4al+pl黏土,其参数见表1. 试验中所用的格栅(图2所示)为实际工程中常用到的双向土工格栅,其参数见表2.

    图  1  道砟试样级配
    Figure  1.  Particle size distribution of ballast samples
    表  1  试验材料参数表
    Table  1.  Physical properties of test materials
    材料堆积密度ρ/(kg•m−3比重 d 孔隙比 e 含水率ω/%液限/%塑限/%
    道砟14322.660.858
    黏土11782.701.29222.242.122.4
    下载: 导出CSV 
    | 显示表格
    图  2  试验土工格栅
    Figure  2.  Diagram of geogrid used in the experiment
    表  2  土工格栅参数表
    Table  2.  Physical and technical properties of geogrid
    材料 类型 孔径/mm 空隙率/% 极限抗拉强度/(kN•m−1
    聚丙烯 双向方形格栅 55 × 55 81 30
    下载: 导出CSV 
    | 显示表格

    本研究根据道砟特点,结合以往学者[13-15]的研究,采用自行研制的大型直剪仪,直剪仪的布置如图3. 剪切盒由固定的600 mm × 600 mm × 250 mm的上剪切盒和可移动的700 mm × 600 mm × 250 mm 的下剪切盒组成,其尺寸约为道砟最大粒径的9倍,可以忽略尺寸效应的影响[16]. 法向压力由固定于反力架横梁上法向的15 t高精度千斤顶提供,水平剪切力由抵在底座反力板上同型号水平向的千斤顶提供. 千斤顶通过与其相连的液压伺服控制设备实时调控压力,以保证法向压力稳定和水平向的匀速加载. 在千斤顶顶头底部安装压力传感器,以实时监测记录法向压力和水平剪切力;同时,由连接在静态电阻应变仪上的法向位移传感器和水平位移传感器实时监测记录法向位移和水平位移.

    图  3  试验装置
    Figure  3.  Diagram of direct shear test apparatus

    为研究不同污染程度对土工格栅加固道砟的影响,试样中要掺入不同量的黏土污染物,所选用的污染指标为Indraratna等[17]提出的VCI (void contamination index, VVCI):

    VVCI=1+efeb×Gs,bGs,f×MfMb×100%,
    (1)

    式中:ef为污染物的孔隙率;eb为道砟的孔隙率;Gs,b为道砟的比重;Gs,f为污染物的比重;Mf为污染物的干质量;Mb为道砟的干质量.

    本次试验分别对表3所示试样进行直剪试验.

    表  3  试验方案
    Table  3.  Experiment schemes
    试验组别试样编号VCI/%法向压力/kPa
    有格栅1 号015、35、55、75
    2 号20
    3 号40
    无格栅4 号015、35、55、75
    5 号20
    6 号40
    下载: 导出CSV 
    | 显示表格

    试验前,对道砟材料进行清洗、晒干、筛分. 在配制试样时,根据VCI指标的计算结果称取相应质量的道砟和黏土进行混合,控制道砟级配和所有试样的总质量保持一致. 将混合均匀的材料分4层装入剪切盒中,每层装填完成后,振动夯实至相同指定高度,使试样密实并保证试样初始状态的一致性. 在装填至上下剪切盒之间时,铺设双向土工格栅,格栅挂在固定于下剪切盒上的螺纹钢筋上来模拟非位移边界[18]. 整个试样装填完成后,用约100 kPa的法向压力进行预压约1 min. 试验开始后,先对试样施加法向压力,待法向压力稳定后,进行水平向加载,使下剪切盒保持3 mm/min的剪切速度. 当下剪切盒达到60 mm的剪切位移(即样本10%的剪切应变)时停止加载,加载过程中每达到1.5 mm的水平位移做记录点,记录法向应力、剪切应力和法向位移等数据.

    所有道砟试样直剪得到的剪切强度与水平位移关系曲线如图4所示,各曲线峰值剪切强度如表4,图中:σn为法向压力. 对比未加土工格栅不同污染程度下的剪切结果发现,在4种法向压力下,道砟试样的剪切强度峰值都随污染程度的增加而呈现下降趋势;且在相同法向压力下,VCI大的试样剪切强度上升趋势更为平缓,说明黏土污染不利于道砟试样剪切强度的形成. 其原因是试样的剪切强度主要来源于道砟颗粒间的咬合力和摩擦力,而黏土污染物填充了道砟颗粒之间的空隙,起到了对道砟颗粒的润滑作用,减少了道砟颗粒的咬合,从而降低了试样的剪切强度. 在各污染程度下,土工格栅加固均可有效增大道砟试样的剪切强度. 1号、2号、3号试样在4种法向压力下的平均峰值剪切强度分别比4号、5号、6号试样提高了17%,24%和22%. 对比1号试样与2号试样、4号试样与5号试样的结果发现,黏土污染使未铺设土工格栅试样的峰值剪切强度下降了13%,而采用了土工格栅加固的试样的峰值剪切强度只下降了8%. 说明土工格栅可以减小黏土污染物对道砟试样剪切强度的削弱作用. 原因是土工格栅的肋与节点可以限制道砟的移动,提供类似道砟颗粒间咬合的作用,减弱由于黏土颗粒填充孔隙对道砟颗粒咬合及摩擦带来的不利影响,从而增加了试样的剪切强度.

    图  4  剪切强度与水平位移关系曲线
    Figure  4.  Shear stress-shear displacement curves of ballast samples
    表  4  峰值剪切强度数据
    Table  4.  Results of peak shear stress
    试样编号σn=15 kPaσn=35 kPaσn=55 kPaσn=75 kPa
    1 号80.28130.28155.44178.19
    2 号74.17115.28141.53169.03
    3 号58.1994.03128.89153.47
    4 号70.4294.10140.28164.93
    5 号54.4484.44136.04149.72
    6 号46.8167.64111.39141.94
    下载: 导出CSV 
    | 显示表格

    由于道砟散体的特性,线性的摩尔库伦强度准则不能很好的表征出试样剪切强度的变化规律[16]. 因此本文采用Indraratna等[19]提出的非线性强度准则对道砟强度进行分析,如式(2).

    τpσc=a(σnσc)b,
    (2)

    式中:τp为试样的最大剪切强度;σc为道砟母岩的单轴抗压强度;ab为无量纲的拟合参数,后文分为有格栅(a1b1)和无格栅(a2b2)两种.

    表4中数据采用式(2)拟合,结果如图5所示(σc=130 MPa). 可以看出:随着污染程度的增加,道砟剪切强度包络线向下移动,道砟剪切强度降低.

    图  5  不同污染程度下道砟强度包络图
    Figure  5.  Strength envelops of ballast at different VCIs

    图5的拟合结果可以看出:拟合参数ab随污染程度的改变而变化,且在试验涉及的污染范围内(VCI = 0~40%)与VCI有较高的相关性. 采用指数函数对参数ab与VCI进行拟合,结果如图6所示. 参数ab增加会降低试样的峰值剪切强度,从图6中可以看到:与无污染试样相比,添加黏土污染试样的参数均会增大,并且参数增加的速率随着VCI的增加而逐渐增大,从而说明污染物含量VCI的增加对道砟试样剪切强度形成的不利影响是逐渐加剧的. 根据图6的拟合结果,在实际工程中,可以通过对VCI的测量,对道砟的剪切强度做出初步的估计,对实际工程具有一定的参考意义. 例如,将通过图6拟合公式计算出的有格栅拟合参数a1b1代入式(2)中,利用得到的强度包络线与洁净道砟的强度包络线进行对比,可以判断土工格栅的加固效果. 而在VCI大于40%时,试样内黏土的含量较多,试样峰值剪切强度下降较大,此时道砟层需要维护更新,不再具有应用的价值.

    图  6  非线性强度准则参数拟合曲线
    Figure  6.  Fitting curves by nonlinear failure criterion

    图7显示了试样剪切位移和法向位移之间的关系曲线,可以看出不同VVCI条件下,所有试样都是先剪缩后剪胀,符合密实材料变形的典型特征. 通常认为最大剪胀量小的试样,颗粒间的咬合更加紧密,更有利于剪切强度的形成. 从图中可以观察到:相较未铺设土工格栅时的结果,土工格栅的加固作用可以减小试样的最大剪胀量,这和以往学者的研究相同[20-21],并且本文在2号试验时,抑制剪胀的平均效果达到了26.6%,此时土工格栅的加固效果最明显. 然而也可以观察到,如在75 kPa法向压力下,5号和6号的最大剪胀量相较于4号,分别下降了40%和54%,即黏土污染会减少试样的最大剪胀量,而且VVCI越大,试样剪缩时间持续得越长,剪胀量减小得越明显. 但是根据周跃峰等[22]的研究,该剪胀量的减小是由于填充于道砟颗粒间的欠密实黏土的剪缩引起的,对道砟的剪切强度并没有提升作用. 相反,根据相关学者的研究[7,18],占据了道砟空隙的黏土颗粒会减少道砟重新排列的空间,不利于道砟重新排列形成相互咬合,这也可以由前文黏土污染降低试样剪切强度的试验结果进行印证.

    图  7  法向位移与水平位移关系曲线
    Figure  7.  Vertical displacement versus horizontal displacement curves

    峰值剪胀角是水平位移和法向位移曲线中最大斜率与水平线的夹角,反映了材料剪切变形的特点. 图8为不同法向压力下,试样的剪胀角随着污染指标变化的关系图. 从图中可以看出:随着VVCI的增加,所有法向压力下试样峰值剪胀角都呈现出减小的趋势;1号、2号、3号试样的峰值剪胀角比4号、5号、6号试样约小0.7°~3.7°,说明土工格栅可以限制道砟颗粒和污染物的移动,起到了抑制试样剪胀的作用.

    图  8  峰值剪胀角与污染程度关系
    Figure  8.  Relationship between the peak friction angle and VCI under various normal stress

    本文利用自主研发的大型直剪仪,采用黏土作为污染物,对不同VVCI的污染指标、不同格栅加固状态下的道砟试样进行了一系列的直剪试验. 分析了直剪条件下,污染物及土工格栅相互作用对道砟体强度及变形的影响,得出以下结论:

    1) 随着黏土污染物含量的增加,道砟体的峰值剪切强度呈现下降的趋势,剪切强度曲线上升趋势更加平缓,削弱了道砟层的剪切力学性能.

    2) 在相同污染程度下,土工格栅加固可以增大试样的峰值剪切强度,减小试样的最大剪胀量和约0.7°~3.7° 的峰值剪胀角;不同污染程度下,土工格栅的加固效果不同,本研究中在VVCI为20%条件下,格栅加固的效果最为明显.

    3) 采用土工格栅加固和未用土工格栅的道砟剪切强度均表现出典型的非线性,土工格栅会影响非线性强度准则中的无量纲参数. 在VVCI为0~40%时,污染指标和无量纲参数具有较高的相关性. 可以利用该相关性,通过对污染指标测量,预估道床的峰值剪切强度,对实际工程中受污染道砟的维护工作具有一定的参考意义.

  • 图 1  道砟颗粒建模

    Figure 1.  Ballast particle modeling

    图 2  有砟轨道离散元整体模型(单位:m)

    Figure 2.  Discrete element model of ballast bed (unit:m)

    图 3  板结道床局部模型

    Figure 3.  Local model of hardening ballast bed

    图 4  道床振动加速度

    Figure 4.  Ballast bed vibration acceleration

    图 5  道砟颗粒振动加速度矢量图

    Figure 5.  Vector diagram of ballast particle vibration acceleration

    图 6  5号轨枕下方道砟颗粒振动加速度时程曲线

    Figure 6.  Time history curves of ballast particle vibration acceleration under No. 5 sleeper

    图 7  道砟颗粒振动加速度幅值

    Figure 7.  Peak vibration acceleration of ballast particles

    图 8  道砟颗粒接触力均值时程曲线

    Figure 8.  Time history curves of mean contact force of ballast particles

    图 9  道砟颗粒平均接触力峰值

    Figure 9.  Peak of average contact force of ballast particles

    图 10  道砟颗粒滑动分数

    Figure 10.  Slip fraction of ballast particles

    图 11  单一时步中道砟颗粒切向滑动位移

    Figure 11.  Sliding displacement of ballast particles in a single time step

    图 12  道砟颗粒运动示意

    Figure 12.  Schematic diagram of the ballast particles movements

    表  1  计算工况

    Table  1.   Calculation conditions

    工况编号道床板结程度/%脏污颗粒数量/粒
    1洁净道床00
    2板结道床2517605
    3板结道床5035213
    4板结道床10070427
    下载: 导出CSV

    表  2  离散元模型计算参数

    Table  2.   Parameters of DE model

    参数取值
    道砟颗粒法向刚度/(N•m−15.0 × 108
    道砟颗粒切向刚度/(N•m−15.0 × 108
    道砟颗粒摩擦系数0.55
    道砟颗粒密度/(kg•m−32 500
    脏污颗粒密度/(kg•m−32500
    钢轨密度/(kg•m−31029.5
    轨枕密度/(kg•m−33025.8
    脏污颗粒黏结强度/MPa0.5
    黏结半径比例系数1.0
    下载: 导出CSV
  • [1] 王青波. 大秦线人工更换桥梁道砟施工技术[J]. 铁道建筑技术,2015(2): 66-69. doi: 10.3969/j.issn.1009-4539.2015.02.016

    WANG Qingbo. Construction technology of artificial replacement of bridge ballast for datong-qinhuangdao railway[J]. Railway Construction Technology, 2015(2): 66-69. doi: 10.3969/j.issn.1009-4539.2015.02.016
    [2] 郭善辛. 道床板结的形成及其分析[J]. 铁道工程学报,1985(2): 156-164.

    GUO Shanxin. Formation and analysis of ballast hardening[J]. Journal of Railway Engineering Society, 1985(2): 156-164.
    [3] 赵成江,王晓荣,夏琼,等. 不同含泥量与含水率下铁路板结道床的冻胀特性[J]. 铁道工程学报,2011,28(9): 36-39,96. doi: 10.3969/j.issn.1006-2106.2011.09.007

    ZHAO Chengjiang, WANG Xiaorong, XIA Qiong, et al. Frost heaving behaviors of hardened ballast beds under different content of clay and water[J]. Journal of Railway Engineering Society, 2011, 28(9): 36-39,96. doi: 10.3969/j.issn.1006-2106.2011.09.007
    [4] FORTUNATO E, PINELO A, FERNANDES M M. Characterization of the fouled ballast layer in the substructure of a 19th century railway track under renewal[J]. Soils and Foundations, 2010, 50(1): 55-62. doi: 10.3208/sandf.50.55
    [5] TANG A M, CUI Y J, DUPLA J C, et al. Mechanical characterisation of the fouled ballast in ancient railway track substructure by large-scale triaxial tests[J]. Soils and Foundations, 2012, 52(3): 511-523. doi: 10.1016/j.sandf.2012.05.009
    [6] PAIVA C, FERREIRA M, FERREIRA A. Ballast drainage in Brazilian railway infrastructures[J]. Construction and Building Materials, 2015, 92: 58-63. doi: 10.1016/j.conbuildmat.2014.06.006
    [7] ISHIKAWA T, FUKU S, NAKAMURA T, et al. Influence of water content on shear behavior of unsaturated fouled ballast[J]. Procedia Engineering, 2016, 143: 268-275. doi: 10.1016/j.proeng.2016.06.034
    [8] 金寿延,谭复兴,王午生. 道床板结对轨道振动的影响[J]. 同济大学学报(自然科学版),2001,29(10): 1162-1168. doi: 10.3321/j.issn:0253-374X.2001.10.006

    JIN Shouyan, TAN Fuxing, WANG Wusheng. Researching for the influence of cement hardening on the rail track vibration[J]. Journal of Tongji University (Natural Science), 2001, 29(10): 1162-1168. doi: 10.3321/j.issn:0253-374X.2001.10.006
    [9] 赵宁. TLV加载参数优化及轨道结构动态特征识别方法研究[D]. 北京: 北京交通大学, 2015.
    [10] 肖军华,张德,王延海,等. 基于DEM-FDM耦合的普通铁路碎石道床-土质基床界面接触应力分析[J]. 工程力学,2018,35(9): 170-179. doi: 10.6052/j.issn.1000-4750.2017.06.0422

    XIAO Junhua, ZHANG De, WANG Yanhai, et al. Study on interface stress between ballast and subgrade for traditional railway based on coupled DEM-FDM.[J]. Engineering Mechanics, 2018, 35(9): 170-179. doi: 10.6052/j.issn.1000-4750.2017.06.0422
    [11] ZHANG X, ZHAO C, ZHAI W. Dynamic behavior analysis of high-speed railway ballast under moving vehicle loads using discrete element method[J]. International Journal of Geomechanics, 2017, 17(7): 04016157.1-04016157.14.
    [12] 井国庆,强伟乐,常锦秀,等. 针片状指数对道砟直剪力学特性的影响[J]. 西南交通大学学报,2020,55(4): 688-694.

    JING Guoqing, QIANG Weile, CHANG Jinxiu, et al. Effect of flakiness-elongation index on shear behavior of railway ballast[J]. Journal of Southwest Jiaotong University, 2020, 55(4): 688-694.
    [13] 中华人民共和国铁道部. 铁路碎石道砟: TB/T 2140—2018[S]. 北京: 中国铁道出版社, 2018.
    [14] 付龙龙,宫全美,周顺华,等. 列车荷载作用下有砟轨道轨面沉降与路基不均匀沉降间的相关关系[J]. 振动与冲击,2013,32(14): 23-28. doi: 10.3969/j.issn.1000-3835.2013.14.004

    FU Longlong, GONG Quanmei, ZHOU Shunhua, et al. Relationships between ballasted track settlement and differential subgrade settlement under train load[J]. Journal of Vibration and Shock, 2013, 32(14): 23-28. doi: 10.3969/j.issn.1000-3835.2013.14.004
    [15] Itasca Consulting Group. Inc. PFC2D (particle flow code in 2 dimensions) theory and background[M]. Minneapolis : Itasca Consulting Group, 2002.
    [16] 伍超. 铁路板结道床捣固作业机理研究[D]. 成都: 西南交通大学, 2016.
    [17] 徐旸,高亮,井国庆,等. 脏污对道床剪切性能影响及评估指标的离散元分析[J]. 工程力学,2015,32(8): 96-102. doi: 10.6052/j.issn.1000-4750.2014.01.0051

    XU Yang, GAO Liang, JING Guoqing, et al. Shear behavior analysis of fouled railroad ballast by dem and its evaluation index[J]. Engineering Mechanics, 2015, 32(8): 96-102. doi: 10.6052/j.issn.1000-4750.2014.01.0051
    [18] 国家铁路局. 轨道检测道床应变、振动测试方法: TB/T 2488—2017[S]. 北京: 中国铁道出版社, 2017.
    [19] 郑 建. 增加有砟轨道弹性的措施及其对轨道力学性能影响[D]. 成都: 西南交通大学, 2017.
  • 加载中
图(12) / 表(2)
计量
  • 文章访问数:  595
  • HTML全文浏览量:  279
  • PDF下载量:  23
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-03-24
  • 修回日期:  2020-06-03
  • 网络出版日期:  2020-07-07
  • 刊出日期:  2020-07-07

目录

/

返回文章
返回