0 引　言

在海洋强国的时代背景下，国家大力引导发展海洋旅游业的经济，作业于近远海的超大型浮式综合旅游体应运而生。而海上不比陆地，其情况复杂多变，结构要承受的环境载荷更加突出，浮体长期主要承受波浪造成的交变载荷，浮体结构极易因疲劳而产生破坏，因此对大型浮体结构疲劳特性的研究具有十分广阔的应用前景 ^[1]。疲劳估价法主要包括时域疲劳估价法和频域疲劳估价法。在实际项目中，主要根据所获得结构应力的形式和难易程度来选择疲劳评估方法。当结构应力的时间过程较容易得到时，通常采用时域疲劳估价法 ^{[1 − 4]}。当结构受力的功率谱密度函数容易算出时，通常采用频域疲劳估价法^{[5 − 7]}。对于一些需要考虑复杂非线性因素的结构，在疲劳分析过程中可以用时域疲劳估价法进行计算，分析精度很高，但计算量却很大。频域疲劳估价法基于线性化假定，忽略非线性因素以及平均应力的影响，可执行程度高，计算量小。

目前，对于船舶这样的浮式结构来说，其应力反应的时间过程很难获得，因此对其疲劳损伤的评估通常采用频域谱分析方法^{[8 − 9]}。张朝阳等 ^[8]运用热应力谱分析方法，计算得出深水半潜式平台薄弱节点处的疲劳寿命，对不同浪向及各短程海况对节点总疲劳损伤度的贡献进行了进一步分析。马网扣等 ^[9]运用程序对某新型深水半潜式钻井平台典型节点疲劳寿命进行分析。谢文会等^[10]采用疲劳谱分析方法，基于Miner准则，计算了薄弱节点的热点应力传递函数，并对平台典型节点的疲劳寿命进行评估。

常规的海洋工程结构物的疲劳研究，行业内有着丰富的经验，常规海洋工程结构物的疲劳节点位置都位于主体受力框架结构上^[11]。而随着浮体结构形式的改变，结构受力分析更加困难，如何准确分析结构的疲劳成为一大难点。在对结构疲劳分析中，疲劳位置传递函数的计算对疲劳分析的准确度起着至关重要的作业。一般有4种传递函数的计算方法：1）理论评估法；2）实部虚部合成法；3）简化理论分析法；4）规则波数值模拟试验法。其中，规则波数值模拟试验法的计算原理最接近真实海况。本文以某海上大尺度异型镂空浮体为对象，采用规则波数值模拟试验法计算了浮体疲劳点位的传递函数，采用2种谱分析方法，即名义应力和热点应力谱分析方法对异形体结构进行计算，对比了2种计算结果的差异性，给出了合理的工程设计建议。

1 结构疲劳分析的基本原理 1.1 谱分析法的基本原理

频谱分析法是研究船舶结构荷载和响应的一种常用方法，其理论是将随机理论中的线性系统做变换^[11]。而船舶结构又是动力系统的一种，系统的输入部分$\eta (t)$为作用在结构上的波浪过程，系统的响应部分$X(t)$为结构受到波浪荷载产生的交变应力。系统的响应过程与输入过程之间存在以下关系^[12]：

$ X(t) = L\left[ {\eta (t)} \right]。$

(1)

式中：$L$为将$\eta (t)$转化为$X(t)$的算子。当$L$为线性算子时，整个系统即为线性系统。

根据线性理论，在线性系统条件下，便可以计算船舶结构在波浪载荷下的疲劳寿命，不妨假设波浪是一个平稳随机过程，则转化后的响应亦是一个平稳的随机过程。由随机理论，2个平稳随机过程的功率谱密度函数关系如下^[12]：

$ {G_X}(\omega ) = {\left| {H(\omega )} \right|^2}\cdot {S_\eta }(\omega )。$

(2)

式中：$\left| {H(\omega )} \right|$为传递函数；${\left| {H(\omega )} \right|^2}$为响应幅值算子RAO；$ {G_X}\left( \omega \right) $为应力的响应谱。

1.2 传递函数计算方法

