一種有限體積懸沙輸運模型在蘇北輻射沙脊群海域的應用

孫大威，張 鑫

（1.內蒙古民族大學，通遼028000；2河北省水運工程規劃設計院，天津300074）

一種有限體積懸沙輸運模型在蘇北輻射沙脊群海域的應用

孫大威1，張 鑫2

（1.內蒙古民族大學，通遼028000；2河北省水運工程規劃設計院，天津300074）

懸沙輸運是近岸水域地形演變的重要因素之一，南黃海輻射沙脊群海域懸沙在波浪和潮流的共同作用下有著巨量的輸運。文章構建了基于無結構三角形網格的Roe-MUSCL有限體積懸沙輸運數值模型，其中心思想是使懸沙濃度求解過程同步且具有高分辨率。基于無結構三角形網格，將每一個獨立的三角形單元作為控制個體，分別在各三角形單元上對控制方程進行Gauss線積分離散，采用二階MUS?CL重構格式和Roe型數值通量形式對控制方程進行空間離散，采用一階向前差分格式的時間離散方法，通過對擴散項的平衡處理，建立了能夠適應復雜地形和水動力條件下的高分辨二維有限體積懸沙輸運數學模型，并根據當地海域的懸沙濃度實測資料進行數值驗證。

懸沙輸運方程；有限體積法；無結構網格

數值模擬是了解懸沙輸運特征的主要手段，在近海大范圍開敞水域，水平尺度要遠大于垂向尺度，忽略物理參數在垂向的變化，采用二維懸沙輸運模型進行數值研究是可行的。

上世紀五十年代開始，泥沙數學模型的研究經了歷從一維到二維再到三維，從只考慮單因素的影響到考慮多種因素共同作用，從非耦合到耦合的過程，自上世紀八十年代開始獲得了迅速發展。Chang［1］等對淺水環境下等坡度海灘懸沙輸運特征受波浪、潮流的影響進行了研究。竇國仁等［2］推導得出了波浪和潮流共同作用下的懸沙輸運控制方程和挾沙能力公式，構建了淺水環境下的二維泥沙數學模型。Pruszak［3］等選擇具有代表性的中值粒徑作為研究對象，進行了示蹤沙床面泥沙濃度的實驗研究。Green［4］等對新西蘭海岸波浪以及連續潮流對泥沙運動影響的研究。辛文杰［5］等將波浪概化為時間平均的波浪流分布場，將波浪引起的輻射應力、波流摩阻力以及波流挾沙力耦合到水流運動方程和懸沙輸運方程中，得到了潮流、波浪聯合作用下的二維懸沙數學模型。Eynde［6］等運用拉格朗日法所構建的二維泥沙輸運模型，考慮了潮流和波浪對底部的共同作用，并對比利時海岸疏浚棄土擴散過程進行了數值模擬。

眾多學者長期以來對淺水復雜水域的懸沙輸運模擬研究，弱化了對精度的要求，更多體現在懸沙輸運趨勢方面，而本文采用的懸沙輸運高精度數值方法基于無結構網格和有限體積法原理的高分辨模擬技術研究是合理且必要的。

1 懸沙輸運控制方程

取笛卡爾坐標系，z軸垂直向上，XOY坐標面取未擾動靜止水平面，得到垂向平均的二維懸沙輸運計算控制方程

其中：懸移質泥沙水平擴散項表示為

Biography：SUN Da?wei（1983-），male，lecturer.

式中：u，v代表水深平均的流速在x，y方向的流速分量；H=h+s代表總水深，其中S代表水位；x,y代表平均海平面以下深度；ρx，ρy分別代表懸沙x，y方向的泥沙紊動擴散系數，C代表含沙量；Fs代表床面沖淤函數，表達式為

式中：ω代表泥沙沉降速度；α3代表泥沙沉降幾率，通常取0.67～0.84；ω0代表水流挾沙力；uc代表泥沙起動流速；uf代表泥沙懸浮流速；p代表挾沙力恢復飽和系數；γ代表含沙量恢復飽和系數。

C*代表波、流共同作用下的挾沙力

式中：α=0.023，β=0.000 4；T代表波的周期；Hw代表波高；C2代表謝才系數，根據曼寧公式獲得。

2 控制方程的離散和求解

