畸形波作用下的浮式结构水弹性响应分析

2022-10-29朱仁庆苑中排

船舶力学 2022年10期

朱仁庆，谢彤，刘一，苑中排

（江苏科技大学船舶与海洋工程学院，江苏镇江 212003）

0 引言

大型浮式海洋结构物（very large floating structure，简称VLFS）是开发海洋资源和利用海洋空间的重大海洋装备，可用作海上机场、岛礁浮式保障平台、海上航天发射场、漂浮海上生态城市等[1-3]。VLFS 体现两个特点：一是其某一维度的尺度要达到数公里；二是在波浪等外载荷作用下的弹性运动特别明显而必须采用水弹性理论对其载荷与响应进行求解和分析。

船舶与海洋工程水弹性理论自上世纪70 年代末提出到应用，已经从二维拓展到了三维、从线性拓展到了非线性、从势流理论向CFD领域扩展[4]。

目前，VLFS 水弹性响应分析大多是在线性波浪理论的框架下进行，少量考虑了波浪的弱非线性效应。当海洋工程向深海发展时，海洋结构物会不可避免遭遇畸形波（freak wave），这是一种波高极大、波峰异常尖瘦、能量相当集中的不规则波[5]，对海洋平台安全性是个重要威胁。在此海洋环境下，必须考虑强非线性波浪下结构的水弹性响应。Takag（i1996）[6]在研究受风暴或海啸引起的非线性长波作用下弹性浮板水弹性响应时发现，波浪的非线性对VLFS 水弹性响应影响不可忽略；Liu（2002）[7]基于BEM-FEM 方法模拟了柔性浮体结构在非线性波浪中的水弹性响应，通过与实验结果相比较，验证了方法的有效性；程勇（2017）[8]基于高阶BEM 方法建立了完全非线性的二维时域数值水池，用于模拟聚焦波浪与弹性浮板的水弹性响应，研究表明，浮体结构在极端波浪情况下的水弹性响应具有极强的非线性特征。以上研究仍是在势流理论框架下进行的，当遇到波浪涌上浮体上部以及波浪冲击浮体后波浪破碎等强非线性效应时，势流理论存在一定的局限性。相比于势流理论，计算流体力学（computational fluid dynamics，简称CFD）在描述流场中速度场与压力场的改变、捕捉自由液面非线性运动等方面更为有效和具有准确性。基于CFD理论，研究者们开发了CFD-FEM耦合方法来研究时域下浮体的水弹性响应。Seng（2012）[9]基于OpenFoam 开发了一种CFD-FEM 的耦合方法，其中将弹性体简化为梁模型进行计算，结果表明，该方法在计算弹性体水弹性响应方面具有较高的精度；Lee（2014）[10]开发了一种基于RANS 方法的CFD 和动态FEM 的单向和双向耦合算法，发现该方法可以很好地预报由规则波和不规则波激励引起的载荷；Qin（2017）[11]基于自开发的CFD-FEM 完全耦合的求解器，对比了弹性板垂直下落到非线性畸形波与规则波上结构响应的不同。

本文将大型浮体简化为带有自由边界的、材料为各向同性的弹性浮板，借助数值波浪水池，基于CFD-FEM方法，研究在畸形波浪作用下弹性浮板的水弹性响应。

1 数值模型

1.1 控制方程

（1）流体运动方程

当忽略流体的压缩性、不计流体表面张力时，描述流体运动的连续性方程与N-S方程为

式中，xi(i= 1,2,3 )为直角坐标系的三个分量；ui为xi方向速度；ρ为流体密度；t为时间；Fbi表示xi方向体积力；p为压强；τij表示粘性应力张量，其表达式为

其中，μ为粘度系数。

粘性流体运动存在湍流运动。目前处理湍流问题的方法主要有直接求解N-S 方程，又称为DNS法；雷诺平均法，即RANS法；大涡模拟法，即LES法。这里采用RANS法，其方程为

现在虽有多种湍流模型可供选择，但针对不同粘流问题需要选择合适的湍流模型。在数值模拟考虑粘流影响的波流问题时，SSTk-ω模型是个常用的选择。其在方程中增加了交叉扩散项，并且在湍流粘性系数中考虑了剪切应力的影响，使其计算稳定性好，计算效率与计算精度较高，故本文采用此模型。SSTk-ω模型的输运方程为

式中，k与ω分别表示湍流动能与湍流耗散率，ui与uj为速度分量，G͂k是指由于平均速度引起的湍动能k产生的项，Gω是单位耗散项，Γk和Γω分别是k和ω的有效扩散项，Yk和Yω分别是k和ω的湍动能耗散项，Dω是交叉扩散项，Sk与Sω是自定义源项。

