浅海风-浪-流-海床耦合场非定常时空演化规律及评价指标

陈 静, 柯世堂, 李文杰, 朱庭瑞, 员亦雯, 任贺贺

(南京航空航天大学 土木与机场工程系; 江苏省机场基础设施安全工程研究中心, 南京 211106)

随着我国海洋强国战略[1]的快速推进,近海岸海洋工程建设需求与日俱增,南海海域经历复杂的地质演化过程形成不同类型的海床地貌,其三维效应对风-浪-流耦合场演化规律具有显著影响,导致近岸风、浪、流场的时空分布及统计特性不同于开阔海域,准确预测复杂海洋环境是海岸工程安全保障的重要前提.目前国内外规范[2-3]关于浅海区波浪的计算方法是基于线性波浪理论假定,不能考虑海床面粗糙度对波流运动的影响,工程设计中海洋环境预测结果仍只能过度依赖实测数据校准和关键参数取值经验.随着计算机技术高速发展,计算流体力学(CFD)数值模拟现在已成为风-浪-流环境耦合场预测的重要手段,但目前模拟采用的边界造波法忽略了初始场未形成已改变风浪流要素的问题,一定程度限制了浅海地区风-浪-流耦合场的预测精度.

现有关于风-浪-流耦合场研究主要集中在数值方法改进[4-6]、流场演化规律[7-10]、结构水动力模型[11-14]等方面,文献[9]中基于Fluent平台二次开发研究发现风-浪-流三者耦合作用下,风会导致波浪传播过程中波长增大、波幅减小;文献[6]中结合理论推导将波浪水质点速度化作单位矢量代入风剖面公式之中,改进了风-浪耦合边界自由液面处风与水质点速度方向的一致性.但上述模拟过程中风-浪-流耦合模型通常假定风速、流速沿高度均匀分布和水深不变,无法考虑浅海地区梯度风、剪切流的垂向速度不均匀和海床地貌不规则对耦合场的影响.针对海床地貌的影响研究,仅有部分学者对浅水波浪的非线性变化[15]或斜坡地形的波浪破碎[16-21]问题进行探索,文献[19]中结合物理模型实验结果分析了礁前坡度、坡顶淹没水深与规则波在岛礁地形上的破碎指标之间的相关性;文献[20]中指出随海床面坡度的增大,孤立波在传播过程中受地形因素导致的能量耗散作用为主,浅化作用为辅;文献[21]中为研究孤立波与斜坡作用,利用新一代CFD软件STAR CCM+以造波边界设置多层流体速度方式构建二维内孤立波数值水槽,模拟波浪振幅与设计振幅误差在0.5%以内.此类研究大都局限于斜坡海床下波浪场的破碎现象,未涉及真实海域典型海床地貌与风-浪-流耦合场的非定常演变物理过程研究,且缺乏基于不同海床对风-浪-流场非定常演变特征的评价标准.

鉴于此,基于STAR-CCM+二次开发考虑了梯度风、剪切流的垂向速度不均匀对耦合场的影响,开展多层三维黏性风-浪-流-海床耦合场仿真模拟,分析4种典型海床地貌下波浪、海流和风场的非定常时空演变规律,对比分析了波高、波长、平均流速、流剖面指数、基本风速及风剖面指数的差异及产生原因,最后结合主成分分析法提出各海床工况下风-浪-流耦合场全生命周期非定常综合评价指标,可为浅海地区典型海床地貌下的风-浪-流场精准预测提供参考.

1 风-浪-流-海床耦合场模拟方法

1.1 数值模拟方法

1.1.1典型海床概况 我国南海海域包含着纷繁复杂的地貌单元,整体发展趋势从四周向中部水深逐渐变大,将陆相-海陆过渡相-海相的发展区域地貌单元依次划分为大陆架、大陆坡、边缘海盆地[22-23].研究区域位于南海海域的东北位置,属于大陆架与大陆坡的交界区域,重点以大陆坡典型海床地貌作为研究对象,其中大陆坡的三级地貌共分为海底平原、海底斜坡和海槽,以3种地貌平均深度较低的四级地貌[23]:海底平原的深海水道、海底斜坡的褶皱型“微凹陷”和海槽地带的“U” 形“微凹陷”作为典型地貌,再以无地貌的平坦海床作为对比工况,以此来研究典型海床地貌下波浪场、风场与海流场演化过程.图1给出了海底平原、海底斜坡和海槽典型地貌的数值模型,床面高度指海床最大高度差,分别为18、5、10 m;褶皱角度(30°)指褶皱顶部角度,U弯角度(15°)指中央拱起角度.

