李大鳴，陽　婷，李楊楊，歐陽錫鈺

（天津大學　水利工程仿真與安全國家重點實驗室，天津　300072）

?

潮流泥沙數學模型在青島港擋沙堤工程的應用

李大鳴，陽婷，李楊楊，歐陽錫鈺

（天津大學水利工程仿真與安全國家重點實驗室，天津300072）

摘要：考慮波浪輻射應力對潮流場和泥沙運動的影響，建立了青島港前灣三期碼頭前沿擋沙堤工程二維潮流泥沙數學模型。在采用實測資料驗證的基礎上，運用模型對本海域在無擋沙堤及不同擋沙堤長的各種方案的流場變化和泥沙回淤情況進行計算研究。結果表明，無擋沙堤時，由于三期工程的建設縮窄了河口至海區(qū)的斷面面積，斷面西側的淺水區(qū)水流速度增大，容易掀起泥沙輸移至斷面東側開挖后的深水區(qū)，使泥沙在碼頭前港池中落淤，碼頭前沿最大淤積強度約為0.818 m/day；而建設擋沙堤后將顯著減小碼頭前沿的泥沙淤積。經過比較，從擋沙堤附近流場與港池航道回淤情況的角度考慮，認為方案二對碼頭前沿攔沙的效果較好。

關鍵詞：波浪輻射應力；潮流；泥沙；數學模型；擋沙堤；青島港

近岸水動力環(huán)境對海岸帶及港工建筑物的影響十分重要，而潮流、波浪和泥沙等相互作用又使得海岸帶附近的水動力環(huán)境非常復雜。在對其的諸多研究方法中，以數學模型最為方便經濟，并成為研究大范圍水動力條件變化的主要手段之一（韓亮，2010）。近岸潮流、波浪和泥沙運動規(guī)律的基本方程大多都是由偏微分方程來描述的。通過數值計算方法求解這些方程，首先要將這些方程進行離散形成代數方程組，然后在計算機上求解（嚴愷，2002）。常使用的方法有有限差分法、有限單元法、有限體積法、有限分析法和邊界元法等；按差分網格形狀的不同可以分為三角形、矩形、四邊形、多邊形、曲線坐標網格以及不同網格的組合等（張滌明，1991）。

在對近岸水動力環(huán)境的大量研究中，河口海岸地區(qū)波浪與潮流共同作用下的泥沙運動及其模擬屬于前沿課題。早期研究中，Swart（1976）通過將波浪、潮流共同作用下的床面剪切應力替代河流泥沙輸運公式中的剪切應力來研究河流輸沙問題；David等（1991）采用波流相互作用模型來計算底部切應力，認為在研究泥沙再懸浮及淺水河口床面泥沙運動時考慮波流相互作用是十分重要的。隨著相關研究的發(fā)展，不少學者采用流場模型，在懸沙擴散方程中采用波流共同作用下的挾沙能力來研究泥沙擴散問題（竇希萍，1992）；竇國仁等（1995）導出了波浪和潮流共同作用下的懸沙輸移方程式和挾沙能力公式，建立了河口海岸平面二維泥沙數學模型。另外，也有許多研究者在潮流方程中引入輻射應力（嚴以新，1988），結合懸沙擴散方程討論波流共同作用下的泥沙擴散（曹祖德等，1993）；辛文杰（1997）將“波浪輻射應力”、“波流摩阻力”以及“波流挾沙能力”三個要素歸納到水流運動方程和懸沙輸送方程中去，建立潮流、波浪綜合作用下的河口二維懸沙數學模型。從其他角度出發(fā)，白玉川等（2000）根據N-S方程和質量傳輸方程，采用Reynolds分解法建立了模擬波浪潮流共同輸沙及海岸沖淤演變模式；吳永勝等（2002）通過引入紊動封閉方程，耦合波浪和水流之間的相互作用因素，建立了統(tǒng)一描述波浪、水流流場在波浪邊界層內、外的變化規(guī)律的數學模型。李大鳴等（2014）考慮了波浪場產生的“波浪輻射應力”和“波流挾沙力”動力要素，采用有限元加權質量法建立了波浪作用的潮流泥沙數學模型，對本文有著很大的理論指導意義。

