孫大威，張 鑫

（1.內(nèi)蒙古民族大學(xué)，通遼028000；2河北省水運(yùn)工程規(guī)劃設(shè)計(jì)院，天津300074）

一種有限體積懸沙輸運(yùn)模型在蘇北輻射沙脊群海域的應(yīng)用

孫大威1，張 鑫2

（1.內(nèi)蒙古民族大學(xué)，通遼028000；2河北省水運(yùn)工程規(guī)劃設(shè)計(jì)院，天津300074）

懸沙輸運(yùn)是近岸水域地形演變的重要因素之一，南黃海輻射沙脊群海域懸沙在波浪和潮流的共同作用下有著巨量的輸運(yùn)。文章構(gòu)建了基于無結(jié)構(gòu)三角形網(wǎng)格的Roe-MUSCL有限體積懸沙輸運(yùn)數(shù)值模型，其中心思想是使懸沙濃度求解過程同步且具有高分辨率。基于無結(jié)構(gòu)三角形網(wǎng)格，將每一個(gè)獨(dú)立的三角形單元作為控制個(gè)體，分別在各三角形單元上對(duì)控制方程進(jìn)行Gauss線積分離散，采用二階MUS?CL重構(gòu)格式和Roe型數(shù)值通量形式對(duì)控制方程進(jìn)行空間離散，采用一階向前差分格式的時(shí)間離散方法，通過對(duì)擴(kuò)散項(xiàng)的平衡處理，建立了能夠適應(yīng)復(fù)雜地形和水動(dòng)力條件下的高分辨二維有限體積懸沙輸運(yùn)數(shù)學(xué)模型，并根據(jù)當(dāng)?shù)睾Ｓ虻膽疑碀舛葘?shí)測(cè)資料進(jìn)行數(shù)值驗(yàn)證。

懸沙輸運(yùn)方程；有限體積法；無結(jié)構(gòu)網(wǎng)格

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

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

眾多學(xué)者長(zhǎng)期以來對(duì)淺水復(fù)雜水域的懸沙輸運(yùn)模擬研究，弱化了對(duì)精度的要求，更多體現(xiàn)在懸沙輸運(yùn)趨勢(shì)方面，而本文采用的懸沙輸運(yùn)高精度數(shù)值方法基于無結(jié)構(gòu)網(wǎng)格和有限體積法原理的高分辨模擬技術(shù)研究是合理且必要的。

1 懸沙輸運(yùn)控制方程

取笛卡爾坐標(biāo)系，z軸垂直向上，XOY坐標(biāo)面取未擾動(dòng)靜止水平面，得到垂向平均的二維懸沙輸運(yùn)計(jì)算控制方程

其中：懸移質(zhì)泥沙水平擴(kuò)散項(xiàng)表示為

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

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

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

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

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

2 控制方程的離散和求解

式（1）是懸沙輸運(yùn)運(yùn)動(dòng)方程的守恒形式，改寫為

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

本文選擇MUSCL—TVD格式對(duì)方程進(jìn)行求解

則所構(gòu)造的MUSCL-TVD格式水動(dòng)力計(jì)算公式如下

當(dāng)線性坡度是較小值時(shí)，限制函數(shù)取Vi+1-Vi；當(dāng)單元邊界兩側(cè)坡度方向相反時(shí)，限制函數(shù)取0值。

時(shí)間項(xiàng)進(jìn)行向前差分處理，

擴(kuò)散項(xiàng)進(jìn)行降階處理的方法

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

根據(jù)格林公式得到

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

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

水體中的溫度、鹽度、懸沙濃度等調(diào)整過程較為緩慢，需要給出合理的初始場(chǎng)，初始場(chǎng)可以取自實(shí)測(cè)資料或均勻場(chǎng)，即

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

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

3 模型應(yīng)用與驗(yàn)證

模型的計(jì)算區(qū)域?yàn)?40 km× 180 km，包括北緯31°74′～35°24′，東經(jīng)119°15′～123°00′，如圖2所示。為適應(yīng)江蘇北部海岸線不規(guī)則的特點(diǎn)，盡量弱化固邊界對(duì)計(jì)算結(jié)果的影響，采用無結(jié)構(gòu)三角形網(wǎng)格單元對(duì)計(jì)算區(qū)域進(jìn)行剖分，如圖3所示，網(wǎng)格邊長(zhǎng)最大尺度約400m，最小尺度約100m。

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

圖2 水深測(cè)站示意Fig.2 Sketch of survey station

圖3 計(jì)算域網(wǎng)格劃分Fig.3 Sketch of triangle mesh

根據(jù)對(duì)此海域波況的觀測(cè)統(tǒng)計(jì)，以NE、E、ESE向波最為頻繁［9］。模型采用年最大風(fēng)速下NE向2.4 m有效波高［10］為波浪場(chǎng)計(jì)算初始值。

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

本文采用《江蘇省近海海洋綜合調(diào)查與評(píng)價(jià)專項(xiàng)》所提供的實(shí)測(cè)資料作為懸沙輸運(yùn)數(shù)學(xué)模型的邊界條件和驗(yàn)證依據(jù)。懸沙濃度測(cè)站位置如圖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 大潮懸沙濃度驗(yàn)證Fig.6 Suspended sediment validation during spring tide

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

由于驗(yàn)證數(shù)據(jù)有限，故大潮時(shí)對(duì)7個(gè)觀測(cè)點(diǎn)進(jìn)行驗(yàn)證，懸沙濃度驗(yàn)證圖6表明，懸沙濃度的模擬值同實(shí)測(cè)值的分布趨勢(shì)基本相同。

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

4 結(jié)論

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

對(duì)江蘇北部輻射沙脊群海域進(jìn)行了懸沙數(shù)值模擬，通過對(duì)計(jì)算結(jié)果和觀測(cè)數(shù)據(jù)的比較，表明該模型能較好地模擬各測(cè)站的懸沙濃度過程，對(duì)復(fù)雜地形具有很好的適應(yīng)性。

［1］竇國(guó)仁,董鳳舞,竇希萍,等.河口海岸泥沙數(shù)學(xué)模型研究［J］.中國(guó)科學(xué)（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］辛文杰.潮流、波浪綜合作用下河口二維懸沙數(shù)學(xué)模型［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］李紹武，鄭建軍.回流區(qū)水流運(yùn)動(dòng)二維模擬［J］.港工技術(shù)，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］楊耀中.南黃海輻射沙脊群懸沙通量數(shù)值研究［D］.南京:河海大學(xué),2010.

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

［7］王穎.黃海陸架輻射沙脊群［M］.北京：中國(guó)環(huán)境科學(xué)出版社，2002.

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

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

［10］王穎.黃海陸架輻射沙脊群［M］.北京：中國(guó)環(huán)境科學(xué)出版社，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

內(nèi)蒙古民族大學(xué)博士科研啟動(dòng)基金資助項(xiàng)目（BS325）

孫大威（1983-），男，遼寧省蓋縣人，講師，主要從事淺水環(huán)境水動(dòng)力數(shù)值模擬研究。