張永蘭蔣 勤*王麗珠

(1.河海大學(xué) 港口海岸與近海工程學(xué)院,江蘇 南京210098;2.南京師范大學(xué) 海洋科學(xué)與工程學(xué)院,江蘇 南京210023)

拋石潛堤具有結(jié)構(gòu)簡(jiǎn)單,易于施工,整體穩(wěn)定性好等優(yōu)點(diǎn),所以它是目前海岸工程中廣泛應(yīng)用的防護(hù)結(jié)構(gòu)型式。當(dāng)波浪行進(jìn)到拋石潛堤時(shí),部分波浪會(huì)在堤前被反射,部分波浪會(huì)越過潛堤,還有一部分波浪水體會(huì)以滲流的形式透過拋石堤的孔隙傳入掩護(hù)區(qū)內(nèi)。由于在拋石潛堤孔隙中的滲流水體受到塊體摩擦阻力的影響,部分入射波能被耗散,因此拋石防波堤可以起到很好的消波作用。因此,研究波浪與拋石防波堤的相互作用對(duì)防波堤設(shè)計(jì)具有重要的工程價(jià)值。

近年來,隨著計(jì)算機(jī)模擬技術(shù)的發(fā)展,數(shù)值模擬方法已成為研究波浪與結(jié)構(gòu)物相互作用的重要手段之一。波浪與拋石潛堤相互作用問題是一個(gè)具有復(fù)雜流-固界面和大自由表面變形,以及多孔介質(zhì)滲流的強(qiáng)非線性紊流運(yùn)動(dòng)問題。鑒于對(duì)拋石體內(nèi)復(fù)雜的紊動(dòng)滲流問題進(jìn)行嚴(yán)格的理論求解存在一定的困難,目前通常將拋石結(jié)構(gòu)簡(jiǎn)化為均質(zhì)多孔介質(zhì),通過建立數(shù)值計(jì)算模型的方法模擬波浪與可滲結(jié)構(gòu)物的相互作用。對(duì)多孔介質(zhì)內(nèi)滲流的研究可追溯到19世紀(jì)初期,Darcy[1]從流量與壓降關(guān)系的線性假定出發(fā),提出了分析多孔介質(zhì)內(nèi)低雷諾數(shù)流體運(yùn)動(dòng)的基本原理。Forchheimer[2]針對(duì)具有較高雷諾數(shù)的滲流問題,在動(dòng)量方程中引入一個(gè)附加非線性阻力項(xiàng),建立了描述非達(dá)西滲流運(yùn)動(dòng)的控制方程。Sollitt和Cross[3]通過在動(dòng)量方程中添加慣性項(xiàng)和線性、非線性阻力項(xiàng),提出了描述多孔介質(zhì)中非達(dá)西流動(dòng)的廣義滲流模型。Gent[4]基于描述流體運(yùn)動(dòng)的雷諾時(shí)均N-S方程,并借助“U”型管內(nèi)一維滲流的實(shí)驗(yàn)結(jié)果來調(diào)整模型中線性、非線性阻力系數(shù)的方法,研究了波浪與多孔介質(zhì)的相互作用問題。Liu等[5]基于流體運(yùn)動(dòng)的N-S方程,推導(dǎo)出描述多孔介質(zhì)內(nèi)滲流運(yùn)動(dòng)的空間平均N-S方程。Hsu等[6]在Liu等[5]的基礎(chǔ)上給出了可滲結(jié)構(gòu)內(nèi)部紊流運(yùn)動(dòng)的體平均/雷諾平均N-S方程,以描述介質(zhì)的孔隙尺寸較大、紊流效應(yīng)無法忽略的紊流滲流問題。迄今為止,波浪與拋石結(jié)構(gòu)物相互作用的數(shù)值模擬大多基于上述多孔介質(zhì)滲流模型,并采用網(wǎng)格法進(jìn)行數(shù)值求解[6-9]。但是,傳統(tǒng)網(wǎng)格法在對(duì)動(dòng)量方程中的對(duì)流項(xiàng)和大變形自由表面運(yùn)動(dòng)求解時(shí)存在著數(shù)值耗散問題,難以準(zhǔn)確再現(xiàn)具有復(fù)雜流-固界面和大變形自由面的強(qiáng)非線性滲流問題。

