目前已运行的海上风电场中，风机基础型式采用单桩基础的比例近80%^[1]，是应用最广泛的基础型式. 海上风机的运行环境恶劣，风机基础结构长期受到循环载荷作用，易发生疲劳破坏^[2-4].

为了准确评价海洋结构物的疲劳损伤，需要进行全时域的动力分析. Jia^[5]对波浪载荷作用下的导管架结构进行了时域疲劳损伤分析，为缩减计算时间对工况进行了合并简化. Kvittem等^[6]在进行时域疲劳分析时，发现当模拟时长大于1 h时，模拟时段内随机性对损伤计算的影响降低.

全时域计算非常耗时，频域疲劳分析则相对简便. ABS规范给出了窄带响应下的疲劳损伤频域计算方法，Du等^[7]利用该方法计算了波浪载荷作用下的单桩基础海洋平台的疲劳损伤，并将其与时域分析结果对比，发现频域法计算的疲劳损伤是时域法结果的3.4倍. Mohammadi等^[8]认为ABS规范频域法计算结果过于保守的原因在于实际的应力幅概率分布并不严格符合ABS规范所假定的Rayleigh分布. 为此，一些学者提出适用于宽带过程的疲劳损伤分析方法^[9-12]. Benasciutti等^[13]研究了频域法的宽带过程疲劳问题，并讨论了带宽参数的取值与组合对疲劳损伤的影响. 双峰过程是宽带过程的1种特例^[14]，可由组合谱法^[14-17]和疲劳损伤组合法^[18-20]计算.

在进行海上风电结构疲劳分析时，可以借鉴海洋平台的有关研究成果，同时面临一些新问题. 与海洋平台相比，海上风电结构除受到波浪荷载外，风荷载及其引起的弯矩作用更加突出，气动阻尼影响更大，导致海上风电结构的动力响应和疲劳损伤机制更加复杂.

当前在海上风电结构疲劳设计时，为了简化载荷输入和减小计算量，通常的做法是分别计算风、浪单独作用时基础的疲劳损伤，再对损伤值进行叠加作为基础的总疲劳损伤. Yamashita等^[21]将海上单桩基础风机简化为多质点体系，计算表明风、浪联合作用的疲劳损伤比单独作用、损伤叠加得到的结果大25%. Dong等^[22]计算了导管架基础的疲劳损伤，结果表明疲劳损伤主要由风载荷引起，且考虑风、浪联合作用比单独作用计算得到的疲劳损伤大12.79%~24.60%. Yeter等^[23]采用频域方法分析了三桩式海上风机基础的疲劳损伤. Tempel^[24]分别采用时域和频域方法计算了海上风机基础疲劳损伤，发现由于2种方法得到的应力谱不同，损伤结果存在明显差异. Ziegler等^[25]采用频域方法分析了海上风机结构的浪致疲劳损伤，认为水深和波浪周期对浪致疲劳载荷具有显著影响. Marino等^[26]发现风机运行状态不同会有不同的气动阻尼，对疲劳损伤计算结果影响较大. 阻尼是一个非常复杂的因素，不同文献和风机运行工况的阻尼比取值差别较大，气动阻尼比取值范围在3%~9%^[27-28]，材料阻尼比取值范围在0.19%~1.5%^[29-30]，基础阻尼比取值范围在0.17%~1.3%^[30-31]. 由于实际风况和海况条件的复杂性，以及各海上风电文献所采用的风、浪载荷和计算参数不同，不同时域和频域计算方法之间难以相互验证和比较，难以评估不同方法之间的误差程度.

在工程设计中，疲劳寿命往往是风机基础尺寸选择的控制性指标，但当前的不同计算模型和方法所得的疲劳损伤差异大，非常不利于进行海上风机基础优化设计，可能会造成结构设计上的浪费和不安全. 本文以海上风机单桩基础为例进行疲劳分析，研究气动阻尼比取值对基础疲劳损伤的影响；建立全时域的动力分析模型，评价风、浪联合作用对基础疲劳损伤的影响；建立海上风机基础频域疲劳分析流程，评价不同应力幅概率分布模型的影响.

