鄒國(guó)良，張慶河

（1. 天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072；2. 南京水利科學(xué)研究院，南京 210024）

?

基于非靜壓方程的波浪與孔隙結(jié)構(gòu)物相互作用模擬

鄒國(guó)良1, 2，張慶河1

（1. 天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072；2. 南京水利科學(xué)研究院，南京 210024）

摘 要：通過將Darcy-Forchheimer運(yùn)動(dòng)方程引入非靜壓方程中，提出了描述孔隙流體運(yùn)動(dòng)的非靜壓波浪模型.將所建立的數(shù)學(xué)模型應(yīng)用到波浪與孔隙結(jié)構(gòu)物相互作用的數(shù)值模擬中，并采用物理實(shí)驗(yàn)的實(shí)驗(yàn)值與數(shù)值解進(jìn)行了對(duì)比，結(jié)果表明，所建立的模型在垂向只需要較少層數(shù)（2～3層）即可合理地描述孔隙結(jié)構(gòu)物前后的波浪變形，為描述波浪與結(jié)構(gòu)物相互作用提供了一種新的解決方案.

關(guān)鍵詞：非靜壓方程；孔隙結(jié)構(gòu)物；部分反射；拋石防波堤

近岸區(qū)由于建筑物的存在，尤其是港口工程中如防波堤等建筑物，使得波浪的模擬變得復(fù)雜.如在模擬港內(nèi)波浪的泊穩(wěn)條件時(shí)，需要考慮防波堤對(duì)波浪的掩護(hù)作用，尤其是在防波堤附近因波浪的反射、繞射以及透射等綜合作用對(duì)港內(nèi)碼頭泊位前泊穩(wěn)條件的影響.因此波浪數(shù)學(xué)模型不可避免地需要解決波浪與建筑物相互作用的模擬.孔隙結(jié)構(gòu)是海岸防護(hù)工程的主要結(jié)構(gòu)形式，這類結(jié)構(gòu)除了受到外部波浪、水流等直接作用外，水體還會(huì)透過孔隙結(jié)構(gòu)直接進(jìn)入其內(nèi)部，形成較為復(fù)雜的內(nèi)部孔隙流動(dòng).因此，了解孔隙介質(zhì)中的流體運(yùn)動(dòng)對(duì)海岸工程的設(shè)計(jì)、預(yù)估等具有十分重要的意義.

早期關(guān)于波浪與孔隙結(jié)構(gòu)相互作用的研究主要采用勢(shì)流理論對(duì)波浪遇到孔隙結(jié)構(gòu)時(shí)的反射和透射系數(shù)進(jìn)行解析解推導(dǎo)[1].采用數(shù)學(xué)模型對(duì)波浪與孔隙結(jié)構(gòu)相互作用進(jìn)行數(shù)值研究的工作主要是結(jié)合Darcy-Forchheimer運(yùn)動(dòng)方程進(jìn)行的.其中，平面二維數(shù)學(xué)模型應(yīng)用較多的為沿水深積分的緩坡方程[2]和Boussinesq方程[3]模型，這類模型在工程中應(yīng)用較為廣泛.但由于方程沿水深進(jìn)行了積分，不能刻畫潛式結(jié)構(gòu)物.三維（垂向二維）數(shù)學(xué)模型主要是借助于VOF法[4-7]、SPH方法[8]求解雷諾時(shí)均N-S方程（RANS方程）來模擬波浪與孔隙介質(zhì)相互作用，這類模型可以描述復(fù)雜的自由表面運(yùn)動(dòng)，但對(duì)計(jì)算的穩(wěn)定性及計(jì)算量的要求仍較高.