青島港前灣港區(qū)地處山東半島膠州灣黃島東南岸線的前灣內，其一﹑二期工程已投入運營。為滿足港口貨運增長的要求，計劃在一﹑二期工程西南面沿岸線，進行三期工程的建設（圖1）。隨著三期工程向西發(fā)展，工程占用了前灣灣頂辛安河口的一部分海區(qū)面積，導致河道縮窄；與此同時，三期工程最西側港池、航道的挖深使得青島前灣納潮量增加，港池邊坡塌落，淺水區(qū)的掀沙隨潮流向深水區(qū)落淤的泥沙將威脅西側港池航道的正常水深。為順利實施青島港三期工程和減少港池回淤，擬在三期工程最西端建一座擋沙堤，故本文建立波浪潮流泥沙數學模型，對該區(qū)的流場動力條件和泥沙回淤情況進行分析研究，為青島港前灣三期工程的實施及今后港口的建設、維護提供科學依據。

圖1　青島港前灣港區(qū)地理位置和模型范圍

1數學模型理論

1.1波浪場數學模型

1.1.1基本方程

采用三角網格剖分的有限海域折繞射聯(lián)合計算，并依據波能衰減的緩坡方程（潘軍寧等，2001）來確定波浪場：式中：φ表示波勢；h=（?/?x，?/?y）；C表示波速；Cg表示群波速；ω為角頻率；k表示波數；F為波能衰減因子；

1.1.2邊界條件

入射邊界條件：

透、反射邊界條件：

式中：γ為反射參數，由下式計算：

式（2）至（6）中，n為入射波向；表示入射波勢；φr表示反射波勢；R為反射系數；ε為相位差；R =1時為全反射邊界，R =0時為全透射邊界。

1.2考慮波浪作用的潮流場數學模型

1.2.1基本方程

采用三角網格剖分的有限元加權集中質量法計算潮流場（吳克田，1990）。依據波浪輻射應力作用的二維淺水環(huán)流方程：

其中：u、v為流速在x、y方向的分量；z為水位；g為重力加速度；c為謝才阻力系數；h為水深；f為柯氏力系數；S為輻射應力；AH為水平渦粘系數；H為波高；k為波數；ρ為水體密度；α為波向角。

1.2.2邊界條件

式中：Z為（x0，y0）點上的潮位過程；Γ1為模型水域水線邊界；θ為模型陸邊界法向角；Γ2為模型水域陸線邊界；Ω為平面計算域。

1.3泥沙場數學模型

依據以下方程采用三角網格剖分的三角差分方法計算懸沙含沙量、底沙量和海床變化。

1.3.1考慮紊動水動力影響的泥沙擴散方程（王尚毅，1990）

其中：Cs為懸沙平均含沙量；Cb為床面層內部底沙含沙量；εx、εy為x、y方向的擴散系數；ωs為懸沙平均沉速；d50s為懸沙中值粒徑；As為與懸沙平均沉速和中值粒徑相關的修正系數。

1.3.2考慮兩相流相互作用的底沙輸移的連續(xù)和運動方程

其中：ωb為底沙平均沉速；d50b為底沙中值粒徑；Ab為與底沙平均沉速和中值粒徑相關的修正系數；v為水流運動粘性系數；ks為廣義滲透率。

1.3.3床面變形方程

其中：γc為底沙干容重；γs為懸沙容重；CA為全沙含沙量。（注：以上平均含沙量均采用含沙相對體積比）。

1.3.4邊界條件

模型水域水線邊界Γ1處網格節(jié)點的含沙量以平衡挾沙理論（曹志先，1997）為基礎，根據該邊界處的水動力條件及附近海域的泥沙中值粒徑計算確定；而模型水域陸線邊界Γ2處網格節(jié)點的含沙量始終為0，即：

海床底部變化邊界條件由（16）式確定。

2模型布置

青島港前灣港區(qū)潮流泥沙數學模型的模型范圍如圖1所示，東起團島鼻與腳子石嘴一線海域，西至辛安河口西岸線（120°17＇24＂E-120°10＇24＂E），東西最大寬度約10.5 km；南起海西灣灣頂，北到團島嘴與黃島東岸一線海域（35°58＇N-36°3＇30＂N），模型長約10.5 km。模型覆蓋海域面積達60 km2。模型中部采用規(guī)則的有限元三角形單元進行剖分，邊界則采用不規(guī)則三角形網格使工程和岸線與單元網格完全擬合。驗證模型以三期工程建設前的地形數據為基礎，空間步長約為160m，平面剖分網格見圖2左部分。為較準確的描述前灣港區(qū)三期工程西端擋沙堤對淺灘泥沙的攔擋作用，計算模型采用如圖2右部分所示的平面剖分網格，在三期工程附近的淺灘和港區(qū)周邊取1.6 km×8 km的矩形加密區(qū)，空間最小步長約為25m。對不同長度的擋沙堤，模型平面剖分網格略有差別，剖分網格圖略。