式（1）是懸沙輸運運動方程的守恒形式，改寫為

令Fz=（m，n），則向量形式（6）簡化為

本文選擇MUSCL—TVD格式對方程進行求解

則所構造的MUSCL-TVD格式水動力計算公式如下

當線性坡度是較小值時，限制函數取Vi+1-Vi；當單元邊界兩側坡度方向相反時，限制函數取0值。

時間項進行向前差分處理，

擴散項進行降階處理的方法

式中：w代表定義在單元形心處的任意函數值

根據格林公式得到

式中：Ω表示的是虛線路徑包圍的面積

式中：xi,yi,wi(i=1,2,3,4)分別代表O1,P,O2,Q的平面坐標和w的值。計算過程中用到的控制體頂點的物理量根據文獻［7］得到。

水體中的溫度、鹽度、懸沙濃度等調整過程較為緩慢，需要給出合理的初始場，初始場可以取自實測資料或均勻場，即

由于計算范圍較大，缺少開邊界潮位觀測資料，本文利用中國海大范圍潮波模型提供水邊界條件，在邊界單元處給定水位過程和懸沙濃度過程［8］

式中：COB（t）代表開邊界懸沙濃度過程，當邊界為入流時由Cin(x,y,t)確定；當邊界為出流時由確定。固邊界單元處采用鏡像法：CR=CL。

3 模型應用與驗證

模型的計算區域為140 km× 180 km，包括北緯31°74′～35°24′，東經119°15′～123°00′，如圖2所示。為適應江蘇北部海岸線不規則的特點，盡量弱化固邊界對計算結果的影響，采用無結構三角形網格單元對計算區域進行剖分，如圖3所示，網格邊長最大尺度約400m，最小尺度約100m。

圖1 積分路線示意圖Fig.1 Sketch of integral path

圖2 水深測站示意Fig.2 Sketch of survey station

圖3 計算域網格劃分Fig.3 Sketch of triangle mesh

根據對此海域波況的觀測統計，以NE、E、ESE向波最為頻繁［9］。模型采用年最大風速下NE向2.4 m有效波高［10］為波浪場計算初始值。

水流模型以淺水方程作為控制方程，離散求解采用與本文懸沙輸運模型相同的二階MUSCL重構格式和Roe型數值通量形式，這里不再重述。

本文采用《江蘇省近海海洋綜合調查與評價專項》所提供的實測資料作為懸沙輸運數學模型的邊界條件和驗證依據。懸沙濃度測站位置如圖2所示。北部海州彎水域懸沙中值粒徑范圍在0.006～0.022 mm之間分布；廢黃河口至小洋口港海域懸沙中值粒徑范圍在0.007～0.01 mm之間分布；洋口港以南海域懸沙中值粒徑范圍在0.009～0.013 mm之間分布［11］。

圖4 大潮漲急平均含沙量Fig.4 Sediment concentration for the maximum flood during spring tide

圖5 大潮落急平均含沙量Fig.5 Sediment concentration for the maximum ebb during spring tide

圖6 大潮懸沙濃度驗證Fig.6 Suspended sediment validation during spring tide

由圖4、5可知，大潮水體含沙量在輻射沙脊群中心處達到最高值，隨著水深的不斷加大而逐漸降低，等值線分布形式與等水深線大體一致（如圖2），高含沙量水體存在于沙脊群水深較淺區域。在漲潮過程中，各方向水流通過沙脊間潮流水道向沙脊群交匯處集中；在落潮過程中，由沙脊群頂點向外呈輻散狀流動。該海域水動力強勁，泥沙易于懸浮，形成了高含沙量水域分布特征。隨著水深的不斷加大，水體的挾沙能力逐漸下降，在此過程中外圍海域水體中的懸移質沙體逐漸落淤，含沙量也隨之下降，形成了與等水深線分布相似的含沙量等值線，這與文獻［8］的分析結果相近。

由于驗證數據有限，故大潮時對7個觀測點進行驗證，懸沙濃度驗證圖6表明，懸沙濃度的模擬值同實測值的分布趨勢基本相同。