基于非靜壓方程的波浪數(shù)學(xué)模型與RANS模型類似都需要求解壓力泊松方程獲得壓力項(xiàng).所不同的是，非靜壓模型通過求解水位方程來追蹤自由面，因此在自由面附近的垂向網(wǎng)格劃分不需要像RANS模型（VOF、MAC等方法求解自由面）那么精細(xì).在描述非線性、色散性較強(qiáng)的非破碎波運(yùn)動(dòng)時(shí)，目前較先進(jìn)的非靜壓波浪模型通過將壓力項(xiàng)定義在網(wǎng)格單元邊中心，并結(jié)合Keller-box離散格式[9-10]，故可精確描述自由面壓力邊界條件，因此在模擬強(qiáng)非線性、色散性波浪傳播變形時(shí)，垂向一般只需要2～3 層[11-12]即可獲得較準(zhǔn)確的模擬結(jié)果，大大提高了模型的計(jì)算效率與適用性，成為近年來研究較多的新型波浪數(shù)學(xué)模型.這類模型基于非線性淺水方程，在動(dòng)量方程中考慮了非靜水壓力梯度項(xiàng)，數(shù)值離散格式相對(duì)于高階Boussinesq方程簡(jiǎn)單，也不受相對(duì)水深、無黏性流體、無旋運(yùn)動(dòng)等條件限制，適用于波浪從深水到淺水傳播變形、波流相互作用等問題的研究.而采用非靜壓波浪模型研究波浪與結(jié)構(gòu)物相互作用的研究仍然較少.為此，本文結(jié)合Darcy-Forchheimer運(yùn)動(dòng)方程，基于垂向二維非靜壓型實(shí)現(xiàn)波浪與孔隙結(jié)構(gòu)物相互作用的模擬，為近岸存在結(jié)構(gòu)物時(shí)波浪的準(zhǔn)確模擬奠定基礎(chǔ).

1　控制方程

假定孔隙介質(zhì)為各向同性、均勻的剛性介質(zhì)，且孔隙介質(zhì)底部邊界為剛性不可滲透邊界.引入Darcy-Forchheimer方程[13]，孔隙介質(zhì)內(nèi)流體運(yùn)動(dòng)的控制方程為

式中：ξ為自由面；g為重力加速度；t為時(shí)間；u、w分別為x和z方向的速度分量；υ為渦黏系數(shù)；q為非靜壓項(xiàng)；對(duì)于波浪的運(yùn)動(dòng)密度可認(rèn)為恒定，設(shè)密度ρ＝1,kg/m3；cr為附加慣性系數(shù)；ɑ和b分別為線性和非線性阻力系數(shù).

通過將Darcy-Forchheimer運(yùn)動(dòng)方程與非靜壓方程相結(jié)合建立的孔隙流體運(yùn)動(dòng)的控制方程，不需要明確地區(qū)分孔隙介質(zhì)與流體的邊界[14].因此，可以同時(shí)求解孔隙介質(zhì)內(nèi)外流體的運(yùn)動(dòng)方程.孔隙介質(zhì)外的流體孔隙率n＝1，附加慣性系數(shù)cr＝1，線性和非線性阻力系數(shù)不予以考慮，即ɑ＝b＝0.孔隙介質(zhì)內(nèi)的流體運(yùn)動(dòng)則通過考慮附加慣性力項(xiàng)、線性和非線性阻力項(xiàng)來求解，此時(shí)方程中的u和w分別表示為x和z方向的孔隙流體速度分量.

將連續(xù)性方程（3）沿水深積分并采用自由面和底部運(yùn)動(dòng)學(xué)邊界條件通過萊布尼茲積分轉(zhuǎn)換可得到非靜壓模型的水位方程

式中d為靜水深.所采用的自由面和底部運(yùn)動(dòng)學(xué)邊界條件分別為

垂向采用σ坐標(biāo)變換，即

則垂向坐標(biāo)變換后的控制方程為

式中ω為σ坐標(biāo)下的垂向流速.ω與垂向流速w之間的關(guān)系為

采用Darcy-Forchheimer運(yùn)動(dòng)方程研究波浪與孔隙結(jié)構(gòu)相互作用時(shí)，其關(guān)鍵在于如何確定線性和非線性阻力系數(shù)ɑ和b.Engelund[15]基于系列物理模型實(shí)驗(yàn)研究針對(duì)海岸工程中的孔隙結(jié)構(gòu)提出了應(yīng)用較為廣泛的線性和非線性阻力系數(shù)公式

式中：n為孔隙率；d50為特征孔隙粒徑；ν為水的運(yùn)動(dòng)黏滯系數(shù)；α和β為Forchheimer常數(shù).Ergun[16]建議α＝150，β＝1.75；Engelund[15]和Madsen[17]建議α＝780～1,500，β＝1.8～3.6；van Gent[18]通過物理實(shí)驗(yàn)建議α＝1,000，β＝1.1.

附加慣性系數(shù)cr[2]可表示為