传递函数是谱分析法中最关键的要素。传递函数的物理意义是在线性动力系统中做简谐振动，其圆频率为$\omega $时，响应过程的波幅比上输入过程的波幅。传递函数是指以波浪运动为输入过程，交变应力为响应过程，结构在圆频率为$\omega $的规则余弦波作用下应力幅值比上波幅。本文采取规则波数值模拟试验法进行结构传递函数的计算，具体关系如图1所示。利用频率响应法，对每组波的不同频率，分别测定各组波的波高及其相应系统的反应，进行一系列的规则波实验，即可求得传递函数。

1.3 波浪谱

本项目浮体作业位置离岸较近，水深较浅，可采用改进的JONSWAP谱进行计算^[13]。其表达形式如下：

$ S(f)=\beta_{J} H_{\frac{1}{3}}^{2} T_{P}^{-4} f^{-5} \exp \left[-\frac{5}{4}\left(T_{P} f\right)^{-4}\right] \gamma^{{{\mathrm{exp}}}\left[-\left(\frac{f}{f_{F}}-1\right)^{2} / 2 \sigma^{x}\right]}。$

(3)

式中：

$ {{B_J} =\dfrac{{0.06238}}{{0.23 + 0.0336\gamma - 0.185{{(1.9 + \gamma )}^{ - 1}}}}[1.094 - 0.01915\ln \gamma ]。}$

$\gamma $为峰高因子，本项目中取峰高因子的平均值3.3；$\sigma $为峰形参数，频率在峰值点左侧时取0.07，频率在峰值点右侧时取0.09给出有效波高与谱峰周期即可确定JONSWAP谱曲线。

根据公式$ S(\omega)\mathrm{d}\omega=S(f)\mathrm{d}f，\omega=2\text{π}f $可将JONSWAP谱换算为自变量为圆频率的形式：

$ \begin{aligned}S(\omega)= & \frac{1}{2\text{π}}\beta_JH_{\frac{1}{3}}^2T_P^{-4}\left(\frac{\omega}{2\text{π}}\right)^{-5}\times \\ & \exp\Bigg[-\frac{5}{4}\left(\frac{T_P\omega}{2\text{π}}\right)^{-4}\Bigg]\gamma^{\exp[-(\frac{\omega}{2\text{π}f_P}-1)^2\mathord{\left/\vphantom{-(\frac{\omega}{2\text{π}f_P}-1)^22\sigma^2}\right.}2\sigma^2]}。\end{aligned} $

(4)

再由式（2）即可得到疲劳位置的应力响应谱。

1.4 异形浮体疲劳累计损伤度的计算

在获得结构应力的响应谱后, 即可计算出：

$ {G_{XX}}\left( {{\omega _{\text{e}}}} \right) = \int_{ - \text{π} /2}^{\text{π} /2} {f\left( \beta \right)} {\left[ {H\left( {{\omega _e},\theta - \beta } \right)} \right]^2} \cdot {G_{\text{x}}}\left( \omega \right){\mathrm{d}}\beta。$

(5)

式中：$ f\left( \beta \right) $为波浪扩散函数；$ \theta $为主浪向与航向之间的夹角。

波浪扩散函数表达式为：

$ f\left( \beta \right) = k{\cos ^n} \beta 。$

(6)

式中：k和n的取值要满足使$ \int_{-\text{π}/2}^{\text{π}/2}f\left(\beta\right)\mathrm{d}\beta=1 $的条件。

本文在计算中根据的ITTC推荐，取$ k = 2/\text{π} ，n = 2 $，则具体的扩散函数表达如下：

$ f\left( \beta \right) = \frac{2}{\text{π} }{\cos ^2} \beta 。$

(7)

交变载荷过程的跨零率${\nu _0}$将用于疲劳计算，即单位时间内以正斜率跨越零均值的平均次数，以便于求得给定时间内的荷载循环次数，其表达式如下：

$ {\nu _0}{\text{ = }}\frac{1}{{2\text{π} }}\sqrt {\frac{{{m_2}}}{{{m_0}}}}。$

(8)

式中：${m_0}、{m_2}$分别为功率密度谱的0阶矩和2阶矩，可按下式计算。

$ {m_n} = \int\limits_0^{ + \infty } {\omega _e^n} \cdot {G_{XX}}({\omega _e}){\mathrm{d}}{\omega _e}。$

(9)

在实际使用时，将窄带平稳随机过程认定为某海况给定航向下的交变应力过程，其均值为0，由随机理论可知，结构的应力峰值服从Rayleigh分布，其概率密度函数为^[12]：

