王寧舸，宋麗佳

（1.南京水利科學(xué)研究院，江蘇 南京 210029；2.中國船舶重工集團公司第七二四研究所，江蘇 南京 211153）

0 引言

在沙質(zhì)海岸建設(shè)港口工程時，擋沙堤工程會引起上游岸灘淤積，需準(zhǔn)確判斷上游淤積防護年限；下游岸灘沖刷后退，需采取防護措施。自20世紀(jì)70—80年代服務(wù)于毛里塔尼亞友誼港以來，我國已積累了豐富的物理模型試驗研究成果和試驗經(jīng)驗[1-2]。近20年來，隨著泥沙運動基本理論的深入研究和計算機運算能力的不斷提高，數(shù)值模擬成為相對更常規(guī)的手段。基于過程的岸灘長歷時演變數(shù)值模擬難度相對較大，計算耗時較長，數(shù)模研究大多從一線理論或岸灘短期演變角度探討泥沙運動的基本規(guī)律[3]。解鳴曉等[4]采用一線模型對西非塞拉利昂碼頭工程所在海岸的沿岸輸沙能力進行了模擬。張海文等[5]對友誼港北側(cè)沉船附近岸灘模擬了60 d的變化。Tang等[6]研究了離岸堤與直型丁壩對岸灘短歷時沖淤的影響。也有少數(shù)沙質(zhì)海岸岸灘長歷時演變模擬的成果，如Stive等[7]介紹了人工養(yǎng)灘新方式“補沙引擎”，預(yù)測了沙源投放后20 a的岸灘演變趨勢。

沙質(zhì)海岸岸灘長歷時演變模擬仍有許多技術(shù)環(huán)節(jié)尚未明確。本項研究通過數(shù)值模擬，以期探討部分技術(shù)環(huán)節(jié)的處理方式，為沙質(zhì)海岸港口工程相關(guān)數(shù)值模擬研究工作提供一定的參考。

1 研究區(qū)域概況

本項研究以毛里塔尼亞某港口工程（以下簡稱“A港”）為例，工程位于友誼港以南175 km的沙質(zhì)海岸，岸線大致呈12°N—192°N走向，港口工程地理位置示意見圖1。A港口工程建設(shè)挖入式港池，采用防浪擋沙堤掩護進出航道。防浪擋沙堤堤長初步設(shè)計為630 m，對上游來沙進行擋沙防護。下游采用150 m直型丁壩，對下游岸灘進行沖刷防護。A港初步設(shè)計方案見圖2。

圖1 毛里塔尼亞A港工程位置Fig.1 Location of Port A in Mauritania

A港海域面向開敞大西洋，沿岸波浪較強，外海大浪方向一般出現(xiàn)在WNW—NNW向范圍內(nèi)，這3個方向的波浪出現(xiàn)頻率合計約為87%。港區(qū)海域波高Hs在1.0～2.0 m之間的頻率為74%，2.0～3.0 m的波高出現(xiàn)頻率為11%。經(jīng)外海波浪推算以及波能加權(quán)平均統(tǒng)計方法計算等，判斷A港代表波方向與岸線法線方向夾角為19°，H1/10代表波高為1.26 m，代表波周期為10.88 s。A港海域潮汐潮流弱，多年平均潮差不足1.0 m，潮流平均流速不足0.1 m/s[8]。

A港海域近岸泥沙中值粒徑為0.23 mm，泥沙運動以波浪作用下的沿岸輸沙為主要特征，沿岸輸沙率約100萬m3/a，方向總體為自北向南。圖2還給出了工程前水深等深線圖，岸線及水下地形較平直，近岸地形較陡，平均坡度約1/50，泥沙運動上界高度和下界水深分別約為+3 m和-12 m[8]。

圖2 A港工程平面布置與工程前水下地形等深線示意圖Fig.2 Layout of Port A and isobath of underwater terrain before the project

2 沙質(zhì)海岸岸灘演變模擬敏感性分析

2.1 模型概況

