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

数字孪生演进模型及其在智能制造中的应用

江海凡 丁国富 肖通 樊孟杰 付建林 张剑

宋启峰, 陈光雄, 董丙杰, 张峻才, 冯晓航. 地铁梯形轨枕轨道钢轨波磨成因研究[J]. 西南交通大学学报. doi: 10.3969/j.issn.0258-2724.20230573
引用本文: 江海凡, 丁国富, 肖通, 樊孟杰, 付建林, 张剑. 数字孪生演进模型及其在智能制造中的应用[J]. 西南交通大学学报, 2022, 57(6): 1386-1394. doi: 10.3969/j.issn.0258-2724.20210202
SONG Qifeng, CHEN Guangxiong, DONG Bingjie, ZHANG Juncai, FENG Xiaohang. Cause of Rail Corrugation on Ladder Sleeper Track[J]. Journal of Southwest Jiaotong University. doi: 10.3969/j.issn.0258-2724.20230573
Citation: JIANG Haifan, DING Guofu, XIAO Tong, FAN Mengjie, FU Jianlin, ZHANG Jian. Digital Twin Evolution Model and Its Applications in Intelligent Manufacturing[J]. Journal of Southwest Jiaotong University, 2022, 57(6): 1386-1394. doi: 10.3969/j.issn.0258-2724.20210202

数字孪生演进模型及其在智能制造中的应用

doi: 10.3969/j.issn.0258-2724.20210202
基金项目: 国家重点研发计划(2020YFB1708001)
详细信息
    作者简介:

    江海凡(1992—),男,博士研究生,研究方向为数字化车间建模与仿真,数字孪生,E-mail:JiangHaifan@my.swjtu.edu.cn

    通讯作者:

    张剑(1972—),女,教授,研究方向为智能制造,生产排程与调度,E-mail:jerrysmail@263.net

  • 中图分类号: TP301.6