$ f(x) = \frac{x}{{{\sigma ^2}}}\exp \left( { - \frac{{{x^2}}}{{2{\sigma ^2}}}} \right)\begin{array}{*{20}{c}} {}&{0 \leqslant x < + \infty } 。\end{array} $

(10)

式中：$x$为结构的应力峰值，$\sigma $为在交变应力下的标准差。则假设应力范围S和应力峰值$x$之间有以下关系：

$ S = 2x。$

(11)

最终，得到其在应力范围下的概率密度函数：

$ {f_S}(S) = \frac{S}{{4{\sigma ^2}}}\exp \left( { - \frac{{{S^2}}}{{8{\sigma ^2}}}} \right)\begin{array}{*{20}{c}}，\end{array}0 \leqslant S < + \infty。$

(12)

假设船舶航行时间为${T_{ij}}$，下标表示为第$i$海况、第$j$航向，并用${D_{ij}}$表示在该期间的累积损伤度，则有

$ D_{ij}=\frac{T_{ij}.\nu_{0ij}}{A}\int_0^{+\infty}S^mf_{Sij}(S)\mathrm{d}S。$

(13)

式中，${\nu _{0ij}}$为该交变应力过程的跨零率。

应力范围S和循环次数N的关系可表示为：

$ N \cdot {S^m} = A 。$

(14)

式中：A、m均为S-N曲线的参数。具体可参照《海洋工程结构物疲劳强度评估技术指南》，根据疲劳点位置的不同选取对应的参数^[14]。

将相应的短期应力范围分布${f_{Sij}}(S)$代入式（13）得：

$ \begin{aligned} {D_{ij}} = &\frac{{{T_{ij}}.{\nu _{0ij}}}}{A}\cdot{(2\sqrt 2 {\sigma _{Xij}})^m}\cdot \Gamma \Bigg(1 + \frac{m}{2}\Bigg) = \\ &\frac{{{T_{ij}}.{\nu _{0ij}}}}{A}\cdot {(2\sqrt {2{m_{0ij}}} )^m}\cdot \Gamma \Bigg(1 + \frac{m}{2}\Bigg)。\end{aligned} $

(15)

对于服役时间为$T$，有$m$个应力长期分布的海况组成，并且有出现概率是${p_i}$的各海况；划分$n$个对应的航向角，各航向出现的概率为${p_j}$，则

$ {T_{ij}} = T\cdot{p_i}\cdot{p_j}。$

(16)

该期间，总损伤应为

$\begin{aligned} D = &\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {{D_{ij}}} } = \frac{T}{A}\cdot \Gamma \Bigg (1 + \frac{m}{2}\Bigg )\cdot \\ &\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {{p_i}.{p_j}.} } {\nu _{0ij}}.{(2\sqrt {2{m_{0ij}}} )^m}。\end{aligned} $

(17)

2 异形浮体疲劳位置的选取 2.1 异形浮体模型的建立

本文采用通用的建模软件Ansys-APDL建立异形浮体的有限元模型。其中，板材和T型材采用shell单元，扶墙材和角钢采均用beam单元模拟建立。浮体的有限元模型如图2所示。

此异形浮体整体呈“心”型，中间有4个镂空部分，结构设计新颖，尺寸偏大（400 m×400 m），该结构尺寸超过了一般船用规范计算标准，属于超大型海上综合式浮体，并且宽度方向上存在2个宽度较窄的长条型连接结构，此部分结构会受到巨大的扭矩和弯矩，强度偏弱，在复杂多变的海洋环境中作业极易产生疲劳破坏。

2.2 异形浮体的水动力分析

从浮体有限元模型提取湿表面生成水动力计算模型，将湿表面模型导入到专业从事水动力计算的软件Aqwa中进行水动力计算的前处理，结构的重心坐标及转动惯量均从模型的静力分析计算得出。图3所示为浮体湿表面模型，共包括26542个单元。将处理好的浮体质量模型导入水动力计算软件Aqwa中进行求解，求解完成后即可得到浮体任意剖面内的剪力和弯矩，即为传递函数。传递函数是单位波高产生的浮体响应函数，与浮体外形型、遭遇频率、浪向角等紧密相关。图4为计算工况下浮体中部垂向弯矩的传递函数。

2.3 疲劳位置的确定