（2）结构运动方程

假设结构为线弹性材料，在波浪等外载荷作用下相对于原平衡位置作刚体运动和变形，则线弹性结构经有限元离散后的总体运动方程为

式中，M为结构质量矩阵，C为结构阻尼矩阵，K为结构刚度矩阵，[x]为节点位移矩阵，[F(t)]为外界各种力合成的等效节点力列阵。

1.2 自由表面捕捉方法

在有网格的自由表面数值模拟方法中，常使用流体体积分数方法（volume of fluid method，简称VOF 法）来捕捉自由液面。交界面的相分布和位置由相体积分数αi来进行描述。网格单元中所有相的体积分数总和为1，意味着该单元内不存在空相，单元被所有相填充满，不存在空隙。相i的体积分数定义如下：

式中，Vi为网格单元中相i的体积，V为网格单元的体积，N为总相数。

根据体积分数的值，可以区分网格单元中是否存在不同相或流体：

包含交界面的网格单元中同时存在两种或者两种以上的物理相，该网格单元可视为混合相，其体积与动力粘度可以通过公式（12）进行计算：

式中，ρi为网格单元中相i的体积，μi是网格单元中相i的动力粘度。

1.3 畸形波浪模型

（1）组合聚焦模型

选取JONSWAP波浪谱来描述不规则海况，其表达式为

由于极限波聚焦模型中存在能量过于集中的问题，因此选取文献[12]中的极限波-随机波组合模型来实现畸形波浪的生成，表达式为

式中，通过指定水池位置xp，指定时刻tp叠加产生畸形波，N为叠加波浪的数量，坐标系同2.2节中的坐标系，当Ep=0 时为极限波浪模型，Ep=1 时为随机波浪模型。可通过改变Ep的数值来调整波谱能量在畸形波模拟中的分配问题，使模拟得到的波浪满足畸形波的特征；Ai、k i和ωi分别为第i个组成波的振幅、波数和圆频率；ϕi为随机相位，是(0,2π)内的随机数。

（2）仿物理造波

采用推板造波的方法来实现畸形波的生成。根据推板造波的原理，推板造波运动方程为

式中，Tr(ωi)为第i个组成波与造波板运动速度之间的运动传递函数，其表达式为

式中，h为数值水池的水深。

将畸形波的波面方程（14）代入公式（15）中，得到畸形波推板运动方程为

1.4 CFD-FEM 耦合方法

本文采用CFD 和FEM 之间的双向耦合方法，图1 为耦合计算流程图，图中t0代表初始时刻，Δt表示时间增量。在初始时刻使用CFD计算出的浮板表面压力以数据映射的方法传递给有限元模块中的浮板模型，在压力载荷的作用下，浮板有限元节点产生速度与加速度的变化将导致流固耦合交界面的变形；之后把变形后的节点数据传递给CFD计算程序，进行交界面的更新。由CFD计算出的压力场和速度场以及FEM 计算出的节点速度和加速度将传递给下一时间步。耦合计算中数据映射至关重要，因为CFD 与FEM 之间的网格离散不同，网格节点不相对应，计算中的映射采用形状函数插值的方法进行数据映射。

2 数值计算

2.1 数值模型验证

为了验证CFD-FEM 方法的有效性，将本文计算的浮体水弹性响应计算结果与文献[13]的试验结果以及文献[14]三维势流方法计算结果进行比较，其中弹性浮板模型参数为：板长L= 3.97 m，板宽B= 0.97 m，板厚D= 0.06 m，水深Dw= 2.3 m，吃水T= 0.037 m，弯曲刚度EI= 50.832 N·m，泊松比σ1= 0.3。弹性浮板的迎浪端距离速度入口造波边界的距离为4 m，水池宽度为4 m，设置浮板位于水池的中心；入射波波高为0.02 m，波长分别取2 m 和4 m；流场采用正六面体网格划分，弹性浮板采用四面体单元划分，并在自由液面上下波动空间和弹性浮板所在区域及附近流域内进行了网格加密。计算时间步为0.005 s，时间离散为2 阶精度，时间步内的迭代次数为10 次；监测弹性浮板迎浪端、中部以及背浪端位置处的垂向位移的变化情况。图2 为弹性浮板在迎浪端垂向位移最大时，浮板中线处的无因次垂向位移，并与试验数据[13]及势流计算结果[14]进行比较。由此可见，本文计算结果与两者的结果吻合较好，说明了采用CFD-FEM 方法模拟非线性波浪作用下超大型浮体水弹性响应的可行性。

