陳漢寶，鐘 生，陳松貴*，時 健

(1.交通運輸部天津水運工程科學(xué)研究所，天津 300456；2.河海大學(xué) 港口海岸與近海工程學(xué)院，南京 210098)

珊瑚礁主要是由珊瑚蟲、鈣化藻、有孔蟲等造礁生物的灰質(zhì)骨骼殘體世代堆積形成。我國的珊瑚礁主要分布于中國大陸南方沿岸、海南島和臺灣島[1]，以及南海的西沙、東沙、中沙和南沙群島。珊瑚礁生態(tài)系統(tǒng)為人類社會提供了多種多樣的服務(wù)，發(fā)揮著重要的作用，如海洋生物資源、漁民生計、建筑材料和海洋藥物以及文化旅游娛樂等。因此，如何在維持原有生態(tài)環(huán)境情況下提高珊瑚礁系統(tǒng)的多樣性，就需要對珊瑚礁地形上的水動力特性有充分的了解和研究[2]。

由于珊瑚礁地形是一種坡度急劇變化的陡變地形[3]，在此類復(fù)雜地形條件下，研究波浪在珊瑚島礁地形上的波高、流速、增減水等變化規(guī)律，對珊瑚礁系統(tǒng)的生態(tài)環(huán)境保護和利用起著至關(guān)重要的作用。目前針對珊瑚島礁地形下的波流研究主要方法包括：現(xiàn)場觀測、物理模型試驗和數(shù)值模擬等。對于現(xiàn)場觀測而言，通常只能針對特定區(qū)域或測站，并且投入時間長、資金多，同時極易受到自然因素的影響。另外，物模試驗的試驗設(shè)計以及試驗開展需要投入較多的人力、物力和財力等，還可能會受到場地限制和試驗材料等因素的影響，并且模擬的試驗工況也會相對受限，特別是極端工況下的水動力模擬非常困難。為了獲取更多類別及工況的試驗數(shù)據(jù)，開展數(shù)值模擬試驗既能解決以上所面臨的問題，還能夠節(jié)約大量的試驗成本。

隨著計算機性能的不斷提高，計算方法的不斷進步，促進了波浪數(shù)值模擬技術(shù)的研究與應(yīng)用的發(fā)展[4]。其中，非靜壓波浪模型得到快速發(fā)展，其中基于非線性淺水方程的非靜壓波浪模型SWASH，采用單值函數(shù)對自由表面進行追蹤[5]，可以對較為復(fù)雜地形下的波面和淺水流進行大范圍、長時間的模擬。并且，Suzuki等[6]通過SWASH模擬結(jié)果與試驗數(shù)據(jù)的對比，驗證SWASH模型在波浪變形和越浪計算的適用性；Guimaraes等[7]通過SWASH對巴西Tramandaí海灘波浪進行了模擬，并對其爬坡模擬值進行了驗證，得出SWASH可以很好地模擬該海域的波浪狀況。白志剛[8]對基于Boussinesq方程的FUNWAVE模型和基于非線性淺水方程的SWASH模型從計算精度和運行速度等方面進行分析對比，得出SWASH模型比FUNWAVE模型更具有實用性；王良才等[9]將SWASH模型模擬結(jié)果與物理模型的實測數(shù)據(jù)進行對比后發(fā)現(xiàn)，SWASH模型可以對近海岸及港口區(qū)域的多種物理現(xiàn)象如淺水變形、繞射、折射、反射、破碎等現(xiàn)象進行有效的模擬。不過，諸多涉及SWASH模型的文章并沒有特別關(guān)注SWASH模型模擬下空間步長和時間步長的變化對試驗結(jié)果的影響，因此本文將采用SWASH模型對島礁進行三維模擬，并選取不同的空間步長和時間步長模型模擬的試驗數(shù)據(jù)進行SWASH模型的時空敏感性研究。

1 非靜壓模型SWASH簡介

1.1 控制方程

