0 引　言

深远海风能具有资源丰富、风力稳定、环境友好等诸多优点，近年来在全球范围内获得了广泛关注和积极发展^[1]。为了降低深远海风能开发成本，发电机组采用漂浮式基础结构已成为共识^[2]。当前，深远海风能开发主要面临机组服役条件十分恶劣和由于离岸较远而造成高昂维修保养费用等挑战^[3-4]。此外，与陆上风能和近海风能开发相比，深远海风能开发中支撑结构的费用急剧增长^[1]。上述因素使得当前浮式风机尚未完全实现商业化，限制了深远海风能的大规模开发。为此，应当采用基于整体可靠度的结构设计方法，实现兼顾浮式风机结构安全性与经济性的定量设计，为深远海风能的大规模开发提供关键技术保障。因此，浮式风机结构的整体可靠度分析具有重要意义。

当前，复杂结构的整体可靠度评估通常采用数值分析方法，主要涉及三个方面的问题：结构响应分析的数值模型、随机荷载的建模、结构可靠度分析方法^[5]。浮式风机结构的精细化数值分析模型是其整体可靠度分析的关键基础^[6]。由于浮式风机结构中同时存在多种复杂耦合效应^[7-8]，例如桨叶-机舱-塔筒-浮体-锚索等部件之间的相互作用、刚体运动与弹性变形耦合、气动-结构-机电-水动耦合、桨叶转动效应和几何非线性效应等，使得浮式风机结构的一体化动力学建模成为一项难题^[9]。虽然针对该建模问题当前已有较多研究，但是在实际工程中仍主要采用分离式设计方法，难以合理考虑上述复杂的耦合效应，也难以进行结构的整体可靠度分析。一般来说，浮式风机结构的一体化建模主要基于多体动力学理论^[9]，采用Lagrange方法或者Kane方法建立结构的整体动力学方程，并可以进一步分为多刚体动力学模型和多柔性体动力学模型^[10]。多刚体动力学模型将浮式风机中的所有部件均考虑为刚体，因而无法获得柔性部件（桨叶、塔筒等）的弹性变形和内力响应，具有较大的局限性^[10]。但是，由于多刚体动力学模型建模简便、分析效率较高，因此多用于浮式风机结构的概念设计和初始设计阶段^[11]。多柔性体模型能够同时考虑系统中的刚体运动和弹性变形响应，其中一般采用有限元方法或者模态叠加方法建立柔性部件的弹性模型^[7]。虽然多柔性体模型能够全面反映系统的力学特性和复杂的耦合效应，但是该类模型建模较为复杂，同时存在分析效率不高的问题，难以应用于需要进行大量数值模拟分析的问题中，例如结构的可靠度评估与优化设计等^[12]。

浮式风机在服役期间的主要环境荷载为风荷载和波浪荷载，且两种荷载具有相关性^[13]。因此，在浮式风机结构的动力分析与可靠度评估中应当合理考虑风和波浪荷载的联合作用^[14]。实测数据表明，10 min内的风速时程和20 min至3～6 h内的波浪时程均可以认为是平稳过程，并可以采用平均风速、有效波高、谱峰周期、风-浪夹角等物理量描述上述风况和海况^[13]。为了刻画结构在服役期内所经历的所有工况的概率大小，可基于结构所在场址的长期观测数据建立上述物理量的联合概率分布，这也是当前海洋结构设计规范中所推荐的方法^[13]。在建立风-浪参数的联合概率分布模型时，可采用条件概率分布建模方法^[15]、Nataf变换方法^[16]和Copula方法^[17-18]。条件概率分布建模方法是当前海洋结构设计规范中推荐使用的方法^[13]。该方法首先确定主导随机变量，并建立其边缘概率分布模型，然后根据主导随机变量的取值将数据划分为一系列区间，进而估计主导随机变量取不同值条件下的其他随机变量的概率分布。但是，当随机变量的个数较多时，该方法需要进行多次条件概率分布建模。此外，在某些条件下（例如高风速条件），观测数据较少，将导致该条件下的概率分布建模结果误差较大^[18]。Nataf变换方法首先基于等概率变换原则，将原始随机变量变换为服从标准高斯分布的随机变量，进而利用多维高斯分布模型建立变换后的随机变量的联合概率分布，最后再利用反变换获得原始随机变量的联合概率分布^[16]。该方法虽然简便易行，但是当随机变量为非线性、非高斯相关时，例如风-浪参数间的相关结构通常呈现这种特点，该方法将引入较大误差。Copula方法是近年来发展的多维随机变量联合概率分布建模的新方法，对于非线性相关问题十分有效^[17]。该方法在联合概率分布建模中，将随机变量的边缘概率分布建模与变量间的概率相关结构建模分开处理，具有简便、灵活、精确等优点。该方法目前已在多个领域获得广泛关注^[19]，但在海洋工程中尚应用较少。

结构可靠度分析方法主要包括一次和二次可靠度方法、基于跨越过程理论的首次超越破坏可靠度方法、随机模拟方法、概率密度演化理论等^{[5, 20]}。一次和二次可靠度方法常用于静力问题的分析，具有较高的效率，然而该方法需要将系统的极限状态函数展开为一阶或二阶形式，对于强非线性系统（例如浮式风机结构）会产生较大误差^[9]。基于跨越过程理论的首次超越破坏可靠度方法适用于动力系统，但通常引入Poisson假定，具有一定的近似性，尤其对于大失效概率问题其误差较大^[5]。随机模拟方法主要包括Monte Carlo模拟及其改进方法以及子集模拟等方法，这类方法不受问题的维度、非线性程度等因素的限制，应用范围较广，但通常具有随机收敛和效率不高等问题^[20]。为了将这类方法应用于浮式风机的可靠度分析中，通常发展浮式风机结构分析的代理模型以提高效率，但同时也引入了新的误差^[21]。概率密度演化理论给出了随机动力系统响应的概率密度函数的演化方程，揭示了系统响应随机性演化的内在物理规律^[5]，该方法在分析系统的可靠度时同样不受问题类型、问题维度和非线性程度等因素的影响，能给出系统响应量的概率密度函数演化过程，同时具有较高的分析效率，因此适用于浮式风机的可靠度分析，但目前在此领域的应用研究较少。

