船舶力学

量纲参数（单位）	值	无量纲参数（单位）	值
杨氏模量E/pa	2.1×10¹¹	附加质量系数C_a	1.0
管道长度L/m	50	固体颗粒阻力系数C_Ds	0.9
管道外径D/m	0.24	线性阻尼系数A₁	27000
管道内径d/m	0.21	非线性阻尼系数A₂	0.24
固体颗粒直径d_s/m	0.1d
管道密度ρ_p/（kg·m⁻³）	7800
内部流体密度ρ_f/（kg·m⁻³）	1020
内部固体密度ρ_s/（kg·m⁻³）	3500
末端集中块质量m_e/kg	1000

量纲参数（单位）	值	无量纲参数（单位）	值
杨氏模量E/pa	2.1×10¹¹	附加质量系数C_a	1.0
管道长度L/m	50	固体颗粒阻力系数C_Ds	0.9
管道外径D/m	0.24	线性阻尼系数A₁	27000
管道内径d/m	0.21	非线性阻尼系数A₂	0.24
固体颗粒直径d_s/m	0.1d
管道密度ρ_p/（kg·m⁻³）	7800
内部流体密度ρ_f/（kg·m⁻³）	1020
内部固体密度ρ_s/（kg·m⁻³）	3500
末端集中块质量m_e/kg	1000

深海粗颗粒矿石垂直提升管道振动稳定性特征研究

PDF下载

高云 ¹^,² , 陈晓东 ¹ , 张磊 ¹ , 齐正杨 ¹

船舶力学 | 流体力学 2026,30(8): 1279-1293

收起

船舶力学 |流体力学 2026 , 30 (8) : 1279 -1293

深海粗颗粒矿石垂直提升管道振动稳定性特征研究

全屏

高云¹^,², 陈晓东¹, 张磊¹, 齐正杨¹

作者信息

¹哈尔滨工业大学（威海）海洋工程学院，山东威海 264209

²上海交通大学海洋工程全国重点实验室，上海 200240

作者简介:

高　云（1985–），男，副教授，博士生导师，通讯作者，E-mail：gaoyun@hit.edu.cn

Vibration stability characteristics of vertical lifting pipelines for deep-sea coarse ore

Yun GAO¹^,², Xiao-dong CHEN¹, Lei ZHANG¹, Zheng-yang QI¹

Affiliations

¹School of Ocean Engineering, Harbin Institute of Technology, Weihai 264209, China

²State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

出版时间: 2026-08-15 doi: 10.3969/j.issn.1007-7294.2026.08.006

文章导航

摘要

收起

本文基于修正的哈密顿原理，建立了深海粗颗粒矿石垂直提升管道系统的非线性动力学控制方程；随后通过Galerkin方法对方程进行离散，将偏微分方程转化为一系列常微分方程构成的方程组；最后利用特征值法以及Newton-Raphson迭代法对管道的稳定性特征进行了数值求解。系统分析了颗粒进料浓度、末端集中质量块、管道质量以及液相速度这四个关键无量纲参数对结构临界流速的影响。研究结果表明，当进料浓度C_v=0.1时，结构仅在前三阶模态发生Hopf分岔，失稳集中于这三阶模态；当进料浓度增加到C_v=0.3时，第四阶模态也发生Hopf分岔，且最小临界流速u_cr显著降低。当内流速度小于临界速度时，系统复特征频率虚部均为正值，结构振动形式表现为衰减振动；当内流速度大于临界速度时，系统的一阶复特征频率虚部会出现负值，结构振动呈现颤振失稳。当管道质量α保持不变时，系统临界流速u_cr随无量纲末端质量κ₅的增大而逐渐上升；当 κ₅保持不变时，随着管道质量α的增大，临界流速显著升高，呈现近似抛物线的增长趋势。

关键词

固液两相流 / 流固耦合 / 稳定性特征 / 悬臂管道

Abstract

收起

First, based on the modified Hamiltonian principle, the nonlinear dynamic control equation of the deep-sea coarse-grained ore vertical lifting pipeline system was established. Then, the equation was discretized using the Galerkin method, converting the partial differential equation into a system of ordinary differential equations. Finally, numerical solutions were obtained for the stability characteristics of the pipeline using the eigenvalue method and Newton-Raphson iteration method. The system analyzed the influence of four key dimensionless parameters, namely the particle feed concentration, terminal lumped mass, gravity coefficient of the pipeline, and fluid velocity, on the critical flow velocity of the structure. Results show that when the feed concentration is C_v is 0.1, Hopf bifurcations occur only in the first three modes. When the feed concentration increases to C_v = 0.3, a Hopf bifurcation also occurs in the fourth mode, and the minimum critical flow velocity u_cr decreases significantly. When the internal flow velocity is below the critical value, the imaginary parts of the complex eigenfrequencies remain positive, indicating decaying structural vibrations. Once the internal flow velocity exceeds the critical threshold, the imaginary part of the first complex eigenfrequency becomes negative, signifying flutter instability. When the gravity coefficient α is held constantly, the critical flow velocity u_cr increases progressively with the dimensionless terminal mass κ₅. Conversely, with fixed κ₅, an increase in the gravity coefficient leads to a significant rise in critical flow velocity, exhibiting an approximately parabolic growth trend.

Key words

solid-liquid two-phase flow / fluid-structure interaction / stability characteristics / cantilevered pipeline

引用本文

高云, 陈晓东, 张磊, 齐正杨. 深海粗颗粒矿石垂直提升管道振动稳定性特征研究. 船舶力学, 2026 , 30 (8) : 1279 -1293 . DOI: 10.3969/j.issn.1007-7294.2026.08.006

Yun GAO, Xiao-dong CHEN, Lei ZHANG, Zheng-yang QI. Vibration stability characteristics of vertical lifting pipelines for deep-sea coarse ore[J]. Journal of Ship Mechanics, 2026 , 30 (8) : 1279 -1293 . DOI: 10.3969/j.issn.1007-7294.2026.08.006

正文

收起

0　引言

收起

开发海洋资源和利用海洋空间是21世纪人类开发海洋的两大核心主题，其中海底蕴含着极其丰富的矿产资源，其有效开发不仅具有重要的战略意义，也蕴含巨大的商业潜力^[1]。这些矿产资源大多富含于数千米水深的深海底部，导致其开采难度极大。水力提升式采矿系统是目前国际上公认的最具有实用价值和发展前景的开采系统，如图1所示，其主要由集矿机、柔性输送软管、中间仓、垂直提升管道以及采矿船等部分组成。在深海矿产资源开采过程中，具有较大长细比的悬臂式垂直输流管道是整个深海采矿系统中最为薄弱、技术难度最大的部分，同时也是最核心的部分。随着海洋工程设计与建造理论技术的飞速发展，非线性动力学问题日益凸显，亟待深入研究与解决。其中包括由外部海流诱导的柔性细长立管结构涡激振动以及内流激励下引起的振动稳定性与振动响应等。针对外流引起的涡激振动，已有广泛研究^[2−5]。由内流导致的振动现象仍面临诸多挑战，其成因主要包括内部流体运动产生的离心力与科里奥利力，这些力在高流速条件下可能引发管道失稳^[6]。因此，如何在高流速环境下确保海底管道系统在输送矿产资源过程中的结构稳定性与运行效率，已成为当前海洋工程领域的重要研究方向。

海洋工程中的输流管道，依据边界条件，可大致分为两类：一类是传统用于油气输运的两端固支的输流管道；另一类是用于矿石输运的一端固支、另一端自由的悬臂式输流管道。由于悬臂输流管道自由端与外界存在质量和动量交换，因此，与两端固支的输流管道相比，悬臂输流管道的动力学响应特征要更为复杂。针对悬臂式输流管道的研究，根据管内介质的相组成不同，又进一步细分为单相内流激励下输流管道动力响应研究以及多相内流激励下悬臂输流管道动力响应研究。

对于单相内流的研究，最早可追溯到上个世纪。1966年，Gregory和 Paidoussis^[7]首先基于牛顿法推导了管道的控制方程，研究了悬臂式管道在输送流体过程中的振动特性，并揭示了流体质量比与管道临界不稳定状态之间的定量关系。1974年，Paidoussis 和 Issid^[8]进一步从理论层面探讨了不同边界条件下输流管道的稳定性。研究发现，两端支撑的管道在较高流速下易发生由发散不稳定和模态耦合引起的颤振现象；而对于悬臂式管道，结构稳定性主要受限于第一阶单模颤振的影响。1980年，Païdoussis等^[9−11]对特殊约束与末端带有质量块的悬臂管道的频率与混沌现象进行了研究，发现系统在低流速下的不稳定现象与倍周期分叉导致混沌现象。步入新世纪后，非线性单相内流的研究模型逐渐得到了更多的重视。2005年，Kuiper和Metrikine^[12]推导了吸液提升管道的振动模型，以分析自由悬挂管道在水下环境中的稳定性。研究表明，负压对结构稳定性的影响较为有限。2007年，Kuiper和Metrikine^[13]基于吸液提升管道的振动模型，进一步分析了引入流体阻力与常规边界条件相结合的影响。通过对浸没式刚性圆柱体的实验数据分析，发现流体阻力会提高管道颤振临界速度。2019年，王乙坤等^[14]开展了悬臂管道在内流激励条件下的参数共振现象研究，指出分布式运动约束力在显著影响管道位移响应方面起到了关键作用。2020年，方孟孟等^[15]探讨了悬臂输流管道在基础激励与脉动内流共同作用下的动力学特性，结果表明相位差与频率比对系统的混沌行为比例具有显著调控作用。同年，柳博瀚等^[16]运用CFD单向、双向流固耦合分析方法，开展了管道内流对海洋弹性管道振动响应影响的研究，研究表明：相对比较小内流速度的对管道的振动频率及振幅影响较小；相对较大的内流速度，对横流向振幅及顺流向振型影响较为明显。2022年，赵志贤等人^[17]分析了该类管道的自由振动特性，研究发现，采用一端具弹性支承的约束方式有助于增强系统在自由振动中的稳定性。随后，漆发辉等^[18]进一步研究了分布载荷作用下悬臂输流管道的稳定性问题，揭示了分布载荷、质量比及黏弹性系数对系统临界失稳流速的影响规律。