式中cm為附加質(zhì)量系數(shù).Van Gent[18]通過物理試驗(yàn)建議取cm＝0.34，Sulisz[19]根據(jù)物理實(shí)驗(yàn)建議取cm＝2；Huang等[6]、任冰等[8]在其RANS模型中將附加質(zhì)量力系數(shù)取0，取得了較好的模擬結(jié)果.

對(duì)于潛式孔隙結(jié)構(gòu)，當(dāng)垂向網(wǎng)格劃分使得某計(jì)算層同時(shí)含有水體和孔隙結(jié)構(gòu)時(shí)，該層整體均被考慮為孔隙結(jié)構(gòu)，但孔隙率需要根據(jù)實(shí)際孔隙材料所占的體積進(jìn)行調(diào)整.如圖1所示，計(jì)算域垂向分為3層，第2層同時(shí)含有水體和結(jié)構(gòu)物，此時(shí)該層的孔隙率則根據(jù)該層結(jié)構(gòu)物的高度與結(jié)構(gòu)孔隙率進(jìn)行插值，即第2層的孔隙率n2為

式中h1和h2分別為第2層網(wǎng)格中水體的高度和結(jié)構(gòu)物的高度.

圖1　潛式孔隙結(jié)構(gòu)的垂向分層Fig.1 Vertical layers for submerged porous structure

2　邊界條件

自由面處垂向流速滿足運(yùn)動(dòng)學(xué)邊界條件式（5），且自由面壓力邊界滿足

底部邊界上的垂向流速滿足底部運(yùn)動(dòng)學(xué)邊界條件式（6），且底部水平流速梯度滿足

固壁邊界采用全反射邊界，即垂直于邊界的速度為零，對(duì)于垂向二維水槽則有

出流邊界采用輻射邊界條件并結(jié)合海綿層消波技術(shù).輻射邊界條件為

波浪入射邊界采用鄒國(guó)良等[20]基于非靜壓模型提出的域內(nèi)質(zhì)量線源造波方法以避免波浪二次反射.數(shù)值消波邊界采用Zou等[21]基于非靜壓模型改進(jìn)的海綿層消波方式進(jìn)行數(shù)值消波.

3　數(shù)值離散和計(jì)算方法

數(shù)值離散時(shí)采用有限差分格式和有限體積法相結(jié)合的方式進(jìn)行空間離散.變量布置采用交錯(cuò)網(wǎng)格布置，其中非靜壓q定義在單元邊中心，如圖2所示.動(dòng)量方程中的對(duì)流項(xiàng)采用Stelling等[22]針對(duì)交錯(cuò)網(wǎng)格提出的動(dòng)量守恒的對(duì)流項(xiàng)格式離散.垂向動(dòng)量方程采用Keller-box格式離散，該格式垂向只需要較少的網(wǎng)格層數(shù)就可以較準(zhǔn)確地描述短波的色散特性[12].對(duì)于時(shí)間項(xiàng)的積分，動(dòng)量方程中的對(duì)流項(xiàng)采用顯格式，靜壓梯度項(xiàng)和非靜壓梯度項(xiàng)采用θ隱格式.

圖2　交錯(cuò)網(wǎng)格系統(tǒng)變量布置Fig.2 Arrangement of the variables in a staggered grid system

沿各層積分后孔隙流體的水平方向動(dòng)量方程為

各層垂向動(dòng)量方程為

各層的連續(xù)性方程為

式中：σk為第k層的σ坐標(biāo)；hk為第k層的高度；uk為第k層的水平流速；wk為第k層z坐標(biāo)下的垂向流速；kω為第k層σ坐標(biāo)下的垂向流速；qk為第k層的非靜壓項(xiàng).

模型在求解上述壓力速度耦合方程時(shí)采用分步法并結(jié)合壓力校正法求解[12].具體離散格式見文獻(xiàn)[12]和文獻(xiàn)[20]，數(shù)值計(jì)算流程見圖3.

圖3　數(shù)值計(jì)算流程Fig.3 Flow chart of numerical calculation

4　計(jì)算結(jié)果與分析

為了對(duì)孔隙流體的非靜壓方程進(jìn)行驗(yàn)證，分別采用行進(jìn)波遇孔隙結(jié)構(gòu)物的部分反射實(shí)驗(yàn)和波浪與拋石防波堤相互作用的實(shí)驗(yàn)來驗(yàn)證.

