对角修形对斜齿轮系统准静态及动态特性的影响研究
袁冰1, 常山1,2, 吴立言1, 刘更1, 陈策3     
1. 西北工业大学 陕西省机电传动与控制工程实验室, 陕西 西安 710072;
2. 中国船舶重工集团 第七○三研究所, 黑龙江 哈尔滨 150078;
3. 陆航研究所 机体发动机室, 北京 101121
摘要: 基于齿面承载接触分析方法,建立了斜齿轮副啮合刚度和轮齿修形耦合非线性激励计算模型,研究了对角修形对齿面载荷分布、综合啮合刚度以及承载传递误差的影响。基于Timoshenko梁理论建立了考虑轴段变形的斜齿轮-转子-轴承系统广义有限元模型,通过分离系统传递误差激振力,将系统参变微分方程组转化为定常微分方程组,实现了快速求解。以传递误差激振力波动量最小为优化目标,采用图解法得到了最佳对角修形参数。研究表明:对于修形齿轮副,啮合刚度波动量并不是影响齿轮传动系统振动的主导因素,传递误差激振力波动量的大小决定了系统振动的剧烈程度;对角修形不仅对轮齿刚度的削弱程度较小,而且可以显著降低动态传递误差均方根值和轴承动载荷波动量,进而降低系统振动噪声。
关键词: 斜齿轮     对角修形     啮合刚度     传递误差     轴承动载荷     傅里叶级数     均方差    

由于制造误差、装配误差和承载变形的存在, 齿轮传动系统将不可避免地产生啮合冲击和参数自激振动, 并通过轴、轴承等传递到箱体, 产生箱体结构噪声和辐射噪声, 降低设备的舒适性和安全性。为了降低齿轮传动系统振动噪声, 齿轮修形已成为重要的设计和工艺手段。

国内外学者和技术人员对直、斜齿轮副的齿面修形技术进行了大量研究。Conry和Seireg采用柔度矩阵法建立了齿面接触方程, 提出了直齿轮和斜齿轮的修形方法[1]。Houser以静态传递误差的啮合频率谐波成分幅值最小为目标优化了直齿轮齿廓修形参数[2]。Lin研究了直线和抛物线齿廓修形对直齿轮系统动态响应的影响, 发现抛物线修形对轮齿载荷的敏感性较低[3]。Maatar研究了齿廓修形和齿向修形对窄斜齿轮副准静态和动态特性的影响[4]。Bonori采用有限元法计算齿轮副传递误差, 并以传递误差波动量最小为目标函数, 采用遗传算法优化得到齿廓修形参数[5]。Chen建立了啮合刚度与轮齿误差耦合非线性激励解析计算模型, 研究了齿廓修形对直齿轮齿面承载接触特性的影响[6]。Ghosh建立了综合考虑时变啮合刚度、轮齿误差、齿面摩擦和齿侧间隙的直齿轮系统动力学模型, 以轴承动载荷波动量和啮合动载荷最大值加权最小构造目标函数, 研究了不同转速下的最佳齿廓修形参数[7]

除此之外, 部分学者开展了斜齿轮对角修形的相关研究。由于斜齿轮接触线与轴线不平行, 齿廓修形和齿向修形会减少斜齿轮实际接触线总长度, 显著降低重合度, 削弱轮齿强度, 而对角修形只对齿轮啮入啮出端沿接触线方向修形, 不仅可以降低齿轮振动噪声, 也可以保证修形前后实际接触线总长度变化较小, 重合度较高。王宪法以降低啮合刚度波动量为前提, 提出了对角修形的最佳修形量计算公式[8]。蒋进科建立了斜齿轮副纯扭转振动模型, 以承载传递误差幅值、啮合冲击力以及啮合线方向振动加速度均方根最小为优化目标, 采用遗传算法确定了对角修形参数[9]

从现有文献来看, 关于斜齿轮对角修形的研究相对较少, 且未详细分析对角修形对系统准静态特性 (即齿面承载接触特性) 的影响规律以及对系统振动激励和响应的影响机理。本文基于齿面承载接触分析方法, 建立了斜齿轮副啮合刚度和轮齿修形耦合非线性激励计算模型, 研究了对角修形对齿面承载接触特性的影响。建立了考虑轴段变形的斜齿轮-转子-轴承系统广义有限元模型, 通过分离系统传递误差激振力, 实现了傅里叶级数法快速求解, 提高了求解效率。以传递误差激振力波动量最小为优化目标, 采用图解法得到了最佳对角修形参数。