风是一种典型的随机过程，一般认为风由平均风和脉动风两部分组成，风速可表示为平均风速和脉动风速之和.

平均风速沿高度方向的分布由风廓线描述，可由指数分布表示；认为平均风速长时间的概率分布满足威布尔分布. 通常假定脉动风是平稳的高斯过程，可以用功率谱的形式描述风在不同频率下的能量分布，常用的风谱包括Kaimal谱、Von Karman谱、Davenport谱等. 本文使用IEC标准^[32]中推荐的Kaimal谱来模拟脉动风速.

风作用在叶轮上的水平推力（风机载荷）由推力系数法确定，推力系数C_T由叶素-动量理论计算^[33]. 为简化计算，本文只考虑顺风向的风载荷：

(1) $F_{\rm{T}} = \frac{1}{2}{C_{ _{\rm{T}}}}\rho A{u^2}. $

式中：F_T为水平推力，ρ为空气密度，A为叶轮扫风面积，u为瞬时风速.

除了叶轮的水平推力外，风载荷还包括风作用在塔架、机舱、轮毂上的载荷，计算方法为

(2) $F = \frac{1}{2}{\mu _{\rm{s}}}{A_{\rm{t}}}\rho {u^2}. $

式中：μ_s为体型系数，A_t为塔架、机舱、轮毂在垂直于风向平面上的投影面积.

海上风机单桩基础的直径与波长之比D/λ<0.2，属于小尺度桩，在计算波浪对桩的作用力时可以不考虑桩对波场的影响，一般使用莫里森公式计算波浪载荷. 莫里森公式认为波浪力是速度分力F_D和拖拽分力F_I之和：

(3) $F = {F_{\rm{D}}} + {F_{\rm{I}}}. $

本文采用莫里森公式结合Pierson-Moscowitz波浪谱的方法计算波浪载荷^[34].

在进行动力分析时，风机基础与塔筒均采用欧拉-伯努利梁模拟，桩土作用由p-y、t-z、Q-z弹簧模拟，分别表示水平土抗力与桩身水平位移、桩身摩阻力与桩身竖向位移、桩端土抗力与桩端竖向位移的关系. 将机舱、叶轮、轮毂简化为一集中质量点添加在梁顶位置. 将气动阻尼添加在集中质量点处，塔筒与基础分别赋予材料阻尼和基础阻尼. 海上风机动力分析模型如图1所示. 图中，k为弹簧刚度.

在频域疲劳分析中，将结构看作是一个动力学系统，把风、浪载荷等随机过程看作是系统的输入，将应力等结构的响应看作系统的输出，系统的输入和输出由功率谱密度函数（power spectral density，PSD）表示，两者之间的变换关系由传递函数决定：

(4) ${S_{yy}}\left( \omega \right) = {\left| {{T_{xy}}\left( \omega \right)} \right|^2}{S_{xx}}\left( \omega \right). $

式中：S_xx（ω）为输入谱，S_yy（ω）为输出谱，T_xy（ω）为传递函数（transfer function，TRF），是系统（结构）本身所具有的特性. 只要明确传递函数，即可根据输入谱求得输出谱.

在得到应力谱后，如何处理应力谱是频域法计算疲劳损伤的一个关键问题. 通常是利用应力幅的概率分布函数计算结构在各应力水平下的循环次数. 对于平稳的高斯随机过程，计算应力幅概率密度函数（probability density function，PDF）应用较多的包括Rayleigh模型、Tunna模型、Benasciutti-Tovo模型、Dirlik模型等.

Dirlik模型采用1个指数分布和2个Rayleigh分布来近似描述应力幅的概率分布，其概率密度函数为

(5) $ \begin{split} p\left( {\Delta \sigma } \right){\rm{ = }}& \frac{{\rm{1}}}{{2\sqrt {{m_0}} }}\left[ {\frac{{{D_1}}}{Q}\exp \left( { - \frac{Z}{Q}} \right) + } \right.\\ & \left. {\frac{{{D_2}Z}}{{{R^2}}}\exp \left( { - \frac{{{Z^2}}}{{2{R^2}}}} \right) + {D_3}Z\exp \left( { - \frac{{{Z^2}}}{2}} \right)} \right], \end{split} $

