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

任喜峰，孫昭晨，梁書秀

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

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

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

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

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

1 數(shù)值計算方法

1.1 控制方程

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

(1)

(2)

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

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

(3)

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

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

(4)

(5)

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

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

(6)

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

1.2 域粒子著色方法

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

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

圖2 薄板端點(diǎn)附近區(qū)域劃分示意圖Fig.2 Sketch of regions near deck end

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

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

表1 薄板端點(diǎn)附近分區(qū)相互作用關(guān)系表Tab.1 Domain relationship of a singular deck end

注：表中帶有’的編號表示區(qū)域關(guān)于結(jié)構(gòu)對稱的鏡像粒子域。

(7)

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

1.3 邊界條件

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

2 數(shù)值水槽的建立

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

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

(8)

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

(9)

式中：k為波數(shù)；d為水深。

圖3 水槽中波腹點(diǎn)(x=10 m)和波節(jié)點(diǎn)(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)

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

2.2 明基床上透空沉箱參數(shù)設(shè)置

數(shù)值水槽長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所示。透空式沉箱的幾何尺寸和壓力測點(diǎn)的位置見圖5，沉箱開孔率ε=0.3。

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

3 數(shù)值計算結(jié)果及分析

3.1 測點(diǎn)壓力過程分析

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

圖6 1#、3#、4#、7#和12#測點(diǎn)的數(shù)值模擬結(jié)果于實驗結(jié)果的比較Fig.6 Comparison of simulated results and experimental results at sensor 1#, 3#, 4#, 7# and 12#

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

3.2 反射系數(shù)的影響因素分析

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

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

圖11 SPH計算結(jié)果與經(jīng)驗公式的反射系數(shù)比較Fig.11 Comparison of reflection coefficient of results from present SPH method and empirical formula

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

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

(10)

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

4 結(jié)論

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