1 对角修形齿面设计

对角修形是指对齿轮啮入啮出端沿接触线方向的修形。由于主动轮由齿根啮入齿顶啮出, 从动轮由齿顶啮入齿根啮出, 对主动轮啮入位置修形等价于对从动轮啮出位置修形, 因此本文分别对主动轮齿顶 (啮出位置) 和从动轮齿顶 (啮入位置) 进行对角修形。

图 1 斜齿轮副啮合面及对角修形曲线

对角修形参数包括修形曲线、修形量和修形长度3个基本要素。修形曲线一般有直线、抛物线等高次曲线。本文采用二次抛物线修形曲线。对角修形参数可用啮合面数据表达, 如图 1所示, 对角修形区域为ΔABC和ΔDEF组成的阴影部分, BC和EF为修形起始位置且与接触线平行。为方便表示, 建立旋转坐标系x2oy2, 则对角修形曲面方程可表示为

(1)

式中, E1E2表示啮入啮出端修形量; xP表示接触点在旋转坐标系x2oy2中的横坐标, LAMLDN表示沿接触线法向的修形长度, E(x, y) 表示坐标为 (x, y) 的接触点的修形量。

2 齿面承载接触特性的计算

修形斜齿轮副的啮合过程可视为两弹性体的准静态承载接触过程, 依据文献[10]提出的圆柱齿轮齿面承载接触分析方法, 将连续的接触线离散为一系列接触点, 并将离散接触点处的承载变形分为宏观弯曲剪切变形和局部赫兹接触变形, 各个离散点处的修形量可视为接触点对之间的初始间隙量, 通过有限元子结构法提取齿面法向弯曲剪切变形柔度矩阵, 采用接触变形解析公式计算赫兹接触变形柔度矩阵, 则对角修形斜齿轮副齿面承载接触方程的矩阵形式可写为

(2)

式中, [λ]Global为接触点的宏观变形柔度矩阵; {F}为接触点的载荷向量; {u}Local为接触点的局部接触变形向量; {d}为接触点的剩余间隙向量; {ε}为接触点的初始间隙向量。求解该矩阵方程可得到修形斜齿轮副齿面载荷分布F和承载传递误差xs

则修形斜齿轮副综合啮合刚度为

(3)

修形斜齿轮副综合啮合误差为

(4)
3 系统动力学建模及求解 3.1 系统动力学建模

集中质量法不能考虑支承系统结构对系统动态响应的影响, 且不能详细计算各轴承支承位置的振动响应, 本节基于有限元法的基本思想, 将斜齿轮-转子-轴承系统划分为啮合单元、轴段单元和轴承单元3个子系统, 分别建立各单元的动力学模型, 通过有限元基本理论组装得到系统整体广义有限元模型。

图 2 斜齿轮副啮合单元

主动轮和从动轮被视为由弹簧-阻尼结构在基圆处相切连接的刚性圆盘。图 2为主动轮左旋且绕z轴正向逆时针旋转时斜齿轮副全自由度动力学模型, 其中, OpOg分别为主动轮和从动轮的回转中心, rprg分别为主动轮和从动轮的基圆半径, α为齿轮副啮合角, ψ为主动轮至从动轮中心连线矢量与x轴正向的夹角。定义啮合单元的节点广义坐标列向量为

(5)

啮合单元动力学方程的矩阵形式可写为

(6)

式中, Mm为啮合单元质量矩阵, Km为啮合单元刚度矩阵, Cm为啮合单元阻尼矩阵, em为齿面修形引入的齿轮副综合啮合误差。

考虑轴段结构的弯扭轴摆综合变形, 轴段单元采用Timoshenko梁单元模拟。轴承单元为轴承所连接的轴节点通过轴承刚度和轴承阻尼与箱体相连接。基于有限元法基本理论, 由各单元动力学模型组装得到系统整体动力学模型, 系统整体运动微分方程的矩阵形式可写为

(7)

式中, M为系统质量矩阵, C为系统阻尼矩阵, K为系统刚度矩阵, x为系统广义坐标列向量, e为系统综合啮合误差列向量, P0为系统外载荷激励列向量。