(6) $Z = {{\Delta \sigma } / {(2\sqrt {{m_0}} )}},$

(7) ${D_1} = {{2\left( {{\chi _{_ {\rm{m}}}} - {\alpha _2}^2} \right)}/{\left( {1 + {\alpha _2}^2} \right)}},$

(8) ${D_2} = {{\left( {1 - {\alpha _2} - {D_1} + {D_1}^2} \right)} / {\left( {1 - R} \right)}},$

(9) ${D_3} = 1 - {D_1} - {D_2},$

(10) ${\alpha _2} = {{{m_2}} / {\sqrt {{m_0}{m_4}} }}$

(11) ${\chi _{_ {\rm{m}}}} = {{{m_1}}}/{{{m_0}}}\sqrt {{{{m_2}}}/{{{m_4}}}} ,$

(12) $R = \frac{{{\alpha _2} - {\chi _{_ {\rm{m}}}} - {D_1}^2}}{{1 - {\alpha _2} - {D_1} + {D_1}^2}},$

(13) $Q = {{1. 25\left( {{\alpha _2} - {D_3} - {D_2}R} \right)} / {{D_1}}}. $

式中：D₁、D₂、D₃分别为1个指数分布和2个Rayleigh分布的权重系数，Z、α₂、χ_m、R、Q为m_i的函数，m_i为应力谱的i阶矩，表示为

(14) ${m_i} = \int_0^{ + \infty } {{f^i}} S\left( f \right){\rm{d}}f. $

得到应力幅的概率密度函数后，使用期内结构在不同应力幅下的应力循环次数由下式计算.

(15) ${n_i}\left( {\Delta \sigma } \right) = \sqrt {{{{m_4}} / {{m_2}}}} {T_{\rm{d}}}{p_i}\left( {\Delta \sigma } \right). $

式中：T_d为结构的设计使用寿命，p_i(∆σ)为不同应力幅的出现频率.

对海上风机而言，风导致的随机过程是宽带过程，波浪导致的随机过程是窄带过程. 概率密度函数模型中Rayleigh模型通常用于窄带过程，其余用于宽带过程.

除了利用概率密度函数计算应力循环次数外，还可以通过对应力谱进行快速傅里叶逆变换（inverse fast Fourier transform，IFFT）得到应力时程，然后采用雨流计数法统计应力循环次数. IFFT过程由谐波叠加法实现：

(16) $\sigma \left( t \right) = \sum\nolimits_{i = 1}^N {{A_i}} \cos\; \left( {{f_i}t + {\theta _i}} \right),$

(17) ${A_i} = \sqrt {2S\left( {{f_i}} \right)\Delta f} ,$

(18) $ {f_i} = {f_{\min }} + {{\left( {2i - 1} \right)\Delta f} / 2}. $

式中：σ为应力，t为时间，f_i为第i个频带的频率，θ_i为随机初始相位角，A_i为随机应力幅，∆f为频率间隔，f_min为第i个频带的频率最小值.

如图2所示为频域法疲劳评价的分析流程，图中，S_stress（f）为应力谱，S_load（f）为载荷谱，D_wind为风致疲劳损伤，D_wave浪致疲劳损伤，P为联合分布概率. 具体步骤如下：1）根据风与波浪的载荷时程，在Matlab软件中分别计算出载荷的功率谱密度函数（PSD），见图中①和②；2）在ABAQUS软件的稳态动力学分析步中进行谐响应分析，输出应力的传递函数（TRF），见图中③；3）分别计算风与波浪作用下的应力功率谱密度函数（PSD），见图中④和⑤；4）对于风致应力谱（宽谱），由Dirlik模型计算应力幅的概率密度函数（PDF），得到各应力水平下的应力循环次数，见图中⑥；5）对于浪致应力谱（窄谱），通过IFFT方法将应力谱转换为应力时程，由雨流计数法统计应力循环次数，见图中⑦.

