波流聯合作用下廣東陽西電廠海域二維泥沙數值模擬研究

譚 超，黃本勝，邱 靜，黃廣靈，劉 達

(1.廣東省水利水電科學研究院，廣州 510630；2.河口水利技術國家地方聯合工程實驗室，廣州 510630；3. 廣東省水動力學應用研究重點實驗室，廣州 510630)

0 引 言

在河口海岸地區, 引起泥沙運動的主要動力因素是潮流和波浪, “波浪掀沙、潮流輸沙”被公認為是泥沙運動機制[1,2]。近岸水動力環境對海岸帶及港工建筑物的影響十分重要，而潮流、波浪和泥沙等相互作用又使得海岸帶附近的水動力環境非常復雜。因此，在復雜動力條件下，海岸工程建設對沿岸泥沙輸運的影響過程研究備受國內外學者關注，也取得了系列研究成果[3-9]。在對其諸多研究方法中，以數學模型最為方便經濟，并成為研究大范圍水動力條件變化的主要手段之一[10]。 電廠的建設改變了岸線形態，導致水動力條件和泥沙運動的變化。同時，華南海岸經常遭受不同程度臺風的侵襲，臺風期間產生的大浪和風暴潮極易掀起海床泥沙，使含沙量劇增，泥沙被水流搬運至近岸處發生驟淤，對近岸港口、航道及取、排水口的安全等造成較大影響[11,12]。

廣東陽西火力發電廠位于廣東省陽江市海陵灣的外灣-青海灣的西端海岸地段。青海灣呈半圓形，為岬角控制的弧形海灣，其東側為海陵島閘坡岬角，西端為白虎嶺岬角。2個岬角在地貌上均屬于侵蝕剝蝕丘陵，由花崗巖基岸及風化殼組成。岬角附近海域分布有礁石和小島嶼，較大的有雙山島。灣口大體向南，東西寬約14 km，灣頂為沙質海灘，其東北端為與海陵灣內灣相連潮汐水道。清海灣灣口水深約10 m，海灣東部大部水深為2～5 m，西部則一般為5～10 m。工程附近海域雙山島至馬村-白虎嶺沿岸水深多變，一般為2～5 m，最深處為白虎嶺南側岬角附近，其10 m等深線距岸僅400 m左右， 研究區域位置見圖1。

圖1 研究區域位置

本文擬建立波流聯合作用下的二維潮流泥沙數學模型，對陽西電廠附近海域的潮流場和泥沙場進行模擬研究，計算取水口及防波堤區域沖淤量，分析大風浪作用下泥沙運動規律，其成果不僅可加深對海岸泥沙輸運物質輸移的關鍵過程的認識，同時可為陽西電廠及類似工程設計及運行維護提供技術支撐。

1 數學模型與工程方案

1.1 基本方程

1.1.1 潮流數學模型

連續方程:

(1)

動量方程:

(2)

(3)

(4)

應力項Ax、Ay為包括水平黏滯應力、表面風應力、底部切應力和波浪輻射應力。模型的空間離散是使用單元中心有限體積法，時間差分格式采用顯式迎風格式。

1.1.2 波浪數學模型

波浪數學模型模擬的物理過程包括：地形和海流空間變化導致的波浪折射作用；地形和海流空間變化導致的淺水變形作用；逆向流造成的障礙和反射作用；障礙物的阻擋或部分傳播作用。其基本方程為：

(5)

N=N(σ,θ,x,y,t)=E(σ,θ,x,y,t)/σ

式中：N為波作用密度譜；σ為相對波頻；θ為波向；Cx和Cy為波浪在x和y方向的傳播速度分量；Cσ和Cθ為σ、θ空間的波浪傳播速度；S是以譜密度表示的源項，包括風能輸入、波與波之間的非線性相互作用和能量耗散。

該模型采用了全隱式有限差分格式對波作用動量守恒能量平衡方程進行離散，然后在4個象限中用迭代的方法進行求解，其計算是無條件穩定的，因而允許較大的時間步長。

1.1.3 泥沙數學模型

(1)基本方程。根據海區泥沙場分析結果可知，本海區泥沙粒徑較小，表層沉積物以粉沙為主，泥沙輸運以懸移為主，用二維懸沙輸運模型來計算泥沙輸運及泥沙驟淤。