粒子法作為一種拉格朗日方法,具有易于處理大自由表面變形和不同介質(zhì)界面以及由于在控制方程中不包含對(duì)流項(xiàng),所以在理論上能夠避免數(shù)值耗散等優(yōu)勢(shì),可以彌補(bǔ)傳統(tǒng)網(wǎng)格法的不足。因此,近年來,采用粒子法模擬多孔介質(zhì)滲流問題的研究備受關(guān)注。Akbari和Namin[10]以及Shao和Loe[11]忽略孔隙介質(zhì)內(nèi)外流體運(yùn)動(dòng)的紊流效應(yīng),分別采用光滑粒子流體動(dòng)力學(xué)法(Smoothed Particle Hydrodynamics,SPH)和不可壓縮光滑粒子流體動(dòng)力學(xué)法(Incompressible SPH,ISPH)求解可滲結(jié)構(gòu)物內(nèi)外的流體運(yùn)動(dòng)問題。同樣,Fu和Jin[12]忽略滲流的紊流效應(yīng),利用移動(dòng)粒子半隱式法(Moving Particle Semi-implicit method,MPS)建立的數(shù)值模型研究了多孔介質(zhì)河床上的水流流動(dòng)問題。此外,考慮到波浪與可滲結(jié)構(gòu)相互作用中常伴有波浪破碎等強(qiáng)紊動(dòng)現(xiàn)象,Ren[13]忽略孔隙介質(zhì)內(nèi)部的紊流運(yùn)動(dòng),將大渦模擬技術(shù)用于再現(xiàn)孔隙介質(zhì)外部的紊流運(yùn)動(dòng),建立了分析波浪與可滲結(jié)構(gòu)物相互作用的弱可壓光滑粒子流體動(dòng)力學(xué)法(Weakly Compressible SPH,WCSPH)數(shù)值計(jì)算模型。Peng等[14]考慮孔隙介質(zhì)內(nèi)外的紊流運(yùn)動(dòng),基于Drew[15]的二相流混合理論,建立了模擬多孔介質(zhì)中流體流動(dòng)的WCSPH 模型。

相較于網(wǎng)格法,無網(wǎng)格粒子法是一種較新的數(shù)值流體計(jì)算方法,尚存在自由表面粒子誤判、壓力求解精度較低以及在模擬波浪傳播時(shí)數(shù)值能量衰減顯著等問題[16]。因此,本文利用改進(jìn)的MPS法流體運(yùn)動(dòng)模型[16,17]構(gòu)建了垂向二維紊流滲流數(shù)值計(jì)算模型,模擬波浪與拋石結(jié)構(gòu)物的相互作用。通過對(duì)“U”型管中多孔介質(zhì)內(nèi)滲流運(yùn)動(dòng)和孤立波在拋石潛堤上傳播過程的數(shù)值模擬及其與理論解和實(shí)測(cè)結(jié)果的對(duì)比分析,對(duì)所建立的MPS法紊流滲流模型在用于波浪與可滲結(jié)構(gòu)物相互作用問題時(shí)的求解精度進(jìn)行驗(yàn)證。

1 MPS紊流滲流模型

1.1 控制方程

MPS(Moving Particle Semi-implicit)法即移動(dòng)粒子半隱式法是由Koshizuka等[18]于1996年提出的一種無網(wǎng)格數(shù)值計(jì)算方法。它采用帶有物理信息的粒子將計(jì)算域進(jìn)行離散,通過追蹤離散粒子的拉格朗日運(yùn)動(dòng),模擬流體的運(yùn)動(dòng)過程。因此,在采用MPS 法建立模擬多孔介質(zhì)內(nèi)外流體運(yùn)動(dòng)的數(shù)值模型時(shí),基于Drew[15]的二相混合流理論,可以得出描述斷面二維多孔介質(zhì)內(nèi)外流體運(yùn)動(dòng)的控制方程為