本項研究采用Delft3D數(shù)值模擬軟件[9]建立了平面二維沿岸輸沙與岸灘演變數(shù)學(xué)模型。水流模塊基于沿水深平均的淺水方程。波浪模塊采用SWAN模式，基于動譜平衡方程。泥沙計算采用無黏性沙模塊，基于Van Rijn（1993）公式。

數(shù)學(xué)模型采用矩形網(wǎng)格，模型范圍覆蓋全部A港港區(qū)、港區(qū)上下游各10 km及離岸6 km。離岸邊界附近網(wǎng)格步長為50 m，近港區(qū)逐步加密至10 m。模型地形采用現(xiàn)場實測地形。

2.2 岸灘演變模擬敏感性參數(shù)分析

本項數(shù)學(xué)模型采用2016年現(xiàn)場實測潮位和流速資料進行了驗證，并耦合了代表波的傳播變形計算，可模擬近岸破波沿岸流。模型中泥沙中值粒徑設(shè)置為0.23 mm，通過對各項輸沙量計算系數(shù)的率定，最終控制模型上游來沙大致穩(wěn)定在100萬m3/a。

2.2.1 地貌加速因子敏感性分析

地貌加速因子（Morphological scale factor）是一個無量綱的系數(shù)，通過對每個時間步長床面的沖淤量進行該系數(shù)的倍乘，起到加速地形變化、減少模型計算時間的作用，在較長時間序列岸灘演變模擬中尤為重要。地貌加速因子取值太小會影響模型計算效率，但取值過大會導(dǎo)致岸灘變化對水動力的急劇影響，從而影響模型計算精度和穩(wěn)定性。本項研究地貌加速因子考慮4、7和10三種情況。

不同加速因子條件下，地形等深線3 a計算結(jié)果對比見圖3。圖中給出了工程前岸線（+2 m）位置，擋沙堤上游岸線和水下等深線前進，表明岸灘淤積，下游反之表明岸灘沖刷。這3種地貌加速因子條件下的岸灘沖淤演變趨勢一致，沖淤幅度和速率比較接近，大致在3 a后擋沙堤堤頭開始出現(xiàn)繞沙。具體差別來看，對于上游淤積，加速因子采用10時岸灘淤積略快，采用7時岸灘淤積速率居中，采用4時岸灘淤積較前兩者略慢，而下游岸灘沖刷結(jié)果與上游淤積呈相反趨勢。因此，隨著地貌加速因子的增大，岸灘演變計算結(jié)果呈上游淤積越快、下游沖刷越慢的趨勢。

圖3 不同加速因子地形等深線計算結(jié)果對比Fig.3 Comparison of terrain isobath with different acceleration factors

為了說明計算結(jié)果的合理性，采用JTS 145—2015《港口與航道水文規(guī)范》推薦的沙質(zhì)海岸突堤式建筑物上游岸線演變預(yù)報計算方法，估算A港擋沙堤完全有效攔沙時間。擋沙堤堤長630 m，垂直于岸線方向的有效長度約570 m，減去上游岸灘淤積后的剖面寬度約375 m（變形高度15×淤積剖面坡比25），則擋沙堤有效攔沙長度約195 m。再根據(jù)波浪入射角19°和岸灘淤積角約13°，得到A港擋沙堤完全有效攔沙時間約2.9 a，數(shù)模計算結(jié)果（接近3 a）與公式預(yù)報結(jié)果接近。

從計算穩(wěn)定性和精度看，這3種因子取值均合適，并未出現(xiàn)因取值變化而導(dǎo)致計算結(jié)果較大差異。從工程安全和實際應(yīng)用角度，如何選取地貌加速因子值得進一步討論。對于下游岸灘演變，因其還涉及到波浪繞射、下游反向輸沙等復(fù)雜水沙問題，數(shù)學(xué)模型重在反映趨勢，工程實際應(yīng)用還需要結(jié)合物理模型進行綜合分析。而對于上游岸灘演變，當(dāng)發(fā)生較大淤積后，堤頭繞沙將逐步發(fā)生，并產(chǎn)生港池航道回淤問題。因此，從工程安全角度，數(shù)模宜重點考慮上游淤積問題，在保證計算精度的情況下，加速因子可取值偏大一些。本項數(shù)模研究也對加速因子20進行了計算，但計算過程中出現(xiàn)發(fā)散而報錯，因此地貌加速因子選擇10對于工程實際應(yīng)用而言是相對合適的。