黏性泥沙輸運模型涉及泥沙在水體中的運動以及泥沙與底床的相互作用。懸移泥沙的輸運一般建立在水動力模型中的對流項中，可用以下方程來描述：

(6)

(2)波流共同作用下的床面剪切應力。海岸地區的泥沙輸移大多是在波浪和潮流共同作用下完成的，波浪主要起掀沙作用，但本身搬運泥沙能力較弱，而泥沙一旦處于懸浮狀態，一般潮流就能搬運，也即前面提到的“波浪掀沙、潮流輸沙”。所以計算時必須考慮波流共同作用下的泥沙運動。在模型中，床面剪切應力可采用以下公式進行計算。

波浪作用下的底床平均剪切力計算公式為：

(7)

式中：fw為波浪摩擦系數，0.047 (Swart 1974)；Ub為波浪水質點底部軌跡運動速度：

(8)

式中：Hs為有效波高；T為周期。

純流作用下的底床平均剪切應力計算公式為：

(9)

(10)

式中：V為平均流速；fc為流阻系數：h為水深；k為底床糙度。

波流相互作用的底床剪切力計算公式(soulsby eta al.1993)為：

(11)

(12)

式中：τc為純流作用底床剪切力；τw為純波作用底床剪切力；τmax為波流相互作用最大底床剪切力；τmean為波流相互作用平均底床剪切力；a、b、m、n、p、q為常數,由波流理論來確定。

1.2 數學模型構建

1.2.1 計算范圍

為了更好的提供水動力模型所需要的邊界條件以保證局部流場計算符合潮流場的整體物理特征，潮流模擬采用嵌套計算。大模型范圍包括整個海陵灣，小模型范圍和網格均與波浪模型的范圍和網格一致，2層模型的計算區域見圖2。

圖2 模型網格布置

根據實測海流的流向及大小，大模型東邊界定在海陵島下塘見村附近，在電廠以東約34 km，西邊界定在茂名市電白縣博賀鎮北山村附近，在電廠以西約38 km，南邊界定在廠址以南約47 km之處，模型包括整個海陵灣水域。模擬原體水域面積約2 560 km2。模型采用能夠很好地貼合岸線和方便對工程區域加密的非結構三角形網格單元，模型區域共剖分了44 191 個網格，共有22 709 個節點，其中最小的網格單元邊長約為50 m，廠址附近海域使用最小網格，可以較好地模擬岸線變化及擬建電廠平面布置。

1.2.2 邊界選取

(1)波浪邊界條件。本次計算選擇本項目工在工程海域測波站2004 年1月-2005年1月的實測波浪統計資料，測波站的水深為18 m，與波浪模型的南邊界水深基本一致，因此本次波浪計算可以直接采用測波站的波浪統計資料作為計算邊界。根據波浪統計結果，本海區的常浪向為ESE、SE、SSE和S向，其年內出現的頻率分別為25.3%、40.73%、13.27%和13.48%，因此進行波流耦合計算時，需分4組工況進行。各組波浪邊界情況見表1。

表1 波浪模型計算邊界

(2)潮流邊界條件。模型計算區域為近岸水域，因此邊界由潮位來給。鑒于本海域的潮汐比較復雜，本次大模型試驗的邊界采用調和常數換算成水位來試算，然后通過與實測潮位的不斷驗證、調整和耦合，最終得出一個能反映實測情況的潮位邊界，最終應用此邊界作為大模型的計算邊界。大模型計算結果為小模型提供水位或者流量邊界。

(3)泥沙邊界條件。本次建立泥沙數學模型的目的是計算工程后該海域的年內沖淤演變情況，因此，計算水動力的選擇必須按年來考慮。在泥沙模型計算過程中，潮流邊界由調和常數計算出的水位提供，計算邊界包含整個全潮過程。波浪計算邊界則以實測統計資料為基礎，以各個主浪向的波要素為邊界條件，分別計算出各個主浪向波浪作用1 a的沖淤量，然后根據它們各自的年出現頻率進行加權統計，最后得到年內的沖淤變化量。

1.2.3 模型參數取值

