基于ECOMSED模型的湛江灣水道三維潮流數(shù)值模擬

趙婉璐 , 郝瑞霞

(1.太原理工大學 水利科學與工程學院, 山西 太原 030024; 2.山西水利職業(yè)技術學院, 山西 運城 044004)

濱海地區(qū)近岸近海工程建設會在不同程度上影響原有工程區(qū)的水動力特性, 勢必也會對該工程海域的自然生態(tài)環(huán)境造成一定的影響。因此, 合理預報工程建筑物對潮流運動特性的影響尤為重要。采用數(shù)值模擬方法對待開發(fā)區(qū)域內(nèi)的海洋水動力進行研究是一個經(jīng)濟、有效的手段。由于計算機技術和數(shù)值求解技術的限制, 工程中常采用平面二維模型對潮流進行數(shù)值模擬[1-2], 這難以體現(xiàn)天然水體的三維特性。近年來, 隨著工程實踐對數(shù)值模擬要求的提高, 以及計算機技術的發(fā)展, 越來越多的工程采用三維數(shù)值模擬進行研究[3-5], 趙群[4]根據(jù)黃驊港外航道淤積問題的具體情況, 采用風浪模型SWAN (Simulation WAve Nearshore)和三維海洋紊流模式ECOMSED(Estuarine, Coastal and Ocean Modeling System with Sediments)對黃驊港海域的波浪、潮流、泥沙進行了模擬, 指出波浪是造成黃驊港外航道泥沙淤積的主要動力因素。郝瑞霞等[5]針對冷卻水工程的實際問題, 對有橫向來流條件下的表面湍浮力射流進行了三維的數(shù)值模擬, 成功地預測了由溫度引起的分層流現(xiàn)象, 拓寬了湍流模型在該領域的應用。本文利用ECOMSED水動力模式對湛江灣上游水道的潮流進行數(shù)值模擬, 并用原體觀測資料進行驗證。

1 三維潮流數(shù)值模型及計算方法

1.1 ECOMSED簡介

ECOMSED[6]模型是由 Blumhberg, Mellor等在POM和ECOM模型的基礎上發(fā)展起來的一個較為成熟的淺海三維水動力學模型。該模型適用于河口及近岸的海洋模擬, 可以綜合考慮潮流, 徑流, 風力等水文氣象因素作用下的水流, 鹽度, 溫度及泥沙等物理量的時空分布。模型采用模態(tài)分離技術, 外模態(tài)用于計算水平對流和擴散的二維變量, 內(nèi)模態(tài)考慮垂向分層, 用于計算密度場為主因的三維變量, 從而有效地提高計算速度。ECOMSED模型在水平方向上采用 Arakawa C交錯網(wǎng)格, 加密網(wǎng)格以更好地擬合岸線邊界; 在垂直方向上采用σ坐標系統(tǒng), 能夠提高淺海海底地形的處理能力。本次模擬即是用水動力模塊來模擬湛江灣灣頂水域的潮流運動。

1.2 控制方程

連續(xù)性方程:

雷諾時均動量方程:

其中,為水平對流速度,W為垂向速度。U為水平x軸向速度,V為水平y(tǒng)軸向速度;ρ0為基準密度,ρ為計算液體的密度;KM為垂向湍流摻混系數(shù), 其大小決定速度的垂直分布;Fx和Fy為湍流擴散項,f為科氏參數(shù)。

采用靜力學假設和 Boussinesq近似下的深度壓強方程:

邊界條件:

運動學邊界條件遵循

在自由水面η(x,y)

在水底H(x,y)

動力學邊界條件要求邊界層的切應力應為速度梯度的函數(shù)。其中的(τox,τoy) 為自由表面摩擦力, (τbx,τby)為底面摩擦力。其他各參數(shù)基本意義詳見文獻[6]。

2 計算細節(jié)

計算區(qū)域為湛江灣灣頂霞山到石門之間的狹長水道, 地理坐標為 110°23′~110°29′E, 21°10′~21°24′N(圖1), 整個水域面積60 km2, 受潮汐汊道發(fā)育影響, 水道地形南深北淺, 縱深相差達20 m。上游石門橋處海面突然束窄, 為了滿足分辨率, 水平網(wǎng)格尺寸為 60m×60m。計算全域平面共有200×400個網(wǎng)格, 垂向采用σ坐標分為6層, 較好地反映了模擬區(qū)域岸線及地形變化。

