解鳴曉，李 姍，張 弛，李 鑫，陽志文

（1.南京水利科學研究院水文水資源與水利工程科學國家重點實驗室，南京210029；2.交通運輸部天津水運工程科學研究所工程泥沙交通行業(yè)重點實驗室，天津300456；3.交通運輸部天津水運工程科學研究所港口水工建筑技術(shù)國家工程實驗室，天津300456；4.中交天津港灣工程研究院有限公司中國交建海岸水動力重點實驗室，天津300222；5.河海大學港口海岸與近海工程學院，南京210098）

沙質(zhì)海岸破波帶內(nèi)底部離岸流及沙壩遷移數(shù)值模擬研究

解鳴曉1，2，3，李 姍4，張 弛5，李 鑫1，2，陽志文2，3

（1.南京水利科學研究院水文水資源與水利工程科學國家重點實驗室，南京210029；2.交通運輸部天津水運工程科學研究所工程泥沙交通行業(yè)重點實驗室，天津300456；3.交通運輸部天津水運工程科學研究所港口水工建筑技術(shù)國家工程實驗室，天津300456；4.中交天津港灣工程研究院有限公司中國交建海岸水動力重點實驗室，天津300222；5.河海大學港口海岸與近海工程學院，南京210098）

建立了破波帶內(nèi)三維水動力泥沙數(shù)學模型體系，模型中同時考慮了波浪、波生底部離岸流、非粘性泥沙運動及地貌演變等動力機制。通過波浪水槽實驗實測水、沙及地形數(shù)據(jù)驗證，所建模式可有效模擬破波帶內(nèi)底部離岸流、泥沙運動過程及離岸沙壩的隨時間遷移規(guī)律。模擬成果表明波浪破碎后引起的波生流系在破波帶內(nèi)形成一個時均垂向環(huán)流結(jié)構(gòu)，其中底部時均流向指向外海，表層則指向近岸，靠近破波點處，底部離岸流的流速逐漸降低。在時均垂向環(huán)流體系的輸送下，破波帶內(nèi)被波浪掀起的床面泥沙在底部離岸流作用下持續(xù)向外海輸送，并在破波點附近堆積形成沙壩。沙壩的形成改變了局部波流場結(jié)構(gòu)，在波浪的持續(xù)作用下逐漸向外海運移直至剖面達到基本平衡。

破波帶；三維水沙數(shù)學模型；底部離岸流；泥沙運動；沙壩遷移

近岸沙灘位于波浪破碎帶內(nèi)，泥沙的沿岸及橫向運動直接控制了沙灘地貌格局和岸線的穩(wěn)定性。近10年來，我國風暴潮災害頻率增加，70％以上的沙質(zhì)岸線遭遇明顯侵蝕。目前，數(shù)值模擬技術(shù)已廣泛應用于海岸工程設計中。為進一步實現(xiàn)合理的海岸防護，發(fā)展合適的數(shù)學模型體系預測岸灘剖面演變十分重要，是近年來國內(nèi)外較為熱點的研究課題。沙灘地貌的主要塑造動力為近岸波浪，床面泥沙被波浪掀起后，可在沿岸流、底部離岸流、裂流等波生流系的作用下運移。對沙灘剖面演變和橫向輸沙而言，破波帶內(nèi)的垂向波生流結(jié)構(gòu)，即Stokes漂流、波生時均垂向回流和破波水滾均可起重要作用，動力機制復雜。大量現(xiàn)場經(jīng)驗和水槽實驗數(shù)據(jù)表明，在破波帶內(nèi)流系的輸送下，破波點附近發(fā)育平行岸線的沙壩，并存在一定的遷移性。沙壩的遷移改變了剖面地貌，進而反饋至波浪傳播和波生流系的運移，是一個動態(tài)平衡過程。