图1 南海东北部海域海床地貌及简化模型图Fig.1 Seabed topography and simplified model map in northeastern South China Sea

1.1.2黏性流体流动数学模型 通常假设浮体周围流域内的流体是不可压缩性的,不考虑流体的表面张力,可结合连续性方程和N-S方程来表达流体的黏性运动[24]:

(1)

(2)

式中:u为三维流体x、y、z方向的速度矢量;ρ为密度;t为流体运动时间;μ为动力黏度系数与涡动黏度系数之和;p为压力;g为重力矢量;fs为消波源项.

海床面的运动边界条件:

(3)

式中:D为海水深度.

1.1.3改进多层风-浪-流耦合模型 线性波浪波面高程方程为

(4)

式中:η为波面高程;H为波高;k为波数;ω为圆频率.

以五阶Stokes波模拟周期性波浪,该波浪波面具有非对称余弦曲线函数状态,可直观观察到波峰与波谷形态,比线性波浪更接近实际海况.其波面高程方程为

(5)

x方向的速度表达式:

ux=

(6)

z方向的速度表达式:

uz=

(7)

式中:d为水深;c=cosh(kd);λ1=λ,λ2=λ2B22+λ4B24,λ3=λ2B33+λ4B35,λ4=λ2B44,λ5=λ2B55;λ、B系列参数详见文献[25].

风场依据《建筑结构荷载规范》[26]采用A类地貌风剖面系数,以指数率梯度风进行定义:

(8)

式中:uw、u10分别为高度z和高度10 m对应的平均风速;α为地面粗糙度指数,A类地貌取值0.12.

自由液面以下的海流流速以1/7定律的指数型分布表示:

(9)

式中:uc、vc分别为海流的流速和平均流速;l为距海床底距离.

实际海域中水流速度在垂向上具有一定的分布差异,从而导致波浪沿不同水深发生变化;同时风在垂向上也呈现梯度分布,导致波面传播会产生变形现象.为研究近海风-浪-流演化问题,需要同时考虑风、流垂向分布的影响.以多层质点速度耦合方法应用于风浪流共存情况,如图2所示.图中:在多层风浪流耦合模型中建立笛卡尔坐标系OXYZ;原点O为海平面位置;Z坐标描述水深;Y坐标描述风-浪-流横向宽度;风-浪-流沿着X正方向传播.对计算域进行风波耦合区、波流耦合区、水深一半以上区和水深一半以下区的分层,不同层的风速、流速以风剖面、剪切流的速度进行分解输入,以此建立新型多层风-浪-流耦合模型.该模型的速度uwj、uci分别为j、i层的风浪耦合速度和浪流耦合速度.风-浪-流在x方向上的速度:

图2 多层风-浪-流耦合场模型示意图Fig.2 Diagram of multilayer wind-wave-current coupling field model

Ux=

(10)

1.2 CFD模拟参数设定

图3给出了风-浪-流-海床三维数值黏性水池模型示意图.计算域由风-浪-流生成区、阻尼消波区和风-浪-流-海床耦合演化区组成.为保证流体充分发展,计算流体域取400 m×200 m×330 m(流向x×展向y×纵向z),水深按照实际海深取值30 m,尾端设置尺寸为66 m×200 m×330 m的消波区.

图3 多层风-浪-流-海床三维数值黏性水池模型示意图Fig.3 3D numerical viscous pool model of multilayer wind-wave-current-seabed

