基于直线电机驱动的插件机振动特性研究

胡意波, 高自成, 李立君, 邹洋, 罗义, 赵凯杰

胡意波, 高自成, 李立君, 邹洋, 罗义, 赵凯杰. 基于直线电机驱动的插件机振动特性研究[J]. 制造技术与机床, 2023, (12): 32-39. DOI: 10.19287/j.mtmt.1005-2402.2023.12.004
引用本文: 胡意波, 高自成, 李立君, 邹洋, 罗义, 赵凯杰. 基于直线电机驱动的插件机振动特性研究[J]. 制造技术与机床, 2023, (12): 32-39. DOI: 10.19287/j.mtmt.1005-2402.2023.12.004
HU Yibo, GAO Zicheng, LI Lijun, ZOU Yang, LUO Yi, ZHAO Kaijie. Research on vibration characteristics of plug-in machine driven by linear motor[J]. Manufacturing Technology & Machine Tool, 2023, (12): 32-39. DOI: 10.19287/j.mtmt.1005-2402.2023.12.004
Citation: HU Yibo, GAO Zicheng, LI Lijun, ZOU Yang, LUO Yi, ZHAO Kaijie. Research on vibration characteristics of plug-in machine driven by linear motor[J]. Manufacturing Technology & Machine Tool, 2023, (12): 32-39. DOI: 10.19287/j.mtmt.1005-2402.2023.12.004

基于直线电机驱动的插件机振动特性研究

基金项目: 湖南省重点研发计划项目(2021NK2023)
详细信息
    作者简介:

    胡意波,1998年生,男,硕士研究生,研究方向为智能技术及装备。E-mail:huyibo527679@163.com

    通讯作者:

    高自成,1968年生,男,硕士生导师,教授,研究方向为智能装备设计与开发、现代林业装备等。E-mail:gzc1968@126.com

  • 中图分类号: TH-39,TH113.1,TN403