在對破波帶內(nèi)水沙動力機制的數(shù)值模式開發(fā)中，歷史研究多數(shù)采用平面二維模型［1-3］，亦有學者發(fā)展了三維模型［4-5］。以往研究在波生流的模擬中多數(shù)采用輻射應力理論，然而輻射應力為平面二維概念，不能刻畫時均剩余動量的垂向分布，以致無法模擬底部離岸流現(xiàn)象。近年來，不同學者亦基于不同理論推導了輻射應力的垂向分布公式［6-10］。基于最新研究進展，Warner等人［11］，Wu和Zhang［12］等人分別將不同表達式嵌入三維水動力模型中，但未嵌入泥沙運動模塊。對于沙壩運動的模擬，尹晶［13］和張弛等人［14］分別基于高階Boussinesq方程和邊界層理論建立了一維模型體系。

筆者［15-16］建立了一個基于時均動力過程的波生流運動三維數(shù)值模式，模型中考慮了輻射應力的垂向分布、破波水滾和波浪紊動效應。本文中基于筆者所建的水動力數(shù)值模式，通過增加非粘性泥沙運動和地貌演變模塊，將其拓展為一個三維水沙動力地貌數(shù)值模式，并采用實驗室內(nèi)沙壩變遷實測數(shù)據(jù)對模型的有效性進行檢驗，并進行相關(guān)討論。

1三維水沙數(shù)值模擬建立

1.1水動力模式

水動力方程的基本形式為雷諾方程，見式（1）～（4）。在對波生流的模擬中，引入波浪時均剩余動量的垂向分布、破波水滾和波浪紊動作為驅(qū)動力。模型網(wǎng)格剖分中采用水平向矩形、垂向貼體坐標的形式。

式中：t為時間；x和y為水平向坐標；η為水位；U和V分別為x和y方向的水平流速分量；ω為σ坐標系下垂向流速；D為水深；M為時均波生剩余動量，可由Lin?Zhang公式求解［7］，表達式見式（5），R為破波水滾動量；KMc和AMc分別為波流共同作用下的垂向和水平紊動系數(shù)；ρw為海水密度。

式中：E為波能；n為波能傳遞率；k為波數(shù)；T為波周期；δ為Kronecker張量標記；i和j分別代表x和y方向。波流共同作用下的底部剪切力τcw采用Soulsby等人［17］的理論。

在波流共同作用的紊動描述中，分別單獨求解水流與波浪引起的紊動系數(shù)AM與KM，并將其疊加［15-16］，水流引起的水平紊動系數(shù)AM（σ）采用Smagorinsky方程求解，形式見式（6），其中Δx和Δy為x和y方向的空間步長；Cs為經(jīng)驗系數(shù)，取值在0.1～0.2左右。波浪引起的水平紊動系數(shù)采用筆者［15-16］推導的公式，見式（7）中所示，其中λ為無因次系數(shù)，取值在0.2～0.5左右［15］。水流引起的垂向紊動系數(shù)KM（σ）采用Mellor?Yamada紊流閉合方程求解，波浪引起的垂向紊動系數(shù)表達式見式（8），其中b為無因次系數(shù)，取值在0.001左右［15］。

筆者［15］推導了波浪破碎引起的水滾能量傳輸方程，見式（9），式中α為水滾能量傳遞率，取值在0.0～1.0之間；AR為水滾面積；ρR為水滾密度；C為波速；KR＝0.375（0.3+2.4s），s為底坡坡度。當波浪參數(shù)確定后，取破波點外AR＝0，然后自破波點向岸迭代求解。水滾能量垂向分布采用Haas和Warner［18］建議形式。

方程求解技術(shù)采取內(nèi)外模式分裂法，平面求解采用顯式差分，垂向求解采用隱式差分。網(wǎng)格配置采用Arakawa-C網(wǎng)格，時均剩余動量和水滾動量項布置在網(wǎng)格中心。動邊界處理采用Oey［19］提出的OGCM法。

1.2波浪模型

波浪模型采用折繞射聯(lián)合REF/DIF數(shù)值模式，模型基于拋物緩坡方程理論，可模擬波浪的淺水變形、折射、繞射和破碎等現(xiàn)象，同時可考慮水流對波浪的反饋影響。

1.3非粘性泥沙運動模型