3.2 系统激振力的分离及求解

由 (7) 式可知, 系统整体运动微分方程为参变微分方程组, 必须采用数值积分法求解, 然而齿轮-转子-轴承系统自由度较多, 求解时间长, 在保证一定求解精度的前提下, 本节对系统运动微分方程进行近似变换处理, 将方程时变部分以系统激振力的形式转化到方程右侧, 从而将参变微分方程组转化为定常微分方程组, 此时不仅可采用频域求解方法快速求解, 大大提高求解效率, 而且可以更清晰地认识到系统振动产生的根本原因。

由系统整体运动微分方程式 (7) 退化可得到系统静力学平衡方程为

(8)

K(t)=K0K, (7) 式可转化为

(9)

xs(t)=xs0xs, 联立 (8) 式, 由于阻尼在非共振区对系统振动的影响较小, 此处忽略由阻尼引起的误差激振力[11], 则有

(10)

将方程右端K0Δxs定义为传递误差激振力波动量。显然, 传递误差激振力波动量由啮合刚度均值与传递误差波动量的乘积决定, 该时变项决定了系统振动的大小。此时, 原始参变微分方程组转化为定常微分方程组, 采用傅里叶级数法求解系统动力学方程可得到系统动态位移x(t)。

动态传递误差可表示为

(11)

轴承动载荷可表示为

(12)
4 算例分析 4.1 算例描述

以某斜齿轮-转子-轴承系统为研究对象, 其结构简图如图 3所示, 各轴段参数如表 1所示。采用十二自由度空间Timoshenko梁单元将系统离散化, 将主动轮所在轴段划分为13个单元, 从动轮所在轴段划分为11个单元, 共26个节点156个自由度, 节点由轴段节点、轴承节点和功率节点3种类型组成。啮合单元由主动轮和从动轮所在轴段节点通过弹簧和阻尼相连接构成。

图 3 系统结构简图
表 1 轴段参数
轴段参数/mm
轴段轴1轴2
外径内径长度外径内径长度
1900110200120260
211006023012030
31200100320120200
41300180200120300
5120010028012030
6110080
790070
表 2 斜齿轮副基本参数
部位主动轮从动轮
齿数z37106
模数mn/mm55
压力角α/(°)2020
螺旋角β/(°)1818
齿顶高系数ha11
顶隙系数c0.250.25
齿宽B/mm6060
4.2 对角修形对齿面承载接触特性的影响 4.2.1 修形量的影响

齿轮副总重合度为2.8297, 在一个啮合周期内, 齿轮副处于双齿与三齿交替啮合状态, 对啮入啮出位置进行对角修形, 采用文献[12]中的修形长度近似计算公式Ln=(εα-1)Pba得到对角修形长度为10.015 7 mm。

当修形长度不变时, 不同修形量对综合啮合刚度和承载传递误差的影响如图 4所示。

图 4 不同修形量对综合啮合刚度和承载传递误差的影响

可以看出, 当修形量小于12 μm时, 综合啮合刚度曲线变化较小, 三齿啮合区承载传递误差随着修形量的增加而增加。当修形量大于12 μm时, 三齿啮合区承载传递误差继续增加, 逐步大于双齿啮合区, 三齿啮合区部分啮合位置综合啮合刚度显著减小, 甚至小于双齿啮合区。这是因为对角修形会使齿轮啮入啮出位置的齿面载荷分布较未修形时发生很大改变。

修形量分别为0 μm、8 μm、12 μm、20 μm时的齿面载荷分布如图 5所示。

图 5 不同修形量时齿面载荷分布

可以看出, 未修形时, 齿轮副啮入啮出位置承担较大载荷, 当修形量较小时, 啮入啮出位置承担载荷明显降低。当修形量增加到一定程度时, 啮入啮出位置开始出现部分接触线不承担载荷。此时, 三齿啮合区部分啮合位置转变为双齿啮合, 然而此时实际接触线总长度减小较少, 有利于系统平稳。当修形量进一步增加时, 不参与接触的接触线长度逐渐增加, 实际接触线总长度减少过多, 将使原始三齿啮合区综合啮合刚度显著降低且小于双齿啮合区, 承载传递误差显著增大, 甚至大于双齿啮合区。

4.2.2 修形长度的影响

取修形量为12μm, 修形长度对承载传递误差和综合啮合刚度的影响如图 6所示。