SWASH (Simulating Waves till Shore)模型由代爾夫特理工大學(xué)的Stelling、Zijlema和Smit等學(xué)者開發(fā)完成的一款淺水非靜壓數(shù)值模型[10-12]。模型考慮了波浪運動中的非靜壓影響，適用于深海至淺水條件下波浪和水流運動的模擬，特別是對于復(fù)雜地形如三維地形下的波浪破碎所引起的波浪增水和波生流有著較好的適用性。特別地，SWASH的控制方程里特別考慮到垂向動量方程和水平動量方程中附加的非靜壓項，在求解時把控制方程分解為靜壓項和非靜壓項兩部分分別求解[13]，其控制方程如下

(1)

(2)

(3)

(4)

式中：t為時間；x、y和z方向的速度分量分別為u、v、w，水平向粘滯系數(shù)為vH，垂向粘滯系數(shù)為νv。ξ代表自由面；g代表重力加速度；q代表非靜水壓力。

1.2 邊界條件

自由面的壓力邊界條件為

q|z=ξ=0

(5)

當(dāng)垂向流速滿足式(5)的邊界條件時水平流速滿足

(6)

在SWASH中關(guān)注入流邊界位置，其邊界條件的設(shè)定，需要明確速度參數(shù)，同時考慮弱反射問題

(7)

式中：ub為入流邊界的流速；邊界入射波的波幅由ξb表示；式(7)中邊界位置是通過正負(fù)號來描述的，+號表示入流邊界為西、南邊界，-號表示入流邊界為東、北邊界。并且采用干濕處理技術(shù)來描述動邊界，模型通過給定最小水深值限制網(wǎng)格點計算過程中出現(xiàn)的負(fù)水深，在計算的單位時間步長內(nèi)對網(wǎng)格點進行干濕判斷，當(dāng)計算的網(wǎng)格點水深低于最小水深時，將最小水深設(shè)為該網(wǎng)格點的水深，避免模型計算中出現(xiàn)負(fù)水深引起模型運行中斷。

1.3 數(shù)值求解

SWASH模型計算采用有限差分格式進行空間離散，垂向采用σ分層網(wǎng)格，變量通過交錯網(wǎng)格方式布置。非靜水壓力q定義在單元邊中心，如圖1所示。

圖1 網(wǎng)格垂向分層和變量布置示意圖Fig.1 Schematic diagram of grid vertical stratification and variable layout

2 物理模型簡介及數(shù)模參數(shù)設(shè)置

2.1 物理模型簡介

圖2 模型尺寸及試驗儀器設(shè)置Fig.2 Model size and test instrument setting

整體物理模型建于交通運輸部天津水運工程科學(xué)研究院綜合試驗廳的波浪水池中，水池長60 m，寬42 m，深1 m，水池一側(cè)的造波系統(tǒng)由9臺推板式造波機和微機控制系統(tǒng)組成。在距造波機34 m處設(shè)置坡度為1:8的斜坡模擬礁前斜面，水平投影長3.2 m，斜面后接長度為14 m的水平平臺模擬礁坪，其頂高程為0.4 m。在兩個礁坪中間存在寬度為6 m的裂口，礁坪后存在寬度為5 m的潟湖，最后在潟湖后設(shè)置坡度為1:3.3的礁后斜坡，其余斜坡坡度均為1:1。物理模型底部用砂石填充，表面進行光滑水泥抹面。試驗中采集波面使用的是LG1型電容式浪高儀，浪高儀量程為0.6 m，絕對誤差小于1 mm，采集單點流速分布使用的是ADVs(超聲多普勒流速儀)和螺旋槳式流速儀。由于該礁坪-潟湖-裂口系統(tǒng)物理模型對稱分布，因此所有的儀器主要布置在裂口的一側(cè)。具體的試驗?zāi)Ｐ统叽绾驮囼灧治鏊玫降膬x器編號位置如圖2所示，其中G3、G5等表示波高儀布設(shè)位置，1#、2#、3#，ADV1、2，LXJ-1、3等表示流速儀布設(shè)位置。