2.2 畸形波生成

选取JONSWAP 谱的波浪参数为：谱峰圆频率ωp=3.491 rad/s（即谱峰周期TP= 1.8 s），叠加波浪的数量为N= 100，聚焦位置xp= 10 m，聚焦时间tp= 15 s。

数值波浪水池参数设置见图4。

在距离造波边界6 m、9 m、10 m与12 m处布置四个浪高仪来监测波面的升高。

在进行数值模拟时，计算初期数值计算的波面升高小于理论值计算结果，在图5（a）中尤为明显；但是在畸形波生成阶段，不同浪高仪位置处检测到的波浪升高与理论值贴合较好，在相位上无明显差异。

由图5（c）可观察到，在tp= 15 s 时，在设定的聚焦位置xp= 10 m 处生成了畸形波，基于畸形波的判定标准，在10 m 处浪高仪检测到的最大波高为0.128 m，与有义波高0.04 m 的比值为3.204，达到畸形波Hmax/H＞2.2 的标准，生成的波浪达到了畸形波的要求。在6 m 位置浪高仪检测到的最大波高为0.116 m，与理论计算的最大波高0.110 m 差异在5.1%；在9 m 位置检测到的最大波高为0.132 m，与理论值波高差异在4.7%。从四个不同位置处浪高仪所测值可知，粘性水池的数值解与理论解具有相当好的一致性。

2.3 畸形波浪中弹性浮板与刚性浮板运动响应对比

选取有义波高为H1/3= 0.06 m、谱峰周期TP= 1.8 s的波谱生成畸形波，畸形波生成位置xp= 10 m，聚焦时刻tp= 15 s，叠加波浪的数量为N= 100，能量集中频率区间为[ωL,ωH]=[2.318,6.891]。浮板模型参数与2.1节的相同，浮板迎浪端距离造波边界6 m，畸形波聚焦位置在浮板背浪端。计算区域设置如图6 所示，所有计算区域的设置以浮板的坐标系为基准，划分为3 个区域：第一个区域是造波区域，生成畸形波；第二个区域为外部流场计算区域，施加流固耦合交界面条件，以求解浮板表面流体压力载荷的变化；第三个区域为结构计算区域，施加结构位移运动条件，来求解在流体压力载荷下浮板运动和弹性变形。

图7为计算流域与弹性浮板结构单元划分图，图8给出了弹性与刚性浮板在迎浪端、中部与背浪端三个节点位置处浮板中线处垂向位移的变化。

从图8可见，刚性浮板垂向位移的幅值相比弹性浮板明显要小，这是因为弹性浮板除了刚体运动垂向位移外，还叠加了弹性垂弯的位移。

如图9 所示，对比甲板上浪的情况可以发现，刚性浮板的上浪情况比弹性浮板更严重，这是因为弹性浮板能够随着畸形波的运动产生垂弯位移，相应地减小了甲板高度与波高的差距，导致弹性浮板甲板上浪的水量少于刚性浮板甲板上浪的水量，所以刚性浮板的甲板上浪情况更严重。在15 s时，畸形波抵达聚焦位置10 m 处，弹性甲板上已没有积水，这是因为波峰离开迎浪端后，波谷到达迎浪端，而浮板会随着畸形波的运动产生垂弯，导致甲板上的积水顺着迎浪端的垂弯流入波谷；刚性浮板由于没有变形运动以减小甲板与波浪的高度差，且甲板上的积水不能及时排出，所以甲板上积水严重。

2.4 畸形波波高对弹性浮板水弹性影响

为了研究不同有义波高下弹性浮板的水弹性响应，这里选取H1/3= 0.05 m（Case-1），H1/3= 0.06 m（Case-2），H1/3= 0.07 m（Case-3）三个不同有义波高，以及迎浪与横浪两个浪向进行模拟。畸形波的其他参数与2.3 节的一致，对于迎浪工况，弹性浮板迎浪端距离造波边界6 m，畸形波聚焦位置在弹性浮板背浪端；对于横浪工况，弹性浮板前侧点距离造波边界9 m，畸形波聚焦位置在弹性浮板的后侧点。迎浪工况计算区域与2.3节的相同，横浪工况依此类推。浮板模型参数与2.1节的相同。

2.4.1 迎浪畸形波对弹性板的作用

