0 引　言

由于海上风电具有诸多陆上风电不可比拟的优势而受到很多国家重视，但海上风机支撑结构成本高、服务周期长、工作环境恶劣，因此其安全性不容忽视^{[1 –2]}。支撑结构属细长柱体结构，长期承受交变随机载荷，属典型的高周疲劳结构，容易产生疲劳损伤，结构疲劳破坏是其主要失效形式^[3]。且一旦发生疲劳破坏，后果将不堪设想。因此，在设计阶段对支撑结构进行疲劳寿命计算尤为重要。

目前，海上风机支撑结构疲劳寿命计算有采用时域法，即先通过动力响应分析得到危险节点的应力时程曲线、试验统计得到材料的S-N曲线，然后对应力时程曲线进行雨流计数得到应力循环次数，最后结合Palmgren-miner线性累计损伤法计算出结构的疲劳寿命^[4]。但随着装机容量越来越大，工作海域越来越深，支撑结构受到的环境载荷也将十分复杂。用时域法计算其疲劳寿命时需处理的数据庞大，耗时长、效率低^{[5 –6]}。频域法通过动力学仿真和有限元分析计算出支撑结构危险节点的功率谱密度（Power spectral density function，PSD），结合Dirlik概率模型计算出结构的疲劳寿命，省去了时域法中统计应力循环次数这一复杂繁琐过程，可大大缩短计算时间，适用于受长期随机交变载荷作用下的海洋风机支撑结构疲劳寿命计算^[7]。为了探索频域法在海上风机支撑结构疲劳寿命计算中的应用，本文首先建立海上三脚架固定式风机有限元模型，利用FAST软件对 5MW风机叶轮附近空气流场进行仿真，计算得到作用于风机叶片上的气动力载荷。并用Stokes五阶波理论和Morison方程计算出浪流联合载荷。然后分别用时域法和频域法计算出支撑结构危险节点的疲劳寿命，并对 2 种方法的计算结果进行对比分析。

1 疲劳载荷

三脚架式支撑结构受力形式如图 1所示。包括作用于塔柱上的风机载荷、风载荷；作用于斜撑及三脚架桩腿结构上的浪流联合载荷等。其中最主要的载荷是风机载荷和浪流联合载荷^[8]。

1.1 风机载荷

当叶轮转动发电时，会产生作用于塔柱顶部 3 个方向上的力和转动力矩。风机载荷作用形式如图 2所示。按照IEC 64100-3 标准，利用FAST软件建立风场模型，采用a级湍流强度的卡曼湍流模型^[9]。轮毂高度为 90 m，风轮附近空气流场网格面积为 145 m × 145 m，叶片为柔性单元，桨距角自动调节。计算得到风机载荷如图 3所示。

1.2 浪流联合载荷

因为风机支撑结构为细长柱体， $D/L \leqslant 0.2$ ，故采用Stokes五阶波理论和Morison方程计算浪流载荷较为合适^[10]。选取的海域为Atlantic region 2 区域，海洋环境相关参数如表 1所示。通过编程计算桩腿和塔柱上的波浪力分布函数，再沿纵向积分可得桩腿和塔柱浪流联合载荷如图 4所示。

2 海上风机有限元模型及其模态分析

本文采用的风机模型是NREL的 5MW海上风机模型^[11]，基本参数如表 2所示。

风机的支撑结构模型为三脚架式结构模型，它由位于中心的塔柱、斜撑及 3 根插入海床的桩腿组成。三脚架式基础结构简单，建造和施工方便，可减少海底冲刷、波浪流载荷，有良好的强度和稳定性，因此适用于水深较大且海床较为坚固的海域^[12]。其有限元模型如图 5所示。

当风机工作时，支撑结构不仅受到风轮旋转所产生的周期性激励，还会受到风浪流载荷的作用，因此，支撑结构会产生振动，当结构的激励频率和固有频率接近时，结构会产生共振，造成极大的破坏。因此，在进行疲劳分析前，需先对其进行模态分析。 图 6为模态分析的前六阶模态振型，对应的自振频率如表 3所示。