根据水动力计算结果，从Aqwa中提取0°浪向 0.07767 Hz、30°浪向0.06 588 Hz、60°浪向 0.04818 Hz、90°浪向 0.08357 Hz、120°浪向 0.07178 Hz、150°浪向 0.06588 Hz、180°浪向 0.6588 Hz的波浪载荷文件。将波浪载荷加载到不含静水压的浮式综合体有限元模型上，并计算得到结果。

以0°浪向 0.07767 Hz浮体应力云图为例，如图5～图6所示。浮体结构分段如图7所示。

经过计算，从各个浪向的应力响应计算结果中选出7个应力较大的位置作为疲劳分析的关注点，7个点的单元编号及位置如表1所示。

根据《海洋工程结构物疲劳强度评估技术指南》（2022）^[14]规范所述，疲劳点位置的有限元模型网格需要细化。以2号疲劳点位置为例，网格细化见图8。

3 传递函数的数值模拟

以0°浪向计算疲劳点的应力响应为例说明规则波数值模拟计算方法。

利用Aqwa软件，计算得到0°浪向角圆频率为0.001 Hz的单位波幅按一定的初相位间隔排列的一系列规则波的浮体运动以及外部水动压力的响应，相位取为0°～360°，间隔45°，每个浪向频率下8个规则波，提取波浪载荷文件。分别加载到结构的有限元模型上，得到不同的应力响应，再分别采取名义应力法和热点应力法计算疲劳位置的应力。最后通过正弦拟合得到疲劳点位在同一频率下所有相位的应力响应。

1）名义应力法：名义应力可以由提取疲劳点计算得到的应力乘以应力集中系数SCF（stress concentration factor）得到，针对不同的位置有不同的系数计算公式，具体的应力集中系数可查阅CCS《海洋工程结构物疲劳强度评估技术指南》3.2～3.4节得到^[14]。采用名义应力法得到的响应曲线，如图9所示。

2）热点应力法：焊缝处的热点应力值是通过插值的方法求出的，具体方法如图10所示。采用热点应力法得到的响应曲线，如图11所示。求焊缝1号点的热点应力，不能直接读取1号点的应力，需要分别计算2条路径（由3到2到1或者由5到4到1）的线性插值，取其中较大值为1号点的热点应力。具体可由下式求得：

$ {\sigma _1} = \max \left\{ \begin{gathered} 1.5{\sigma _2} - 0.5{\sigma _3}，\\ 1.5{\sigma _4} - 0.5{\sigma _5}。\\ \end{gathered} \right. $

(18)

3.1 名义应力法计算的传递函数

分别计算每个浪向所对应的50个频率点，并提取前文所述的7个疲劳关注点处的应力值计算出名义应力，最终标定的传递函数结果如图12～图15所示，此处这里仅给出4个点的传递函数。

3.2 热点应力法计算的传递函数

采用第3节所述热点应力法线性外插的方法计算出各个疲劳关注点的热点应力值，得到异形浮体结构的传递函数如图16～图19所示。

4 浮体疲劳疲劳寿命的计算

在得到疲劳点的传递函数后，选取式（4）给出的JONSWAP谱，按照式（2）即可求出应力的响应谱，将响应谱按照式（5）进行积分可得到结构交变载荷响应的功率谱密度，为了获得给定时间内的应力循环次数，要用到交变载荷过程的跨零率${\nu _0}$，可按式（8）计算得到。本项目浮体设计寿命为25年，最后根据前述疲劳损伤度的计算原理中式（12）～式（17）即可求得各个疲劳点的疲劳累积损伤度及疲劳寿命^[15]。2种方法的计算结果如表2所示。

5 结　语

本文运用规则波数值模拟试验法计算得到海上大尺度异型镂空浮体在波浪载荷下的应力响应后，分别应用热点应力法和名义应力法对该大型异形浮体结构进行疲劳损伤及疲劳寿命评估计算，得到以下结论：

1）浮体的疲劳位置基本都产生在焊缝位置，提高焊接技术有利提高结构的疲劳寿命。

2）当浮体受到150°方向的波浪荷载时，结构的应力值最大，可以基于此对结构进行针对性的加强。

3）采用热点应力法计算的结构疲劳寿命要小于采用名义应力法所计算的结构疲劳寿命，结果偏于保守。在未能准确对焊缝进行建模时，建议采取热点应力法计算结构的疲劳寿命。