港珠澳大橋人工島水下灘槽演變的數(shù)值模擬與工程檢驗(yàn)

閆 禹，何 杰，徐貝貝，謝至正

(1.港珠澳大橋管理局，廣東 珠海 519060; 2.南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029; 3.廣東港珠澳大橋材料腐蝕與工程安全國家野外科學(xué)觀測(cè)研究站，廣東 珠海 519060)

港珠澳大橋跨越珠江口伶仃洋海域，是連接香港特別行政區(qū)、廣東省珠海市、澳門特別行政區(qū)的大型跨海通道。大橋采用“橋—島—隧”相結(jié)合的方式跨越伶仃洋，工程項(xiàng)目主體之一為主通航區(qū)所在的東、西兩個(gè)大型人工島工程。伶仃洋作為珠江的重要入?？冢碌匦谓倌陙硎冀K保持著“三灘兩槽”的分布格局。而在伶仃洋海域中間建設(shè)大型人工島必將會(huì)改變?nèi)斯u水域的水動(dòng)力條件，進(jìn)而改變?nèi)斯u附近水域的水下地形分布形態(tài)，并在潮汐、波浪等多種動(dòng)力的長(zhǎng)期作用下形成新的水下灘槽格局。作為廣州港重要出海航道的伶仃航道，從東、西人工島中間水域穿過。大型人工島工程的建設(shè)，是否會(huì)對(duì)伶仃航道的水沙環(huán)境和水下地形沖淤形態(tài)產(chǎn)生影響，對(duì)航道今后的發(fā)展具有重要意義。因此，開展大型人工島工程對(duì)灘槽演變的影響研究是十分必要的。

關(guān)于相關(guān)海域人工島對(duì)周邊水沙環(huán)境的影響預(yù)測(cè)已開展了大量研究[1-12]。人工島的存在，改變了局部海域的水下地形分布形態(tài)，同時(shí)改變了局部水域的水動(dòng)力結(jié)構(gòu)，繼而人工島及其周邊的地貌演變產(chǎn)生相應(yīng)的變化。在預(yù)測(cè)海中人工島對(duì)周邊環(huán)境的影響研究中，數(shù)學(xué)模型是一種被廣泛采用的技術(shù)手段[5-10, 13]。例如，盛天航等[9]使用平面二維數(shù)值模型模擬了秦皇島湯河河口人工島建設(shè)后的流場(chǎng)，表明人工島在實(shí)施河道清淤的情況下，有利于泥沙的沖刷；李松喆[7]通過潮流、波浪數(shù)值模型研究了海南省紅塘灣海域人工島的平面布置型式與岸灘演變之間的制衡關(guān)系；何杰等[13-15]使用平面二維數(shù)值模型，探究了港珠澳大橋人工島建設(shè)對(duì)珠江口水動(dòng)力的影響程度，對(duì)沉管隧道基槽的泥沙回淤和東人工島島隧結(jié)合部水動(dòng)力條件變化進(jìn)行了預(yù)測(cè)；陳亮鴻等[4]采用COAWST模型探究海南鋪前灣人工島建設(shè)對(duì)岬角渦旋及海灣地形沖淤的影響。在人工島工程設(shè)計(jì)階段，數(shù)值模型作為預(yù)測(cè)工具，往往采用水動(dòng)力過程[4, 10, 13]、實(shí)測(cè)含沙量[9]、工程前的地形演變過程[9]等來進(jìn)行模型的沖淤驗(yàn)證。模型預(yù)測(cè)的結(jié)果很少能與工程建設(shè)后實(shí)際產(chǎn)生的影響效應(yīng)進(jìn)行直接對(duì)比。因此通過跟蹤人工島工程實(shí)際產(chǎn)生的工程影響效應(yīng)并對(duì)當(dāng)初設(shè)計(jì)階段預(yù)測(cè)的地形演變結(jié)果進(jìn)行檢驗(yàn)將是一項(xiàng)十分有意義的研究。

