基于SPH方法的明基床開孔沉箱數值模擬研究

任喜峰，孫昭晨，梁書秀

(大連理工大學 海岸及近海工程國家重點實驗室，大連 116024)

海洋環境復雜多變，防波堤作為消減、抵御外海傳播到近岸的波浪的基礎設施，是設計港口、碼頭必不可少的一部分。傳統的防波堤包括斜坡堤、重力式沉箱防波堤等，隨著港口需求的增加、建港條件的日益復雜，傳統防波堤在經濟性和適用性等方面均不能滿足新的需求。隨著科學工作者對波浪與結構物之間相互作用認識的逐步深入，多種新型防波堤不斷出現。例如：開孔沉箱式防波堤、圓形防波堤和梳式防波堤等。

自1961年Jarlan[1]提出開孔沉箱以來，國內外大量采用了這種形式的防波堤來減小波浪的反射及其對建筑物的作用力。對于波浪與開孔沉箱相互作用的研究，目前主要有物模試驗和數值模擬兩種途徑。陳雪峰等[2]利用物理模型試驗研究了開孔沉箱在規則波作用下的反射率、總水平力與相對水深、相對消浪室寬度等主要影響因素之間的關系；行天強等[3]利用物理模型試驗研究了明基床上開孔沉箱在規則波作用下的反射系數情況；朱大同[4]得出了波浪與開孔沉箱相互作用的反射率的嚴格封閉解析公式；陳雪峰等[5]采用改進的VOF方法結合k-ε模型建立二維數值波浪水槽，完成了線性規則波和不規則波與單層開孔沉箱相互作用的數值模擬。

但是，由于波浪與開孔沉箱結構相互作用過程中涉及到波浪破碎、湍流等多種復雜的水動力問題，現有的理論分析和實驗分析尚不能完全揭示其作用機理。相較于傳統的VOF方法，SPH方法是一種純拉格朗日性質的無網格方法，對于處理大變形、自由表面追蹤等復雜非線性問題具有一定的優勢。姜峰等[6]基于SPH方法建立了波浪與平底開孔沉箱作用的二維數值模型。任喜峰[7]等采用域粒子著色法(Color Domain Particle，CDP)改進了SPH方法對于計算流體與開孔沉箱等薄壁結構的邊界處理方法。對于帶有明基床的開孔沉箱，由于明基床的存在導致波浪在與開孔沉箱作用前已經發生了非線性變形，使得波浪與開孔式沉箱的作用更加復雜。因此，有必要對明基床上開孔式沉箱與波浪的相互作用問題進行進一步的研究。

本文基于微可壓縮光滑粒子流體動力學方法(WCSPH)，通過域粒子著色(CDP)技術，對方法進行加強。應用加強的CDP-SPH方法，建立了模擬波浪與明基床上開孔沉箱相互作用的數值模型，通過與實驗數據對比驗證了方法的適用性和可靠性。最后應用建立的數值模型分析了規則波作用下帶有明基床的開孔式沉箱的反射率和水平力與消浪室寬度的關系。

1 數值計算方法

1.1 控制方程

基本的控制方程為連續性方程和拉格朗日形式的Navier-Stokes方程

(1)

(2)

式中：ρ為流體密度；u為速度矢量；t為時間；p為壓力；f為體積力矢量；為向量微分算子。

SPH方法的基本思想是插值。對于任意的函數A(r)，其在i點的值，可以通過對其周圍的N個點插值得到，即

(3)

式中：i,j分別為不同位置的下標；W(rij,h)為核函數；r是位置矢量；rij為從位置i到位置j的相對位置矢量；h為光滑長度；m質點質量。

因此，離散形式的控制方程可以表示為

(4)

(5)

式中：Fij為粒子關系函數，將在章節1.2給出；I為為了降低數值計算中壓力波動所引入的數值耗散項；II為為了數值穩定引入的人工粘性項；c0為參考聲速；其取值一般為10倍的最大粒子速度，以保證粒子的密度變化率小于1%；c為聲速；cij=max(ci,cj)；nij為由粒子i指向粒子j的單位方向矢量；α為人工粘性項系數，由于本文主要研究波浪問題，取值過大容易造成波浪在傳播的過程中衰減過大，因此在計算中取值為0.02；πij=(uj-ui)·rji/|rij|2。

對于微可壓縮SPH方法，粒子的壓力由狀態方程給出，即