式中,u為速度矢量;ρ為密度;t為時(shí)間;P為壓強(qiáng);g為重力加速度;ν為運(yùn)動(dòng)黏度;φ為孔隙介質(zhì)內(nèi)流體所占體積的百分比,即：φ=VF/V,VF為液相體積,V為多孔介質(zhì)的體積;R為多孔介質(zhì)對(duì)孔隙水流的阻力;νt為紊流黏性系數(shù),可采用Smagorinsky紊流模型[19]計(jì)算,即νt=(CsΔl)2(2SijSij)1/2,其中Cs為Smagorinsky常數(shù),在本模型中取Cs的值為0.1,Δl為粒子初始間距,Sij為應(yīng)變張量,i=1,2,j=1,2。

在上述基于二相混合流理論[15]的方程中,包含了流體和固體(多孔介質(zhì)骨架)兩相。本文假定多孔介質(zhì)骨架不隨外力變形,因此無需求解固相運(yùn)動(dòng)。同時(shí),假定多孔介質(zhì)內(nèi)流體處于飽和狀態(tài),流體穿過多孔介質(zhì)區(qū)所受到的阻力為R,根據(jù)Peng等[14]的計(jì)算方法,R的計(jì)算公式為

式中,μ為動(dòng)力黏性系數(shù);kp為滲透性系數(shù);φ為體積分?jǐn)?shù);Dc為拋石塊體的平均粒徑,本文取Dc=d50。

當(dāng)采用MPS法進(jìn)行數(shù)值求解時(shí),首先利用核函數(shù)(或稱權(quán)重函數(shù))建立粒子間相互作用模型。然后,對(duì)控制方程進(jìn)行離散求解,以實(shí)現(xiàn)對(duì)流體運(yùn)動(dòng)過程的數(shù)值模擬。MPS法以目標(biāo)粒子與其相鄰粒子間的距離遠(yuǎn)近來描述粒子間相互影響的程度,最常用的核函數(shù)為Koshizuka等[18]于1996年提出的標(biāo)準(zhǔn)核函數(shù)w(rij),具體公式為

式中,re為目標(biāo)粒子的影響域半徑;rij為粒子i和粒子j之間的距離,即rij=|rj-ri|,ri和rj分別為i,j粒子的位移矢量;。參考以往的研究經(jīng)驗(yàn),對(duì)于梯度模型,取re=2.1l0;對(duì)拉普拉斯模型,取re=3.1l0,在此l0表示初始粒子分布的粒子間距離,即離散粒子的直徑。

通過構(gòu)建梯度(?φ)算子、拉普拉斯(?2φ)算子和散度)的離散模型對(duì)控制微分方程進(jìn)行離散,實(shí)現(xiàn)對(duì)流體運(yùn)動(dòng)過程的數(shù)值求解。常用的算子的離散模型為

式中,φ為任意的一個(gè)標(biāo)量為任意的一個(gè)矢量;Ds為數(shù)值計(jì)算模型的空間維度,對(duì)二維數(shù)值計(jì)算模型Ds取值為2;n0為初始粒子數(shù)密度;λ為一個(gè)模型參數(shù),定義為

粒子數(shù)密度n是流體密度在 MPS 法中的一個(gè)衍生變量,在物理含義上兩者是對(duì)等的。MPS法通過在計(jì)算過程中保持流體的粒子數(shù)密度不變來滿足流體的不可壓縮條件。粒子數(shù)密度n是一個(gè)無量綱參數(shù),定義為目標(biāo)粒子i與其影響域內(nèi)所有相鄰粒子j之間的核函數(shù)值的累加值,具體表示為

1.2 壓力泊松方程