(1)水動力計算參數。模型海域的初始潮位取各邊界潮位的平均值，初始流速取0。經反復調試，大、小模型中的渦黏系數取0.28，曼寧糙率取0.025，模型計算時間步長60 s，最小時間步長取0.01 s，CFL數取0.8。

(2)泥沙計算參數。

①海床床面泥沙級配。通過對沉積物樣品的粒度分析表明，工程海域沉積物中值粒徑分布總體上表現為東細西粗、南細北粗的特征，近岸及雙山島附近沉積物中值粒徑大部分小于2 Ф(即d50>0.25 mm)，泥沙顆粒較粗，在5 m等深線以外海域中值粒徑一般大于4 Ф(即d50<0.063 mm)。 沉積物類型以沙、黏土質粉沙為主，其他類型為礫石質沙、粉沙質沙、粉沙質黏土。各類型沉積物空間分布見圖3。

圖3 沉積物類型空間分布圖

②泥沙參數。懸沙中值粒徑d50=0.004 mm，由于泥沙顆粒較小，需考慮泥沙的絮凝作用，在考慮泥沙絮凝作用下的懸沙中值粒徑d50=0.024 mm；絮團泥沙沉速w=0.000 45 m/s；泥沙沉降幾率根據經驗取值a=0.45；懸沙干容重取值r0=720 kg/m3。

1.3 數學模型驗證

潮流數學模型采用2013年夏季的水文資料進行驗證，2013年夏季小潮期實測時間為7月16日12∶00至7月17日13∶00；中潮期實測時間為7月19日12∶00至7月20日13∶00。大潮期實測時間為7月22日14∶00至7月23日15∶00。2013年夏夏季小、中、大潮潮位驗證相對誤差分別為0.044、0.053、0.058 m(見表2)；流速驗證平均絕對誤差控制在0.12 m/s以內(見表3)；從含沙量驗證情況看(見圖4)，計算含沙量與實測含沙量在總體趨勢和峰值的波動變化與實測含沙量過程較為同步且吻合較好。驗證結果表明該模型可用于本工程海域水動力與泥沙模擬。

表2 潮位誤差統計 m

表3 流速誤差統計 m/s

圖4 夏季中潮懸沙驗證曲線

1.4 工程方案

本次擬建工程包括電廠二期取排水工程、碼頭擴建工程。二期工程新建5～8號機組取、排水口，每臺機組取、排水量均為78 m3/s。二期工程同時對碼頭進行擴建，碼頭的擴建方案為在一期卸煤碼頭的基礎上向外延伸310 m，增加一個7 萬t的泊位。

2 計算結果分析

2.1 工程擴建前、后流場變化

2.1.1 流態變化

工程附近海區的影響主要體現在3個方面：防波堤延長改變了原有的海岸線，使流場發生變化；防波堤西南側新加的7、8號取水口改變附近流向，增大流速；電廠北側排水口亦改變附近流向、增大流速，見圖5。

圖5 工程擴建前后流場對比

延長防波堤對海岸線改變比較大，因而對流態的影響也最為顯著，防波堤兩側流速減小，堤頭外由于挑流作用流速增大。

2.1.2 流速變化

圖6為夏季大潮漲急、落急時刻的工程擴建前后流速變化等值線。由圖6可見，工程建設以后，防波堤兩側流速明顯減小，而防波堤端口東側由于挑流作用，流速增大；取排水口出流速明顯增大。工程實施后，港池內平均流速有小幅度降低，范圍為3.1%～9.9%；受到防波堤對漲落潮流的遮擋作用，防波堤北面與南面流速大幅降低，范圍為25.8%～60%；排水口處，受電廠排水影響，流速增大可達167%；排水口附近，但離排水口有一定距離，流速變化不大，為0.8%～10.9%，可見電廠排水對海域流速大小的影響范圍較?。稽c9位于延長防波堤東側，由于挑流作用，流速增大幅度可達28.7%；取水口處，受電廠取水影響，流速增大，增幅最大可達58%。

圖6 工程擴建前后流速變化等值線

2.2 年回淤分析