为了分析浮式风机结构的整体可靠度，首先引入结合多刚体动力学理论和有限元方法建立起来的浮式风机结构一体化分析模型（StoDRAFOWT模型）。该模型不仅能同时获得浮式风机的刚体运动和弹性变形响应，而且具有建模简便、易于求解等优点^[9]。然后，基于我国南海某场址长期风-浪数据并引入Copula方法，建立风-浪联合作用的多参数非线性相关的联合概率分布模型，进而采用概率密度演化理论分析浮式风机结构的整体可靠度。由此，建立了风-浪联合作用下浮式风机结构整体可靠度分析的基本框架。针对风机结构设计中需要考虑的额定风速和极端风剪不利工况，采用上述方法进行了单柱式浮式风机结构的整体可靠度分析。

1 浮式风机一体化动力学建模 1.1 基本建模思想与坐标系统

如前所述，浮式风机的多刚体动力学模型虽然分析效率较高，但是无法分析结构的弹性变形。而多柔性体动力学模型虽然能够全面分析结构的动力响应，但是建模复杂、分析效率不高。本文首先介绍将多刚体动力学模型与有限元方法相结合的建模策略，不仅能够全面分析结构响应，而且建模较为简便，分析效率较高^[9]。已有研究表明^[10-11]，对于当前兆瓦级的浮式风机系统，结构的弹性变形对系统刚体运动响应的影响可以忽略，因此可以首先采用多刚体动力学模型分析系统的刚体运动响应。而系统的刚体运动对结构的弹性变形具有重要影响。为此，采用Kane方法^[22-23]推导了考虑系统刚体运动效应的结构有限元模型，从而能够全面获得结构的动力响应。

图1为单柱式浮式风机的示意图和坐标系的定义。为描述浮式风机系统的状态，定义三组笛卡尔坐标系，分别为惯性坐标系 $O{x_0}{y_0}{z_0}$ 、浮体-塔筒的局部坐标系 $Oxyz$ 、桨叶的局部坐标系 $Ox'y'z'$ 。其中 $O{x_0}{y_0}{z_0}$ 坐标系的原点位于海平面位置， ${x_0}$ 轴与风速方向一致， ${z_0}$ 轴竖直向上。在初始时刻浮式风机未发生位移时，坐标系 $Oxyz$ 与 $O{x_0}{y_0}{z_0}$ 重合。坐标系 $Ox'y'z'$ 的原点位于轮毂中心， $x'$ 轴与 $x$ 轴平行， $z'$ 轴由叶根指向叶尖。

1.2 一体化结构分析模型

为了建立浮式风机系统的多刚体动力学模型，通常基于Lagrange方程得到考虑转子转动和浮体六自由度刚体运动的动力方程^[10]，即包含七个自由度。由于现代风力机通常配备桨距控制系统，为此，将变桨距控制方程与上述七自由度动力方程耦合，得到考虑变速、变桨距控制的浮式风机多刚体动力学模型^{[9, 24]}：

$ \left[ {\begin{array}{*{20}{c}} {{{\left[ {{{\boldsymbol{M}}_0}} \right]}_{6 \times 6}}}& {\begin{array}{*{20}{c}} 0&0&0 \\ 0&0&0 \\ 0&0&0 \\ {{I_{rx}}}&0&0 \\ 0&0&0 \\ 0&0&0 \\\end{array}}\\ {\begin{array}{*{20}{c}} 0&0&0&{{I_{rx}}}&0&0 \\ 0&0&0&0&0&0 \\ 0&0&0&0&0&0 \end{array}}& {\begin{array}{*{20}{c}} {{I_{rx}}}&0&0 \\ 0&1&0 \\ 0&0&1 \end{array}} \end{array}} \right] \left\{ \begin{gathered} {{{\dot {\boldsymbol{v}}}}_p} \\ {{\dot {\boldsymbol{\omega}} }_p} \\ {{\dot \omega }_b} \\ {\dot \varphi } \\ {\dot \beta } \\ \end{gathered} \right\} +$

$ \begin{split}& \left[ {\begin{array}{*{20}{c}} {{{\left[ {{{\boldsymbol{C}}_0}} \right]}_{6 \times 6}}}&{\begin{array}{*{20}{c}} \;\;\;\;\;\;0&\;\;\;\;\;\;\;\;\;\;0&\;\;\;\;\;0 \\ \;\;\;\;\;\;0&\;\;\;\;\;\;\;\;\;\;0&\;\;\;\;\;0 \\ \;\;\;\;\;\;0&\;\;\;\;\;\;\;\;\;\;0&\;\;\;\;\;0 \\ \;\;\;\;\;\;0&\;\;\;\;\;\;\;\;\;\;0&\;\;\;\;\;0 \\ {{I_{rx}}{\omega _z}}&\;\;\;\;\;\;\;\;\;\;0&\;\;\;\;\;0 \\ { - {I_{rx}}{\omega _y}}&\;\;\;\;\;\;\;\;\;\;0&\;\;\;\;\;0 \end{array}} \\ {\begin{array}{*{20}{c}} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{array}} & {\begin{array}{*{20}{c}}\\ {{c_0}} & 0 & 0 \\ { - 1} & 0 & 0 \\ {\dfrac{G}{{{\tau _p}}}\Bigg(\dfrac{{{\tau _d}{c_0}}}{{{I_{rx}}}} - 1\Bigg)} & {\dfrac{{ - G}}{{{\tau _p}{\tau _i}}}} & {\dfrac{1}{{{\tau _p}}}} \end{array}} \end{array}} \right]\left\{ \begin{gathered} {{{{{\boldsymbol{v}}}}}_p} \\ {{{{{\boldsymbol{\omega}} }}}_p} \\ {\omega _b} \\ \varphi \\ \beta \\ \end{gathered} \right\} \\& = \left\{ {\begin{array}{*{20}{c}} {{\boldsymbol {F}_G} + {\boldsymbol {F}_D}} \\ {{\boldsymbol {M}_G} + {\boldsymbol {M}_D}} \\ {{F_M}} \\ 0 \\ { - \dfrac{{G{\omega _0}}}{{{\tau _p}}}\Bigg(1 + \dfrac{t}{{{\tau _i}}}\Bigg) + \dfrac{{G{\tau _d}{F_M}}}{{{\tau _p}{I_{rx}}}}} \end{array}} \right\} \\[-45pt] \end{split} $