與SMAC法類似,MPS法采用映射法對(duì)離散后的控制方程進(jìn)行求解。針對(duì)某一時(shí)步,模型的求解過程分為兩步：第一步為預(yù)測(cè)步,對(duì)動(dòng)量方程中除壓力梯度項(xiàng)以外的所有項(xiàng)進(jìn)行積分,得到一個(gè)臨時(shí)的速度場(chǎng)和位移場(chǎng);第二部為校正步,考慮流體的不可壓縮條件,對(duì)包括預(yù)測(cè)步中未使用的壓力梯度項(xiàng)的動(dòng)量方程進(jìn)行積分,推導(dǎo)出壓力泊松方程,通過求解壓力泊松方程得到新的壓力場(chǎng),再采用壓力梯度項(xiàng)對(duì)預(yù)測(cè)步得到的臨時(shí)速度場(chǎng)和位移場(chǎng)進(jìn)行修正,求出下一時(shí)步物理量的計(jì)算結(jié)果。

具體地,在采用MPS法進(jìn)行數(shù)值求解時(shí),目標(biāo)粒子i在k+1時(shí)步的流速可以拆分為中間流速和校正流速)之和,即

式中,i代表目標(biāo)粒子。根據(jù)由Drew[15]的二相流混合理論得到的動(dòng)量方程式(2),在預(yù)測(cè)步中,目標(biāo)粒子i的中間流速可以采用上一時(shí)步k得到的重力、黏性力以及阻力項(xiàng)的值顯式求解,即

式中,Δt為計(jì)算時(shí)間步長(zhǎng)。在MPS法中,預(yù)測(cè)步對(duì)應(yīng)于虛擬時(shí)間步k+1/2,用“*”表示。在校正步中,目標(biāo)粒子i在k+1時(shí)步的流速為中間速度矢量和源于壓力梯度的修正流速矢量之和,即

對(duì)式(11)兩邊同時(shí)求散度,可得

結(jié)合式(1)可得

最后,綜合式(13)和式(14),可以得到可滲結(jié)構(gòu)物內(nèi)的壓力泊松方程為

式中,ρ為流體密度;n0為初始粒子數(shù)密度為k時(shí)步的粒子數(shù)密度為k時(shí)步粒子i的體積分?jǐn)?shù)其中,nw為可滲結(jié)構(gòu)物孔隙率,Ωi為目標(biāo)粒子i的鄰域搜索范圍,其直徑等于可滲區(qū)域過渡區(qū)厚度,w(rij)為表征目標(biāo)粒子i與其相鄰粒子j之間相互關(guān)系的核函數(shù)。

1.3 壓力梯度模型

為了避免因原始?jí)毫μ荻饶Ｐ陀?jì)算精度低而引起的數(shù)值能量耗散問題,本研究采用Wang等[16]基于泰勒級(jí)數(shù)展開提出的高精度壓力梯度模型,其表達(dá)式為

1.4 自由表面粒子的識(shí)別

與網(wǎng)格法中需要采用特殊的自由表面追蹤方法相比,無網(wǎng)格粒子法的最大優(yōu)勢(shì)之一是可以自然、簡(jiǎn)單地實(shí)現(xiàn)自由表面的追蹤模擬。在原始MPS 法中,自由表面粒子的識(shí)別依據(jù)自由表面處粒子的影響域被自由表面邊界截?cái)?其粒子數(shù)密度小于內(nèi)部流體粒子的粒子數(shù)密度的特點(diǎn),給出的自由表面粒子的判別條件為：

式中,β為經(jīng)驗(yàn)系數(shù),一般取值為0.80～0.97,本研究β=0.95。

式(17)中的自由表面粒子的判定條件簡(jiǎn)單易行,但是,當(dāng)內(nèi)部流體粒子的分布不夠均勻、混亂無序時(shí),內(nèi)部的流體粒子常常被誤判為自由表面粒子,從而引起嚴(yán)重的非物理壓力震蕩問題。因此,本研究采用Wang等[19]提出的壓力判別法來提高自由表面的識(shí)別精度,即：在式(17)的自由表面粒子判別條件的基礎(chǔ)上,引入以壓力值為參數(shù)的附加判別標(biāo)準(zhǔn)：