在工程設(shè)計(jì)階段采用平面二維潮流—懸沙數(shù)學(xué)模型探究了港珠澳大橋人工島對(duì)工程周邊灘槽的演變影響。在模型回淤驗(yàn)證時(shí)對(duì)人工試挖槽的回淤過程進(jìn)行了很好地驗(yàn)證，對(duì)人工島周邊灘槽的演變趨勢(shì)進(jìn)行了預(yù)測(cè)。人工島建成10年后，根據(jù)人工島附近水下地形的監(jiān)測(cè)結(jié)果，對(duì)比人工島工程設(shè)計(jì)階段預(yù)測(cè)的地形與建設(shè)后的實(shí)測(cè)地形差異，分析人工島水下灘槽演變的工程效應(yīng)，同時(shí)為港珠澳大橋工程設(shè)計(jì)階段所采用數(shù)學(xué)模型的預(yù)測(cè)精度和適用性進(jìn)行了佐證。

1 伶仃洋海域水沙環(huán)境

伶仃洋是珠江口最大也是最重要的河口灣，有珠江口東四口門(虎門、蕉門、洪奇瀝和橫門)注入，灣型呈喇叭狀，走向接近NNW-SSE方向，灣頂寬約4 km(虎門口)，灣口寬約30 km(澳門至香港大濠島之間)，縱向長(zhǎng)達(dá)72 km，水域面積約為2 110 km2。伶仃洋水下地形具有西部淺、東部深的橫向分布和灣頂窄深、灣腰寬淺、灣口寬深的縱向分布特點(diǎn)，水下地形呈“三灘兩槽”的基本格局。潮汐類型屬不規(guī)則半日混合潮型，潮差較小，平均潮差為0.86～1.69 m，最大潮差為2.29～3.36 m。潮差由東向西逐漸遞減，由灣口向?yàn)稠斨饾u遞增。伶仃洋潮流流向基本上為較穩(wěn)定的南北向往復(fù)流，東、西部動(dòng)力差異較為明顯。伶仃洋的懸移質(zhì)含沙量具有深槽小、淺灘大，東部低、西部高、枯季清、汛期渾等主要分布特征，多年平均含沙量介于0.1～0.2 kg/m3之間。懸沙運(yùn)動(dòng)主要由灣內(nèi)向外海輸送，凈輸沙量平均為570 kg/m3·d。伶仃洋河床質(zhì)的粒徑具有中灘粗、邊灘細(xì)，灣頂附近較粗、灣口一帶較細(xì)的分布特點(diǎn)，航槽和下游深水區(qū)的床沙中值粒徑一般為0.005～0.010 mm，西灘和東灘的底質(zhì)中值粒徑大部分在0.010 mm以下，但在蕉門口、橫門口以及交椅沙、公沙、攔江沙等處河床質(zhì)明顯粗化，最大中值粒徑為0.1～0.5 mm。

2 港珠澳大橋工程概況

港珠澳大橋東、西兩個(gè)大型人工島分別處在伶仃洋大濠深槽兩側(cè)，實(shí)現(xiàn)了港珠澳大橋橋隧結(jié)構(gòu)的轉(zhuǎn)換。東、西兩個(gè)人工島均于2009年建成。兩人工島形狀一致，呈現(xiàn)鵝卵石形狀，如圖1中所示，西人工島為東西走向，東人工島則呈斜向布置，大頭側(cè)指向橋梁，小頭側(cè)指向隧道，迎水面長(zhǎng)度均為625 m。西島外形略鼓，順?biāo)较虼箢^端180 m、小頭端100 m，面積93 688 m2；東島外形略扁，順?biāo)较虼箢^端225 m、小頭端115 m，面積稍大，為102 462 m2。東、西人工島間距為5 573 m，人工島岸壁按1∶2設(shè)置護(hù)坡。

圖1 東西人工島平面示意

3 平面二維潮流懸沙數(shù)學(xué)模型

3.1 控制方程