綜合上述分析，以A港工程為例，沙質(zhì)海岸岸灘演變模擬中地貌加速因子可初步選擇10。

2.2.2 模型邊界條件敏感性分析

為了防止邊界條件對模型計算結(jié)果的影響，需根據(jù)工況選擇合適的邊界條件或邊界遠(yuǎn)離工程區(qū)域。這涉及到兩個問題，其一是應(yīng)該選擇什么樣的邊界條件，以保證邊界的穩(wěn)定并盡可能符合現(xiàn)場實際情況，這是本小節(jié)討論的重點；其二是邊界應(yīng)該遠(yuǎn)離工程多遠(yuǎn)，既保證工程影響范圍不會觸及邊界，又使得模型范圍不會太大以致降低模型計算效率，這將在下一小節(jié)重點討論。

采用地貌加速因子10，探討不同邊界條件設(shè)置對沙質(zhì)海岸岸灘演變模擬的影響。數(shù)學(xué)模型常用邊界條件如潮位-流速組合邊界條件，除此之外，Delft3D數(shù)值模擬軟件還提供了紐曼邊界（Neumann boundary conditions），可與潮位邊界組合。

紐曼邊界設(shè)置示意見圖4中的AA′和BB′，位于垂直于岸線的兩側(cè)邊界，設(shè)置為水位梯度條件，在風(fēng)場和波浪場較為穩(wěn)定、沿岸潮波變化不大的情況下，水位梯度可概化為0。離岸方向開邊界設(shè)置潮位，見圖4中的AB邊界。紐曼邊界可用以專門解決相對平直海岸的橫向流速、潮位分布問題，防止邊界擾動的發(fā)展。沙質(zhì)海岸多呈順直或順直微彎走向，根據(jù)上述介紹，紐曼邊界類型可適用。

圖4 紐曼邊界設(shè)置示意Fig.4 Neumann boundary setting

分別采用潮位-流速組合和潮位-紐曼組合邊界條件，進行A港工程實施后沿岸輸沙模擬，并統(tǒng)計沿岸輸沙率分布，結(jié)果見表1、表2。表明兩種邊界條件下，距離工程較遠(yuǎn)海域的沿岸輸沙率均能達(dá)到100萬m3/a，并呈現(xiàn)上游沿岸輸沙受阻擋而逐步減小、下游沿岸輸沙逐漸恢復(fù)的趨勢。

表1 防浪擋沙堤上游不同邊界條件模型沿岸輸沙率分布Table 1 Longshore sediment transport rate in the upstream of breakwater under different boundary condition models 萬m3/a

表2 防護丁壩下游不同邊界條件模型沿岸輸沙率分布Table 2 Longshore sediment transport rate in the downstream of groyne under different boundary condition models 萬m3/a

兩種邊界條件的沿岸輸沙率分布差異主要集中在模型上下游邊界處，或者說是100萬m3/a沿岸輸沙率的實現(xiàn)方式不同。潮位-流速組合邊界條件下，模型上下游邊界的沿岸輸沙梯度很大，上游100萬m3/a輸沙率是通過邊界附近泥沙沿程迅速補充得到的，下游則表現(xiàn)為向邊界附近沿岸輸沙急劇衰減，這顯然并不合理。盡管在工程主要影響范圍內(nèi)，潮位-流速組合邊界條件的沿岸輸沙分布尚合理，但在模型邊界附近體現(xiàn)不出現(xiàn)場沿岸輸沙平衡狀態(tài)。相比之下，潮位-紐曼組合邊界條件的處理結(jié)果明顯改善，邊界附近沿岸輸沙大致平衡，能夠反映上游平衡來沙和下游輸沙平衡過境的狀態(tài)。