2.2 SWASH模型參數(shù)設(shè)置

根據(jù)鄒國良等[14]的研究，在模擬波浪傳播時，非靜壓模型的水平網(wǎng)格選取ΔX=L/30便可較為精確地模擬波面、流速等結(jié)果。本文選取物理模型一組特征工況為：礁坪淹沒水深hr=0.04 m，外海入射波高H0=0.04 m，周期T0=2 s，通過計算可得到波長L=3.847 m，因此選用模型網(wǎng)格尺寸為ΔX=0.1 m，ΔY=0.1 m。波浪入射為垂直左邊界入射，在右邊界處設(shè)置5 m寬的海綿層，模型垂向取4層，糙率通過設(shè)置曼寧粗糙系數(shù)實現(xiàn)。考慮到計算網(wǎng)格尺寸會影響到模型的空間分辨率，因此本文以模型網(wǎng)格尺寸為上限，選取精度介于Δx=L/100～L/30間的七組不同計算網(wǎng)格尺寸0.04 m、0.05 m、0.06 m、0.07 m、0.08 m、0.09 m、0.1 m進行模擬研究，并將模擬數(shù)據(jù)進行分析對比。此外，SWASH使用手冊提到SWASH中的時間步長應(yīng)足夠小，以解決計算出波場本身的時間變化，因此選取介于0.005～0.03 s的0.005 s、0.01 s、0.02 s、0.03 s四組不同計算時間步長進行模擬研究，以此來分析模型值隨計算網(wǎng)格尺寸和計算時間步長的敏感性變化。

3 敏感性分析及方法簡介

本文主要選取與物模流速測點和礁坪中線對應(yīng)位置的平均流速和有效波高進行敏感性分析，現(xiàn)針對本文數(shù)學(xué)模型所用的敏感度分析法做如下簡介：

3.1 流速隨計算網(wǎng)格尺寸的敏感性變化

圖3展示了不同計算網(wǎng)格尺寸、計算時間步長均為△t=0.01s下裂口中線和礁坪中線在x方向各測點數(shù)值模擬和物理模型的向岸平均流速，平均流速取造波穩(wěn)定后30個周期內(nèi)的平均值。通過對比可以發(fā)現(xiàn)，當(dāng)△x在0.04 m附近時，這三個測點位置向岸平均流速吻合度最高。對比礁坪上測點的向岸平均流速也可以發(fā)現(xiàn)，各個測點的數(shù)值模擬值與實測值吻合度也較好，在礁坪上數(shù)值模擬值與實測值有微小差距，分析原因在于，波浪傳播至礁緣處時，波浪發(fā)生破碎，波浪形態(tài)和流速發(fā)生改變，并且實際物理模型與數(shù)值模擬的地形條件不可能完全相同，因此模擬結(jié)果也會存在相應(yīng)的差異，不過其最大誤差仍在允許范圍內(nèi)。

3-a 裂口中線測點平均流速 3-b 礁坪中線測點平均流速圖3 不同△x下裂口和礁坪中線測點平均流速Fig.3 Average flow velocity at the measuring points in the middle line of crack and reef flat under different △x

表1 不同△x下各測點向岸平均流速Tab.1 Average velocity of each measuring point towards the shore under different △x m/s

表2 不同△x下各測點向岸平均流速敏感度Tab.2 Sensitivity of shore average velocity of each measuring point under different △x %

3.2 波高隨計算網(wǎng)格尺寸的敏感性變化