1.5 流-固界面的耦合

為了保證流體運(yùn)動(dòng)的連續(xù)性,需要對(duì)流體與多孔介質(zhì)間的流-固界面進(jìn)行特殊的邊界處理。在傳統(tǒng)歐拉網(wǎng)格法中,通常依據(jù)連續(xù)方程和動(dòng)量守恒方程,在流-固界面網(wǎng)格上施加相應(yīng)的流速或壓力等的連續(xù)條件。而在MPS粒子法中,由于離散流體粒子具有拉格朗日特性,粒子的位置坐標(biāo)是隨流體運(yùn)動(dòng)不斷變化的,因此無法像傳統(tǒng)網(wǎng)格法那樣在離散粒子上施加相應(yīng)的邊界條件。因此,本研究借鑒Akbari和Namin[10]的邊界處理方法,在流體域和多孔介質(zhì)域之間設(shè)置一個(gè)孔隙率漸變的過渡區(qū)域(圖1),使得流體粒子進(jìn)出多孔介質(zhì)時(shí)受到一個(gè)連續(xù)變化的阻力作用,以實(shí)現(xiàn)多孔介質(zhì)內(nèi)外流體運(yùn)動(dòng)的連續(xù)性。由于實(shí)際問題中可滲結(jié)構(gòu)物界面一般具有不規(guī)則的幾何形狀,難以精準(zhǔn)再現(xiàn)可滲界面,因此,可滲界面通常被視為具有一定范圍的一個(gè)流-固過渡區(qū)域,這一過渡區(qū)域的大小可根據(jù)組成可滲結(jié)構(gòu)物的材質(zhì)的平均直徑來確定。前人的研究結(jié)果[19]表明：當(dāng)選取的過渡區(qū)域的厚度δ小于可滲結(jié)構(gòu)物材料的平均粒徑(δ＜d50)時(shí),將導(dǎo)致數(shù)值計(jì)算的不穩(wěn)定;當(dāng)選取的過渡區(qū)域厚度δ大于可滲結(jié)構(gòu)物材料的平均粒徑(δ＞d50)時(shí),所選取的過渡區(qū)域厚度越小,計(jì)算結(jié)果的精度卻越高。通過模型模擬試算,本文將可滲結(jié)構(gòu)物材料中值粒徑值(δ=d50)作為流-固過渡區(qū)域的計(jì)算厚度。

圖1 流-固界面耦合示意圖Fig.1 A sketch map showing the coupling at the fluid-solid boundary

2 模型驗(yàn)證與結(jié)果分析

2.1 “U”型管內(nèi)多孔介質(zhì)滲流

對(duì)“U”型管中多孔介質(zhì)內(nèi)滲流問題進(jìn)行了數(shù)值實(shí)驗(yàn)研究,分析了 由“U”型管兩端液面差引起的多孔介質(zhì)內(nèi)不同時(shí)刻的滲流速度和液面變化過程。數(shù)值實(shí)驗(yàn)的計(jì)算條件為“U”型管長(zhǎng)B=4 m,高H=3.8 m,在“U”型管底部中央?yún)^(qū)設(shè)置了寬Lp=1 m,高Hp=1 m的多孔介質(zhì)結(jié)構(gòu)物;“U”型管左端初始水柱高HL=3.35 m,右端水柱高HR=2 m,左右端水位差ΔH0=1.35 m(圖2)。由于“U”型管兩端水柱高度不同,在重力作用下,左側(cè)水體會(huì)穿過多孔介質(zhì)滲流到“U”型管的右側(cè),當(dāng)管內(nèi)兩端液面高度持平時(shí)達(dá)到最終的平衡狀態(tài)。

圖2 “U”型管內(nèi)滲流概化圖(m)Fig.2 A sketch map of seepage flow motion within a U-shaped pipe(m)