兩種邊界條件對邊界附近沿岸輸沙梯度的影響最終體現(xiàn)在地形沖淤上，計算3 a后水下地形等深線對比見圖5。潮位-流速組合邊界條件下，由于邊界沿岸輸沙不平衡，導(dǎo)致上、下游邊界水下地形等深線分別呈沖刷后退和淤積前進的趨勢，岸灘不穩(wěn)定。相比之下，潮位-紐曼組合邊界條件下邊界附近等深線總體沖淤穩(wěn)定并保持平順，基本體現(xiàn)了自然岸灘的穩(wěn)定狀態(tài)，更加符合實際。

圖5 不同邊界條件地形計算結(jié)果對比Fig.5 Comparison of terrain isobath with different boundary conditions

據(jù)此可進一步推測，當(dāng)模擬更長歷時的岸灘沖淤變化時，潮位-流速邊界條件引起的邊界不穩(wěn)定會影響模型計算的準(zhǔn)確性，如上游沿岸輸沙需通過邊界岸灘的持續(xù)沖刷進行泥沙供給，但這種供給量是有限的，隨著時間的推移，當(dāng)邊界岸灘沖刷殆盡導(dǎo)致泥沙供給不足后，沙源需求及相應(yīng)的沖刷范圍將會向工程區(qū)延伸，最終影響工程區(qū)附近的計算精度。

綜合上述分析，對于自然狀態(tài)下沿岸輸沙基本平衡、岸灘基本沖淤穩(wěn)定的沙質(zhì)海岸，其岸灘演變模擬需合理選擇邊界條件。采用Delft3D數(shù)值模擬軟件計算時，宜采用潮位-紐曼組合邊界條件。

2.2.3 模型計算范圍敏感性分析

港口工程對沙質(zhì)海岸上下游沿岸輸沙和岸灘沖淤具有一定的影響范圍，因此岸灘演變模擬需要考慮合適的計算范圍。根據(jù)以上敏感性分析結(jié)果，模型地貌加速因子采用10，邊界采用潮位-紐曼組合邊界條件，分別計算模型上、下游各5 km和10 km范圍對沙質(zhì)海岸岸灘演變模擬的影響。

兩種模型計算范圍下，港口工程上下游沿岸輸沙分布見表3、表4。從中可見，在港口工程上下游各10 km計算范圍下，在上下游各6 km之外，沿岸輸沙大致呈不受工程影響的自然平衡狀態(tài)。自6 km位置起向港口方向，沿岸輸沙開始逐步減弱，表明港口工程影響范圍大致達(dá)到上下游各6 km，模型10 km計算范圍是合適的。相比之下，在港口工程上下游各5 km計算范圍下，盡管模型邊界能夠大致達(dá)到100萬m3/a的沿岸輸沙強度，但由于工程實際影響范圍已經(jīng)達(dá)到上下游各6 km，因此模型范圍采用上下游各5 km是偏小的，邊界已經(jīng)受到工程影響。

表3 防浪擋沙堤上游不同計算范圍模型沿岸輸沙率分布Table 3 Longshore sediment transport rate in the upstream of breakwater under different calculation range models 萬m3/a

表4 防護丁壩下游不同計算范圍模型沿岸輸沙率分布Table 4 Longshore sediment transport rate in the downstream of groyne under different calculation range models 萬m3/a

兩種模型范圍下岸灘沖淤演變6 a的計算結(jié)果對比見圖6。從中可以看出，由于港口工程對沿岸輸沙影響范圍大致在上下游各6 km，因此在邊界處輸沙量大致均在100萬m3/a的情況下，模型上下游各5 km計算范圍導(dǎo)致泥沙淤積和沖刷的空間范圍由實際的6 km被強行壓縮至5 km，從而導(dǎo)致上游岸灘淤積前進、下游岸灘沖刷后退的速率偏快。如圖6顯示，模型上下游各5 km范圍條件下，上游水下地形等深線淤積前進幅度略大一些，并且上游泥沙更早越過堤頭并淤積在航道內(nèi)，下游等深線沖刷后退幅度也相對大一些。

