沿程非均勻分布植物區(qū)對海嘯波消減的數(shù)值研究

蔣昌波，徐 進，鄧 斌，陳 杰，屈 科

(1.長沙理工大學 水利工程學院，湖南 長沙 410114; 2.水沙科學與水災害防治湖南省重點實驗室，湖南 長沙 410114; 3.洞庭湖水環(huán)境治理與生態(tài)修復湖南省重點實驗室，湖南 長沙 410114)

海嘯引起的極端波況容易造成巨額的財產損失和眾多的人員傷亡[1-2]。2004年的印尼海嘯直接導致了超過23萬人的傷亡[3]，12萬棟房屋被摧毀[4]，經濟損失高達140億美元[5]。2011年日本海嘯的最大波爬高接近40 m，也造成1.8萬人傷亡[6]，導致300座橋梁和16.2萬棟房屋的損壞[7]。此外，隨著南中國海的大規(guī)模開發(fā)，南中國海海嘯的研究也成為熱點[8-9]。為了防范海嘯的破壞作用，人們在實踐中采用了多種工程措施，如建設防波堤和海墻。這些海岸保護措施，一定程度上起到了防浪作用[10]，但也會直接或者間接破壞海岸輸沙平衡，致使岸灘演變加劇[11]，影響海岸的美觀。而海岸植物作為一種自然屏障，既有利于提高海灣的綠化，也可以有效的減少海嘯波的強度，降低海嘯波對海岸區(qū)域的破壞[12-13]。隨著“藍色海灣整治行動”的持續(xù)推進，工程與生態(tài)措施協(xié)調發(fā)展成為關注的熱點。因此，學者們日益關注沿海植物對海岸區(qū)域的保護機理[14]。

過去幾十年，學者們通過理論分析[15-16]、物理試驗[17-19]和數(shù)值模擬[20-22]開展了大量植物消波作用的研究工作，結果表明植物的消波特性主要受到水動力條件和植物分布特征的影響。以往學者們多將植物概化成為均勻分布的圓柱群，研究不同的入射波要素條件下均勻植物的消波特性。事實上，海岸帶植物分布具有非均勻性，這種植物密度分布的非均勻性對植物消波特性影響方面的研究較少，非均勻植物區(qū)中的消減機理研究還不夠。有學者通過現(xiàn)場觀測[23-26]研究植物分布的非均勻性對植物消波特性的影響，但其高成本和復雜的現(xiàn)場條件在很大程度上制約研究的深入開展[27]。植物區(qū)的不均勻性主要分為沿植物區(qū)寬度方向(水平方向)的不均勻性和植物區(qū)高度方向(豎直方向)的不均勻性。Iimura等[28]在物理模型試驗中組合兩種不同水平分布密度的植物(密度大的在前或者密度大的在后)，發(fā)現(xiàn)孤立波經過不同分布的植物區(qū)后波高和流速接近。陳杰等[18]在試驗中將紅樹林概化成根、莖、葉三層，研究豎直方向上不同分布對植物消波的影響。但是受到水槽尺度和試驗模型布置等因素的限制，試驗研究中的植物分布密度變化范圍有限，難以系統(tǒng)揭示植物密度分布非均勻性對植物消波特性的影響。

近些年雖然有學者認為孤立波不能準確地模擬海嘯波，但是其在傳播中能夠維持穩(wěn)定的波形以及與長時間傳播的海嘯波有相似的形態(tài)[29]，因此目前仍然有許多學者用孤立波模擬海嘯波。本文在研究非均勻分布剛性植物對海嘯波的消減作用時使用孤立波代替海嘯波。采用非靜壓單相流模型NHWAVE計算不同水平密度分布植物與孤立波的相互作用，研究它們消波特性的差異性。對比分析不同植物密度分布情況下反射波、波高和最大流速沿程衰減，并與均勻分布植物進行對比。

1 理論模型

1.1 控制方程

NHWAVE控制方程為σ坐標系下的連續(xù)性方程和動量方程，如式(1)和(2)所示：

(1)

(2)

采用k-ε湍流模型，其控制方程為：

(3)

(4)

1.2 拖曳力系數(shù)

由于研究非均勻分布植物對海嘯的消波特性，需要根據不同的植物分布密度選取相應的拖曳力系數(shù)。根據Huang等[31]的試驗研究，雷諾數(shù)Re的范圍1.72×103～5.74×103，與文中計算工況的雷諾數(shù)處于同一區(qū)間，可將Huang等[31]得到的拖曳力系數(shù)經驗公式添加到模型中，其具體表達式為：

CD=1.245+4.587φi/(1-φi)

(5)

2 模型驗證

