0 引　言

为了保证海上风机的结构完整性，国际电工委员会在设计标准IEC 61400-3^[1]中，对支撑结构和叶轮-机舱组合件的设计提出了极限强度分析的要求，其中设计载荷 ${F_d}$ 是分析中的必要参数，求解方法通常为IEC 61400-1^[2]中介绍的基本方法，直接积分法。虽然直接积分法便于理解与扩展，但是由于算法本身的缺陷，当考虑的环境变量中不止有风还有波浪等时，所需的仿真次数也急剧增加，因此给海上风机极限载荷的求解造成了困难。

为了减少仿真数量，反向一阶可靠度法IFORM（Inverse First-Order Reliability Method）逐渐被运用到考虑多维变量时极限载荷的求解上。P.J. Moriarty等^[3]曾采用直接积分法求解WP_Baselline1.5MW风机叶片根部面外弯矩在1年重现周期下的极端载荷，而Korn Saranyasoontorn等^[4]采用2维IFORM法（环境等值线法）以更少的仿真数据便得到了同样良好的结果。K. Saranyasoontorn等^[5]曾根据IFORM的3种形式建立了求解风机名义载荷的基本模型，求解了某600 kW陆上风机叶片的极端载荷；D. Karmakar等ADDIN EN.CITE.DATA ^{[6 – 8]}利用IFORM法求解了3种不同浮式基础NREL5MW海上风机的极限载荷，并对计算结果与浮式基础的关系进行比较；Puneet Agarwal等^[9]在样本经验分布函数的基础上采用环境等值线法求取了工作在20 m浅水区域NREL5MW海上风机的长期载荷。

但是在上述研究中，一般都是等间距地在搜索区域上选取搜寻点，外推结果的精度由间隔搜寻角的大小直接决定，间隔角越小外推结果就越精准，但同时仿真数量也会增加。而采用IFORM求解极端载荷实际上也是一个最优化过程，因此本文将环境等值线法与1维最优化方法二分法结合起来，形成一种具有搜寻策略的算法来求解极端载荷，通过多级搜寻逐步缩小搜寻区间，减少不必要的搜寻工作，提高外推精度。同时与传统的直接积分法进行对比验证。

1 计算原理 1.1 IFORM原理简介

直接积分法的计算公式如下

$\begin{aligned}{P_T} = {Pr} \left\{ {L > {l_T}} \right\} = \int\nolimits_{{X_1}} {\int\nolimits_{{X_2}} { \cdots \int\nolimits_{{X_n}} {Pr \{ {L > {l_T}\left| {{X_1}} \right. ={x_1} ,{X_2} =}}}}\\{{{{ {x_2}, \ldots ,{X_n} \!=\! {x_n}} \}{f_{{X_1}{X_2}, \ldots ,{X_n}}}\left( {{x_1},{x_2}, \ldots ,{x_n}} \right){\rm d}{x_1}{\rm d}{x_2}, \ldots ,{\rm d}{x_n}} } } \text{。} \end{aligned}$

(1)

其中： $T$ 为重现周期； $L$ 为短期仿真中的载荷极值； ${l_T}$ 和 ${P_T}$ 为对应的极限载荷和超越概率； ${X_1}, {X_2},\ldots ,{X_n}$ 为环境变量，如平均风速、有义波高、谱峰周期等； ${f_{{X_1}{X_2}，...，{X_n}}}$ 为环境变量分布的联合概率密度函数。

从反向可靠性的角度分析积分法求解极限载荷 ${l_T}$ 的过程。现假设 ${l_T}$ 事先已知并且作为结构设计抗力 ${l_{dec}}$ ，则事件 $L > {l_{dec}}$ 为失效事件，概率 ${Pr} \left\{ {L > {l_{dec}}} \right\}$ 为失效概率，记为 ${P_f}$ ，式(1)左半部分变为如下

${P_f} = {Pr} \left\{ {L > {l_{dec}}} \right\}\text{。}$

(2)

为了简化问题，忽略给定环境下短期载荷极值 $L$ 的不确定性，将 $L$ 当作环境变量 $X$ 的确定性函数。再将相互依赖的环境变量 $X$ 转换为独立且服从正态分布的标准正态变量 $U$ ，则短期载荷极值 $L$ 也成为标准正态变量 $U$ 的确定性函数

