基于ECOMSED模型的湛江湾水道三维潮流数值模拟

赵婉璐 , 郝瑞霞

(1.太原理工大学 水利科学与工程学院, 山西 太原 030024; 2.山西水利职业技术学院, 山西 运城 044004)

滨海地区近岸近海工程建设会在不同程度上影响原有工程区的水动力特性, 势必也会对该工程海域的自然生态环境造成一定的影响。因此, 合理预报工程建筑物对潮流运动特性的影响尤为重要。采用数值模拟方法对待开发区域内的海洋水动力进行研究是一个经济、有效的手段。由于计算机技术和数值求解技术的限制, 工程中常采用平面二维模型对潮流进行数值模拟[1-2], 这难以体现天然水体的三维特性。近年来, 随着工程实践对数值模拟要求的提高, 以及计算机技术的发展, 越来越多的工程采用三维数值模拟进行研究[3-5], 赵群[4]根据黄骅港外航道淤积问题的具体情况, 采用风浪模型SWAN (Simulation WAve Nearshore)和三维海洋紊流模式ECOMSED(Estuarine, Coastal and Ocean Modeling System with Sediments)对黄骅港海域的波浪、潮流、泥沙进行了模拟, 指出波浪是造成黄骅港外航道泥沙淤积的主要动力因素。郝瑞霞等[5]针对冷却水工程的实际问题, 对有横向来流条件下的表面湍浮力射流进行了三维的数值模拟, 成功地预测了由温度引起的分层流现象, 拓宽了湍流模型在该领域的应用。本文利用ECOMSED水动力模式对湛江湾上游水道的潮流进行数值模拟, 并用原体观测资料进行验证。

1 三维潮流数值模型及计算方法

1.1 ECOMSED简介

ECOMSED[6]模型是由 Blumhberg, Mellor等在POM和ECOM模型的基础上发展起来的一个较为成熟的浅海三维水动力学模型。该模型适用于河口及近岸的海洋模拟, 可以综合考虑潮流, 径流, 风力等水文气象因素作用下的水流, 盐度, 温度及泥沙等物理量的时空分布。模型采用模态分离技术, 外模态用于计算水平对流和扩散的二维变量, 内模态考虑垂向分层, 用于计算密度场为主因的三维变量, 从而有效地提高计算速度。ECOMSED模型在水平方向上采用 Arakawa C交错网格, 加密网格以更好地拟合岸线边界; 在垂直方向上采用σ坐标系统, 能够提高浅海海底地形的处理能力。本次模拟即是用水动力模块来模拟湛江湾湾顶水域的潮流运动。

1.2 控制方程

连续性方程:

雷诺时均动量方程:

其中,为水平对流速度,W为垂向速度。U为水平x轴向速度,V为水平y轴向速度;ρ0为基准密度,ρ为计算液体的密度;KM为垂向湍流掺混系数, 其大小决定速度的垂直分布;Fx和Fy为湍流扩散项,f为科氏参数。

采用静力学假设和 Boussinesq近似下的深度压强方程:

边界条件:

运动学边界条件遵循

在自由水面η(x,y)

在水底H(x,y)

动力学边界条件要求边界层的切应力应为速度梯度的函数。其中的(τox,τoy) 为自由表面摩擦力, (τbx,τby)为底面摩擦力。其他各参数基本意义详见文献[6]。

2 计算细节

计算区域为湛江湾湾顶霞山到石门之间的狭长水道, 地理坐标为 110°23′~110°29′E, 21°10′~21°24′N(图1), 整个水域面积60 km2, 受潮汐汊道发育影响, 水道地形南深北浅, 纵深相差达20 m。上游石门桥处海面突然束窄, 为了满足分辨率, 水平网格尺寸为 60m×60m。计算全域平面共有200×400个网格, 垂向采用σ坐标分为6层, 较好地反映了模拟区域岸线及地形变化。