式中,t為時(shí)間;ΔH0為“U”型管左右兩側(cè)水柱的初始水位差;LP為滲流路徑的長(zhǎng)度;kh為導(dǎo)水率,即kh=ρgkp/μ。其中g(shù)為重力加速度;kp為滲透性系數(shù);μ為流體的動(dòng)力黏性系數(shù);ρ為流體密度。

同理,多孔介質(zhì)內(nèi)滲流速度隨時(shí)間變化的理論公式為

采用線性達(dá)西滲流的假設(shè)[1],忽略式(4)中第二項(xiàng)的非線性阻力項(xiàng),可得滲流水頭ΔH0隨時(shí)間變化的理論公式為：

式中,t為時(shí)間;ΔH0初始水位差;LP為滲流路徑的長(zhǎng)度;u為多孔介質(zhì)內(nèi)的達(dá)西滲流速度。

為了與理論值進(jìn)行比較,采用所建立的MPS法滲流模型,選取離散粒子直徑l0=0.05 m,粒子總數(shù)N=5 000,流體密度ρ取水的密度1 000 kg/m3,流體的運(yùn)動(dòng)黏滯系數(shù)ν為10-6m2/s,對(duì)“U”型管中多孔介質(zhì)內(nèi)滲流進(jìn)行了數(shù)值模擬。對(duì)比分析當(dāng)kh為0.05和0.025時(shí)計(jì)算得出的“U”型管兩端水柱水頭差隨時(shí)間變化的數(shù)值解與理論值(圖3)可見：在兩端水頭差較大的初期階段(t≤10 s),由于兩端液面高度差較大,重力作用較強(qiáng),數(shù)值解與理論值均給出了液面差急速下降的變化趨勢(shì),二者間的誤差低于10%;此后,隨著兩端液面差逐漸變小,數(shù)值解與理論值給出的液面差變化速度變緩,最終達(dá)到平衡,計(jì)算值與理論值間的誤差低于20%。相應(yīng)地,對(duì)比分析管內(nèi)可滲區(qū)中心斷面平均滲流速度隨時(shí)間變化的數(shù)值解與理論值(圖4)可見：模型模擬結(jié)果與實(shí)驗(yàn)得到的滲流速度的變化趨勢(shì)與液面水頭差變化趨勢(shì)基本一致,即在初期階段(t≤10 s),平均滲流速度隨時(shí)間急速減弱;此后平均滲流速度緩慢下降。數(shù)值計(jì)算得到的平均滲流速度與理論值間的誤差低于15%。

圖3 不同時(shí)刻液面差變化計(jì)算結(jié)果與理論值比較Fig.3 Comparison between the calculated results and the theoretical values of the variations in water level difference at different time

圖4 不同時(shí)刻平均滲流速度的計(jì)算結(jié)果與理論值比較Fig.4 Comparison between the calculated results and the theoretical values of the averaged seepage velocity at different time

2.2 拋石潛堤上孤立波傳播

采用本研究建立的MPS滲流模型對(duì)孤立波在拋石潛堤上的傳播過程進(jìn)行了模擬研究,并與Wu和Hsiao[21]的物理模型實(shí)驗(yàn)結(jié)果進(jìn)行了比較。依據(jù)物理模型實(shí)驗(yàn),計(jì)算域設(shè)置如圖5所示,水槽長(zhǎng)L=8 m,水深d=0.11 m,在水槽的左端采用推板式造波,入射波高為H=0.047 7 m。在水槽中部距水槽左端4 m 處設(shè)置了長(zhǎng)B=13 cm,高h(yuǎn)=6.5 cm 的矩形拋石潛堤。組成拋石潛堤的碎石平均粒徑d50=0.015 m,孔隙率nw=0.52。計(jì)算中,將拋石堤前趾堤腳作為坐標(biāo)原點(diǎn),分別在堤前1.8 m(x=-1.8 m)處和堤后1.8 m(x=+1.8 m)處設(shè)置P1和P2兩個(gè)測(cè)波點(diǎn),記錄孤立波的波面形態(tài);離散粒子直徑取l0=0.005 m,離散粒子總數(shù)N=27 000個(gè);流體密度ρ=1 000 kg/m3,運(yùn)動(dòng)黏滯系數(shù)ν=10-6m2/s。