$L = l\left( U \right) \text{。} $

(3)

在“正向”的一阶可靠度问题中，设计抗力 ${l_{dec}}$ 事先已知，要求的可靠性指标 $\beta $ 由以下最优化问题得到^[10]。

$\begin{array}{l}{\text{优化变量}}:{\text{标准正态变量}}U\text{，}\\{\text{约束条件}}:g\left( U \right) = {l_{dec}} - l\left( U \right) = 0\text{，}\\{\text{目标函数}}:\beta = \min \left\{ {\left| U \right|} \right\}\text{。}\end{array}$

而在“反向”的一阶可靠度问题中，事先已知的是可靠性指标 $\beta = {\Phi ^{ - 1}}\left( {1 - {P_f}} \right)$ ，要求的变量是设计抗力 ${l_{dec}}$ ，问题转化为如下形式

$\begin{array}{l}{\text{优化变量}}:{\text{标准正态变量}}U\text{，}\\{\text{约束条件}}:\left| U \right| = \beta\text{，} \\{\text{目标函数}}:{l_{dec}} = \max \left\{ {l\left( U \right)} \right\}\text{。}\end{array}$

对比可知，在对结构进行可靠性分析中，搜寻的目标，即可靠性指标 $\beta $ 为标准正态空间内失效面上与坐标原点距离最近的点。而在求解设计抗力 ${l_{dec}}$ 的问题中，搜寻的目标为标准正态空间内半径为 $\beta $ 的球面上短期载荷极值 $L = l\left( U \right)$ 最大的点。

1.2 IFORM的3种形式

对于海上风机，通常主要考虑环境变量为风速 $V$ 和有义波高 $Hs$ 。记 $V$ 的概率密度函数和分布函数分别为 ${f_V}\left( v \right)$ 和 ${F_V}\left( v \right)$ ；给定 $V$ 下 $Hs$ 的密度函数和分布函数分别为 ${f_{Hs\left| V \right.}}\left( {hs\left| v \right.} \right)$ 和 ${f_{Hs\left| V \right.}}\left( {hs\left| v \right.} \right)$ ；给定 $V$ ， $Hs$ 下短期载荷极值 $L$ 的分布函数为 ${F_{L\left| {V，Hs} \right.}}\left( {l\left| {v，hs} \right.} \right)$ 。则式（2）可改写为

$\begin{aligned}{P_f}\!\!= \! \!\! &\Pr \left\{ {L \!>\! {l_{dec}}} \right\} \!\!=\!\! \int\nolimits_{{V_{in}}}^{{V_{out}}} {\int\nolimits_0^\infty \!\!{{Pr} \left\{ {L\! >\! {l_{dec}}\left| {V \!\!=\!\! v\text{，}\!\!\!Hs \!\!= \!\!hs} \right.} \right\}}} \text{，} \\ &{{{f_{Hs\left| V \right.}}\left( {hs\left| v \right.} \right){f_V}\left( v \right){\rm d}v{\rm d}hs} } \text{。} \end{aligned}$

(4)

其中 ${V_{in}}$ 和 ${V_{out}}$ 分别为切入风速和切出风速。采用Rosenblatt变换^[11]将物理变量 $V$ ， $Hs$ ， $L$ 转换为标准正态变量 ${U_1}$ ， ${U_2}$ ， ${U_3}$ 。

$\varPhi \left( {{u_1}} \right) = {F_V}\left( v \right) \text{，} $

(5)

$\varPhi \left( {{u_2}} \right) = {F_{Hs\left| V \right.}}\left( {hs\left| v \right.} \right) \text{，} $

(6)

$\varPhi \left( {{u_3}} \right) = {F_{L\left| {V，Hs} \right.}}\left( {l\left| {v，hs} \right.} \right) \text{，} $

(7)

$u_1^2 + u_2^2 + u_3^2 = {\beta ^2} \text{。} $

(8)

由此便可在标准正态空间内搜寻引起短期载荷极值 $L$ 最大的风浪组合点。根据简化程度的不同，通常分为3种情况，如图 1所示。