计算水域地形采用1∶40000的湛江港海图, 经数字化获得网格水深后利用内插方法进行计算。模型以石门和霞山两侧作为整个模型的控制边界, 采用典型中潮时期的逐时潮位作为水流控制边界条件,潮位资料由珠江委水文水资源勘测中心提供。在冷启动模式下采用运行了 3个周期后的潮流数据进行分析。根据前人经验确定模型的几个关键参数为: 最小底摩擦系数BFRIC为0.0025; 底粗糙系数ZOB为0.02 m; 水平对流扩散参数HORCON取0.1; 垂向紊动参数UMOL取1×10–6m2/s。内模时间步长为10 s,外模步长为1 s。

图1 计算区域及网格示意图Fig.1 Computational domain and grid view

3 计算成果分析

本文对2002年8月21日12:00~8月22日13:00 (中潮)进行了三维数值模拟计算。在冷启动模式下, 将中潮期1个周期内的湛江边界的实测潮位值重复为6个周期的潮位系列进行潮流模拟, 所建模型的潮流场在运行3个潮周期后趋于稳定, 选择第5个周期的潮流数据进行分析。并与计算域内3个全潮观测点(V1~V3)的实测数据进行对比分析。3个潮流测点V1~V3见图1。

3.1 潮位计算结果

石门桥下游附近设有一个潮位测站(图1), 图2为模型计算潮位与该测站实测潮位的对比图。从图中可以看出二者基本吻合, 同时可以看出, 本海域一天中出现两次高潮和两次低潮, 半日周期相邻两潮期的高潮或低潮高度明显不相等, 且涨潮时间与落潮时间也不相等, 表现出典型的不规则半日潮性质。

3.2 潮流计算结果

湛江湾顶水道内的潮流特性主要受到外海潮流的动力作用, 具体到本模拟区域, 模型中海潮的涨落模拟主要受到边界点的潮位控制。

图2 潮位验证Fig.2 Verification of tidal level

图3~图5为V1~V3测站在计算时段内模拟的流速、流向与实测值的比较, 从图中可以看出计算潮流与实测潮流有良好的一致性, 表层、0.6H(H为测点水深)、底层流速与流向大小均较为接近。同时可以看出, 平均涨潮流速小于落潮流速, 涨潮平均历时大于落潮平均历时, 同一时刻, 表层流速大于底层流速, 表明模拟结果能够较好地反映湛江湾顶水道的潮流特性。

图3 V1测站的流速、流向实测与模拟数据对比Fig.3 Current speed and current direction comparison between observation and model at V1

图4 V2测站的流速、流向实测与模拟数据对比Fig.4 Current speed and current direction comparison between observation and model at V2

图5 V3测站的流速、流向实测与模拟数据对比Fig.5 Current speed and current direction comparison between observation and model at V3

图6是中潮涨急时海流的深度平均流场, 由流场图可以看出, 模拟海域的潮流流向受岸线和深槽走向控制, 多为南北向或近似为南北向。

图6 中潮期涨急时流场图Fig.6 Current field of the maximum flood of the medium tide

4 结论

本文采用 ECOMSED数值模型, 建立湛江湾顶海域的三维潮流数值模拟, 通过与原体资料的对比分析, 潮流场的计算结果与观测资料吻合良好, 计算出的潮流结果能够体现本研究海域的实际潮流状况。本次数值模拟说明, ECOMSED模型可以很好地复演湛江湾上游海域的潮流运动情况。运用该模型对本流域内的温度变化、盐度分布及泥沙运动情况的模拟效果是否满足要求有待进一步的研究。

[1]吴江航, 韩庆书.计算流体力学、方法和应用[M].北京:科学出版社, 1988: 1-3.

[2]曾平, 段杰辉, 黄柱崇, 等.二维流冰消融数学模型[J].水利学报, 1997, 5: 15-22.

[3]吴水波, 尹翠芳, 张乾, 等.近海三维数值模型简介[J].污染防治技术, 2010, 23(5): 17-29.

[4]赵群.基于SWAN和ECOMSED模式的大风作用下黄骅港波浪、潮流、泥沙的三维数值模拟[J].泥沙研究,2007, 4: 17-26.

[5]郝瑞霞, 周力行, 陈惠泉.冷却水工程中湍浮力射流的三维数值模拟[J].水动力学研究与进展(A辑), 1999,14(4): 484-492.

[6]HydroQual Inc.A Primer for ECOMSED (version 1.3)[R].NJ: HydroQual, Inc., 2002.