圖5 孤立波在可滲潛堤上傳播概化圖Fig.5 A sketch map of solitary wave propagation over submerged permeable breakwater

不同時(shí)刻(t分別為1.45,1.85,2.05和2.25 s)潛堤附近壓力場(chǎng)與流速場(chǎng)的模擬結(jié)果如圖6所示。需要說明的是本研究將孤立波波峰到達(dá)x=-1.8 m 處的瞬間定義t=0 s。由圖6a可見,在流-固界面處粒子分布均勻,且壓力分布的連續(xù)性較好,表明本研究所采用的流-固界面耦合的邊界處理方法是合理可行的。由圖6b可見：當(dāng)孤立波傳播到堤前趾附近時(shí),由于受到拋石潛堤的阻礙作用,波浪水質(zhì)點(diǎn)速度明顯減弱。在堤的后方由于波浪水流速度的垂向差異,形成了明顯的紊流渦旋,且渦旋的尺度和強(qiáng)度隨時(shí)間逐漸增大。在t=2.05 s時(shí),渦旋基本發(fā)育成熟,到t=2.25 s時(shí),渦旋充分成長(zhǎng),渦旋外圍延伸到自由液面附近。由圖6a的壓力場(chǎng)和圖6b的流速場(chǎng)的模擬結(jié)果可見：本數(shù)值模型有效地保證了壓力場(chǎng)的均勻變化,無明顯數(shù)值壓力震蕩,合理地再現(xiàn)了潛堤后方紊動(dòng)渦旋的形成發(fā)展過程。

圖6 不同時(shí)刻(t=1.45 s,t=1.85 s,t=2.05 s,t=2.25 s)潛堤周圍壓力場(chǎng)和流速場(chǎng)模擬結(jié)果Fig.6 Simulated results of the pressure fields and flow speed fields around the submerged permeable breakwater at different time(t=1.45 s,t=1.85 s,t=2.05 s,t=2.25 s)

對(duì)比由MPS數(shù)值模型得到的拋石潛堤前后孤立波波面歷時(shí)曲線與Wu和Hsiao[21]的實(shí)驗(yàn)結(jié)果和Gui等[22]ISPH 模擬結(jié)果(圖7)可知：在x= -1.8和1.8 m 處,本數(shù)值模型得到的孤立波波面形態(tài)變化與Wu和Hsiao[21]的實(shí)驗(yàn)結(jié)果基本一致,同時(shí),相較于Gui等[22]采用ISPH 法得計(jì)算結(jié)果,本模型的模擬結(jié)果與實(shí)驗(yàn)值更為接近。

圖7 拋石潛堤前后孤立波波面形態(tài)模擬結(jié)果與實(shí)測(cè)值比較Fig.7 Comparison of the calculated and the measured profiles of solitary wave in front and back of the submerged rubble-mound breakwater

圖8 潛堤周圍典型斷面水平速度垂向分布對(duì)比圖Fig.8 Comparison of vertical distributions of horizontal velocity along typical sections around permeable structures