计算域边界条件设置如下:左边界为风-浪-流速度入口,右边界设为压力出口,前后边界设为对称边界,顶部边界设为滑移壁面,底部边界设为无滑移壁面,其中为了流体得以更快地形成稳定场,速度入口和压力出口设置成周期性接触,形成充分发展性边界.雷诺平均N-S方程(RANS)非定常计算模型采用更适合于求解雷诺应力分量输运方程的SSTk-ω湍流模型,使用流体体积方法对自由液面进行追踪处理.通过自定义场函数接口,将多层风浪流耦合场模型施加在入口边界上,初始化时停用第一波前设置,可生成对应的风浪流初始场.为防止出口风浪流反射,形成壁面效应,在水池尾部添加动量源阻尼项.根据DNV-RP-C205规范[3]中统计的各海区风浪流长期分布参数,选择研究区域属于40号海区,选取海区作业海况作为输入条件,波高设置为2.1 m,周期为5.3 s,设10 m处风速为11.4 m/s,平均流速为0.68 m/s,通过波浪速度分解x、y和z向,与自定义风剖面、剪切流叠加计算.

在计算域全域沿高度定义波峰以上、波谷以下及波浪发展范围内的风-浪-流解耦速度(式(10)),初始时刻计算域全域即实现了风-浪-流场解耦,解决现阶段数值造波方法存在的初始场形成前,风、浪、流要素已改变问题,同时提高计算效率和精度.图4给出初始化流场分布图.

图4 风-浪-流水池初始化速度场云图Fig.4 Initial velocity field cloud image of wind-wave-flow tank

图5给出了风-浪-流-海床耦合数值水池网格划分示意图.整体采用混合网格离散形式,以切割型六面体网格加密程度区分液面加密域、液面过渡域和其他区域.其中对液面加密域沿波浪传播方向、波高方向加密网格,液面过渡区相对于液面加密区两倍的网格尺寸进行加密.y+值为无量纲壁面距离,关乎近壁面湍流的发展状态,近壁面区域使用全y+壁面处理函数方法,近壁面第1层网格厚度为0.041 m,y+范围在150～250内,满足壁面函数法要求.

图5 风-浪-流-海床三维数值粘性水池网格划分示意图Fig.5 Diagram of wind-wave-stream-seabed three-dimensional numerical viscous pool grid

表1给出了不同水线加密工况下参数对比,其中L为波浪波长.由表可知:随着网格数量的增加,波高误差逐渐减小趋势,波高衰减情况改善.工况4和工况5的波高误差与规范结果相差较小,且两者之间也无明显差异,综合考虑计算精度与效率的平衡,本文选取网格工况4为数值模拟网格方案,时间步长按照T/(2.4r)取值.其中:T为波浪周期;r为一个波长的网格数量.时间步长为0.022 s.

表1 网格划分方案工况Tab.1 Working condition of grid partitioning scheme

1.3 有效性验证

风-浪-流参数模拟精确度与试验结果有很大的关联性,为验证数值水池模拟试验结果,分别对波浪单独作用的波浪场及消波效果、对初始时刻风-浪-流耦合作用下的风场及海流场进行分析,通过距入口、出口20 m处监测的波面高程、纵向分布风速及流速的模拟值与理论值进行对比,来验证试验模拟的准确性,其中风-浪-流规范值参数如表2所示.

表2 风-浪-流参数Tab.2 Wind-wave-current parameter

图6(a)为波浪单独作用时距入口20 m处波浪模拟值与理论值对比曲线,在模拟的60～90 s阶段,波浪的模拟值与理论值波形基本吻合,且无衰减的趋势;图6(b)为距出口20 m消波区波面高程时程曲线,由图可见相较于未消波和消波后的理论值,波浪在两个周期以后已基本消失,验证了阻尼消波的可行性.

图6 波高模拟值与理论值对比时程曲线Fig.6 Comparison of time-history of wave height simulation value with theoretical value

图7给出了风-浪-流耦合作用下初始时刻计算域距入口20 m处沿高度风速和流速模拟结果对比曲线.结果表明入口处风剖面和剪切流剖面与理论解吻合,初始条件下风场和海流场达到稳定.

图7 距入口20 m处风速及流速剖面模拟值与理论值对比曲线Fig.7 Comparison of simulated and theoretical values of wind speed and current speed at 20 m from inlet

2 耦合场时空演化规律

2.1 波浪场时空演化