4.1波浪遇孔隙結(jié)構(gòu)的部分反射

Hirayama[23]通過在Boussinesq方程模型的動(dòng)量方程中引入層流線性阻力項(xiàng)和紊流非線性阻力項(xiàng)，實(shí)現(xiàn)了波浪遇建筑物的部分反射，并將數(shù)值模型計(jì)算結(jié)果與其實(shí)驗(yàn)進(jìn)行了對(duì)比.實(shí)驗(yàn)水槽長(zhǎng)35,m，水深為0.279,m，孔隙結(jié)構(gòu)設(shè)置在水槽右端一高0.004,m、坡度為1∶30的斜坡上，如圖4所示.孔隙結(jié)構(gòu)由不規(guī)則的四角塊體組成，實(shí)驗(yàn)測(cè)定結(jié)構(gòu)的孔隙率為0.45，平均粒徑為0.045,8,m，孔隙結(jié)構(gòu)的寬度為0.5,m.本文在計(jì)算時(shí)采用與物理實(shí)驗(yàn)材料一致的孔隙率以及孔隙粒徑，F(xiàn)orchheimer常數(shù)α和β按照Hirayama[23]建議值分別取為2,100和2.2，附加質(zhì)量系數(shù)為0. 波浪部分反射實(shí)驗(yàn)條件如表1所示.對(duì)于規(guī)則波，表中波高和周期分別為規(guī)則波波高H和周期T；對(duì)于不規(guī)則波表中波高和周期分別指有效波高H1/3和周期T1/3，頻譜采用B-M譜[23].

圖4　Hirayama實(shí)驗(yàn)示意Fig.4 Experiment by Hirayama

表1　波浪部分反射實(shí)驗(yàn)條件Tab.1 Experimental conditions of wave partial reflections

數(shù)值計(jì)算時(shí)水槽尺寸同實(shí)驗(yàn)水槽，水槽左端設(shè)置海綿層消波，水平向網(wǎng)格設(shè)置為Δx＝L/50，垂向分為2層，時(shí)間步長(zhǎng)為T/100.為了獲得穩(wěn)定的反射波，規(guī)則波計(jì)算總時(shí)間為120,s，不規(guī)則波計(jì)算總時(shí)間為600,s.這里的數(shù)值模擬采用相同的反射系數(shù)計(jì)算方法.圖5和圖6分別為規(guī)則波和不規(guī)則波作用下，采用Boussinesq方程模型計(jì)算的反射系數(shù)[23]和非靜壓模型計(jì)算出的反射系數(shù)與實(shí)驗(yàn)給出的反射系數(shù)對(duì)比結(jié)果.根據(jù)反射系數(shù)的計(jì)算結(jié)果可看出，非靜壓模型計(jì)算出的反射系數(shù)隨波高以及周期變化的趨勢(shì)與實(shí)驗(yàn)結(jié)果趨勢(shì)接近.從反射系數(shù)的大小來看，采用Boussinesq方程模型計(jì)算出的反射系數(shù)整體比實(shí)驗(yàn)結(jié)果偏大，而采用非靜壓模型計(jì)算出的反射系數(shù)整體上接近實(shí)驗(yàn)結(jié)果.

圖7為相同計(jì)算參數(shù)條件下水槽右側(cè)分別設(shè)置海綿層消波以及孔隙結(jié)構(gòu)部分反射時(shí)，沿水槽波高計(jì)算結(jié)果.根據(jù)計(jì)算出的波高結(jié)果可明顯看出，水槽右側(cè)邊界在采用海綿層消波時(shí)，消波段內(nèi)波高逐漸衰減，左側(cè)計(jì)算出的波高穩(wěn)定，波浪無反射；而水槽右側(cè)設(shè)置孔隙結(jié)構(gòu)時(shí)，波浪在孔隙結(jié)構(gòu)物內(nèi)產(chǎn)生反射波，與入射波疊加形成立波.由此也可反映出模型較好地描述了波浪消波邊界和部分反射邊界.

圖5　規(guī)則波作用下反射系數(shù)計(jì)算值與實(shí)測(cè)值對(duì)比Fig.5 Comparison of measured and computed reflection coefficients under regular waves

圖6　不規(guī)則波作用下反射系數(shù)計(jì)算值與實(shí)測(cè)值對(duì)比Fig.6 Comparison of measured and computed reflection coefficients under irregular waves