研究水域?yàn)槊娣e較大的開闊水域，其平面尺度遠(yuǎn)大于垂直尺度，故可忽略垂向上的物理過程，使用平面二維模型。平面二維潮流數(shù)值模型的水流控制方程為：

(1)

式中：U=(H,Hu,Hv,Hs)T，H為全水深，H=h+η，h為水平面以下水深，η為水面波動(dòng)。

E=(F,G)

(2)

其中，u，v表示x，y方向的流速；s表示水體含沙量。

水流運(yùn)動(dòng)方程中的紊動(dòng)擴(kuò)散項(xiàng)可表示為：

Ed=(Fd,Gd)

(3)

源項(xiàng)S可表示如下：

(4)

其中，S0x、S0y分別是x、y方向的傾斜效應(yīng)項(xiàng)即河床底部高程變化，S0x=-?zb/x，S0y=-?zb/y，zb為河床底面高程。Sfx，Sfy是在x、y方向的底摩擦效應(yīng)項(xiàng)，可表示為：

其中，n為曼寧系數(shù)，f為柯氏系數(shù)，f=2ωsinφ,ω表示地轉(zhuǎn)速度，φ為當(dāng)?shù)氐乩砭暥取?/p>

Fs為床面沖淤函數(shù)，可用下式表示：

Fs=-αω(β1·s*-β2·s)

(5)

式中：α為泥沙的沉降機(jī)率；ω為泥沙沉速；s*為水流挾沙率；uc為泥沙起動(dòng)流速；uf為泥沙懸浮流速。

3.2 數(shù)值求解

采用三角形單元對(duì)計(jì)算區(qū)域進(jìn)行離散，物理變量放置在每個(gè)單元的中心。對(duì)式(1)在計(jì)算域內(nèi)積分，利用Green公式將面積分化為線積分，可得：

(6)

對(duì)第i個(gè)單元積分結(jié)果可表示為：

(7)

其中，Ai表示第i個(gè)單元的面積，l為三角形的邊長(zhǎng)，對(duì)應(yīng)三角形有：

(8)

式(7)的求解可分為三個(gè)部分，對(duì)流項(xiàng)的數(shù)值通量求解，紊動(dòng)項(xiàng)的求解和源項(xiàng)中底坡項(xiàng)的處理。對(duì)流項(xiàng)基面數(shù)值通量的求解采用Roe格式的近似Riemann解[16]。紊動(dòng)黏性項(xiàng)采用單元交界面的平均值進(jìn)行估算[17]。模型的詳細(xì)求解過程可見文獻(xiàn)[18]。

3.3 底坡源項(xiàng)處理

底坡源項(xiàng)采用類似對(duì)流項(xiàng)特征向量分解形式離散,克服了簡(jiǎn)單處理底坡源項(xiàng)出現(xiàn)的方程不平衡問題，有效避免了靜水虛假流動(dòng)現(xiàn)象的產(chǎn)生。

為了保證底坡項(xiàng)和壓力項(xiàng)的平衡，將底坡源項(xiàng)改寫成下式(h表示靜水深)：

通過三角單元的每條邊參照對(duì)流項(xiàng)特征矩陣分解格式，令：

4 網(wǎng)格剖分及計(jì)算參數(shù)確定

4.1 網(wǎng)格剖分

數(shù)學(xué)模型計(jì)算域北起虎門口，南至大萬山島以南5 km處，西邊界止于珠海的炮臺(tái)山，東邊界止于香港汲水門。模型的空間尺度為東西向51 km，南北向108 km，控制面積達(dá)3 877 km2。計(jì)算域內(nèi)共包含有10萬余個(gè)三角單元，對(duì)于水面寬闊的外伶仃洋海域，采用大尺度網(wǎng)格剖分，而對(duì)各主要航道途徑水域和人工島附近水域進(jìn)行網(wǎng)格加密。為了提高工程前后水動(dòng)力變化數(shù)值模擬的精確度，對(duì)橋墩和人工島也進(jìn)行了網(wǎng)格剖分，在現(xiàn)狀條件下，所有剖分的網(wǎng)格都參與模型計(jì)算，橋墩和人工島邊壁作為固體邊界進(jìn)行處理(剖分效果見圖2)。