图8给出了计算域中间(x=0 m)4种海床地貌下波浪场波形变化时程曲线.对比分析发现风-浪-流场与平坦地貌、海底斜坡和海槽的耦合作用使得波高大幅度减小、波长变长;平坦地貌、海底斜坡和海槽工况呈现衰减演化过程,海底平原工况呈现波浪破碎过程.平坦地貌、海底斜坡和海槽3种地貌下的衰减演化过程均可分为波面激增、衰减和稳定3个阶段,其中衰减阶段平坦地貌和海底斜坡工况波浪场处于持续衰减,而海槽地貌由于“浅-深-浅”的地形特征和色散定律,波浪场形成波群性衰减;海底平原工况先经历外破波和内破波两个阶段的波浪破碎,呈现波高突变和波形衰减的非线性变化,最后在爬波阶段受到触底回流的冲击,形成最后一次高达6 m的“直立”形波浪破碎现象.

图8 不同海床地貌下波形高程-时程变化曲线Fig.8 Waveform time-history curves in different seabed topographies

图9～12给出了演化阶段下各海床工况全域空间波形变化示意图.由图可知:平坦地貌、海底斜坡和海槽3个工况经历过激增阶段的波面陡立后,都呈现不同程度的衰减现象,其中海槽地形也由于色散定理和水深特征,波浪的前缘速度和后缘速度不同,波面的横向不对称性相比其他工况较明显;海底平原工况外破波阶段波浪破碎形态多为激破波与卷破波组合的破碎形式,内破波阶段形成出口端波浪折射现象和部分前沿陡立、后坡平坦的段波,与爬坡阶段的触底回流碰撞形成崩破波.

图9 平坦地貌不同阶段三维空间波形演变云图Fig.9 Three-dimensional spatial waveform variation diagram under flat geomorphic conditions

图10 海底平原不同阶段三维空间波形变化图Fig.10 Three-dimensional spatial waveform variation diagram under submarine plain working conditions

图11 海底斜坡不同阶段三维空间波形变化图Fig.11 Three-dimensional spatial waveform variation diagram under submarine slope conditions

图12 海槽不同阶段三维空间波形变化图Fig.12 Three-dimensional spatial waveform variation diagram under trough working conditions

2.2 海流场时空演化

图13为各工况的剪切流模拟结果对比,其中每幅下图对应上图的最后时刻,β值为流剖面指数,初始为1/7.对比可见:与平坦地貌相反,海底平原、海底斜坡和海槽工况对流剖面要素均存在一定影响,其中平均流速皆呈现随时间逐渐减小的趋势;海底平原工况由于触底回流的原因,演化终期垂向流速基本为负值.分别采用指数型、多项式对一半水深以下、一半水深以上和波流耦合区3个区域进行拟合,并给出如下4个工况拟合后的垂向流速解析式.

教师主要以课堂教学的知识和理论基础为基础，与直接获得实际和实践能力的生产和科研实践相分离。因而，通过实践教育基地的建设，给教师搭建了产学研合作的平台，不仅扩宽了教师的视野，实现知识的转化，而且还能提高教师的科研水平；同时，企业也获得了人才、技术的有力支持，对于提高企业新技术、新产品的开发有了进一步的保障。

图13 不同海床地貌下剪切流模拟结果对比Fig.13 Comparison of shear flow simulation results in different seabed geomorphologies

平坦地貌:

(11)

海底平原:

(12)

海底斜坡:

(13)

海槽:

(14)

图14～17给出各工况下速度流线空间分布图.分析可得:不同于传统的风-浪-流耦合场波峰下部“V”形、波谷下部“U”形间隔涡区的研究结果,耦合场由于垂向压差,形成了不同尺度的周期性涡旋;平坦地貌工况地形对涡旋发展无阻碍,形成“多涡积聚”,海底斜坡工况“褶皱”地形的阻碍导致流场形成“多涡共存”的现象;而海底平原工况由于卷破波形成时部分空气卷入,加剧低空处风和波浪速度,波浪场加剧变化使得全域波形呈现“坡形”,风场与其碰撞形成图15(b)中风波耦合区的涡旋,爬坡阶段涡旋发展过程遇到触底回流,海流场涡旋体系被破坏;对于海槽工况,海流场发展空间较小,但经过长时间的演化,形成的小涡旋逐渐破坏初始“U+V”涡区形态.

图14 平坦地貌工况流场时空演化图Fig.14 Space-time evolution diagram of flow field under flat landform conditions

图15 海底平原工况流场时空演化图Fig.15 Space-time evolution diagram of flow field in submarine plain