圖7　消波邊界與部分反射邊界計(jì)算波高Fig.7 Computed wave heights on absorbing and partial reflective boundaries

4.2波浪與拋石防波堤的相互作用

Garcia等[24]以及Lara等[25]通過在RANS方程中引入Forchheimer運(yùn)動(dòng)方程中的線性和非線性阻力項(xiàng)，結(jié)合VOF法求解自由面進(jìn)行了規(guī)則波和不規(guī)則波與拋石防波堤的相互作用的研究，并與其水槽中的波浪與拋石防波堤相互作用的物理模型實(shí)驗(yàn)進(jìn)行了對(duì)比.上述兩篇文獻(xiàn)分別給出了規(guī)則波和不規(guī)則波作用下拋石防波堤前后的物理實(shí)驗(yàn)波高等相關(guān)實(shí)測(cè)數(shù)據(jù).實(shí)驗(yàn)水槽長(zhǎng)24.05,m、寬0.6,m、深0.8,m.拋石防波堤堤頂寬1,m，堤底寬2,m，堤身高0.25,m，堤兩側(cè)坡度為1∶2.防波堤由護(hù)面塊石以及堤心石兩層組成，其中護(hù)面塊石的平均粒徑為3.94,cm，孔隙率為0.53，堤心石平均粒徑為1.18,cm，孔隙率為0.49.防波堤設(shè)置在水槽右端高度為0.1,m、坡度為1∶20的斜坡上，堤后同樣設(shè)有1∶20斜坡進(jìn)行消波.波浪與拋石防波堤相互作用水槽實(shí)驗(yàn)如圖8所示.

圖8　波浪與拋石防波堤相互作用水槽實(shí)驗(yàn)（單位：m）Fig.8 Experiment on wave interaction with rubble mound breakwater（unit：m）

數(shù)值波浪水槽的長(zhǎng)度為35,m，水槽兩端各設(shè)置3,m寬的海綿層，網(wǎng)格尺寸為Δx＝0.02,m，垂向分為3層.Garcia等[24]在數(shù)值計(jì)算時(shí)將防波堤按照原型中護(hù)面塊石和堤心石分別給定不同的參數(shù)，并分別給出了Forchheimer常數(shù)，護(hù)面層α和β 分別取1,000和0.8，防波堤內(nèi)部則分別取為1,000和1.2.本文為便于計(jì)算將防波堤設(shè)置為統(tǒng)一材料屬性，α 和β分別取1,000和1.0，孔隙率和平均粒徑均取兩種材料的平均值分別為0.51和2.56,cm，附加質(zhì)量系數(shù)為0.34.

規(guī)則波計(jì)算時(shí)間為30,s，取最后10,s進(jìn)行波高計(jì)算，不規(guī)則波總計(jì)算時(shí)間為150,s，測(cè)點(diǎn)采樣頻率為10,Hz.圖9分別為波高計(jì)算值與實(shí)測(cè)值對(duì)比，其中不規(guī)則波波高為均方根波高.圖10為不規(guī)則波作用下防波堤前后主要測(cè)點(diǎn)的波能譜計(jì)算值與實(shí)驗(yàn)值對(duì)比.根據(jù)對(duì)比結(jié)果可知，規(guī)則波作用下，防波堤前A1～A3測(cè)點(diǎn)以及堤后的A7～A11測(cè)點(diǎn)的計(jì)算波高基本與實(shí)測(cè)波高一致，而防波堤前A5測(cè)點(diǎn)的計(jì)算波高略微偏小.不規(guī)則波作用下，堤身A6測(cè)點(diǎn)以及堤前A5測(cè)點(diǎn)計(jì)算波高略微偏小，而堤前、堤后的計(jì)算波高與實(shí)測(cè)波高基本一致.堤前計(jì)算波高偏小原因一方面是非靜壓模型所取的Forchheimer常數(shù)是根據(jù)兩種不同孔隙結(jié)構(gòu)進(jìn)行的平均，不能反映兩種不同結(jié)構(gòu)的真實(shí)滲透屬性，因此參數(shù)需要進(jìn)一步的率定；另一方面是非靜壓模型計(jì)算的垂向網(wǎng)格數(shù)并沒有如RANS模型所需要的密，從而對(duì)防波堤形狀的刻畫不如后者.但計(jì)算波高總體與實(shí)驗(yàn)實(shí)測(cè)波高接近.