图 6 修形长度对综合啮合刚度和承载传递误差的影响

与修形量的影响相类似, 当修形长度在一定范围内时, 静态啮合刚度的变化较小。当修形长度取值过大时, 将使三齿啮合区部分啮合位置实际接触线总长度减少过多, 导致该啮合位置综合啮合刚度减小, 承载传递误差增大。

4.2.3 修形量和修形长度的综合影响

从系统整体运动微分方程式 (10) 可知, 系统振动是由齿面承载接触特性参数承载传递误差和综合啮合刚度的时变性引起的, 为方便描述承载传递误差和综合啮合刚度的变化情况, 定义齿面承载接触特性参数变化率为

(13)

式中, ν表示齿面承载接触特性参数变化率; pn表示未修形齿面承载接触特性参数值; pm表示修形后齿面承载接触特性参数值。

当修形长度和修形量同时变化时, 啮合刚度波动量和承载传递误差波动量的变化如图 7图 8所示。可以看出, 当修形量和修形长度较小时, 啮合刚度波动量和承载传递误差波动量变化较小, 随着修形量和修形长度的增加, 啮合刚度波动量增幅由5%逐渐增大, 直至大于40%, 且不存在使啮合刚度波动量减小的区域。然而, 承载传递误差波动量随着修形量和修形长度的增加呈现先减小后增加的趋势, 存在一明显的窄带区域使承载传递误差波动量较小, 且降幅可达80%以上, 远远小于未修形, 此窄带区域即为最佳修形参数选择区域。此区域任意一组修形量和修形长度均可以使承载传递误差波动量较小。同时, 修形量越小, 则要求选取的修形长度越大, 反之亦然。

图 7 修形参数对啮合刚度波动量的影响 (单位:%)
图 8 修形参数对承载传递误差波动量的影响 (单位:%)

对角修形参数对啮合刚度均值的影响如图 9所示。可以看出, 随着修形量和修形长度的增加, 引入齿面的分布式误差量逐渐增大, 由于轮齿误差和啮合刚度存在非线性耦合, 啮合刚度均值总体上呈现非线性减小趋势。然而, 由于对角修形是在斜齿轮副啮入啮出位置沿接触线方向进行修形, 当修形量在一定范围内时, 修形前后接触线总长度变化较小, 因此啮合刚度均值的变化也较小, 且降幅由1%逐渐增加到4%。这也从侧面说明, 合适的对角修形对轮齿刚度和强度的削弱较小。

图 9 修形参数对啮合刚度均值的影响 (单位:%)
4.3 对角修形参数的优化及系统动态响应分析

从 (10) 式及以上分析可知, 斜齿轮-转子-轴承系统的振动由传递误差激振力决定, 本节以对角修形修形量和修形长度为设计变量, 以传递误差激振力波动量最小为目标函数对对角修形参数进行优化选取。由于传递误差激振力由啮合刚度均值和承载传递误差波动量的乘积构成, 从图 9可知, 对角修形对啮合刚度均值的影响很小, 因此传递误差激振力波动量最小可视为承载传递误差波动量最小, 采用图解法由图 8可得到最佳修形长度为10.5 mm, 最佳修形量为11 μm。

系统输入转速为1 500 r/min, 负载扭矩为25 00 N·m。由fm=n1z1/60可得到输入转速下齿轮副啮合频率为925 Hz。修形前后动态传递误差如图 10所示。由于对角修形引入齿面主动误差, 因此齿轮副动态传递误差增大, 然而其波动量却显著减小, 啮合频率一倍频和二倍频成分均显著下降。修形前后轴承1处沿y轴方向的轴承动载荷如图 11所示, 由于对角修形使传递误差激振力波动降低, 因此轴承动载荷波动量明显下降, 同时轴承动载荷的啮合频率一倍频和二倍频成分大幅降低。计算结果进一步说明, 对于对角修形齿轮副, 啮合刚度波动量并不是影响系统振动的主导因素, 承载传递误差波动量的大小直接决定了系统振动的剧烈程度。

图 10 对角修形前后动态传递误差
图 11 对角修形前后轴承1处沿y轴方向轴承动载荷

