船舶力学

基于双自由度系统的转子特征值耦合特性分析

PDF下载

周兆华 ¹ , 雷晓薇 ¹ , 张雪冰 ¹ , 王刚伟 ¹ , 杨俊 ¹ , 何鹏 ²

船舶力学 | 流体与结构声学 2025,29(3): 506-516

收起

船舶力学 | 流体与结构声学 2025, 29(3): 506-516

基于双自由度系统的转子特征值耦合特性分析

全屏

周兆华¹, 雷晓薇¹, 张雪冰¹, 王刚伟¹, 杨俊¹, 何鹏²

作者信息

^1.武汉第二船舶设计研究所，武汉　430064

^2.哈尔滨工业大学，哈尔滨　150000

周兆华（1997-），男，硕士，工程师，通讯作者，E-mail：zhouxuexiooo@163.com

雷晓薇（1997-），女，硕士，工程师。

通讯作者:

通讯作者，E-mail：zhouxuexiooo@163.com

Coupling characteristics analysis of rotor eigenvalues based on two degrees of freedom system

Zhao-hua ZHOU¹, Xiao-wei LEI¹, Xue-bing ZHANG¹, Gang-wei WANG¹, Jun YANG¹, Peng HE²

Affiliations

^1.Wuhan Second Ship Design and Research Institute, Wuhan 430064, China

^2.Harbin Institute of Technology, Harbin 150000, China

出版时间: 2025-03-20 doi: 10.3969/j.issn.1007-7294.2025.03.015

文章导航

摘要

收起

随着大流量涡轮泵的迅速发展，其叶轮尺寸和转速日益增大，导致涡轮泵本身的耦合振动现象凸显。为保持涡轮泵的平稳运行并降低其振动噪声，本文研究了涡轮泵转子系统的耦合动力学特性。首先，基于简化的涡轮泵转子有限元模型进行分析，发现当转子不同阶频率在Campbell图中相互趋近时，接近重合的区域会发生转向和合并两种耦合特征现象。进一步分析发生耦合特征现象时的转子模态，发现互相耦合的不同阶特征模态会相互影响。接着，采用双自由度系统模拟转子耦合特征现象，分析耦合度、阻尼以及陀螺力矩对其特性的影响。最后，分析双自由度系统耦合动力学特性，发现系统发生转向及合并现象时，其动力响应会发生突变，导致系统耦合振动幅值增大，甚至系统失稳。研究成果可为转子系统安全稳定运行提供理论指导。

关键词

特征值耦合 / 动力学特性 / 转子系统 / 坎贝尔图

Abstract

收起

With the development of high flow turbopumps, the size and speed of impellers have increased rapidly, resulting in the coupling vibration phenomenon of the turbopump. In order to keep the turbopump running smoothly and reduce its vibration and noise, the coupling dynamics characteristics of the turbopump rotor system were studied in this paper. Firstly, based on the simplified finite element model of the turbopump, it was found that when the rotor frequencies of different orders approach each other in the Campbell chart, two coupling characteristics phenomena of curve turning and merging would occur. Further analysis of the rotor modes conducted for the coupling characteristic phenomenon showed the different orders of the coupling characteristic modes would affect each other. Secondly, a two-degree-of-freedom system was used to simulate the coupling characteristics of the rotor, and then the effects of coupling degree, damping and gyroscopic torque on the characteristics were analyzed. Finally, investigation on the coupling dynamic characteristics of the two-degree-of-freedom system revealed that the dynamic response would change abruptly when the curve turning and merging occurred, leading to the increase of the coupling vibration amplitude and even the instability of the system. The research results could provide theoretical guidance for the safe and stable operation of rotor systems.

Key words

eigenvalue coupling / dynamic characteristics / rotor system / Campbell chart

引用本文

周兆华, 雷晓薇, 张雪冰, 王刚伟, 杨俊, 何鹏. 基于双自由度系统的转子特征值耦合特性分析. 船舶力学, 2025 , 29 (3) : 506 -516 . DOI: 10.3969/j.issn.1007-7294.2025.03.015