懸沙運動采用三維對流擴散方程，見式（10），其中C為含沙量；DH、DV分別為泥沙水平及垂向紊動擴散系數(shù)；ω為泥沙沉速，由張瑞瑾公式計算，見式（11）；D50為泥沙中值粒徑；ν為水體運動粘性系數(shù)；ρs為泥沙顆粒密度。泥沙懸浮在水體中，將會改變海水密度，從而引起密度斜壓效應，見式（12），其中Cvol為體積含沙量，模擬時根據(jù)每一時步所求解的水體含沙量數(shù)值修正含沙水體密度，并代回水動力方程中求解。懸浮泥沙與床面的交換形式采用“參考高度”概念，參考點含沙量采用式（13）～（15）求解［4］，其中ks為Nikuradse糙率；在懸沙與床面的交換項處理中，采用式（16）～（17），其中kmx代表垂向網(wǎng)格中高于參考高度a的第一層網(wǎng)格中心位置。在低于參考高度至床面的區(qū)間，認為泥沙以推移質(zhì)形式運動，輸沙率表達式見式（18）～（22），其中q=qbc+qbw為總輸沙率；u*′為摩阻流速；ub為底層網(wǎng)格中心處的流速數(shù)值，vr為垂線平均流速；Uw為波浪底部最大質(zhì)點流速。床面變形方程見式（23），其中zbed為床面高程。在式（17）中的波浪輸沙率計算中，同時考慮了非線性波浪水質(zhì)點運動不對稱的影響［4］。

2沙壩變遷模擬及驗證

為驗證模型的有效性，采用Arcilla等人［20］于Delft水工所開展的LIP11D大型波浪水槽實驗作為算例。實驗所用水槽的長度、寬度和深度分別為233 m、5 m和7 m。在實驗1 b中，坡腳前波高達1.5 m，波周期5.0 m、床面泥沙中值粒徑0.22 mm，接近天然沙灘前緣的波浪環(huán)境，實驗中持續(xù)17 h的窄譜不規(guī)則波作用于原始沙壩地形上，并同步測量了多個站點不同時刻的波高、增減水、時均流速、含沙量和地形變化等關(guān)鍵參數(shù)。圖1中示意了實驗地形和觀測站位置信息，表1中給出了數(shù)學模型中所采用的主要參數(shù)情況。

為量化比較模擬值與實測值間的誤差情況，采用2個誤差指標進行分析，分別為均方差RMS和相關(guān)系數(shù)COR，其中前者反映模擬值和實測值間的偏差，后者則反映兩套數(shù)值的相關(guān)度。

圖2和圖3中分別給出了T=1 h，T=9 h和T=17 h時刻的模擬與實測波高、波浪增減水的對比情況。基于對比結(jié)果，模型可較好地反映出不同時刻波浪沿水槽的傳播規(guī)律和增減水沿程分布情況，模擬誤差均很小?；谀M結(jié)果，波浪在實驗的沙壩地形上存在3次明顯的破碎過程，分別發(fā)生在（1）x≈50 m處，此處波浪初傳播至陡坡上方，在淺水變形作用下發(fā)生破碎；（2）x≈140 m附近，此處位于外側(cè)沙壩頂端；（3）x≈160 m附近，此處位于內(nèi)側(cè)沙壩頂端。在內(nèi)、外兩個初始沙壩間的溝槽內(nèi)（140 m＜x＜160 m），波高略有增大。至于增減水的沿程分布，則呈現(xiàn)出破波區(qū)內(nèi)增水、淺水變形區(qū)內(nèi)減水的整體趨勢。

圖1 LIP11D 1b實驗地形及測站位置示意Fig.1BathymetryandobservationsectionlocationsoftheLIP11D1bcase

表1 LIP11D 1b實驗數(shù)值模擬采用主要參數(shù)情況Tab.1 Input parameters in the numerical model for LIP11D experiment case 1b

圖2 LIP11D 1b實驗不同時刻實測與模擬波高分布對比Fig.2 Comparisons between the simulated and observed wave height for LIP11D 1b case

圖3 LIP11D 1b實驗不同時刻實測與模擬波浪增減水分布對比Fig.3 Comparisons between the simulated and observed wave setup for LIP11D 1b case