輻射沙脊群海域含沙量等值線與水深等值線分布相似，水體含沙濃度隨著水深的減小而變大，在沙脊群中心處達到極大值。漲潮過程中水流由各方向潮流通道匯聚于此，落潮時沿著相反方向分散，所以輻射沙脊群海域的流場具有輻射狀特征。由于泥沙在強勁流場作用下呈懸移質狀態，是造成此區域水體的高含沙量的原因。水體協沙能力隨著水深的增加而變弱，在外圍海域逐漸落淤，形成推移質，故水體含沙量變小，形成了與水深等值線相似的含沙量等值線分布。

4 結論

本文構建了基于無結構三角形網格的Roe-MUSCL有限體積懸沙數值模型，其中心思想是使水動力求解與懸沙濃度求解過程同步且具有高分辨率。基于無結構三角形網格，將每一個獨立的三角形單元作為控制個體，分別在各三角形單元上對控制方程進行Gauss線積分離散，采用二階MUSCL重構格式和Roe型數值通量形式對控制方程進行空間離散，采用一階向前差分格式的時間離散方法，建立了能夠適應復雜地形和水動力條件下的高分辨二維有限體積懸沙數學模型。

對江蘇北部輻射沙脊群海域進行了懸沙數值模擬，通過對計算結果和觀測數據的比較，表明該模型能較好地模擬各測站的懸沙濃度過程，對復雜地形具有很好的適應性。

［1］竇國仁,董鳳舞,竇希萍,等.河口海岸泥沙數學模型研究［J］.中國科學（A輯），1995,25（9）：995-1 001.

［2］Green Malcolm O.,Black Kerry P.And Amos Carl L.Control of estuarine sediment dynamics by interactions between currents and waves at.Several scales［J］.Marine Geology,1997，20（1）：24-32.

［3］辛文杰.潮流、波浪綜合作用下河口二維懸沙數學模型［J］.海洋工程,1997,15（1）：30-47. XIN W J.Numerical Model of 2D Estuarial Suspended Sediment Motion Under the Interaction of Tide Flow and Waves［J］.The Ocean Engineering,1997（1）：30-47.

［4］李紹武，鄭建軍.回流區水流運動二維模擬［J］.港工技術，2006（4）：8-11. LI S W,ZHENG J J.2D Numerical Modeling for Circulation in a Flume with Dead Zone［J］.Port Enginering Technology,2006 （4）：8-11.

［5］楊耀中.南黃海輻射沙脊群懸沙通量數值研究［D］.南京:河海大學,2010.

［6］任美鍔.江蘇海岸帶和海涂資源綜合調查報告［M］.北京：海洋出版社,1986.

［7］王穎.黃海陸架輻射沙脊群［M］.北京：中國環境科學出版社，2002.

［8］楊耀中.南黃海輻射沙脊群懸沙通量數值研究［D］.南京：河海大學，2010.

［9］任美鍔.江蘇海岸帶和海涂資源綜合調查報告［M］.北京：海洋出版社，1986.

［10］王穎.黃海陸架輻射沙脊群［M］.北京：中國環境科學出版社，2002.

A sediment transport numerical model with finite volume method using in the radial sand ridge field of South Yellow Sea

SUN Da?wei1,ZHANG Xin2
（1.Inner Mongolia University for the Nationalities,Tongliao 028000,China;2.Hebei Port and Waterway Engineering Planning and Designing Institute,Tianjin 300074,China）

Suspended sediment transport is an important factor for the evolution of the offshore area.A large amount of suspended sediment is transported under the wave and current action in the radial sand ridge field of the South Yellow Sea.Based on triangular grid,each single triangle element was seemed as the control volume to make Gauss line integral and discrete process.By the format of second order MUSCL Reconstruction and Roe flux,the space discretization of the governing equations was made.With the method of time discretization which based on the format of one?order forward difference and the equilibrium treatments of diffusion term，the two?dimensional water and sediment mathematical model with finite volume and high precision was established for complex terrain and hy?drodynamic condition,and the water flow of the radial sand ridge field in the northern part of Jiangsu province was numerically simulated.

suspended sediment transport equation;finite volume method;unstructured mesh

TV 139.16

A

1005-8443（2016）03-0237-05

2015-08-10；

2015-08-28

內蒙古民族大學博士科研啟動基金資助項目（BS325）

孫大威（1983-），男，遼寧省蓋縣人，講師，主要從事淺水環境水動力數值模擬研究。