由于齿轮动态传递误差和轴承动载荷与相应节点的振动位移直接相关。因此, 不同转速下的动态传递误差和轴承动载荷与该转速下的系统模态关系密切。图 12图 13分别为不同转速下的齿轮副动态传递误差均方根历程和轴承1处沿y轴方向轴承动载荷波动量历程, 可以看出, 当输入转速为840 r/min和2 600 r/min时, 齿轮啮合频率为518 Hz和1 603 Hz, 分别与系统第9阶固有频率和第15阶固有频率相重合, 此时系统发生主共振, 当啮合频率等于第9阶固有频率和第15阶固有频率的1/i (i=1, 2, …, n) 倍时, 啮合频率的倍频成分将使系统发生超谐波共振。由于对角修形对啮合刚度均值的影响较小, 因此对角修形后系统并未发生类似齿向修形引起的共振峰偏移现象[13]。在宽转速范围内, 对角修形后齿轮副动态传递误差均方根值和轴承动载荷波动量与未修形时相比均显著减小, 可见对角修形的减振效果也适用于系统的不同转速工况。

图 12 动态传递误差均方根历程
图 13 轴承1处沿y轴方向轴承动载荷波动量历程
5 结论

1) 基于齿面承载接触分析方法建立了对角修形斜齿轮啮合刚度和齿面修形非线性耦合激励计算模型, 研究了对角修形参数对斜齿轮副综合啮合刚度、承载传递误差和齿面载荷分布的影响。当修形量在一定范围内时, 对角修形对轮齿刚度的削弱程度较小, 但可以显著减小承载传递误差波动量, 随着修形量和修形长度的增加, 啮合刚度波动量呈非线性增加趋势。

2) 基于有限元法基本思想, 建立了考虑轴段变形的斜齿轮-转子-轴承系统广义有限元模型, 通过分离系统传递误差激振力, 发现系统振动主要由啮合刚度均值和承载传递误差波动量的乘积决定, 将系统参变微分方程组转化为定常微分方程组, 实现了傅里叶级数法求解, 提高了求解效率。

3) 以传递误差激振力波动量最小为目标函数, 采用图解法得到了最佳对角修形参数。通过研究修形前后系统的时域和频域响应, 发现所设计的对角修形可有效减小动态传递误差均方根值和轴承动载荷波动量, 降低系统振动噪声, 同时反映了对于对角修形齿轮副, 啮合刚度波动量的大小并不是影响系统振动大小的主导因素, 承载传递误差波动量的大小直接决定了系统振动的剧烈程度。