圖4中分別給出了不同觀測時刻、不同斷面位置處的模擬與實測時均流速垂向分布的對比情況。經(jīng)對比，在外側(cè)沙壩頂端附近（130 m＜x＜150 m），模型計算所得時均流速數(shù)值較實測值略低。經(jīng)分析認為，這一誤差的原因在于所模擬床面為動床狀態(tài)，特別是在波浪和時均流速的作用下外側(cè)沙壩有明顯遷移，從而導致模擬所得地形和實測地形存在一定差異，進而影響了時均流速的計算精度。盡管如此，總體來說所建模型可較為有效的反映出底部離岸流的垂向分布趨勢，且隨時間的變化規(guī)律亦與實測值較為接近。

從時均流速的垂向分布趨勢來看，最大離岸流速出現(xiàn)在近床面附近，并向表層逐漸衰減；從沿程分布來看，最大底部離岸流出現(xiàn)在沙壩附近區(qū)域。至表層水體，離岸流逐漸反向，流向轉(zhuǎn)為指向岸線。這一現(xiàn)象反映出在破波帶內(nèi)存在一個時均垂向回流結(jié)構(gòu)，其表層向岸、底層離岸，是沙壩向海側(cè)運移的主要驅(qū)動力。

圖5中分別給出了不同觀測時刻、不同斷面位置處的模擬與實測時均含沙量垂向分布的對比情況。與圖4中所反映出的規(guī)律類似，模擬結(jié)果表明在地形變化較為平緩的離岸區(qū)內(nèi)（60 m＜X＜130 m），模擬值與實測值的吻合度較佳，最大誤差出現(xiàn)在外側(cè)沙壩附近（130 m＜X＜140）。

經(jīng)分析，這一誤差的成因主要有三：（1）根據(jù)圖4中結(jié)果所示，在該區(qū)域的模擬所得底部離岸流速小于實測值，導致掀沙能力較實際略有不足；（2）實驗中的實際地形形態(tài)存在沙紋現(xiàn)象，即床面并非嚴格平整，而這些沙紋將引起床面局部的小尺度渦旋，但由于本文中所采用的模擬理論為基于時均模型，難以反映一個亞波周期內(nèi)的小尺度紊動結(jié)構(gòu)；（3）在模型中采用的底部摩阻高度取為常值，而在實際水槽中的沙紋結(jié)構(gòu)并非沿程均勻，導致了一定的不匹配。實際上，由于模型計算中將動力過程和地貌過程耦合，從而地形模擬的精度誤差亦可反饋至水沙運動信息，亦將引起一定的模擬誤差。

圖4 LIP11D 1b實驗不同時刻實測與模擬時均流速分布對比Fig.4 Comparisons between the simulated and observed undertow velocity for LIP11D 1b case

圖5 LIP11D 1b實驗不同時刻實測與模擬時均含沙量分布對比Fig.5 Comparisons between the simulated and observed suspended SSC for LIP11D 1b case

圖6中給出了T=17 h時刻的模擬與實測輸沙率的對比情況?？傮w來說，模型可較好地刻畫輸沙率的沿程分布規(guī)律和量級。在沙壩的外海側(cè)，泥沙輸移方向為指向近岸，而在沙壩的向岸側(cè)，泥沙輸移方向則為指向外海，最大輸沙率出現(xiàn)在沙壩頂端的波浪強破碎區(qū)。圖7中給出了T=17 h時刻的模擬與實測沙灘剖面地形的對比情況。經(jīng)統(tǒng)計，模擬所得的沙壩高度略低于實測高度，且沙壩間的溝槽深度略大于實測深度，誤差的成因與模擬所得底部離岸流強度和含沙量的誤差所致。

盡管如此，綜合以上對不同參數(shù)的誤差分析，從整體上來看，模擬所得的底部離岸流流速、含沙量、輸沙率和剖面地形與實測值均達到了較佳的相關(guān)度，在絕大多數(shù)斷面處，相關(guān)系數(shù)COR均超過了80％，且均方差RMS不大。因此，認為所建模型可有效地復演破波帶內(nèi)沙灘剖面在強波浪作用下的水動力特征、泥沙輸移規(guī)律和沙壩遷移、演變特征，模型精度較高。