以某海上风电场3.3 MW风机为例，该风机单桩基础直径D_P=5.5 m，壁厚t_P=75 mm，桩长L_P=83 m. 泥面以下桩长为55 m，泥面至水面距离为17 m，水面至基础顶的距离为11 m. 塔筒为一变截面结构，长度L_t=78.67 m，塔底直径D_tb=5.5 m，t_tb=30 mm，塔顶直径D_tt=3.316 m，壁厚t_tt=38 mm. 轮毂高度为92 m. 机舱、叶轮、轮毂的质量为218.28 t. 风机基础与塔筒选用DH36型钢材，弹性模量E=210 GPa，密度ρ_s=7 850 kg/m³，屈服强度为f_y=355 MPa. 基础上每间隔1 m采用一组弹簧来模拟桩土相互作用.

风机的运行风速为2.5~20 m/s，共划分为9种风况（9个风速区间Δu），每种风况（W_i）选取一个风速代表值u参与计算. 依据风电场的风资源特征参数，威布尔形状因子l=2.29，尺度因子c=8. 46，确定平均风速的概率分布函数，得到各风况划分结果如表1所示. 利用2 a的波浪实测资料，得到有效波高H_s与平均周期T的联合分布概率P，由此划分的海况（S_i）如表1所示. 组合风况与海况，得到共计9×9=81个风、浪联合作用下的工况，各工况出现的概率如图3所示.

频域法的有限元模型与时域法相同，在稳态动力学分析步中进行谐响应分析，气动阻尼比ξ_aero=3%，材料阻尼比ξ_m=2%，基础阻尼比ξ_fdn=3%^[30-31]，分别得到风、浪载荷单独作用下计算点（本文计算点取基础−31 m截面处，下同）的应力传递函数（TRF）. 在Matlab软件中计算各工况下风、浪载荷的功率谱密度函数（PSD）. 由式（4）得到各工况下计算点的应力功率谱密度函数（PSD）. 以风况W₆和海况S₈为例对比由频域法和时域法得到的应力谱，如图4所示.

图4中时域法与频域法所得的应力谱均在0.27 Hz处出现尖峰，即结构的一阶固有频率处，虽然尖峰处的幅值存在偏差，但整体吻合良好，说明频域法具有合理性.

以W₄和S₉为例，分别采用时域法和频域法讨论不同气动阻尼比取值对基础疲劳损伤的影响. 本例中气动阻尼比ξ_aero取值为0、1%、2%、3%、4%、5%，材料阻尼比ξ_m与基础阻尼比ξ_fdn保持不变，取值与3.2节相同. 分析时考虑疲劳设计因子（design fatigue factor，DFF），DFF在空气中取1，水面以下取3. 疲劳损伤D以对数坐标表示，如图5和6所示.

时域法与频域法结果均表明，基础疲劳损伤对气动阻尼比的取值十分敏感，不同阻尼比取值下的基础疲劳损伤存在很大差异. 例如当气动阻尼比取2%时，时域法计算得到的基础风致和浪致疲劳损伤分别是气动阻尼比取3%时的1.9、1.5倍；频域法计算得到的ξ_aero=2%时基础风致和浪致疲劳损伤分别是ξ_aero=3%时的2.9、1.9倍. 图5中由时域法和频域法计算得到的风致疲劳的差值随气动阻尼比的增大而减小；图6中时域法和频域法计算得到的浪致疲劳的差值受气动阻尼比取值的影响不大. 气动阻尼比对风致动力响应的影响更加突出.

仍以W₄和S₉为例，叠加3.3节中由频域法得到的不同气动阻尼比取值的风致疲劳损伤和浪致疲劳损伤作为总疲劳损伤，并由时域法计算W₄和S₉联合作用下的基础总疲劳损伤，阻尼比取值与3.3节相同，结果以对数坐标表示，如图7所示. 由图7可见，随着气动阻尼比的增大，2种方法计算得到的疲劳损伤都在降低，当阻尼比大于1%时，频域法风、浪疲劳损伤先分算再叠加的结果降低的幅度更大一些，使得2种方法所得结果之间的误差越来越大.

采用时域法计算全部工况的风、浪单独作用和风、浪联合作用时基础在25 a内的疲劳损伤，阻尼比取值与3.2节相同. 风、浪单独作用时结果如表2所示，疲劳损伤贡献率由K表示，风致疲劳损伤为0.461，浪致疲劳损伤为8.50×10⁻³，两者叠加得到总疲劳损伤为0.470；风、浪联合作用时各工况下的结果如图8所示，总疲劳损伤为0.829.