图16 海底斜坡工况流场时空演化图Fig.16 Temporal and spatial evolution diagram of flow field under submarine slope conditions

图17 海槽工况流场时空演化图Fig.17 Temporal and spatial evolution diagram of flow field in trough conditions

2.3 风场时空演化

图18为不同工况下风剖面拟合结果对比,其中每幅下图对应上图的最后时刻.对比发现4种海床地貌引起的波浪场变化均会造成风剖面指数的增大和基本风速的减小,且参考表1,发现海床床面高度与风剖面指数成正比关系,与A类地貌取值0.12[26]相比,平坦地貌、海底平原、海底斜坡和海槽的演化终期风剖面指数分别为0.164、0.197、0.173和 0.176;海底平原工况由于外破波阶段波高全域性突变至7～8 m,海表面的粗糙度大幅度提高,此时风场也在波面抬升中损失大量能量,导致图18(b)中在t=7 s附近基本风速降到极小值5 m/s,风剖面指数提高到极大值0.4,在后面的演化阶段,风场恢复至稳定,但其基本风速的稳定值小于其他3个工况.

图18 不同海床地貌下风剖面模拟结果对比Fig.18 Comparison of simulation results of wind profile in different seabed topographies

3 耦合场非定常效应评价指标

(15)

结合主成分分析方法,根据风-浪-流演化过程中每0.02 s模拟出的6个非定常环境因子作为样本,提取最大贡献主成分.图19分别为4个海床工况的非定常环境因子在主成分1×主成分2 上的投影.图中:圆点的大小与演化时间成正比;箭头终点的坐标值代表各非定常因子对主成分的贡献值.表3为不同海床地貌下主成分最大贡献数量及累计贡献率.其中平坦地貌、海底斜坡及海槽工况前3个主成分的累计贡献率分别达到91.7%、85.6%和94.6%;海底平原工况所需4个主成分,累计贡献率达到90.0%,已满足主成分分析法降维精度要求[27].

表3 不同工况主成分贡献数量及累计贡献率Tab.3 Number and cumulative contribution rate of principal components under different working conditions

图19 非定常环境因子和样本数据对主成分贡献数值图Fig.19 Unsteady environmental factors and sample data contribute to numerical diagram of principal components

据此可建立风-浪-流全生命周期非定常综合评价函数:

F=

(16)

式中:X1～X6依次为波高、波长、平均流速、流剖面指数、基本风速与风剖面指数的非定常环境因子.基于此取演化终期的非定常环境因子数值,建立各海床地貌工况下风-浪-流耦合非定常评价指标参考值.得到平坦地貌、海底平原、海底斜坡和海槽的非定常评价指标分别为0.268、4.612、0.672和0.926.

4 结论

系统研究了典型海床地貌下风-浪-流耦合过程非定常演化规律及评价指标,尤其在多层风-浪-流耦合模型、风-浪-流场阶段性发展形态演化和非定常综合评价指标方面取得了原始创新,具体如下:

(1) 新型多层风-浪-流耦合模型能精确模拟出垂向风速、流速不均匀对风-浪-流耦合场造成的影响.以风-浪-流解耦模型运用于CFD初始条件中,有效提高计算效率和精度.

(2) 各海况的风-浪-流-海床演化过程:波浪场演变分为三阶段,其中平坦地貌、海底斜坡和海槽工况分为波面激增、衰减和稳定阶段;海底平原工况共分为外破波、内破波和爬坡阶段.海流场随演变过程在空间上呈多段式分布,一半水深以下的指数率函数、一半水深以上及波流耦合区的多次项函数.

(3) 海流速度场由于垂向压差形成涡旋,各工况根据海床地貌不同形成多涡积聚或多涡共存的现象;风-浪-流-海床耦合演化对风剖面指数产生放大作用,海床床面高度与风剖面系数成正比关系.

(4) 建立了基于主成分分析法的耦合场全生命周期非定常评价指标, 并提出演化终期各典型海床下非定常评价指标数值:平坦地貌、海底平原、海底斜坡和海槽分别为0.268、4.612、0.672和0.926.非定常评价指标可精准预测浅海地区典型海床地貌下风、浪、流时空分布变异性特征, 可为我国海岸工程安全建设中海洋环境预测提供参考.