Digital Twin Evolution Model and Its Applications in Intelligent Manufacturing

  • 摘要:

    数字孪生作为实现智能制造信息物理融合的关键使能技术受到广泛关注,而如何构建数字孪生模型成为当前研究的热点. 目前,数字孪生模型多聚焦于概念抽象或具体工程应用,而较少从构建方法和过程的角度考虑如何分阶段、有步骤地构建和应用数字孪生模型. 因此,本文提出数字孪生演进模型的概念,将数字孪生构建与应用过程分为数字模型、数字投影以及数字孪生3个演进阶段,给出各演进阶段的应用方法、关键技术与工具平台,并探讨数字孪生演进模型在智能装备、智能生产、智能运维方面的典型应用. 研究结果表明:所提模型为数字孪生在智能制造中的分步实施提供了可行的技术路线与有益的应用参考.

     

  • 钢轨波磨是铁路行业尚未解决的难题,在地铁线路上尤为严重. 任何一条地铁线路上,只要曲线轨道半径低于350 m,则该轨道内轨上出现波磨的几率接近100%;然而在大半径曲线轨道(R≥650~800 m)和直线轨道上,很少出现钢轨波磨:这是地铁线路上钢轨波磨的普遍规律[1-2]. 当列车通过波磨轨道时,轮轨界面上的剧烈振动加剧了车辆和轨道部件的疲劳断裂,甚至危害行车安全[3-4].

    目前钢轨波磨的形成机理可分为3类:1) 轮轨表面初始不平顺激励下,轮轨界面的瞬时动力作用引起摩擦功波动导致钢轨波磨[5-10];2) 轮轨系统的粘-滑振动引起了钢轨波磨[11-13];3) 饱和轮轨蠕滑力激发了轮轨摩擦耦合系统的自激振动导致钢轨波磨[14-19]. 第一类理论模型为它激振动模型,该理论可以较好地解释地铁线路上的部分波磨现象,但不足以解释地铁线路上波磨主要出现在小半径曲线内轨,在大半径曲线和直线轨道上波磨发生率较低这一普遍现象;第二类理论模型预测得出的波磨频率约为20~80 Hz[12-13],远小于地铁线路上绝大多数波磨的通过频率;第三类波磨理论由陈光雄教授提出,该理论认为当车辆通过小半径曲线轨道时,导向轮对上的轮轨蠕滑力饱和,引发了轮轨系统的摩擦自激振动,致使钢轨表面出现波磨,该理论可以合理地解释地铁线路上规律性的波磨现象,正逐渐被学者们认可[20-23].

    梯形轨枕因其具有良好的减振降噪效果被铺设在减振要求高的地段[24]. 然而,在小半径曲线梯形轨枕轨道上出现了严重的钢轨波磨[25],目前对于梯形轨枕轨道波磨成因的研究还不充分. 本文基于轮轨系统摩擦自激振动导致钢轨波磨的观点对小半径曲线上的梯形轨枕轨道波磨展开研究,建立了导向轮对-梯形轨枕轨道系统有限元模型. 该模型采用实体单元对扣件系统进行建模,考虑了弹条的扣压作用以及轨距挡块的横向限位作用,建模方式更加符合实际工况. 应用复特征值分析和瞬时动态分析分别求解了该摩擦耦合系统的运动稳定性和时域动态响应,揭示了小半径曲线上梯形轨枕轨道波磨的成因.

    当车辆通过小半径曲线梯形轨枕轨道时,导向轮对与钢轨的接触位置如图1所示. 外轮的轮缘与高轨的轨距面接触,内轮的踏面与低轨的轨顶面接触. 图1中:δLδR分别为外轮、内轮与钢轨间的接触角,FLNL分别为外轮与高轨接触界面上的横向蠕滑力和法向接触力,FRNR分别为内轮与低轨接触界面上的横向蠕滑力和法向接触力,FSVLFSVR为作用在导向轮对左、右轴端的垂向悬挂力,FSLLFSLR为作用在导向轮对左、右轴端的横向悬挂力,KC为弹条刚度,KDLCDL分别为横向缓冲垫的刚度和阻尼,KDVCDV分别为垂向减振垫的刚度和阻尼.

    图  1  导向轮对-梯形轨枕轨道系统接触模型简图
    Figure  1.  Contact model of leading wheelset‒ladder sleeper track system

    扣件系统由轨下垫板、弹条、轨距挡块、铁垫板等组成,轨距挡块约束着钢轨的横向运动,弹条产生的扣压力将钢轨紧紧地压在轨道上,钢轨的纵向运动通过钢轨与轨下垫板间的摩擦力来约束,螺旋道钉将铁垫板固定在梯形轨枕上. 为了简化模型,本文仅对扣件系统中的轨下垫板和铁垫板进行实体单元建模.

    依托于图1所示的轮轨系统接触模型,使用ABAQUS有限元分析软件建立了导向轮对-梯形轨枕轨道系统有限元模型,如图2所示. 车轮廓形为LM磨耗型踏面,轮对滚动圆半径为0.42 m;钢轨廓形为60 kg/m标准钢轨截面,模型中单根钢轨的长度为36 m;轨下垫板离散地分布在钢轨底部,相邻轨下垫板/铁垫板间的距离为0.625 m;梯形轨枕由横向钢管和预应力混凝土纵梁组合而成,相邻横向钢管间的距离为2.5 m;导向轮对、钢轨、轨下垫板、铁垫板和梯形轨枕通过实体单元建模,整个有限元模型采用C3D8I实体单元划分,模型中共包含422968个单元,606942个节点. 轮轨间的接触以及钢轨和轨下垫板间的接触采用面-面接触算法和有限滑移准则进行建模. 轮轨间的摩擦系数为0.40,钢轨与轨下垫板间的摩擦系数为0.75[26]. 通过垂向压缩弹簧来模拟弹条的扣压作用,通过“Equation”约束模拟轨距挡块的横向限位作用,通过“Tie”约束模拟轨下垫板和铁垫板以及铁垫板和梯形轨枕的连接关系. 模型中采用接地弹簧和阻尼单元模拟侧向缓冲垫和垂向减振垫对梯形轨枕的缓冲减振作用,相邻缓冲减振垫间的距离为1.25 m. 缓冲减振垫的刚度为2 × 107 N/m,阻尼为7 × 104 N·s/m[27]. 各部件的材料参数见表1.

    图  2  导向轮对-梯形轨枕轨道系统有限元模型
    Figure  2.  Finite element model of leading wheelset‒ladder sleeper track system
    表  1  材料参数
    Table  1.  Material parameters
    部件 密度/
    (×103 kg·m−3
    弹性模量/MPa 泊松比
    钢轨/轮对 7.80 210000 0.30
    轨下垫板 1.19 12 0.45
    铁垫板 7.80 173000 0.30
    横向钢管 7.80 210000 0.30
    混凝土纵梁 2.50 36000 0.20
    下载: 导出CSV 
    | 显示表格

    采用复特征值分析来求解导向轮对-梯形轨枕轨道系统的运动稳定性. 在小半径曲线轨道段,导向轮对上的轮轨蠕滑力达到饱和,约等于滑动摩擦力. 在此情形下,轮轨摩擦耦合系统的运动方程为

    M¨u+C˙u+Ku=0, (1)

    式中:MCK分别为系统的质量、阻尼和刚度矩阵,u为节点位移向量.

    式(1)的通解为[28]

    u(t)=zi=1φieλit=zi=1φie(αi+jωi)t, (2)

    式中:z为拟提取的特征值阶数;${\lambda _i} = {\alpha _i} + {\mathrm{j}}{\omega _i}$,为第i阶特征值;${{\boldsymbol{\varphi}} _i} $为特征向量;t为时间.

    由式(2)可以看出,当特征值实部为正时,节点位移随时间指数式增长,这意味着轮轨系统发生了不稳定振动. 通常采用等效阻尼比$ \xi $来判断系统发生自激振动的趋势,如式(3).

    ξ=2αi|ωi|. (3)

    当$ \xi <$0时,轮轨摩擦耦合系统发生了自激振动. 而且,$ \xi $越小,轮轨系统发生自激振动的趋势越强,钢轨波磨更易发生.

    复特征值分析是一种频域分析方法,而瞬时动态分析是一种时域非线性分析方法,可以获取摩擦自激振动发生时系统的动态响应. 本文使用ABAQUS/Standard隐式求解器计算了轮轨系统的动态响应. 该求解器采用Newmark方法计算每个增量步结束时刻系统的平衡方程,进而获取节点的位移和速度[29],如式(4)、(5).

    u|t+Δt=u|t+Δt˙u|t+Δt2((12β)¨u|t+β¨u|t+Δt), (4)
    ˙u|t+Δt=˙u|t+Δt((1γ)¨u|t+γ¨u|t+Δt), (5)

    式中:Δt为积分步长,β = (1–α2)/4,γ = 0.5–α,−0.5≤α≤0.

    复特征值分析的步骤为:

    步骤1 在导向轮对轴端施加垂向悬挂力和横向悬挂力;

    步骤2 引入轮轨间的摩擦耦合作用;

    步骤3 提取轮轨系统的固有频率;

    步骤4 进行复特征值分析,提取轮轨系统的摩擦自激振动.

    通过定义轮对的横向滑动速度引入导向轮对和钢轨间的摩擦耦合作用. 应用复特征值分析提取了50~1200 Hz内导向轮对-梯形轨枕轨道系统的摩擦自激振动,如图3所示. 可以看到,在该频率范围内系统共存在7个自激振动,对应频率分别为150.35、163.07、422.16、447.81、471.05、493.97 Hz和519.94 Hz.

    图  3  导向轮对-梯形轨枕轨道系统的自激振动频率分布
    Figure  3.  Distribution of self-excited vibration frequencies of leading wheelset‒ladder sleeper track system

    图4为7个自激振动的振动模态,自激振动主要发生在低轨和内轮上,与钢轨波磨发生在小半径曲线梯形轨枕轨道内轨上相吻合. 根据等效阻尼比越小,自激振动发生趋势越强的原则,可知频率为150.35 Hz和493.97 Hz的自激振动更易发生.

    图  4  导向轮对-梯形轨枕轨道系统的自激振动振型
    Figure  4.  Mode shapes of self-excited vibrations of leading wheelset‒ladder sleeper track system

    图5为曲线半径为350 m的梯形轨枕轨道上的内轨波磨,其波长$\lambda $为60~80 mm. 列车通过该轨道的平均速度$v $为37 km/h,由$\lambda = v/f$可得波磨通过频率为130~180 Hz[27]. 综上可得:当车辆通过小半径曲线段梯形轨枕轨道时,发生在导向轮对内侧轮轨界面上频率为150.35 Hz的自激振动致使内轨上出现了波长为60~80 mm的波磨.

    图  5  梯形轨枕轨道内轨波磨[27]
    Figure  5.  Inner rail corrugation on ladder sleeper track[27]

    应用瞬时动态分析求解轮轨系统的动态响应分2个步骤完成:第一步,在导向轮对的轴端施加垂向悬挂力和横向悬挂力;第二步,设置轮对运行速度为40 km/h,运行时间为0.15 s,采样间隔为5 × 10–5 s. 由于未施加任何外部激励且钢轨表面无初始不平顺,因此,计算所得的轮轨系统的响应为自激振动响应. 图6所示为轮轨界面上的法向接触力的时间历程曲线. 从图中可以看出,内侧轮轨界面上的法向接触力波动幅值要大于外侧轮轨界面上的,其中外轮与高轨间法向接触力的波动幅值约为20 kN,而内轮与低轨间法向接触力的波动幅值约为30 kN.

    图  6  轮轨法向接触力时间历程曲线
    Figure  6.  Time history of wheel-rail normal contact force

    图7为轮轨法向接触力的功率谱密度分析结果. 可以看到,内外侧轮轨界面上法向接触力的振动主频为160 Hz. 由$\lambda = v/f$可得波长约为69 mm,该预测结果与线路上梯形轨枕轨道波磨的波长(60~80 mm)接近.

    图  7  轮轨法向接触力的功率谱密度分析
    Figure  7.  Power spectral density analysis of wheel-rail normal contact force

    不同侧向缓冲垫阻尼CDL,轮轨系统自激振动频率的分布如图8所示. 从图中可以看出:轮轨系统自激振动的数量不随侧向缓冲垫阻尼变化,轮轨系统产生了7个自激振动;侧向缓冲垫阻尼对400~600 Hz内的自激振动的影响很小,而对150~200 Hz内的2个自激振动有较大影响,2个自激振动的等效阻尼比随侧向缓冲垫阻尼值的增大而增大,这说明系统发生自激振动的趋势减弱. 因此,增大侧向缓冲垫阻尼可以在一定程度上抑制梯形轨枕轨道上的内轨波磨.

    图  8  不同侧向缓冲垫阻尼下轮轨系统自激振动频率的分布
    Figure  8.  Distribution of self-excited vibration frequencies of wheel-rail system under different damping levels of lateral cushioning pad

    图9所示为不同垂向减振垫刚度KDV,轮轨系统自激振动频率的分布. 从图9可以看出,当垂向KDV=10 kN/mm时,轮轨系统有6个自激振动,最小等效阻尼比为–0.00789;当KDV=20 kN/mm时,轮轨系统有7个自激振动,最小等效阻尼比为−0.00506;当KDV=30 kN/mm时,轮轨系统有6个自激振动,最小等效阻尼比为−0.00679;当KDV=40 kN/mm时,轮轨系统有6个自激振动,最小等效阻尼比为−0.00798. 由以上分析可知,当KDV=20 kN/mm时,最小等效阻尼比高于其余3种工况,这意味着系统发生自激振动的趋势减弱. 过小或过大的垂向减振垫刚度均不利于抑制梯形轨枕轨道上的内轨波磨.

    图  9  不同垂向减振垫刚度下轮轨系统自激振动频率的分布
    Figure  9.  Distribution of self-excited vibration frequencies of wheel-rail system under different stiffness levels of vertical damping pad

    图10所示为4种不同结构的梯形轨枕. 每块轨枕上横向钢管的数量N=2根时,对应结构的梯形轨枕相邻横向钢管间的距离为5.000 m;N=3根时,对应结构的梯形轨枕相邻横向钢管间的距离为2.500 m;N=5根时,对应结构的梯形轨枕相邻横向钢管间的距离为1.250 m;N=9根时,对应结构的梯形轨枕相邻横向钢管间的距离为0.625 m. 需要说明的是,横向钢管间距的设计值为2.500 m.

    图  10  不同结构的梯形轨枕
    Figure  10.  Ladder sleepers with different structures

    图11所示为4个不同结构的导向轮对-梯形轨枕轨道系统自激振动频率的分布. 从图中可以看出:梯形轨枕结构对轮轨系统的摩擦自激振动有显著影响,并非每块梯形轨枕上的横向钢管数量越多,轮轨系统就越稳定. 如蓝色椭圆圈内所示,当每块梯形轨枕上安装9根横向钢管时,轮轨系统产生自激振动的等效阻尼比反而更小,轮轨系统更不稳定,钢轨波磨更易出现. 红色椭圆圈内所示为低频段的2个自激振动,当每块梯形轨枕上安装2根横向钢管时,等效阻尼比更低,轮轨系统更易趋于不稳定. 当每块梯形轨枕上安装5根横向钢管时,等效阻尼比最大,轮轨系统发生自激振动的趋势减弱. 因此,在地铁小半径曲线地段,铺设横向钢管间距为1.250 m的梯形轨枕可以一定程度抑制梯形轨枕轨道波磨.

    图  11  不同结构的轮轨系统的自激振动频率分布
    Figure  11.  Distribution of self-excited vibration frequencies of wheel-rail systems with different structures

    本文研究了地铁小半径曲线梯形轨枕轨道上的内轨波磨,建立了导向轮对-梯形轨枕轨道系统有限元模型,模型中采用实体单元对扣件系统进行建模,应用复特征值分析和瞬时动态分析分别求解了轮轨系统的运动稳定性和时域动态响应. 研究了缓冲减振垫参数及梯形轨枕结构对轮轨系统摩擦自激振动的影响,得到以下结论:

    1) 饱和轮轨蠕滑力引起的频率为150 Hz的轮轨系统共振振动是地铁线路小半径曲线段梯形轨枕轨道上波长为60~80 mm内轨波磨的成因.

    2) 增大侧向缓冲垫阻尼可以一定程度上抑制梯形轨枕轨道上的内轨波磨,过小或过大的垂向减振垫刚度均不利于抑制钢轨波磨.

    3) 改变单块梯形轨枕上横向钢管的数量实际上是改变了导向轮对–梯形轨枕轨道系统的刚度矩阵,进而影响了轮轨系统发生自激振动的趋势以及钢轨波磨的发展趋势. 铺设横向钢管间距为1.250 m结构的梯形轨枕可以抑制梯形轨枕轨道上的钢轨波磨.

  • 图 1  数字孪生演进模型及其内涵

    Figure 1.  Digital twin evolution model and its connotations

    图 2  数字孪生演进模型应用方法

    Figure 2.  Application method of digital twin evolution model

    图 3  数字孪生不同应用阶段使能技术与工具/平台

    Figure 3.  Enabling technologies and tools/platforms for digital twin at different stages of application

    图 4  数字孪生演进模型在智能设备中的应用

    ①:虚拟机床;②: 控制指令和数控程序;③: 验证和优化后的控制指令和数控程序;④:实时数据、仿真数据归集形成历史数据集;⑤ :数据决策与控制.

    Figure 4.  Application of digital twin evolution model in smart devices

    图 5  数字孪生演进模型在智能生产中的应用

    SCADA:数据采集与监视控制系统;HDFC:hadoop分布式文件系统;ERP:企业资源计划;MES:制造执行系统;①: 虚拟车间;②: 生产指令和生产计划/调度方案;③: 验证和优化后的生产指令和生产计划/调度方案;④: 实时数据、仿真数据归集形成历史数据集;⑤:数据决策与控制.

    Figure 5.  Application of digital twin evolution model in smart production

    图 6  数字孪生演进模型在智能运维中的应用

    ①:虚拟列车;②: 车辆运行采集数据和环境数据;③: 仿真评估后的车辆及其部件状态性能(提前预判) ;④:实时数据、仿真数据归集形成历史数据集;⑤: 云边协同模式下的数据决策与反馈优化.

    Figure 6.  Application of digital twin evolution model in smart management and maintenance

  • [1] TAO F, ZHANG M. Digital twin shop-floor: a new shop-floor paradigm towards smart manufacturing[J]. IEEE Access, 2017, 5: 20418-20427. doi: 10.1109/ACCESS.2017.2756069
    [2] DING K, CHAN F T S, ZHANG X D, et al. Defining a digital twin-based cyber-physical production system for autonomous manufacturing in smart shop floors[J]. International Journal of Production Research, 2019, 57(20): 6315-6334. doi: 10.1080/00207543.2019.1566661
    [3] TAO F, ZHANG M, LIU Y S, et al. Digital twin driven prognostics and health management for complex equipment[J]. CIRP Annals, 2018, 67(1): 169-172. doi: 10.1016/j.cirp.2018.04.055
    [4] 陶飞,刘蔚然,张萌,等. 数字孪生五维模型及十大领域应用[J]. 计算机集成制造系统,2019,25(1): 1-18.

    TAO Fei, LIU Weiran, ZHANG Meng, et al. Five-dimension digital twin model and its ten applications[J]. Computer Integrated Manufacturing Systems, 2019, 25(1): 1-18.
    [5] BAO J S, GUO D S, LI J, et al. The modelling and operations for the digital twin in the context of manufacturing[J]. Enterprise Information Systems, 2019, 13(4): 534-556. doi: 10.1080/17517575.2018.1526324
    [6] CECIL J, ALBUHAMOOD S, CECIL-XAVIER A, et al. An advanced cyber physical framework for micro devices assembly[J]. IEEE Transactions on Systems, Man, and Cybernetics: Systems, 2019, 49(1): 92-106. doi: 10.1109/TSMC.2017.2733542
    [7] JIANG H F, QIN S F, FU J L, et al. How to model and implement connections between physical and virtual models for digital twin application[J]. Journal of Manufacturing Systems, 2021, 58: 36-51. doi: 10.1016/j.jmsy.2020.05.012
    [8] KRITZINGER W, KARNER M, TRAAR G, et al. Digital twin in manufacturing: a categorical literature review and classification[J]. IFAC-PapersOnLine, 2018, 51(11): 1016-1022. doi: 10.1016/j.ifacol.2018.08.474
    [9] 江海凡,丁国富,张剑. 数字孪生车间演化机理及运行机制[J]. 中国机械工程,2020,31(7): 824-832, 841. doi: 10.3969/j.issn.1004-132X.2020.07.008

    JIANG Haifan, DING Guofu, ZHANG Jian. Evolution and operation mechanism of digital twin shopfloors[J]. China Mechanical Engineering, 2020, 31(7): 824-832, 841. doi: 10.3969/j.issn.1004-132X.2020.07.008
    [10] QI Q L, TAO F, ZUO Y, et al. Digital twin service towards smart manufacturing[J]. Procedia CIRP, 2018, 72: 237-242. doi: 10.1016/j.procir.2018.03.103
    [11] 赵颖,侯俊杰,于成龙,等. 面向生产管控的工业大数据研究及应用[J]. 计算机科学,2019,46(增1): 45-51.

    ZHAO Ying, HOU Junjie, YU Chenglong, et al. Study and application of industrial big data in production management and control[J]. Computer Science, 2019, 46(S1): 45-51.
    [12] 陈建. 通用五轴数控加工仿真系统研发[D]. 成都: 西南交通大学, 2014.
    [13] 肖通,江海凡,丁国富,等. 五轴磨床数字孪生建模与监控研究[J]. 系统仿真学报,2021,33(12): 2880-2890.

    XIAO Tong, JIANG Haifan, DING Guofu, et al. Research on digital twin-based modeling and monitoring of five-axis grinder[J]. Journal of System Simulation, 2021, 33(12): 2880-2890.
    [14] 骆伟超. 基于Digital Twin的数控机床预测性维护关键技术研究[D]. 济南: 山东大学, 2020.
    [15] 丁国富, 江海凡, 罗樟圳, 等. 一种任务车间生产计划验证方法: CN110989527A[P]. 2020-04-10.
    [16] 罗樟圳,江海凡,付建林,等. 基于组合赋权的离散车间生产计划综合评价[J]. 系统仿真学报,2021,33(8): 1856-1865.

    LUO Zhangzhen, JIANG Haifan, FU Jianlin, et al. Combination weighting-based comprehensive evaluation for discrete workshop production plan[J]. Journal of System Simulation, 2021, 33(8): 1856-1865.
  • 加载中
图(6)
计量
  • 文章访问数:  807
  • HTML全文浏览量:  842
  • PDF下载量:  173
  • 被引次数: 0
出版历程
  • 收稿日期:  2021-03-23
  • 修回日期:  2021-06-21
  • 网络出版日期:  2022-10-20
  • 刊出日期:  2021-07-06

目录

/

返回文章
返回