圖2　驗證模型剖分網格（左）及擋沙堤建設前模型加密網格（右）圖

3模型驗證

3.1觀測站潮流驗證

模型驗證計算采用2001年6月5日8時-6日9時大潮過程、6月10日8時-11日9時中潮過程和6月17日9時-18日10時小潮過程對模型中四個驗潮站進行了驗證，驗潮站分別為L1、L2、L3 和L4，位置如圖3所示。通過對四點全潮逐時的潮位、流速、流向驗證資料進行調試驗證，可以看到各點潮位吻合良好，流速、流向趨勢基本一致。限于篇幅，這里只給出各測點大潮潮位、流速、流向逐時過程驗證結果，見圖4。

圖3　驗潮站位置示意圖

3.2潮流場驗證

模型還分別驗證了整個模型區(qū)域在大、中、小潮過程的漲落潮流場。這里只給出大潮過程的漲落潮流場圖，見圖5。從圖中可以看出，漲潮主流由外海經灣口流入，一股強流向北流入膠州灣，一股較弱的分流流向西南的前灣，向南流入海西灣的水流流量最?。宦涑狈謩e由膠州灣、前灣和海西灣匯流后流向外海，流態(tài)與實際觀測基本一致。

圖4　大潮潮位、流速、水位過程驗證

圖5　驗證模型大潮落潮（左）、漲潮（右）流場圖

3.3觀測站含沙量驗證

含沙量的驗證計算采用中、大、小潮同期資料。其中各個測點大潮驗證結果見圖6，與實測過程基本一致；整個模型區(qū)域大潮落、漲急情況的含沙量分布見圖7，可見含沙量分布基本合理，模型基本能反映原型的物理特征。

圖6　大潮含沙量過程驗證

圖7　模型驗證大潮含沙量分布圖

4擋沙堤方案計算

4.1擋沙堤布置方案

青島港前灣港區(qū)三期工程擬建的擋沙堤由碼頭前沿向東南延伸。為比較不同長度擋沙堤的攔沙效果，分別設計長度為180m、200m和220m，將這三種情況分別稱為第一、第二和第三方案，未建擋沙堤前的情況稱為工程前方案。為比較各方案擋沙堤附近流場與港池航道回淤情況，在擋沙堤擬建位置附近沿三期工程泊位岸線，選擇5個特征點，用于定量比較工程前后主要物理量的變化，特征點距擋沙堤距離分別為0.0m、100m、200m、300m 和400m。擋沙堤及各個特征點位置如圖8所示。

圖8　擋沙堤附近特征點位置示意圖

4.2各方案計算

4.2.1潮流場計算

考慮上述四種方案下對模型區(qū)域在大、中、小潮的漲落潮進行潮流場計算。其中，在青島港前灣港區(qū)三期工程建設完成、港池航道開挖后無擋沙堤情況下的工程前方案中，由于工程占用海區(qū)，三期工程西端形成河口到海區(qū)縮窄斷面，斷面西側為淺水區(qū)，水流速度較大，斷面東側為開挖后的深水區(qū)，水流速度較小，形成了一個大小流速交換區(qū)，如圖9左邊部分所示的大潮落潮過程。而擋沙堤工程建設后，縮窄斷面處繼續(xù)縮窄，使淺水區(qū)與深水區(qū)縮窄斷面離開碼頭前沿，大流速與小流速交換區(qū)離開碼頭前沿。工程第二方案，即擋沙堤為200m時，三期工程附近大潮期落潮流場如圖9右部分所示；第一和第三方案潮流場與第二方案相似，不再以圖形方式給出。表1給出各方案特征點上大潮落潮的極值情況。

圖9　擋沙堤工程前（左）及第二方案（右）大潮落潮流場

表1　各方案特征點大潮落潮流速極值

4.2.2日淤積強度計算

