肖天葆

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

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

肖天葆1，2

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

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

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

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

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

1 自然條件

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

本海域以風浪為主，常浪為偏東北向，多年出現(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。

連云港海域地處淤泥質海岸，水淺坡緩，海床泥沙中值粒徑0.002～0.004 mm，干密度600～640 kg/m3，近岸多年平均含沙量0.2～0.25 kg/m3。“波浪掀沙，潮流輸沙”是連云港海域泥沙運動的控制機制。

2 基本理論

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

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

對于近岸淺水海域，有

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

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

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

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

3.1 潮流基本方程

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

3.2 波浪基本方程

在直角坐標系下，動譜平衡方程形式如下

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

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

3.3 泥沙基本方程

懸移質輸移擴散方程

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

床面沖淤方程

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

底部沖淤函數(shù)Fs

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

波流共同作用下的床面剪切應力

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

3.4 網格劃分、邊界條件及求解方法

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

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

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

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

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

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

3.5 潮流模型的率定及驗證

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

3.6 波浪模型的率定及驗證

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

表1 大西山測站常風天波浪參數(shù)驗證結果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 海床臨界沖刷切應力的確定

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

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

圖2 技術路線圖Fig.2 Technology roadmap

4.1 多年平均含沙量驗證

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

4.2 含沙量分布驗證

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

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

圖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)與采用遙感技術定量反演出的2005年6月2日中潮某一時刻表層懸沙含量分布圖［11］（詳見圖4）很相似。

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

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

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

4.3 海床沖淤平衡驗證

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

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

4.4 臨界沖刷切應力場

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

5 結論

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

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

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

［1］鄧伯強.淤泥臨界沖刷切應力研究［J］.水動力學研究與進展，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］.水科學進展，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］龐啟秀，白玉川，楊華，等.淤泥質淺灘泥沙臨界起動切應力剖面確定［J］.水科學進展，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］鄭俊.近岸水沙界面通量與挾沙力關系及其應用研究［D］.南京：河海大學，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］解鳴曉，張瑋，張庭榮.淤泥質海岸泥沙運動模擬及進港航道大風天回淤特性研究［J］.應用基礎與工程科學學報，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，海岸與河口潮流泥沙模擬技術規(guī)程［S］.

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

［11］左書華，龐啟秀，楊華，等.海州灣海域懸沙分布特征及運動規(guī)律分析［J］.山東科技大學學報（自然科學版），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.