图 1（c）为IFORM法最基本的形式，即三维IFORM法。图 1（b）为IFORM法的二维简化形式。设 ${l_{0.5}}$ 为短期载荷极值分布 ${F_{L\left| {V，Hs} \right.}}\left( {l\left| {v，hs} \right.} \right)$ 的中位数，满足条件

$\Pr \left\{ {L \geqslant {F_{L\left| {V，Hs} \right.}}\left( {{l_{0.5}}\left| {v，hs} \right.} \right)} \right\} = 0.5\text{。}$

(9)

则当式（7）等号右边 $l$ 取 ${l_{0.5}}$ 时，有 ${u_3} \equiv 0$ 。此时标准正态空间内的三维球面简化为 $O{U_1}{U_2}$ 平面上的二维圆，圆上 ${l_{0.5}}$ 最大的点便成为要求的设计点。由于该圆是 $\beta $ 的等值线，所以其在物理空间环境变量平面 $OVHs$ 上对应的曲线称为环境等值线（Environmental Contour），该方法也称为环境等值线法或EC法。

需要注意的是，在EC法的搜寻过程中，只比较中位数 ${l_{0.5}}$ 的大小而忽略短期载荷极值 $L$ 的不确定性必然有所误差。所以需要采取额外的修正措施来弥补EC法的上述不足，例如采用更高的分位数（大于0.5）或者更大的重现周期等。本文采取如下修正方法。

将 $L$ 的不确定性分为2部分：

$L = {L^*} \cdot \varepsilon \text{，} $

(10)

其中： ${L^*}$ 为由风浪条件变化引起的不确定性； $\varepsilon $ 为同一风浪条件下短期载荷极值本身的不确定性，即仿真中随机数的影响。这2个变量都可以用对数正态分布来拟合，它们取自然对数后的标准差为

${\sigma _{\ln {L^*}}} = \frac{{\ln \left( {L_{T_1}^*/L_{T_2}^*} \right)}}{{{\beta _{T_1}} - {\beta _{T_2}}}} \text{，} $

(11)

${\sigma _{\ln \varepsilon }} = \frac{{{\varepsilon _{p2}} - {\varepsilon _{p1}}}}{{{\varPhi ^{ - 1}}\left( {p_2} \right) - {\varPhi ^{ - 1}}\left( {p_1} \right)}} \text{。} $

(12)

其中： $T_1$ 为要求的重现周期； $T_2$ 为比 $T_1$ 稍短的重现周期； ${\beta _{T1}}$ 和 ${\beta _{T2}}$ 为对应的可靠性指标； $p_1$ 为初始的分位数水平0.5； $p_2$ 是比 $p_1$ 稍高的分位数水平。在得到 ${\sigma _{\ln {L^*}}}$ 和 ${\sigma _{\ln \varepsilon }} $ 后由下式计算修正因子 ${R_T}$

${R_T} = \frac{{{L_{T_1}}}}{{L_{T_1}^*}} = \exp \left[ {\left( {{\sigma _{\ln L}} - {\sigma _{\ln {L^*}}}} \right){\beta _{T_1}}} \right] \text{，} $

(13)

其中 $\sigma _{\ln L}^2 = \sigma _{\ln {L^*}}^2 + \sigma _{\ln \varepsilon }^2$ 。

如果在式（10）中取 $hs$ 为 $h{s_{0.5}}$ ，则 ${u_2}$ 也恒等于0，此时三维IFORM便简化为一维，唯一的搜寻点即当作设计点。

2 数值算例 2.1 Spar型海上浮式风机

NREL海上5MW风机是美国国家可再生能源实验室（National Renewable Energy Laboratory）开发的一款风机^[12]，根据子结构和基础形式的不同分为多种类型。其中，OC3Hywind是一种Spar型浮式基础的海上风机，如图 2所示。基本参数见表 1。

2.2 计算工况

计算工况以IEC 61400-3中正常发电工况的DLC 1.1为例，风机等级IECⅠ-A。参考Korn Saranyasoontorn等^[13]中丹麦Horns Rev处的海洋气象条件，调整相应的参数使风机运行在一个较合理的环境下。假定平均风速服从韦伯分布