從上述潮流場的計算可以看出，三期工程西端形成的河口到海區(qū)縮窄斷面西側為水流速度較大的淺水區(qū)，容易帶動淺水區(qū)掀起的泥沙隨水流運動；斷面東側為水流速度較小的開挖后深水區(qū)，泥沙不容易起動，水流挾沙能力較弱；在落潮流作用下，淺水區(qū)泥沙隨水流經縮窄斷面進入開挖后的深水區(qū)，水流挾沙能力迅速減弱，大部分泥沙快速沉積，少部分泥沙向前灣口方向擴散，但在漲潮流作用下，少部分泥沙的一部分會隨水流作用回到縮窄斷面或擋沙堤東側落淤，增加了斷面處深水側的淤積厚度。根據潮流場計算結果得到開挖完成后各方案大、中、小潮作用24 h后的含沙量分布與港池航道淤積圖，以第二方案大潮過程淤積情況及含沙量分布為例，如圖10所示；各方案特征點上大潮24 h累積淤積量以表2給出。

表2　大潮不同方案擋沙堤東側不同部位日淤積強度

從以上圖表中可以看出，不加擋沙堤的工程前方案，碼頭前沿在大潮情況下淤積強度最大，最大淤積強度發(fā)生在淺水區(qū)與深水開挖區(qū)的交界斷面處，最大數值為0.818m/day；建設擋沙堤后碼頭前沿淤積厚度減小，最大淤積強度仍發(fā)生在大潮情況下，淤積部位為擋沙堤東側堤腳，第一方案擋沙堤長度為180m時，大潮最大淤積強度為0.222m/day；第二方案擋沙堤長度為200m時，大潮最大淤積強度為0.090 m/day；第三方案擋沙堤長度為220 m時，大潮最大淤積強度為0.172m/day，第二方案擋沙堤攔沙作用最明顯，可將碼頭前沿的最大淤積強度減少約89％。

圖10　擋沙堤工程第二方案大潮淤積分布（左）及含沙量分布（右）

4.2.3泥沙年淤積強度計算

青島港前灣港區(qū)三期工程建設完成、港池航道開挖后，初期會由于邊坡塌落和淺灘來沙出現較高強度的回淤，但隨著回淤厚度的增加回淤強度會逐漸減弱?？紤]大、中、小潮作用，模擬港池航道開挖后一年回淤情況，如表3所示，比較工程前方案與擋沙堤第二方案特征點上的年淤積強度，工程前方案與擋沙堤第二方案年淤積強度分布見圖11。

表3　第二方案與工程前方案年淤積強度比較

圖11　擋沙堤工程前（左）及第二方案（右）年淤積強度分布

從圖表中可以看出，建設擋沙堤主要對三期工程西側，有明顯的減淤效果。在表3中，特征點上工程前年最大淤積強度為0.851m/a；建設200m擋沙堤后，碼頭前沿的年最大淤積強度為0.212m/a，最大淤積強度位于第5特征點處，與工程前方案相比約減少75％。通過對比表2和表3可以發(fā)現，工程第二方案最大日淤積強度與最大年淤積強度位置不同，前者在第2點位處最大而第5點位處最小，后者則相反；因為前者計算的是港池開挖后某一特征潮作用下24 h內的淤積高度，而后者計算的是全年的淤積高度。港池開挖后，在一次大潮、中潮或小潮獨立作用的時間內，泥沙會在開挖后港池邊緣擋沙堤堤頭處的地方先落淤，受此時淤積帶的影響，擋沙堤東側第2點位處淤積強度相對較大，而此時第5點位處還未受到淤積帶影響；而在全年潮流的作用下，港池內泥沙淤積帶向第五點位方向發(fā)展，最大淤積強度點也向港池東側推移。從圖11中還可以看出擋沙堤能夠引起三期工程港池航道內淤積分布的改變，但港池航道整體的淤積總量變化不大，泥沙淤積空間分布也逐漸趨于均勻。

5結論

本文建立了考慮波浪作用的二維潮流數學模擬計算模型，并使用實測海流資料驗證所建立的潮流數學模型的合理性。通過對擋沙堤工程各方案進行典型潮的模擬潮流計算，分析工程和流態(tài)的相互影響，并從擋沙堤工程布置方案上對流場變化、泥沙淤積和運移角度加以分析，得到以下結論：（1）青島港前灣港區(qū)三期工程建設完成、港池航道開挖后，在無擋沙堤的工程前方案中，碼頭前沿在大潮情況下淤積強度最大，最大數值為0.818 m/day。（2）青島港前灣港區(qū)三期工程建設200m擋沙堤時，對碼頭前沿攔沙效果最好，大潮時碼頭前沿最大淤積強度減少為0.090m/day，比未建設擋沙堤的情況，碼頭前沿的最大淤積強度減少約89 ％，擋沙堤對減小三期工程碼頭前沿淤積的作用明顯。（3）建設200m擋沙堤后，碼頭前沿的年最大淤積強度為0.212m/a，約減少75％，擋沙堤對減小三期工程碼頭前沿年淤積強度的作用也較明顯。