相比于单相内流激励下悬臂输流管道的研究，针对多相内流激励下悬臂输流管道动力响应的研究起步较晚，在新世纪逐渐得到重视。2004年，Monette and Pettigrew^[19]首先提出了一种修正的气−液（G-L）两相流的滑移速度模型，该模型为后续的研究提供了关键的滑移速度公式。2010年，Zhang 和Xu^[20]通过实验手段分析了内部气流对管道振动特性影响，结果表明，系统的动态响应主要受气体体积分数控制，其变化显著影响管道的稳定性。2015年，An和Su^[21]基于Monette 和 Pettigrew提出的两相流振动模型，进一步优化滑移比因子模型，更准确地表征气−液（G-L）两相流的动态特性。研究表明，随着气体体积分数和流速的增加，管道的振动幅值显著上升。2019年，Ebrahimi-Mamaghani等^[22]在建模中引入耗散力的影响，提出了一种改进的数学模型，其分析结果同样强调了气体体积分数在影响管道稳定性方面的重要作用。2021年，王予琪等^[23]使用标准k-ω模型，运用颗粒相离散元法（CFD-DEM）耦合方法，构建相同体积浓度、不同颗粒粒径的流动模型来研究粒径对于泵外特性的影响，研究表明，混合输流道内，颗粒在入口处发生低速堆积现象，颗粒速度随径向距离增大而增大。2023年，Gao等^[24−25]分别构建了固液两相流的简支管道和悬臂模型，分析了对屈曲临界速度和结构非线性响应特性的影响。研究结果表明，无量纲长径比、无量纲液相体积和无量纲管道质量参数的变化将导致临界速度以线性或非线性形式变化。2024年，Su等^[26]以尺度垂直升力管道模型为研究对象，采用液相计算流体动力学CFD-DEM耦合方法，模拟垂直管道中的固液两相流。研究表明：颗粒的提升受输送参数和结构因素的影响。同年，Zhu等^[27]构建了柔性立管在内部液固两相流和外部剪切流激励下的耦合响应的实验模型，实验结果表明：当立管浸入静水中时，只有内部液固流激发了平面内响应。随着内部流动密度或速度的增长，平面内响应增强，导致压差和顶部张力的波动增加，当存在外部流动的情况下，平面外响应也会被激发。

由以上研究可以发现，针对内流激励下悬臂输流管道动力响应特征研究存在以下几个明显不足：第一、目前绝大多数研究，均是针对单相内流展开，而针对多相内流研究仍以气液两相流为主，固液两相流的研究多基于CFD方法对于管道内部固液两相流模型中的颗粒分布、内流速度场、以及管壁剪切应力等研究，未涉及到固液两相流作用于管道结构控制方程，也未分析其对管道稳定性的影响。此外，在有限的针对多相内流的研究中，大多均是针对气液两相内流展开，而针对固液两相内流的研究则十分少见。第二、在目前的研究中，绝大多数均假设管道质量分布具备均质特征，没有考虑附着在管道上的集中质量块引起的非均质特征。第三、目前绝大多数研究均是针对排流管道展开，而针对吸流管道的研究十分少见。因此，可以看出：针对用于深海采矿作业的垂直提升管道（该管道具备粗颗粒两相内流、下端弱约束、内流为吸流以及端部布置集中质量块等很多独有特征）在内流激励下的动力响应特征研究非常罕见，还有待深入研究。该研究对垂直输流管道早期设计及服役期的安全工作具有重要的理论和工程价值。

基于以上研究不足，本文从真实垂直提升管道内流流动特征、管道质量分布特征、以及边界约束特征出发，建立了粗颗粒两相内流激励下末端带集中质量块的悬臂式吸流管道动力响应理论分析和数值预报模型。本文结构为：在第1节和第2节中，分别建立内流激励下管道振动控制数学方程以及求解方程的数值分析方法；在第3节中，依次对不同关键无量纲参数下管道稳定性特征展开深入分析。

1　问题描述

收起

如图2所示，假设一悬臂垂直管道长度为L、外径为D、内径为d，考虑其在内流激励下引起的平面内振动响应问题。当悬臂管道为采矿管道时，其内部流体为固液两相流，且流体流动方向自下往上。取管道顶端为原点o；z方向与管道未变形时轴线方向重合，且定义向下为正方向；x方向与管道振动方向重合，且定义向右为正方向；x和z两个坐标方向构成了结构模型振动所在的二维平面直角坐标系。管道顶端采用固接边界条件，末端采用自由边界条件，且末端沿管道轴线方向附着一质量为m_e的集中质量块。

由于悬臂管道自由端与外界存在质量与动量交换，经典的Hamilton原理不适用于本文所研究的悬臂管道模型，因此需要通过修正的Hamilton原理建立管道振动方程。当自由端与外界存在质量与动量交换时，修正的哈密顿原理可表示为^[28−29]

$ \delta \int \nolimits_{{t}_{1}}^{{t}_{2}}\mathcal{L}\mathrm{d}t+\int \nolimits_{{t}_{1}}^{{t}_{2}}\delta W\mathrm{d}t=\int \nolimits_{{t}_{1}}^{{t}_{2}}MU\left(\frac{\partial {\boldsymbol{r}}_\text{L}}{\partial t}+U{\boldsymbol{\tau }}_\text{L}\right)\cdot \delta {\boldsymbol{r}}_\text{L}\mathrm{d}t $

(1)

式中：

$ \mathcal{L} $

为系统拉格朗日函数，W为非保守力做功。式右侧，M为管道内流总质量，U为管道内流速度，r_L为管道自由端位置向量，τ_L为管道自由端空间切向量。当仅考虑管道在xoz平面内振动时，r_L和τ_L可表示为

$ {\boldsymbol{r}}_\text{L}=x\boldsymbol{i}+z\boldsymbol{k},\quad {\boldsymbol{\tau }}_\text{L}=\frac{\partial x}{\partial s}\boldsymbol{i}+\frac{\partial z}{\partial s}\boldsymbol{k} $

(2)

式中：s是沿管道轴线的曲线坐标。值得注意的是，式（1）表示的是常见的内流为单相内流的修正Hamilton原理表达式。当内部流体为固液两相流、且假设固相和液相流动的速度大小分别定义为U_s和U_f时，式（1）可以改写为

$ \delta \int \nolimits_{{t}_{1}}^{{t}_{2}}\mathcal{L}\mathrm{d}t+\int \nolimits_{{t}_{1}}^{{t}_{2}}\delta W\mathrm{d}t=-\int \nolimits_{{t}_{1}}^{{t}_{2}}{M}_\text{f}{U}_\text{f}\left(\frac{\partial {\boldsymbol{r}}_\text{L}}{\partial t}-{U}_\text{f}{\boldsymbol{\tau }}_\text{L}\right)\cdot \delta {\boldsymbol{r}}_\text{L}\mathrm{d}t-\int \nolimits_{{t}_{1}}^{{t}_{2}}{M}_\text{s}{U}_\text{s}\left(\frac{\partial {\boldsymbol{r}}_\text{L}}{\partial t}-{U}_\text{s}{\boldsymbol{\tau }}_\text{L}\right)\cdot \delta {\boldsymbol{r}}_\text{L}\mathrm{d}t $

(3)

值得注意得是，公式（3）中固相速度和液相速度分别取为−U_s以及−U_f，这是由于本文研究的管道为吸流管道，管道内固相和液相流动方向均向上，与图1中定义的坐标轴z正向相反。式（3）中，M_f为管道内部液相质量，M_s为管道内部固相质量。

$ \mathcal{L} $

为系统拉格朗日函数，可表示为系统动能T与系统势能V之差，为

$ \mathcal{L}=T-V,\quad T={T}_{{\mathrm{p}}}+{T}_\text{s}+{T}_\text{f},\quad V={V}_{{\mathrm{p}}}+{V}_\text{s}+{V}_\text{f} $

(4)

式中：系统总动能T包括三部分，分别为管道结构动能T_p、管道内部固相动能T_s以及管道内部液相动能T_f。系统总势能V同样包括三部分，分别为管道结构势能V_p、管道内部固相势能V_s以及管道内部液相势能V_f。