由表2可见，本例中风致疲劳损伤占主导作用. 风速越大，风致疲劳损伤越大，例如W₅、W₆、W₇的出现概率为23%，K=89.87%；波高越大，浪致疲劳损伤越大，例如S₈、S₉所示，波高相同时，周期越小，浪致疲劳损伤越大，例如S₆、S₇所示.

对比表2和图8结果可知，风、浪单独作用比联合作用的基础疲劳损伤小43.3%，表明风、浪联合作用对基础疲劳损伤具有较大影响，分析时不能忽视. 以S₂与各风况组合的情况为例，分别计算风、浪单独作用后叠加，以及风、浪联合作用的损伤值，如图9所示. 结果表明，基础疲劳损伤均随风速的增大而增大，与单独作用相比，风、浪联合计算总是产生更大的疲劳损伤. 以大风速的W₇与各海况组合的情况为例，同理得到风、浪单独作用后叠加，以及风、浪联合作用的损伤值，如图10所示. 结果表明，大风速时的基础疲劳损伤主要由风况与海况组合后出现概率较大的工况贡献，例如当W₇与S₂、S₃、S₅组合时（出现概率如图3所示）.

另外，如果只将图8中的W₄、W₅、W₆、W₇与全部海况进行组合，可将原先9×9工况矩阵缩减为4×9矩阵，则近似地可以减少55.56%的计算量，得到的基础疲劳损伤为7.79×10⁻¹，与原计算损伤值的计算误差为6.03%，这将十分有助于减少风、浪联合作用时的计算量.

频域法中应力幅概率分布模型的选取是影响疲劳损伤计算准确性的一个重要因素. 本节讨论Rayleigh模型、Tunna模型、Benasciutti-Tovo（BT）模型、Dirlik模型4种应力幅概率分布模型以及IFFT方法对疲劳损伤结果的影响. 阻尼比取值与3.2节相同，得到风、浪单独作用下基础在25 a内的疲劳损伤，并将其与风、浪单独作用的时域结果对比，结果如表3所示.

与风、浪单独作用时的时域法疲劳损伤结果相比，计算风致疲劳损伤时Tunna模型、BT模型、IFFT方法的结果过于保守，Dirlik模型更接近时域结果，误差为13.48%；计算浪致疲劳损伤时Rayleigh模型过于保守，IFFT方法准确性较好，误差为14.03%. 虽然计算风致和浪致疲劳损伤时均可以使用IFFT方法，但本例结果显示基于谐波叠加法的IFFT方法并不适用于风致疲劳损伤（应力谱为典型宽谱）的计算. 叠加分别由Dirlik模型和IFFT方法计算的风致和浪致疲劳损伤，得到频域法的基础总疲劳损伤为5.27×10⁻¹，小于时域法风、浪联合作用的基础疲劳损伤.

（1）基础疲劳损伤对气动阻尼比的取值十分敏感，气动阻尼比对风致动力响应的影响比浪致响应更加显著，气动阻尼比越大，时域法和频域法分别计算得到的风致疲劳损伤的差值越小，在进行疲劳分析时应当合理、准确地选取阻尼比.

（2）时域方法中风、浪单独作用，然后再叠加得到的基础疲劳损伤较风、浪联合作用的结果小43.3%，且这一差异主要由大载荷的工况导致，表明在疲劳计算时考虑风、浪联合作用是十分必要的，简单叠加风致疲劳和浪致疲劳的方法会导致计算结果偏小. 本算例中，风机基础疲劳损伤主要由大风速、大波高且小周期、出现概率大的工况所产生，这对于简化风、浪联合计算的工况数量具有指导意义.

（3）应力幅概率分布模型对频域法计算的疲劳损伤影响显著，本例中风致疲劳损伤由Dirlik模型计算最接近时域结果，浪致疲劳损伤由IFFT方法计算最接近时域结果. 采用频域法计算的风、浪损伤值叠加后得到的总疲劳损伤小于风、浪联合作用时的计算结果.