Zhao-hua ZHOU, Xiao-wei LEI, Xue-bing ZHANG, Gang-wei WANG, Jun YANG, Peng HE. Coupling characteristics analysis of rotor eigenvalues based on two degrees of freedom system[J]. Journal of Ship Mechanics, 2025 , 29 (3) : 506 -516 . DOI: 10.3969/j.issn.1007-7294.2025.03.015

正文

收起

0　引言

收起

求解系统运动微分方程组是分析转子系统动力学特性问题最简单且有效的方法，由于直接求解多自由度系统方程组的解析解较为困难，因此目前学者们大多采用不同的数值方法进行计算。但无论是采用有限梁单元方法、全三维有限元方法或是其他常用的数值方法求解，在绘制系统的坎贝尔（Campbell）图时，以上方法均会由于计算量的限制（仅取少量的转速点来描述特征值随转速的变化趋势）而导致忽略了由耦合等因素引起的系统本质特征现象。

随着国内外学者对系统耦合方向研究的深入，Giannini^[1]基于动力学方程得到了考虑各影响因素的理论模型，并提出当系统考虑耦合等因素作用时，特征值不再是互不影响的线性变化曲线，而是发生交叉、转向或者合并的非线性变化；Pastora^[2]提出了一种衡量模态间耦合程度的置信度（MAC）比较方法，以归一化耦合作用；Chun等^[3]着重地考虑了科氏力、陀螺力矩以及离心力对系统动力学特性的影响；Anegawa等^[4-5]分析了叶片内部与轴平行运动的耦合，基于弹簧质量系统动力学方程，分析轴与叶片系统出现不稳定区域时的现象，与试验结果对应较好；周兆华^[6]采用CUF方法分析了叶片-轮盘-转子系统的耦合现象；牛宏伟等^[7]采用三维有限元软件对叶片模态进行分析，探究叶片耦合导致的频率转向和不同振型间的转换现象，得到了叶片动频曲线的转向趋势以及耦合时的振型规律。

综上所述，目前学者们针对转子系统Campbell图中的特征值耦合特征现象的研究较少，并且缺少对特征现象发生时的动力学特性的分析。本文首先对简化的涡轮泵转子的有限元模型进行特征值耦合现象分析，发现当特征值发生转向和合并时，转子模态将相互影响并发生变化。由于采用三维有限元的分析方法受到计算量过大以及影响因素复杂的制约，为了从机理上探究转子系统耦合特征现象的影响因素与动力学特性，本文基于双自由度系统动力学公式，模拟转子系统特征值接近时可能产生的转向以及合并现象，讨论在不同耦合度、阻尼以及陀螺力矩的取值下两自由度特征值的变化情况。最后，对双自由度系统施加简谐激励，分析系统在发生交叉、转向或合并时的时域响应，可为转子系统安全稳定运行提供理论指导。

1　转子耦合特性分析

收起

涡轮泵转子常用于较高转速的场合，在不考虑其他外部激励的情况下，涡轮泵的噪声主要是由强线谱引起的。在涡轮泵升速或降速的过程中，部分特殊频率的线谱特征非常明显，这部分特殊频率并非涡轮泵转子的倍频，而是涡轮泵转子自身的耦合频率。本文分析的原始涡轮泵转子模型如图1所示，其中图1（a）为涡轮盘的三维模型，图1（b）为将涡轮盘等效后的涡轮泵转子有限元模型（六面体单元数量16万个），转子两端由滑动轴承支撑，最高转速为25 000 r/min，涡轮泵在实际运行过程中会受到陀螺力矩、阻尼、离心力以及支撑刚度等因素的影响。

1.1　转子特征值耦合现象分析

由于原始涡轮泵转子有限元模型的计算单元数量较多，导致仿真计算量过大，为研究转子耦合现象，现将涡轮泵转子模型进行简化处理，简化为偏盘转子模型。简化后尺寸为：左轴长100 mm，直径20 mm；右轴长400 mm，直径20 mm；圆盘半径300 mm，厚度20 mm。转轴两端为两个轴承，约束轴向自由度，径向刚度为1×10⁸ N/m，同时考虑阻尼、陀螺力矩以及离心力的影响。模型采用八节点六面体单元，材料的杨氏模量为210 GPa，泊松比为0.3，密度为7850 kg/m³。