管道内部固相绝对流动速度大小v_s以及液相绝对流动速度大小v_f，分别与管道振动速度以及各相相对于管道的流动速度有关，假设结构发生小变形振动，可设s ≈ z，当仅考虑管道横向振动时，v_s和v_f可表示为

$ {v}_\text{s}=\frac{\partial x}{\partial t}-{U}_\text{s}\frac{\partial x}{\partial s}\approx \frac{\partial x}{\partial t}-{U}_\text{s}\frac{\partial x}{\partial z},\quad {v}_\text{f}=\frac{\partial x}{\partial t}-{U}_\text{f}\frac{\partial x}{\partial s}\approx \frac{\partial x}{\partial t}-{U}_\text{f}\frac{\partial x}{\partial z} $

(5)

管道内部固相动能T_s以及管道内部液相动能T_f可表示为沿管道轴线方向的积分形式，并联立式（5），得到

$ {T}_\text{s}=\frac{1}{2}{m}_\text{s}\int \nolimits_{0}^{L}v_\text{s}^{2}\mathrm{d}z=\frac{1}{2}{m}_\text{s}\int \nolimits_{0}^{L}{\left(\dot{x}-{U}_\text{s}{x}^{\prime}\right)}^{2}\mathrm{d}z,\quad {T}_\text{f}=\frac{1}{2}{m}_\text{f}\int \nolimits_{0}^{L}v_\text{f}^{2}\mathrm{d}z=\frac{1}{2}{m}_\text{f}\int \nolimits_{0}^{L}{\left(\dot{x}-{U}_\text{f}{x}^{\prime}\right)}^{2}\mathrm{d}z $

(6)

式中：m_s以及m_f分别为单位长度管道内部的固相质量以及液相质量。值得注意的是，为了表述方便，在式（6）中，用“⋅”表示对时间t求偏导；用“′”表示对空间坐标z求偏导，后续为了表达方便均采用这种简化表达形式。如图2所示，当管道结构末端附着有集中质量块时，管道结构动能可表示为

$ {T}_\text{p}=\frac{1}{2}\left[{m}_\text{p}+{m}_\text{e}\delta \left(L-z\right)\right]\int \nolimits_{0}^{L}{\dot{x}}^{2}\mathrm{d}s $

(7)

式中：m_p为单位长度管道结构的质量，m_e为末端集中质量块的质量，δ（L − z）为Dirac函数，单位为L⁻¹，其具体表达式为