圖2 人工島和橋墩網(wǎng)格剖分示意

4.2 泥沙參數(shù)選取

根據(jù)珠江口伶仃洋海域泥沙運(yùn)動(dòng)的特性以及珠江口多項(xiàng)建港條件水沙環(huán)境論證的工程經(jīng)驗(yàn)，模型相關(guān)的泥沙計(jì)算公式及參數(shù)選取采用如下：

水流挾沙率：

(9)

泥沙沉降速率：

(10)

泥沙起動(dòng)流速：

(11)

泥沙懸浮流速：

uf=0.812d0.4ω0.2H0.2

(12)

泥沙干容重：

(13)

泥沙沉降機(jī)率：

a=0.6

(14)

式中：ω0為泥沙在清水中的沉速；d為泥沙粒徑，d50為中值粒徑；γs和γ分別為海水和清水密度。

5 模型驗(yàn)證

5.1 水動(dòng)力模型驗(yàn)證

模型采用洪季大潮水情(2007年8月13日至14日)和枯季大潮水情(2009年3月27日至28日)兩組水文實(shí)測(cè)資料對(duì)數(shù)值模型進(jìn)行了驗(yàn)證計(jì)算。汛期大潮水情組合有9站潮位和15條垂線的流速資料可供模型驗(yàn)證，枯季大潮水情組合則有11站潮位和11條垂線的流速資料可供模型驗(yàn)證。圖3至圖5分別為部分潮位站和橋區(qū)水域水文垂線的潮位和流速、流向過程驗(yàn)證。從這些圖中點(diǎn)、實(shí)曲線的線型以及峰谷位置比較來看，模型所模擬的各站潮位過程與實(shí)測(cè)過程吻合，模型計(jì)算出的各條垂線的流速大小、流向變化過程形態(tài)與實(shí)測(cè)結(jié)果基本趨于一致。計(jì)算的高、最低潮位和實(shí)測(cè)值偏差在10 cm以內(nèi)，漲潮和落潮時(shí)段的平均流速偏差在±10%以內(nèi)，相位誤差小于0.5 h。驗(yàn)證結(jié)果總體上符合模擬規(guī)程的精度要求。

圖3 潮位過程驗(yàn)證

圖4 流速過程驗(yàn)證

圖5 流向過程驗(yàn)證

5.2 沖淤模型驗(yàn)證

懸沙數(shù)學(xué)模型在對(duì)海床沖淤變化驗(yàn)證時(shí)選取港珠澳大橋的人工試挖槽作為驗(yàn)證對(duì)象。港珠澳大橋工程的試挖槽選擇在西人工島東側(cè)540 m處的水域，海床平均標(biāo)高-9.5 m(理基)，槽底寬21 m，槽底長(zhǎng)100 m，槽型為東西走向，挖深至底標(biāo)高-21 m，相對(duì)開挖深度平均為11.5 m。試挖槽邊坡比分別取1∶5、1∶6(南側(cè))和1∶8、1∶10(北側(cè))，東西兩端縱向邊坡比則均取1∶10，形狀如圖6所示。

圖6 港珠澳大橋試挖槽平面示意

從試挖槽成型后的底槽平均水深變化來看，從2009年2月6日到5月8日三個(gè)月的枯季水情期間，試挖槽底槽平均淤厚為0.40 m；而在5月8日到10月13日五個(gè)月的洪季水情期間，試挖槽底槽平均淤厚為0.83 m，見圖7。

圖7 試挖槽實(shí)測(cè)底槽平均淤積厚度變化