Research on vibration characteristics of plug-in machine driven by linear motor

  • 摘要: 针对由直线电机引起的机械振动影响异形插件机精度的问题,对异形插件机的振动特性进行研究,并根据振动特性进行结构优化。首先对异形插件机上的直线电机进行推力和频率测试;其次对异形插件机的整机模型进行模态分析,根据模态分析结果对零部件进行结构优化,提升异形插件机的低阶固有频率;然后对优化后的整机模型进行谐响应分析,计算出取样点的振幅-频率响应数据;最后利用视觉相机对异形插件机进行重复定位精度测试实验。结果表明,优化后的异形插件机的前六阶固有频率有不同程度的提升,采样点在0~50 Hz的直线电机激励下,振动位移在0.03 mm左右。实验测试重复定位精度稳定在±0.03 mm内,验证了结构优化的有效性和仿真结果的可靠性。
    Abstract: Aiming at the problem that the mechanical vibration caused by the linear motor affects the accuracy of the special-shaped plug-in machine, the vibration characteristics of the special-shaped plug-in machine were studied, and the structure was optimized according to the vibration characteristics. Firstly, the thrust and frequency tests of the linear motor on the special-shaped plug-in machine were carried out. Secondly, the modal analysis of the whole machine model of the special-shaped plug-in machine was carried out, and the structural optimization of the parts was carried out according to the modal analysis results, so as to improve the low-order natural frequency of the special-shaped plug-in machine. Then, the harmonic response analysis of the optimized whole machine model was carried out, and the amplitude-frequency response data of the sampling point is calculated. Finally, the vision camera was used to test the repeatability of the special-shaped plug-in machine. The results show that the first six natural frequencies of the optimized special-shaped plug-in machine have different degrees of improvement, and the vibration displacement of the sampling point was about 0.03 mm under the excitation of the linear motor of 0~50 Hz. The experimental test was that the repeatability of positioning was stable within ±0.03 mm, which verifies the effectiveness of structural optimization and the reliability of simulation results.
  • 随着工业自动化、智能化的推进,“灯塔工厂”等一系列先进工厂对电子元器件组装生产线上贴插设备的要求越来越高。异形插件机是电子元器件组装生产线上不可或缺的自动化贴插设备,主要负责异形元器件的插装工作,其插装速度和插件精度影响着整个生产线的生产效率和质量[1]

    贴插设备在高速高加速度运动的状态下,容易产生机械振动,所以在贴插设备上高速和高精度是很难兼顾的。提高贴插设备精度的方法主要包括速度控制研究、图像处理以及机械结构研究:Hu G Q等[2]提出了一种基于S形速度规划、效率最优的实时运动控制方法,可以在提高运动速度的同时有效地保证精度。武洪恩等[3]提出了一种基于Hu矩和递进Hough变换的SOT元件识别算法,对视觉相机采集的图像进行处理,可以提高贴片的定位精度。Huang Y H等[4]提出了一种不仅能适应灰度线性和灰度非线性变化,而且支持尺度和视角变化的几何特征识别定位算法,可以有效提升定位精度。孙继斌等[5]根据模态分析结果对贴片机的横梁结构进行优化设计,提高横梁的抗振性。王莹[6]通过对高速贴片机的关键部件进行结构优化,提升了整机模型的固有频率。贴插设备结构动态特性是影响精度的主要因素,对零部件进行模态分析和结构优化,需要验证零部件之间是否会产生共振。

    选用高速高加速度的直线电机作为异形插件机的驱动部件,在直线电机的极限加速度运动状态下,对异形插件机整机模型进行振动特性研究,通过零部件的结构优化提升整机结构的固有频率,避开共振,提升精度。实验确定了异形插件机上安装的直线电机所产生的随机频率和推力范围,以此作为判断是否会发生共振的依据和谐响应分析的边界条件。利用模态分析和谐响应分析得到整机结构的固有频率、振型和振动幅值。重复定位精度的测试结果验证了优化结果和仿真结果的可靠性,确定了优化后的异形插件机不会发生共振现象。

    异形插件机的整机结构如图1所示。所研究的异形插件机是双横梁结构,前后都可以用供料器进行供料,横梁部分的一侧由固定在大理石梁上的直线电机驱动,另一侧通过导轨和滑块支撑在支撑梁上。插件头固定在横梁上,由横梁上的直线电机进行驱动。供料器和PCB板运输机构安装在底座上。大理石梁采用大理石制成,具有刚性好、硬度高、耐磨性强以及受温度影响变形小等特点[7-8]

    图  1  异形插件机的整机结构

    异形插件机的工作对象为异形元器件等轻质量物料,受静力作用影响比较小,因此目前设计的异形插件机的结构静强度和静刚度基本能满足要求。然而,异形插件机配置直线电机具有高速高加速度等特性,因此更需要考虑结构的动态特性。

    在传统的拱架式贴插设备的设计上,大多数都由电机+滚珠丝杆的方式进行驱动[5-69-10],直到近些年才出现以直线电机作为驱动方式的插件机。直线电机直接驱动的方式也被称为直驱技术,由电机直接和运动执行部分结合,没有中间的机械传动环节,结构大幅度简化,具有高速高精度、高效低耗、高刚度快响应和低噪音等优点,被广泛运用在工业领域[1112]。但是随着研究的深入,发现了直线电机中存在一些问题。其中最受关注的是直线电机的推力波动现象,它会引起动子速度和频率波动,产生振动和噪声,甚至与零部件发生共振现象,影响插件精度,甚至损坏零件。

    对直线电机在极限加速度下的推力波动和随机频率范围进行实验测试,实验平台如图2所示。

    图  2  测试平台

    将大理石梁上的直线电机和横梁上的直线电机分别安装在实验平台上进行测试。实验仪器主要有力传感器、示波器、变送器等。反馈元件采用海德汉的封闭式高精度直线光栅尺,分别采集不同运动距离下的直线电机的频率和推力数据。利用博世力士乐NYCe-4000系统的NYCeScope软件监测运动过程中的推力数据。

    直线电机参数设置:峰值电流27.6 A、持续电流19.7 A、力矩常数152 N/A、极距32 mm、电感12 mH、电阻2.5 Ω。大理石梁上的直线电机的频率测试界面如图3所示。

    图  3  频率测试界面

    将数据导出,用Origin2019绘制频率波峰瀑布图,两种直线电机的频幅图分别如图4图5所示。

    图  4  大理石梁直线电机频幅图

    图4图5可以看出,大理石梁和横梁上的直线电机的主要电流幅值对应的频率集中分布在30~50 Hz频率带内,其他频率对应的电流波峰幅值相对较低。

    大理石梁和横梁上直线电机的推力如图6图7所示。

    图  5  横梁直线电机频幅图
    图  6  大理石梁推力图
    图  7  横梁直线电机推力图

    图6图7可知,大理石梁上的直线电机产生的推力在780 N左右,横梁上的直线电机产生的推力在390 N左右。

    模态分析及其相关技术被广泛应用于精密机械和龙门机床的研发和优化设计中[1315]。对高速高加速度的异形插件机来说,结构的动态特性对插件精度有着十分重大的影响。模态分析可以计算出固有频率和振型,找到整机结构的薄弱部分,是评判异形插件机机械性能和优化结构的重要手段。

    将异形插件机视为由多个弹簧组成的一个自由度为n的线性系统,它的模态理论方程为

    $$ {\boldsymbol{M}}\ddot{X}+{\boldsymbol{C}}\dot{X}+{\boldsymbol{K}}X={\boldsymbol{F}}\left(t\right) $$ (1)

    式中:$ \ddot{X} $、$ \dot{X} $和X分别为加速度、速度和位移向量;M为质量矩阵;K为刚度矩阵;C为阻尼矩阵;F(t)为外部激励力矩阵。

    前面介绍了异形插件机的移动装置的传动方式为直线电机直接驱动,直线电机产生的摩擦较小,忽视整个设备的阻尼,即C=0,并假设结构为自由振动,则有

    $$ M\ddot{X}+KX=F\left(t\right) $$ (2)

    对于线性系统,MK取常数,自由振动为简谐运动。

    $$ u={\varphi }_{i}\mathrm{c}\mathrm{o}\mathrm{s}\left({\omega }_{i}t\right) $$ (3)

    典型的无阻尼模态分析求解的基本方程是经典的特征值问题。

    $$ K{\varphi }_{i}={\omega }_{i}^{2}M{\phi }_{i} $$ (4)

    式中:$ {\varphi }_{i} $为特征向量;$ {\omega }_{i}^{} $为i阶模态的固有频率(特征值)。

    理论上对该方程求解可以得到m个解,即得到结构自由固有频率为

    $$ {\omega }_{1}\leqslant \cdots \leqslant {\omega }_{i}=\sqrt{\frac{{K}_{i}}{{M}_{i}}}{\omega }_{1}\leqslant \cdots \leqslant {\omega }_{m} $$ (5)

    但是对于运动起主导作用的主要是前面的几阶模态,对于异形插件机,计算其前六阶固有频率即可。

    对整机结构进行模态分析,研究整机结构的前六阶固有频率和振型。模态分析步骤如下。

    (1)简化模型:将插件头替换成同等体积的箱体。去除PCB运输底板、拖链和线缆,删除地脚、圆角、倒角以及螺栓孔等。

    (2)编辑材料:大理石梁的材料为大理石,横梁为ZL104,直线电机与横梁间的连接板为6061(铝),其余部分使用Q235。

    (3)网格划分:利用Hypermesh对横梁部分的网格进行细化处理。

    (4)添加边界条件:在底座与地面接触处添加固定约束,螺栓连接处用固定约束代替,横梁与支撑梁上滑块的约束设置为位移约束(Y=0)。

    (5)求解:模态分析的结果如图8所示。

    图  8  整机前六阶模态振型图

    各阶模态振型描述见表1。由前六阶整机固有频率可知,一阶固有频率为52.28 Hz,而直线电机产生的随机频率范围为0~50 Hz,因此异形插件机与直线电机在运动过程中产生的随机频率可能会发生共振现象,影响插件精度。所以需要提升整机的低阶固有频率,避开直线电机产生的随机频率范围。

    表  1  优化前整机固有频率与振型描述 Hz
    阶数 固有频率 主要振型描述
    1 52.28 底座沿Z方向摆动,横梁扭转,插件头位置振型变化较大
    2 75.12 横梁跟随大理石梁绕对称轴在X-Z平面转动,横梁末端振型变化较大
    3 76.80 支撑梁沿Z方向摆动,支撑梁振型变化较大
    4 98.86 横梁沿X方向摆动,横梁末端振型变化较大
    5 105.47 整机绕中轴线在X-Z平面摆动,横梁末端振型变化较大
    6 109.60 横梁扭转,插件头位置变化较大
    下载: 导出CSV 
    | 显示表格

    通过零部件优化来提升整体的固有频率。由各阶振型可知:底座Z方向振型变化较大,可以通过适当添加加强筋使底座强度提升;支撑梁固定在底座上,可以通过优化内部结构,减少重量等来提升整体的固有频率;对横梁内部结构进行优化,提升横梁结构的固有频率。

    横梁使用拓扑优化的方法进行结构优化,详见文献[16]。底座的优化采用焊接钢板和加强筋的方式,优化后的底座结构如图9所示。优化后的整机模态分析结果如图10所示。优化前后的整机固有频率对比见表2

    图  9  优化后的底座结构
    图  10  优化后的整机模态

    优化后的整机一阶主要振型为底座沿Z方向摆动,二阶主要振型为大理石梁和横梁绕对称轴在X-Z平面转动,三阶主要振型主要为横梁沿X方向摆动,四阶主要振型为支撑梁沿Z方向摆动,五阶主要振型为横梁扭转,六阶振型为整机绕中轴线在X-Z平面转动。

    表  2  优化前后的整机固有频率对比表 Hz
    阶数 1 2 3 4 5 6
    优化前 52.28 75.12 76.80 98.86 105.47 109.60
    优化后 68.01 86.26 96.24 108.08 119.76 124.60
    下载: 导出CSV 
    | 显示表格

    优化前后的前六阶模态主要振型特点基本上相同,优化后的各阶固有频率分别提升30.1%、14.8%、25.3%、9.3%、13.5%和13.7%。可见优化效果较明显。一阶固有频率由52.28 Hz提升为68.01 Hz,避开了直线电机的频率范围。

    谐响应分析是一种线性分析,仅对稳态受迫振动进行分析,不考虑瞬态振动的影响[17]。谐响应分析是计算在已知频率的载荷作用下结构响应的技术,对结构设计、优化及避免产生共振等具有重要意义[18-19]。在模态分析的基础上进行谐响应分析,两台直线电机位置分别添加780 N和390 N的持续推力,求解频率范围为0~150 Hz,间隔为5。优化后的整机谐响应分析结果如图11所示。

    图  11  谐响应分析求解结果

    选取横梁和箱体上的点为采样点,求解振幅与频率的响应数据。4个取样点的位置如图12所示。4个取样点各方向振幅-频率响应曲线如图13所示。

    图  12  4个取样点的位置示意图
    图  13  四个取样点各方向上的振幅-频率响应曲线

    由采样点的位移-频率响应曲线可知:

    (1)采样点各方向上的振幅–频率响应曲线的整体趋势一致。在X方向,采样点的振幅变化随着频率的增大而缓慢增大,在0~50 Hz,振幅在0.03 mm左右;当频率为65 Hz时,振幅达到第一个峰值,该峰值在0.06 mm左右;当频率为85 Hz时,振幅达到第二个峰值,该峰值在0.18 mm左右;当频率为95 Hz时,振幅达到第三个峰值,该峰值在0.27 mm左右,为最大峰值;当频率为125 Hz时,振幅达到第四个峰值,该峰值在0.05 mm左右。在Y方向,采样点的振幅变化随着频率的增大而缓慢增大,采样点达到各峰值的频率与X方向基本一致,不同的是,采样点3和采样点4在65 Hz时的振幅最大,而采样点1和采样点2在120 Hz时的振幅最大。在Z方向,采样点的振幅变化随着频率的增大而缓慢增大,在0~50 Hz,振幅在0.02 mm以内、在65 Hz时为最大振幅,该值在0.045 mm左右;此后出现峰值的频率与X方向基本一致。因此,采样点在65 Hz、85 Hz、95 Hz和125 Hz容易发生共振,且在65 Hz时,Y方向和Z方向的振动位移较大,在95 Hz时,X方向的振幅较大。

    (2)为了提高插件机的精度,延长使用寿命,插件机的外部激励应避免上述共振频率区域。结合直线电机的产生的随机频率,可以确定优化后的异形插件机不会与之产生共振现象。

    以横梁为采样对象,计算其应力-频率响应曲线,结果如图14所示。

    图  14  横梁的应力-频率响应曲线

    横梁上,XYZ这3个方向上应力幅值对应的频率为65 Hz,85 Hz、95 Hz、110 Hz、125 Hz。其中Z方向上的应力幅值相对较大,在65 Hz时为0.23 MPa。XY方向的应力幅值在95 Hz时较大,分别为0.03 MPa和0.06 MPa。XYZ方向上应力响应在0.25 MPa以内,远小于材料许用应力,即使发生共振,也不会导致横梁损坏,满足工作要求。

    测试系统使用博世力士乐NYCe4000运动控制系统,是集成驱动传动技术的一款功能齐全的运动控制系统。

    实验过程为插件头在水平面内重复运动,不进行取料和插件的升降动作。PCB测试板固定在运输轨道上,PCB测试板上有很多专门用于测试的点。选择PCB测试板上的4个点,使用相机进行对焦,并设定为运动轨迹坐标。控制插件头在PCB测试板上方的水平面内按设定的坐标运动,运动路径为正方形,即插件头将在水平面按“1”→“2”→“3”→“4”→“3”→“2”→“1”的顺序循环运动。将“1”设置为Mark点,第一次运动到Mark点的坐标设为(0,0)。此后插件头每次运动到Mark点,相机进行一次拍照,经过图像处理换算后得到Mark点在运动过程中的坐标。图15所示为测试示意图。Mark点的坐标变化如图16所示。

    图  15  测试示意图
    图  16  Mark点坐标变化图

    图16可知,优化前的Mark点坐标在X方向和Z方向的变化范围基本上在±0.04 mm以内,有时会出现0.05 mm,说明了运动过程中结构不稳定,可能存在共振现象。优化后的异形插件机的Mark点坐标在XZ方向的变化范围稳定在±0.03 mm以内,该数值与谐响应分析结果在0~50 Hz时的振动幅值较为一致,验证了谐响应分析的可靠性,也说明了优化后的异形插件机在直线电机的驱动下不会产生共振现象。

    文章针对一种由直线电机驱动的异形插件机进行振动特性研究,并且以提升异形插件机的结构稳定性和提高重复精度为目的进行结构优化。通过实验和仿真结合的方法验证了仿真结果的可靠性和结构优化的有效性,得到以下结论:

    (1)所研究的直线电机产生的随机频率为0~50 Hz,在30~50 Hz频率带电流幅值较大。优化前的整机模型第一阶固有频率为52.28 Hz,与直线电机产生的随机频率容易发生共振现象。

    (2)通过模态分析确定了整机结构的薄弱位置,对薄弱位置的零部件进行优化设计可以有效地提高整机结构的稳定性。优化后的整机模型第一阶固有频率为68.01 Hz,较优化前提升30.1%。重复精度由±0.04 mm提升为±0.03 mm,验证了结构优化的有效性。

    (3)以直线电机的推力为边界条件,在模态分析基础上对优化后的整机结构进行谐响应分析,计算异形插件机的振动幅值。频率为0~50 Hz的激励下,采样点的振动位移在0.03 mm左右,结合重复精度的测试结果可以确定基于直线电机驱动的异形插件机不会发生共振现象。

    (4)通过仿真分析结合实验验证的方法来研究结构的振动特性,未来将通过其他方法来研究异形插件机的振动特性以及其他因素对振动特性的影响。

  • 图  1   异形插件机的整机结构

    图  2   测试平台

    图  3   频率测试界面

    图  4   大理石梁直线电机频幅图

    图  5   横梁直线电机频幅图

    图  6   大理石梁推力图

    图  7   横梁直线电机推力图

    图  8   整机前六阶模态振型图

    图  9   优化后的底座结构

    图  10   优化后的整机模态

    图  11   谐响应分析求解结果

    图  12   4个取样点的位置示意图

    图  13   四个取样点各方向上的振幅-频率响应曲线

    图  14   横梁的应力-频率响应曲线

    图  15   测试示意图

    图  16   Mark点坐标变化图

    表  1   优化前整机固有频率与振型描述 Hz

    阶数 固有频率 主要振型描述
    1 52.28 底座沿Z方向摆动,横梁扭转,插件头位置振型变化较大
    2 75.12 横梁跟随大理石梁绕对称轴在X-Z平面转动,横梁末端振型变化较大
    3 76.80 支撑梁沿Z方向摆动,支撑梁振型变化较大
    4 98.86 横梁沿X方向摆动,横梁末端振型变化较大
    5 105.47 整机绕中轴线在X-Z平面摆动,横梁末端振型变化较大
    6 109.60 横梁扭转,插件头位置变化较大
    下载: 导出CSV

    表  2   优化前后的整机固有频率对比表 Hz

    阶数 1 2 3 4 5 6
    优化前 52.28 75.12 76.80 98.86 105.47 109.60
    优化后 68.01 86.26 96.24 108.08 119.76 124.60
    下载: 导出CSV
  • [1] 王雪岐. 异形电子元件插件机视觉检测与定位技术研究[D]. 长沙:湖南大学,2019.
    [2]

    Hu G Q,Ma J W,Zuo Y M,et al. A flexible velocity planning algorithm for high-speed mounter with both efficiency and precision[J]. Journal of the Brazilian Society of Mechanical Sciences and Engineering,2021,43(12):1-18.

    [3] 武洪恩,倪良月,王凯,等. 基于Hu矩和递进Hough变换的SOT元件识别定位算法[J]. 组合机床与自动化加工技术,2019(8):4-7.
    [4]

    Huang Y H ,Li D B. Research on geometric positioning algorithm of vision system in high speed and high precision chip mounter[J]. Applied Mechanics and Materials,2014,519-520:667-671.

    [5] 孙继斌,骆敏舟,姜涛,等. 高精度贴片机横梁结构动态优化设计[J]. 机械设计与制造,2014(8):45-47,51.
    [6] 王莹. 高速贴片机关键部件动态模型及优化设计[D]. 青岛:山东科技大学,2020.
    [7] 刘舒颖,王福增,郭子宇,等. 基于有限元/离散元耦合的大理石高速划擦过程仿真[J]. 金刚石与磨料磨具工程,2019,39(1):95-100.
    [8] 徐宏伟,张广鹏,黄玉美,等. 大理石机床床身对主轴箱温度场分布的影响研究[J]. 中国机械工程,2011,22(11):1274-1278.
    [9] 姜凤鹏,莫锦秋. 基于ANSYS的贴片机床身结构优化设计[J]. 机电工程,2008(10):95-97.
    [10] 姜涛,唐火红,冯宝林. 拱架式贴片机横梁刚柔耦合动态仿真优化设计[J]. 机械设计与制造,2014(6):31-34.
    [11] 吕刚. 直线电机在轨道交通中的应用与关键技术综述[J]. 中国电机工程学报,2020,40(17):5665-5675.
    [12]

    Liang J T ,Ming Z F ,Li P D. System construction of a four-side primary permanent-magnet linear motor driv e mechanical press[J]. Frontiers of Mechanical Engineering,2020,15(4):600-609.

    [13] 漆力子,卢礼华,高强,等. 精密电动反射镜架的结构设计与动态特性分析[J]. 中国激光,2020,47(3):82-92.
    [14]

    Lu X,Chen Y W,Sheng X. Analysis of the influence of structural parameters on the modes of a gantry machine tool crossbeam based on analytical method[J]. Applied Sciences,2022,12(21):11083.

    [15]

    Chen P,Zhuang H,Chang Y,et al. Modal analysis of an aerostatic spindle system for ultra-precision machine tools at different spin speeds[J]. Advances in Mechanical Engineering,2018,10(9):168781401879852.

    [16] 胡意波,高自成,李立君,等. 插件机横梁结构优化与动态特性分析[J]. 机电工程,2023,40(1):144-151.
    [17] 郑彬,钟发,张敬东,等. 鼓式制动器模态及谐响应分析[J]. 中国农机化学报,2020,41(4):117-122.
    [18]

    Cecen F,Aktas B. Modal and harmonic response analysis of new CFRP laminate reinforced concrete railway sleepers[J]. Engineering Failure Analysis,2021,127:105471. DOI: 10.1016/j.engfailanal.2021.105471

    [19]

    Su W,Luo W H,Dong X S,et al. Harmonic response analysis for ball grid array package using computer finite element simulation[J]. Journal of Physics:Conference Series,2021,2033(1).

  • 期刊类型引用(0)

    其他类型引用(4)

图(16)  /  表(2)
计量
  • 文章访问数:  79
  • HTML全文浏览量:  11
  • PDF下载量:  32
  • 被引次数: 4
出版历程
  • 修回日期:  2023-04-29
  • 网络出版日期:  2023-12-03
  • 刊出日期:  2023-12-01

目录

/

返回文章
返回