参考文献
[1] Conry T F, Seireg A. A Mathematical Programming Technique for the Evaluation of Load Distribution and Optimal Modifications for Gear Systems[J]. Journal of Engineering for Industry, 1973, 95(4): 1115–1122. DOI:10.1115/1.3438259
[2] Tavakoli M S, Houser D R. Optimum Profile Modifications for the Minimization of Static Transmission Errors of Spur Gears[J]. Journal of Mechanisms, Transmissions, and Automation in Design, 1986, 108(1): 86–94. DOI:10.1115/1.3260791
[3] Lin H H, Oswald F B, Townsend D P. Dynamic Loading of Spur Gears with Linear or Parabolic Tooth Profile Modifications[J]. Mechanism and Machine Theory, 1994, 29(8): 1115–1129. DOI:10.1016/0094-114X(94)90003-5
[4] Maatar M, Velex P. Quasi-Static and Dynamic Analysis of Narrow-Faced Helical Gears with Profile and Lead Modifications[J]. Journal of Mechanical Design, 1997, 119(4): 474–480. DOI:10.1115/1.2826392
[5] Bonori G, Barbieri M, Pellicano F. Optimum Profile Modifications of Spur Gears By Means of Genetic Algorithms[J]. Journal of Sound and Vibration, 2008, 313(3): 603–616.
[6] Chen Z, Shao Y. Mesh Stiffness Calculation of a Spur Gear Pair with Tooth Profile Modification and Tooth Root Crack[J]. Mechanism and Machine Theory, 2013, 62: 63–74. DOI:10.1016/j.mechmachtheory.2012.10.012
[7] Ghosh S S, Chakraborty G. On Optimal Tooth Profile Modification for Reduction of Vibration and Noise in Spur Gear Pairs[J]. Mechanism and Machine Theory, 2016, 105: 145–163. DOI:10.1016/j.mechmachtheory.2016.06.008
[8] 王宪法, 李华敏. 窄斜齿轮对角修形的最佳修形量计算及修形区域的选择[J]. 哈尔滨工业大学学报, 1993, 25 (5): 62–66.
Wang Xianfa, Li Huamin. The Calculation of Optimal Diagonal Modification Values and Selection of Modification Area of the Narrow Helical Gears[J]. Journal of Harbin Institute of Technology, 1993, 25(5): 62–66. (in Chinese)
[9] 蒋进科, 方宗德, 王峰. 降低斜齿轮噪声的对角修形优化设计[J]. 振动与冲击, 2014, 33 (7): 63–67.
Jiang Jinke, Fang Zongde. Optimal Design with Diagonal Modification for Reducing Helical Gear Noise[J]. Journal of Vibration and Shock, 2014, 33(7): 63–67. (in Chinese)
[10] 常乐浩, 刘更, 吴立言. 齿轮综合啮合误差计算方法及对系统振动的影响[J]. 机械工程学报, 2015, 51 (1): 123–130.
Chang Lehao, Liu Geng, Wu Liyan. Determination of Composite Meshing Errors and Its Influence on the Vibration of Gear System[J]. Journal of Mechanical Engineering, 2015, 51(1): 123–130. (in Chinese)
[11] Özgüven H N, Houser D R. Dynamic Analysis of High Speed Gears By Using Loaded Static Transmission Error[J]. Journal of Sound and Vibration, 1988, 125(1): 71–83. DOI:10.1016/0022-460X(88)90416-6
[12] Velex P, Maatar M. A Mathematical Model for Analyzing the Influence of Shape Deviations and Mounting Errors on Gear Dynamic Behaviour[J]. Journal of Sound and Vibration, 1996, 191(5): 629–660. DOI:10.1006/jsvi.1996.0148
[13] 王奇斌, 张义民. 齿向修形直齿轮系统动力学特性分析[J]. 振动工程学报, 2016, 29 (1): 61–68.
Wang Qibin, Zhang Yimin. Dynamic Characteristics Analysis of a Spur Gear System with the Lead Crown Relief[J]. Journal of Vibration Engineering, 2016, 29(1): 61–68. (in Chinese)
Effect of Cross Modification on the Quasi-Static and Dynamic Characteristics of Helical Gear System
Yuan Bing1, Chang Shan1,2, Wu Liyan1, Liu Geng1, Chen Ce3     
1. Shaanxi Engineering Laboratory for Transmissions and Controls, Northwestern Polytechnical University, Xi'an 710072, China;
2. China Shipbuilding Industry Corporation 703 Institute, Harbin 150078, China;
3. Helicopter and Engine Department, Army Aviation Institution, Beijing 101121, China
Abstract: A nonlinear excitation model coupling mesh stiffness and tooth modification is established based on loaded tooth contact analysis method. The effects of cross modification on load distribution, mesh stiffness and loaded transmission error of helical gear are investigated. Considering shaft deflection, a generalized finite element model of helical-rotor-bearing system is developed based on Timoshenko beam theory. By separating transmission error exciting force, the differential equation with variable coefficients is transformed into ordinary differential equation, so the equation is solved fast. The best cross modification parameters are obtained based on graphical method with the objective of minimum of fluctuation of transmission error exciting force. The results show that mesh stiffness fluctuation is not the dominant factor affecting the vibration of gear transmission with modification, the system vibration is determined by the fluctuation of transmission error exciting force, and the cross modification not only weakens the gear mesh stiffness less, but also decreases the mean square error of dynamic transmission error and the fluctuation of dynamic bearing load significantly, thereby reduces gear transmission vibration and noise.
Key words: helical gear     cross modification     mesh stiffness     transmission error     dynamic bearing load     Fourier series     mean square error    
西北工业大学主办。
0

文章信息

袁冰, 常山, 吴立言, 刘更, 陈策
Yuan Bing, Chang Shan, Wu Liyan, Liu Geng, Chen Ce
对角修形对斜齿轮系统准静态及动态特性的影响研究
Effect of Cross Modification on the Quasi-Static and Dynamic Characteristics of Helical Gear System
西北工业大学学报, 2017, 35(2): 232-239.
Journal of Northwestern Polytechnical University, 2017, 35(2): 232-239.

文章历史

收稿日期: 2016-09-27

相关文章

工作空间