為了驗證數(shù)值模型計算植物消波的準確性以及拖曳力系數(shù)經驗公式在模型中的適用性，計算了孤立波和非淹沒剛性植物的相互作用，并且與Huang等[31]的試驗數(shù)據進行了詳細對比分析，由于文獻中流速數(shù)據的缺失，因此本文只對自由液面進行了驗證。Huang等[31]的試驗在長32 m、寬0.55 m的試驗水槽中進行，植物區(qū)布置在水槽的中間。剛性植物由有機玻璃管制成，單株植物外徑為0.01 m，試驗中設置了A、B、C三種排列方式，如圖1(a)所示，這三種分布方式所對應的植物分布密度φ分別為0.175、0.087、0.044，與此對應的單位面積植物株數(shù)N分別為2 228、1 108和560。在試驗中共設置了三種不同的植物區(qū)長度W，分別為0.545 m、1.090 m和1.635 m，具體的試驗布置以及波高儀分布位置如圖1(b)所示。

圖1 植物區(qū)三種排列方式和波高儀安放的位置

文中共驗證5組不同算例，所有工況的水深h=0.15 m，波高H分別為0.03 m、0.04 m和0.041 7 m，計算工況如表1所示。建立二維數(shù)值水槽，計算區(qū)域長度為15 m，水平網格大小0.01 m，所有的工況水深均為0.15 m。為了確定豎直網格的數(shù)量，對豎直網格的收斂性進行驗證，如圖2所示。圖2分別展示了A1工況下，測點G1和G5處分別在豎直方向設置10、15和20層網格時的水位對比，可認為豎直網格在20層時已收斂。驗證結果如圖3所示。圖3(a)中G1數(shù)據列為G1波高儀測得水位的時間序列，其第一個峰值為入射波波高，第二個峰值為植物區(qū)的反射波，而G5為波高儀G5測得孤立波經過植物區(qū)的透射波數(shù)據列。由于工況B2中植物區(qū)密度較小，所以反射波也較小。比對圖3(a)和(b)中的反射波，圖3(c)植物區(qū)長度分別為0.545 m、1.090 m和1.635 m時的相對波高沿程變化，可看出計算結果與試驗測量結果整體吻合度較高，表明數(shù)值模型可以有效模擬植物區(qū)的波浪傳播變形。通過上述5組算例的驗證可以看出Huang等[31]提出的經驗公式對于文中數(shù)值模型是適用的，使用該經驗公式算出來的拖曳力系數(shù)可以準確地模擬植物區(qū)對孤立波的反射、透射以及孤立波波高在植物區(qū)中沿程波高的衰減趨勢，因此，使用Huang等[31]提出的經驗公式(5)計算拖曳力系數(shù)CD是可靠的，可以用于研究不同密度分布植物對孤立波的消波作用。

表1 驗證的工況

圖2 網格收斂性驗證

圖3 工況A1、B2自由液面驗證和C1、C2、C3沿程相對波高驗證

3 工況設置

為研究不同密度分布植物對孤立波的消波特性，文中建立了二維數(shù)值水槽，計算區(qū)域網格設置與驗證算例相同，所有的工況水深均為0.15 m，左側為造波邊界，右側設置有長度為6 m的消波區(qū)，波高H設置為0.02～0.05 m。為滿足計算穩(wěn)定性的要求，計算均采用自適應時間步長，最大CFL數(shù)設置為0.2。

剛性植物區(qū)布置在X0=6.0 m和Xn(Xn=X0+W，W為植物區(qū)的長度)之間，其中植物區(qū)長度W設置為0.545 m、1.090 m和1.635 m，測點布置如圖1(b)所示。在植物區(qū)設置了四種不同的密度線性變化方式，分別是：1)L1，植物密度從X0處的N0增大到Xn處的Nn；2)L2，植物密度從X0處的Nn減小到Xn處的N0；3)L3，植物密度在X0和X0+W/2之間由N0增大到Nn，在X0+W/2和Xn之間由Nn減小到N0；4)L4，植物密度在X0和X0+W/2之間由Nn減小到N0，在X0+W/2和Xn之間由N0增大到Nn。其中N0=560 m-2，Nn=2 228 m-2。圖4(a)以W=1.090 m為例分別展示了L1、L2、L3、L4植物密度N沿程變化，圖4(b)為沿程拖曳力系數(shù)，橫坐標x為沿水平正方向與G3測點的距離，x=0 m即為G3所在的位置，所有計算工況如表2所示。

圖4 W=1.090 m時植物區(qū)密度和拖曳力系數(shù)沿程變化

表2 計算工況表

4 結果分析

4.1 反射波對植物分布密度的響應機制

