波流聯(lián)合作用下廣東陽西電廠海域二維泥沙數(shù)值模擬研究

譚 超，黃本勝，邱 靜，黃廣靈，劉 達(dá)

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

0 引 言

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

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

圖1 研究區(qū)域位置

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

1 數(shù)學(xué)模型與工程方案

1.1 基本方程

1.1.1 潮流數(shù)學(xué)模型

連續(xù)方程:

(1)

動量方程:

(2)

(3)

(4)

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

1.1.2 波浪數(shù)學(xué)模型

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

(5)

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

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

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

1.1.3 泥沙數(shù)學(xué)模型

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

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

(6)

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

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

(7)

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

(8)

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

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

(9)

(10)

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

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

(11)

(12)

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

1.2 數(shù)學(xué)模型構(gòu)建

1.2.1 計算范圍

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

圖2 模型網(wǎng)格布置

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

1.2.2 邊界選取

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

表1 波浪模型計算邊界

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

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

1.2.3 模型參數(shù)取值

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

(2)泥沙計算參數(shù)。

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

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

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

1.3 數(shù)學(xué)模型驗證

潮流數(shù)學(xué)模型采用2013年夏季的水文資料進(jìn)行驗證，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以內(nèi)(見表3)；從含沙量驗證情況看(見圖4)，計算含沙量與實測含沙量在總體趨勢和峰值的波動變化與實測含沙量過程較為同步且吻合較好。驗證結(jié)果表明該模型可用于本工程海域水動力與泥沙模擬。

表2 潮位誤差統(tǒng)計 m

表3 流速誤差統(tǒng)計 m/s

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

1.4 工程方案

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

2 計算結(jié)果分析

2.1 工程擴(kuò)建前、后流場變化

2.1.1 流態(tài)變化

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

圖5 工程擴(kuò)建前后流場對比

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

2.1.2 流速變化

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

圖6 工程擴(kuò)建前后流速變化等值線

2.2 年回淤分析

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

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

圖7 波流共同作用下地形沖淤結(jié)果

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

表4 工程特征部位的沖淤強(qiáng)度

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

2.3 驟淤計算

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

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

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

圖8 9713號臺風(fēng)過后地形沖淤狀況

3 結(jié) 論

本文建立工程海區(qū)波流共同作用下的平面二維泥沙數(shù)學(xué)模型，在利用工程海區(qū)實測水沙資料對模型進(jìn)行驗證的基礎(chǔ)上，計算分析了波流聯(lián)合作用下廣東陽西電廠海域潮流泥沙過程及暴風(fēng)驟淤特征，主要得出以下結(jié)論。

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

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

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