采用有限元分析软件对偏盘转子模型在不同转速下的模态进行分析。当数据点较少时，特征值曲线可能会将耦合区域略过，图形曲线几乎呈线性变化，无法体现出局部特征。因此，为观察特征值随转速的变化情况，需加密转速取值。综合考虑计算量和效果显示等因素后，选择计算转子前20阶模态，转速范围为0~27 000 r/min，每个转速的间隔为50 r/min，共计算540个数据点，得到如图2所示的Campbell图。由图可以看到，此时的各阶特征值是非线性变化的曲线，且特征值曲线会在特征值重合的区域发生突变。

对特征值曲线突变的区域进行局部放大分析，发现当两个或多个特征值曲线逼近时，特征值曲线会发生转向的现象，如图3（a）所示。特征值曲线转向后偏离原变化方向。除此以外，部分特征随曲线逼近时会发生“合并”的现象，如图3（b）所示。在“合并”的区域内，两特征值大小和走向完全一致，但不完全与原来的两个特征值方向相同，处于一种聚合的状态。上述两种现象，在特征值转向或合并后，随着转速增加特征值并不沿原轨迹方向前进，而是变化为与之相互影响的其他特征值方向。

1.2　转子耦合模态分析

单个轮盘的模态振型是不同阶的节圆、节径及其组合的模态，单个转轴的模态振型是不同阶的弯曲、扭转及其组合的模态。考虑阻尼、陀螺力矩以及离心力等因素的影响，对涡轮泵偏盘转子模型进行有限元模态分析，其模态振型主要跟轮盘和转轴的耦合作用大小有关，不同阶的模态振型分别是轮盘模态主导、转轴模态主导以及两者耦合作用下的复合模态。

由图2和图3（a）的Campbell图可以发现，转子模型第9阶特征值与第12阶特征值在转速升高的过程中会发生转向现象，转向后特征值沿原方向变化。由此分析第9阶模态和第12阶模态在转向前、转向时和转向后的模态变化情况，如图4所示。

图4（a）、图4（b）所示分别为第9阶和第12阶模态转向前后的模态云图，图4（c）、图4（d）分别为第9和12阶模态转向时的模态云图。由上图可以发现，不同阶模态之间形状并不相同，且在特征值转向前和转向后模态均无变化，而在转向时的模态不同于第9阶和第12阶模态中的任何一个。针对两阶模态分别分析会发现，当Campbell图中特征值发生转向时，其模态开始变化，在转向曲率最大时，模态的变化最为显著。此时的模态属于相互耦合模态的混合状态，然后随着转向曲率减小，其模态逐渐恢复原样。

图2、图3（b）的Campbell图可以发现，模态转子模型第5阶特征值与第7阶特征值在转速升高的过程中会发生“合并”现象，合并后的特征值沿原方向变化。由此分析第5阶模态和第7阶模态在合并前、合并时和合并后的模态变化情况，如图5所示。

图5（a）、图5（b）所示分别为第5和第7阶模态合并前后的模态云图，图5（c）、图5（d）分别为第5和7阶模态合并时的模态云图。不同阶模态之间形状并不相同，在合并前和合并后两者均无变化，在合并时的模态不同于两者中的任何一个。针对两阶模态分别分析时发现，当Campbell图中特征值刚发生合并现象时，其模态开始发生变化，在完全合并区域模态不再变化，此时的模态属于相互耦合模态的混合状态，然后随着合并状态结束，曲线再度分离，其模态变化逐渐恢复原样。

转子系统的耦合作用是复杂的，各自由度间均会有不同程度的影响，因此某处发生耦合现象的原因并非完全由相互靠近的某两个自由度引起，更可能是多个自由度在多种因素影响下共同作用的结果，因此多自由度系统在进行有限元分析时出现耦合现象的种类和频率是难以控制的。不仅如此，在进行多自由度有限元仿真时，在某一转速下施加激励是困难的，并且大多数情况下只能计算动态激励下的稳态结果，分析模型瞬态变化的计算量巨大。