圖5展示了植物區(qū)長度W=1.635 m時植物區(qū)前的反射系數(shù)隨波高的變化，圖中H1為G1所測波高，即入射波高，h為水深，CR為反射系數(shù)。圖5反映出不同分布密度，植物區(qū)波浪的反射效應并不相同，反射系數(shù)從大到小依次是L2、L4、均勻分布、L3和L1，說明植物越集中在前面部分反射波越大。其中L1與L2的差值最大，H1/h=0.13時L2的反射系數(shù)比L1大80.2%，H1/h=0.33時L2的反射系數(shù)比L1大78.3%。從圖5中還看出所有分布的反射系數(shù)均與入射波高線性相關，隨著入射波高的增大而增大。

圖5 W=1.635 m時反射系數(shù)隨波高變化

圖6 W=1.635 m、H1/h=0.33時L1和L2反射波隨時間的變化

4.2 波高衰減對植物分布密度的響應機制

4.2.1 入射波高的影響

圖7展示了W=1.090 m，入射波高H1/h分別為0.13、0.20、0.27、0.33時孤立波在植物區(qū)沿程波高衰減。橫坐標x為沿水平正方向與G3測點的距離，x=0 m即為G3所在的位置。對于不同的波高，孤立波在植物區(qū)中沿程消減規(guī)律相近，L1與L2的曲線組成類“0”字形，L3和L4的曲線組成類“8”字形，“8”字形的交點均在x/h=4.10～4.20之間，而均勻分布的曲線在“0”字和“8”字的中間。造成這種現(xiàn)象的原因是，L1和L2在植物區(qū)密度分布相反，L1的密度線性增大，在植物區(qū)最前端密度最小，所以植物區(qū)前端波高衰減的強度弱，而L2的密度線性減小，在植物區(qū)最前端密度最大，所以植物區(qū)前端波高衰減的強度大，但是隨著密度沿程的變化L1波高衰減的強度逐漸大于L2，最終在植物區(qū)尾端波高之間的差小于1.5%。L3在植物區(qū)中點處密度最大，向兩端線性減小，消波效率先增大后減小，而L4與此相反，在中點處密度最小，向兩端線性增大，所以在植物區(qū)中點后會得到相同的波高即兩條曲線的交點，最終L3、L4、L1、L2、均勻分布的透射波高十分接近，其之間的差距小于2.5%。

圖8展示了W=0.545 m、1.090 m和1.635 m時透射系數(shù)隨波高的變化，從圖中可看出透射系數(shù)隨著入射波高的增大而減小。而L1、L2、L3、L4和均勻分布的透射系數(shù)隨波高變化的趨勢是一致的，因此在植物區(qū)長度和波高相同的情況下5種分布的透射系數(shù)相同。

圖7 W=1.090 m時波高沿程變化

圖8 透射系數(shù)隨波高的變化

4.2.2 植物區(qū)長度的影響

圖9(a)、圖7(d)、圖9(b)分布展示了H1/h=0.33，植物區(qū)長度分別為0.545 m、1.090 m、1.635 m時植物區(qū)中孤立波沿程波高。從圖中可看出隨著植物帶長度的增大，植物區(qū)中波高沿程消減的機制沒有明顯的變化，L1、L2、L3和L4分布的沿程曲線仍然呈現(xiàn)“0”字形和“8”字形，而均勻分布的曲線位于它們中間?？梢园l(fā)現(xiàn)W=0.545 m和W=1.090 m時L3和L4的交點分別在x/h=2.08和x/h=6.18處，皆在植物區(qū)中部偏后的位置。

圖10展示了H1/h=0.33時，透射系數(shù)隨植物區(qū)長度的變化，圖中橫坐標W為植物區(qū)長度。從圖中可看出透射系數(shù)均隨著植物區(qū)長度的增大而減小，但是4種非均勻分布和均勻分布植物的透射系數(shù)隨植物區(qū)長度變化的規(guī)律相同，因此出現(xiàn)了當植物區(qū)長度相同時5種分布具有相近透射系數(shù)的情況。

圖9 H1/h=0.33時波高沿程變化

圖10 H1/h=0.33時透射系數(shù)隨植物區(qū)長度的變化

4.3 植物區(qū)中流速對植物分布密度的響應機制

4.3.1 入射波高的影響

圖11分別展示了W=0.545 m、1.090 m和1.635 m時流速透射系數(shù)UT隨波高的變化。4種非均勻分布和均勻分布植物的流速透射系數(shù)UT均隨著入射波高的增大而減小，例如W=1.090 m時，H1/h=0.33的流速透射系數(shù)只有H1/h=0.13的64.1%。然而5種分布的流速透射系數(shù)在隨著波高變化的過程中保持一致性，即UT隨波高變化呈現(xiàn)出一樣的規(guī)律，在植物區(qū)長度和波高相同時具有相近的流速透射系數(shù)。