$ \delta \left(L-z\right)=\left\{\begin{aligned}& 0,& z\neq L\\& 1,& z=L\end{aligned}\right. $

(8)

式（4）中，管道内部固相重力势能V_s以及管道内部液相重力势能V_f可表示为^[30]

$ {V}_\text{s}=\frac{1}{2}{m}_\text{s}g\int \nolimits_{0}^{L}\left(L-z\right){{{x}^{\prime}}}^{2}\text{d}z,\quad {V}_\text{f}=\frac{1}{2}{m}_\text{f}g\int \nolimits_{0}^{L}\left(L-z\right){{{x}^{\prime}}}^{2}\text{d}z $

(9)

式（4）中，管道结构势能V_p包括两部分，分别为管道和集中质量块引起的重力势能V_p1以及由管道发生弯曲变形引起的弹性势能V_p2，表示为^[30]

$ {V}_{{\mathrm{p}}1}=\frac{1}{2}\left({m}_\text{p}g-{\rho }_\text{f}{A}_\text{e}g\right)\int \nolimits_{0}^{L}\left(L-z\right){{{x}^{\prime}}}^{2}\text{d}s+\frac{1}{2}{m}_\text{e}g\int \nolimits_{0}^{L}{{{x}^{\prime}}}^{2}\text{d}z $

(10)

$ {V}_{{\mathrm{p}}2}=\frac{1}{2}EI\int \nolimits_{0}^{L}{\kappa }^{2}\text{d}z+\frac{1}{2}\left({P}_\text{f}\left(z\right){A}_\text{e}-{P}_{{\mathrm{i}}}\left(z\right){A}_{{\mathrm{i}}}\right)\int \nolimits_{0}^{L}{{{x}^{\prime}}}^{2}\text{d}z $

(11)

式（10）和（11）中：ρ_f为外部流体密度；A_e和A_i分别为管道外横截面以及内横截面的面积；g为重力加速度；EI为管道弯曲刚度；κ为管道曲率，在小变形振动假设条件下，κ ≈ x″；P_f（z）和P_i（z）分别表示管道外部以及管道内部的流体压强。当管道处在静水中、且为吸流管道时，结构末端会存在负压，进一步导致结构外部和内部存在一定的压力差，但这种压差对于管道动力响应特征影响很小^[12]，因此，为计算方便，这里忽略管道结构的内外压差，将式（11）简化为

$ {V}_{{\mathrm{p}}2}=\frac{1}{2}EI\int \nolimits_{0}^{L}{\kappa }^{2}\text{d}z $

(12)

对于静水中的悬臂输流管道，该系统的非保守力f（z,t）包含两部分，分别为管道周围流体引起的附加质量力f_in （z,t）以及流体阻尼力f_d（z,t），表示为^[13]

$ f\left(z,t\right)={f}_{{\mathrm{in}}}\left(z,t\right)+{f}_{{\mathrm{d}}}\left(z,t\right)={m}_{{\mathrm{a}}}\ddot{x}+\mu {A}_{1}\dot{x}+\frac{1}{2}{A}_{2}{\rho }_\text{f}D\left| \dot{x}\right| \dot{x} $

(13)

式中：m_a为单位长度管道周围的附加流体质量，可表示为：m_a = C_aρ_fA_e，其中，C_a为附加质量系数，通常可取1。μ为周围海水的动力粘性系数，A₁为线性阻尼系数，A₂为非线性阻尼系数，A₁和A₂是根据实验数据拟合得到。由于式（13）中的非保守力（包括：附加质量力以及流体阻尼力）在管道振动过程中均做负功，因此，其所做的功W可表示为

$ W=-\int \nolimits_{0}^{L}f\left(z,t\right)x\text{d}z=-\int \nolimits_{0}^{L}\left({m}_{{\mathrm{a}}}\ddot{x}+\mu {A}_{1}\dot{x}+\frac{1}{2}{A}_{2}{\rho }_\text{f}D\left| \dot{x}\right| \dot{x}\right)x\text{d}z $

(14)

将式（6）~（14）代入式（3），并进行一系列变分简化操作，得到如下表达式

$ \begin{split}&\left[{m}_\text{p}+{m}_\text{e}\delta \left(L-z\right)+{m}_\text{s}+{m}_\text{f}+{m}_\text{a}\right]\ddot{x}-\left(2{m}_\text{s}{U}_\text{s}+2{m}_\text{f}{U}_\text{f}\right)\dot{{x}^{\prime}}+\left({m}_\text{s}U_\text{s}^{2}+{m}_\text{f}U_\text{f}^{2}\right){x}^{\prime\prime}+EI{{{x}^{\prime\prime\prime\prime}}}-\\&\qquad \left({m}_\text{s}g+{m}_\text{f}g+{m}_\text{p}g-{\rho }_\text{f}{A}_\text{e}g\right)\left(L-z\right){x}^{\prime\prime}+\left[{m}_\text{s}g+{m}_\text{f}g+{m}_\text{p}g+{m}_\text{e}g\delta \left(L-z\right)-{\rho }_\text{f}{A}_\text{e}g\right]{x}^{\prime}-\\&\qquad {m}_\text{e}g{x}^{\prime\prime}+\mu {A}_{1}\dot{x}+\frac{1}{2}{A}_{2}{\rho }_\text{f}D\left| \dot{x}\right| \dot{x}=0\end{split} $

(15)

如图2所示，对于上端固支、下端自由的悬臂管道，其边界条件可写作如下形式

$ x(0,t)={x}^{\prime}(0,t)={x}^{\prime\prime}(L,t)={x}^{\prime\prime\prime}(L,t)=0 $

(16)

为合理描述输流管道的内部流动特征，本文对管内颗粒与流体之间的耦合作用模型展开详细讨论。为方便研究，假设固体颗粒为理想的球体。基于上述假设，管内固体颗粒运动速度与管内流体流动速度之间的关系可采用滑移速度模型加以描述。这里，基于Xia等^[31]的实验分析结果对管内多相流体流动特征加以描述，表示为

$ {U}_{{\mathrm{slip}}}=\left[1-{\left({d}_\text{s}/d\right)}^{2}\right]{\omega }_\text{s}{{\mathrm{e}}}^{\left(-2.65{C}_\text{v}+3.32C_\text{v}^{2.2}\right)} $

(17)

式中：U_slip为滑移速度，可表示为：U_slip = U_f − U_s，d_s表示固相颗粒直径，ω_s表示无界流体中单个颗粒的沉降速度，根据斯托克斯方程，可表示为

$ {\omega }_\text{s}=\sqrt{\frac{4g{d}_\text{s}({\rho }_\text{s}-{\rho }_\text{f})}{3{C}_{{\mathrm{Ds}}}{\rho }_\text{f}}} $

(18)

式中：ρ_s表示固相颗粒的密度，C_Ds表示固相颗粒的阻力系数。式（17）中，C_v表示初始进料时管内进料浓度，可表示为

$ {C}_{{\mathrm{v}}}=\frac{{V}_{{\mathrm{so}}}}{{V}_{{\mathrm{so}}}+{V}_{{\mathrm{fo}}}} $

(19)

式中：V_so和V_fo分别表示初始进料时管内固相和液相的所占体积。由于固相颗粒在垂直向上运输过程中，通常会出现堆积现象^[32]，会导致稳定输送时管内的局部浓度C_vl略高于初始进料时管内的进料浓度C_v。为了更准确地描述稳定输送时管内流体的流动特征，这里采用局部浓度对其加以描述，依据Xia等^[31]的实验数据拟合结果，可得局部浓度C_vl和进料浓度C_v之间存在如下关系

$ {C}_\text{vl}=\frac{1}{2}\left(1-\frac{{U}_\text{m}}{{U}_\text{slip}}\right)+{\left[\frac{1}{4}{\left(\frac{{U}_\text{m}}{{U}_\text{slip}}-1\right)}^{2}+{C}_\text{v}\frac{{U}_\text{m}}{{U}_\text{slip}}\right]}^{\frac{1}{2}} $

(20)

式中：U_m表示稳态输送时管内固液两相流的平均速度，可基于据动量守恒定理表示为

$ {U}_\text{m}=\frac{{m}_\text{so}{U}_\text{s}+{m}_\text{fo}{U}_\text{f}}{{m}_\text{so}+{m}_\text{fo}} $

(21)

式中：m_so和m_fo分别表示：初始时刻进入管道内部固相质量以及液相质量。为方便表述固相速度U_s与液相速度U_f之间的关系，引入滑移因子β，表示为

$ \beta =\frac{{U}_\text{s}}{{U}_\text{f}}=\frac{{U}_\text{f}-{U}_\text{slip}}{{U}_\text{f}}=1-\frac{{U}_\text{slip}}{{U}_\text{f}} $

(22)

将式（19）、（21）以及（22）代入式（20）可进一步得到局部浓度C_vl的表达式为

$ {C}_{{\mathrm{vl}}}=\dfrac{1}{2}\left[1-\dfrac{\left(\dfrac{{\rho }_\text{s}}{{\rho }_\text{f}}{C}_\text{v}\beta +1-{C}_\text{v}\right)}{\left(\dfrac{{\rho }_\text{s}}{{\rho }_\text{f}}{C}_\text{v}+1-{C}_\text{v}\right)\left(1-\beta \right)}\right]+{\left[\dfrac{1}{4}{\left(\dfrac{\left(\dfrac{{\rho }_\text{s}}{{\rho }_\text{f}}{C}_\text{v}\beta +1-{C}_\text{v}\right)}{\left(\dfrac{{\rho }_\text{s}}{{\rho }_\text{f}}{C}_\text{v}+1-{C}_\text{v}\right)\left(1-\beta \right)}-1\right)}^{2}+{C}_\text{v}\dfrac{\left(\dfrac{{\rho }_\text{s}}{{\rho }_\text{f}}{C}_\text{v}\beta +1-{C}_\text{v}\right)}{\left(\dfrac{{\rho }_\text{s}}{{\rho }_\text{f}}{C}_\text{v}+1-{C}_\text{v}\right)\left(1-\beta \right)}\right]}^{\frac{1}{2}} $

(23)

基于式（23）求得管内局部浓度C_vl后，便可进一步得到稳定输送时单位长度管内的固相质量和液相质量，表示为

$ {m}_\text{s}=\frac{\text{π} {d}^{2}}{4}{C}_\text{vl}{\rho }_\text{s},\quad {m}_\text{f}=\frac{\text{π} {d}^{2}}{4}\left(1-{C}_\text{vl}\right){\rho }_\text{f} $

(24)

为了更好地观察管道模型的动力学响应，引入如下无量纲变量，对方程（15）进行无量纲化

$ \eta =\frac{x}{D},\xi =\frac{z}{L},\tau =\sqrt{\frac{EI}{{m}_\text{p}}}\frac{t}{{L}^{2}} $

(25)

式中：η为无量纲位移，ξ为无量纲坐标位置，τ为无量纲时间。将公式（25）代入到公式（15），得到如下的无量纲方程

$ \begin{split}&\left[1+\frac{{\kappa }_{5}}{\alpha \chi }\delta \left(1-\xi \right)+\frac{{\kappa }_{1}}{\alpha }+\frac{{\kappa }_{2}}{\alpha }+\frac{{\kappa }_{4}}{\alpha }\right]\frac{{\partial }^{2}\eta }{\partial {\tau }^{2}}-2\left({\kappa }_{1}\beta +{\kappa }_{2}\right)\frac{\chi {u}_\text{f}}{\sqrt{\alpha }}\frac{\partial \eta }{\partial \xi \partial \tau }+\left({\kappa }_{1}{\beta }^{2}+{\kappa }_{2}\right){\chi }^{2}u_\text{f}^{2}\frac{{\partial }^{2}\eta }{\partial {\xi }^{2}}+\\&\qquad \frac{{\partial }^{4}\eta }{\partial {\xi }^{4}}-\left({\kappa }_{1}+{\kappa }_{2}+\alpha -{\kappa }_{3}\right){\chi }^{3}\left(1-\xi \right)\frac{{\partial }^{2}\eta }{\partial {\xi }^{2}}+\left({\kappa }_{1}+{\kappa }_{2}+\alpha -{\kappa }_{3}\right){\chi }^{3}\frac{\partial \eta }{\partial \xi }+\\&\qquad {\kappa }_{5}\delta \left(1-\xi \right){\chi }^{2}\frac{\partial \eta }{\partial \xi }-{\kappa }_{5}{\chi }^{2}\frac{{\partial }^{2}\eta }{\partial {\xi }^{2}}+{A}_{1}\varsigma \frac{\partial \eta }{\partial \tau }+{A}_{2}\gamma \left| \frac{\partial \eta }{\partial \tau }\right| \frac{\partial \eta }{\partial \tau }=0\end{split}$

(26)

式（26）中，各无量纲系数可表示为

$ \begin{split}&{\kappa }_{1}={m}_\text{s}g\frac{{D}^{3}}{EI},\quad {\kappa }_{2}={m}_\text{f}g\frac{{D}^{3}}{EI},\quad {\kappa }_{3}={\rho }_\text{f}{A}_\text{e}g\frac{{D}^{3}}{EI},\quad {\kappa }_{4}={m}_\text{a}g\frac{{D}^{3}}{EI},\quad {\kappa }_{5}={m}_\text{e}g\frac{{D}^{2}}{EI},\quad \alpha ={m}_\text{p}g\frac{{D}^{3}}{EI}\\&{u}_\text{f}=\frac{{U}_\text{f}}{\sqrt{gD}},\quad \chi =\frac{L}{D},\quad \varsigma =\frac{\mu {L}^{2}}{\sqrt{{m}_\text{p}EI}},\quad \gamma =\frac{{\rho }_\text{f}{D}^{2}}{2{m}_\text{p}}\end{split} $

(27)

式中：κ₁为无量纲固相质量，κ₂为无量纲液相质量，κ₃为无量纲液相浮重，κ₄为无量纲附加质量，κ₅为无量纲末端质量，α为无量纲管道质量，u_f为无量纲液相流速，χ为长细比，ς为无量线性阻尼系数，γ为无量纲非线性阻尼系数。

同理，将式（16）所给出的边界条件也写成无量纲形式为

$ \eta (0,\tau )={\eta }^{\prime}(0,\tau )={\eta }^{\prime\prime}(1,\tau )={\eta }^{\prime\prime\prime}(1,\tau )=0 $

(28)

2　数值方法

收起

本节将对悬臂输流管道的稳定性特征以及振动响应特征数值分析方法展开详细介绍。基于Galerkin法，对结构振动无量纲位移η进行离散，并截取前N阶模态振型，表示为

$ \eta \left(\xi ,\tau \right)=\sum \limits_{i=1}^{N}{\phi }_{i}\left(\xi \right){\overline{\eta }}_{i}\left(\tau \right) $

(29)

式中：φ_i（ξ）为悬臂边界条件下输流管道的第i阶模态振型，

$ {\overline{\eta }}_{i}\left(\tau \right) $

是第i阶广义坐标。对于悬臂管道结构，其模态振型可表示为

$ {\phi }_{i}\left(\xi \right)=\cosh \left({\beta }_{i}\xi \right)-\cos \left({\beta }_{i}\xi \right)-{\sigma }_{i}\left[\sinh \left({\beta }_{i}\xi \right)-\sin \left({\beta }_{i}\xi \right)\right],\quad {\sigma }_{{i}}=\frac{\cosh \left({\beta }_{i}\right)+\cos \left({\beta }_{i}\right)}{\sinh \left({\beta }_{i}\right)+\sin \left({\beta }_{i}\right)},{i}=1,\cdots ,N $

(30)

式（30）中，β_i满足如下特征方程

$ \cos \left({\beta }_{i}\right)\cosh \left({\beta }_{i}\right)+1=0 $

(31)

通过求解式（31），可非常方便地得到前8阶模态振型对应的β_i值，分别为1.8751，4.6941，7.8548，10.9955，14.1372，17.2788，20.4204，23.5619。将式（29）代入式（26），方程左右两端同时乘以振型函数φ_i（ξ），并在区间[0, 1]上进行定积分，可以得到如下表达式

$ \begin{split}&\int \nolimits_{0}^{1}{\phi }_{i}\left[1+\frac{{\kappa }_{5}}{\alpha \chi }\delta \left(1-\xi \right)+\frac{{\kappa }_{1}}{\alpha }+\frac{{\kappa }_{2}}{\alpha }+\frac{{\kappa }_{4}}{\alpha }\right]{\phi }_{j}{\ddot{\overline{\eta }}}_{j}\mathrm{d}\xi -\int \nolimits_{0}^{1}{\phi }_{i}2\left({\kappa }_{1}\beta +{\kappa }_{2}\right)\frac{\chi {u}_\text{f}}{\sqrt{\alpha }}{\phi }_{j}^{\prime}{\dot{\overline{\eta }}}_{j}\mathrm{d}\xi+ \\&\qquad \int \nolimits_{0}^{1}{\phi }_{i}\left({\kappa }_{1}{\beta }^{2}+{\kappa }_{2}\right){\chi }^{2}u_\text{f}^{2}{\phi }_{j}^{\prime\prime}{\overline{\eta }}_{j}\mathrm{d}\xi -\left({\kappa }_{1}+{\kappa }_{2}+\alpha -{\kappa }_{3}\right){\chi }^{3}\int \nolimits_{0}^{1}{\phi }_{i}\left(1-\xi \right){\phi }_{j}^{\prime\prime}{\overline{\eta }}_{j}\mathrm{d}\xi+ \\&\qquad \left({\kappa }_{1}+{\kappa }_{2}+\alpha -{\kappa }_{3}\right){\chi }^{3}\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}^{\prime}{\overline{\eta }}_{j}\mathrm{d}\xi +\int \nolimits_{0}^{1}{\phi }_{i}{\kappa }_{5}\delta \left(1-\xi \right){\chi }^{2}{\phi }_{j}^{\prime}{\overline{\eta }}_{j}\mathrm{d}\xi -\\&\qquad \int \nolimits_{0}^{1}{\kappa }_{5}{\chi }^{2}{\phi }_{i}{\phi }_{j}^{\prime\prime}{\overline{\eta }}_{j}\mathrm{d}\xi +\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}{{^{\prime\prime}}}^{\prime\prime}{\overline{\eta }}_{j}\mathrm{d}\xi +\int \nolimits_{0}^{1}{\phi }_{i}{A}_{1}\varsigma {\phi }_{j}{\dot{\overline{\eta }}}_{j}\mathrm{d}\xi +\int \nolimits_{0}^{1}{\phi }_{i}{A}_{2}\gamma \left| {\phi }_{j}{\dot{\overline{\eta }}}_{j}\right| {\phi }_{k}{\dot{\overline{\eta }}}_{k}\mathrm{d}\xi =0\end{split} $