需指出的是，在沙壩的橫向遷移中，同時可受指向岸側(cè)的Stokes漂流和指向海側(cè)的底部離岸流雙重影響。在本文采用的算例中，重點討論了沙壩在強底部離岸流條件下的向海遷移。然而，在弱浪環(huán)境下，由于波浪破碎強度低，底部離岸流速較弱，而Stokes漂流占優(yōu)，可能導致沙壩的向岸側(cè)遷移，因此在后續(xù)工作中，應進一步針對這一現(xiàn)象進行相關(guān)研究工作。

圖6 LIP11D1b實驗T=17h時刻實測與模擬輸沙率沿程分布對比Fig.6Comparisonsbetweenthesimulatedandobservedsediment transportrateforLIP11D1bcase(T=17h)

圖7 LIP11D 1b實驗T=17 h時刻實測與模擬地形沿程分布對比Fig.7 Comparisons between the simulated and observed bed elevation for LIP11D 1b case(T=17 h)

3結(jié)論

本文建立了一個基于動力過程的破波帶內(nèi)水動力及沙灘剖面地貌演變的三維數(shù)值模式，模型中同時考慮了波生流的垂向分布、懸沙和底沙運動以及剖面地形的耦合反饋。采用基于大型波浪水槽的實驗數(shù)據(jù)對模型進行了詳細驗證，并評價了模型表現(xiàn)力。研究結(jié)果表明所建數(shù)值模式可較為有效地復演破波帶內(nèi)的水沙動力機制和沙壩發(fā)育和遷移規(guī)律，模擬精度較高。主要結(jié)論如下：

（1）波浪破碎后，在Stokes漂流、波生時均剩余動量和破波水滾動量傳輸?shù)木C合作用下，在破波點岸側(cè)形成一個垂向時均環(huán)流結(jié)構(gòu)，即底部時均流速指向外海，而表層時均流速則指向岸側(cè)，最大底部離岸流速出現(xiàn)在破波點內(nèi)側(cè)。

（2）在破波帶內(nèi)底部離岸流的輸送下，波浪起動的沙灘泥沙可以懸沙和底沙的形式持續(xù)向海側(cè)輸送。由于底部離岸流速向海側(cè)逐漸降低，從而泥沙在破波點附近堆積形成沙壩。沙壩的發(fā)育改變了沙灘剖面地形，進一步改變了波浪傳播和底部離岸流的分布格局，進而使得沙壩有向海遷移的趨勢，直至動力環(huán)境和地貌形態(tài)達到基本穩(wěn)定狀態(tài)。

［1］Nam P T,Larson M,Hanson H,et al.A numerical model of beach morphological evolution due to waves and currents in the vicini?ty of coastal structures［J］.Coastal Engineering,2011,58:863-876.

［2］Nam P T,Larson M.Model of nearshore waves and wave?Induced currents around a detached breakwater［J］.Journal of Waterway, Port,Coastal and Ocean Engineering,2010,6:156-176.

［3］Ding Y,Wang S S Y.Development and application of coastal and estuarine morphological process modeling system［J］.Journal of Coastal Research,Special Issue,2008,52:127-140.

［4］Lesser G R,Roelvink J A,van Kester J A T M,et al.Development and validation of a three?dimensional morphological model［J］. Coastal Engineering,2004,51:883-915.

［5］Roelvink J A.Coastal morphodynamic evolution techniques［J］.Coastal Engineering,2006,53:277-287.

［6］Ardhuin F,Rascle N,Belibassakis K A.Explicit wave?averaged primitive equations using a generalized Lagrangian mean［J］. Ocean Modeling,2008,20:35-60.

［7］Lin P Z,Zhang D.The depth?dependent radiation stresses and their effect on coastal currents［C］//ICHD.Proceedings of the 6th In?ternational Conference of Hydrodynamics:Hydrodynamics VI Theory and Applications,2004:247-253.

［8］Mellor G L.Some consequences of the three dimensional current and surface wave equations［J］.Journal of Physical Oceanogra?phy,2005,33:1 978-1 989.

［9］McWilliams J C,Restrepo J M,Lane E M.An asymptotic theory for the interaction of waves and currents in coastal waters［J］.Jour?nal of Fluid Mechanics,2004,511:135-178.

［10］Xia H Y,Xia Z W,Zhu L S.Vertical variation in radiation stress and wave?induced current［J］.Coastal Engineering,2004,51: 309-321.