年回淤量計算采用的波浪邊界為實測統計值(見表1)，以各個主浪向的波要素為邊界條件，分別計算出各個主浪向波浪作用1 a的沖淤量，然后根據它們各自的年出現頻率進行加權統計，最后得到年內的沖淤變化量。計算中采用的潮流動力為典型的全潮過程。

泥沙數學模型計算的波浪、潮流共同作用下工程附近海域地形沖淤變化見圖7[其中圖7(a)為工程擴建前的沖淤結果，圖7(b)為工程擴建后的沖淤結果]，圖7中正值表示淤積，負值表示沖刷。

圖7 波流共同作用下地形沖淤結果

由圖7及表4可見，工程擴建前后工程海域的總體沖淤趨勢基本一致，防波堤南、北兩側均表現為淤積，防波堤東側為沖刷。工程后由于防波堤加長，東側沖刷區相應的東移，沖刷的最大深度有所增加(增加約0.02 m/a)，而防波堤南、北側的淤積區由于防波堤的加長而淤積強度稍有增強(約0.01 m/a)。

表4 工程特征部位的沖淤強度

工程擴建后，取水口東側、南側均為泥沙淤積區，取水口附近年沖淤厚度范圍為0.1～0.2 m/a；加長碼頭防波堤后，堤頭東側水域則會有所沖刷，最大沖刷為0.2 m/a左右；而碼頭港池、泊位、回旋水域則稍有淤積，年淤積厚度為0.04～0.10 m/a。

2.3 驟淤計算

驟淤是惡劣天氣條件(臺風、風暴潮、熱帶氣旋)下短時間內產生的泥沙淤積。根據工程海區的地質資料，廠址附近淺海區泥沙較粗，為沙粒，水深在10 m以上的外海區泥沙粒徑較小為粉沙及黏土等，這種泥沙在風暴潮時很容易起動，水體含沙量將顯著增大，有可能造成取水口及港池的大量淤積，因此必要進行驟淤分析。

本次驟淤計算本試驗選取200 a一遇高潮位4.16 m、遭遇9713號臺風(特大臺風，S向，H4%=6.6 m(H有效=5.5 m)，且連續作用2 d作為臺風驟淤的試驗條件。

從圖8及表5可見，工程擴建后，取水口東南側均為泥沙淤積區，最大淤積厚度可達0.5 m，距取水口位置約1.6 km，取水口附近淤積厚度為0.20～0.25 m；加長碼頭防波堤后，堤頭東側會有所沖刷，最大沖刷可達0.5 m左右；而碼頭港池、泊位、回旋水域則稍有淤積，淤積厚度為0.2～0.5 m。

圖8 9713號臺風過后地形沖淤狀況

3 結 論

本文建立工程海區波流共同作用下的平面二維泥沙數學模型，在利用工程海區實測水沙資料對模型進行驗證的基礎上，計算分析了波流聯合作用下廣東陽西電廠海域潮流泥沙過程及暴風驟淤特征，主要得出以下結論。

(1)陽西電廠二期工程方案實施后，流態變化比較大的區域主要是取排水口附近、防波堤兩側以及防波堤東側水域。排水口處流速明顯增大，增幅可達167%， 取水口處與防波堤東側水域流速也增大，增幅最大達到58%，防波堤兩側流速則明顯減小，降幅為25.8%～60.0%。

(2)陽西電廠二期工程方案實施后，取水口東側、南側均為泥沙淤積區，取水口附近年沖淤厚度范圍為0.1～0.2 m/a；加長碼頭防波堤后，堤頭東側水域則會有所沖刷，最大沖刷為0.2 m/a左右；而碼頭港池、泊位、回旋水域則稍有淤積，年淤積厚度為0.04～0.10 m/a。

(3)臺風作用下，海域含沙量明顯增大，進入港池和取水口的含沙量也隨之增大，因此可造成港池和取水口的較強淤積。驟淤條件下電廠取水口含沙量為0.08～0.10 kg/m3。取水口東南側均為泥沙淤積區，最大淤積厚度可達0.5 m距取水口位置約1.6 km，取水口附近淤積厚度為0.20～0.25 m；加長碼頭防波堤后，堤頭東側會有所沖刷，最大沖刷可達0.5 m左右；而碼頭港池、泊位、回旋水域則稍有淤積，淤積厚度為0.2～0.5 m。