(32)

为了考虑末端集中质量块对结构振型的影响，将式（32）中基于Dirac函数描述的集中质量块函数进一步改写为^[32]

$ \int \nolimits_{0}^{1}\delta \left(1-\xi \right){\phi }_{i}{\phi }_{j}\mathrm{d}\xi ={\phi }_{i}\left(1\right){\phi }_{j}\left(1\right),\quad \int \nolimits_{0}^{1}\delta \left(1-\xi \right){\phi }_{i}{{{\phi }^{\prime}}}_{j}\mathrm{d}\xi ={\phi }_{i}\left(1\right){{{\phi }^{\prime}}}_{j}\left(1\right) $

(33)

将式（33）代入式（32）得到

$ \begin{split}&\left(1+\frac{{\kappa }_{1}}{\alpha }+\frac{{\kappa }_{2}}{\alpha }+\frac{{\kappa }_{4}}{\alpha }\right)\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}{\ddot{\overline{\eta }}}_{j}\mathrm{d}\xi +\frac{{\kappa }_{5}}{\alpha \chi }{\phi }_{i}\left(1\right){\phi }_{j}\left(1\right)\ddot{\overline{\eta }}-2\left({\kappa }_{1}\beta +{\kappa }_{2}\right)\frac{\chi {u}_\text{f}}{\sqrt{\alpha }}\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}^{\prime}{\dot{\overline{\eta }}}_{j}\mathrm{d}\xi -{\kappa }_{5}{\chi }^{2}\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}^{\prime\prime}{\overline{\eta }}_{j}\mathrm{d}\xi+ \\&\qquad \left({\kappa }_{1}{\beta }^{2}+{\kappa }_{2}\right){\chi }^{2}u_\text{f}^{2}\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}^{\prime\prime}{\overline{\eta }}_{j}\mathrm{d}\xi +\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}^{{{\prime\prime\prime\prime}}}{\overline{\eta }}_{j}\mathrm{d}\xi -\left({\kappa }_{1}+{\kappa }_{2}+\alpha -{\kappa }_{3}\right){\chi }^{3}\int \nolimits_{0}^{1}{\phi }_{i}\left(1-\xi \right){\phi }_{j}^{\prime\prime}{\overline{\eta }}_{j}\mathrm{d}\xi+ \\&\qquad \left({\kappa }_{1}+{\kappa }_{2}+\alpha -{\kappa }_{3}\right){\chi }^{3}\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}^{\prime}{\overline{\eta }}_{j}\mathrm{d}\xi +{\kappa }_{5}{\chi }^{2}{\phi }_{i}\left(1\right){\phi }_{j}^{\prime}\left(1\right){\overline{\eta }}_{j}+{A}_{1}{\varsigma }_{1}\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}{\dot{\overline{\eta }}}_{j}\mathrm{d}\xi +{A}_{2}\gamma \int \nolimits_{0}^{1}{\phi }_{i}\left| {\phi }_{j}{\dot{\overline{\eta }}}_{j}\right| {\phi }_{k}{\dot{\overline{\eta }}}_{k}\mathrm{d}\xi =0\end{split} $

(34)

在分析静水中管道的稳定性特征时，非线性阻尼系数A₂对管道颤振后动力影响不可忽视，所以在研究稳定性时，可略去式（34）中的非线性阻尼项^[13]，并将其改写为矩阵形式

$ {\boldsymbol{M}}_{N\times N}{{\ddot{\overline{\boldsymbol{\eta }}}}}_{N\times 1}\left(\tau \right)\boldsymbol+{\boldsymbol{C}}_{N\times N}{{\dot{\overline{\boldsymbol{\eta }}}}}_{N\times 1}\left(\tau \right)+{\boldsymbol{K}}_{N\times N}{{\overline{\boldsymbol{\eta }}}}_{N\times 1}\left(\tau \right)={\boldsymbol{0}}_{N\times 1} $

(35)

式中：质量矩阵M、阻尼矩阵C以及刚度矩阵K中对应的各矩阵元素M_ij、C_ij以及K_ij可表示如下：

$\begin{split}&{M}_{ij}=\left(1+\frac{{\kappa }_{1}}{\alpha }+\frac{{\kappa }_{2}}{\alpha }+\frac{{\kappa }_{4}}{\alpha }\right)\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}\mathrm{d}\xi +\frac{{\kappa }_{5}}{\alpha \chi }{\phi }_{i}\left(1\right){\phi }_{j}\left(1\right),\quad {C}_{ij}=-2\left({\kappa }_{1}\beta +{\kappa }_{2}\right)\frac{\chi {u}_\text{f}}{\sqrt{\alpha }}\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}^{\prime}\mathrm{d}\xi +{A}_{1}\varsigma \int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}\mathrm{d}\xi \\&{K}_{ij}=\left({\kappa }_{1}{\beta }^{2}+{\kappa }_{2}\right){\chi }^{2}u_\text{f}^{2}\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}^{\prime\prime}\mathrm{d}\xi +\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}^{\prime\prime\prime\prime}\mathrm{d}\xi -\left({\kappa }_{1}+{\kappa }_{2}+\alpha -{\kappa }_{3}\right){\chi }^{3}\int \nolimits_{0}^{1}\left(1-\xi \right){\phi }_{i}{\phi }_{j}^{\prime\prime}\mathrm{d}\xi +\\&\qquad \left({\kappa }_{1}+{\kappa }_{2}+\alpha -{\kappa }_{3}\right){\chi }^{3}\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}^{\prime}\mathrm{d}\xi +{\kappa }_{5}{\chi }^{2}{\phi }_{i}\left(1\right){\phi }_{j}^{\prime}\left(1\right)-{\kappa }_{5}{\chi }^{2}\int \nolimits_{0}^{1}{\phi }_{i}{\phi }_{j}^{\prime\prime}\mathrm{d}\xi \end{split}$

(36)

将式（35）转化为一阶微分方程，可表示为

$ {{\dot{\boldsymbol{Z}}}}_{2N\times 1}\left(\tau \right)={\boldsymbol{Y}}_{2N\times 2N}{\boldsymbol{Z}}_{2N\times 1}\left(\tau \right) $

(37)

式中：

