島礁三維SWASH模型的時(shí)空敏感性分析

陳漢寶，鐘 生，陳松貴*，時(shí) 健

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

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

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

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

1 非靜壓模型SWASH簡(jiǎn)介

1.1 控制方程

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

(1)

(2)

(3)

(4)

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

1.2 邊界條件

自由面的壓力邊界條件為

q|z=ξ=0

(5)

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

(6)

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

(7)

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

1.3 數(shù)值求解

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

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

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

2.1 物理模型簡(jiǎn)介

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

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

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

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

3 敏感性分析及方法簡(jiǎn)介

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

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

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

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

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

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

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

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

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

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

表3 不同△x下各測(cè)點(diǎn)有效波高值Tab.3 Effective wave height of each measuring point under different △x m

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

3.3 流速隨計(jì)算時(shí)間步長(zhǎng)的敏感性變化

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

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

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

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

3.4 波高隨計(jì)算時(shí)間步長(zhǎng)的敏感性變化

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

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

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

表7 不同Δt下各測(cè)點(diǎn)有效波高值Tab.7 Effective wave height of each measuring point under different Δt m

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

4 結(jié)論

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