圖4展示了不同計算網(wǎng)格尺寸、計算時間步長均為△t=0.01 s下裂口中心線和礁坪中心線在x方向各測點數(shù)值模擬和物理模型的有效波高，波高選取造波穩(wěn)定后30個周期內(nèi)的波高值。通過對比可以發(fā)現(xiàn)，在裂口通道G39及以前各測點數(shù)值模擬的有效波高值與實測值吻合度較好，在潟湖內(nèi)和靠近礁后斜坡處G40、G41測點的有效波高值隨計算網(wǎng)格尺寸的改變也發(fā)生變化，并且其影響不可忽略。此外，通過對比礁坪上測點的有效波高值也可以發(fā)現(xiàn)，自外海向潟湖內(nèi)部，各個測點的有效波高值與實測值吻合度也較好，礁坪中部G12測點數(shù)值模擬值與實測值有微小差距，在靠近礁后斜坡G15測點的有效波高值隨著計算尺寸的改變也發(fā)生變化，并且其影響也不可忽略，分析原因在于波浪發(fā)生破碎后，波流形態(tài)發(fā)生紊動，并且礁后斜坡處存在強烈的反射現(xiàn)象，因此物理模型和數(shù)值模擬的流速和波高在靠近礁后斜坡處附近存在相應(yīng)的差異，不過可以發(fā)現(xiàn)，當(dāng)△x在0.05 m附近時，各測點的有效波高值與實際物理模型值吻合度較好。

4-a 裂口中線測點有效波高 4-b 礁坪中線測點有效波高圖4 不同△x下裂口和礁坪中線測點有效波高Fig.4 Effective wave height at the measuring points in the middle line of the crack and reef flat under different △x

表3是不同計算網(wǎng)格尺寸、計算時間步長均為△t=0.01 s時各測點有效波高值，波高選取造波穩(wěn)定后30個周期內(nèi)的波高值。以最高精度計算網(wǎng)格尺寸△x=0.04 m為基準(zhǔn)組，應(yīng)用敏感性分析方法對不同計算網(wǎng)格精度下所有測點的有效波高值敏感度S(Hsi,xi)=(△Hsi/△xi)/(Hsi,xi)進行計算，其計算結(jié)果如表4所示，其中：△Hsi表示每組計算網(wǎng)格尺寸下的有效波高相對于基準(zhǔn)組△x=0.04 m的有效波高的差值；△xi表示每組計算網(wǎng)格尺寸與基準(zhǔn)組△x=0.04 m的計算網(wǎng)格尺寸的差值；Hsi表示每組計算網(wǎng)格尺寸下對應(yīng)測點各自的有效波高值；xi表示基準(zhǔn)組參數(shù)，此處為基準(zhǔn)計算網(wǎng)格尺寸△x=0.04 m。由表4可以發(fā)現(xiàn)，裂口中線測點G40以前的各測點位置的敏感度都小于10%，潟湖內(nèi)測點G40和靠近礁后斜坡的測點G41的有效波高值對計算網(wǎng)格尺寸的變化表現(xiàn)得比較敏感，再對比礁坪中線上有效波高值也可以發(fā)現(xiàn)，測點G12以前的各測點位置的敏感度幾乎都小于10%，測點G12至測點G15，有效波高值也對計算網(wǎng)格尺寸變得比較敏感。分析原因，在復(fù)雜三維地形下，低精度網(wǎng)格尺寸在輸出時因讀取不到相鄰網(wǎng)格點間某些較大的數(shù)值而導(dǎo)致結(jié)果偏小，因此利用非靜壓模型SWASH模擬近岸復(fù)雜地形時，計算網(wǎng)格精度是一個不可忽略的重要因素，高精度的計算網(wǎng)格對模型的準(zhǔn)確性起著至關(guān)重要的作用。

表3 不同△x下各測點有效波高值Tab.3 Effective wave height of each measuring point under different △x m

表4 不同Δx下各測點有效波高值敏感度Tab.4 Sensitivity of effective wave height value of each measuring point under different Δx %

3.3 流速隨計算時間步長的敏感性變化

圖5展示了不同時間步長、計算網(wǎng)格尺寸均為Δx=0.04 m下裂口中心線和礁坪中心線在x方向各測點數(shù)值模型的向岸平均流速。通過對比可以發(fā)現(xiàn)，改變時間步長，數(shù)值模型的流速值沒有明顯的變化，其影響可忽略不計。