$ \boldsymbol{Y}=-{\boldsymbol{B}}^{-1}\boldsymbol{E},\boldsymbol{B}=\left[\begin{matrix}\boldsymbol{0} & \boldsymbol{M}\\\boldsymbol{M} & \boldsymbol{C}\end{matrix}\right],\quad \boldsymbol{E}=\left[\begin{matrix}-\boldsymbol{M} & \boldsymbol{0}\\\boldsymbol{0} & \boldsymbol{K}\end{matrix}\right],\quad \boldsymbol{Z}\left(\tau \right)\boldsymbol=\left[\begin{array}{c}\boldsymbol{\dot{\overline{\boldsymbol{\eta }}}}\left(\tau \right)\\\boldsymbol{\overline{\boldsymbol{\eta }}}\left(\tau \right)\end{array}\right] $

(38)

假设Z（τ）的表达式写为Z（τ） = Ae^λτ，并将其代入式（37），得到

$ \left(\lambda \boldsymbol{I}-\boldsymbol{Y}\right)\boldsymbol{A}=0 $

(39)

式中：I为单位矩阵，由式（39）可进一步得到

$ \det \left(\boldsymbol{Y}-\lambda \boldsymbol{I}\right)=0 $

(40)

由式（40）可看出，通过求解矩阵Y的特征值便可以确定内流对结构稳定性的影响。复特征值λ是一系列复数，复特征值λ与复特征频率ω之间存在如下恒定关系，即：λ ≡ iω。对于特征值问题，复特征频率的实部（或复特征值的虚部）反映系统的固有频率；而复特征频率的虚部（或复特征值的实部）反映系统的振动稳定性特征。若存在复特征频率虚部小于0（或复特征值实部大于0）的情况，则认为结构发生颤振失稳。

基于四阶Houbolt差分格式对方程（35）进行离散，并通过Newton–Raphson迭代法求解所得非线性方程组，从而获得结构的动力响应特征。具体数值方法可参考Gao等^[24]的研究。

3　分析与讨论

收起

3.1　数值方法验证

本节对第2章介绍的内流引起的结构稳定性特征数值计算方法进行验证。值得注意的是，由于缺乏管内为固液两相流的末端带集中质量块的悬臂吸流管道模型的相关数据，因此，无法直接对文中提出的数值模型展开验证，本小节仅对管内为单相流体的输流管道模型展开验证。为充分验证本文所提数值方法的可靠性，本小节依次进行两次对比验证：第一次是将本文的数值模型简化为不带集中质量块的悬臂吸流管道模型，并与Kuiper等的结果^[13]进行对比；第二次是将本文的数值模型简化为带集中质量块的悬臂排流管道模型，并与ElNajjar等的结果^[32]进行对比。

这里首先介绍第一次对比验证。为方便对比，略去方程（15）中的固相质量以及末端集中质量，将其退化为管内为单相流体的悬臂吸流管道模型，表示为

$ \begin{split}&\left({m}_\text{p}+{m}_\text{f}+{m}_\text{a}\right)\ddot{x}-2{m}_\text{f}{U}_\text{f}\dot{{x}^{\prime}}+{m}_\text{f}U_\text{f}^{2}{x}^{\prime\prime}+EI{{{x}^{\prime\prime\prime\prime}}}-\left({m}_\text{f}g+{m}_\text{p}g-{\rho }_\text{f}{A}_\text{e}g\right)\left(L-z\right){x}^{\prime\prime}+\\&\qquad \left({m}_\text{f}g+{m}_\text{p}g-{\rho }_\text{f}{A}_\text{e}g\right){x}^{\prime}+\mu {A}_{1}\dot{x}+\frac{1}{2}{A}_{2}{\rho }_\text{f}D\left| \dot{x}\right| \dot{x}=0\end{split} $

(41)

将方程（41）转化为无量纲形式，令

$ \tilde{\eta }=\frac{x}{L},\quad \xi =\frac{z}{L},\quad \tilde{\tau }=\sqrt{\frac{EI}{M}}\frac{t}{{L}^{2}} $

(42)

式（42）中，为与Kuiper等的^[13]研究保持一致，引入了新的无量纲位移

$ \tilde{\eta } $

以及无量纲时间

$ \tilde{\tau } $

，表达式

$ \tilde{\tau } $

中M为单位长度系统总质量，包括3部分：单位长度管道结构质量m_p、单位长度管道内部流体质量m_f以及单位长度管道外部附加流体质量m_a，即M = m_p+m_f+m_a。将式（42）代入式（41），可得到内流为单相内流时悬臂吸流管道振动无量纲方程，表示为

$ \frac{{\partial }^{4}\tilde{\eta }}{\partial {\xi }^{4}}-2\tilde{\beta }U\frac{{\partial }^{2}\tilde{\eta }}{\partial \xi \partial \tilde{\tau }}+{U}^{2}\frac{{\partial }^{2}\tilde{\eta }}{\partial {\xi }^{2}}+\frac{{\partial }^{2}\tilde{\eta }}{\partial {\tilde{\tau }}^{2}}-\tilde{\alpha }\left(1-\xi \right)\frac{{\partial }^{2}\tilde{\eta }}{\partial {\xi }^{2}}+\tilde{\alpha }\frac{\partial \tilde{\eta }}{\partial \xi }+\tilde{\varsigma }{A}_{1}\frac{\partial \tilde{\eta }}{\partial \tilde{\tau }}+\tilde{\gamma }{A}_{2}\left| \frac{\partial \tilde{\eta }}{\partial \tilde{\tau }}\right| \frac{\partial \tilde{\eta }}{\partial \tilde{\tau }}=0 $

(43)

式中：

$ \tilde{\alpha } $

为无量纲重力系数，

$ \tilde{\beta } $

为无量纲质量比，U为无量纲内流速度，

$ \xi $

为无量纲线性阻尼，

$ \tilde{\gamma } $

为无量纲非线性流体阻尼，这些无量纲系数可表示为

$ \tilde{\alpha }=\frac{\left({m}_\text{p}g+{m}_\text{f}g-{\rho }_\text{f}{A}_\text{e}g\right){L}^{3}}{EI},\quad \tilde{\beta }=\sqrt{\frac{{m}_\text{f}}{M}},\quad U={u}_\text{f}L\sqrt{\frac{{m}_\text{f}}{EI}},\quad \tilde{\varsigma }=\frac{\mu {L}^{2}}{\sqrt{MEI}},\quad \tilde{\gamma }=\frac{{\rho }_\text{f}DL}{2M} $

(44)

基于3.1节提出的结构稳定性数值计算方法，对公式（43）进行数值求解。但Kuiper等的^[13]研究表明，在计算系统稳定性时，可忽略式（43）中的非线性流体阻尼项，即A₂取0。图3展示了基于本文数值计算方法得到的复特征值Argand图，并将其与Kuiper等^[13]研究得到的复特征值Argand图进行了对比。

如图3所示，当内部流速U = 0时，系统的复特征值λ的实部均小于零，表明结构处于稳定状态；随着内部流速U的增加，当U达到临界流速U_cr时，复特征值λ的实部达到零点，结构发生颤振失稳。由图3（a）可看出，当截断模态N取1时，基于本文数值方法计算得到的U_cr与Kuiper得到的U_cr分别为8.91和8.74；由图3（b）可看出，当截断模态N取2时，基于本文数值方法计算得到的U_cr与Kuiper得到的U_cr分别为4.48和4.45。基于以上分析可看出：基于本文数值方法计算得到的内流稳定性特征与文献[13]基本一致。同时可以看出：N = 2时的临界流速仅约为N = 1时临界流速的一半，表明至少需要2阶及以上模态才能可靠地反映其振动稳定性特征。

下面介绍第二次对比验证，将式（15）进行退化并进行无量纲化处理，得到如下表达式

$ \frac{{\partial }^{2}\tilde{\eta }}{\partial {\tau }^{2}}+2\overline{\beta }U\frac{{\partial }^{2}\tilde{\eta }}{\partial \xi \partial \tau }+\left[1+\mu \delta \left(1-{\xi }_\text{p}\right)\right]\frac{{\partial }^{4}\tilde{\eta }}{\partial {\xi }^{4}}+\left[{U}^{2}-\overset{\frown }{\gamma }(1-\xi )\right]\frac{{\partial }^{2}\tilde{\eta }}{\partial {\xi }^{2}}+\overset{\frown }{\gamma }\frac{\partial \tilde{\eta }}{\partial \xi }=0 $

(45)

式（45）中的无量纲变量可表示为

$ \tilde{\eta }=\frac{x}{L},\;\;\;\xi =\frac{z}{L},\;\;\;\overset{\frown }{\tau }=\sqrt{\frac{EI}{{m}_\text{p}+{m}_\text{f}}}\frac{t}{{L}^{2}},\;\;\;U=\sqrt{\frac{{m}_\text{f}}{EI}}L{u}_\text{f},\;\;\;\overline{\beta }=\sqrt{\frac{{m}_\text{f}}{{m}_\text{p}+{m}_\text{f}}},\;\;\;\overset{\frown }{\gamma }=\left(\frac{{m}_\text{p}+{m}_\text{f}}{EI}\right)\frac{{L}^{3}}{g},\;\;\;\mu =\frac{{m}_\text{e}}{\left({m}_\text{p}+{m}_\text{f}\right)L} $

