波浪與拋石潛堤相互作用的MPS法數值模型研究

張永蘭蔣 勤*王麗珠

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

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

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

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

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

1 MPS紊流滲流模型

1.1 控制方程

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

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

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

式中,μ為動力黏性系數;kp為滲透性系數;φ為體積分數;Dc為拋石塊體的平均粒徑,本文取Dc=d50。

當采用MPS法進行數值求解時,首先利用核函數(或稱權重函數)建立粒子間相互作用模型。然后,對控制方程進行離散求解,以實現對流體運動過程的數值模擬。MPS法以目標粒子與其相鄰粒子間的距離遠近來描述粒子間相互影響的程度,最常用的核函數為Koshizuka等[18]于1996年提出的標準核函數w(rij),具體公式為

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

通過構建梯度(?φ)算子、拉普拉斯(?2φ)算子和散度)的離散模型對控制微分方程進行離散,實現對流體運動過程的數值求解。常用的算子的離散模型為

式中,φ為任意的一個標量為任意的一個矢量;Ds為數值計算模型的空間維度,對二維數值計算模型Ds取值為2;n0為初始粒子數密度;λ為一個模型參數,定義為

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

1.2 壓力泊松方程

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

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

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

式中,Δt為計算時間步長。在MPS法中,預測步對應于虛擬時間步k+1/2,用“*”表示。在校正步中,目標粒子i在k+1時步的流速為中間速度矢量和源于壓力梯度的修正流速矢量之和,即

對式(11)兩邊同時求散度,可得

結合式(1)可得

最后,綜合式(13)和式(14),可以得到可滲結構物內的壓力泊松方程為

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

1.3 壓力梯度模型

為了避免因原始壓力梯度模型計算精度低而引起的數值能量耗散問題,本研究采用Wang等[16]基于泰勒級數展開提出的高精度壓力梯度模型,其表達式為

1.4 自由表面粒子的識別

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

式中,β為經驗系數,一般取值為0.80～0.97,本研究β=0.95。

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

1.5 流-固界面的耦合

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

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

2 模型驗證與結果分析

2.1 “U”型管內多孔介質滲流

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

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

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

同理,多孔介質內滲流速度隨時間變化的理論公式為

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

式中,t為時間;ΔH0初始水位差;LP為滲流路徑的長度;u為多孔介質內的達西滲流速度。

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

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

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

2.2 拋石潛堤上孤立波傳播

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

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

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

圖6 不同時刻(t=1.45 s,t=1.85 s,t=2.05 s,t=2.25 s)潛堤周圍壓力場和流速場模擬結果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)

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

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

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

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

表1 流速測量點x 方向坐標表Table 1 The x coordinates of the flow velocity measuring points

2.3 模型收斂性驗證

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

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

3 結 論

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

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

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

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