(1)

式中， $ {{\boldsymbol{M}}_0} $ 和 $ {{\boldsymbol{C}}_0} $ 分别为浮式风机的质量矩阵和阻尼矩阵，其表达式可参见文献[24]； $ {I_{rx}} $ 为转子相对于 $x'$ 轴的转动惯量； ${{{\boldsymbol{v}}}_p}$ 和 ${{\boldsymbol{\omega}} _p} = {\left( {{\omega _x},{\omega _y},{\omega _z}} \right)^{\text{T}}}$ 分别为浮体在 $Oxyz$ 坐标系中的平动和转动速度； $ {\omega _b} $ 和 $ \varphi $ 分别为桨叶的转动速度和转动角度； $ \beta $ 为桨叶的桨距角； $ {c_0} $ 、 $ G $ 、 $ {\tau _d} $ 和 $ {\tau _p} $ 是与风力机桨距控制和转矩控制有关的参数； $ {\omega _0} $ 为低转速轴的额定转速； ${{{\boldsymbol{F}}}_G}$ 为重力在 $Oxyz$ 坐标系下的效应； ${{{\boldsymbol{M}}}_G}$ 为在 $Oxyz$ 坐标系下重力关于坐标系原点的力矩； ${{\boldsymbol{F}}_D}$ 为外部荷载在 $Oxyz$ 坐标系下的效应； ${{{\boldsymbol{M}}}_D}$ 为在 $Oxyz$ 坐标系下外部荷载关于坐标系原点的力矩； $ {F_M} $ 为低转速轴上的净扭矩。

通过求解方程式（1）即可获得浮式风机中浮体的六自由度刚体运动响应、桨叶转速和桨叶的桨距角。为了分析浮式风机中柔性部件的弹性变形，基于Kane方程分别推导得到塔筒单元和桨叶单元的有限元表达式。

在坐标系 $Oxyz$ 中，可以得到考虑刚体运动效应的塔筒单元的动力方程^[9]：

$ {{\boldsymbol{M}}^t}\ddot {\boldsymbol{\mu}} + {{\boldsymbol{C}}^t}{\dot {\boldsymbol{\mu}} } + {{\boldsymbol{K}}^t}{{\boldsymbol{\mu }}} = {{{\boldsymbol{F}}}^t} $

(2)

式中： ${{\boldsymbol{M}}^t}$ 和 $ \;{{\boldsymbol{\mu}} } $ 分别为单元的质量矩阵和结点位移向量； ${{\boldsymbol{C}}^t}$ 、 ${{\boldsymbol{K}}^t}$ 和 $ {{{\boldsymbol{F}}}^t} $ 分别为单元的刚度矩阵、阻尼矩阵和荷载向量，表达式为：

$ {{\boldsymbol{C}}^t} = {\boldsymbol{C}}_0^t{\text{ + }}{\boldsymbol{C}}_d^t $

(3)

$ {{\boldsymbol{K}}^t} = {\boldsymbol{K}}_0^t{\text{ + }}{\boldsymbol{K}}_d^t $

(4)

$ {{{\boldsymbol{F}}}^t} = {{\boldsymbol{F}}}_0^t + {{\boldsymbol{F}}}_d^t $

(5)

其中： ${\boldsymbol{C}}_0^t$ 、 ${\boldsymbol{K}}_0^t$ 和 $ {{\boldsymbol{F}}}_0^t $ 分别为三维常规梁单元的阻尼矩阵、刚度矩阵和荷载向量； ${\boldsymbol{C}}_d^t$ 、 ${\boldsymbol{K}}_d^t$ 和 $ {{\boldsymbol{F}}}_d^t $ 分别为由于塔筒的刚体运动效应产生的附加阻尼矩阵、附加刚度矩阵和附加荷载向量，其表达式为：

$ {\boldsymbol{C}}_d^t = \int_V {2\rho {{{\boldsymbol{\tilde N}}}^{\text{T}}}{\boldsymbol{\tilde {\boldsymbol{\omega }}}}} {\boldsymbol{\tilde N}}{\text{d}}V $

(6)

$ {\boldsymbol{K}}_d^t = \int_V {\rho {{{\boldsymbol{\tilde N}}}^{\text{T}}}\left( {{\boldsymbol{\dot {\tilde {{\boldsymbol{\omega}}}} }} + {\boldsymbol{\tilde {{\boldsymbol{\omega}}} \tilde {{\boldsymbol{\omega}} }}}} \right)} {\boldsymbol{\tilde N}}{\text{d}}V $

(7)

$ \boldsymbol {F}_d^t = - \int_V {\rho {{{\boldsymbol{\tilde N}}}^{\text{T}}}\left[ {\left( {{\boldsymbol{\dot {\tilde {\boldsymbol{\omega}}} }} + {\boldsymbol{\tilde{\boldsymbol{ \omega}} \tilde {\boldsymbol{\omega}} }}} \right){{{\boldsymbol{r}}}_k} + {{\boldsymbol{R}}^{\text{T}}}{{{\ddot {\boldsymbol{r}}}}_o}} \right]} \;{\text{d}}V $

(8)

式中： $ \;\rho $ 为单元的材料密度； $ V $ 为单元体积； $ {\boldsymbol{\tilde N}} $ 为形函数矩阵；上标“T”表示转置矩阵； ${{\tilde {\boldsymbol{\omega}} }}$ 为浮体转动速度 ${{{\boldsymbol{\omega}} }_p}$ 对应的反对称矩阵^[24]； $ {\boldsymbol{R}} $ 为坐标转换矩阵，将坐标系 $Oxyz$ 中的矢量转换至坐标系 $O{x_0}{y_0}{z_0}$ 中； $ {{{\boldsymbol{r}}}_o} $ 为 $Oxyz$ 的坐标系原点在坐标系 $O{x_0}{y_0}{z_0}$ 中的位矢； $ {{{\boldsymbol{r}}}_k} $ 为单元中任意一点在坐标系 $Oxyz$ 中的位矢。