(46)

基于3.1节提出的结构稳定性数值计算方法对公式（45）进行数值求解，并通过特征频率得到管道发生颤振失稳时的临界速度。图4给出了当μ = 0.2、ξ_p = 0.5以及

$ \overset{\frown }{\gamma } $

= 10时，基于本文数值计算方法得到的临界速度与ElNajjar等^[32]的验证结果对比。由图4可看出：临界速度u_cr随

$ \overline{\beta } $

的变化趋势与文献基本一致。

上述两次对比验证结果表明，本文关于“末端带集中质量块时内流激励下悬臂吸流管道”振动稳定性数值计算理论和模型具有可靠性。

3.2　结构振动稳定性特征分析

本节结构参数与Gao等^[24]一致。固相颗粒密度取自Dreiseitl ^[33]对海底结核的物理特征分析结果，为3500 kg/m³。流体阻尼系数A₁与A₂的取值与Kuiper等^[13]研究中的取值保持一致。固相颗粒阻力系数C_Ds参考Zhou等的研究取值^[34]。本文数值计算所采用的主要物理参数如表1所示。

为充分认识粗颗粒固液两相内流流动引起的悬臂管道动力学响应特征，本小节首先对粗颗粒固液两相流在管内的流动特征展开研究。基于前面部分提出的滑移模型，本节集中讨论了处于稳态输送阶段时管道无量纲内部流速以及进料浓度参数对滑移因子以及局部浓度分布的影响。本节所采用的全部数值计算参数均基于表1中得到，图5（a）给出了处于稳定输送阶段时不同进料浓度C_v下滑移因子β随无量纲内流速度u_f的变化关系。当内流速度保持不变时，滑移因子β随进料浓度C_v的增加而整体呈上升趋势，但当C_v超过0.3后，该增长趋势趋于缓和。当进料浓度保持不变时，滑移因子β随u_f的增大而单调上升。需要特别指出的是，在u_f ≤0.67的低速阶段，滑移因子会出现负值，表明固相颗粒未能随液相上升，反而在管道内发生沉降。随着u_f的增加，当u_f>0.67时，滑移因子迅速增大，并逐渐趋近于1，这表明液相对固相的携带能力显著增强。

图5（b）展示了处于稳定输送阶段时不同进料浓度C_v下管内局部浓度C_vl随无量纲内流速度u_f的变化趋势。结果表明，在低速阶段（0 < u_f ≤ 1），局部浓度随u_f的减小而迅速升高。这一现象可解释为：当液相流速较低时，其携带固相颗粒的能力减弱，颗粒更易在管内沉积和聚集，从而导致局部浓度显著升高。随着u_f的进一步增加（u_f>1），局部浓度C_vl逐渐接近进料浓度C_v，表明固相颗粒运动趋于均匀分布，输送过程趋于稳定。

图6给出了进料浓度C_v 分别为0.1和0.3时，随无量纲内流速度变化的管道结构复特征频率Argand图。图中横轴和纵轴分别表示复特征频率的实部和虚部。实部表示管道固有频率；虚部表示振动时管道的阻尼特征。系统阻尼比ζ通常可定义为：ζ= Im（ω）/Re（ω），这意味着：当Argand图的其中一个分支穿过横轴（即Im（ω）=0）时，系统的复特征频率虚部Im（ω）从正值越过0点变为负值，此时结构系统阻尼比小于0，结构将从内流中持续获得能量，进而发生颤振失稳。这种系统阻尼比从正值跨过0点达到负值引起的颤振失稳模式，通常被称为Hopf分岔失稳，发生Hopf分岔的最小流速，记为临界流速u_cr。

从图6可以看到，在进料浓度C_v= 0.1时，管道结构前三阶模态呈现出Hopf分岔现象，且每阶模态对应的颤振失稳流速u_f不同。图6 （a）表明：当C_v = 0.1时，管道前三阶模态对应的临界流速u_f分别3.43、6.51、12.39。值得注意的是，图6 （b）中，当进料浓度C_v增加到0.3时，第四阶模态呈现出Hopf分岔现象，前四阶模态对应的临界流速u_f分别为2.58、4.43、7.78、15.38。可看出：随进料浓度C_v的提高，会导致高阶模态更易被激发，从而在更高阶数的模态中也出现了颤振失稳。且对于两种不同的进料浓度C_v条件下，管道结构最先在第一阶模态下发生颤振失稳，说明一阶模态的动态特性对系统稳定性具有主导作用。此外，当C_v = 0.1以及0.3时，管道最先开始失稳临界流速u_cr分别为3.43和2.58，可见，随着进料浓度C_v的增加，临界流速u_cr呈明显下降趋势。这一现象主要是因为：进料浓度升高所引起的管内固相密度与总质量的增加，进而增强了系统内部流动所产生的离心力与科里奥利力。其中，离心力增大会削弱管道的结构刚度，而科里奥利力作为非保守激励源，在吸流条件下会增强系统的能量输入，导致结构更易进入颤振状态^[13]。

为了确保结构稳定性特征计算方法的可靠性，基于Newton-Raphson 迭代法，图7分别对结构内流速度处于临界速度u_cr前的亚临界内流速度和临界速度u_cr后的超临界内流速度对应的结构动力学特征展开了分析。图7分别给出了图6（b）中进料浓度C_v=0.3时，管道处于亚临界内流速度（u_f=2.57）以及临界内流速度（u_cr=2.58）对应的结构末端处振动位移。由图7可以看出，当内流速度处于亚临界区间时，结构末端位移时程曲线均随时间而逐渐衰减，此时结构振动形式为衰减振动；当内流速度到达临界流速时，结构末端位移时程曲线均呈现出随时间增加而逐渐放大的特征，此时结构振动呈现出颤振失稳形式。

图8给出了四种不同的无量纲端部质量块κ₅下系统临界流速u_cr与管道质量α的变化关系图。由图可以看出，当α保持不变时，随着κ₅的增大，临界流速呈现逐渐上升的趋势，这是因为末端质量的增大在提高系统整体惯性的同时，也增强了管道结构的有效刚度，从而提升了结构对流动诱导失稳的抵抗能力。当κ₅保持不变时，随着α增大，临界流速亦随之升高，呈现近似抛物线的增长趋势。此外，当α较大时，κ₅对u_cr的影响显著增大。这一现象表明，较大管道质量提高了系统的刚度，显著削弱科里奥利力的影响，从而提升管道的临界流速。

4　结论

收起

本文基于修正的哈密顿原理建立了多相内流激励下垂直吸流管道非线性振动控制方程，使用四阶Galerkin法对方程进行了时空离散，得到了常微分方程组。将常微分方程组写成矩阵形式，基于特征值法研究了管道的稳定性特征，并利用Newton – Raphson迭代法验证了临界流速。本文针对不同无量纲液相内部流速u_f、进料浓度C_v、末端质量块ĸ₅以及管道质量α下管道振动稳定性特征进行了系统分析。基于以上分析结果，得到如下结论：

（1）当内流速度处于亚临界区间时，系统特征频率虚部均为正值，结构无法从流体中获取能量，其振动形式为衰减振动；当内流速度处于超临界区间时，系统的一阶特征频率虚部会出现负值，结构会从流体中持续不断地获取能量，结构发生颤振失稳。

（2）当进料浓度较低时，结构仅在前三阶模态发生Hopf分岔，系统的临界流速较高，失稳主要集中在前三阶模态；当进料浓度升高至C_v= 0.3时，第四阶模态也发生Hopf分岔，且最小临界流速u_cr显著降低；进料浓度升高导致管内固相质量与总质量增加，加强了科里奥利力的作用，使系统更易从内流中获取能量，进而进入颤振状态。

（3）当管道质量 α保持不变时，系统临界流速u_cr随无量纲末端质量κ₅的增大而逐渐上升，表明末端质量增强了系统惯性与结构刚度；当 κ₅保持不变时，随着管道质量α 的增大，临界流速显著升高，呈现近似抛物线的增长趋势，较大管道质量可有效削弱科里奥利力的非保守激励效应，从而提升系统稳定性。

由前面部分的分析可以看出，固液两相流滑移模型可以很好地反映细长柔性采矿管道的稳定性。现简要阐述固液两相流滑移模型的优缺点。与CFD方法和实验方法相比，固液两相流滑移模型最大的优点在于计算时间短、消耗成本低，本文所有工况的计算时间均在一个小时内完成。因此固液两相流滑移模型更适用于计算实际海洋工程中细长柔性采矿管道的稳定性特性的多工况问题。但是鉴于模型自身的简化性，固液两相流滑移模型同样存在不足之处，其中最大的不足便是计算精度问题，固液两相流滑移模型是一种基于试验参数的半经验模型方法，其计算精度受到很多因素制约，如固相颗粒阻力系数的选取、固相颗粒形状、固相颗粒与管壁的碰撞效应等。如何提高固液两相流滑移模型的计算精度使其更好地适用于实际海洋工程中，这是一个具有实际工程价值的问题，后续将对该问题进行更深入的研究，并开展相应的实验以验证该数值模型的准确性。

基金

收起

