目前研究浮式风机动力响应的方法主要有模型试验和数值模拟两种，模型试验又可以细分为在海岸上进行的较大尺寸模型试验^[1]及在波浪水池中进行的较小尺寸模型试验^[2-4]2种类型。文献[5-8]主要基于势流理论研究了浮式风机的动力响应，并利用Morison方程^[9]或者试验数据来考虑粘性影响，其中Morison方程是一种半经验公式，难于准确反映粘性影响；并且试验数据十分稀少，存在尺度效应，无法真实地反映粘性影响。有学者借助计算流体力学方法基于粘性流理论对风机动力响应问题进行数值模拟^[10-13]，但是这种方法也存在耗时巨大、效率低下、难于模拟随机风浪的缺点。

本文主要基于势流理论来构造浮式风机的时域运动方程。在此基础上，采用计算流体力学(computational fluid dynamics，CFD)方法考虑粘性影响，采用准静定系泊模型来计算系泊力，采用叶素动量理论来计算风力，对DeepCwind半潜式平台风机进行了动力响应分析，并将所得结果与试验结果以及文献[5, 7, 10]的数值模拟结果进行了对比分析。

1 浮式风机动力响应分析基本原理 1.1 时域运动方程及粘性阻尼计算原理概述

浮式风机时域运动方程为

$ \begin{array}{l} \sum\limits_{J = 1}^6 {\left[ {\left( {{M_{ij}} + {\mu _{ij}}} \right){{\ddot X}_j}\left( t \right) + \int_0^t {{K_{ij}}\left( {t - \tau } \right){{\dot X}_j}\left( \tau \right){\rm{d}}\tau {\rm{ + }}} } \right.} \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\left. {{b_{ij}}\left| {{{\dot X}_j}\left( t \right)} \right|{{\dot X}_j}\left( t \right) + {C_{ij}}{X_j}} \right] = \\ \;\;\;\;\;\;\;\;\;{F_{{\rm{W}}\mathit{i}}}\left( t \right) + {F_{{\rm{M}}\mathit{i}}}\left( t \right),\;\;\;\;i = 1,2, \cdot \cdot \cdot ,6 \end{array} $

(1)

式中:μ_ij为频率无穷大时的附加质量；K_ij为时延函数；F_Wi为时域广义波浪力；b_ij、F_Mi分别为非线性阻尼系数以及广义系泊力。

对于半潜式平台浮式风机而言，主要由界层分离所产生的粘性阻尼在总阻尼中占很大成分^[5]，本文采用CFD方法来考虑这一主要成分为非线性的粘性阻尼，且仅计及非线性阻尼系数矩阵的主对角线成分，忽略耦合项影响。

1.2 系泊力计算原理 1.2.1 本文所采用系泊模型

本文所采用不计系泊动力效应、忽略系缆弯矩和扭矩的准静定模型^[14]。考虑如图 1所示的单根系缆，图中H_A、V_A和H、V分别表示系缆在锚处以及导缆孔处所受张力的水平、垂向分量，l和h则表示系缆的水平、垂向投影长度，可表示为:

$ \left\{ \begin{array}{l} l = \frac{H}{W}\left[ {{{\sinh }^{ - 1}}\left( {\frac{V}{H}} \right)} \right] - {\sinh ^{ - 1}}\left( {\frac{{V - WL}}{H}} \right) + \frac{{HL}}{E}\\ h = \frac{H}{W}\left[ {\sqrt {1 + {{\left( {\frac{V}{H}} \right)}^2}} - \sqrt {1 + {{\left( {\frac{{V - WL}}{H}} \right)}^2}} } \right] + \\ \;\;\;\;\;\;\frac{1}{E}\left( {VL - \frac{{W{L^2}}}{2}} \right) \end{array} \right. $

(2)

式中:L为系缆未拉伸时的原长；W为海水中单位长度系缆的重量；E系缆的轴向刚度。

当V>WL时系缆不会产生躺底段，式(2)成立；而当V≤WL时，系缆将会产生躺底段(如图 1(a)中实线所示，图中L_B即为躺底段长度，且L_B=L-V/W)，式(2)转化为

$ \left\{ \begin{array}{l} l = {L_B} + \frac{H}{W}{\sinh ^{ - 1}}\left( {\frac{V}{H}} \right) + \frac{{HL}}{E} + \\ \;\;\;\;\;\frac{{{C_B}W}}{{2E}}\left[ {\mu \left( {L - \frac{V}{W} - \frac{H}{{{C_B}W}}} \right) - {{\left( {L - \frac{V}{W}} \right)}^2}} \right]\\ h = \frac{H}{W}\left[ {\sqrt {1 + {{\left( {\frac{V}{H}} \right)}^2}} - 1} \right] + \frac{{{V^2}}}{{2EW}} \end{array} \right. $

(3)

式中C_B为海底摩擦力系数。

1.2.2 StarCCM+系泊模型

StarCCM+的Catenary Coupling计算模块所采用的系泊模型为经典的弹性系缆在两固定点之间悬垂的模型，当V≤WL时，系缆会如图 1(b)所示呈凹形分布而不会产生躺底段，图中u₁ (x₁ , y₁)、u₂ (x₂ , y₂)分别表示参数u在锚处和导缆孔处的值，则在图中系缆上任取一点(x, y)，可得如下参数方程:

$ \left\{ \begin{array}{l} x = \frac{H}{W}u + \frac{{{V^2}}}{{EW}}\sinh \left( u \right) + \alpha \\ y = \frac{H}{W}\cos \left( u \right) + \frac{{{V^2}}}{{2EW}}{\sinh ^2}\left( u \right) + \beta \end{array} \right. $

(4)

式中:参数u满足sinh(u)=dy/dx；α、β为常量。

需要指出的是，在不考虑躺底段的情形下，该系泊模型与本文所采用的系泊模型是等价的。

1.3 风力计算原理

定常风风速的空间分布按下式进行计算:

$ V\left( Z \right) = {V_{{\rm{hub}}}}{\left( {Z/{Z_{{\rm{hub}}}}} \right)^\alpha } $

(5)

式中:Z_hub为风机轮毂距离水面高度；V_hub为风速；a为幂律指数。

随机风风速频域分布遵循NPD风谱^[5]，利用逆离散傅里叶变换方法^[15]，由NPD风谱得到风机轮毂处的风速时历，进而运用式(5)得到风速空间分布。本文主要采用增加了Prandtl叶尖损失修正^[15]以及Grauert修正^[15]的叶素动量理论^[15]来计算风机叶片的风力。

对于塔柱所受风力，本文采用FAST软件中使用的方法进行简化求解。将塔柱划分为若干面元，对面元积分以计算塔柱风力:

$ F=\iint\limits_{S}{0.5{{\rho }_{\text{air}}}{{C}_{\text{d}}}V_{\text{n}}^{2}\text{d}\mathit{s}} $

(6)

式中:ρ_air为空气密度；C_d为阻力系数；V_n为面元法向相对风速。

2 浮式风机动力响应数值计算 2.1 水动力系数及时延函数

本文的研究对象为DeepCwind半潜式平台浮式风机，图 2以侧视图与俯视图的方式展示了其主要的结构尺寸、系泊系统布置情况、以及空间固定坐标系位置、来波来风方向等研究条件和参数。关于该浮式风机系统更详实的数据可参阅文献[3-5]。

水动力系数以及时延函数的结果对比如图 3所示，图中A_ij、B_ij分别表示频域附加质量以及阻尼系数，K_ij则表示时延函数，本文的计算结果与文献[5, 16]的结果吻合得很好。

2.2 系泊系统特性

结合系泊系统的布置情况(图 2)以及系缆张力的计算原理，可以推断出，当风机不发生较大位移的运动时，系缆均会有躺底段存在，因此本文假设存在躺底段，以式(3)来求解各系缆张力，再根据所得结果验证假设是否成立，若假设不成立，则再利用不考虑躺底段的式(2)来计算系缆张力。

图 4为纵荡与横荡2个模态下系泊回复力随风机位移变化的曲线，其中StarCCM+曲线的系泊回复力根据式(2)计算得出的。可以看出，StarCCM+曲线较试验值有所偏小，而本文所采用的系泊模型所得的结果则与试验值吻合的很好。

2.3 粘性阻尼 2.3.1 边界条件设置及网格划分

本文CFD数值模拟采用的是SST k-ω湍流模型，对于自由液面，采用液体体积方法(fluid of volume, VOF)来进行捕捉，对于风机运动，采用重叠网格法来进行处理。图 5(a)显示了背景区域(400 m×400 m×300 m)的边界条件，其中，顶部、底部及左侧为速度入口，右侧为压力出口，水深(即静水面到底部的垂向距离)为200 m，而对于重叠区域(140 m×140 m×90 m)，模型表面为物面条件，重叠区域与背景区域的交界面则为overset边界条件；图 5(b)是网格剖面图，显示了网格划分情况，对于自由表面区域以及重叠区域外围的背景区域，本文均做了网格加密，以更好的捕捉自由液面形状和风机运动。

在对浮式风机进行自由衰减数值模拟时，本文考虑了纵荡、垂荡以及纵摇3种模态，并均计及了系泊系统。对于垂荡、纵摇运动而言，回复力主要由平台自身产生，系泊回复力对运动的影响较小，因此本文选择采用Catenary Coupling模块来模拟系泊力。而对于纵荡运动而言，其回复力主要由系泊系统来提供，而StarCCM+软件的Catenary Coupling模块所得的纵荡系泊回复力相对于试验值偏小，无法准确模拟风机的纵荡自由衰减运动，故本文选择利用Linear Spring Coupling模块来模拟系泊力，此模块中的弹簧刚度系数由本文所采用的系泊模型计算得出。需要指出的是，根据文献[5]，为了考虑模型试验时的线缆对平台纵荡运动的影响，在进行纵荡模拟时，本文附加了一个刚度系数为7.39 kN/m的约束弹簧。

2.3.2 网格数量及时间步长选取

如图 6所示，以纵摇为例，在不同的网格数量以及不同的时间步长情况下，纵荡衰减曲线的模拟结果整体差距很小，说明数值模拟结果稳定可靠；同时，注意到中等网格的结果与精细网格的结果差距甚微，而模拟耗时则少得多，因此本文采用中等网格方案，同理，本文也采用中等时间步长方案。

3种模态的网格数量及时间步长选取结果如表 1所示。

2.3.3 阻尼系数计算结果

表 2为非线性阻尼系数的计算结果对比，图 7为阻尼比的计算结果对比。

如图 7所示，对于纵摇模态，FAST以及CFD的阻尼比计算结果均与试验值差别不大，而且可以看到CFD的结果相较于FAST的结果更接近于试验值；对于垂荡模态，从整体上而言，本文的CFD结果偏小于试验值，而文献[10]的CFD结果则是偏大于试验值，对于较小的初始垂荡位移，FAST的计算结果更接近于试验值，而对于较大的初始垂荡位移，则是CFD的模拟结果更接近于试验值；试验值与数值模拟结果差别最大的是纵荡模态，在大初始纵荡位移情况下FAST以及CFD结果均明显偏小于试验值，原因可能在于本文及文献[10]都是对风机进行实尺度数值模拟，而试验值则是由模型尺度换算得来，这种尺度的不匹配性造成雷诺数的不匹配性，从而导致所受阻尼存在差异，在大初始位移的情况下这一差异被显著放大。3种模态总体来看，在大尺度初始位移的情况下，本文CFD结果与试验值的偏差均小于FAST，这在一定程度上说明CFD方法在处理较强非线性问题方面具有优势性，可以得出比FAST更为合理的结果。表 2也揭示出一致的结论，即本文垂荡、纵荡结果偏小于试验值，而纵摇结果则与试验值相吻合。

2.4 风力计算

本文利用风力计算原理，计算了表 3所列的5种定常风工况的风机叶片推力。每一种工况，风机整体均固定不动，风机转速在作业过程中均保持不变，叶片偏转角也固定为6.4°。所得推力由图 8给出，可以看到本文所得结果与试验值以及FAST所得结果吻合的很好。

3 动力响应分析 3.1 规则波

本文以表 4给示的6种规则波工况作为激励，研究风机在仅有规则波作用下的动力响应。为了使响应结果达到稳定状态，本文时域程序对于每一种工况均模拟了足够长时间(1 600 s)。

图 9为风机的纵荡、垂荡运动以及系缆张力RAO的数值模拟结果与试验值的对比。由图 9可知，对于纵荡运动RAO，6种工况下的各数值模拟结果与试验值均吻合得很好；对于垂荡运动RAO，工况1和2的各数值模拟结果也与试验值相吻合，而在工况3~6的情况下，FAST+OrcaFlex(A)、AQWA-C以及AQWA-D的结果与试验值较为接近，它们的共同点在于都是利用Morison方程来计及粘性影响，而利用CFD数值模拟或试验数据来计及粘性影响的本文、FAST以及FAST+OrcaFlex(B)的结果均明显偏小于试验值，而工况3~6所对应的波浪周期在风机的垂荡固有周期附近，可能是产生了共振而使垂荡RAO偏小；对于系缆张力RAO，各数值模拟结果基本都明显偏小于试验值，而考虑了系泊系统动力效应的AQWA-D、FAST+OrcaFlex(A)以及FAST+OrcaFlex(B)的结果虽然相对于未考虑系泊动力效应的其他数值模拟结果更接近于试验值，但相对于试验值仍是明显偏小的。

3.2 不规则波

本节以有义波高为H_s=11.3 m(超过了缅因湾的百年一遇有义波高^[5])的宽带白噪声谱作为激励，研究风机在仅有不规则波作用下的动力响应，模拟时长为3 h。相关模拟结果如图 10所示。

由图 10可知，本文的波浪PSD曲线与文献[5]是高度吻合的。在0.04~0.20 Hz波频范围内，本文纵荡响应结果与试验值吻合的很好；垂荡响应结果则较试验值则略有偏小；在纵荡固有频率0.009 3 Hz附近，本文结果明显偏小于试验值，原因在于忽略了二阶波浪力；而系缆2张力响应则与规则波作用下的结果一致，大幅偏小于试验值。

3.3 不规则波和风

本节讨论风浪联合作用下的风机动力响应，其中的不规则波激励与3.2节一致，而对于风则考虑了定常风(表 3的工况3)以及随机风(表 3的工况6)2种情况，模拟时长仍取为3 h，相关结果由图 11给出。

由图 11可知，本文所模拟的NPD风谱与理论值是相当吻合的。在0.04~0.20 Hz波频范围内，无论是定常风还是随机风，纵荡以及垂荡响应结果与无风时几乎没有差别，这说明风对于波频范围内的运动响应影响不大，而这一现象与文献[3, 4]的试验数据是吻合的。风会减小纵荡固有频率0.009 3 Hz附近的风机纵荡响应，虽然本文未考虑二阶波浪力的影响，此处的响应在无风状态下就已大幅小于试验值，但这一减小的趋势在定常风作用下时还是体现的很明显的，而对于随机风状态，此趋势则不太明显。而风对于系缆2张力响应的影响则是相当明显的，风力的存在相当于给系缆2附加了一个拉伸力，因此其张力响应会有所增大，而且风速越大增大的越明显。

4 结论

1) 本文所采用的势流理论为主，CFD数值模拟为辅的方法所得的动力响应结果大体上是与试验值相吻合的，从而其在总体上适用于风机的动力响应研究。

2) 风对于风机动力响应的影响主要体现在减小风机在纵荡固有频率附近的运动响应，而对于风机运动的波频响应影响不大；同时，风会影响系缆张力响应，而且风速越大影响效果越明显。

本文所采用的未考虑动力效应的系泊模型会低估风机的系缆张力响应，故其不适用于对系泊系统的研究设计。