參考文獻

Schoellhamer David H,Levesque Victor A,1991.W ind generated wave resuspension of sediment in Old Tampa Bay.Proceedings-National Conferenceon Hydraulic Engineering,(29):85-90.

Swart D H,1976.Coastal sediment transport.Computation of Longshore Transport,ReportR968.

白玉川，顧元棪，蔣昌波，2000.潮流波浪聯(lián)合輸沙及海床沖淤演變的理論體系與其數學模擬.海洋與湖沼，31（2）：186-196.

曹志先，1997.泥沙數學模型近底邊界條件I：平衡輸沙.水利學報，（1）：11-19.

曹祖德，王桂芬，1993.波浪掀沙、潮流輸沙的數值模擬.海洋學報，15（1）：107-118.

竇國仁，董鳳舞，竇希萍，李禔來，1995.河口海岸泥沙數學模型研究.中國科學，25（9）：995-1001.

竇希萍，羅肇森，1992.波浪潮流共同作用下的二維泥沙數學模型.水利水運科學研究，（4）：331-338.

韓亮，2010.二維潮流波浪數學模型及工程應用.天津大學.

李大鳴，歐陽錫鈺，潘番，等，2014.考慮波浪輻射應力的潮流泥沙數學模型.海洋通報，33（6）：703-711.

潘軍寧，左其華，洪廣文，2001.一種推廣的緩坡方程.海洋工程，19（1）：24-31.

王尚毅，顧元蒳，郭傳鎮(zhèn)，1990.河口工程泥沙數學模型.海洋出版社，120-192.

吳克田，1990.二維淺水波問題Galerkin解法的質量集中法.計算物理，7（1）：1-6.

吳永勝，練繼建，王兆印，等，2002.波浪-水流相互作用模型.水利學報，（4）：13-17.

辛文杰，1997.潮流、波浪綜合作用下河口二維懸沙數學模型.海洋工程，15（1）：30-47.

嚴愷，2002.海岸工程.海洋出版社，622-625.

嚴以新，1988.潮汐河口地區(qū)波流相互作用的數學模型.海洋工程，7（3）：47-56.

張滌明，1991.計算流體力學.中山大學出版社，181-187.

（本文編輯：岳心陽）

App lication of thenum ericalm odelof tidal currentsand sedim ents in the dike engineering of Qingdao Port

LIDa-ming,YANG Ting,LIYang-yang,OUYANG Xi-yu
（State Key LaboratoryofHydraulic EngineeringSimulation and Safety,Tianjin University,Tianjin 300072,China）

Abstract：In this paper,a 2-D tidal currentand sedimentnumericalmodel for the dike engineering at the former bay of Qingdao Port Phase 3 is established in view of＂wave radiation stress＂generated by the wave field on the flow field and sedimentmovement.Based on the observed data,themodel isapplied to study the change of currentsand back silting in the portareawithout the dike and with the dike of different length.The results show thatwhen there is no dike here,the cross section area between the estuary and sea narrows due to the construction of the Phase 3,and itcauses the increasing ofwater flow rate at the shallow water on the west side of cross section,which makes sediments transport from the west to the east deepwater side easilywith the sedimentation in theharborbasin.In this case,themaximum silting intensity ofwharfapron is about0.818m/day.W hereas the dikebuildingwillsignificantly reduce the sedimentation of thewharfapron.According to the comparison,case2 isbetter in view of the sedimentdeduction.

Keywords：wave radiation stress;tidalcurrents;sediments;numericalmodeling;dike;Qingdao Port

通訊作者：陽婷，碩士研究。電子郵箱：544052703@qq.com。

作者簡介：李大鳴（1957-），男，博士，教授，主要從事水力學與河流動力學研究工作。電子郵箱：lidaming@tju.edu.cn。

收稿日期：2014-12-29；

修訂日期：2015-04-30

Doi：10.11840/j.issn.1001-6392.2016.01.014

中圖分類號：TV148

文獻標識碼：A

文章編號：1001-6932（2016）01-0103-09