數(shù)值模擬中海床臨界沖刷切應(yīng)力的計算方法

肖天葆

（1.廣東茂名濱海新區(qū)管理委員會，茂名525000；2.河海大學(xué)港口海岸與近海工程學(xué)院，南京210098）

數(shù)值模擬中海床臨界沖刷切應(yīng)力的計算方法

肖天葆1，2

（1.廣東茂名濱海新區(qū)管理委員會，茂名525000；2.河海大學(xué)港口海岸與近海工程學(xué)院，南京210098）

針對目前計算海床臨界沖刷切應(yīng)力公式適用性不強的問題，文章從實測資料出發(fā)，借助于波流共同作用下的泥沙數(shù)學(xué)模型，提出了一種經(jīng)過驗證的計算海床臨界沖刷切應(yīng)力的通用性方法。該方法不需要大范圍的海域泥沙實測資料，其適用性相對較強，且能在一定程度上很好的保證海床臨界沖刷切應(yīng)力的連續(xù)性，減少人為隨意的手動調(diào)試。

臨界沖刷切應(yīng)力；海床；數(shù)值模擬

目前，在以切應(yīng)力法處理水沙界面通量的近岸泥沙運動的數(shù)值模擬中，海床臨界沖刷切應(yīng)力的計算是一大難點。雖然人們對此的研究比較多［1-3］，并有計算臨界沖刷切應(yīng)力的相應(yīng)公式，且一般認為其是泥沙粒徑、干密度及水深等參數(shù)的函數(shù)，如張瑞謹公式、Van Rijn公式、Nicholson和O′Connor公式等［4］。但具體到某一特定海域，若實測資料缺乏，再加上前人得出的經(jīng)驗公式一般適用性有限，那么在這一特定海域就不一定適用。因此，提出一種計算海床臨界沖刷切應(yīng)力的通用性方法顯得十分必要。

本文擬根據(jù)連云港歷年實測資料及懸沙遙感定量反演圖，借助于波流共同作用下的泥沙數(shù)學(xué)模型，從基本理論出發(fā)，提出了一種經(jīng)過驗證的計算海床臨界沖刷切應(yīng)力的通用性方法，為泥沙運動的研究提供參考。

1 自然條件

連云港地處江蘇省北部黃海海州灣西南岸，潮汐運動受制于南黃海駐波系統(tǒng)，為非正規(guī)半日潮型，多年平均潮差3.39 m，平均漲潮歷時5 h38 min，平均落潮歷時6 h48 min，由于受近岸地形影響，潮流從外海向岸邊逐漸由逆時針旋轉(zhuǎn)流向往復(fù)流過渡。

本海域以風(fēng)浪為主，常浪為偏東北向，多年出現(xiàn)頻率26.41%；強浪為偏北向，多年出現(xiàn)頻率18.40%。據(jù)1962～2003年實測波浪資料統(tǒng)計，累年有效波高H1/3為0.4～0.6 m，累年平均波高H1/10為0.5 m。

連云港海域地處淤泥質(zhì)海岸，水淺坡緩，海床泥沙中值粒徑0.002～0.004 mm，干密度600～640 kg/m3，近岸多年平均含沙量0.2～0.25 kg/m3?！安ɡ讼粕?，潮流輸沙”是連云港海域泥沙運動的控制機制。

2 基本理論

海床的臨界沖刷切應(yīng)力是泥沙數(shù)學(xué)模型計算的重要參數(shù)，一般認為其與計算海域的泥沙特性有關(guān)。根據(jù)Sanford和Maa［5］的研究，淤泥質(zhì)海床的臨界沖刷切應(yīng)力隨海床泥沙密度垂向的增大而迅速增大，若海床沖刷1 cm，則臨界沖刷切應(yīng)力可增加0.8 N/m2以上。此外，龐啟秀［3］在室內(nèi)水槽和環(huán)形槽試驗中也得出類似結(jié)論，即臨界沖刷切應(yīng)力隨水深及泥沙密度的增大而增大，直至達到某一固定值。另據(jù)盧惠泉［6］的研究：水動力對沉積物的機械分選作用顯著，使沉積物的分布具有較明顯的規(guī)律性。這也就是說水動力機械分選沉積物，進而影響海床臨界沖刷切應(yīng)力，即海床臨界沖刷切應(yīng)力會隨水動力的變化而變化，直至達到某一種平衡，水動力強，其相應(yīng)的臨界沖刷切應(yīng)力也就大，水動力弱，其相應(yīng)的臨界沖刷切應(yīng)力也就小。綜上所述，可以這樣認為：水動力在塑造海床的同時，也塑造了相應(yīng)的臨界沖刷切應(yīng)力，在海床沖淤平衡的條件下，床面切應(yīng)力與臨界沖刷切應(yīng)力正相關(guān)。