${F_V}\left( v \right) = 1 - \exp \left[ { - {{\left( {\frac{v}{a}} \right)}^k}} \right],\;\;\;a = 12.258\;3,k = 1.8 \text{，} $

(14)

并且只考虑在工作风速区间内失效的情况，因此对切入风速以下和切出风速以上进行截断，将平均风速的实际分布函数改写为

${F_V}\left( v \right) = \frac{{G\left( {{v_{in}}} \right) - G\left( v \right)}}{{G\left( {{v_{in}}} \right) - G\left( {{v_{out}}} \right)}} \text{，} $

(15)

$G\left( v \right) = \exp \left[ { - {{\left( {\frac{v}{a}} \right)}^k}} \right] \text{。} $

(16)

假定有义波高服从正态分布。

$\begin{aligned}& {F_{Hs\left| V \right.}}\left( {hs\left| v \right.} \right) \!=\! \Phi \left[ {\frac{{hs - {\mu _{Hs}}\left( v \right)}}{{{\sigma _{Hs}}}}} \right],\\\;\;\;& {\mu _{Hs}}\left( v \right) \!=\! 0.13v,{\sigma _{Hs}} \!=\! 0.24\text{。}\end{aligned}$

(17)

为了简化问题，将谱峰周期 ${T_p}$ 作为有义波高的确定性函数。参考IEC61400-3中确定性设计波浪对应周期范围的公式

$11.1\sqrt {{H_{S,NSS}}\left( V \right)/g} \leqslant T \leqslant 14.3\sqrt {{H_{S,NSS}}\left( V \right)/g} \text{，} $

(18)

对谱峰周期取值如下：

${T_p} = \frac{{\left( {11.1 + 14.3} \right)}}{2}\sqrt {\frac{{Hs\left( V \right)}}{g}}\text{。} $

(19)

式中 $g$ 为重力加速度，取9.8 m/s。

3 极限载荷求解 3.1 搜索算法

采用环境等值线法求解该风机1年重现周期和20年重现周期下的长期载荷，求解对象为叶片根部面外弯矩OOPB（Out-of Plane Blade Root Moment）以及塔筒基底首尾弯矩TBM（Fore-aft Tower Base Moment）。基本流程如图 3所示。

3.2 1年重现周期结果

设定收敛条件 $\delta $ 等于3°。1年重现周期下OOPB的极端载荷见表 2和表 3。在进行Rosenblatt变换时，环境等值线上风速较低的2个搜寻点出现了有义波高为负的情况，这是有义波高分布函数本身的定义造成的。而最大短期载荷极值多发生在海况较恶劣的风浪点，因此对设计点的搜寻不造成影响。经过试算，发现同等风速下有义波高较小时中位数 ${l_{0.5}}$ 也较小，因此相同风速的海况中有义波高较小的一组没有进行仿真试验。

最终获得1年重现周期下OOPB极限载荷 $14\;703 \times$ $1.027\;6 = 15\;109\;{\rm{kN}} \cdot {\rm{m}}$ 。标准正态空间中的OOPB等值线与物理空间中的OOPB等值线如图 4所示。图中标注的数字对应表 1中的搜寻顺序。

由于计算中发现TBM短期极值的不确定性比OOPB更大，因此设定收敛条件为 $\delta $ =2°。1年重现周期下TBM的极端载荷见表 4和表 5。

最终获得1年重现周期下的TBM极端载荷 $145\;129 \times$ $1.048\;8 = 152\;211\;{\rm{kN}} \cdot {\rm{m}}$ 。标准正态空间中的TBM等值线与物理空间中的TBM等值线如图 5所示。

3.3 20年重现周期结果

为了之后与直接积分法进行更多的对比验证，同时还计算了20年重现周期时的极限载荷。搜寻过程与1年重现周期类似，在此不予赘述。最终获得20年重现周期下OOPB极限载荷为 $15\;167 \times 1.048\;9 = $ $15\;909\;{\rm{kN}} \cdot {\rm{m}}$ 。20年重现周期下TBM极限载荷为 $145\;472 \times 1.119\;7 = 162\;885\;{\rm{kN}} \cdot {\rm{m}}$ 。

4 积分法对计算结果的验证