本文使用的风机叶轮为三叶片式，每转 1 周塔柱受到 3 次激励振动，所以共振的激励源为 1 P和 3 P，3 P为作用在塔柱上的轴向推力频率^[10]。塔柱固有频率必须避开激励振动频率。风机模型正常工作时的叶轮转速为 12.1 r/min，则 $1{\rm{P}} = 12.1/60 = 0.2017\;{\rm{Hz}}$ ， $3{\rm{P}} = 0.6051\;{\rm{Hz}}$ ，国标GB/T 19072-2003 中要求取 ± 10% 的裕度，则风机固有频率应避开的频率区间为[0.1815 Hz，0.6656 Hz]，模态分析的固有频率最小值为1阶固有频率 0.7926 Hz，故塔柱工作频率不在共振频率禁区内，理论上结构不会发生共振现象。

3 时域法计算支撑结构疲劳寿命

目前，海上风机支撑结构的疲劳寿命分析多采用时域法，即通过动力响应分析得到的应力时程曲线和试验统计得到的S-N曲线，应用Palmgren-miner线性累计损伤法计算结构的疲劳寿命。

3.1 Palmgren-miner法则

支撑结构在多级循环应力作用下的疲劳损伤是一个累积的过程，每一级应力水平都会产生一组疲劳损伤，从而消耗掉一定分量的疲劳寿命。当结构总的损伤达到某一值时，结构就会发生疲劳破坏。目前最常用的疲劳损伤理论是Palmgren-miner线性累计损伤理论。这一理论认为结构受多级循环应力作用时，各级循环应力造成的损伤可以线性累加^[10]。对于单一循环应力 s _i 作用产生的损伤度为

${D_i} = \frac{{{n_i}}}{{{N_i}}}{\text{，}}$

(1)

式中： n _i 为应力范围 s _i 的实际循环次数； N _i 为结构在应力范围 s _i 作用下产生疲劳破坏所需的循环次数。

若有 k 级循环应力，那么总损伤度为

$D = \sum\limits_{i = 1}^k {{D_i}} = \sum\limits_{i = 1}^k {\frac{{{n_i}}}{{{N_i}}}}{\text{。}} $

(2)

当 D = 1 时，结构便产生疲劳破坏。

3.2 应力时程曲线处理及雨流计数

用Ansys对有限元模型做动力响应分析，可以得到所有单元节点的应力时间历程，其中最大等效应力时程曲线如图 7所示。

从应力时程曲线可看出，应力时程存在一些不会产生疲劳损伤的小应力循环，为了减少处理的数据，提高计算效率，需要对其进行小载荷处理，删除不必要的小应力循环。本文建立一个取舍标准，处理后的应力时程曲线如图 8所示。

因为分析处理得到的载荷是随机不规则的，不能直接用于疲劳寿命的计算，需先对其应力循环进行计数，本文采用的是工程上应用较为广泛的雨流计数法。根据此方法，可以得到应力幅值、应力均值和对应的循环次数，支撑结构最危险节点的应力循环统计图如图 9所示。

3.3 基于时域法的疲劳寿命计算

根据统计处理得到的应力循环统计图和S-N曲线，利用Palmgren-miner线性累计损伤理论，可以求出风机支撑结构最危险节点的疲劳寿命。

4 频域法计算支撑结构疲劳寿命

在使用频域法计算疲劳寿命时，需先对有限元模型加载求得最危险节点的动应力响应，由动应力响应经过傅里叶变换可得到应力的功率谱密度函数 $G\left( f \right)$ ，进而求得峰值概率密度函数 $p(s)$ ，由概率密度函数结合Palmgren-miner线性累计损伤法则，即可进行海上风机支撑结构的疲劳寿命计算^[13]。

4.1 Dirlik概率模型

在频域疲劳分析方法中，得到应力功率谱密度后，可进而求各阶矩 m _n ^[14]。

${m_n} = \int {{f^n}} \cdot G\left( f \right){\rm d}f {\text{，}}$

(3)

式中 $G\left( f \right)$ 为当频率为 f 时的PSD值。

因为风电支撑结构的应力属于宽带随机过程，故可选择Dirlik公式作为评估模型，对支撑结构进行疲劳分析。Dirlik通过运用门特卡罗技术做了大量计算机模拟，得出峰值概率密度函数^[15]。其表达式为：

$p\left( S \right) = \frac{{\frac{{{D_1}}}{Q} \cdot {g^{ - \frac{Z}{Q}}} + \frac{{{D_1} \cdot Z}}{{{R^2}}}{e^{ - \frac{{{Z^2}}}{{2{R^2}}}}} + {D_3} \cdot Z \cdot {e^{ - \frac{{{Z^2}}}{2}}}}}{{2\sqrt {{m_0}} }}{\text{，}}$

(4)

式中：