由此，提出以下計算海床臨界沖刷切應(yīng)力的公式

對于近岸淺水海域，有

對于外海深水海域，由于波流的掀沙作用已不明顯，因此可以認為該海域海床不可沖刷，即

式中：τce表示海床臨界沖刷切應(yīng)力，是判定泥沙起懸的一個物理量，若τce小于床面切應(yīng)力，則床面泥沙將起懸；τmax、τmean分別表示床面切應(yīng)力的最大值和平均值，其值可以從數(shù)模結(jié)果中提取出來。k值表征海床沖刷的難易程度，應(yīng)通過模型率定給出，其值越大，表示海床越難沖刷。一般來說，從近岸淺水海域到外海深水海域，為保證海床沖刷難易的連續(xù)性，k值應(yīng)逐漸增大。

3 數(shù)學(xué)模型的建立

根據(jù)基本理論的分析，海床的臨界沖刷切應(yīng)力在達到某一種平衡之前是隨水動力的變化而變化的，因此，若要計算出一個不隨時間而變化的臨界沖刷切應(yīng)力場，就需要考慮水動力的平均情況，即選取代表水動力。根據(jù)解鳴曉［7］對連云港海域的研究，代表潮可選取2005年9月的實際中潮，因為該實際潮平均潮差3.82 m，高于多年平均潮差12.6%，符合Latteux［8］對代表潮的研究。由于連云港海域常浪向偏NE向，為簡化模型計算，本次研究僅選取NE向為代表波向，其相應(yīng)的風(fēng)向取45°。

3.1 潮流基本方程

式中：x、y為直角坐標(biāo)系坐標(biāo)；t為時間變量；Y為平均水深；ζ為相對于平均海平面的潮位；Ux、Uy為x、y方向上的垂線平均速度；r為水流密度；g為重力加速度；Nx、Ny為x、y方向的水平紊動粘性系數(shù)；f為科氏參數(shù)（f= 2wsinf，w為地球旋轉(zhuǎn)角速度，f為緯度）；τx、τy為床面剪切應(yīng)力在x、y方向的分量。

3.2 波浪基本方程

在直角坐標(biāo)系下，動譜平衡方程形式如下

式中：N為動譜密度；σ為相對波浪頻率（當(dāng)坐標(biāo)系隨水流運動時觀測到的頻率）；θ為波向，Cx、Cy、Cσ、Cθ為波浪沿x、y、σ、θ方向傳播的速度。S為源匯項，如下表示

式中：Sin為風(fēng)能輸入項；Snl為非線性波-波相互作用的能量傳輸；Sds為波浪白帽耗散造成的能量損失；Sbot為波浪底部摩阻所造成的能量損失；Ssurf為波浪破碎所造成的能量損失。

3.3 泥沙基本方程

懸移質(zhì)輸移擴散方程

式中：S為垂線平均含沙量；Dx、Dy為x、y方向的泥沙擴散系數(shù)；Fs為底部沖淤函數(shù)。

床面沖淤方程

式中：γd為床沙干容重；ηb為海底床面的沖淤變化量。

底部沖淤函數(shù)Fs

式中：τ為水流底部剪切應(yīng)力；τcd為臨界淤積切應(yīng)力；τce為臨界沖刷切應(yīng)力；α為淤積概率；M為沖刷系數(shù)，kg/（m2·s）；ω為泥沙絮凝沉速。