圖8為不同時(shí)刻(t=1.45,1.85,2.05和2.25 s)典型斷面(如表1所示)的水平流速垂線分布的MPS法模型模擬結(jié)果與Wu和Hsiao[21]的模型實(shí)驗(yàn)結(jié)果和Gui等[21]的ISPH 法模擬結(jié)果的比較。相較于Gui等[22]采用ISPH 法得到的模擬結(jié)果,本模型給出的各斷面的流速分布總體上與實(shí)驗(yàn)結(jié)果更加接近。此外,從圖8b～8d可見,在斷面x=0.16 m 處,斷面流速分布的數(shù)值計(jì)算值與實(shí)驗(yàn)值間的誤差明顯大于其他6個(gè)斷面。結(jié)合圖6b可以發(fā)現(xiàn),在x=0.16 m 處的流速變化比較復(fù)雜,并伴有渦的產(chǎn)生和消失,由于目前的數(shù)值模型是將拋石結(jié)構(gòu)物假設(shè)為均勻的多孔介質(zhì)結(jié)構(gòu),不能再現(xiàn)實(shí)際拋石結(jié)構(gòu)物中復(fù)雜的孔隙結(jié)構(gòu),因此導(dǎo)致在流速變化復(fù)雜的區(qū)域模型得到的流速值與實(shí)驗(yàn)值尚存在一定的差異。

表1 流速測(cè)量點(diǎn)x 方向坐標(biāo)表Table 1 The x coordinates of the flow velocity measuring points

2.3 模型收斂性驗(yàn)證

為了驗(yàn)證所建立的數(shù)值計(jì)算模型的收斂性,采用不同的離散粒子解像度對(duì)計(jì)算結(jié)果的模擬精度進(jìn)行了分析。圖9為粒子直徑分別為l0=0.005 m,l0=0.006 m 和l0=0.008 m 的條件下得到的t=2.25 s時(shí)刻,不同斷面(x= -0.04 m,x= +0.04 m,x= +0.12 m)上水平流速垂向分布計(jì)算結(jié)果的對(duì)比圖。由圖9可見,3種不同離散粒子尺度下的計(jì)算結(jié)果差異很小,并與實(shí)驗(yàn)值基本吻合,表明所建立的MPS滲流模型是收斂的。

圖9 不同粒子解像度條件下水平流速計(jì)算結(jié)果比較Fig.9 Calculated horizontal flow velocities under the conditions of different resolutions

3 結(jié) 論

本文以Wang等[16,19]提出的改進(jìn)MPS法流體運(yùn)動(dòng)模型為基礎(chǔ),基于Drew[15]的多孔介質(zhì)二相流運(yùn)動(dòng)方程,采用Gotoh等[23]提出的亞粒子尺度(SPS)紊流模型,構(gòu)建了垂向二維MPS法紊流滲流模型,用以研究波浪與拋石潛堤相互作用。通過模擬“U”型管中多孔介質(zhì)內(nèi)滲流運(yùn)動(dòng)和孤立波在拋石潛堤上的傳播過程,驗(yàn)證了本模型的合理性和模擬精度,具體結(jié)論為：

1)對(duì)“U”型管中多孔介質(zhì)內(nèi)滲流問題的模擬研究表明：模型得到的不同時(shí)刻“U”型管兩側(cè)水柱液面差和管內(nèi)滲流區(qū)中心斷面水平速度平均值變化與理論值基本吻合,變化趨勢(shì)一致,證明本文提出的MPS法紊流滲流模型可已很好地模擬線性滲流問題。

2)由孤立波在拋石潛堤上傳播過程的模擬結(jié)果可見：由本文建立的數(shù)值模型得到的孤立波波面形態(tài)、滲流區(qū)內(nèi)外典型斷面上垂向流速分布的計(jì)算結(jié)果與實(shí)測(cè)結(jié)果吻合較好;計(jì)算得到的壓力場(chǎng)分布均勻、無明顯的數(shù)值壓力震蕩;流速場(chǎng)的計(jì)算結(jié)果清晰地再現(xiàn)出波浪與拋石潛堤相互作用過程中的渦流現(xiàn)象及其運(yùn)動(dòng)過程,表明本模型在模擬具有復(fù)雜介質(zhì)界面和大自由表面變形的非線性滲流問題時(shí)具有顯著優(yōu)勢(shì)。

3)通過對(duì)不同解像度條件下滲流區(qū)典型斷面水平流速垂向分布的模擬結(jié)果與實(shí)測(cè)值的比較,驗(yàn)證了本數(shù)值模型具有良好的收斂性。