船舶力学

网格数	$ {C}_{\mathrm{d}} $（Ave）	St	St 误差^[24]（与文献[23]比较）
73500	1.31	0.231	13.8%
104084	1.49	0.212	4.4%
147200	1.51	0.209	2.9%

网格数	$ {C}_{\mathrm{d}} $（Ave）	St	St 误差^[24]（与文献[23]比较）
73500	1.31	0.231	13.8%
104084	1.49	0.212	4.4%
147200	1.51	0.209	2.9%

涡激振动冻结边界法研究

PDF下载

袁国辉 , 洪方文 , 许影博

船舶力学 | 流体力学 2026,30(4): 520-531

收起

船舶力学 |流体力学 2026 , 30 (4) : 520 -531

涡激振动冻结边界法研究

全屏

[Author(id=1281204008085865238, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, orderNo=0, firstName=null, middleName=null, lastName=null, nameCn=null, orcid=null, stid=null, country=null, authorPic=null, dead=0, email=1470307651@qq.com, emailSecond=null, emailThird=null, correspondingAuthor=1, authorType=1, ext={EN=AuthorExt(id=1281204008379466520, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, authorId=1281204008085865238, language=EN, stringName=Guo-hui YUAN, firstName=Guo-hui, middleName=null, lastName=YUAN, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=null, address=^1.China Ship Scientific Research Center, Wuxi 214082, China, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null), CN=AuthorExt(id=1281204008471741209, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, authorId=1281204008085865238, language=CN, stringName=袁国辉, firstName=null, middleName=null, lastName=null, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=null, address=^1.中国船舶科学研究中心，江苏无锡 214082, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null)}, companyList=[AuthorCompany(id=1281204007997784850, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, xref=1., ext=[AuthorCompanyExt(id=1281204008006173459, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, companyId=1281204007997784850, language=EN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=^1.China Ship Scientific Research Center, Wuxi 214082, China), AuthorCompanyExt(id=1281204008018756372, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, companyId=1281204007997784850, language=CN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=^1.中国船舶科学研究中心，江苏无锡 214082)])]), Author(id=1281204008534655771, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, orderNo=1, firstName=null, middleName=null, lastName=null, nameCn=null, orcid=null, stid=null, country=null, authorPic=null, dead=0, email=null, emailSecond=null, emailThird=null, correspondingAuthor=0, authorType=1, ext={EN=AuthorExt(id=1281204008643707677, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, authorId=1281204008534655771, language=EN, stringName=Fang-wen HONG, firstName=Fang-wen, middleName=null, lastName=HONG, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=null, address=^1.China Ship Scientific Research Center, Wuxi 214082, China, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null), CN=AuthorExt(id=1281204008807285534, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, authorId=1281204008534655771, language=CN, stringName=洪方文, firstName=null, middleName=null, lastName=null, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=null, address=^1.中国船舶科学研究中心，江苏无锡 214082, bio={"content":"

洪方文（1971–），男，研究员

"}, bioImg=null, bioContent=

洪方文（1971–），男，研究员

, aboutCorrespAuthor=null)}, companyList=[AuthorCompany(id=1281204007997784850, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, xref=1., ext=[AuthorCompanyExt(id=1281204008006173459, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, companyId=1281204007997784850, language=EN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=^1.China Ship Scientific Research Center, Wuxi 214082, China), AuthorCompanyExt(id=1281204008018756372, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, companyId=1281204007997784850, language=CN, country=null, province=null, city=null, postcode=null, companyName=null, departmentName=null, remark=^1.中国船舶科学研究中心，江苏无锡 214082)])]), Author(id=1281204008899560224, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, orderNo=2, firstName=null, middleName=null, lastName=null, nameCn=null, orcid=null, stid=null, country=null, authorPic=null, dead=0, email=null, emailSecond=null, emailThird=null, correspondingAuthor=0, authorType=1, ext={EN=AuthorExt(id=1281204009008612130, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, authorId=1281204008899560224, language=EN, stringName=Ying-bo XU, firstName=Ying-bo, middleName=null, lastName=XU, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=null, address=^1.China Ship Scientific Research Center, Wuxi 214082, China, bio=null, bioImg=null, bioContent=null, aboutCorrespAuthor=null), CN=AuthorExt(id=1281204009071526691, tenantId=1146029695717560320, journalId=1240685776644648972, articleId=1281203945318105703, authorId=1281204008899560224, language=CN, stringName=许影博, firstName=null, middleName=null, lastName=null, prefix=null, suffix=null, authorComment=null, nameInitials=null, affiliation=null, department=null, xref=null, address=^1.中国船舶科学研究中心，江苏无锡 214082, bio={"content":"

许影博（1985–），男，研究员

"}, bioImg=null, bioContent=

许影博（1985–），男，研究员

袁国辉, 洪方文, 许影博

作者信息

^1.中国船舶科学研究中心，江苏无锡 214082

通讯作者:

袁国辉（2000–），男，硕士研究生，通讯作者，E-mail：1470307651@qq.com

作者简介:

洪方文（1971–），男，研究员

许影博（1985–），男，研究员

Research on the frozen boundary method for vortex-induced vibration

Guo-hui YUAN, Fang-wen HONG, Ying-bo XU

Affiliations

^1.China Ship Scientific Research Center, Wuxi 214082, China

出版时间: 2026-04-15 doi: 10.3969/j.issn.1007-7294.2026.04.002

文章导航

摘要

收起

自然界和工程实际中普遍存在着振动幅度小、流固耦合效应强的流激振动现象。本文提出了一种边界固定不动，利用质量源和动量源表达边界振动效应的冻结边界法，并利用该方法实现了圆柱单自由度强迫振动和涡激振动计算。结果表明，在强迫振动中，冻结边界法相较于动网格方法，在保证了计算精度的前提下提高了计算速率，验证了该方法的可靠性。对于单自由度圆柱涡激振动成功计算出了锁频现象，由于锁频区共振的影响，圆柱尾流场横向波动较大，尾涡锁频是振动涡系和脱落涡系竞争的结果，当振动涡系主导时表现为频率的锁定，此时发生共振。冻结边界法为计算流固耦合问题提供了一种新的思路和实现途径。

关键词

涡激振动 / 锁频 / 流固耦合 / 数值模拟 / 尾涡脱落

Abstract

收起

The phenomenon of flow-induced vibrations with small amplitude and strong fluid-structure interaction is commonly found in both nature and engineering practice. This paper proposes a frozen boundary method, which keeps the boundary fixed and represents the boundary vibration effects using mass sources and momentum sources. This method is applied to calculate forced vibrations and vortex-induced vibrations of a circular cylinder with a single degree of freedom. The results show that in forced vibrations, the frozen boundary method improves computation speed compared to the dynamic mesh method while ensuring calculation accuracy, thus validating the reliability of the method. For the single degree of freedom vortex-induced vibration of the circular cylinder, the phenomenon of lock-in was successfully computed. Due to the resonance effects within the lock-in region, the lateral fluctuations of the cylinder's wake field are considerable. The wake-vortex lock-in results from the competition between the vibrating vortex system and the detached vortex system. When the vibrating vortex system dominates, it manifests as frequency locking, leading to resonance. The frozen boundary method provides a new perspective and implementation approach for calculating fluid-structure interaction problems.

Key words

vortex-induced vibration / lock-in / fluid-structure interaction / numerical simulation / wake vortex shedding

引用本文

袁国辉, 洪方文, 许影博. 涡激振动冻结边界法研究. 船舶力学, 2026 , 30 (4) : 520 -531 . DOI: 10.3969/j.issn.1007-7294.2026.04.002

Guo-hui YUAN, Fang-wen HONG, Ying-bo XU. Research on the frozen boundary method for vortex-induced vibration[J]. Journal of Ship Mechanics, 2026 , 30 (4) : 520 -531 . DOI: 10.3969/j.issn.1007-7294.2026.04.002

正文

收起

0　引　言

收起

涡激振动（Vortex-inducted vibration，VIV）现象广泛存在于海洋工程^[1]、道桥工程^[2]和船舶工程^[3]中，这种振动的激励源来自于物体后缘交替脱落的涡流。一方面，如果漩涡脱落的频率接近结构的固有频率，则可能会导致共振，造成结构疲劳破坏或引发强烈的噪声问题^[4]。另一方面，近年来相关学者将涡激振动现象应用到了海流能发电领域，将流过装置的水动能转化为振子机械能，进而转换为可收集的电能^[5]。因此，无论是从预防涡激振动带来的危害，还是从加深对其理解进而推广其应用的角度来看，对涡激振动现象特性的研究都尤为重要。

随着计算机技术的不断发展，数值模拟已经成为重要的研究手段。目前用于涡激振动计算的数值模拟方法主要是分区式流固耦合法^[6]，即将流体域和固体域分区进行迭代计算，在流固交界面上进行位移和力等计算信息的交换。从流固耦合界面的处理方式来看，又分为有网格方法^[7−8]和无网格方法^[9−10]。无网格方法是基于拉格朗日思想，将计算域用一系列分布的节点进行离散，在涡激振动计算方面，比较典型的无网格方法有离散涡方法^[11]和光滑粒子流体动力学（Smoothed Particle Hydrodynamic, SPH）方法^[12]。离散涡方法通过在流场中布置离散的涡元来模拟尾流结构的主要特征，并利用涡元的诱导作用获得流场中任意位置处的物理量，进而模拟流体的运动过程。离散涡方法已经发展得较为成熟，它的优点在于计算结果从物理上易于理解和解释，但由于该方法的复杂性，很多基于二维的离散涡模型都不适用，并且因为涉及多体相互作用带来的庞大计算量，三维涡方法的应用也较少。在二维和三维计算中，由于涡元数量随时间不断增加，计算量也随之大幅增加^[13]。SPH方法则使用一组粒子离散和代表所模拟的介质（流体或固体），并且基于粒子体系近似和估算介质运动的控制方程，每个粒子代表宏观连续区域的介质团，伴随介质的运动以当地速度移动，SPH方法最主要的特点为拉格朗日特性和无网格粒子法特性，这使得它在自由面流动^[14−16]计算中具有显著优势，而在非自由面流动问题中其应用相对较少^[17]，但由于采用移动的粒子系统代表所模拟的连续介质，也给边界的处理和边界条件的实施带来了困难。在考虑边界层影响的流动计算时，为精确计算边界层内的粘性力，需要在边界布置高分辨率粒子，这会导致计算量急剧上升。此外，由于边界层内速度梯度较大，在使用SPH方法计算的过程中易出现粒子分布不规则带来的数值噪声甚至计算失败^[18]。

从计算问题的适用性和准确性以及方法推广应用的泛化性来说，目前涡激振动计算无网格方法应用得较少，采用有网格方法比较普遍。有网格方法基于欧拉描述，能够准确描述移动边界，通过网格变形重构，保持与固体边界的贴合，能够精确捕捉边界层，进而有效模拟物体周围流场特征，目前已经发展得较为成熟，并且在涡激振动计算方面得到了广泛的应用^[19−21]。通常使用动网格^[7]或重叠网格^[8]来计算涡激振动问题，即保证流体网格始终与运动中的固体界面匹配重合，通过交界面完成流体物理量与固体物理量之间的数据传递。在网格运动的过程中，网格需要变形或与其他网格进行信息插值，而在涡激振动问题中为了保证计算的精确性，在壁面处需要进行网格加密，这导致在其变形过程中网格容易畸变甚至产生负体积网格，同时网格的变形计算也消耗大量的计算资源。针对网格变形带来的问题，还有学者采用浸没边界法对于涡激振动问题进行计算^[22]，浸没边界法将流体与固体在接触边界上的作用等效为额外的体积力项施加在动量方程中，从而避免了计算网格的变形，但该方法的不足之处在于对耦合界面捕捉不够精确和通过插值施加力时容易带来误差^[6]。

对于螺旋桨涡激振动问题，虽然振幅只有微米量级，但它会导致唱音的发生^[4,23]，严重影响舰船隐蔽性。此外，螺旋桨工作过程中其表面雷诺数较高，因此将动网格方法应用在螺旋桨的涡激振动上时，会出现计算不稳定和效率低的问题。为了避免动网格方法中网格移动带来的网格畸变与计算量增大，以及浸没边界法中固定网格带来的流固耦合界面捕捉以及边界层流动计算不准确的问题，本文提出一种新方法：在固定欧拉网格体系下，通过在流固耦合界面区域构建动态源项，将固体边界的运动效应及流场的非定常动态特征，直接嵌入动量方程与连续性方程的源项修正中。在计算过程中根据物体的运动状态调整源项，实现固体边界条件的等效施加，完成双向流固耦合计算。对于小振幅振动场景（振动幅值远小于特征长度尺度），无需进行网格拓扑调整或界面重构，即可精准表征流固双向强耦合效应。相较于传统固定网格方法，该源项建模策略避免了界面插值带来的精度损失，并且通过边界层网格加密，显著提升了近壁区流动细节的分辨率，同时保留了固定网格的计算高效性，规避了动网格的畸变风险与额外计算成本。为验证该方法的有效性，本文选取圆柱单自由度强迫振动与涡激振动锁频现象作为基准算例，系统模拟了不同折合速度下的流场结构演化、结构振动响应及锁频区间特征，验证了该方法在流场分离涡捕捉和锁频机制复现方面的准确性与可靠性。

1　数理问题描述

收起

涡激振动指的是流动中的物体表面发生流动分离产生涡脱落，脱落的涡在物体上产生非定常作用激励物体振动，条件适当的情况下，物体振动会反过来影响涡脱落，从而形成流固耦合作用，如图1所示。

本文中计算的是单自由度刚体圆柱位于来流速度为U的流场中，其表面边界层在尾流场中会形成交替脱落的涡流，从而产生周期性的力作用在结构上引起振动。

对于涡激振动系统，我们采用分区式流固耦合方法进行描述，则有以下函数关系

$ \left\{\begin{aligned}& {\mathrm{F}}\left(\boldsymbol{u},\boldsymbol{p},t;{S}_{{\mathrm{m}}},{S}_{{\mathrm{M}}}\right)=0\\ &{\mathrm{S}}\left(\boldsymbol{x},\mathit{\Theta },t\right)=0\\&{{\mathrm{F}}}_{\mathit{\Omega }}={{\mathrm{S}}}_{\mathit{\Omega }}\end{aligned}\right. $

(1)

其中，F表示流体域，S_m、S_M分别为用于修正流体域的质量源项和动量源项，S表示固体域，Ω表示流体域和固体域的边界，两者保持速度、压力位移条件的匹配。在使用动网格方法对该问题进行计算时，源项为0，需要在每一个时间步调整流固耦合界面的网格使其满足位移匹配条件，对于大振幅问题这是必须考虑的。但对于小振幅流固耦合问题，如流激噪声，其振动幅度很小，仅在微米量级，此时在模拟过程中仍然选择动网格方法进行计算往往得不偿失，将有限的计算资源大量消耗在了微幅的变形模拟上，同时也可能对计算过程中的网格质量产生影响。因此在冻结边界法中，我们在考虑流固耦合界面满足的关系时，仅考虑速度和压力匹配条件，即通过速度和压力匹配表示流固耦合界面的物理量关系，其表达式如下

$ {\boldsymbol{U}}_{{\mathrm{f}}}\left(t+\mathrm{d}t\right)={\boldsymbol{U}}_{{\mathrm{S}}}\left(t+\mathrm{d}t\right) $

(2)

$ {p}_{{\mathrm{f}}}\left(t+\mathrm{d}t\right)={p}_{{\mathrm{S}}}\left(t+\mathrm{d}t\right) $

(3)

将这两个匹配条件直接转化为流体控制方程中的动量源项与质量源项，嵌入到式（1）中。对于流固界面的速度匹配条件式（2），冻结边界法中网格被“冻结”，因此需通过动量源项

$ {{S}}_{\mathrm{M}} $

补偿“流体未随固体运动”的动量差，使流体速度强制满足式（2）的匹配约束。定义流固界面处流体速度的目标值（固体速度）与计算值的偏差为

$ \Delta u={u}_{{\mathrm{s}}}\left(t+\mathrm{d}t\right)-u_{{\mathrm{f}}}^{*}\left(t+\mathrm{d}t\right) $

(4)

其中，

$ u_{{\mathrm{f}}}^{*} $

为无源性项时流体的预测速度，则动量源项满足

${{S}}_{\mathrm{M}}=\rho \frac{\Delta u}{\mathrm{d}t},\left(x\in \mathit{\Omega }\right) $

(5)

S_M的物理意义是在固定网格下，为了使流体在界面处的速度强制跟随固体运动，需要向流体单元注入的动量率，直接体现了式（2）的速度匹配约束，并嵌入流体动量方程中。

对于流固界面的压力匹配条件式（3），其本质是流固界面处质量守恒的宏观体现。在冻结边界法中，固体边界的微幅运动会引发界面处流体的质量重分布，而固定网格无法通过几何变形来适应这种质量变化，因此需要引入质量源项S_m来补偿该效应，使流体压力场满足式（3）的匹配约束。定义流固界面处流体压力的目标值（固体界面压力）与无源性项时的预测值偏差为

$ \Delta p={p}_{{\mathrm{s}}}\left(t+{\mathrm{d}}t\right)-p_{{\mathrm{f}}}^{*}\left(t+{\mathrm{d}}t\right) $

(6)

其中，

$ p_{{\mathrm{f}}}^{*} $

为无源性项时流体的预测压力，结合不可压缩流体的连续方程，质量源项满足

$ {S}_{\mathrm{m}}=\nabla \cdot \left(\rho \left({\boldsymbol{U}}_{{\mathrm{s}}}-\boldsymbol{U}_{{\mathrm{f}}}^{*}\right)\right),\left(x\in \mathit{\Omega }\right) $

(7)

S_m的物理意义是在固定网格下，为了使流体在界面处的压力连续，需要向流体单元注入的质量率，直接体现了式（3）的压力匹配约束，并嵌入流体连续方程中。

流体域控制方程采用惯性参考系下流场的连续性方程和动量方程，其表达式为

$ \frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \overrightarrow{\nu }\right)={S}_{{\mathrm{m}}} $

(8)

$ \frac{\partial }{\partial t}\left(\rho \overrightarrow{\nu }\right)+\nabla \cdot \left(\rho \overrightarrow{\nu }\overrightarrow{\nu }\right)=-\nabla p+\nabla \cdot \left(\overline{\overline{\tau }}\right)+\rho \overrightarrow{g}+{S}_{{\mathrm{M}}} $

(9)

其中，

$ {S}_{\mathrm{m}} $

和

$ {{S}}_{\mathrm{M}} $

分别为添加的质量源项和动量源项，仅在流固界面Ω上非零，直接将式（2）、（3）的匹配约束转化为流场的源项驱动，使固定网格下的流体场能够精准响应固体边界的微幅运动。

$ \overline{\overline{\tau }} $

为应力张量。

对于涡激振动问题，边界层流动一般是湍流流动。目前对于湍流的主要计算方法有直接数值模拟、雷诺应力平均及大涡模拟。本文采用

$ \text{SST}\;k\text{-}\omega $

湍流模型进行流场计算，

$ \text{SST}\;k\text{-}\omega $

模型的湍动能

$ k $

和比耗散率

$ \omega $

的输运方程如下

$ \frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial {x}_{i}}\left(\rho k{u}_{i}\right)=\frac{\partial }{\partial {x}_{j}}\left[{\mathit{\Gamma }}_{k}\frac{\partial k}{\partial {x}_{j}}\right]+{G}_{k}-{Y}_{k}+{S}_{k} $

(10)

$ \frac{\partial }{\partial t}\left(\rho \omega \right)+\frac{\partial }{\partial {x}_{i}}\left(\rho \omega {u}_{i}\right)=\frac{\partial }{\partial {x}_{j}}\left[{\mathit{\Gamma }}_{\omega }\frac{\partial \omega }{\partial {x}_{j}}\right]+{G}_{\omega }-{Y}_{\omega }+{S}_{\omega }$

(11)

式中：

$ {G}_{k} $

为由于平均速度梯度而产生的湍动能，

$ {G}_{\omega } $

为比耗散率产生的；

$ {\mathit{\Gamma }}_{k} $

为湍动能的有效扩散系数；

$ {\mathit{\Gamma }}_{\omega } $

为比耗散率的有效扩散系数；

$ {S}_{k} $

和

$ {S}_{\omega } $

为用户定义源项。

为了简化研究，侧重方法的建立，本文计算中结构为刚体，根据结构动力学相关知识，考虑圆柱单自由度振动，其控制方程如下

$ m\frac{{{\mathrm{d}}}^{2}y}{{\mathrm{d}}{t}^{2}}+c\frac{{\mathrm{d}}y}{{\mathrm{d}}t}+ky={F}_{{{}_{y}}}\left(t\right) $

(12)

其中，

$ m $

为圆柱的等效质量，

$ c $

为阻尼系数，

$ k $

为弹簧刚度系数，

$ y $

为圆柱在横流向的位移，F_y为圆柱所受升力，主要通过流场计算离散的N-S方程后使用用户自定义函数获得。

2　数值方法

收起

2.1　冻结边界法

冻结边界法的主要思想在于添加源项来替代网格运动效应，对于动网格方法，其边界移动的任意控制体积V上的一般标量

$ \phi $

的守恒方程的积分形式如下

$ \frac{{\mathrm{d}}}{{\mathrm{d}}t}{\int }_{V}\,\rho \phi {\mathrm{d}}V+{\int }_{\partial V}\,\rho \phi \left(\overrightarrow{u}-{\overrightarrow{u}}_{g}\right)\cdot {\mathrm{d}}\overrightarrow{A}={\int }_{\partial V}\,\mathit{\Gamma }\nabla \phi \cdot {\mathrm{d}}\overrightarrow{A}+{\int }_{V}\,{S}_{\phi }{\mathrm{d}}V $

(13)

其中，

$ \rho $

为流体密度，

$ \overrightarrow{u} $

为流速矢量，

$ {\overrightarrow{u}}_{g} $

为动网格的移动速度，

$ \mathit{\Gamma } $

为扩散系数，

$ {S}_{\phi } $

为

$ \phi $

的源项，

$ \partial V $

用于表示控制体积V的边界。

将代表网格运动的一项分离出来，有

$ \frac{{\mathrm{d}}}{{\mathrm{d}}t}{\int }_{V}\,\rho \phi {\mathrm{d}}V+{\int }_{\partial V}\,\rho \phi \overrightarrow{u}\cdot {\mathrm{d}}\overrightarrow{A}-{\int }_{\partial V}\,\rho \phi {\overrightarrow{u}}_{g}\cdot {\mathrm{d}}\overrightarrow{A}={\int }_{\partial V}\,\mathit{\Gamma }\nabla \phi \cdot {\mathrm{d}}\overrightarrow{A}+{\int }_{V}\,{S}_{\phi }{\mathrm{d}}V $

(14)

移项到方程的右边得

$ \frac{{\mathrm{d}}}{{\mathrm{d}}t}{\int }_{V}\,\rho \phi {\mathrm{d}}V+{\int }_{\partial V}\,\rho \phi \overrightarrow{u}\cdot {\mathrm{d}}\overrightarrow{A}={\int }_{\partial V}\,\mathit{\Gamma }\nabla \phi \cdot {\mathrm{d}}\overrightarrow{A}+{\int }_{V}\,{S}_{\phi }{\mathrm{d}}V+{\int }_{\partial V}\,\rho \phi {\overrightarrow{u}}_{g}\cdot {\mathrm{d}}\overrightarrow{A} $

(15)

通过高斯散度定理将网格运动项并入源项中，有

$ \frac{{\mathrm{d}}}{{\mathrm{d}}t}{\int }_{V}\,\rho \phi {\mathrm{d}}V+{\int }_{\partial V}\,\rho \overrightarrow{u}\cdot {\mathrm{d}}\overrightarrow{A}={\int }_{\partial V}\,\mathit{\Gamma }\cdot {\mathrm{d}}\overrightarrow{A}+{\int }_{V}\,{{{S}_{\phi }}}^{*}{\mathrm{d}}V $

(16)

其中，

$ {{{S}_{\phi }}}^{\mathrm{*}} $

为移项合并后得到的新的源项，具体为

$ {{{S}_{\phi }}}^{*}={S}_{\phi }+\nabla \cdot (\rho \phi {\overrightarrow{u}}_{g}) $

，因此可以通过添加源项的方式来表达物体边界的运动，从而在保持物体静止时实现其运动的模拟计算。

当物体在流场中运动时，因其表面不可穿透且考虑不可压缩流动，可以认为周围流体在物体的影响下也以相同的速度在运动。以圆心为坐标原点建立直角坐标系如图2，圆柱半径取r，则其表面任意一点的坐标表示为（rcosα, rsinα）。

以上表面一段微元ds为例进行分析，当圆柱以速度v向上运动时，对于其表面ds段的质量流量为

$ {\mathrm{d}}m=\rho \boldsymbol{v}\cdot {\mathrm{d}}\boldsymbol{s} $

(17)

其中，

$ \boldsymbol{v} $

和

$ \mathrm{d}\boldsymbol{s} $

为速度和面积微元矢量，在计算中以网格为最小单元，则ds为壁面网格面积，对该网格单元的质量源和动量源项有

$ {\mathrm{d}}{S}_{{\mathrm{m}}}=\frac{\rho \boldsymbol{v}\cdot {\mathrm{d}}\boldsymbol{s}}{{{\mathrm{Volume}}}_{{\mathrm{cell}}}}，{\mathrm{d}}{S}_{{\mathrm{M}}}=\frac{\rho \boldsymbol{v}\cdot {\mathrm{d}}\boldsymbol{s}}{{{\mathrm{Volume}}}_{{\mathrm{cell}}}} \boldsymbol{v} $

(18)

本文使用Fluent软件进行流动计算，通过Fluent用户自定义函数进行流固耦合计算。主要使用到的宏为DEFINE_EXECUTE_AT_END，该宏可以在稳态迭代步结束时调用，也可以在瞬态计算时间步迭代完毕后调用，此处主要用于设定在每个时间步获取物体受力和进行运动方程求解的部分，从而实现时间步内的信息交互和对物理场参数的周期性操作。在本文中，该宏被配置为在每个瞬态时间步收敛后执行以下关键任务：首先通过集成压力、粘性力的空间积分运算实时获取流场对固体边界的流体作用力与力矩，其次使用四阶龙格-库塔算法对结构运动方程进行数值求解，获得物体在该时间步内的位移、速度等运动学参数，最终将这些更新后的运动参数传递至流场求解模块，完成双向耦合系统的跨物理场数据同步。添加源项时使用宏DEFINE_SOURCE，此处主要用于在流固耦合界面添加质量和动量源项，实现了质量、动量守恒方程的源项修正。在质量守恒方程中，基于运动边界速度场分布建立连续性源项，补偿因虚拟边界位移造成的质量输运偏差。在动量守恒方程中，则将边界的加速度效应转化为等效动量源项，按照公式在圆柱周围添加对应的质量源和动量源，从而实现圆柱运动模拟。

2.2　计算对象

圆柱直径取0.033 m，为避免边界对计算结果造成影响，计算域的入口和上下边界都至少距离圆柱圆心10D，出口距离圆心30D。其计算域如图3所示，在网格无关性验证中去除阻尼和弹簧作用，定义计算域左下角为坐标原点，x方向为流向方向，横流向为y方向，流场介质为水。入口设置为速度入口，出口为压力出口，圆柱表面和计算域上下边界为无滑移壁面。

对于弹簧阻尼系统，设置圆柱质量m为20.73 kg，弹簧刚度k为715.5 N/m，阻尼为0.074 N·s/m，阻尼比

$ \xi $

为0.00231，圆频率

$ w $

为5.875，固有频率f_n为0.935。

2.3　计算方法

流动的数值求解使用商用软件Fluent，该软件使用有限体积法对控制方程进行离散。压力速度耦合求解方法选用SIMPLE算法，湍动能、比耗散率和动量均采用二阶迎风格式离散，流场中物理量梯度计算使用基于单元的最小二乘法。压力计算松弛因子取0.3，速度计算松弛因子取0.7。考虑四阶龙格−库塔法精度较高且算法步骤明了易于编程，所以使用该方法离散求解振动方程。流固耦合计算的流程图如图4所示。

2.4　网格无关性验证

对圆柱计算域网格采用全结构化网格划分，并在圆柱周围5D×5D的范围内设置网格加密区，其中采用O型网格划分以保证网格良好的正交性。为了精确捕捉圆柱表面涡的结构，对圆柱近壁面网格进行加密，取无量纲壁面距离y+为1，第一层网格高度为1.2E-4 m，向外网格尺寸逐渐增大，其增长率设置为1.1，加密区内网格最大尺寸为1.3 mm，加密区外最大网格尺寸为3 mm，网格总数为104084，加密区网格划分图如图5所示，采用瞬态计算，时间步长取0.005 s。

在雷诺数为3900的来流条件下，计算静态圆柱绕流，以进行网格无关性验证，利用阻力系数和斯特劳哈尔数（St）对计算结果进行量化。斯特劳哈尔数（St）、阻力系数

$ {C}_{\mathrm{d}} $

定义如下

$ {C}_{{\mathrm{d}}}=\frac{2{F}_{{\mathrm{d}}}}{\rho U_{\infty }^{2}A} $

(19)

$ St=\frac{{f}_{}D}{{U}_{\infty }} $

(20)

其中，

$ {F}_{\mathrm{d}} $

是顺流向的总力，f为涡脱落频率。

选取三组不同分辨率的网格，其数量分别为73500、104084和147200，按照

$ \sqrt{2} $

的比例依次增加，不同网格的结果如表1所示。可以看出，随着网格数增加，结果越来越接近，中等数量网格计算的结果与细网格较为接近，为减小计算量并保证计算精度，下文选择第二套网格进行计算。

2.5　圆柱单自由度强迫振动的冻结边界法验证计算

为了验证本文使用方法对物体运动计算的准确性，通过冻结边界法对二维圆柱的强迫振动进行计算，并与动网格方法计算结果进行对比。计算域设置和网格划分与上一部分相同，设置圆柱运动方程为y方向速度V=0.005sin（10t）。出口设置为压力出口，圆柱表面设置为壁面，动网格方法采用扩散光顺，湍流模型选择

$ \text{SST}\;k\text{-}\omega $

模型，设置时间步长为0.005 s，设置相同数量的线程进行两种方法计算速度的对比。

将计算稳定后圆柱y方向一周期受力进行对比，如图6所示，从图中可以看出圆柱受力均以正弦的形式发生变化，二者的周期均为0.63 s，符合对速度的设置，动网格方法计算圆柱受力最大值为0.0454 N，冻结边界法计算受力最大值为0.0438 N，冻结边界法相对动网格方法计算的偏差为3.5%，而使用相同的计算资源进行一周期的计算，动网格方法用时为13分钟，冻结边界法为9分钟，冻结边界法耗时仅为动网格方法的69%，可以看出冻结边界法在保证了计算精度的前提下很大程度上缩短了计算时间。

3　计算结果与分析

收起

3.1　圆柱单自由度涡激振动的冻结边界法计算

在入口设置不同来流速度，来流速度采用无量纲速度，公式为

$ {U}_{{\mathrm{r}}}=\frac{{U}_{\infty }}{{f}_{{\mathrm{n}}}D}$

(21)

其中，

$ {U}_{\mathrm{\infty }} $

为来流速度，

$ {f}_{\mathrm{n}} $

为系统固有频率，

$ D $

为圆柱直径，对约化速度U_r=1-10的工况进行计算分析。

在流场中，圆柱首先保持静止，当水流流经圆柱时，由于其表面的粘性作用以及逆压梯度的作用，圆柱表面发生边界层分离，从而产生周期性的横向力使圆柱振动。图7展示了一个典型速度下（

$ {U}_{\mathrm{r}} $

=7）圆柱表面升力系数和y方向模拟位移随时间变化的情况。在10 s前，圆柱表面的升力系数出现振荡，其幅度逐渐增大；10 s后，升力系数和圆柱模拟位移产生稳定的周期性振荡，此时，漩涡从圆柱表面以一定规律脱落。为计算涡脱落频率，选取升力系数稳定部分进行快速傅里叶变换，其分析结果如图8所示，在1.1 Hz处存在峰值，所以此时的涡脱落频率为1.1 Hz。

同理，改变来流速度对圆柱涡激振动进行计算，取圆柱稳定后y方向模拟位移最大值无量纲化，其计算公式如下

$ {A}^{*}=\frac{y}{D} $

(22)

式中：

$ y $

为圆柱运动的幅值，

$ D $

为圆柱直径。

分别对各工况进行傅里叶变换获取涡脱落频率，使用系统固有频率对涡脱落频率进行无量纲化。不同来流速度下的频率比与无量纲振幅如图9所示，随着来流速度的增大，涡脱落频率在约化速度3-5外呈现线性增大趋势，这符合斯特劳哈尔定律，其斜率即为St数，值为0.2，但振幅却表现为在约化速度小于3时增幅很小，约化速度为3时急剧增大。进一步分析频率图可知，在约化速度3-5时频率比为1，此时的涡脱落频率和涡激振动系统的固有频率相等，涡脱落频率被“锁定”在结构固有频率上，对应这段速度区间的y方向位移出现了显著的增大，并在约化速度为5时达到峰值，表明在这段速度区间内发生了共振，即涡脱落频率等于系统固有频率的现象。之所以在约化速度为5时达到最大振幅，是因为在这个来流速度下，圆柱自然涡脱落频率等于系统固有频率。当速度大于5时，涡脱落频率不再锁定在结构固有频率上，而是继续随着流速的增大而增大，从振幅图中也可以看出，此时的振幅快速下降。通过冻结边界法的计算，得到了圆柱涡激振动锁频的过程。

3.2　不同来流速度下尾涡发展对比

从上文的计算中可以看出，在锁频区和非锁频区的圆柱振动显著不同，而其激励源是脱落的尾涡，因此为了具体观察尾涡的不同，我们取约化速度2和约化速度5的一周期尾涡发展流线图，如图10所示。在约化速度为2时，0.25T时刻上表面涡旋发展至圆心高度处，此时下表面涡旋也在发展，二者开始相互剪切挤压；在0.5T时，下表面涡旋继续发展变大而导致上表面涡旋脱落，这可以定义为自然脱落涡。在约化速度为5时，上下表面的涡旋均在对称涡旋尚未发展时就已经脱落，结合约化速度为5时的振幅分析，可以认为此时圆柱的振动促进了涡的脱落，从而定义此时的涡主要为下泄涡。因此，可以认为圆柱尾涡的涡系由脱落涡和下泄涡两部分组成。

为了观察圆柱尾流区涡量的分布，分别取约化速度为2、5、8三个工况一周期内圆柱尾涡云图如图11~13所示。可以看出，涡流出现在圆柱的上下表面，并以相反的旋转方向脱落。约化速度为2时，来流速度较小，结合上一部分可知此时圆柱振动幅度也较低，其尾涡脱落情况和静止时基本一致，整体来看此工况下涡量强度较小，分析其具体脱落过程可知，在T=0时，刚好有一负涡完整地从圆柱表面脱落，下一个负涡正在圆柱尾部发展并且正涡在脱落的过程中，随着边界层的分离以及剪切层的相互作用，两种涡从圆柱表面脱落进入尾流场形成排列在圆柱正后方（y方向）长度较长的涡街。约化速度为5时，涡量强度增大，T=0时有一正涡刚好脱落，此时下一个正在形成的正涡卷起的高度越过了圆柱中心，较上一情况更细长，并且也提前与正在发展的负涡剪切导致其涡上的“鱼尾”状结构更长，在涡发展脱落后逐渐排列至圆柱上下两侧呈现x方向长度较长的涡街。约化速度为8时可以看出，涡间距明显增大，T=0时有一正涡刚好脱落，而此时下一个正在形成的负涡已经离圆柱表面较远，在剪切层的相互作用下其“鱼尾”状结构最长，在发展脱落后主要形成在圆柱后方排列的涡街。从尾涡的排列方式来看，约化速度为2和8时，正负涡带都排列在圆柱正后方；而约化速度为5时，则排列在圆柱后方的上下两侧。结合上一部分的振幅可以看出，在约化速度为5时，圆柱振幅较大，尾涡是在振动的影响下脱落形成的；而在振幅较小时，尾涡由剪切层相互作用自然脱落，因此排列在圆柱正后方。

由于振幅较大的涡系与圆柱振幅较小时显著不同，且此时发生了频率锁定现象。分析认为，在物体振动的过程中，其表面会额外生成振动涡系，与自然脱落涡系相互竞争影响，共同组成尾涡。随着振动幅度的增大，振动涡系强度也增大，从而占据主导地位控制尾涡脱落频率，表征为涡脱落频率锁定在结构固有频率上。来流速度进一步增大后，自然脱落涡强度增大，占据主导地位控制尾涡，脱离锁频阶段。因此，尾涡锁频是两种涡系相互竞争的结果。

4　结　论

收起

本文通过冻结边界法表达边界的振动效应，对圆柱强迫振动、单自由度涡激振动进行了数值模拟，并对不同约化速度下的涡脱落频率、圆柱振幅和尾涡发展进行了分析，表明了冻结边界法具有应用于小振幅涡激振动问题（如流激噪声）的潜力，获得以下结论：

（1）在圆柱强迫振动的计算中，冻结边界法计算结果与动网格方法一致，并且其计算一周期运动的速度提升了30%，表明本文的计算方法在保证计算精度的前提下有效提高了计算速度。

（2）通过本文的计算方法对圆柱单自由度涡激振动进行计算，成功计算出了锁频现象，在锁频区间内涡脱落频率被锁定在系统固有频率上，振幅显著增大，锁频区间外涡脱落频率与来流速度呈线性关系。

（3）锁频区间内由于共振的影响，尾涡由前后排列改变为上下排列的正负涡带，尾流涡量场在横流向的波动增大。

（4）尾涡锁频的机制主要在于振动下泄涡与自然脱落涡相互竞争，振幅较小时振动涡系强度也较小，此时脱落涡占据主导，而随着来流速度增大，振动频率接近固有频率，导致振幅逐渐增大，当振动涡系强度达到一定值时，即主导尾涡，从而表现为频率锁定和振幅显著增大。

（5）通过理论推导与典型算例验证，本文提出的冻结边界法为小振幅、强耦合流固耦合问题提供了兼具高效性与高精度的求解框架。船舶螺旋桨在非均匀伴流场中产生的小振幅（微米至毫米级）、强耦合、高频次流固耦合响应，正是该方法极具潜力的核心应用方向。

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

刘晖, 刘蓉, 吉柏锋, 等. 腐蚀海洋立管的涡激振动响应分析[J]. 振动与冲击, 2024, 43(16): 33‒41.

Liu

, Liu

, Ji

B F

, et al. Analysis of vortex-induced vibration response of corroded marine risers[J]. Vibration and Shock, 2024, 43(16): 33‒41. (in Chinese)

李玲瑶, 余志武, 何旭辉, 等. 基于可靠度理论的桥梁涡激振动概率性评价[J]. 工程力学, 2015, 32(10): 154‒160+168.

L Y

, Yu

Z W

, He

X H

, et al. Probabilistic evaluation of vortex-induced vibration of bridges based on reliability theory[J]. Engineering Mechanics, 2015, 32(10): 154‒160+168. (in Chinese)

杨振宇, 姚慧岚, 张怀新. 螺旋桨流激噪声与泄涡频率关系分析[J]. 船舶与海洋工程, 2020, 36(4): 1‒7+43.

Yang

Z Y

, Yao

H L

, Zhang

H X

. Analysis of the relationship between propeller-induced noise and wake vortex frequency[J]. Shipbuilding and Ocean Engineering, 2020, 36(4): 1‒7+43. (in Chinese)

Kim

, Hur

, Lee

. Numerical and experimental analysis of a singing propeller having blunt trailing edges[J]. Journal of Ship Research, 2020, 64(3): 234‒249.

王家正, 孙洪源, 林海花, 等. 考虑质量比影响的双圆柱涡激振动发电装置数值模拟研究[J]. 振动与冲击, 2024, 43(14): 8‒17.

Wang

J Z

, Sun

H Y

, Lin

H H

, et al. Numerical simulation study of vortex-induced vibration energy generation device with dual cylinders considering mass ratio effects[J]. Vibration and Shock, 2024, 43(14): 8‒17. (in Chinese)

张桂勇, 王双强, 孙哲, 等. 船舶与海洋工程流固耦合数值方法研究进展[J]. 中国舰船研究, 2022, 17(5): 52‒73.

Zhang

G Y

, Wang

S Q

, Sun

, et al. Research progress on fluid-structure interaction numerical methods in ship and ocean engineering[J]. Chinese Journal of Ship Research, 2022, 17(5): 52‒73. (in Chinese)

Han

, Huang

Q G

, Pan

, et al. Optimal energy harvesting efficiency from vortex-induced vibration of a circular cylinder[J]. Ocean Engineering, 2023, 282: 114869.

Jia

L Y

, Sang

, Shi

, et al. Investigation on numerical simulation of VIV of deep-sea flexible risers[J]. Applied Sciences, 2023, 13(14): 8096.

Wen

H J

, Zhao

Y M

, Ren

. SPH modelling of the vortex-induced vibration of a near-wall cylinder[J]. China Ocean Engineering, 2023, 37(3): 355‒368.

Jin

G Q

, et al. Numerical analysis of vortex-induced vibration on a flexible cantilever riser for deep-sea mining system[J]. Marine Structures, 2023, 87: 103334.

Jin

G Q

, et al. Analysis of flow-induced vibration characteristics for flexible riser attached with non-rotating water-drop fairing[J]. Ocean Engineering, 2022, 262: 112166.

Sun

P N

, Le

Touze D

, Oger

, et al. An accurate FSI-SPH modeling of challenging fluid-structure interaction problems in two and three dimensions[J]. Ocean Engineering, 2021, 221: 108552.

徐嘉. 基于离散涡模型的涡激振动与系泊动力学耦合分析方法[D]. 北京: 中国舰船研究院, 2024.

. Coupled analysis method of vortex-induced vibration and mooring dynamics based on the discrete vortex model[D]. Beijing: China Shipbuilding Research Institute, 2024. (in Chinese)

Sun

X S

, Sakai

, Yamada

. Three-dimensional simulation of a solid–liquid flow by the DEM–SPH method[J]. Journal of Computational Physics, 2013, 248: 147‒176.

Oger

, Doring

, Alessandrini

, et al. Two-dimensional SPH simulations of wedge water entries[J]. Journal of Computational Physics, 2006, 213(2): 803‒822.

Zheng

, Duan

W Y

, Ma

Q W

. A new scheme for identifying free surface particles in improved SPH[J]. Science China Physics, Mechanics and Astronomy, 2012, 55: 1454‒1463.

黄灿. 光滑粒子动力学方法若干关键数值问题的研究[D]. 北京: 北京理工大学, 2016.

Huang

. Research on several key numerical issues of the smooth particle hydrodynamics method[D]. Beijing: Beijing Institute of Technology, 2016. (in Chinese)

Antuono

, Bouscasse

, Colagrossi

, et al. A measure of spatial disorder in particle methods[J]. Computer Physics Communications, 2014, 185(10): 2609‒2621.

Zhao

, Kaja

, Xiang

, et al. Vortex-induced vibration (VIV) of a circular cylinder in combined steady and oscillatory flow[J]. Ocean Engineering, 2013, 73: 83‒95.

陈东阳, 肖清, 顾超杰, 等. 柱体结构涡激振动数值计算[J]. 振动与冲击, 2020, 39(19): 7‒12+47.

Chen

D Y

, Xiao

, Gu

C J

, et al. Numerical computation of vortex-induced vibration of cylindrical structures[J]. Vibration and Shock, 2020, 39(19): 7‒12+47. (in Chinese)

Zhang

J K

, Huang

, Zhang

M J

, et al. Investigation of vortex-induced vibration characteristics around the NACA0009 hydrofoil with a focus on the lock-in phenomenon[J]. Ocean Engineering, 2024, 312: 119082.

郭涛, 张纹惠, 王文全, 等. 基于IBM法的低雷诺数下涡激振动高质量比效应的研究[J]. 工程力学, 2022, 39(3): 222‒232.

Guo

, Zhang

H W

, Wang

W Q

, et al. Effects of high mass and damping ratio on VIV of a circular cylinder with low Reynolds number based on IBM[J]. Engineering Mechanics, 2022, 39(3): 222‒232. (in Chinese)

Wang

, Cao

, Liang

, et al. Experimental investigation of the effect of propeller characteristic parameters on propeller singing[J]. Ocean Engineering, 2022, 256: 111538.

陈静涛. 圆柱绕流的二维数值模拟和尾迹分析[J]. 计算机辅助工程, 2013, 22(6): 1‒6+27.

Chen

J T

. Two-dimensional numerical simulation and wake analysis of flow around a cylinder[J]. Computer-Aided Engineering, 2013, 22(6): 1‒6+27. (in Chinese)

2026年第30卷第4期

PDF下载

140

引用本文

BibTeX

文章信息

doi: 10.3969/j.issn.1007-7294.2026.04.002

接收时间：2025-08-09
首发时间：2026-07-07
出版时间：2026-04-15

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2025-08-09

基金

作者信息

^1.中国船舶科学研究中心，江苏无锡 214082

通讯作者:

袁国辉（2000–），男，硕士研究生，通讯作者，E-mail：1470307651@qq.com

参考文献

分享链接

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

分享至

全文二维码

扫描看全文

引用本文

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