国家重点研发计划资助项目（2024YFC2814803）；国家自然科学基金资助项目（42476206）；山东省自然科学基金资助项目（ZR2025MS759）；上海交通大学海洋工程全国重点实验室开放基金资助项目（GKZD010089）.

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

Sharma R. Deep-sea mining: Current status and future considerations[M]. Cham: Springer International Publishing, 2017.

徐万海, 谢武德, 高喜峰, 等. 海底多段悬跨管道涡激振动特性分析[J]. 船舶力学, 2017, 166(8): 1025‒1034.

Xu W H, Xie W D, Gao X F, et al. Study on vortex-induced vibrations (VIV) of multi-spans pipelines[J]. Journal of Ship Mechanics, 2017, 166(8): 1025‒1034. (in Chinese)

高云, 潘港辉, 刘磊, 等. 变张力细长柔性立管涡激振动响应特性数值研究[J]. 船舶力学, 2023, 27(9): 1398‒1412.

Gao Y, Pang G H, Liu L, et al. Numerical study on responses of vortex-induced vibration of a long flexible riser with an axial time-varying tension[J]. Journal of Ship Mechanics, 2023, 27(9): 1398‒1412. (in Chinese)

Gao Y, Zhang Z Z, Pan G H, et al. Three-dimensional vortex-induced vibrations of a circular cylinder predicted using a wake oscillator model[J]. Marine Structures, 2021, 80(1038): 103078.

陈正寿, 王帅. 内外流耦合作用下弹性管多模态振动响应研究[J]. 船舶力学, 2025, 29(2): 219‒231.

Chen Z S, Wang S. Multi-mode vibration response of flexible pipe under coupled effect of external current and internal flow[J]. Journal of Ship Mechanics, 2025, 29(2): 219‒231. (in Chinese)

Meng S, Song S D, Che C D, et al. Internal flow effect on the parametric instability of deepwater drilling risers[J]. Ocean Engineering, 2018, 149: 305-312.

Gregory R W, Paidoussis M P. Unstable oscillation of tubular cantilevers conveying fluid I. Theory[J]. Mathematical and Physical Sciences, 1966, 293(1435): 528-542.

Païdoussis M P, Issid N T. Dynamic stability of pipes conveying fluid[J]. Journal of Sound and Vibration, 1974, 33(3): 267-294.

Noah S T, Hopkins G R. Dynamic stability of elastically supported pipes conveying pulsating fluid[J]. Journal of Sound and Vibration, 1980, 71(1): 103-116.

Paidoussis M P, Luu T P. Dynamics of a pipe aspirating fluid such as might be used in ocean mining[J]. Journal of Energy Resources Technology, 1985, 107: 250-255.

Païdoussis M P, Li G X, Moon F C. Chaotic oscillations of the autonomous system of a constrained pipe conveying fluid[J]. Journal of Sound and Vibration, 1989, 135(1): 1-19.

Kuiper G L, Metrikine A V. Dynamic stability of a submerged, free-hanging riser conveying fluid[J]. Journal of Sound and Vibration, 2005, 280(3-5): 1051-1065.

Kuiper G L, Metrikine A V, Battjes J A. A new time-domain drag description and its influence on the dynamic behaviour of a cantilever pipe conveying fluid[J]. Journal of Fluids and Structures, 2007, 23(3): 429-445.

王乙坤, 王琳. 分布式运动约束下悬臂输液管的参数共振研究[J]. 力学学报, 2019, 51(2): 558‒568.

Wang Y K, Wang L. Parametric resonance study of cantilever infusion tube under distributed motion constraints[J]. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(2): 558‒568. (in Chinese)

方孟孟, 郭长青. 悬臂输流管道在基础激励与脉动内流联合作用下的参激振动[J]. 应用力学学报, 2020, 37(2): 653‒660+933−934.

Fang M M, Guo C Q. Parametric vibration of cantilever pipeline under the combined action of foundation excitation and pulsating internal flow[J]. Chinese Journal of Applied Mechanics, 2020, 37(2): 653‒660+933−934. (in Chinese)

柳博瀚, 陈正寿, 鲍健, 等. 管道内流对海洋弹性管振动影响的数值仿真研究[J]. 振动与冲击, 2020, 39(17): 177‒185+202.

Liu B H, Chen Z S, Bao J, et al. Numerical simulation for effects of pipeline internal flow on vibration of flexible marine pipe[J]. Journal of Vibration and Shock, 2020, 39(17): 177‒185+202. (in Chinese)

赵志贤, 随岁寒, 杨昌锦, 等. 基于三阶剪切变形理论的悬臂输流管道自由振动[J]. 力学季刊, 2022, 43(1): 132‒140.

Zhao Z X, Sui S H, Yang C J, et al. Free vibration of cantilever pipeline based on third-order shear deformation theory[J]. Chinese Quarterly of Mechanics, 2022, 43(1): 132‒140. (in Chinese)

漆发辉, 王天林, 郭长青. 分布随从力作用下黏弹性悬臂输流管道的稳定性分析[J]. 南华大学学报(自然科学版), 2022, 36(6): 1‒7.

Qi F H, Wang T L, Guo C Q. Stability analysis of viscoelastic cantilever pipeline under distributed follower forces[J]. Journal of University of South China(Science and Technology), 2022, 36(6): 1‒7. (in Chinese)

Monette C, Pettigrew M J. Fluidelastic instability of flexible tubes subjected to two-phase internal flow[J]. Journal of Fluids and Structures, 2004, 19(7): 943-956.

Zhang M M, Xu J Z. Effect of internal bubbly flow on pipe vibrations[J]. Science China Technological Sciences, 2010, 53: 423-428.

An C, Su J. Dynamic behavior of pipes conveying gas–liquid two-phase flow[J]. Nuclear Engineering and Design, 2015, 292: 204-212.

Ebrahimi-mamaghani A, Sotudeh-gharebagh R, Zarghami R, et al. Dynamics of two-phase flow in vertical pipes[J]. Journal of Fluids and Structures, 2019, 87: 150-173.

王予琪, 宿向辉, 朱祖超. 深海采矿混输泵内流场及粗颗粒运动特性[J]. 排灌机械工程学报, 2022, 40(8): 800‒806.

Wang Y Q, Shu X H, Zhu Z C. Characteristics of flow field and coarse particle motion in multiphase pump for deep-sea mining[J]. Journal of Drainage and Irrigation Machinery Engineering, 2022, 40(8): 800‒806. (in Chinese)

Gao Y, Liu L, Fu, S X, et al. Nonlinear dynamics of a vertical pipe subjected to a two-phase, solid-liquid internal flow[J]. Applied Mathematical Modelling, 2023, 120: 651-666.

高云, 陈晓东, 程玮, 等. 固−液两相内流激励下悬臂输流管道稳定性特征研究[J]. 力学学报, 2024, 56(3): 577‒585.

Gao Y, Chen X D, Cheng W, et al. Study on stability characteristics of cantilever pipeline subjected to a solid-liquid two-phase internal flow[J]. Chinese Journal of Theoretical and Applied Mechanics, 2024, 56(3): 577‒585. (in Chinese)

Su T, Chen S, Yuan H. Research on the characteristics of solid–liquid two-phase flow in the lifting pipeline of seabed mining[J]. Journal of Marine Science and Engineering, 2024, 12(8): 1409.

Zhu H, Zhou M, Gao Y, et al. The vibration of deep-sea mining riser excited by both internal and external flows[J]. Physics of Fluids, 2024, 36(9): 095127.

Benjamin T B. Dynamics of a system of articulated pipes conveying fluid-II. Experiments[J]. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 1961, 261: 487-499.

Benjamin T B. Dynamics of a system of articulated pipes conveying fluid - I. Theory[J]. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 1962, 261: 457-486.

Païdoussis M. P, Pipes conveying fluid: A model dynamicall problem[J]. Journal of Fluids and Structures, 1993, 7(2): 137-204.

Xia J X, Ni J R, Mendoza C. Hydraulic lifting of manganese nodules through a riser[J]. Journal of Offshore Mechanics and Arctic Engineering, 2004, 126(1): 72-77.

Elnajjar J, Daneshmand F. Stability of horizontal and vertical pipes conveying fluid under the effects of additional point masses and springs[J]. Ocean Engineering, 2020, 206: 106943.

Dreiseitl I. About geotechnical properties of the deep seabed polymetallic nodules[C]//18th International Conference on Transport and Sedimentation of Solid Particles, 2017.

Zhou D C, Ma Y Q, Li S Y, et al. Vibration stability of a lifting pipe transporting solid–liquid two phase flow[J]. Ocean Engineering, 2023, 288: 116170.

2026年第30卷第8期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.3969/j.issn.1007-7294.2026.08.006

接收时间：2026-02-25
首发时间：2026-11-03
出版时间：2026-08-15

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2026-02-25

基金

国家重点研发计划资助项目（2024YFC2814803）；国家自然科学基金资助项目（42476206）；山东省自然科学基金资助项目（ZR2025MS759）；上海交通大学海洋工程全国重点实验室开放基金资助项目（GKZD010089）.

作者信息

¹哈尔滨工业大学（威海）海洋工程学院，山东威海 264209

²上海交通大学海洋工程全国重点实验室，上海 200240

参考文献

分享链接

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

分享至

全文二维码

扫描看全文

引用本文

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