圖9　波高計(jì)算值與實(shí)測(cè)值對(duì)比Fig.9 Comparison of measured and computed wave heights

圖10　波能譜計(jì)算值與實(shí)測(cè)值對(duì)比Fig.10 Comparison of measured and computed wave energy spectra

表2為波浪與拋石防波堤相互作用的實(shí)驗(yàn)條件，表中不規(guī)則波采用JONSWAP譜，譜峰頻率為3.3,Hz，波高為有效波高，周期為譜峰周期.

表2　波浪與拋石防波堤相互作用實(shí)驗(yàn)條件Tab.2 Experimental wave conditions of wave interaction with rubble mound breakwater

4.3孤立波與孔隙結(jié)構(gòu)物相互作用

為了進(jìn)一步驗(yàn)證模型在模擬強(qiáng)非線性波與孔隙結(jié)構(gòu)物相互作用時(shí)的精度，采用Lynett等[3]曾在水槽中進(jìn)行的孤立波與直立式孔結(jié)構(gòu)物相互作用的實(shí)驗(yàn)進(jìn)行驗(yàn)證.波浪水槽靜水深為0.1,m，孤立波波高為0.025,m（H/d＝0.25），孔隙介質(zhì)由平均直徑厚度為0.02,m的沙礫構(gòu)成，外部由鐵絲網(wǎng)固定，有效孔隙率為0.5.孔隙結(jié)構(gòu)寬為0.3,m，高為0.14,m（出水結(jié)構(gòu)），在結(jié)構(gòu)物前、后各1,m處分別布置測(cè)點(diǎn)PWG1和PWG2，如圖11所示.

圖11　孤立波與孔隙結(jié)構(gòu)相互作用示意（單位：cm）Fig.11 Sketch of solitary wave interaction with porous structure（unit：cm）

數(shù)值波浪水槽長(zhǎng)30,m，網(wǎng)格尺寸為Δx＝0.02,m，垂向分為3層，孔隙結(jié)構(gòu)物位于x＝10.0～10.3,m. Lynett等[3]利用Boussinesq方程采用3組不同的層流、紊流阻力系數(shù)對(duì)實(shí)驗(yàn)進(jìn)行數(shù)值驗(yàn)證，其結(jié)果表明，當(dāng)α＝1,100、β＝0.81時(shí)Boussinesq方程數(shù)值結(jié)果與實(shí)驗(yàn)結(jié)果較吻合.本文在數(shù)值計(jì)算時(shí)，同樣地采用3組不同阻力系數(shù)進(jìn)行數(shù)值計(jì)算并與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比，附加質(zhì)量系數(shù)為0，孔隙率和平均粒徑均采用與實(shí)驗(yàn)一致的參數(shù).

圖12給出了采用非靜壓模型，結(jié)合不同阻力系數(shù)模擬獲得的測(cè)點(diǎn)PWG1和PWG2位置處波面計(jì)算值與實(shí)測(cè)值對(duì)比，其中將孤立波波峰經(jīng)過測(cè)點(diǎn)PWG1的時(shí)刻作為“0”時(shí)刻.根據(jù)對(duì)比，在3組不同阻力系數(shù)條件下，孔隙結(jié)構(gòu)物前測(cè)點(diǎn)PWG1處入射波和反射波波形以及達(dá)到時(shí)間的計(jì)算值與實(shí)驗(yàn)值基本接近；結(jié)構(gòu)物后測(cè)點(diǎn)PWG2處透射波到達(dá)時(shí)間的計(jì)算值與實(shí)驗(yàn)值接近，但波面存在一定的差別.根據(jù)圖12（b）可知，隨著層流、紊流阻力系數(shù)的增大，波浪傳播至孔隙結(jié)構(gòu)物后透射波波面呈減小趨勢(shì)，其中α＝1800、β＝1.1時(shí)非靜壓模型計(jì)算出的孔隙結(jié)構(gòu)物后的透射波波面與實(shí)測(cè)波面基本吻合.

圖12　波面高程計(jì)算值與實(shí)測(cè)值對(duì)比Fig.12 Comparison of computed and measured wave surface elevation

5　結(jié)論