圖6 不同計算范圍地形計算結(jié)果對比Fig.6 Comparison of terrain isobath with different model ranges

綜合上述分析，沙質(zhì)海岸岸灘演變模擬需考慮工程對沿岸輸沙和岸灘沖淤的影響范圍。根據(jù)A港海域計算分析，計算范圍可初步采用上下游各10 km。當(dāng)然，范圍超過10 km也是可行的，但考慮到模型的計算效率，10 km已基本滿足需要。

3 友誼港上下游岸灘演變驗證

A港海域?qū)崪y地形資料相對缺乏，為了進一步證實以上敏感性分析結(jié)論的可靠性，建立友誼港海域岸灘演變數(shù)學(xué)模型，采用友誼港實測地形對模型進行驗證。

3.1 友誼港概況

友誼港海域的泥沙及水動力環(huán)境與A港相似，港區(qū)近岸泥沙中值粒徑為0.25 mm，岸灘坡度為1/30，近岸自北向南沿岸輸沙量約100萬m3/a。友誼港-10 m水深的H1/10代表波高為1.10 m，波周期為10.8 s，波浪方向與岸線法線方向的夾角為26°。

3.2 岸灘演變驗證

3.2.1 模型設(shè)置

在A港海域模型的基礎(chǔ)上，將港口工程、水深地形、波浪要素、泥沙粒徑等參數(shù)替換為友誼港海域的相關(guān)參數(shù)，地貌加速因子采用10，模型邊界采用潮位-紐曼組合邊界條件，模型計算范圍為港區(qū)上下游各10 km。

3.2.2 岸灘沖淤驗證

模型以1986年12月友誼港上下游實測地形為初始地形，采用1990年7月實測地形作為驗證地形。1990年7月友誼港上下游水下地形驗證結(jié)果見圖7，上游最大淤積點R1和下游最大沖刷點L2的岸灘剖面計算結(jié)果與實測對比見圖8。結(jié)果表明，模型上下游水下地形等深線分布的計算結(jié)果與實測接近，上游最大淤積點和下游最大沖刷點的剖面變化也與實測吻合良好。

圖7 1990年7月友誼港上下游岸灘地形驗證結(jié)果Fig.7 Validation of upstream and downstream morphology of Friendship Port in July 1990

圖8 岸灘剖面計算與實測對比Fig.8 Comparison between calculation and field measurement beach profile

綜上所述，友誼港海域模型基本復(fù)演了1986—1990年友誼港上下游岸灘沖淤變化，采用經(jīng)敏感性分析確定的各項設(shè)置條件是大體合適的，進行沙質(zhì)海岸岸灘演變模擬時可初步考慮這種設(shè)置方式。

4 結(jié)語

通過敏感性分析，探討了沙質(zhì)海岸岸灘演變模擬中幾個重要的數(shù)學(xué)模型設(shè)置條件。研究認(rèn)為，綜合考慮模型計算精度與計算效率，地貌加速因子可以取10。采用Delft3D數(shù)值模擬軟件進行計算時，宜采用潮位-紐曼組合邊界條件。綜合考慮港口工程對沿岸輸沙、岸灘沖淤的影響范圍以及模型計算效率，計算范圍可以采用港口工程上下游各10 km。

沙質(zhì)海岸岸灘演變模擬的模型影響參數(shù)較多，未來還可探討泥沙擴散系數(shù)、泥沙輸運公式等設(shè)置條件，以及考慮垂向分層的三維模型與平面二維數(shù)學(xué)模型計算結(jié)果的差異。從動力角度，還可對潮位變化、波要素分級設(shè)置等進行研究。通過以上研究工作，以期進一步增強對沙質(zhì)海岸岸灘演變模擬的認(rèn)識，提高精細(xì)化模擬水平。