［11］Warner J C,Sherwood C R,Signell R P,et al.Development of a three?dimensional,regional,coupled wave,current,and sediment?transport model［J］.Computers&Geosciences,2008,34:1284-1306.

［12］Wu X Z,Zhang Q H.A three?dimensional nearshore hydrodynamic model with depth?dependent radiation stresses［J］.China Ocean Engineering,2009,23:291-302.

［13］尹晶.海岸沙壩運動的實驗與數(shù)值模擬研究［D］.大連:大連理工大學,2012.

［14］張弛,鄭金海,王義剛.波浪作用下沙壩剖面形成過程的數(shù)值模擬［J］.水科學進展,2012,23（1）:104-109. ZHANG C,ZHENG J H,WANG Y G.Numerical simulation of wave?induced sandbar formation［J］.Advances in Water Science, 2012,23（1）:104-109.

［15］Xie M X.Establishment,validation and discussions of a three dimensional wave?induced current model［J］.Ocean Modelling, 2011,38:230-243.

［16］Xie M X.Three?dimensional numerical modeling of the wave?induced rip currents under irregular bathymetry［J］.Journal of Hy?drodynamics,2012,6:864-872.

［17］Soulsby R L,Hamm L,Klopman G,et al.Wave?current interaction within and outside the bottom boundary layer［J］.Coastal En?gineering,1993,21:41-69.

［18］Haas K A,Warner J C.Comparing a quasi?3D to a full 3D nearshore circulation model:SHORECIRC and ROMS［J］.Ocean Mod?eling,2009,26:91-103.

［19］Oey L Y.An OCGM with movable land?sea boundaries［J］.Ocean Modelling,2006,13:176-195.

［20］Arcilla A S,Roevink J A,O'Connor B A.The Delta Flume'93 experiments［C］//ASCE.Proceedings of Coastal Dynamics,2010: 488-502.

Numerical modeling of the undertow and sandbar migration process in the surfzone

XIE Ming?xiao1,2,3,LI Shan4,ZHANG Chi5,LI Xin1,2,YANG Zhi?wen2,3
(1.State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering,Nanjing Hydraulic Research Institute,Nanjing 210029,China;2.Key Laboratory of Engineering Sediment,Tianjin Research Institute for Water Transport Engineering,M.O.T.,Tianjin 300456,China;3.National Engineering Laboratory for Port Hydraulic Construction Technology,Tianjin Research Institute for Water Transport Engineering,M.O.T.,Tianjin 300456,China; 4.Key Laboratory of Coastal Engineering Hydrodynamics of CCCC,Tianjin Port Engineering Institute,CCCC, Tianjin 300222,China;5.College of Harbor,Coastal and Offshore Engineering,Hohai University,Nanjing 210098, China)

A process?based three?dimensional numerical model for surfzone hydrodynamics and beach profile evolution was developed.The model incorporates the coupling process of waves,wave?induced undertow,sediment transport and morphology response.Using the established model,the wave?induced undertow structure and the sand?bar migration characteristics were modeled,and a series experimental datasets were applied to validate the model. The results show that the model can satisfactorily describe the surfzone hydrodynamic and the sandbar migration characteristics in the surfzone with acceptable precision.Once the wave breaks,the vertical distribution of the un?dertow indicates a phase?averaged circulation structure.Outside the breaking zone,the undertow rapidly diminish?es,and the undertow direction points to the offshore near the water surface.For the sediment transport characteris?tics along the beach profile,the sediments entrained by waves continuously move to the breaking point,and deposit near the breaking point to form a sandbar where the undertow velocity decreases.

surfzone;process?based three?dimensional numerical model;wave?induced undertow;sediment transport;sandbar evolution

TV 142；O 242.1

A

1005-8443（2016）04-0349-07

2016-01-27；

2016-02-25

國家自然科學基金青年科學基金（41306033）；交通運輸部應用基礎研究項目（2014329224330）；水文水資源與水利工程科學國家重點實驗室開放項目（2014492211）；天津市自然科學基金青年項目（16JCQNJC06900）

解鳴曉（1982-），男，山東省青島人，副研究員，博士，主要從事海岸動力學研究。

Biography：XIE Ming?xiao(1982-),male,associate professor.