在迎浪端、中部与背浪端三个节点位置处监测弹性浮板中线处垂向位移的变化情况。从图10可以看出，当聚焦位置和聚焦时间不变，聚焦波的有义波高增加时，弹性浮板垂向位移的变化幅值增大，弹性浮板的水弹性响应运动加剧,这是由于当聚焦波波高增大时，浮板各点处受到的波浪载荷增加，浮板的水弹性响应加剧，随之影响浮板各位置处的垂向位移数值增加。图中三种工况下垂向位移的幅值虽有差异，但是在整体的变形形态上是保持一致的，这是因为入射波浪的组成波的幅值虽然不同，但是组成波的相位与频率是相同的，使畸形波在波浪形态上存在一定程度的相似，波浪载荷作用于浮板上使之出现相似的变形形态。

2.4.2 横浪畸形波对弹性板的作用

进一步研究在横浪时弹性浮板的垂向位移变化情况。计算中选取浮板前侧点（P1）、中部点（P2）与后侧点（P3）监测垂向位移的变化情况，具体结果如图11所示。

改变入射波的有义波高，计算结果如图12 所示。在聚焦时刻15 s 附近时，聚焦位置处P3 监测点的垂向位移的峰值分别为0.086 m、0.099 m 与0.107 m，垂向位移随着入射波波高的增加而增加。

相比于迎浪，横浪下浮板的刚体运动占总体运动的较大成分。从图13 中Hw= 0.06 m 的波高下的侧视图中可以看出，浮板侧面的各点基本保持在一条直线上，弹性变形较小。这是由于在横浪时弹性浮板宽度方向尺度较小，浮板侧面的弯曲刚度相比于正面的弯曲刚度要大，从而导致其弹性变形量减小，刚体运动成为浮板运动的主要成分。

2.4.3 上浪情况研究

浪向角为0°与90°时都出现波浪涌上板面的情况（见图14），入射波高越大，上浪的水量越多。虽然浮板在遭遇波浪时存在一定的抬升，但是波浪的高度要大于浮板的抬升高度，从而出现波浪涌上浮板上部的现象。在浮板迎浪端抬升时，由于浮板上部存在一定的水体，其抬升高度将受到一定的影响，其峰值高度将略小于无上浪时的峰值变化。

2.5 畸形波作用位置对弹性浮板水弹性影响

为了研究畸形波浪位置对浮板水弹性作用的影响，选取了3种不同计算工况：浮板的迎浪端位于聚焦位置（Condition-1），浮板的中部位于聚焦位置（Condition-2）以及浮板背浪端位于聚焦位置（Condition-3）。图15 显示了三种不同工况下弹性浮板迎浪端与造波边界的距离。采用有义波高H1/3=0.04 m 的畸形波进行计算，畸形波其他参数与2.3节的一致。浮板模型参数与2.1 节的相同。从3 种工况可以看出，弹性浮板的水弹性响应峰值出现时刻依次逐渐提前。这是由于弹性浮板离造波边界的距离从Condition-1 到Condition-3 逐渐减小，遭遇到畸形波浪的时间逐渐提前。图16 为Condition-1 和Condition-3中垂向位移各峰值的对比，可以看出，畸形波聚焦位置处都出现了垂向位移的峰值，同时注意到在Condition-2 中，聚焦位置在浮板中部，但是浮板水弹性响应位移最大值却出现在迎浪端位置处，背浪端的位移最大值比迎浪端的略小，这是因为浮板中部位置是浮板的内部节点，当浮板出现变形时，该节点处会受到与之相邻的其他节点的影响，在一定程度上使其运动状态减弱。

2.6 不同谱宽畸形波对弹性浮板水弹性影响

为研究不同频带宽度下畸形波与弹性浮板的相互作用，选取中心频率为fc=0.56 Hz，即谱峰周期为TP=1.8 s，模拟频谱宽度为(a)Δf=0.3 Hz，(b)Δf=0.6 Hz，(c)Δf=0.9 Hz三种工况下浮板的水弹性响应。参照文献[12]选取频率区间，即以中心频率fc=0.56 Hz为基准加上或减去0.5Δf来获取频率区间的上下限，有义波高Hs=0.05 m，畸形波其他参数与2.3节的一致。浮板模型参数与2.1节的相同。

图17中为有义波高Hs= 0.05 m 时不同频带宽度条件下弹性浮板迎浪端、中部和背浪端的垂向位移的时历曲线。在图中可见，当选取的频带宽度较窄时，浮板的垂向位移大于频带较宽的情况，这是由于频带宽度是能量集中的外在表现。当频带越窄时，单位频带所含有的能量越高，从而导致波浪波幅、波浪载荷增加，使浮体的水弹性响应加剧，在同一时刻同一位置处的位移增大。