本发明公开了一种刚性圆柱体横流与顺流方向涡激振动耦合响应预测方法，主要包括以下步骤：1)建立横流方向与顺流方向互为耦合的振动控制方程；2)基于有限差分法对耦合的振动控制方程进行求解；3)基于分析数据，对实例进行计算分析。本发明针对刚性圆柱体横流与顺流方向耦合振动响应进行了研究。建立了完整的横流、顺流方向的耦合振动模型，用于分析同时考虑横流方向振动、顺流方向振动以及结构几何非线性的刚性圆柱体涡激振动响应预测问题。该模型以很好地模拟出结构的锁定以及位移突变现象。
刚性圆柱体横流与顺流方向涡激振动耦合响应预测方法
技术领域
[0001] 本发明涉及圆柱体涡激振动响应领域，具体地说，特别涉及一种刚性圆柱体横流与顺流方向涡激振动耦合响应预测方法。
背景技术
[0002] 圆柱体广泛存在于各种实际工程中，圆柱体在一定的来流中，会在结构两侧形成交替脱落的漩涡，周期性的漩涡脱落会在结构上产生周期性的拖曳力以及升力。若结构为弹性支撑，周期性的拖曳力和升力会引起结构在顺流以及横流方向产生周期性的振动，统称为“涡激振动”(Gao,Y,Fu,SX,Xiong,YM,Zhao,Y,and Liu,LM,Experimental study on response performance of vortex‑induced vibration of a flexible cylinder,Ships and Offshore Structures,2017,12(1):116‑134)。当漩涡脱落频率出现在结构固有频率附近时，便会发生锁定现象。锁定区域内，圆柱体会发生大幅的、危险的涡激振动响应，这种响应会对结构带来很大的疲劳损伤。
[0003] 现有的研究圆柱体涡激振动响应的方法主要可分为实验方法以及数值方法。与数值方法相比，实验方法得到的数据更为可靠、得到的现象更为直观，但实验方法通常研究成本较高。数值方法通常主要包括计算流体动力学(Computational Fluid Dynamics,CFD)方法(Gao,Y,Zong,Z,Zou,L,Takagi,S,Jiang,ZY,Numerical simulation of vortex‑induced vibration of a circular cylinder with different surface roughnesses,Marine Structures,2018,57:165‑179)以及经验模型方法。与经验模型方法相比，CFD方法的计算精度要高，但CFD方法的数值计算成本要比经验模型方法高很多。为了适应海洋工程中的多工况问题，有必要建立一种能够快速预报结构涡激振动响应主要特性的经验模型方法，而尾流振子模型方法则正是近些年广泛应用于涡激振动响应预报中的一种经验模型方法。
[0004] 目前基于尾流振子模型方法对刚性圆柱体涡激振动响应的研究主要集中只考虑结构在横流方向的振动，而忽略了结构在顺流方向的振动。而实际工程中横流方向与顺流方向振动互为耦合、互相影响。在少量的同时考虑横流以及顺流方向振动响应的研究中，针对结构几何非线性对涡激振动响应的影响研究则更少。对于这种同时考虑结构几何非线性以及CF方向与IL方向互为耦合的圆柱体模型，其涡激振动响应特性有待进一步展开研究。
发明内容
[0005] 为了解决现有技术的问题，本发明实施例提供了一种刚性圆柱体横流与顺流方向涡激振动耦合响应预测方法。所述技术方案如下：
[0006] 本发明提供了一种刚性圆柱体横流与顺流方向涡激振动耦合响应预测方法，主要包括以下步骤：
[0007] 1)建立横流方向与顺流方向互为耦合的振动控制方程；
[0008] 2)基于有限差分法对耦合的振动控制方程进行求解；
[0009] 3)基于分析数据，对实例进行计算分析。
[0010] 进一步的，所述步骤1)具体包括：
[0011] a、建立均匀来流下刚性圆柱体二自由度涡激振动耦合模型；
[0012] b、在所述均匀来流下刚性圆柱体二自由度涡激振动耦合模型中，考虑一单位长度、直径为D的刚性圆柱体在均匀来流U作用下引起的CF方向以及IL方向互为耦合的涡激振动响应问题，其中X方向为IL方向，Y方向为CF方向。
[0013] 其中，在分析过程中假设系统在X方向和Y方向呈对称分布，因此在X方向的质量、刚度、阻尼以及固有频率均与Y方向对应的参数相等；
[0014] 所述均匀来流下刚性圆柱体二自由度涡激振动耦合模型中设有4弹簧刚性圆柱体系统。
[0015] 进一步的，所述4弹簧刚性圆柱体系统在X和Y方向的振动方程描述如下：
[0016]
[0017] 式(1)d表示对变量求全导；表示对变量求偏导；T为时间； 为拉格朗日函数，表示为： 其中K为动能，P为势能；u1以及u2分别为X方向和Y方向位移；和 分别为X方向和Y方向速度；F1和F2分别表示作用在结构上的X方向以及Y方向水动力载荷Fx以及Fy；r为阻尼系数，包括结构阻尼系数rs以及流体阻尼系数rf两项，流体阻尼系数rf表示为：rf＝2 2
ωfγρD＝(2πStU/D)γρD ，其中ωf为漩涡脱落频率；St为斯脱哈尔数，这里取为0.2；ρ为流体密度；γ为粘滞力系数，取为0.8。
[0018] 进一步的，假设与圆柱体相连的4根弹簧初始长度均为S，弹簧弹性系数均为k，并假设圆柱体T时刻在X方向以及Y方向的结构位移分别为X(T)以及Y(T)，便得到此时刻圆柱体的动能K以及势能P，表示如下：
[0019]
[0020]
[0021] 式(2)中m为单位长度系统总质量，由结构质量ms以及流体附加质量mf两部分构成，2
流体附加质量表示为：mf＝CMρDπ/4，其中CM为附加质量系数，对于圆柱体取1.0；将式(2)和式(3)代入式(1)中，得到圆柱体在X方向和Y方向的振动方程为:
[0022]
[0023]
[0024] 方程式(4)和(5)中的非线性项体现了结构振动的几何非线性特征。为了求解式(4)和式(5)，需要对二式进行一定的简化处理；函数f(X,Y)在(X0,Y0)的二元泰勒级数展开式表示如下：
[0025]
[0026] 将式(4)和式(5)中的非线性项按照式(6)在(0,0)处进行泰勒级数展开，取前三阶得到：
[0027]
[0028]
[0029] 将式(7)代入式(4)、式(8)代入(5)中得到：
[0030]
[0031] 式(9)中h＝2k， 为了方便表达，记： 式(9)进一步写作：
[0032]
[0033] 式(10)中，h为系统刚度，已知系统总质量m以及系统刚度h，便可得到系统在X方向3 3
以及Y方向的固有频率，表示为： 式(10)中，X以及Y前面的系数 以
2 2
及 反应了结构在轴向变形的几何非线性特性；XY以及YX 前面的系数 以及 反应了CF以及IL方向振动互为耦合的的几何非线性特性；本发明将这 以及 这4个
参数统称为结构的几何非线性系数。
[0034] 进一步的，对于静止圆柱体，作用在圆柱体上的Fx以及Fy的方向分别与拖曳力FD以及升力FL方向一致；但对于振荡圆柱体，由于结构的振动，导致结构与流体之间的相对速度方向不再沿X方向，使得：Fx将不再与FD方向一致，Fy也将不再与FL方向一致。
[0035] 进一步的，假设结构的振动速度V方向与流速U方向的夹角为θ，通常情况下振动速度V与流速U比起来可以看作是小量，因此夹角θ也为小量；此时Fx以及Fy表示如下：
[0036]
[0037] 式(11)中 为平均拖曳力，FD为振荡拖曳力，FL为升力，分别表示如下：
[0038]
[0039] 式(12)中， CD以及CL分别表示平均拖曳力系数、振荡拖曳力系数以及升力系数，这里平均拖曳力系数 取为1.2，CD和CL分别写作：
[0040]
[0041] 式(13)中，CD0与CL0分别表示静止圆柱体上的振荡拖曳力系数以及升力系数，分别取为0.2以及0.3；p(T)以及q(T)分别表示在X以及Y方向的尾流振子的运动；
[0042] 采用Van der Pol方程来描述尾流振子特性，表示如下：
[0043]
[0044] 式(14)中εx、εy、Ax以及Ay分别为经验参数，Ax和Ay取为12，εx取为0.3，而εy则参照数*值与实验的对比结果，对预设数据进行拟合，得到εy关于m的表示式如下：
[0045]
[0046] 将方程(10)和(14)转换为无量纲形式，引入IL以及CF方向的无量纲位移x、y以及无量纲时间t，如下：
[0047] x＝X/D,y＝Y/D,t＝T·ωn (16)
[0048] 将式(16)带入方程式(10)和(14)中，整理得到X方向以及Y方向结构和尾流振子的无量纲方程为：
[0049]
[0050] 式(17)中，αx、βx、αy以及βy为结构的几何非线性无量纲系数，ξ为结构阻尼比，与结构阻尼系数rs之间的关系为：ξ＝rs/2mωn；式(17)中频率比Ω、质量比μ、以及系统无量纲参数 MD以及ML分别表示如下：
[0051]
[0052] 式(17)以及式(18)中Ur为折合速度：Ur＝2πU/ωnD，在分析中使用的质量比为m*，*即：结构质量与被结构排开的液体质量的比值，其与μ之间的关系为：m＝4μ/π‑CM。
[0053] 进一步的，所述步骤2)具体如下：
[0054] 在时间上采用二阶精度中心差分格式对方程式(17)进行先离散后迭代求解；将计算无量纲总时间ttotal划分为N段，得到每个时间步长Δt＝ttotal/N，被离散后的N+1个时刻依次记为：t＝tj(j＝0,1,2,..,N)；假设tn时刻对应的无量纲参数x、y、p以及q分别表示为xn、yn、pn以及qn；那么tn时刻式(17)中各导数项的二阶精度中心差分格式表示为：
[0055]
[0056] 将式(19)代入式(17)整理得到：
[0057]
[0058] x和y的初始条件(t＝t0)设为位移以及速度均为0，即：x0＝y0＝0以及dx/dt＝dy/dt＝0；p和q的初始条件(t＝t0)为：p0以及q0均取一微小的振幅以及dp/dt＝dq/dt＝0；由dx/dt＝dy/dt＝dp/dt＝dq/dt＝0得到：
[0059] x‑1＝x1,y‑1＝y1,p‑1＝p1,q‑1＝q1 (21)
[0060] 将式(21)代入式(20)得到t1时刻x、y、p以及q的值，表示如下：
[0061]
[0062] 依据式(22)得到y1、q1、x1以及p1；至此得到了t0以及t1时刻的x、y、p以及q的值；在此基础上，依据式(20)通过迭代依次求解得到t2,t3,…,tN时刻的x、y、p以及q的值。
[0063] 进一步的，所述步骤3)具体如下：
[0064] 对尾流振子数值模型进行数值验证；在验证过程中，选取的计算参数m*＝2.36,ζ＝0.006，计算总时间ttotal取为1000个无量纲时间，时间步长Δt取为0.1，折合速度Ur区间取为1.0‑14.0，区间间隔ΔUr取为0.1。
[0065] 本发明实施例提供的技术方案带来的有益效果是：
[0066] 本发明针对刚性圆柱体横流与顺流方向耦合振动响应进行了研究。建立了完整的横流、顺流方向的耦合振动模型，用于分析同时考虑横流方向振动、顺流方向振动以及结构几何非线性的刚性圆柱体涡激振动响应预测问题。该模型可以很好地模拟出结构的锁定以及位移突变现象。
附图说明
[0067] 为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0068] 图1是本发明实施例的均匀来流下刚性圆柱体二自由度涡激振动耦合模型的示意图；
[0069] 图2是本发明实施例的振动圆柱体受力示意图；
[0070] 图3是本发明实施例的数值模型CF与IL方向振动幅值与Stappenbelt等人的实验结果对比示意图。
具体实施方式
[0071] 为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。
[0072] 本发明提供了一种刚性圆柱体横流与顺流方向涡激振动耦合响应预测方法，主要包括以下步骤：
[0073] 1)建立横流方向与顺流方向互为耦合的振动控制方程
[0074] 2)基于有限差分法对耦合的振动控制方程进行求解；
[0075] 3)基于分析数据，对实例进行计算分析。
[0076] 步骤1：建立结构和流场互为耦合的振动控制方程，具体如下：
[0077] 如图1所示，考虑一单位长度、直径为D的刚性圆柱体在均匀来流U作用下引起的CF方向以及IL方向互为耦合的涡激振动响应问题，其中X方向为IL方向，Y方向为CF方向。本发明在分析过程中假设系统在X方向和Y方向呈对称分布，因此在X方向的质量、刚度、阻尼以及固有频率均与Y方向对应的参数相等。
[0078] 依据Hamilton原理，图1所示的4弹簧刚性圆柱体系统在X和Y方向的振动方程可描述如下：
[0079]
[0080] 式(1)中T为时间；为拉格朗日函数，可表示为： 其中K为动能，P为势能；u1以及u2分别为X方向和Y方向位移；和 分别为X方向和Y方向速度；F1和F2分别表示作用在结构上的X方向以及Y方向水动力载荷Fx以及Fy；r为阻尼系数，包括结构阻尼系数rs以及
2 2
流体阻尼系数rf两项，流体阻尼系数rf可表示为：rf＝ωfγρD ＝(2πStU/D)γρD ，其中ωf为漩涡脱落频率；St为斯脱哈尔数，这里取为0.2；ρ为流体密度；γ为粘滞力系数，可取为
0.8(Facchinetti,ML,de Langre,E,Biolley,F,Coupling of structure and wake oscillators in vortex‑induced vibrations,Journal of Fluids and Structures,
2004,19:123‑140)。如图1所示，假设与圆柱体相连的4根弹簧初始长度均为S，弹簧弹性系数均为k，并假设圆柱体T时刻在X方向以及Y方向的结构位移分别为X(T)以及Y(T)，便可得到此时刻圆柱体的动能K以及势能P，表示如下：
[0081]
[0082]
[0083] 式(2)中m为单位长度系统总质量，由结构质量ms以及流体附加质量mf两部分构成，2
流体附加质量可表示为：mf＝CMρDπ/4，其中CM为附加质量系数，对于圆柱体通常取1.0(Govardhan,R,Williamson,CHK,Modes of vortex formation and frequency response of a freely vibrating cylinder,Journal of Fluid Mechanics,2000,420:85‑130)。将式(2)和式(3)代入式(1)中，可以得到圆柱体在X方向和Y方向的振动方程为:
[0084]
[0085]
[0086] 方程式(4)和(5)中的非线性项体现了结构振动的几何非线性特征。为了求解式(4)和式(5)，需要对二式进行一定的简化处理。函数f(X,Y)在(X0,Y0)的二元泰勒级数展开式可表示如下：
[0087]
[0088] 将式(4)和式(5)中的非线性项按照式(6)在(0,0)处进行泰勒级数展开，取前三阶得到：
[0089]
[0090]
[0091] 将式(7)代入式(4)、式(8)代入(5)中得到：
[0092]
[0093] 式(9)中h＝2k， 为了方便表达，记： 式(9)可进一步写作：
[0094]
[0095] 式(10)中，h为系统刚度，已知系统总质量m以及系统刚度h，便可得到系统在X方向3 3
以及Y方向的固有频率，可表示为： 式(10)中，X以及Y前面的系数
2 2
以及 反应了结构在轴向变形的几何非线性特性；XY以及YX前面的系数 以及 反应了CF以及IL方向振动互为耦合的的几何非线性特性；本发明将这 以及 这4个
参数统称为结构的几何非线性系数。
[0096] 对于静止圆柱体，作用在圆柱体上的Fx以及Fy的方向分别与拖曳力FD以及升力FL方向一致。但对于振荡圆柱体，由于结构的振动，导致结构与流体之间的相对速度方向不再沿X方向，使得：Fx将不再与FD方向一致，Fy也将不再与FL方向一致。
[0097] 如图2所示，假设结构的振动速度V方向与流速U方向的夹角为θ，通常情况下振动速度V与流速U比起来可以看作是小量，因此夹角θ也为小量。此时Fx以及Fy可表示如下(Wang,XQ,So,RMC,Chan,KT,A non‑linear fluid force model for vortex‑induced vibration of an elastic cylinder,Journal of Sound and Vibration,2003,260:287‑305)：
[0098]
[0099] 式(11)中 为平均拖曳力，FD为振荡拖曳力，FL为升力，分别表示如下：
[0100]
[0101] 式(12)中， CD以及CL分别表示平均拖曳力系数、振荡拖曳力系数以及升力系数，这里平均拖曳力系数 取为1.2(Ge,F,Long,X,Wang,L,Hong,Y,Flow‑induced vibration of a long circular cylinders modelled by coupled nonlinear oscillators,Science in China Series G:Physics,Mechanics and Astronomy,2009,52(7):1086‑1093)，CD和CL分别写作(Facchinetti,ML,de Langre,E,Biolley,F,Coupling of structure and wake oscillators in vortex‑induced vibrations,Journal of Fluids and Structures,2004,19:123‑140)：
[0102]
[0103] 式(13)中，CD0与CL0分别表示静止圆柱体上的振荡拖曳力系数以及升力系数，分别取为0.2以及0.3(Srinil,N.,Zanganeh,H.Modelling of coupled cross‑flow/in‑line vortex‑induced vibrations using double Duffing and van der Pol oscillators.Ocean Engineering,2012；53:83‑97)。p(T)以及q(T)分别表示在X以及Y方向的尾流振子的运动。
[0104] 采用Van der Pol方程来描述尾流振子特性，可表示如下(Nayfeh,N.H.Introduction to Perturbation Techniques.New York:Wiley,1993)：
[0105]
[0106] 式(14)中εx、εy、Ax以及Ay分别为经验参数，Ax和Ay可取为12，εx可取为0.3(Facchinetti,ML,de Langre,E,Biolley,F,Coupling of structure and wake oscillators in vortex‑induced vibrations,Journal of Fluids and Structures,2004,19:123‑140)，而εy则参考文献中(Stappenbelt,B.,Laiji,F.,Tan,G.Low mass th
ratio vortex‑induced motion.In:the 16 Australian Fluid Mechanics Conference,Gold,Coast,Australia,2007,pp.1491‑1497)表1(b)中数值与实验的对比结果，对表中数*
据进行拟合，可得到εy关于m的表示式如下：
[0107]
[0108] 将方程(10)和(14)转换为无量纲形式，引入IL以及CF方向的无量纲位移x、y以及无量纲时间t，如下：
[0109] x＝X/D,y＝Y/D,t＝T·ωn (16)
[0110] 将式(16)带入方程式(10)和(14)中，整理得到X方向以及Y方向结构和尾流振子的无量纲方程为：
[0111]
[0112] 式(17)中，αx、βx、αy以及βy为结构的几何非线性无量纲系数，ξ为结构阻尼比，与结构阻尼系数rs之间的关系为：ξ＝rs/2mωn。式(17)中频率比Ω、质量比μ、以及系统无量纲参数 MD以及ML分别表示如下：
[0113]
[0114] 式(17)以及式(18)中Ur为折合速度：Ur＝2πU/ωnD，这里值得注意的是通常在分*析中大多时使用的质量比为m ，即：结构质量与被结构排开的液体质量的比值，其与μ之间*
的关系为：m＝4μ/π‑CM。
[0115] 步骤2：基于有限差分法对耦合的振动控制方程进行求解，具体如下：
[0116] 在时间上采用二阶精度中心差分格式对方程式(17)进行先离散后迭代求解。将计算无量纲总时间ttotal划分为N段，得到每个时间步长Δt＝ttotal/N，被离散后的N+1个时刻依次记为：t＝tj(j＝0,1,2,..,N)。假设tn时刻对应的无量纲参数x、y、p以及q分别表示为xn、yn、pn以及qn。那么tn时刻式(17)中各导数项的二阶精度中心差分格式可表示为：
[0117]
[0118] 将式(19)代入式(17)整理得到：
[0119]
[0120] x和y的初始条件(t＝t0)设为位移以及速度均为0，即：x0＝y0＝0以及dx/dt＝dy/dt＝0；p和q的初始条件(t＝t0)为：p0以及q0均取一微小的振幅以及dp/dt＝dq/dt＝0。由dx/dt＝dy/dt＝dp/dt＝dq/dt＝0得到：
[0121] x‑1＝x1,y‑1＝y1,p‑1＝p1,q‑1＝q1 (21)
[0122] 将式(21)代入式(20)可得到t1时刻x、y、p以及q的值，表示如下：
[0123]
[0124] 依据式(22)可以此得到y1、q1、x1以及p1。至此得到了t0以及t1时刻的x、y、p以及q的值。在此基础上，依据式(20)可通过迭代依次求解得到t2,t3,…,tN时刻的x、y、p以及q的值。
[0125] 步骤3：基于分析数据，对实例进行计算分析。
[0126] 为了验证本发明中尾流振子数值模型的可靠性，这里对其进行了数值验证。在验证过程中，选取的计算参数与Stappenbelt等(2007)的实验参数相同(m*＝2.36,ξ＝0.006)(Stappenbelt,B.,Laiji,F.,Tan,G.Low mass ratio vortex‑induced motion.In:the th16 Australasian Fluid Mechanics Conference,Gold,Coast,Australia,2007,pp.1491‑1497)，几何非线性无量纲系数αx、βx、αy以及βy均取为0.7(Srinil,N.,Zanganeh,H.Modelling of coupled cross‑flow/in‑line vortex‑induced vibrations using double Duffing and van der Pol oscillators.Ocean Engineering,2012；53:83‑97)，计算总时间ttotal取为1000个无量纲时间，时间步长Δt取为0.1，折合速度Ur区间取为1.0‑
14.0，区间间隔ΔUr取为0.1。
[0127] 图3给出了Stappenbelt等人实验方法以及本发明中尾流振子模型方法得到的CF以及IL方向的振动幅值。这里值得注意的是：由于本发明中数值模型得到的Ax/D包括两部分，即平均拖曳力 引起的平均位移以及振荡拖曳力FD引起的振荡位移。因此，在后续分析中需要将 引起的平均位移加以消除，只留下FD引起的振荡位移。由图3可以看出：本发明中数值方法以及他人实验方法得到的Ay/D以及Ax/D随折合速度Ur的变化趋势一致：随着Ur的增加，振动幅值明显增大，响应区间从初始分支(initial branch)进入上分支(upper branch)；随后Ur的进一步增加，振动幅值出现突降，进入低分支(lower branch)。实验方法以及本发明中数值方法得到的不同折合速度下Ay/D以及Ax/D的最大值均非常接近：Ay/D最大值均在1.5附近；Ax/D最大值均在0.25附近。且两种方法得到的振动幅值最大值均出现在Ur＝8.5附近。总体上来说，本发明中尾流振子模型方法的数值计算结果与Stappenbelt等人的实验结果吻合较好，从而验证了本发明中数值模型的可靠性。
[0128] 本发明实施例提供的技术方案带来的有益效果是：
[0129] 本发明针对刚性圆柱体横流与顺流方向耦合振动响应进行了研究。建立了完整的横流、顺流方向的耦合振动模型，用于分析同时考虑横流方向振动、顺流方向振动以及结构几何非线性的刚性圆柱体涡激振动响应预测问题。该模型可以很好地模拟出结构的锁定以及位移突变现象。
[0130] 以上仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