2　双自由度系统耦合状态分析

收起

由上述分析可知，采用有限元数值方法分析多自由度转子系统时，存在无法区分两自由度间耦合作用的实际影响效果和难以分析耦合现象发生时的系统时域振动特性的问题。因此，本章针对双自由度系统模型进行分析，根据双自由度系统动力学公式，考虑实际涡轮泵转子系统的影响因素，如陀螺力矩、阻尼、离心力以及耦合作用，模拟转子系统特征值接近时可能产生的转向以及合并现象，并对耦合现象进行分析。

在采用惯性框架法推导转子动能时，考虑了离心力以及陀螺效应产生的动能，但由于离心力做功相互抵消，因此得到的运动方程仅增加陀螺效应项。有阻尼转子系统运动微分方程组为

(1)

式中，Ω为转速，M为质量矩阵，G为陀螺矩阵，C为阻尼矩阵，K为刚度矩阵，F为外部激励列向量。

为探究转子系统运动微分方程的耦合项、陀螺项以及阻尼项对系统特征值变化的影响规律，本文基于上述方程采用两自由度模型进行分析，如图6所示。其中各自由度质量分别为m、m₁，刚度为k、k₁，其间连接刚度为ε。

2.1　双自由度系统自由模态耦合分析

首先考虑两自由度无耦合、无阻尼以及无陀螺力矩作用的系统，即仅有刚度k和k₁的运动方程如下：

(2)

假设各自由度质量均为1 kg（统一单位制），则系统特征值为

(3)

由上述公式可知，两个系统特征值互不影响，当刚度k₁线性变化时，k的特征值无变化，k₁的特征值随k₁线性变化，发生如图7所示的交叉现象。

由图7可知，当两自由度间无耦合等影响因素时，两自由度特征值互不影响。

当考虑系统耦合作用时，两自由度间存在耦合项ε，则系统运动方程^[1]如下：

(4)

假设各自由度质量均为1 kg（统一单位制），则系统特征值为

(5)

当考虑系统耦合作用时，两自由度间存在耦合项ε，假设初始刚度k等于1 N/m，当ε的值等于0.01 k和0.02 k时，即耦合刚度等于0.01 N/m以及0.02 N/m，两自由度特征值随刚度k₁（0.9~1.1 N/m）的变化如图8所示，图中k₀和k₁₀为无耦合时候的特征值变化曲线。

由此可见，刚度k₁线性变化时，耦合项的存在会导致原本交叉的两条特征值曲线发生转向现象。存在耦合项时的特征值均大于原本该自由度下的特征值，耦合刚度越大，则转向区域的范围越大，特征值较耦合前增大得也越多。当耦合项为0时，两特征值曲线呈交叉状态（无转向）。转向是更加普遍的现象，只有当两个自由度完全不耦合时才会出现交叉现象。

当考虑系统耦合项与非对称滞后阻尼项时，滞后阻尼与位移项有关，将滞后阻尼公式代入方程得到的系统运动方程为

(6)

对于各自由度质量均为1 kg（统一单位制）的系统特征值为

(7)

假设初始刚度k等于1 N/m，耦合项等于0.01 k不变，增加除了耦合项外的非对称滞后阻尼项jη与jη₁，改变阻尼项的数值，两自由度特征值随刚度k₁（0.9~1.1 N/m）的变化如图9所示。

由图9（a）可知，当两个阻尼系数的差值小于2倍的耦合项ε时，无论阻尼系数为何值，均只发生转向现象，且差值越小转向范围越大。当差值为0时，转向范围只取决于耦合项。由图9（b）可知，此时无论阻尼取值多少，特征值大小均无变化，即当差值大于耦合项的2倍时，特征值只发生交叉现象。分析上述结果可知，阻尼项的存在可能抵消耦合项的作用。

当考虑系统陀螺项以及耦合项时，系统运动方程^[1]为

(8)

对于各自由度质量均为1 kg（统一单位制）的系统特征值为

(9)

当增加除了耦合项外的陀螺项α时，取耦合项为0.01 k不变，两自由度特征值随刚度k₁（0.9~1.1 N/m）的变化如图10所示。