計算水域地形采用1∶40000的湛江港海圖, 經(jīng)數(shù)字化獲得網(wǎng)格水深后利用內(nèi)插方法進行計算。模型以石門和霞山兩側作為整個模型的控制邊界, 采用典型中潮時期的逐時潮位作為水流控制邊界條件,潮位資料由珠江委水文水資源勘測中心提供。在冷啟動模式下采用運行了 3個周期后的潮流數(shù)據(jù)進行分析。根據(jù)前人經(jīng)驗確定模型的幾個關鍵參數(shù)為: 最小底摩擦系數(shù)BFRIC為0.0025; 底粗糙系數(shù)ZOB為0.02 m; 水平對流擴散參數(shù)HORCON取0.1; 垂向紊動參數(shù)UMOL取1×10–6m2/s。內(nèi)模時間步長為10 s,外模步長為1 s。

圖1 計算區(qū)域及網(wǎng)格示意圖Fig.1 Computational domain and grid view

3 計算成果分析

本文對2002年8月21日12:00~8月22日13:00 (中潮)進行了三維數(shù)值模擬計算。在冷啟動模式下, 將中潮期1個周期內(nèi)的湛江邊界的實測潮位值重復為6個周期的潮位系列進行潮流模擬, 所建模型的潮流場在運行3個潮周期后趨于穩(wěn)定, 選擇第5個周期的潮流數(shù)據(jù)進行分析。并與計算域內(nèi)3個全潮觀測點(V1~V3)的實測數(shù)據(jù)進行對比分析。3個潮流測點V1~V3見圖1。

3.1 潮位計算結果

石門橋下游附近設有一個潮位測站(圖1), 圖2為模型計算潮位與該測站實測潮位的對比圖。從圖中可以看出二者基本吻合, 同時可以看出, 本海域一天中出現(xiàn)兩次高潮和兩次低潮, 半日周期相鄰兩潮期的高潮或低潮高度明顯不相等, 且漲潮時間與落潮時間也不相等, 表現(xiàn)出典型的不規(guī)則半日潮性質。

3.2 潮流計算結果

湛江灣頂水道內(nèi)的潮流特性主要受到外海潮流的動力作用, 具體到本模擬區(qū)域, 模型中海潮的漲落模擬主要受到邊界點的潮位控制。

圖2 潮位驗證Fig.2 Verification of tidal level

圖3~圖5為V1~V3測站在計算時段內(nèi)模擬的流速、流向與實測值的比較, 從圖中可以看出計算潮流與實測潮流有良好的一致性, 表層、0.6H(H為測點水深)、底層流速與流向大小均較為接近。同時可以看出, 平均漲潮流速小于落潮流速, 漲潮平均歷時大于落潮平均歷時, 同一時刻, 表層流速大于底層流速, 表明模擬結果能夠較好地反映湛江灣頂水道的潮流特性。

圖3 V1測站的流速、流向實測與模擬數(shù)據(jù)對比Fig.3 Current speed and current direction comparison between observation and model at V1

圖4 V2測站的流速、流向實測與模擬數(shù)據(jù)對比Fig.4 Current speed and current direction comparison between observation and model at V2

圖5 V3測站的流速、流向實測與模擬數(shù)據(jù)對比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數(shù)值模型, 建立湛江灣頂海域的三維潮流數(shù)值模擬, 通過與原體資料的對比分析, 潮流場的計算結果與觀測資料吻合良好, 計算出的潮流結果能夠體現(xiàn)本研究海域的實際潮流狀況。本次數(shù)值模擬說明, ECOMSED模型可以很好地復演湛江灣上游海域的潮流運動情況。運用該模型對本流域內(nèi)的溫度變化、鹽度分布及泥沙運動情況的模擬效果是否滿足要求有待進一步的研究。

[1]吳江航, 韓慶書.計算流體力學、方法和應用[M].北京:科學出版社, 1988: 1-3.

[2]曾平, 段杰輝, 黃柱崇, 等.二維流冰消融數(shù)學模型[J].水利學報, 1997, 5: 15-22.

[3]吳水波, 尹翠芳, 張乾, 等.近海三維數(shù)值模型簡介[J].污染防治技術, 2010, 23(5): 17-29.

[4]趙群.基于SWAN和ECOMSED模式的大風作用下黃驊港波浪、潮流、泥沙的三維數(shù)值模擬[J].泥沙研究,2007, 4: 17-26.

[5]郝瑞霞, 周力行, 陳惠泉.冷卻水工程中湍浮力射流的三維數(shù)值模擬[J].水動力學研究與進展(A輯), 1999,14(4): 484-492.

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