类似地，在坐标系 $Ox'y'z'$ 中，可以得到考虑刚体运动和旋转效应的桨叶单元的动力方程^[9]：

$ {{\boldsymbol{M}}^b}{\ddot {\boldsymbol{\mu}} } + {{\boldsymbol{C}}^b}{\dot {\boldsymbol{\mu}} } + {{\boldsymbol{K}}^b}{{\boldsymbol{\mu}} } = {{{\boldsymbol{F}}}^b} $

(9)

式中： $ {{\boldsymbol{M}}^b} $ 为单元的质量矩阵； $ {{\boldsymbol{C}}^b} $ 、 $ {{\boldsymbol{K}}^b} $ 和 $ {{{\boldsymbol{F}}}^b} $ 分别为单元的刚度矩阵、阻尼矩阵和荷载向量，其表达式与式（3）、式（4）、式（5）类似。其中，由于刚体运动和旋转效应产生的桨叶单元的附加阻尼矩阵、附加刚度矩阵和附加荷载向量的表达式为：

$ {\boldsymbol{C}}_d^b\; = \int_V {2\rho {{{\boldsymbol{\tilde N}}}^{\text{T}}}\left( {{{\boldsymbol{T}}^{\text{T}}}{\boldsymbol{\tilde {\boldsymbol{\omega}} T}} + {\boldsymbol{\tilde \varOmega }}} \right)} {\boldsymbol{\tilde N}}\;{\text{d}}V $

(10)

$ {\boldsymbol{K}}_d^b\; = \int_V {\rho {{{\boldsymbol{\tilde N}}}^{\text{T}}}\left[ {{{\boldsymbol{T}}^{\text{T}}}\left( {{\boldsymbol{\tilde {{\boldsymbol{\omega}}} \tilde {{\boldsymbol{\omega}}} }} + {\boldsymbol{\dot {\tilde {\boldsymbol{\omega}}} }}} \right){\boldsymbol{T}} + 2{{\boldsymbol{T}}^{\text{T}}}{\boldsymbol{\tilde {{\boldsymbol{\omega}}} T\tilde {{\boldsymbol{\omega}}} }} + \left( {{\boldsymbol{\tilde {\boldsymbol{\omega}} \tilde {\boldsymbol{\omega}} }} + {\boldsymbol{\dot {\tilde {\boldsymbol{\omega}}} }}} \right)} \right]{\boldsymbol{\tilde N}}} {\text{d}}V $

(11)