随着陀螺项的值从小变大，观察图10（a）~（c）中两个自由度特征值的变化情况。陀螺项小于耦合项时，两特征值曲线仍是转向状态；当陀螺项等于耦合项时，两特征值曲线处于交叉状态；当陀螺项大于耦合项时，两自由度的特征值曲线会发生合并现象。

分析图10可知，在转向状态时，由于陀螺项的存在，耦合项的作用会减小，也使转向的范围减小。两特征值曲线处于交叉状态时，虽然此时陀螺项的存在抵消了耦合项的转向作用，但是特征值的大小相比无耦合状态时有了明显的增大，两个自由度增大的幅值均匀且相同。当曲线处于合并状态时，陀螺项数值越大合并区域越宽，两自由度特征值在合并区域完全一致，并呈现聚合状态。虽然由于耦合和陀螺力的存在，合并后的特征值整体仍大于原无耦合特征值，但是呈现聚合状态的合并区域会使得局域特征值小于原无耦合状态的特征值。

当耦合系统考虑非对称滞后阻尼作用时，由于滞后阻尼与位移项有关，将滞后阻尼公式代入方程得到的系统运动方程^[1]为

(10)

对于各自由度质量均为1 kg（统一单位制）的系统特征值为

(11)

当增加除了耦合项外的非对称滞后阻尼项jη与jη₁时，耦合项为0.01 k不变，改变阻尼系数的大小，发现当两个阻尼系数的差值小于耦合项的2倍时，无论阻尼系数为何值，均只发生转向现象，如图11（a）所示。当差值大于耦合项的二倍时，特征值发生交叉现象，如图11（b）所示。因此可以推断出，阻尼项的存在可能抵消耦合项的作用使其发生交叉现象。

当陀螺耦合系统考虑阻尼力时，系统运动方程为

(12)

对于各自由度质量均为1 kg（统一单位制）的系统特征值为

(13)

当系统包含耦合、阻尼和陀螺力时，各项不同的取值会带来不同的特征值耦合图，但是整体均为交叉、转向和合并的组合，此时特征值耦合曲线示意如图12所示。

2.2　双自由度系统强迫振动分析

为探究两自由度系统特征频率交叉、转向及合并前后的时域振动情况，分别对上述双自由度系统（m和m₁）中的m₁自由度施加简谐激振力f=Asin（ωt+φ），其中A为激振力幅值、ω为周期频率（与系统模态频率相近），并求解非齐次线性微分方程组，得到两自由度振动位移随时间的变化，分析其振动规律。

在交叉前、交叉时以及交叉后的两自由度时域振动情况如图13所示，图中位移为相对值，其大小与激励力幅值相关。由时域振型图可以发现，两自由度间无论何时都是互不干扰，只有受激励的自由度m₁发生拍振现象。

在转向前和转向后的双自由度系统时域振动情况如图14（a）~（b）、图14（e）~（f）所示，由图可以发现受到直接激励的自由度m₁的振幅较大，与m₁耦合的自由度m的振幅极小（与耦合等影响因素大小有关）。在转向时的双自由度时域振动情况如图14（c）~（d）所示，由图可以发现此时自由度m的振动幅值突然增大，其幅值略小于直接受到激励的自由度m₁，且两个自由度振动幅值大小发生交替变化。

在合并前和合并后的双自由度时域振动情况如图15（a）~（b）、图15（e）~（f）所示，由图可以发现受到直接激励的自由度m₁的振幅较大，与m₁耦合的自由度m₂的振幅较小（与耦合等影响因素大小有关）。在合并时的双自由度时域振动情况如图15（c）~（d）所示，自由度m和m₁的振动幅值均随着时间以指数形式增加，且不受直接激励的自由度m的振幅增长速率更快，双自由度系统处于不收敛的状态。

3　结论

收起

为探究涡轮泵转子耦合动力学特性，本文采用基于双自由度模型的解析方法，解决了有限元等常规方法难以分析耦合现象发生时的系统时域振动特性的问题和无法区分各自由度间实际耦合效果的问题。