波流共同作用下的床面剪切應(yīng)力

式中：τcw為波流共同作用下的切應(yīng)力；τc為考慮純流作用下的切應(yīng)力；τw為考慮純波浪作用的切應(yīng)力；b、p、q為綜合表達式。

3.4 網(wǎng)格劃分、邊界條件及求解方法

模型范圍北起日照（35°22′30″N，119°33′E），南至廢黃河口附近（34°17′00″N，120°17′E），東西寬約99.7 km，南北長約119.3 km，模型內(nèi)水域面積約8 648 km2。（詳見圖1）由于三角形非結(jié)構(gòu)網(wǎng)格能很好的擬合計算域的復(fù)雜邊界，網(wǎng)格易大易小，適應(yīng)地形好。因此本次研究采用三角形非結(jié)構(gòu)網(wǎng)格對模型區(qū)域進行離散。

潮流模型的西邊界、南邊界為陸邊界，北邊界、東邊界為水邊界。在進行數(shù)值計算時，對于開邊界，由東中國海潮波數(shù)學(xué)模型提供；對于閉邊界，根據(jù)不可入原理取法向流速為0。此外，還采用了干濕判別技術(shù)進行動邊界處理。

波浪模型的陸邊界采用全吸收邊界；灌河采用閉邊界；外海開邊界，通過東中國海大模型計算給出，且通過區(qū)分風(fēng)浪和涌浪的形式給出（風(fēng)浪與涌浪的固定分界周期取8 s）。東中國海大模型的波浪邊界條件對關(guān)注區(qū)域（即連云港海域）影響甚微，可以任意選取，這里取浪透邊界。

泥沙模型的水邊界按實測資料及經(jīng)驗方法給垂線平均含沙量邊界。

在對潮流方程的求解中，采用有限差分的ADI法；在對波浪方程的求解中，地理空間和譜空間離散采用單元中心有限體積法，時間離散采用分布法求解；時間步長Δt=30 s。

圖1 連云港海域地形Fig.1 Topography in Lianyungang sea area

3.5 潮流模型的率定及驗證

采用2005年9月中潮的實測潮位、流速和流向，對模型進行率定和驗證，該潮型可代表多年平均情況，率定出曼寧糙率取0.012～0.016，驗證結(jié)果符合規(guī)范要求［9］，限于篇幅，驗證結(jié)果不列出（具體可參見文獻［10］）。

3.6 波浪模型的率定及驗證

由于實測資料缺乏，連云港海域僅大西山一個波浪測站。鑒于本次研究的目的，僅驗證具有統(tǒng)計意義的常風(fēng)天（NE向）代表波高、波周期及波向。驗證結(jié)果（詳見表1）表明：該模型計算出的波浪場可以作為代表波浪場。

表1 大西山測站常風(fēng)天波浪參數(shù)驗證結(jié)果Tab.1 Verification result of wave parameters of Daxishan station in common windy day

表2 大西山等測站多年平均含沙量驗證Tab.2 Validation of annual average sediment concentration at Daxishan station,etc.

4 海床臨界沖刷切應(yīng)力的確定

本次研究在考慮海床沖淤平衡的條件下，通過對連云港海域?qū)崪y多年平均含沙量及含沙量分布的驗證來確定海床臨界沖刷切應(yīng)力。由于采用的是實測多年平均的資料進行模型驗證，那么相應(yīng)的所計算出來的海床臨界沖刷切應(yīng)力也是多年平均的情況。

在此，給出研究技術(shù)路線：通過經(jīng)過驗證的波流數(shù)學(xué)模型計算出床面切應(yīng)力，進而統(tǒng)計出床面切應(yīng)力的最大值和平均值，再結(jié)合公式（1）、（2），推導(dǎo)出初始臨界沖刷切應(yīng)力τce及相應(yīng)的k值，并把其代入波流共同作用下的泥沙數(shù)學(xué)模型，然后進行多年平均含沙量、含沙量分布及海床沖淤平衡這三方面的反復(fù)驗證，直至模型驗證合理，最終率定出相應(yīng)的臨界沖刷切應(yīng)力（詳見圖2）。

圖2 技術(shù)路線圖Fig.2 Technology roadmap

4.1 多年平均含沙量驗證

從大西山等測站多年平均含沙量驗證結(jié)果（詳見表2）可知：年平均含沙量的計算值與實測值偏差未超過30%，符合規(guī)范要求［9］。

4.2 含沙量分布驗證

從模型計算的年平均含沙量場（詳見圖3），可以看出：

（1）連云港海域含沙量分布表現(xiàn)為“近岸高，外海低”，高含沙水體在近岸運動，尤其在近岸破波帶內(nèi)，且含沙量等值線大致與等深線平行。

圖3 年平均含沙量場Fig.3 Annual average sediment concentration field

圖4 表層懸沙含量分布圖Fig.4 Distribution of surface suspended sediment concentration

（2）存在兩個明顯的高含沙量區(qū)：一個位于灌河口外，年平均含沙量在0.6～0.7 kg/m3；另一個位于臨洪河口至西墅灣一帶，年平均含沙量約0.3～0.5 kg/m3。

（3）模型計算出含沙量場分布形態(tài)與采用遙感技術(shù)定量反演出的2005年6月2日中潮某一時刻表層懸沙含量分布圖［11］（詳見圖4）很相似。

以上對比說明：含沙量分布情況與以往研究成果基本一致。

圖5 海床沖淤平衡驗證Fig.5 Validation of seabed erosion and deposition balance

圖6 海床臨界沖刷切應(yīng)力場Fig.6 Seabed critical shear stress field for erosion

4.3 海床沖淤平衡驗證

海床沖淤平衡作為本次研究的限定性條件，但如果把其看成是模型的驗證資料，那么能在一定程度上解決因?qū)崪y資料缺失而引起的模型驗證不準(zhǔn)確性問題。海床沖淤平衡從理論上解釋就是在一段時間內(nèi)同一地點的海床沖淤相抵。從數(shù)模上表現(xiàn)就是最初的地形經(jīng)過一段時間的沖刷與淤積后所得到的地形與最初的地形相同。因此，通過數(shù)模計算出最初的地形與最終的地形的變化即可驗證海床是否沖淤平衡，其變化越小說明越接近平衡。但在自然條件中，一段時間內(nèi)的海床沖淤變化基本上不可能為0，因此，可以根據(jù)海床的歷年實測沖淤變化或相關(guān)經(jīng)驗來設(shè)定一合理數(shù)值進而衡量海床是否沖淤平衡。

從海床沖淤平衡的驗證結(jié)果（詳見圖5），可知：除部分邊界及島嶼附近，海床沖淤變化年均小于0.1 m/a，說明海床沖淤已達平衡。

4.4 臨界沖刷切應(yīng)力場

根據(jù)對多年平均含沙量、含沙量分布及海床沖淤平衡的驗證，可知所確定的海床臨界沖刷切應(yīng)力場合理可信（詳見圖6）。該計算方法具有如下優(yōu)點：一是不需要大范圍的海域泥沙實測資料，其適用性相對較強。二是能在一定程度上很好的保證海床臨界沖刷切應(yīng)力的連續(xù)性，減少人為隨意的手動調(diào)試。

5 結(jié)論

（1）基于海床垂向物理特性隨沖淤變化劇烈的特征，提出了一種計算海床臨界沖刷切應(yīng)力的新方法：水動力在塑造海床的同時，也塑造了相應(yīng)的臨界沖刷切應(yīng)力，在海床沖淤平衡的條件下，床面切應(yīng)力與臨界沖刷切應(yīng)力正相關(guān)。

（2）計算海床臨界沖刷切應(yīng)力的新方法具有如下優(yōu)點：一是不需要大范圍的海域泥沙實測資料，其適用性相對較強。二是能在一定程度上很好的保證海床臨界沖刷切應(yīng)力的連續(xù)性，減少人為隨意的手動調(diào)試。

（3）在進行泥沙模型驗證時，除采用實測多年平均含沙量及含沙量分布進行驗證外，還提出了采用海床沖淤平衡這一條件作為驗證資料的新理論，該理論能在一定程度上解決因?qū)崪y資料缺失而引起的模型驗證不準(zhǔn)確性問題。

［1］鄧伯強.淤泥臨界沖刷切應(yīng)力研究［J］.水動力學(xué)研究與進展，1991（3）：68-75. DENG B Q.A study for critical shear stress on mud erosion［J］.Journal of Hydrodynamics，1991（3）：68-75.

［2］王虎，劉紅軍，王秀海.考慮滲流力的海床臨界沖刷機理及計算方法［J］.水科學(xué)進展，2014（1）：115-121. WANG H，LIU H J，WANG X H.Mechanism of seabed scour and its critical condition estimation by considering seepage forces ［J］.Advances in Water Science，2014（1）：115-121.

［3］龐啟秀，白玉川，楊華，等.淤泥質(zhì)淺灘泥沙臨界起動切應(yīng)力剖面確定［J］.水科學(xué)進展，2012（2）：249-255. PANG Q X，BAI Y C，YANG H，et al.Critical stress profile for incipient sediment motion on muddy shoals［J］.Advances in Water Science，2012（2）：249-255.

［4］鄭俊.近岸水沙界面通量與挾沙力關(guān)系及其應(yīng)用研究［D］.南京：河海大學(xué)，2012.

［5］Sanford L P，Maa P Y.A unified erosion formulation for fine sediments［J］.Marine Geology，2001，179：9-23.

［6］盧惠泉，蔡鋒，孫全.福建海壇海峽峽道動力地貌研究［J］.臺灣海峽，2009（3）：417-424. LU H Q，CAI F，SUN Q.Study on the channel dynamic geomorphology of Haitan Strait，F(xiàn)ujian［J］.Journal of Oceanography in Tai?wan Strait，2009（3）：417-424.

［7］解鳴曉，張瑋，張庭榮.淤泥質(zhì)海岸泥沙運動模擬及進港航道大風(fēng)天回淤特性研究［J］.應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報，2010（2）：262-272. XIE M X，ZHANG W，ZHANG T R.Numerical modeling of sediment transport on muddy coast and siltation feature in approach channel under the impact of strong wind［J］.Journal of Basic Science and Engineering，2010（2）：262-272.

［8］Latteux B.Techniques for long?term morphological simulation under tidal action［J］.Marine Geology，1995，126：129-141.

［9］JTJ/T 231-2-2010，海岸與河口潮流泥沙模擬技術(shù)規(guī)程［S］.

［10］肖天葆.連云港區(qū)不同等級航道的回淤研究［D］.南京：河海大學(xué)，2015.

［11］左書華，龐啟秀，楊華，等.海州灣海域懸沙分布特征及運動規(guī)律分析［J］.山東科技大學(xué)學(xué)報（自然科學(xué)版），2013（1）：10-17. ZUO S H，PANG Q X，YANG H，et al.Analysis on the distribution and movement of suspended sediment in HaiZhou bay sea area ［J］.Journal of Shandong University of Science and Technology（Natural Science），2013（1）：10-17.

Calculation method of seabed critical shear stress for erosion in numerical simulation

XIAO Tian?bao1,2
（1.Guangdong Maoming Binhai New Area Administrative Committee,Maoming 525000,China;2.College of Harbor,Coastal and Offshore Engineering,Hohai University,Nanjing 210098,China）

For the present problem that the applicability of the formula of calculating the seabed critical shear stress for erosion is not strong,starting from the measured data,a proven general method of calculating the seabed critical shear stress for erosion with the help of the sediment mathematical model under wave and tidal flow was put forward in this paper.This method does not need a wide range of the observed data of marine sediment,and its appli?cability is relatively strong.And to some extent,this method can guarantee the continuity of the seabed critical shear stress for erosion very well and reduce artificial optional manual debugging.

critical shear stress for erosion;seabed;numerical simulation

TV 142；O 242.1

A

1005-8443（2016）03-0231-06

2015-11-02；

2015-12-14

肖天葆（1989-），男，湖南邵陽人，助理工程師，主要從事港口航道工程研究。

Biography：XIAO Tian?bao（1989-），male，assistant engineer.