$ \begin{split} \boldsymbol {F}_d^b\; =& - \int_V {\rho {{{\boldsymbol{\tilde N}}}^{\text{T}}}\left\{ {{{\boldsymbol{T}}^{\text{T}}}{{\boldsymbol{R}}^{\text{T}}}{{\boldsymbol{\ddot r}}_o} + {{\boldsymbol{T}}^{\text{T}}}\left( {{\boldsymbol{\tilde {\boldsymbol{\omega}} \tilde {\boldsymbol{\omega}} }} + {\boldsymbol{\dot {\tilde {\boldsymbol{\omega}}} }}} \right){{\boldsymbol r}_{o'}}} \right\}\;{\text{d}}V} - \\& \int_V {\rho {{{\boldsymbol{\tilde N}}}^{\text{T}}}} \left[ {{{\boldsymbol{T}}^{\text{T}}}\left( {{\boldsymbol{\tilde {\boldsymbol{\omega}} \tilde {\boldsymbol{\omega}} }} + {\boldsymbol{\dot {\tilde {\boldsymbol{\omega}}} }}} \right){\boldsymbol{T}} + 2{{\boldsymbol{T}}^{\text{T}}}{\boldsymbol{\tilde {\boldsymbol{\omega}} T\tilde {\boldsymbol{\omega}} }} + \left( {{\boldsymbol{\tilde {\boldsymbol{\omega}} \tilde {\boldsymbol{\omega}} }} + {\boldsymbol{\dot {\tilde {\boldsymbol{\omega}}} }}} \right)} \right]{\boldsymbol{r}_q}{\text{d}}V \end{split} $

(12)

其中： $ {\boldsymbol{T}} $ 为坐标转换矩阵，将坐标系 $Ox'y'z'$ 中的矢量转换至坐标系 $Oxyz$ 中； $ {\boldsymbol{\tilde \varOmega }} $ 为桨叶转动速度 $ {\left( {{\omega _b},\;0,\;0} \right)^{\text{T}}} $ 对应的反对称矩阵； $ {{\ddot {\boldsymbol{r}}}_o} $ 为 $Oxyz$ 的坐标系原点在坐标系 $O{x_0}{y_0}{z_0}$ 中的加速度； $ {{{\boldsymbol{r}}}_{o'}} $ 为 $Ox'y'z'$ 的坐标系原点在坐标系 $Oxyz$ 中的位矢； $ {{{\boldsymbol{r}}}_q} $ 为单元中任意一点在坐标系 $Ox'y'z'$ 中的位矢。

值得注意的是，桨叶截面为非对称形状，并且桨叶截面存在沿翼展方向变化的预扭转角，同时还存在桨距控制。因此，桨叶在风轮平面内、外两个方向的振动相互耦合，在建模过程中应当合理考虑。此外，桨叶在风机正常运转工况下的变形较大，应当在桨叶的三维梁单元模型中考虑几何非线性效应。这是与塔筒有限元模型的不同之处。

基于塔筒和桨叶的有限元模型，并结合前述浮式风机的刚体运动响应，可建立浮式风机柔性构件的整体有限元模型，进而分析浮式风机结构弹性变形响应。

1.3 荷载计算

在浮式风机结构的动力响应分析中，主要应考虑桨叶和塔筒的气动荷载、浮体的水动力荷载和锚泊荷载。本文采用叶素动量理论方法^[25]计算桨叶的气动荷载，并考虑了Prandtl叶尖损失因子和高推力系数条件下的Glauert修正^[25]。塔筒的气动荷载可基于准定常假定计算^[25]。

在计算大型浮体的水动力荷载时，应当合理考虑辐射效应和绕射效应^[26]。本文关注的是单柱式浮式风机，浮体的截面尺寸相对较小，因此可以简化考虑这两种效应。在实际中可采用相对形式的Morison公式计算单柱式浮体的水动力荷载^{[10, 26]}，此时，单位长度浮体的水动力荷载为^[27]：

$ \begin{split} {\boldsymbol {f}_h} =& \frac{1}{2}{\rho _w}{d_p}{C_d}\left( {{{{\boldsymbol{v}}}_w} - {{{\boldsymbol{v}}}_f}} \right)\left\| {{{{\boldsymbol{v}}}_w} - {{{\boldsymbol{v}}}_f}} \right\| + \\& \frac{1}{4}{\rho _w}{\text{π }}d_p^2\left( {1 + {C_A}} \right){{{\boldsymbol{a}}}_w} - \frac{1}{4}{\rho _w}{\text{π }}d_p^2{C_A}{{{\boldsymbol{a}}}_f} \end{split} $

(13)

式中： $ {\rho _w} $ 为海水密度； $ {d_p} $ 为浮体外径； $ {C_d} $ 和 $ {C_A} $ 分别为拖曳力系数和附加质量系数； $ {{{\boldsymbol{v}}}_f} $ 和 $ {{{\boldsymbol{a}}}_f} $ 分别为该浮体微元的速度和加速度向量； $ {{{\boldsymbol{v}}}_w} $ 和 $ {{{\boldsymbol{a}}}_w} $ 分别为水质点的速度和加速度，可采用Airy波浪理论、Stokes波浪理论或其他波浪理论计算^{[13, 25]}； $ \left\| \cdot \right\| $ 表示向量的幅值。应指出，求解结构的动力方程（式1）时，通常将式（13）中的 $ - \left( {{1 \mathord{\left/ {\vphantom {1 4}} \right. } 4}} \right){\rho _w}{\text{π }}d_p^2{C_A}{{{\boldsymbol{a}}}_f} $ 移至动力方程的左边，由此产生浮体的附加质量矩阵。

在计算浮式风机的锚泊荷载时，需要建立锚索的力学分析模型，主要包括锚索的静力模型、拟静力模型和动力模型^[28]。由于静力和拟静力模型在某些工况下会低估锚索的疲劳荷载，因此推荐采用锚索的动力模型^[28]。本文采用锚索的集中质量动力模型^[28]计算锚泊荷载。该模型不仅考虑锚索的受拉弹性变形，而且将锚索的动力效应和水动力效应以及与海底的接触效应等也考虑在内。锚索底部固定在海床上，锚索顶部的边界条件由浮体的位置确定。通过求解锚索的动力方程，可以获得锚索对浮体的荷载效应。

1.4 模型验证

为了检验上述模型的正确性，分别利用该模型和当前风力机结构分析的常用软件FAST，分析风-浪联合作用下的海上单柱式浮式风机结构的动力响应，并进行对比。浮式风机结构参数采用美国国家可再生能源实验室（NREL）给出的5 MW海上风力发电机组参数^[29]和OC3项目中海上单柱式浮式风机基础参数^[30]，关键参数如表1所示。

本文仅给出一种工况条件下的分析对比结果，其他工况条件下的对比结果见文献[24]。此处的工况条件为：定常风速场，且轮毂高度处的风速为15 m/s；规则波浪，波高为6 m，波浪周期为10 s。分析结果如图2所示。结果对比表明，本文模型的分析结果与FAST的分析结果吻合良好，尤其是在响应进入平稳阶段以后。因此，可以利用本文模型进行浮式风机结构的可靠度分析。

2 风-浪联合作用的概率模型

由于浮式风机在服役期间同时受到风荷载和波浪荷载作用，因此在进行浮式风机的动力响应与可靠度分析时，应当合理考虑风-浪联合作用。如前所述，当前海洋结构的设计规范中建议采用风-浪参数的联合概率分布模型考虑风-浪联合作用^[13]。本文采用Copula方法建立风-浪参数的联合概率分布模型。

在实际工程中，通常需要考虑平均风速 ${V_w}$ 、有效波高 ${H_s}$ 、波浪谱峰周期 ${T_p}$ 以及风-浪夹角等物理量。本文主要关注额定风速条件下发生极端风剪时浮式风机结构的可靠度。根据风力机结构的设计规范，在分析该工况下风力机结构的动力响应时，取风-浪夹角为0°^[31]。因此，需要建立额定风速条件下 ${H_s}$ 和 ${T_p}$ 的条件联合概率分布模型。

在建立风-浪参数的联合概率分布模型时，应当基于场址的长期风-浪观测数据。本文采用欧洲中尺度天气预报中心ERA-Interim数据库^[32]提供的风-浪参数再分析数据，获取了我国南海场址（21°N，113°E）处1979年至2018年（共计40年）的风-浪再分析数据。这些数据均为30 min内观测数据的平均值，时间分辨率为6 h。通过对数据预处理，可得到轮毂处平均风速、有效波高和谱峰周期的长期数据。

为了建立额定风速条件下 ${H_s}$ 和 ${T_p}$ 的条件联合概率分布模型，在上述风-浪数据的基础上选取额定风速附近的数据进行概率建模。本文中，风力机的额定风速为11.4 m/s，因此取风速在11.3 m/s ～11.5 m/s区间内的风-浪数据进行 ${H_s}$ 和 ${T_p}$ 的联合概率分布建模，共有681个数据点。

根据Copula联合概率分布建模的基本理论^[17]，对于任意两个随机变量 ${X_1}$ 和 ${X_2}$ ，其联合概率密度函数的表达式为：

$ f\left( {{x_1},{x_2}} \right) = {c_{12}}\left( {{F_1}\left( {{x_1}} \right),{F_2}\left( {{x_2}} \right)} \right) \; {f_1}\left( {{x_1}} \right) \; {f_2}\left( {{x_2}} \right) $

(14)

式中： $ {F_1}\left( {{x_1}} \right) $ 和 $ {F_2}\left( {{x_2}} \right) $ 分别为 ${X_1}$ 和 ${X_2}$ 的边缘概率分布函数； $ {f_1}\left( {{x_1}} \right) $ 和 $ {f_2}\left( {{x_2}} \right) $ 为相应的边缘概率密度函数； $ {c_{12}}\left( { \cdot , \cdot } \right) $ 为二元Copula密度函数，可从常用的Copula函数中选用，例如高斯Copula、t Copula、Clayton Copula、Gumbel Copula、Frank Copula等^[17]。

从式（14）中可以看到，在利用Copula方法建立多维随机变量的联合概率分布模型时，需要分别建立各随机变量的边缘概率分布模型及变量间的相关结构模型。由于边缘概率分布模型及变量间的相关结构模型可以分开处理，因此该方法具有灵活、便捷、精确等优点。

基于上述681个额定风速条件下的数据点，首先采用常用的参数概率分布模型，建立 ${H_s}$ 和 ${T_p}$ 的边缘概率分布模型，并进行假设检验，进而基于AIC准则^[19]确定最优的边缘概率分布模型。结果表明， ${H_s}$ 和 ${T_p}$ 的最优边缘分布均为Gamma分布，结果见表2和图3。表中， $ \varGamma \left( \cdot \right) $ 为Gamma函数。

进一步，采用常用的Copula函数建立 ${H_s}$ 和 ${T_p}$ 的相关结构模型。备选的Copula函数除了上述常用的Copula函数外，还包括按照如下方式构建的非对称Copula函数^[33]：

$ \bar C\left( {u,v} \right) = C\left( {{u^{{{\alpha} _0}}},{v^{{\beta _0}}};\lambda } \right) \; {u^{1 - {{\alpha} _0}}} \; {v^{1 - {\beta _0}}} $

(15)

式中： $ \bar C\left( { \cdot , \cdot } \right) $ 为非对称Copula函数； $ C\left( { \cdot , \cdot } \right) $ 为上述常用的Copula函数； $\lambda $ 、 ${{\alpha} _0}$ 和 $\;{\beta _0}$ 为待定参数，且有 $0 \leqslant {{\alpha} _0} \leqslant 1$ 、 $0 \leqslant {\beta _0} \leqslant 1$ 。

在确定备选的Copula函数后，采用极大似然方法估算Copula函数中的待定参数，然后根据AIC准则选取最优Copula函数并进行拟合优度检验^[17，19]。对于上述额定风速条件下的风浪数据，识别的最优结果为按照式（15）构造的非对称高斯Copula模型，见表2。表中， $ C\left( {u,v} \right) $ 表示高斯Copula， $ {\varPhi _\lambda }\left( \cdot \right) $ 为相关系数为 $ \lambda $ ( $ - 1 \leqslant \lambda \leqslant 1 $ )的二维标准正态分布函数， $ \varPhi \left( \cdot \right) $ 为一维标准正态分布函数， $ {\varPhi ^{ - 1}}\left( \cdot \right) $ 为 $ \varPhi \left( \cdot \right) $ 的反函数。进而，基于式（14），可以得到额定风速条件下 ${H_s}$ 和 ${T_p}$ 的联合概率密度函数，如图4所示。

应指出，本文采用Copula方法建立了二维随机环境变量的条件联合概率分布模型。对于高维随机环境变量，可采用C-vine Copula方法便捷地建立其联合概率分布模型^[18-19]。

3 浮式风机结构整体可靠性分析 3.1 概率密度演化理论简介

本文采用概率密度演化理论进行额定风速下发生极端风剪时浮式风机的可靠度分析。为此，先简要介绍概率密度演化理论。

不失一般性，对于多自由度动力系统，考虑系统的随机因素后，动力方程可写为：

$ \dot {\boldsymbol{Y}}\left( t \right) = {\boldsymbol{G}}\left( {{\boldsymbol{Y}}\left( t \right),\;{\boldsymbol{\varTheta}} ,\;t} \right),\;\;\;\;{\boldsymbol{Y}}\left( {{t_0}} \right) = {{\boldsymbol{Y}}_0} $

(16)

式中： $ {\boldsymbol{Y}}\left( t \right) $ 为系统的状态响应量； $ {{\boldsymbol{Y}}_0} $ 为初始时刻的系统状态； ${\boldsymbol{ \varTheta}} = \left( {{\varTheta _1},\;...,\;{\varTheta _s}} \right)$ 为s维随机向量，刻画了系统中结构参数和外部激励的随机性。 ${\boldsymbol{ \varTheta}}$ 的联合概率密度函数记为 ${p_{_{\boldsymbol{ \varTheta}}} }\left( {{\theta}} \right)$ 。

式(16)的解可写为 ${\boldsymbol{Y}}\left( t \right) = {\boldsymbol{H}}\left( {{\boldsymbol{\varTheta}} ,\;t} \right)$ ，该系统中其他物理量的响应 $ {\boldsymbol{Z}}\left( t \right) = \left( {{Z_1}\left( t \right),\;...,\;{Z_m}\left( t \right)} \right) $ 均可通过 $ {\boldsymbol{Y}}\left( t \right) $ 和 $ \dot {\boldsymbol{Y}}\left( t \right) $ 获得。基于概率守恒原理，李杰和陈建兵^[5]建立了广义概率密度演化方程。当只关注一个物理响应量 $ Z\left( t \right) $ 时，广义概率密度演化方程具有如下形式^[5]：

$ \frac{{\partial {p_{_{Z{{ \varTheta}} }}}\left( {z,\;{{\theta}} ,\;t} \right)}}{{\partial t}} + \dot Z\left( {{\boldsymbol{\varTheta}} ,\;t} \right)\frac{{\partial {p_{_{Z{{ \varTheta}} }}}\left( {z,\;{{\theta}} ,\;t} \right)}}{{\partial z}} = 0 $

(17)

且具有如下的初值条件：

$ {p_{_{Z{{\varTheta}}} }}\left( {z,\;{{\theta}} ,\;t} \right)\left| {_{t = 0}} \right. = \delta \left( {z - {z_0}} \right) {p_{_{{\varTheta}} }}\left( {{\theta}} \right) $

(18)

式中： ${p_{_{Z{{\varTheta}}} }}\left( {z,\;{{\theta}} ,\;t} \right)$ 为 $\left( {{\boldsymbol{\varTheta}} ,\;Z\left( t \right)} \right)$ 的联合概率密度函数； $ \delta \left( \cdot \right) $ 为Dirac函数； $ {z_0} $ 为 $ Z\left( t \right) $ 的初值。

因此， $ Z\left( t \right) $ 的概率密度函数为：

$ {p_{_Z}}\left( {z,\;t} \right) = \int_{{\varOmega _{{ \varTheta}} }} {{p_{_{Z{{ \varTheta}}} }}\left( {z,\;{{\theta}} ,\;t} \right){\text{d}}{{\theta}} } $

(19)

式中， ${\varOmega _{{\varTheta}} }$ 为随机向量 ${\boldsymbol{\varTheta}}$ 的分布域。

为了求解广义概率密度演化方程，已发展了点演化和群演化两类数值求解方法。本文采用点演化方法，具体的求解步骤可参考文献[5]。

基于概率密度演化理论发展了两种结构首次超越破坏可靠度分析方法，即，基于吸收边界条件的可靠度分析方法和基于极值分布的可靠度分析方法^[5]。本文采用基于极值分布的可靠度分析方法。首先构造随机过程 $Z\left( {{\boldsymbol{\varTheta}} ,\;t} \right)$ 的等价极值事件：

$ {Z_{{\text{ext}}}}\left( {\boldsymbol{\varTheta}} \right) = \mathop {\max }\limits_{t \in \left[ {0,T} \right]} \Big| {Z\left( {{\boldsymbol{\varTheta}} ,\;t} \right)} \Big| $

(20)

式中， $ T $ 为结构响应的分析时间。

进而构造如下的虚拟随机过程^[5]：

$ {Z_{{\text{vir}}}}\left( {{\boldsymbol{ \varTheta}} ,\;t} \right) = {Z_{{\text{ext}}}}\left( {\boldsymbol{ \varTheta}} \right)\sin \left( {\frac{{\text{π }}}{2}t} \right),\;\;t \in \left[ {0,1} \right] $

(21)

由于 $Z\left( {{\boldsymbol{\varTheta}} ,\;t} \right)$ 的极值分布等价于 $t = 1$ 时刻 ${Z_{{\text{vir}}}}\left( {{\boldsymbol{\varTheta}} ,\;t} \right)$ 的概率分布，因此可先通过求解广义概率密度演化方程获得 ${Z_{{\text{vir}}}}\left( {{\boldsymbol{\varTheta}} ,\;t} \right)$ 的概率密度演化过程，然后获得 $Z\left( {{\boldsymbol{\varTheta}} ,\;t} \right)$ 的极值分布，即：

$ {p_{_{Z_{\max }}}}\left( z \right) = {p_{_{Z_{{\text{vir}}}}}}\left( {z,t} \right)\left| {_{t = 1}} \right. $

(22)

由此，可得到结构的可靠度为：

$ {F_s} = \int_{ - \infty }^b {{p_{_{Z_{\max }}}}\left( z \right){\text{d}}z} $

(23)

式中， $b$ 为 $Z\left( {{\boldsymbol{\varTheta}} ,\;t} \right)$ 的失效阈值。

3.2 算例

针对风力机结构设计规范^[31]中定义的额定风速条件下发生极端风剪的工况，进行了单柱式浮式风机结构的整体可靠度分析。此处仅考虑竖直方向的正极端风剪，此时风剖面为^[31]：

$ V\left( {z,t} \right) = \left\{ \begin{array}{l} {V_{{\text{hub}}}}{\left( {\dfrac{z}{{{z_{{\text{hub}}}}}}} \right)^{\alpha} } + \left[ {1 - \cos \left( {\dfrac{{2{\text{π}}t}}{{{T_0}}}} \right)} \right] \left( {\dfrac{{z - {z_{{\text{hub}}}}}}{D}} \right) \cdot \\ \qquad \left[ {2.5 + 0.2\beta {\sigma _1}{{\left( {\dfrac{D}{{{\varLambda _1}}}} \right)}^{\frac{1}{4}}}} \right],\;\;0 \leqslant t \leqslant {T_0} \\ {V_{{\text{hub}}}}{\left( {\dfrac{z}{{{z_{{\text{hub}}}}}}} \right)^{\alpha} }\;\;\;\quad\qquad\qquad,\;\;{\text{otherwise}} \end{array} \right. $

(24)

式中： ${V_{{\text{hub}}}}$ 为轮毂高度处的平均风速； $z$ 和 ${z_{{\text{hub}}}}$ 分别为高度和轮毂高度； ${\alpha} = 0.2$ 为风剖面指数律的幂指数； ${T_0} = 12\;{\text{s}}$ 为极端风剪的持续时间； $t$ 为时间； $D$ 为风轮直径； $\;\beta = 6.4$ ； ${\sigma _1}$ 为脉动风速标准差的代表值， ${\sigma _1} = {I_{{\text{ref}}}}\left( {0.75{V_{{\text{hub}}}} + 5.6} \right)$ ，其中 ${I_{{\text{ref}}}}$ 为湍流强度； ${\varLambda _1}$ 为湍流积分尺度； ${\sigma _1}$ 和 ${\varLambda _1}$ 的取值可根据规范^[31]确定。当 ${V_{{\text{hub}}}} = 15\;{\text{m/s}}$ 、 ${I_{{\text{ref}}}} = 0.12$ 时，5 MW风力机在 $t = {{{T_0}} \mathord{\left/ {\vphantom {{{T_0}} 2}} \right. } 2}$ 时的极端风剖面和正常风剖面对比图见图5。

在分析结构的可靠度时，应当充分、合理考虑系统中的随机源。本文采用的概率密度演化理论在分析结构可靠度时能够同时全面考虑结构自身和外部激励的随机性^[5]。作为示例，本文不仅将额定风速条件下的有效波高 ${H_s}$ 和谱峰周期 ${T_p}$ 作为随机变量，而且考虑了脉动风速、波面高程以及塔筒阻尼比 ${\xi _T}$ 的随机性。本文基于谱表达方法并分别采用波数-频率联合功率谱^[34-35]和JONSWAP谱^[13]模拟随机脉动风场和随机波浪场。为了降低随机场模拟中的随机变量（相位）个数，采用演化相位谱方法^[36]，由此分别引入4个和1个相位演化时间，分别记为 ${T_{e1}}\sim{T_{e4}}$ 和 ${T_{e5}}$ 。相位演化时间为随机变量，基于实测风速数据和波浪数据可以识别其概率分布^[36]。相位演化时间和塔筒阻尼比的概率分布信息如表3所示^[9]。 ${H_s}$ 和 ${T_p}$ 的概率分布信息由第2节方法确定。

在求解广义概率密度演化方程时，首先需要将概率空间剖分为 ${N_q}$ 个子域，并在每个子域内选取一个代表点^[5]。然后对每个代表点进行确定性动力响应分析以获得概率密度演化方程（偏微分方程）中的时变系数，进而求解概率密度演化方程，并获得响应量的概率密度演化结果^[5]。本文共考虑8个随机变量，采用基于点集GF-偏差最小的点集优选策略^[37]在八维概率空间中选取了329个代表点。图6中的红点展示的是 ${H_s}$ 和 ${T_p}$ 的代表点，图中的灰线为 ${H_s}$ 和 ${T_p}$ 联合概率密度函数的等值线。

在选定代表点之后，需要模拟每种情况下的随机风场和波浪场，并分析浮式风机结构的动力响应。作为示例，图7给出了 ${H_s} = 2\;{\text{m}}$ 、 ${T_p} = 8.4\;{\text{s}}$ 工况下浮式风机结构动力响应的典型样本。同时分析了不考虑极端风剪时的结构响应并进行对比。动力响应的时间分析长度为1500 s，并从第500 s开始施加极端风剪条件。为了避免初始瞬态响应的影响，仅展示后1200 s的动力激励和结构动力响应结果。图中，“EWS”表示极端风剪。

从图7（a）可以看出，在施加极端风剪条件后， $z = 30\;{\text{m}}$ 处和 $z = 150\;{\text{m}}$ 处的风速差异明显增大，由此导致塔筒和桨叶的位移响应发生变化（如图7d、图7f所示）。从图7（d）还可以看出，考虑极端风剪后，该样本中塔顶位移响应的极值有所增大，因此会影响结构的整体可靠度结果。图7（c）和图7（f）显示，大约在700～750 s时间段内，结构响应有突然的“跌落”，这是由于启动了变桨控制，降低了气动荷载（如图7（j）所示）。图7（g）～图7（j）表明，极端风剪对浮式风机的纵荡响应、纵摇响应、桨叶转速、桨叶的桨距角变化均影响较小。

采用概率密度演化理论可以获得结构响应极值的概率分布。以塔顶顺风向位移为例，通过求解广义概率密度演化方程，得到其响应极值的概率密度函数如图8所示。在给定的失效阈值条件下，当塔顶位移响应超过该阈值时，即认为发生失效。由此可以得到不同失效阈值条件下的失效概率，如图9所示。

作为示例，本文在分析浮式风机结构的可靠度时，仅考虑了塔顶位移。在实际中，应考虑多个指标，例如塔顶加速度、塔底应力以及浮式风机的纵摇响应等。概率密度演化理论能够同时考虑多个失效指标^[5]。当同时采用多个失效指标分析时，可靠度结果可能发生变化。

4 结　论

浮式风机结构的整体可靠性分析是进行基于可靠度优化设计的关键基础，对于大规模深远海风能的开发具有重要意义。本文针对海上单柱式浮式风机结构，首先介绍了一体化耦合动力学分析模型，然后建立了环境参数的联合概率分布模型以合理考虑风-浪联合作用，进而分析浮式风机结构在极端工况条件下的动力可靠度。主要研究结论如下：

1）对于当前兆瓦级浮式风机，柔性部件的弹性变形对浮式风机的刚体运动响应几乎没有影响，因此在建模过程中可以忽略这一因素，从而使建模方法与分析过程更为简便。本文基于这一思想，将多刚体动力学理论与有限元方法结合起来考察了浮式风机结构的一体化模型，不仅能够弥补一般多刚体动力学模型的不足，而且利用Kane方法便捷地考虑了系统中刚体运动对柔性部件变形的影响。虽然从理论上看，全耦合多柔性体模型（例如FAST）的建模方法考虑得更为全面，但是从实际工程应用的角度来看，本文所提的方法更为便捷。两种模型的对比分析结果表明了本文模型具有良好的精度。

2）在浮式风机结构的动力响应分析中，应当合理考虑风-浪联合作用。基于场址的长期风-浪观测数据并采用Copula方法，能够精确、便捷地建立风-浪多参数非线性相关的联合概率分布模型，从而合理地考虑风和波浪对结构的荷载效应。

3）基于概率密度演化理论的可靠度分析方法适用于海上大型浮式风机结构的整体可靠性分析。该方法不仅能够同时考虑外部激励和结构自身的随机性，而且能够获得结构响应的概率密度演化过程，同时具有很高的分析效率。

本文研究建立了海上大型浮式风机结构整体可靠性分析的基本框架，为进行基于整体可靠度的结构优化设计奠定了基础。由于我国东南海域的风电场常年遭受台风侵袭^[38]，因此在未来工作中需要进一步分析在台风等极端环境条件下的浮式风机结构整体可靠度，并开展基于整体可靠度的结构优化设计等研究。