數(shù)學(xué)模型通過枯季和汛期兩種水情模擬試挖槽的泥沙回淤過程。首先采用枯季大潮過程作為計(jì)算水情條件模擬枯季三個(gè)月的泥沙回淤過程，根據(jù)枯季期間的泥沙回淤厚度對(duì)試挖槽水下地形進(jìn)行修正，而后采用洪季大潮過程作為計(jì)算水情條件模擬洪季五個(gè)月的泥沙回淤過程。根據(jù)模型預(yù)測(cè)的試挖槽泥沙回淤結(jié)果(圖8)來看，試挖槽淤積形態(tài)表現(xiàn)為槽底大、邊坡小的分布特點(diǎn)，邊坡淤積又呈現(xiàn)出北邊坡大于南邊破、西邊坡大于東邊坡，這與試挖槽實(shí)際淤積形態(tài)的平面分布比較相似。模型計(jì)算結(jié)果表明，試挖槽在枯季三個(gè)月的平均淤厚為0.38 m，汛期五個(gè)月平均淤厚為0.78 m，與實(shí)際淤厚誤差在10%以內(nèi)，符合泥沙沖淤數(shù)值模擬的精度要求。

圖8 試挖槽淤積分布模擬

6 人工島水域沖淤預(yù)測(cè)結(jié)果分析

根據(jù)圖9顯示的人工島工程實(shí)施一年后的海床沖淤變化數(shù)值模擬結(jié)果，人工島上、下游均有梭狀淤積體形成。對(duì)比人工島上、下游的淤積體分布形態(tài)可以看出，人工島對(duì)伶仃洋水域的水沙環(huán)境影響集中在人工島上下游各5 km水域，呈現(xiàn)出人工島兩端沖刷、上下游形成以島為中心的帶狀淤積體。島南側(cè)的帶狀淤積體較長(zhǎng)、范圍較大，島北側(cè)的淤積體范圍較小，淤積體范圍的大小和人工島南、北兩側(cè)的回流范圍有一定關(guān)系。人工島北側(cè)淤積體的淤積厚度相對(duì)島南側(cè)要大一些，淤積強(qiáng)最大超過2.0 m/a。東、西兩人工島相比較而言，西島南、北兩側(cè)形成的淤積體無論范圍還是強(qiáng)度都要比東島大一些。兩人工島南北側(cè)形成的帶狀淤積體范圍并沒有波及到通航區(qū)的伶仃航道、銅鼓航道和榕樹頭航道。人工島兩側(cè)的挑流作用使得島兩側(cè)均出現(xiàn)不同程度的沖刷，西島的西側(cè)和東島的兩側(cè)形成的沖刷范圍較大，沖刷強(qiáng)度超過0.80 m/a。兩人工島的束水作用使得通航區(qū)的潮流動(dòng)力增強(qiáng)，銅鼓航道(西線)的部分航道出現(xiàn)沖刷，伶仃航道穿過主通航區(qū)一段航道的航槽淤積呈減小趨勢(shì)。

圖9 人工島工程實(shí)施一年后，人工島水域海床沖淤平面分布

模型計(jì)算的設(shè)計(jì)方案考慮了人工島兩端的防撞墩，未考慮人工島周邊的護(hù)坡和島隧結(jié)合部的護(hù)坦。而在工程建設(shè)過程中則取消了人工島兩端的防撞墩，增加了護(hù)坡和護(hù)坦。因此人工島周邊的實(shí)際沖淤情況與模型預(yù)測(cè)會(huì)有一些差異，但人工島對(duì)伶仃洋海床沖淤的總體趨勢(shì)基本一致。

人工島上、下游梭狀淤積體的形成，主要是由于人工島改變了原有的流場(chǎng)，如圖10所示，人工島上下游成為流速減小區(qū)，兩側(cè)成為流速增加區(qū)；人工島背水面的回流區(qū)成為流速減小的主要水域，迎水面一側(cè)也有較小范圍的流速減小區(qū)。人工島的流速減小區(qū)范圍和形狀同人工島的回流范圍和形狀基本相似，西人工島漲急時(shí)刻流速減小區(qū)正北向、落急時(shí)刻南偏東向；東人工島漲急時(shí)刻流速減小區(qū)北偏東，落潮時(shí)刻南偏西。人工島迎水面和背水面形成的弱流區(qū)為泥沙回淤提供了良好的環(huán)境，因此人工島形成的泥沙沖淤部位和強(qiáng)度與人工島周邊的水流條件密切相關(guān)。