5-a 裂口中線測點平均流速 5-b 礁坪中線測點平均流速圖5 不同Δt下裂口和礁坪中線測點平均流速Fig.5 Average flow velocity at the measuring points in the middle line of crack and reef flat under different Δt

表5 不同Δt下各測點向岸平均流速Tab.5 Average velocity of each measuring point towards the shore under different Δt m/s

表6 不同Δt下各測點向岸平均流速敏感度Tab.6 Sensitivity of shore average velocity of each measuring point under different Δt %

3.4 波高隨計算時間步長的敏感性變化

圖6展示了不同時間步長、計算網(wǎng)格尺寸均為Δx=0.04 m下裂口中心線和礁坪中心線在x方向各測點數(shù)值模型的有效波高值。通過對比可以發(fā)現(xiàn)，時間步長對有效波高值的影響程度很低，在計算時間步長為0.005～0.03 s，改變時間步長，對模型輸出的有效波高值幾乎沒有影響。

6-a 裂口中線測點有效波高 6-b 礁坪中線測點有效波高圖6 不同Δt下裂口和礁坪中線測點有效波高Fig.6 Effective wave height at the measuring points in the middle line of the crack and reef flat under different Δt

表7是不同時間步長、計算網(wǎng)格尺寸均為Δx=0.04 m時各測點有效波高值，波高選取造波穩(wěn)定后30個周期內(nèi)的波高值。以最高精度時間步長Δt=0.005 s為基準(zhǔn)組，應(yīng)用上述方法對不同計算網(wǎng)格精度下裂口和礁坪中線上測點的有效波高值敏感度進行計算S(Hsi，xi)=(ΔHsi/Δxi)/(Hsi/xi)，其計算結(jié)果如表8所示，其中：ΔHsi表示每組計算時間步長下的有效波高相對于基準(zhǔn)組Δt=0.005 s的有效波高的差值；Δxi表示每組計算時間步長與基準(zhǔn)組Δt=0.005 s的計算時間步長的差值；Hsi表示每組計算時間步長下對應(yīng)測點各自的有效波高值；xi表示基準(zhǔn)組參數(shù)，此處為基準(zhǔn)計算時間步長Δt=0.005 s。由表8可以發(fā)現(xiàn)，裂口中線測點以及礁坪中線上測點的有效波高值敏感度均小于10%，說明在計算時間步長0.005～0.03 s對模型的有效波高值影響并不顯著，其影響可以忽略不計。

表7 不同Δt下各測點有效波高值Tab.7 Effective wave height of each measuring point under different Δt m

表8 不同Δt下各測點有效波高值敏感度Tab.8 Sensitivity of effective wave height value of each measuring point under different Δt %

4 結(jié)論

采用非靜壓時域波浪模型SWASH，研究了特征規(guī)則波工況入射條件下，珊瑚島礁水域內(nèi)的平均流速和有效波高隨模型中的計算網(wǎng)格尺寸和計算時間步長的變化情況，也可稱之為時空敏感性分析。綜合以上分析表明：SWASH可以較好地模擬近岸復(fù)雜地形下的波流特性；數(shù)值模擬結(jié)果中流速隨計算網(wǎng)格尺寸的變化敏感度較低，有效波高值隨計算網(wǎng)格尺寸的變化表現(xiàn)得較為敏感，特別是當(dāng)波浪發(fā)生破碎變形之后，其敏感程度較高；數(shù)值模擬結(jié)果中的流速和波高隨計算時間步長的變化表現(xiàn)得并不敏感，其影響可以忽略不計。因此，采用SWASH模型模擬珊瑚礁水動力時，在高精度計算時間步長范圍內(nèi)，只需選取合適的計算網(wǎng)格尺寸便可以較好地模擬復(fù)雜地形下的波流傳播特性。