黃飛揚,戴文鴻,姚 毓,陳羿名
(河海大學(xué) 水利水電學(xué)院,江蘇 南京 210098)
波浪溢流[1]是指沿海地區(qū)極端高水位和熱帶氣旋增加產(chǎn)生的風暴潮導(dǎo)致海堤堤前水位超過堤頂,海堤受到越浪與溢流聯(lián)合作用的現(xiàn)象[2]。在全球氣候變暖的背景下,極端氣候事件發(fā)生的強度和頻率都在不斷增加[3],這使得發(fā)生波浪溢流現(xiàn)象的概率也進一步提高。波浪溢流一旦發(fā)生,將會極易引起大規(guī)模的潰堤災(zāi)難[4]。2005年卡特里娜颶風風暴潮期間,美國墨西哥灣沿岸由于波浪溢流產(chǎn)生的海堤后坡侵蝕破壞就多達43處[5],大量的村莊和農(nóng)田被海水沖毀。自此波浪溢流現(xiàn)象受到了海岸工程界學(xué)者的高度重視,波浪溢流也成為了海堤設(shè)計和安全評估中必須考慮的因素。
自2005年卡特里娜颶風以來,波浪溢流產(chǎn)生的潰壩潰堤災(zāi)害受到了整個海岸工程界的高度關(guān)注,國內(nèi)外學(xué)者也對此做了大量的研究。如歐洲的越浪手冊(第一版)[6]、Reeve等[7]、Hughes等[8]和Li等[9]都圍繞著波浪溢流越堤流量參數(shù)對其水動力特征進行了分析計算,然而上述研究僅僅較為宏觀地反映出海堤對波浪溢流水動力特征的影響,對越堤水流的形態(tài)、空間變化的相關(guān)描述卻極為缺乏。為了進一步研究波浪溢流對海堤不同時間、空間分區(qū)的影響,可以通過使用計算流體軟件對海堤各重要位置的水動力特征進行模擬。隨著計算機技術(shù)發(fā)展,基于FLUENT軟件模擬達到工程實用水平的數(shù)值波浪水槽已經(jīng)成為了國際海岸與近海工程界關(guān)注的熱點。如周勤俊等[10]利用FLUENT軟件將模擬的入射波作為人工分布源項加入動量方程之中,發(fā)展得到了適用于流體體積(VOF)方法的源項造波法;黃華等[11]針對波浪及其影響效應(yīng)的研究需求,采用了推波板造波方法對典型二維Stokes波產(chǎn)生、發(fā)展和消亡過程進行了準確的仿真研究;徐剛等[12]通過在入口處添加速度場來模擬波浪在受數(shù)值因素影響下所造成的衰減,并得到數(shù)值仿真精度更高的模擬方案。通過上述研究,針對波浪溢流水動力變化過程的特點,將采用UDF速度邊界造波方法對規(guī)則波中的二階Stokes波進行模擬。該方法主要通過對入口水質(zhì)點輸入自定義速度獲得相應(yīng)的波浪波形,經(jīng)驗證輸出的波形較其他方法十分準確,同時模型計算量急劇減小,使得數(shù)值仿真精度和效率都得到了一定的提高[13-14]。阻尼消波法主要對波浪遇建筑物時的消波效果較好[15],同實驗室的消波方法原理相同,加入阻尼源項后,便于對阻尼消波區(qū)域進行控制,使得參數(shù)設(shè)置更加方便。因此,為了防止波形受到海堤反射的作用,可采用阻尼消波方法消除反射波對波浪溢流水動力特性的影響。
綜上所述,將使用UDF速度邊界造波和阻尼消波法對不同超高和波浪要素的波浪溢流現(xiàn)象進行模擬。將模擬結(jié)果與Hughes等[8]和Li等[9]物理模型試驗相關(guān)研究成果進行對比驗證分析,并對波浪溢流中越浪和溢流的主導(dǎo)性關(guān)系進行探討。再進一步通過各測點水動力變化的時間序列研究海堤后坡水流厚度的動態(tài)過程,同時對海堤后坡穩(wěn)定水深公式進行擬合并驗證。
在對波浪溢流進行模擬時,研究的二維波浪屬于不可壓縮黏性流體,即將流動視為包含水和空氣的兩相流進行處理,空氣和水的交界面即為模擬波浪的波面,且二者之間不可滲透。這里取水平方向為x軸(向右為正),取垂直方向為y軸(向上為正)。因此得到的N-S(Navier-Stokes)控制方程可簡化為:
連續(xù)方程的表達式為:
(1)
動量方程的表達式為:
(2)
式中:u和v分別表示為在兩個方向(x、y方向)上速度的分量;P為流體微元體上的壓力;ρ為流體密度;ν為流體動力黏度系數(shù)。
當波浪接近海堤時,由于海堤坡度的影響使得水流產(chǎn)生較大的紊動甚至導(dǎo)致波浪破碎,因此在波浪溢流數(shù)值水槽建立之前應(yīng)對湍流模型進行考慮。RNG k-ε模型是由標準的k-ε模型改進而來,最早由Yakhot等[16]提出,該模型考慮了大面積流體運動和空間位置變化的同時還可以捕捉到流線彎曲程度較大的流動,十分適合對波浪溢流和海堤相互作用的模擬。其湍流模型數(shù)學(xué)表達式為:
(3)
(4)
其中,
μeff=μ+μt
(5)
(6)
式中:Gk由平均速度梯度引起的湍動能產(chǎn)生;αk和αε分別為湍流能和耗散率有效普朗特數(shù)的倒數(shù);C1ε和C2ε為經(jīng)驗參數(shù),取默認值分別為1.42和1.68;由于文中流體運動的雷諾數(shù)較大(Re>30 000),故Cμ取為0.084 5;μeff為湍流流動時的名義黏性系數(shù);μ為水體的分子黏性系數(shù);μt為湍流黏性系數(shù),即名義黏性系數(shù)。模型糙率的選擇將參照Hughes等[8]物理模型,因Hughes模型中海堤模型表面為碾壓混凝土,其曼寧系數(shù)n為0.018,通過FLUENT手冊對模型糙率參數(shù)進行換算,在FLUENT模型中需要輸入的粗糙度常數(shù)參數(shù)Cs取為0.53。
在基于N-S方程建立的數(shù)值水槽中,能夠準確捕捉自由表面的大變形運動是關(guān)鍵。VOF方法適用于有清晰相界面的水體流動,對于處理具有復(fù)雜自由表面的波浪溢流問題十分理想,可以清晰地分析出波浪溢流作用于海堤時的沖擊壓力和流場特征。
在VOF方法中引入了體積函數(shù)aq用于定義第q相流體所占體積與單元體積的比值,則體積函數(shù)aq的值存在以下3種可能:若aq=0時,單元中沒有第q相流體;若aq=1時,單元中充滿了第q相流體;若0 (7) 運用速度邊界造波法對規(guī)則波中的二階Stokes波進行數(shù)值模擬,速度邊界造波法是將入口邊界水質(zhì)點賦予已知的波水平和垂向周期性振蕩速度,促使水質(zhì)點運動產(chǎn)生波浪,并通過波面方程最終確定輸入波浪的波形。此方法可以通過FLUENT UDF進行實現(xiàn)。UDF(用戶自定義函數(shù))通過用戶編制C語言程序連接到FLUENT所提供的預(yù)定義宏來實現(xiàn)特有功能,從而達到對軟件進行二次開發(fā)的目的。利用FLUENT中可以直接訪問單元中流體變量的宏DEFINE PROFILE(name,t,i)來定義入口邊界水質(zhì)點的速度方程和波面方程。對于二階Stokes速度入口邊界處的流速分布方程和波面方程可表示為: 波面方程: (8) x方向速度: (9) (10) 式中:ω為波頻率;k為波數(shù);H為波高;h為靜水位。 波浪和海堤相互作用后,產(chǎn)生的反射波會導(dǎo)致水槽內(nèi)的波浪特性發(fā)生改變,為了消除海堤對波浪的二次反射作用,需要對數(shù)值水槽進行消波處理。阻尼消波法是根據(jù)實驗室水槽消波的原理,在特定的計算區(qū)域間加入阻尼源項,起到對該區(qū)域間反射波減弱和消除的作用。加入的阻尼源項表達式為: (11) 式中:S為y方向動量方程的源項;C1為線性阻尼常數(shù);C2為二次阻尼常數(shù);V為沿y方向的速度;f(y)為沿y方向的阻尼函數(shù);f(x)為沿x方向的阻尼函數(shù)。其中, (12) (13) 式中:xs、xe分別是阻尼層在x方向上的起始位置和終點位置;yfs、ye分別是阻尼層在y方向上的起始位置和終點位置,即自由表面高程和底部高程。 以Hughes等[8]和Li等[9]的大水槽試驗為基礎(chǔ),建立了一個與Hughes和Li物理水槽完全相同的二維數(shù)值波浪水槽模型。其中模型水槽長49.3 m,高2.80 m;海堤模型高0.65 m,前坡坡度為1∶4.25,后坡坡度為1∶3。在海堤模型前還設(shè)置了一個坡度為1∶10、長度為2 m的陡坡,數(shù)值波浪水槽和海堤模型如圖1所示。 圖1 數(shù)值波浪溢流水槽示意Fig.1 Diagram of numerical simulation of combined wave and surge overtopping 數(shù)值水槽的劃分可利用Workbench中的Meshing模塊將模型劃分為四邊形結(jié)構(gòu)性網(wǎng)格,文中主要研究波浪溢流與海堤的作用,需要重點觀察海堤表面的水動力特征,則要在海堤模型上加入邊界層網(wǎng)格,劃分后的全局網(wǎng)格尺寸為0.025 m,邊界層網(wǎng)格尺寸為0.010 m。數(shù)值水槽的上方設(shè)為壓力出口邊界;底部設(shè)為壁面邊界;左端為UDF函數(shù)輸入的速度入口邊界;右端有自由出流邊界;在海堤前坡位置設(shè)置了阻尼消波區(qū)?;谏鲜瞿P偷挠嬎憔W(wǎng)格和邊界條件,將采用有限體積法(FVM)對N-S控制方程和湍流模型進行離散化分析計算,其中選定的時間步長為0.005 s。對于方程中的擴散項采用中心差分格式;對于方程中的對流項,采用二階迎風差分格式;對于壓力和速度的耦合計算,采用壓力隱式算子分割(PISO)法進行迭代求解。 文中依據(jù)Stokes波浪理論和速度方程對二階Stokes波進行仿真模擬。在Stokes波浪二階理論方程中,相對波高H/h(波高與靜水深的比值)的使用范圍在0.05~0.50之間。因此為了提高模擬的準確性,輸入的Stokes波采用了波高H=0.1 m、周期T=1.5 s、超高Rc=0.06 m(即水深為h=0.71 m)的算例,模擬得到的波浪溢流模擬效果如圖2(a)所示。為了得到模擬波浪的波形,具體操作如下:在x=10 m處建立體積分數(shù)為0.5的監(jiān)測點,隨著波浪的運動,監(jiān)測其y軸坐標隨時間的變化,由其輸出的y坐標值便可得到數(shù)值波浪水槽造波波形圖。計算得到的Stokes波形如圖2(b)所示。 圖2 Stokes波波形模擬與理論值對比Fig.2 Comparison of experimental and theoretical values of Stokes wave waveform 生成的Stokes波在比較了8個周期數(shù)值波浪模擬產(chǎn)生的波形后,發(fā)現(xiàn)基本與理論波形吻合,都呈現(xiàn)了“上尖下平”的特點,波浪要素試驗值和理論值在各個相位點之間相對誤差都小于4%。由以上結(jié)果表明該數(shù)值波浪水槽對二階Stokes波的模擬十分準確。 一般檢修舉措往往會耗費大量的人力、物力、財力,且故障排查不及時,不能消除潛在的安全隱患,往往會造成重大安全事故。而安全有效的帶電檢測不僅可以大大降低設(shè)備損壞的可能性,延長開關(guān)柜的使用壽命,而且還能使整個供電系統(tǒng)擁有長期穩(wěn)定的保障。 波浪溢流可視為越浪和溢流過程的非線性疊加,由此可分析在不同的波浪要素前提下,堤前水位超高對波浪溢流量的影響。具體操作為:在海堤模型上建立10個垂直于水平面的測量面用于測量平均波浪溢流量,其中沿著海堤堤頂?shù)染嚯x建立5個,再沿后坡每隔10 cm建立1個,共建5個,分別命名為Q1~Q10。具體觀測面布置如圖3所示。 圖3 數(shù)值模擬流量測量面分布Fig.3 Numerical simulation of flow measurement surface distribution 為了驗證數(shù)值模擬結(jié)果的正確性,選取波浪溢流的基本參數(shù)平均波浪溢流量對模型進行驗證。分別選取了不同超高Rc和不同波高H的波浪溢流現(xiàn)象進行模擬,其結(jié)果與Hughes等[8]和Li等[9]同尺寸物理試驗?zāi)P偷钠骄ɡ艘缌髁窟M行比較。結(jié)果如圖4所示,圖中橫坐標為無量綱的相對超高,縱坐標為無量綱波浪溢流量。 圖4 數(shù)值模擬和物理模型實測越堤流量結(jié)果比較Fig.4 Comparison of numerical simulation and physical model measured discharge 由圖4可知,波浪溢流量模擬結(jié)果與Hughes等[8]和Li等[9]的物理試驗結(jié)果在總體趨勢上十分接近,無量綱波浪溢流量隨著相對超高的增大不斷減小。但在超高相同的條件下數(shù)值模擬的平均波浪溢流量都相較于物理試驗結(jié)果偏小,造成這個現(xiàn)象的主要原因是物理試驗中流量測量方法和模型試驗條件的不同,從而導(dǎo)致了波浪溢流量的差異。在Hughes等[8]試驗?zāi)P椭胁ɡ艘缌髁恐饕ㄟ^ADV流速儀和波高儀對海堤后坡的流速與水流厚度進行測量,再通過流量和流速、水流厚度的數(shù)值關(guān)系計算波浪溢流量。但當小超高的條件下波浪溢流摻混著空氣導(dǎo)致水流厚度變大,從而導(dǎo)致波浪溢流量的測量值相較于真實值偏大。而Li等[9]在測量波浪溢流量時雖然改進了上述方法,采用了稱重法對波浪溢流量進行測量,但其海堤模型使用普通混凝土進行制作,表面糙率較小,所以波浪溢流量也與數(shù)值波浪模擬值相比偏大。但是總體上數(shù)值模擬結(jié)果和二者試驗數(shù)據(jù)相差不大,通過數(shù)值模擬結(jié)果與Hughes波浪溢流計算公式進行比較,發(fā)現(xiàn)數(shù)值模擬結(jié)果數(shù)據(jù)的決定系數(shù)R2為0.988,均方根誤差RMSE為0.044,由此可以認為波浪溢流數(shù)值水槽模型可以準確地模擬出波浪溢流現(xiàn)象。 對圖4進行觀察可以發(fā)現(xiàn)無量綱波浪溢流量隨相對超高增大而減小,呈現(xiàn)單調(diào)性關(guān)系,為了探究無量綱波浪溢流量與相對超高之間的關(guān)系,通過Hughes等[8]波浪溢流物理試驗數(shù)據(jù)和Pan等[17]后坡護坡為高性能加筋草的波浪溢流物理試驗數(shù)據(jù)建立了如圖5(a)所示的指數(shù)關(guān)系曲線。其中由Hughes和Pan二者計算波浪溢流量曲線的表達式可以發(fā)現(xiàn)在無量綱波浪溢流量與相對超高在Rc/H≤-0.3的范圍內(nèi)試驗數(shù)據(jù)與擬合曲線存在較好的相關(guān)性,而在-0.3 為了進一步研究穩(wěn)定溢流量與波浪溢流量之間的一般規(guī)律,首先使用Kindsvater純溢流量計算公式計算與波浪溢流數(shù)值模擬試驗組次同等海堤超高Rc下的穩(wěn)定溢流量qs的數(shù)值(計算時摩阻經(jīng)驗參數(shù)Cf取為0.544 5,同數(shù)值模型的摩阻經(jīng)驗參數(shù)相同),然后以數(shù)值模擬的波浪溢流量qws和同等海堤超高Rc下對應(yīng)的穩(wěn)定溢流量qs的比值qws/qs為縱坐標,以海堤相對超高Rc/H為橫坐標繪制圖5(b)。 圖5 波浪溢流量中波浪與溢流主導(dǎo)性關(guān)系Fig.5 Dominant relationship between wave and overflow in combined wave and surge overtopping 由圖5(b)可知,當海堤相對超高Rc/H絕對值較大時,比值qws/qs趨近于1,即波浪溢流量接近對應(yīng)的穩(wěn)定溢流量,此時波浪溢流過程中溢流占據(jù)主導(dǎo)作用,波浪溢流與純溢流相似,相對超高與波浪溢流量之間主要呈現(xiàn)為指數(shù)關(guān)系;當海堤相對超高Rc/H絕對值趨近于0時,比值qws/qs迅速增大,波浪溢流量達到對應(yīng)穩(wěn)定溢流量的數(shù)倍,此時波浪溢流過程中越浪逐漸占據(jù)主導(dǎo)作用,由于越浪量主要滿足韋伯分布,打破了波浪溢流量與相對超高之間較為穩(wěn)定的指數(shù)關(guān)系,從而使得波浪溢流量隨著相對超高的增加逐漸表現(xiàn)為隨機性,使用Hughes和Pan的指數(shù)函數(shù)公式計算平均波浪溢流量將不再準確。這種誤差以Rc/H=-0.3為明顯的分界點,這解釋了圖5(a)所示的波浪溢流量在Rc/H=-0.3兩側(cè)表現(xiàn)出不同的分布趨勢的原因。因此,可以以Rc/H=-0.3為界,將Rc/H≤-0.3的波浪溢流稱為溢流主導(dǎo)的波浪溢流,將-0.3 在對波浪溢流主導(dǎo)性作用研究中發(fā)現(xiàn),波浪溢流越過堤頂后,當?shù)啄ψ?、紊動等產(chǎn)生的摩阻比降與坡度相同時,后坡水深從脈動運動逐漸達到相對穩(wěn)定的狀態(tài)。為了進一步研究波浪溢流水力學(xué)參數(shù)后坡穩(wěn)定水深變化情況,可以利用FLUENT中監(jiān)測功能沿海堤模型表面定義10個流速測量點(D1~D10),各個測點與流量測量面Q1~Q10處于同一水平位置。通過流量—流速—水流厚度之間的關(guān)系近似得到該測點的水流厚度。選擇3種不同超高的模擬工況對后坡水流厚度進行研究,具體模擬組次安排如表1所示。 表1 模擬試驗組次Tab.1 Test groups 由水力學(xué)基本理論可知波浪溢流發(fā)生后,波浪溢流由堤頂流向后坡的過程中當摩阻比降和坡度相等時,水流流速和厚度不再變化并達到后坡平均穩(wěn)定水流厚度dsm。為了準確了解后坡水流狀態(tài)隨時間的變化情況,選取上述Rc=-0.06 m的模擬組次WS7,對后坡D4、D6、D9、D10處水流厚度結(jié)果進行分析,其中D4位于堤頂,D6測點位于后坡上段,D9和D10測點位于后坡下段(具體位置分布見圖3)。其模擬結(jié)果如圖6所示。 如圖6所示,當水流在D4、D6處時,水流厚度呈現(xiàn)大波動形態(tài),尤其D6處于堤頂與后坡的交界位置,受結(jié)構(gòu)突變的影響,水流厚度變化十分劇烈。隨著時間推移波浪溢流非恒定水流峰值也在不斷的向后移動,產(chǎn)生如此現(xiàn)象的主要原因是受到了波浪運動的影響導(dǎo)致堤頂和后坡上端水流產(chǎn)生巨大紊動,峰值在后堤向下周期性運動,水流狀態(tài)十分不穩(wěn)定;隨著水流運動后移到達后坡下段D9、D10處時,相同時間節(jié)點下(t=42 s)二者的平均水流厚度相同,此時水流受后坡阻力的影響,水位漸漸處于穩(wěn)定,呈現(xiàn)為多峰小波動狀態(tài)。其中后坡流速不斷提高,導(dǎo)致后坡下端水流厚度小于堤前水位,最終平均水流厚度逐漸趨向穩(wěn)定,達到后坡平均穩(wěn)定水流厚度dsm。 圖6 試驗4個測點處Rc=-0.06 m時水流厚度d隨時間變化Fig.6 Variation of flow thickness d with time at Rc=-0.06 m at four test points 根據(jù)表1提供的不同模擬工況,測量計算各個測點(D1~D10)水流厚度平均值dm沿著堤頂后坡的變化情況,并與相同超高條件下的純溢流平均水流厚度進行比較,比較結(jié)果如圖7所示。 由圖7可知,當超高Rc較大時,波浪溢流作用下海堤堤頂和后坡的平均水流厚度dm相較于純溢流偏大,隨著超高的減小,平均水流厚度dm將趨近純溢流作用下沿海堤堤頂和后坡的水流厚度。這主要是因為當超高較大時,波浪溢流越堤流量受波浪影響較大,波浪對波浪溢流量的影響占據(jù)主導(dǎo)地位,產(chǎn)生了漲水效果并有一定的越浪量,因此平均水流厚度較大。隨著超高Rc的減小,波浪影響減弱,平均波浪溢流越堤流量逐漸接近純溢流量。當相同超高條件下,不同組次的平均水流厚度沿著后坡逐漸減小并接近,這主要是因為波浪作用范圍主要在堤頂和后坡上端,而后坡下端受波浪影響較小,從而水流狀態(tài)逐漸穩(wěn)定類似于純溢流。越堤水流在由堤頂流向后坡的過程中水流厚度逐漸減小,水流流速逐漸增大,當因底摩阻、紊動等產(chǎn)生的摩阻比降和坡度相等時,平均水流厚度和水流流速不再發(fā)生變化,此時的水流厚度達到后坡平均穩(wěn)定水流厚度dsm。 圖7 平均水流水深dm沿堤頂后坡的變化Fig.7 Variation of average flow velocity along the dike 為了探究平均穩(wěn)定水流厚度dsm估算方法,經(jīng)過與不同特征參數(shù)進行組合分析,建立無量綱化后坡平均穩(wěn)定水流厚度dsm/H與無量綱化平均波浪溢流越堤流量qws/(gH3)1/2之間的相關(guān)關(guān)系,結(jié)果如圖8所示。 圖8 無量綱化后坡平均穩(wěn)定水流厚度dsm/H與平均波浪溢流越堤流量qws/(gH3)1/2關(guān)系Fig.8 Relationship between average steady flow thickness and average combined wave and surge overtopping discharge 圖8中的擬合曲線可表示為: (14) 式中:dsm為波浪溢流平均穩(wěn)定水流厚度,qws為平均波浪溢流越堤流量。上述公式的適用范圍為-0.9 在運用FLUENT軟件的基礎(chǔ)上,建立了自定義的UDF速度邊界函數(shù)進行造波,以VOF方法為基礎(chǔ)加入了RANS控制方程進行離散求解。驗證了數(shù)值波浪溢流水槽造波和消波的正確性,完整的模擬了波浪溢流的生成、傳播以及與海堤相互作用的過程。根據(jù)此模型研究了波浪溢流發(fā)生后波浪溢流主導(dǎo)性問題和海堤后坡水流厚度等內(nèi)容,得出了如下結(jié)論: 1)對波浪溢流的主導(dǎo)性進行定義,以Rc/H=-0.3為界,將Rc/H≤-0.3的波浪溢流稱為溢流主導(dǎo)的波浪溢流,將-0.3 2)由于受到了波浪運動的影響,導(dǎo)致堤頂和后坡上端水流產(chǎn)生巨大紊動,水體厚度峰值在后堤向下周期性運動,導(dǎo)致波浪溢流非恒定水流峰值隨著時間不斷后移,當?shù)竭_后坡下段D9、D10處時,相同時間節(jié)點下平均水流厚度趨近一致,最終達到后坡平均穩(wěn)定水流厚度dsm。 3)以結(jié)論2)為基礎(chǔ),對不同特征參數(shù)進行組合分析,建立無量綱化后坡平均穩(wěn)定水流厚度dsm/H與無量綱化平均波浪溢流越堤流量qws/(gH3)1/2之間的函數(shù)關(guān)系,得到了波浪溢流后坡穩(wěn)定水流厚度的計算公式。1.3 波浪水槽的造波和消波
2 波浪溢流水槽的建立與驗證
2.1 數(shù)值波浪溢流模型和參數(shù)設(shè)置
2.2 數(shù)值波浪水槽的驗證
3 模擬結(jié)果與分析
3.1 波浪溢流量分析
3.2 波浪溢流后坡水流厚度分析
4 結(jié) 語