（1）將Darcy-Forchheimer運(yùn)動(dòng)方程引入到垂向二維非靜壓方程中，并結(jié)合非靜壓模型數(shù)值求解框架實(shí)現(xiàn)了波浪與孔隙結(jié)構(gòu)物相互作用的模擬.

（2）將孔隙流體的非靜壓方程應(yīng)用到行進(jìn)波遇孔隙結(jié)構(gòu)部分反射的數(shù)值模擬中，數(shù)值計(jì)算與實(shí)驗(yàn)結(jié)果對(duì)比表明，模型可較好地描述孔隙結(jié)構(gòu)物前波浪的部分反射；波浪水槽右端設(shè)置消波邊界的計(jì)算波高與部分反射計(jì)算的波高對(duì)比表明，模型的消波邊界具有較高的消波效率.

（3）孔隙流體非靜壓方程的波浪與拋石防波堤相互作用的數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比表明，模型在垂向僅需要較少的網(wǎng)格層數(shù)即可合理地描述拋石防波堤前的反射、堤后的透射.

（4）孤立波與孔隙結(jié)構(gòu)物的非靜壓計(jì)算結(jié)果表明，模型可合理地描述強(qiáng)非線性波浪在孔隙結(jié)構(gòu)物前后的反射和透射作用.

由結(jié)論（2）～（4）可知，文中所建立的描述波浪與孔隙結(jié)構(gòu)相互作用的非靜壓模型在垂向僅需要較少的分層條件下即可合理地描述孔隙結(jié)構(gòu)物前后的波浪部分反射、透射作用，具有進(jìn)一步在工程波浪數(shù)值模擬中推廣應(yīng)用的前景.

參考文獻(xiàn)：

［1］Dalrymple R A，Losada M A，Martin P A. Reflection and transmission fron porous structures under oblique wave tank[J]. Journɑl of Fluid Mechɑnics，1991，224：625-644.

［2］Hsu T W，Chang J Y，Lan Y J，et al. A parabolic equation for wave propagation over porous structures[J]. Coɑstɑl Engineering，2008，55（12）：1148-1158.

［3］Lynett P，Liu P L F，Losada I J，et al. Solitary wave interaction with porous breakwaters[J]. Journɑl of Wɑterwɑy，Port，Coɑstɑl，ɑnd Oceɑn Engineering，2000，126（6）：314-322.

［4］Liu P L F，Lin P，Chang K A，et al. Numerical modeling of wave interaction with porous structures[J]. Journɑl of Wɑterwɑy，Port，Coɑstɑl，ɑnd Oceɑn Engineering，1999，125（6）：322-330.

［5］Hsu T J，Sakakiyama T，Liu P L F. A numerical model for wave motions and turbulence flows in front of a composite breakwater[J]. Coɑstɑl Engineering，2002，46（1）：25-50.

［6］Huang C J，Chang H H，Hwung H H. Structural permeability effects on the interaction of a solitary wave and a submerged breakwater[J]. Coɑstɑl Engineering，2003，49（1/2）：1-24.

［7］Akbari H，Namin M M. Moving particle method for modeling wave interaction with porous structures[J]. Coɑstɑl Engineering，2013，74：59-73.

［8］任 冰，葉曉文，高 睿，等. 波浪與多空介質(zhì)結(jié)構(gòu)物相互作用SPH模擬[J]. 海洋工程，2012，30（2）：46-53. Ren Bing，Ye Xiaowen，Gao Rui，et al. SPH modelling of wave interaction with porous structures[J]. The Oceɑn Engineering，2012，30（2）：46-53（in Chinese）.

［9］Zijlema M，Stelling G S. Further experiences with computing non-hydrostatic free-surface flows involving water waves[J]. Internɑtionɑl Journɑl for Numericɑl Methods in Fluids，2005，48（2）：169-197.

［10］Ma G，Shi F，Kirby J T. Shock-capturing nonhydrostatic model for fully dispersive surface wave processes[J]. Oceɑn Modelling，2012，43/44：22-35.

［11］Ai C，Jin S，Lü B. A new fully non-hydrostatic 3D free surface flow model for water wave motions[J]. Internɑtionɑl Journɑl for Numericɑl Methods in Fluids，2011，66（11）：1354-1370.

［12］Zijlema M，Stelling G S，Smit P. SWASH：An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters[J]. Coɑstɑl Engineering，2011，58（10）：992-1012.