(6)

式中：γ為常數，對于水一般可取7。

1.2 域粒子著色方法

圖1 薄板導致的支持域分割Fig.1 Thin structure divides particle support domain

SPH方法基于插值的概念，一個位置處的函數值通過其周圍粒子的插值得到。包含周圍粒子的域稱為粒子的支持域。在數值計算中，支持域通常取為φ(|rij|

圖2 薄板端點附近區域劃分示意圖Fig.2 Sketch of regions near deck end

當進行流體與結構物相互作用的模擬時，如果結構物為薄壁機構且結構物兩側都存在流體，例如透空式沉箱的前面板，雙層板防波堤等，為了避免物理上不存在的結構物兩側粒子的直接相互作用(圖1)，需根據選擇邊界條件的不同要求結構物的厚度 ?>2R(虛粒子法)或?

為了避免由于薄板導致的不必要的粒子細分，這里我們采用域粒子著色方法(CDP)，即先對計算域進行分區，然后對于邊界區域內的粒子進行著色處理，使粒子攜帶它能夠與之作用的其他區域信息。在計算粒子之間的關系時檢查兩個粒子來自的域，并判斷兩個粒子是否能夠發生相互作用,因此，公式(4)和公式(5)中的函數Fij可以取為如下的形式

表1 薄板端點附近分區相互作用關系表Tab.1 Domain relationship of a singular deck end

注：表中帶有’的編號表示區域關于結構對稱的鏡像粒子域。

(7)

圖2給出了薄板一端的域劃分示意圖，表1給出了區域之間的關系表，表中X′為與X關于薄板對稱的固定鏡像粒子[8]域。通過這樣的計算域劃分和域關系設定，能夠避免區域1和區域3的直接相互作用。需要指出的是，這樣的設定并不唯一，對于區域復雜的情況，可以根據需要合理的劃分區域和設定相互作用關系表。

1.3 邊界條件

文中除特別說明外，所有的固壁邊界均設置為自由滑移邊界條件。邊界處理采用固定鏡像粒子法[8]。

2 數值水槽的建立

2.1 主動吸收式波浪數值水槽的建立與驗證

數值水槽的長度直接影響計算量的大小，為了減少計算量并消除二次反射，數值水槽的造波邊界設置為主動吸收式推板造波邊界，根據高睿等[9]的研究，造波板的水平運動速度可以寫成

(8)

式中：X0為造波板的位置；U為造波板的水平速度；η0和η分別為造波板處的理論波面和實際波面；ω為波的圓頻率；φ為傳遞函數，使用一階理論，它的形式可以通過下式給出

(9)

式中：k為波數；d為水深。

圖3 水槽中波腹點(x=10 m)和波節點(x=9.230 m) 的波面時間序列(H=10 cm)Fig.3 Time series of free surface at node(x=10 m) and antinode (x=9.230 m) in the wave tank (H=10 cm)

為了驗證數值水槽的造波和消波性能，取計算域為平底水槽，水槽長度L=10 m，水深d=0.5 m，波高H=10 cm，周期T=1.6 s。設置水槽的右壁面為光滑直立壁面，波浪在右側壁面產生完全反射，入射波和反射波疊加后會在水槽內形成駐波。圖給出了駐波波腹點和波節點歷時曲線，結果顯示在腹點處，波浪的振幅能夠達到2倍的入射波波高，波浪節點處水面的振幅很小。但是由于波高的增大，波浪的非線性導致波谷變緩，波峰變陡。

2.2 明基床上透空沉箱參數設置

數值水槽長11.2 m，水深d=0.40 m，水槽末端為基床，高度hm=0.10 m，坡度1:2，透空式沉箱座于基床上，沉箱距離造波板9.10 m，沉箱距離基床頂面前緣0.25 cm，沉箱高度0.60 cm, 消浪室寬度Bc=0.30 m，如圖4所示。透空式沉箱的幾何尺寸和壓力測點的位置見圖5，沉箱開孔率ε=0.3。

圖4 二維明基床透空式沉箱數值水槽布置圖Fig.4 Sketch of perforated caisson sitting on rubble mound foundation in two-dimensional numerical flume圖5 開孔沉箱幾何尺寸和壓力測點布置圖Fig.5 Geometry of perforated caisson and layout of pressure sensors

3 數值計算結果及分析

3.1 測點壓力過程分析

透空式沉箱消浪室內外的水體僅能通過開孔進行交換，沉箱內外的水壓差導致水體高速通過開孔，進而導致開孔附近產生較強的湍混，水流結構復雜。這里僅以周期T=1.2 s, 波高H=6 cm的波浪與開孔沉箱的作用過程進行分析。圖6給出了1#、3#、4#、7#和12#測點的數值模擬結果與實驗結果的比較。從圖中可以看出，相較于VOF計算結果，采用本文SPH模型的數值計算結果與實驗結果更加吻合，尤其是對于靜水面以上的壓力測點(4#、12#)，SPH的計算結果與實驗值吻合非常好，較VOF計算結果優勢明顯。

圖6 1#、3#、4#、7#和12#測點的數值模擬結果于實驗結果的比較Fig.6 Comparison of simulated results and experimental results at sensor 1#, 3#, 4#, 7# and 12#

圖7給出了t=13.8 s時開孔消浪室附近自由面位置對比。從圖中可以看出，SPH模擬結果的自由面較不規整；在消浪室后墻處，自由位置高出VOF模擬結果較多，這與12#測點SPH方法的點壓力計算結果大于VOF的點壓力計算結果相吻合。這也說明了，SPH方法對于處理流體與開孔沉箱的作用過程中發生的自由面大變形問題相較于傳統的VOF具有一定的優勢。

3.2 反射系數的影響因素分析

反射系數是評價開孔沉箱性能的重要指標之一，通過本文建立的SPH模型，對波高分別為0.06 m和0.08 m，周期為1.0 s、1.2 s、1.4 s和1.6 s, 基床高度為0.1 m的工況組合進行數值模擬，提取1#～5#位置處的波面歷時數據，采用Goda等[10]的兩點法分離輸入波和反射波并計算反射系數。圖8給出了T=1.2 s、H=0.06 m時的分離的入射波與反射波能量譜。從圖中可以看出，入射波為單色波，頻率為fp=0.384 Hz，與造波板生成的入射波一致。反射波的能量主要集中在一倍頻和二倍頻處，且二倍頻的能量更顯著。

圖9 反射系數與相對消浪室寬度的關系Fig.9 Reflection coefficient, Kr vs. Bc/L圖10 反射系數與相對基床高度的關系Fig.10 Reflection coefficient, Kr vs. hm/L

圖11 SPH計算結果與經驗公式的反射系數比較Fig.11 Comparison of reflection coefficient of results from present SPH method and empirical formula

圖9、圖10給出了反射系數與消浪室寬度(Bc/L)、相對基床高度(hm/L)的關系。從圖中可以看出，透空式沉箱的反射系數隨相對消浪室寬度的增加先減小后增大；隨相對基床高度的增加先減小后增大。從圖10中還可以看出，相對基床高度較大時，大波高的反射系數更大，原因可以歸結于基床淺化作用導致波浪變形，大波高的波浪的非線性更強。

行天強等[3]通過實驗研究認為透空式沉箱的反射系數與相對消浪室寬度、相對基床高度、相對水深(d/L)以及開孔率(ε)有關，擬合了規則波作用下明基床上開孔沉箱的反射系數計算公式

(10)

圖11給出了數值計算結果與上述經驗公式的比較，可以看出，SPH方法的數值計算結果與實驗公式基本吻合，大多數點在y=(1±10%)x的包絡線內。所以，本文所建立的SPH模型計算結果可以為分析開孔沉箱反射率提供參考依據。

4 結論

本文基于光滑粒子流體動力學(SPH)方法，并應用域粒子著色(CDP)法，建立了二維數值水槽。進而對波浪與明基床上開孔沉箱的相互作用進行數值模擬。數值計算結果在各壓力測點與實驗值吻合良好，相較于VOF計算結果，靜水面以上的壓力測點優勢更加明顯。通過對比同一時刻VOF計算結果與SPH結果的自由表面，可以發現SPH方法模擬結果對大自由面變形適應更好。最后進行了多組合工況的數值模擬，分析得出開孔防波堤的反射系數隨相對基床高度、相對消浪室寬度的增加呈現先減小后增大的趨勢。通過與實驗得到的經驗公式比對，可知數值計算結果與實驗結果吻合良好。綜上，本文應用改進的CDP-SPH方法，解決了開孔沉箱與流體相互作用的問題，為開孔沉箱的數值模擬提供了新的手段。