圖10 西人工島工程前后附近水域流態(tài)

7 人工島附水域工程后實(shí)測(cè)淤積速率

人工島工程在2009年建成，在2019年對(duì)人工島水域進(jìn)行了水下地形測(cè)量。通過2009年至2019年間的水下地形實(shí)測(cè)年平均沖淤速率分布情況(圖11)可以看出，東、西人工島和隧道區(qū)在大橋沿程產(chǎn)生的水下地形沖淤變化相對(duì)明顯。隧道沿程水深出現(xiàn)了出現(xiàn)了較大幅度的增加，這是隧道沉管完成后未將基槽完全回填而形成的。西人工島形成以島為中心、南北走向的淤積帶，南、北兩側(cè)的淤積帶長(zhǎng)度分別為5.2 km和4.5 km，淤積體東西向最大寬度為1.5 km，最大淤厚區(qū)域的淤積速率在0.3 m/a左右；東人工島同樣形成以島為中心的淤積體，南、北兩側(cè)的淤積帶則呈NNE-SSW向，長(zhǎng)度分別為5.5 km和4.2 km，淤積體的最大寬度略小于西人工島。上述兩個(gè)人工島形成的淤積區(qū)主要是由島對(duì)漲、落潮流的阻流作用產(chǎn)生，淤積帶的方向與島所在水域的潮流運(yùn)動(dòng)方向一致，這與數(shù)學(xué)模型預(yù)測(cè)的人工島淤積分布形態(tài)是一致的。人工島建設(shè)初期，島體對(duì)周邊海床產(chǎn)生的工程效應(yīng)會(huì)相對(duì)較大，但隨著工程周邊水沙動(dòng)力環(huán)境和灘槽格局的逐步適應(yīng)，工程產(chǎn)生的水下地形沖淤速率會(huì)隨時(shí)間逐漸較小[5]。2009—2019年10年間人工島淤積體的實(shí)測(cè)年平均淤積速率約為0.3 m/a，數(shù)值模型預(yù)測(cè)工程建設(shè)后一年的最大淤積速率2.0 m，該淤積速率會(huì)隨著時(shí)間推移幅度逐步減小?？傮w而言，數(shù)學(xué)模型預(yù)測(cè)結(jié)果無論從淤積形態(tài)還是幅度上都與實(shí)際結(jié)果比較相符，這為工程設(shè)計(jì)階段所采用的數(shù)學(xué)模型模擬精度和模型適用性都提供了良好的佐證。

圖11 工程水域2009—2019年海床年平均沖淤速率分布情況

8 結(jié) 語

通過平面二維潮流懸沙數(shù)值模型模擬港珠澳大橋人工島建成后對(duì)其附近水域地形沖淤變化的影響。在工程建設(shè)前，通過水沙資料和試挖槽實(shí)測(cè)地形資料對(duì)模型進(jìn)行了很好的驗(yàn)證，在此基礎(chǔ)上預(yù)測(cè)了人工島工程建設(shè)后的地形沖淤演變，為工程建設(shè)決策提供了相關(guān)依據(jù)。采用工程建設(shè)10年前后的實(shí)測(cè)水下地形沖淤變化情況對(duì)設(shè)計(jì)階段數(shù)學(xué)模型的預(yù)測(cè)結(jié)果進(jìn)行了檢驗(yàn)，人工島建設(shè)引起的海床沖淤變化趨勢(shì)與數(shù)學(xué)模型預(yù)測(cè)結(jié)果基本一致，此為工程設(shè)計(jì)階段所采用的數(shù)學(xué)模型模擬精度和模型適用性提供了良好的佐證，也為該數(shù)學(xué)模型能為類似跨海通道工程在水沙運(yùn)動(dòng)和水下灘槽演變預(yù)測(cè)方面的應(yīng)用提供了有力證明。