［13］Sollitt C K，Cross R H. Wave transmission through permeable breakwaters[C]// Proceedings of the 13th Internɑtionɑl Conference on Coɑstɑl Engineering. New York，USA，1972：1827-1846.

［14］Beavers G S，Joseph D D. Boundary condition at a naturally permeable wall[J]. Journɑl of Fluid Mechɑnics，1967，30：197-207.

［15］Engelund F. On the laminar and turbulent flows of ground water through homogeneous sand[J]. Trɑnsɑctions of the Dɑnish Acɑdemy of Technicɑl Sciences，1953，3（4）：7-105.

［16］Ergun S. Fluid flow through packed columns[J]. Chemicɑl Engineering Progress，1952，48（2）：89-94.

［17］Madsen O S. Wave transmission through porous structures[J]. Journɑl of Wɑterwɑy，Hɑrbors ɑnd Coɑstɑl Engineering，1974，100（3）：169-188.

［18］Van Gent M R A. Porous flow through rubble-mound material[J]. Journɑl of Wɑterwɑy，Port，Coɑstɑl，ɑnd Oceɑn Engineering，1995，121：176-181.

［19］Sulisz W. Wave reflection and transmission at permeable breakwaters of arbitrary cross-section[J]. Coɑstɑl Engineering，1985，9：371-386.

［20］鄒國(guó)良，張慶河. 非靜壓波浪模型的無反射造波[J].海洋工程，2012，30（4）：55-61. Zou Guoliang，Zhang Qinghe. Wave generation without re-reflection for non-hydrostatic wave mode[J]. The Oceɑn Engineering，2012，30（4）：55-61（in Chinese）.

［21］Zou Guoliang，Zhang Qinghe. Improvement of absorbing boundary conditions for non-hydrostatic wave-flow model SWASH[J]. Applied Mechɑnics ɑnd Mɑteriɑls，2013，353/356：2676-2682.

［22］Stelling G S，Duinmeijer S P A. A staggered conservative scheme for every Froude number in rapidly varied shallow water flows[J]. Internɑtionɑl Journɑl for Numericɑl Methods in Fluids，2003，43：1-23.

［23］Hirayama K. Characteristics of partial reflection boundary with porous layer on the Boussinesq model[C]// Proceedings of the 20th Internɑtionɑl Offshore ɑnd Polɑr Engineering Conference. Kitakyushu，Janpan，2002：628-633.

［24］Garcia N，Lara J L，Losada I J. 2-D numerical analysis of near-field flow at low-crested permeable breakwaters[J]. Coɑstɑl Engineering，2004，51（10）：991-1020.

［25］Lara J L，Garcia N，Losada I J. RANS modelling applied to random wave interaction with submerged permeable structures[J]. Coɑstɑl Engineering，2006，53（5/6）：395-417.

（責(zé)任編輯：樊素英）

Simulation of Wave Interaction with Porous Structures Based on Non-Hydrostatic Equation

Zou Guoliang1, 2，Zhang Qinghe1
（1. State Key Laboratory of Hydraulic Engineering Simulation and Safety，Tianjin University，Tianjin 300072，China；2. Nanjing Hydraulic Research Institute，Nanjing 210024，China）

Abstract：A mathematical non-hydrostatic wave model describing pore fluid motion is developed by introducing Darcy-Forchheimer equation into the non-hydrostatic equations. The proposed model is applied to simulating wave interaction with porous structures. Then the model is verified by using experimental data of physical experiments. Results show that the model can describe wave deformation before and after porous structures reasonably with as few as two to three vertical layers. The method provides a new solution to describing wave interaction with structures.

Keywords：non-hydrostatic equation；porous structure；partial reflection；rubble mound breakwater

通訊作者：張慶河，qhzhang@tju.edu.cn.

作者簡(jiǎn)介：鄒國(guó)良（1983— ），男，博士，guoliangzou@tju.edu.cn.

基金項(xiàng)目：國(guó)家高技術(shù)研究發(fā)展計(jì)劃（863計(jì)劃）資助項(xiàng)目（2012AA051709）；天津市自然科學(xué)基金重點(diǎn)資助項(xiàng)目（12JCZDJC30200）.

收稿日期：2014-06-16；修回日期：2015-02-11.

中圖分類號(hào)：TV139.2

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：0493-2137（2016）01-0065-08

DOI:10.11784/tdxbz201406045