${D_1} = \displaystyle\frac{{2 \cdot ({x_m} - {r^2})}}{{1 + {r^2}}}\text{；}$

${D_2} = \displaystyle\frac{{1 - r - {D_1} + D_1^2}}{{1 - R}}\text{；}$

${D_3} = 1 - {D_1} - {D_2}$ ；

$Z = \displaystyle\frac{s}{{2\sqrt {{m_0}} }}$ ；

$Q = \displaystyle\frac{{1.25(r - {D_3} - {D_2} \cdot R)}}{{{D_1}}}$ 。

其中， $r = \displaystyle\frac{{{m_2}}}{{\sqrt {{m_0} \cdot {m_4}} }}$ ； ${x_m} = \displaystyle\frac{{{m_1}}}{{{m_0}}}\sqrt {\displaystyle\frac{{{m_2}}}{{{m_4}}}} $ ； m ₂ 和 m ₄ 分别为功率谱密度函数的2阶矩和4阶矩。

结构最危险节点应力循环计算式为

$N{S^m} = C{\text{。}}$

(5)

式中 C 和 m 为材料 ${\rm S} - {\rm N}$ 曲线方程中的参数。

由式（4）和式（5）可推导出宽带随机载荷作用下疲劳寿命的计算公式为：

$\begin{split}\\[-12pt]{N_p} = & C/{\left\{ {(2\sqrt {{m_0}} )} \right.^m}\left[ {{D_1}} \right.{Q^m}{e^{ - \frac{{{S_0}}}{{2Q\sqrt {{m_0}} }}}} \cdot \left( {{{\left( {\frac{{{S_0}}}{{2Q\sqrt {{m_0}} }}} \right)}^m} + } \right.\\& \left. {\sum\limits_{k = 1}^m {\frac{{m!}}{{(m - k)!}}{{\left( {\frac{{{S_0}}}{{2Q\sqrt {{m_0}} }}} \right)}^{m - k}}} } \right) + 2\frac{{m - 2}}{2}\Gamma \left( {\frac{{m + 2}}{2}} \right) \times \\& \left( {{R^m}{D_2}Q_\chi ^2\left( {\frac{{{S_0}}}{{4{R^2}{m_0}}},m + 2} \right) + } \right.\left. {{D_3}Q_\chi ^2\left( {\frac{{{S_0}}}{{4{m_0}}},m + 2} \right)} \right){\text{。}}\end{split}$

(6)

式中， $Q_\chi ^2\left( {\chi _a^2,v} \right)$ 为 ${\chi ^2}$ 分布概率函数； $v$ 为 ${\chi ^2}$ 分布自由度； $\varGamma \left( w \right)$ 为Gamma函数， $\varGamma \left( w \right) = 2\int_0^\infty {{y^{2w - 1}}{e^{ - {y^2}}}} {\rm{d}}y$ 。

4.2 频域法疲劳寿命的计算

由应力时程经傅氏变换可得功率谱密度函数如图 10所示。

将功率谱密度代入式（3）、式（4）和式（6）可得到风机支撑结构最危险节点的疲劳寿命。

5 计算结果对比分析

采用雨流循环计数的时域法和Dirlik载荷概率密度函数的频域法分别计算海上风机支撑结构的疲劳寿命，计算结果如表 4所示。

从计算结果可看出，采用Dirlik载荷概率密度函数的频域法计算结果比较准确，与时域法计算结果相对误差在可接受范围内。因为频域法省去了时域法中应力循环计数这一过程，故相对于时域法，频域法在计算结构疲劳寿命过程中具有耗时短的优势，对于受长期交变载荷的复杂结构，这一优势更加明显。但需要注意的是，频域法适用范围较窄，目前仅适用于线性系统，对于非线性系统，需采用时域法进行疲劳寿命的计算。

6 结　语

通过对海上风机支撑结构动力响应分析得到危险节点应力时程曲线，进而求出功率谱密度函数，结合Dirlik概率模型，计算出支撑结构的疲劳寿命，并与时域法计算结果作对比分析。得到基本结论如下：

1） 采用Dirlik概率模型的频域法能够较为准确地计算出复杂交变载荷载荷作用下海上风机支撑结构的疲劳寿命，频域法和时域法计算结果相对误差在可接受范围内。

2）采用Dirlik载荷概率密度函数的频域法替代雨流计数的时域法，可省去统计应力循环次数这一繁琐过程，缩短计算时间，在受长期随机交变载荷作用下的结构疲劳寿命的计算过程中尤为明显。