圖11 流速透射系數(shù)隨波高的變化

圖12展示了植物區(qū)長度W=1.090 m，入射波高H1/h分別為0.13、0.20、0.27、0.33時，植物區(qū)沿程最大流速的變化，圖中的Umax,1為入射波的最大流速，取波峰經過斷面時該斷面上的最大流速。從圖中可看出沿程最大流速的曲線也呈現(xiàn)出類似的“0”字形和“8”字形，只是其在最前端出現(xiàn)開口，這與上述的沿程波高曲線不同。出現(xiàn)不同的原因是圖7和圖9中的波高在x=0處由于反射等原因會使得瞬時波高大于入射波高，而且x=0處的密度越大瞬時波高越大，因此會在波高衰減的過程中出現(xiàn)兩兩交叉的情況。但是最大流速則不同，在x=0處的植物密度越大會使得最大流速越小，所以在消減過程中并不會出現(xiàn)交叉。通過圖12(a)～(d)的對比可發(fā)現(xiàn)入射波高的變化，對最大流速的沿程變化規(guī)律沒有明顯的影響。

圖12 W=1.090 m時沿程最大流速變化

4.3.2 植物區(qū)長度的影響

圖13(a)、圖12(d)、圖13(b)分別展示了H1/h=0.33，植物區(qū)長度分別為0.545 m、1.090 m、1.635 m時植物區(qū)沿程最大流速的變化。從圖中可看出對于不同的植物區(qū)長度，最大流速沿程消減的機制沒有明顯的變化。W=0.545 m、1.090 m和1.635 m的情況下，L1和L2的曲線呈現(xiàn)出開口的“0”字形，而L3和L4的曲線組成開口的“8”字形，均勻分布位于它們中間，而且最終經過植物區(qū)后最大流速相近，最大差為5.4%。

圖14展示了H1/h=0.33時，UT隨植物區(qū)長度的變化。從圖中可看出，UT均隨著植物區(qū)長度的增大而減小，例如W=1.635 m時的流速透射系數(shù)只有W=0.545 m時的51.2%。然而4種非均勻分布和均勻分布的流速透射系數(shù)隨植物區(qū)長度的變化規(guī)律是一致的，在植物區(qū)長度相同的情況下5種分布具有相近的流速透射系數(shù)。

圖13 H1/h=0.33時最大流速沿程變化

圖14 H1/h=0.33時，流速透射系數(shù)隨植物區(qū)長度的變化

5 結 語

采用非靜壓數(shù)值模型開展了非均勻剛性非淹沒植物對海嘯波消減的數(shù)值模擬工作。通過與Huang等[31]的試驗數(shù)據進行對比驗證數(shù)值模型的計算可靠性。文中研究了5種不同密度分布植物與孤立波的相互作用，分析了不同密度分布植物對海嘯的消波特性。通過對計算結果進行分析，主要結論如下：

1)在植物總株數(shù)相同的情況下，植物區(qū)的非均勻性對孤立波的消減作用幾乎沒有影響，只改變了沿程的消波強度分布。對于平均密度相同的5種分布，孤立波波高和流速在不同分布密度的植物中消減過程不同，L1與L2沿程波高和最大流速曲線組成類“0”字形，L3和L4沿程波高和最大流速曲線呈類“8”字形。不同分布密度的植物區(qū)沿程消波強度分布不同，但最終透射波的波高和流速幾乎相同，不同分布的透射波高之差小于3%，最大流速的差小于5.4%。在所有的工況下，孤立波的透射波高和流速的差皆很小，說明透射波的能量相差較小，可以說在植物總株數(shù)相同的情況下，不同分布密度的植物對孤立波的消減作用是相同的。

2)孤立波的入射波高和植物區(qū)的長度變化，不改變波浪在非均勻分布植物區(qū)沿程消減機制，波高和流速沿程衰減曲線皆呈現(xiàn)出類“0”字形和類“8”字形。而且孤立波經過不同分布植物區(qū)的透射系數(shù)和流速透射系數(shù)隨著入射波高和植物區(qū)長度變化的規(guī)律是一致的。

3)通過研究可見，在植物總株數(shù)相同的情況下，植物區(qū)前端密度較大時會引起海嘯波反射增強，但卻對最終的海嘯波消減效率影響有限。因此在實際工程實踐中，在植物總分布密度接近的情況下，可使植物帶密度隨與海岸的距離逐漸減小。這樣不但可以有效的消減海嘯波的破壞作用，并且可以進一步減少海嘯波反射引起的次生災害。