本文考虑了实际涡轮泵受到的影响因素，模拟了耦合现象发生时的系统特征。采用双自由度系统模型可明确在无外力作用的情况下，某一自由度的变化完全是由另一个自由度的耦合作用所引起。分析发生耦合现象的双自由度系统后发现如下规律：（1）系统发生转向和合并现象的基础是考虑各自由度间存在的耦合效应；（2）阻尼能够降低各自度间耦合作用的效果；（3）离心作用仅改变系统的刚度大小，对耦合现象的本身无影响；（4）陀螺力矩对系统耦合效果具有一定的影响，并且陀螺力矩的存在是出现合并现象的必要条件。

进一步分析双自由度系统发生耦合现象时的时域振动特性，还可以发现如下规律：（1）非受激自由度的振动幅值在系统转向现象发生时会明显增大，但其幅值仍小于受激自由度的振动幅值，且两自由度的幅值不断交替变化；（2）两个自由度的振动幅值在合并时均会以指数速率增大，两者均不收敛，且非受激自由度振动幅值增大的速率要大于受激自由度振动幅值增大的速率。

涡轮泵转子系统在实际的运行过程中存在耦合、陀螺力矩以及阻尼等多种影响因素，系统各自由度间耦合作用复杂。因此，由上述分析结果可知：涡轮泵转子运行的转速应首先避开合并现象发生的区域，因为此时的涡轮泵转子系统处于失稳的状态；其次涡轮泵转子运行转速应尽量避开转向现象发生的区域，因为此时涡轮泵转子系统受到的激励会激起与耦合相关的多种振型，其幅值处于交替变化的状态，会导致涡轮泵转子系统的振动频率更为复杂。综上，为使涡轮泵转子系统能够安全稳定地运行，其运行转速应避开转向和合并范围所对应的转速区域。

基金

收起

工信部重大专项(2017-IV-0008-0045)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

Giannini

. Experimental characterization of veering crossing and lock-in in simple mechanical systems[J]. Mechanical Systems and Signal Processing, 2016, 72-73:846-864.

[2]

Pastora

, Binda

, et al. Modal assurance criterion[J]. Procedia Engineering, 2012, 48:543-548.

[3]

Chun

S B

, Lee

C W

. Vibration analysis of shaft-bladed disc system by using substructure synthesis and assumed modes method[J]. Journal of Sound and Vibration, 1996(189): 587-608.

[4]

Anegawa

, Fujiwara

, Matsushita

. Blade-shaft coupled resonance vibration by using active magnetic bearing excitation[C]//ASME Turbo Expo 2008: Power for Land, Sea, and Air, 2008.

[5]

Anegawa

, Fujiwara

, et al. Resonance and instability of blade-shaft coupled bending vibrations with in-plane blade vibration[J]. International Journal of Fluid Machinery & Systems, 2008, 1(1): 169-180.

[6]

周兆华. 基于扩展梁的叶片-轮盘-转子系统动力学耦合特性研究[D]. 哈尔滨：哈尔滨工业大学，2021.

Zhou

Z H

. Study on dynamic coupling characteristics of blade-route-rotor system based on extension beam[D]. Harbin: Harbin Institute of Technology, 2021. (in Chinese)

[7]

牛宏伟，郭海东，文敏. 航空发动机压气机转子叶片频率转向特性研究[J]. 机械研究与应用，2019, 32(3):37-39.

Niu

Hongwei

, Guo

Haidong

, Wen

Min

. Research on frequency steering characteristics of aero-engine compressor rotor blades[J]. Mechanical Research and Application, 2019, 32(3): 37-39. (in Chinese)

2025年第29卷第3期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.3969/j.issn.1007-7294.2025.03.015

接收时间：2024-09-18
首发时间：2026-03-24
出版时间：2025-03-20

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2024-09-18

基金

工信部重大专项(2017-IV-0008-0045)

作者信息

^1.武汉第二船舶设计研究所，武汉　430064

^2.哈尔滨工业大学，哈尔滨　150000

通讯作者:

通讯作者，E-mail：zhouxuexiooo@163.com

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/cblx/CN/10.3969/j.issn.1007-7294.2025.03.015

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT