沙質(zhì)海岸港口工程上下游岸灘演變模擬敏感性參數(shù)研究

王寧舸，宋麗佳

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

0 引言

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

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

1 研究區(qū)域概況

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

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

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

A港海域近岸泥沙中值粒徑為0.23 mm，泥沙運(yùn)動(dòng)以波浪作用下的沿岸輸沙為主要特征，沿岸輸沙率約100萬m3/a，方向總體為自北向南。圖2還給出了工程前水深等深線圖，岸線及水下地形較平直，近岸地形較陡，平均坡度約1/50，泥沙運(yùn)動(dòng)上界高度和下界水深分別約為+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 模型概況

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2.3 模型計(jì)算范圍敏感性分析

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

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

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

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

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

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

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

3 友誼港上下游岸灘演變驗(yàn)證

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

3.1 友誼港概況

友誼港海域的泥沙及水動(dòng)力環(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 岸灘演變驗(yàn)證

3.2.1 模型設(shè)置

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

3.2.2 岸灘沖淤驗(yàn)證

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

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

圖8 岸灘剖面計(jì)算與實(shí)測(cè)對(duì)比Fig.8 Comparison between calculation and field measurement beach profile

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

4 結(jié)語

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

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