为了对比验证搜索算法的外推结果，根据式1采用直接积分法求解OOPB和TBM的长期分布。首先确定整个2重积分的积分域。根据前文定义的计算工况和平均风速的分布函数，取平均风速的积分区间为5 m/s~25 m/s；根据有义波高的分布函数取积分区间0 m~5 m。然后将各自的积分区间划分成多个子区间。为此将平均风速按照2 m/s的分辨率划分成11个子区间，各子区间平均风速为5 m/s，7 m/s，···，25 m/s。将有义波高按照1 m的分辨率划分成5个子区间，各子区间平均有义波高分别为0.5 m，1.5 m，···，4.5 m。综上所述，总积分域被划分成55个子区间。

确定好各子区间后根据其风浪参数，对风机在该环境条件下的运行进行数值仿真，以获取短期极值点 $L$ 。本文对每个风浪点均进行至少6次随机仿真。首先由湍流风软件TurbSim^[14]产生湍流风数据，然后导入风机仿真软件FAST^[15]进行10 min短期仿真，最后从输出结果中提取时域响应的极值点。采用常见的分块法选取极值点，根据夏一青等^[16]的建议取分块数30，即从每次仿真中选取30个极值点。以平均风速13 m/s，有义波高4.5 m的某次随机仿真为例，选取OOPB短期极值点如图 6所示。

在获得短期极值点后，需要对其分布进行拟合，计算 ${F_{L\left| {V{\text{，}}Hs} \right.}}\left( {l\left| {v{\text{，}}hs} \right.} \right)$ 。本文采用常见的3参数韦伯分布进行拟合，以线性矩法估计分布参数，采用Q-Q图对拟合分布进行检验。同样以上述风浪点为例，拟合分布的概率密度函数和累积分布函数如图 7和图 8所示。拟合分布与样本经验分布的分位数对比如图 9所示。

注意以上求解 ${F_{L\left| {V，Hs} \right.}}\left( {l\left| {v，hs} \right.} \right)$ 的过程在搜索算法计算中位数 ${l_{0.5}}$ 时也会用到。而在积分法中，在式（1）右侧取不同的 ${l_T}$ 可得到不同的超越概率 $Pr \left\{ {L > {l_T}} \right\}$ 。由此计算出OOPB和TBM的长期超越概率曲线如图 10和图 11所示。在1年的时长内，有365×24×6=52 560个10 min，因此1年重现周期对应的超越概率为1/52 560=1.902 6E-5，相应地20年重现周期对应的超越概率为9.512 9E-7。

1年和20年重现周期下搜索算法与积分法的对比如表 6所示。其中OOPB1和TBM1为搜索算法直接搜寻的结果，OOPB2和TBM2为经修正点修正之后的结果。

由表中数据可见，在4项对比中，搜索算法预估的失效点与积分法预估的设计点均较一致，平均风速的误差不超过1 m/s，这表明搜索算法具有良好的识别能力，能够有效地分析出风机正常发电时最有可能出现失效的情况。

对于不确定性较小的叶片根部首尾弯矩OOPB，搜索算法外推结果与积分法更接近，修正因子 ${R_T}$ 的数值也更小。塔筒基底首尾弯矩TBM较OOPB的不确定性更大，在和积分法直接比较时，存在一定误差，并且重现周期越大误差也越明显。然而较大的不确定性也导致了较大的修正因子 ${R_T}$ ，所以在一定程度上弥补了搜索算法的不足。经过修正因子修正之后，外推结果与积分法更为接近，误差在5%左右。

5 结　语

将EC法与二分法结合起来形成一种新的搜索算法，应用于求解Spar型浮式风机叶片根部面外弯矩和塔筒基底首尾弯矩在1年和20年重现周期下的极端载荷，主要可得以下结论：

1）相较于积分法55组的仿真数量，搜索算法仅利用17~19组仿真便获得了相应的极限载荷与设计点，因此在计算成本上具有较大的优势。

2）经过与传统直接积分法的外推结果对比验证，发现在极限载荷数值和引发失效的设计点上均较一致，说明了搜索算法同时也具有较高的可靠性。

3）作为参考，EC法还可结合其他最优化方法以进一步减少仿真需要的数量，提高求解精度。