李效民　張　林　郭海燕　姜　?！⊥酢?/p>

(中國海洋大學(xué)工程學(xué)院　青島　266100)

經(jīng)過長期的開采,陸上油氣資源日益減少,人們把目光轉(zhuǎn)向了資源豐富的海洋; 而復(fù)雜的海洋環(huán)境存在許多制約海洋開發(fā)的因素,內(nèi)波即是其中之一。內(nèi)波會對海洋立管、海洋平臺等海洋結(jié)構(gòu)物產(chǎn)生較強(qiáng)的破壞作用,給海洋油氣的開發(fā)帶來極大的挑戰(zhàn)。

海洋內(nèi)波是發(fā)生在密度穩(wěn)定層化的海水內(nèi)部的一種波動。海洋內(nèi)波的頻率和波速遠(yuǎn)小于表面波,但波長和振幅遠(yuǎn)大于表面波。據(jù)記載,最大垂向振幅高達(dá)180m(方欣華等,2005)。內(nèi)孤立波是一種強(qiáng)非線性潮成內(nèi)波,主要通過非線性效應(yīng)與色散效應(yīng)達(dá)到一定程度的平衡,實現(xiàn)穩(wěn)定傳播(Osborneet al,1980)。

目前廣泛使用的內(nèi)孤立波理論模型有 KdV理論、eKdV理論、mKdV理論、MCC理論等。研究內(nèi)孤立波的手段有理論分析、實地觀測及實測圖像分析、物理實驗、數(shù)值模擬,其中,物理實驗和數(shù)值模擬實施簡便、結(jié)果直觀可靠、適用性較強(qiáng),深受研究者的青睞,而數(shù)值模擬更是憑借其成本低、可重復(fù)性好、快速的優(yōu)勢成為最常用的內(nèi)孤立波研究方法。Wei等(2009)采用質(zhì)量源數(shù)值造波方法模擬了兩層流體中內(nèi)孤立波的生成與傳播。付東明等(2009)結(jié)合KdV和mKdV理論發(fā)展了雙推板數(shù)值造波方法。陳鈺等(2009)利用 Fluent軟件,采用速度入射邊界法,成功建立了模擬弱非線性內(nèi)孤立波的分層數(shù)值水槽。王旭等(2014)結(jié)合內(nèi)孤立波理論進(jìn)一步研究了質(zhì)量源數(shù)值造波方法,與理論結(jié)果較符合。Zhang等(2012)使用 CFD方法建立了三維數(shù)值水槽,所得結(jié)果與KdV理論結(jié)果吻合較好。高原雪等(2012)基于 MCC理論,采用速度入射邊界法建立數(shù)值水槽,實現(xiàn)了振幅可控的內(nèi)孤立波數(shù)值模擬。目前國內(nèi)外主要將數(shù)值模擬結(jié)果與理論結(jié)果進(jìn)行對比分析,而將理論、實驗、數(shù)值模擬三者進(jìn)行對比分析的研究較少。此外,對各種數(shù)值造波方法優(yōu)劣的比較也較為缺乏。

本文針對三種內(nèi)孤立波理論模型與數(shù)值造波方法,基于KdV、mKdV和eKdV三種理論模型與雙推板、平板拍擊和速度入射邊界三種數(shù)值造波方法所組成的一系列工況,利用Fluent軟件進(jìn)行計算,并將數(shù)值模擬結(jié)果與理論公式及實驗結(jié)果進(jìn)行對比分析,驗證了數(shù)值造波的合理性,評估了這三種數(shù)值造波方法。

1　內(nèi)孤立波理論模型

采用有限水深兩層流體簡化模型,假設(shè)上下層流體為無旋且不可壓縮的理想流體,上層流體密度為ρ1,深度為h1,下層流體密度為ρ2,深度為h2,總水深為h=h1+h2。如圖1,建立直角坐標(biāo)系oxyz。其中oxy平面為未擾動的兩層流體界面,ox軸正向為內(nèi)孤立波傳播方向,oz軸垂直于oxy面向上。令η(x,t)表示內(nèi)孤立波的波面函數(shù);為振幅大小,符號滿足上凸波為正,下凹波為負(fù);λ為內(nèi)孤立波特征波長;c為內(nèi)孤立波相速度。

圖1　內(nèi)孤立波及其參數(shù)Fig.1　Schema of internal solitary wave and the parameterization

1.1　KdV理論

波面方程的KdV理論解為(Choiet al,1999):

其中:

式中cKdV為 KdV理論下內(nèi)孤立波的相速度;lKdV為KdV理論下內(nèi)孤立波的特征長度。

KdV理論不存在極限振幅。但實驗研究發(fā)現(xiàn),振幅與水深比小于0.1時,KdV理論完全適用,而對于大振幅內(nèi)孤立波不再適用。

1.2　eKdV理論

波面方程的eKdV理論解為(Helfrichet al,2006):

其中:

式中ceKdV為eKdV理論下內(nèi)孤立波的相速度。

eKdV理論適用于振幅與水深比超過0.1的廣泛振幅范圍,且存在極限振幅,其大小為:

1.3　mKdV理論

波面方程的 mKdV理論解為(Michalletet al,1998):

其中:

式中cmKdV為mKdV理論下內(nèi)孤立波的相速度;為極限振幅大小。

mKdV理論適用于振幅與水深比超過0.1且上下層水深比接近臨界水深比的情況。

2　數(shù)值造波方法

目前,數(shù)值造波方法可大致歸為兩類: 一類是仿物理造波法,該方法源于實驗室造波技術(shù),即仿造物理造波機(jī)的工作原理。如: 推板法、平板拍擊法等。另一類是純數(shù)值造波法,如: 速度入射邊界法等。

2.1　雙推板法

圖2為雙推板法的示意圖,水槽上、下邊界及右邊界設(shè)置為固壁邊界,左側(cè)上、下推板設(shè)置為移動邊界。其造波原理為: 選取合適的內(nèi)孤立波理論,采用動網(wǎng)格技術(shù),通過用戶自定義函數(shù)(UDF)的DEFINE_CG_MOTION宏控制兩塊推板的水平反向運(yùn)動實現(xiàn)兩層流體界面處內(nèi)孤立波的模擬。

圖2　雙推板法造波示意圖Fig.2　Wave generation by double push-pedals method

令上板高度為d1,下板高度為d2,上下板的水平運(yùn)動速度分別為:U1、U2。根據(jù)上推板運(yùn)動排開的流體體積與兩層流體界面處水位在水平方向的積分相等,即:

可得上推板運(yùn)動速度為(徐鑫哲,2012):

再根據(jù)流體質(zhì)量守恒(徐鑫哲,2012),得下推板運(yùn)動速度為:

式中c為內(nèi)孤立波的相速度,η(x,t)為內(nèi)孤立波波面函數(shù)(上凸波為正,下凹波為負(fù))。

實際設(shè)置時,應(yīng)使d1略大于以防上下層流體摻混(付東明等,2009)。

2.2　平板拍擊法

圖3為平板拍擊法的示意圖,水槽上、下、左、右邊界及隔板均設(shè)置為固壁邊界,造波板設(shè)置為移動邊界。其造波原理為: 選取合適的內(nèi)孤立波理論,采用動網(wǎng)格技術(shù),通過用戶自定義函數(shù)(UDF)的DEFINE_CG_MOTION宏控制造波板做上下運(yùn)動實現(xiàn)兩層流體界面處內(nèi)孤立波的模擬。

圖3　平板拍擊法造波示意圖Fig.3　Wave generation by moving pedal method

令造波板的長度為D,運(yùn)動速度為U。根據(jù)造波板移動所排開的流體體積與造波板右端兩層流體界面處水位在水平方向的積分相等,即:

可得造波板運(yùn)動速度為(徐鑫哲,2012):

2.3　速度入射邊界法

圖4為速度入射邊界法的造波示意圖,水槽上、下邊界及右邊界設(shè)置為固壁邊界,左側(cè)邊界設(shè)置為速度入射邊界。該方法是以邊界為擾動源,但邊界是固定不運(yùn)動的。其造波原理為: 選取合適的內(nèi)孤立波理論,通過用戶自定義函數(shù)(UDF)的 DEFINE_PROFILE宏在入射邊界上給定該理論下的水流速度實現(xiàn)兩層流體界面處內(nèi)孤立波的模擬。

圖4　速度入射邊界法造波示意圖Fig.4　Wave generation by velocity-inlet boundary method

入射邊界處水流速度設(shè)為:

3　數(shù)值模擬與結(jié)果分析

3.1　內(nèi)孤立波物理實驗

王飛(2015)在中國海洋大學(xué)物理海洋實驗室分層流體水槽中進(jìn)行了內(nèi)孤立波物理實驗。水槽長15m,寬0.35m,高0.7m,額定水深0.6m。采用Oster雙缸法制取密度分層流體,上層流體高度 9.5cm,密度1035kg/m3,下層流體高度48.5cm,密度1054kg/m3。采用重力塌陷法制造不同振幅的內(nèi)孤立波,使用兩個斜面在水槽末端進(jìn)行消波,利用PIV技術(shù)對內(nèi)孤立波致流場進(jìn)行測量,使用高速相機(jī)進(jìn)行圖像采集,利用圖像處理軟件進(jìn)行數(shù)據(jù)分析。物理實驗結(jié)果如圖6、圖8所示。

3.2　建立數(shù)值水槽

為便于與實驗結(jié)果對比,本文參照內(nèi)波實驗建立數(shù)值水槽模型。采用兩層流體簡化模型,利用Fluent進(jìn)行內(nèi)孤立波數(shù)值模擬。采用結(jié)構(gòu)化網(wǎng)格離散,在兩層流體交界面附近對網(wǎng)格劃分進(jìn)行加密。工作區(qū)網(wǎng)格劃分如圖5(a)所示: 底部向上38cm至50cm高度范圍內(nèi),網(wǎng)格尺寸為0.5cm,其余為1cm,x軸方向網(wǎng)格尺寸為2cm。選用二維不可壓縮Navier-Stokes方程作為流體控制方程,使用有限體積法對控制方程進(jìn)行離散。選用κ-ε湍流模型,采用VOF方法追蹤兩層流體的交界面,兩層流體界面的構(gòu)造采用幾何重構(gòu)法,壓力速度耦合選用 PISO算法,壓力插值采用體積力加權(quán)法,對流相及輸運(yùn)方程的離散采用一階迎風(fēng)格式。

由于水槽的右邊界為固壁邊界,易產(chǎn)生反射波,本文采用阻尼消波法在水槽末端進(jìn)行消波。消波段長度取為 1—2倍波長。同時,在消波區(qū)沿x軸方向采用逐漸增大的漸變網(wǎng)格來耗散波浪的能量。消波區(qū)網(wǎng)格劃分如圖 5(b)所示:z軸方向網(wǎng)格劃分與工作區(qū)相同,x軸方向首層網(wǎng)格尺寸為 2cm,后續(xù)網(wǎng)格尺寸按1.04的比例逐漸變大。

圖5　數(shù)值水槽網(wǎng)格劃分示意圖Fig.5　Grids of the numerical flume

根據(jù)內(nèi)孤立波理論對振幅水深比的使用條件(黃文昊等,2013),按照不同的振幅,選取合適的理論模型,采用不同的數(shù)值造波方法進(jìn)行數(shù)值模擬。具體工況如表1所示。

設(shè)置各項參數(shù)后進(jìn)行流場的迭代計算,迭代時間步長取 0.01s,每個時間步長的最大迭代次數(shù)取 20次。數(shù)值計算結(jié)果如圖6、圖8所示。

3.3　計算結(jié)果與分析

圖 6為各工況下數(shù)值模擬所得波形與相應(yīng)理論及實驗結(jié)果的對比。結(jié)果表明: 三種數(shù)值造波方法均能生成內(nèi)孤立波,但模擬效果有差異。由工況 A1、A2、A3(B1、B2、B3 和 C1、C2、C3)可知,對于同一種振幅、同一種內(nèi)孤立波理論,速度入射邊界法所模擬波形與理論及實驗波形吻合最好。平板拍擊法所模擬波形與理論及實驗波形吻合較好,只是在波形尾端存在較小的尾波,但在傳播過程中尾波會逐漸衰減,不會對內(nèi)波產(chǎn)生影響。雙推板法所得模擬波形與理論及實驗波形趨勢相同,但最大振幅達(dá)不到設(shè)計振幅。

表1　數(shù)值模擬工況設(shè)置Tab.1　Condition setting of numerical simulation

由工況 A1、B1、C1(A2、B2、C2和 A3、B3、C3)可知,隨著振幅的增加,速度入射邊界法和平板拍擊法所得波形與理論及實驗波形吻合仍然較好。雙推板所得波形與理論及實驗波形趨勢仍然一致,但最大振幅與設(shè)計振幅的差值越來越大,其原因為流體阻礙了推板的運(yùn)動造成耗散。由圖7(a)可知A1工況下加大推板的運(yùn)動幅值后所得波形與理論及實驗波形吻合較好。因此,可通過加大推板的運(yùn)動幅值使所生成內(nèi)孤立波的幅值滿足要求。此外,雙推板法數(shù)值造波所得幅值還受上下層流體深度比的影響。實驗研究發(fā)現(xiàn): 實測振幅均比設(shè)計振幅小,且兩者之間存在線性關(guān)系(黃文昊等,2013),因此,實際使用時需擬合相同條件下不同設(shè)計振幅與所得振幅之間的關(guān)系,然后根據(jù)實際所需振幅反算出設(shè)計振幅。

圖 6顯示三種數(shù)值造波方法結(jié)果與實驗結(jié)果存在差異,尤其隨著波高的增大,兩者的波形在末端差異越大。其原因有: (1)內(nèi)波在傳播過程中粘滯效應(yīng)影響會導(dǎo)致能量損耗。(2)實驗造波方法為重力塌陷法,波高與塌陷高度差呈正相關(guān),塌陷高度差越大,對流體產(chǎn)生的激蕩也越強(qiáng),引發(fā)的尾波越明顯。(3)實驗的測量設(shè)備、邊界條件等的限制導(dǎo)致實驗結(jié)果與數(shù)值模擬結(jié)果不能完全一致。

圖6　數(shù)值造波波形與理論及實驗結(jié)果對比Fig.6　Comparison in numerical waveform between experimental and theoretical results

圖7　A1工況下修正后的數(shù)值造波結(jié)果與理論及實驗結(jié)果的對比Fig.7　Comparison of calibrated numerical result with experimental and theoretical ones for Case A1

由工況B1、B2、B3可知,mKdV理論波形明顯要寬于實驗波形與數(shù)值模擬波形。這一偏差說明mKdV理論對這種波高水深比情況下的內(nèi)孤立波的描述存在一定偏差,但基于該理論的數(shù)值造波仍然能給出比較接近實際情況的波形。

圖 8為各工況下數(shù)值模擬所得波谷經(jīng)過斷面處波致水平流速沿垂向的分布,并與理論及實驗結(jié)果進(jìn)行了對比。結(jié)果表明: 三種數(shù)值造波方法所得波致水平流速沿垂向分布趨勢基本一致,但有所差異。由工況 A1、A2、A3(B1、B2、B3)可知,平板拍擊法和速度入射邊界法所得波致水平流速均與理論及實驗結(jié)果吻合較好。雙推板法所得波致水平流速相對于理論及實驗結(jié)果明顯偏小。而對于工況C1、C2、C3,由于實驗時上層流體中波致水平流速較大,致使實驗監(jiān)測數(shù)據(jù)失效,實驗結(jié)果失真。但仍然可以看出速度入射邊界法和平板拍擊法所得波致水平流速與理論結(jié)果吻合較好,雙推板法與理論結(jié)果存在一定差異。

圖8　內(nèi)孤立波波致水平流速沿垂向分布Fig.8　Horizontal velocity magnitude induced by internal solitary waves along vertical section

由工況 A1、B1、C1(A2、B2、C2和 A3、B3、C3)可知,隨著振幅的增加,速度入射邊界法和平板拍擊法所得結(jié)果與理論及實驗結(jié)果吻合仍然較好。雙推板法所得結(jié)果與理論及實驗結(jié)果趨勢相同,但上下層流體水平速度與理論及實驗結(jié)果的差值越來越大。加大推板的運(yùn)動幅值,所得內(nèi)孤立波的上下層流體水平速度與理論及實驗結(jié)果吻合較好,如圖7(b)所示。

圖8顯示實驗過程中通過PIV測量得到的上層流體水平流速變化情況比較混亂,沒有規(guī)律,與數(shù)值模擬結(jié)果存在差異,其原因為上層流體中波致水平流速較大,超出了PIV系統(tǒng)的測量范圍,故測得的數(shù)據(jù)存在一定失真。

內(nèi)界面與波面之間區(qū)域的誤差主要原因是實驗與數(shù)值模擬存在密度躍層,密度躍層內(nèi)會有流速變化,而理論模型是簡單地分為兩層處理,沒有中間過渡階段。

綜合上述分析,可以發(fā)現(xiàn)三種造波方法均能生成內(nèi)孤立波,但雙推板法需擬合設(shè)計振幅與所需振幅之間的關(guān)系,通過加大推板的運(yùn)動幅值來生成所需振幅的內(nèi)孤立波,無法直接有效實現(xiàn)對內(nèi)孤立波的數(shù)值模擬。而其余兩種數(shù)值造波方法均能直接有效模擬內(nèi)孤立波的生成與傳播,且所得結(jié)果與理論及實驗結(jié)果吻合較好。

表2為三種數(shù)值造波方法模擬效率的對比。從中可以看出: 在相同電腦配置(windows 7 64位,16G內(nèi)存,i7處理器)、相同振幅、相同理論模型、相同模擬時間的情況下,平板拍擊法與雙推板法所用時間相當(dāng),速度入射邊界法所用時間明顯少于前面兩種方法。其原因是仿物理造波方法使用動網(wǎng)格,需占用更多的計算資源。

表2　三種數(shù)值造波方法模擬效率的比較Tab.2　Efficiency of three numerical wave-generating methods

4　結(jié)論

利用 Fluent軟件計算了基于 KdV、mKdV和eKdV三種理論模型與雙推板、平板拍擊和速度入口三種數(shù)值造波方法所組成的九種工況,并將模擬結(jié)果與理論及實驗結(jié)果進(jìn)行對比和分析,得到如下結(jié)論:

(1) 三種數(shù)值造波方法均能實現(xiàn)對內(nèi)孤立波的數(shù)值模擬,除雙推板法所得結(jié)果與理論及實驗結(jié)果存在一定差異,需通過加大推板的運(yùn)動振幅來實現(xiàn)所需的造波效果外,其余兩種數(shù)值造波方法均能直接有效模擬內(nèi)孤立波的生成與傳播,且所得結(jié)果與理論及實驗結(jié)果吻合較好。

(2) 雙推板法與平板拍擊法均屬于仿物理造波法,但平板拍擊法所得結(jié)果明顯優(yōu)于雙推板法,因此實際使用時應(yīng)優(yōu)先選用平板拍擊法。

(3) 同工況同配置條件下,速度入射邊界法模擬效率最高,所得結(jié)果更準(zhǔn)確,因此實際使用時應(yīng)優(yōu)先選用速度入射邊界法。

王飛,2015. 內(nèi)孤立波作用下小尺度豎直圓柱體的水動力特性研究. 青島: 中國海洋大學(xué)博士學(xué)位論文,21—32

王旭,林忠義,尤云祥,2014. 兩層流體中內(nèi)孤立波質(zhì)量源數(shù)值造波方法. 上海交通大學(xué)學(xué)報,48(6): 850—855

方欣華,杜濤,2005. 海洋內(nèi)波基礎(chǔ)和中國海內(nèi)波. 青島:中國海洋大學(xué)出版社,1—2

付東明,尤云祥,李巍,2009. 兩層流體中內(nèi)孤立波與潛體相互作用數(shù)值模擬. 海洋工程,27(3): 38—44

陳鈺,朱良生,2009. 基于Fluent的海洋內(nèi)孤立波數(shù)值水槽模擬. 海洋技術(shù),28(4): 72—75,100

高原雪,尤云祥,王旭等,2012. 基于MCC理論的內(nèi)孤立波數(shù)值模擬. 海洋工程,30(4): 29—36

徐鑫哲,2012. 內(nèi)波生成機(jī)理及二維內(nèi)波數(shù)值水槽模型研究.哈爾濱: 哈爾濱工程大學(xué)碩士學(xué)位論文,35—36

黃文昊,尤云祥,王旭等,2013. 有限深兩層流體中內(nèi)孤立波造波實驗及其理論模型. 物理學(xué)報,62(8): 084705

Choi W,Camassa R,1999. Fully nonlinear internal waves in a two-fluid system. Journal of Fluid Mechanics,396(3): 1—36

Helfrich K R,Melville W K,2006. Long nonlinear internal waves. Annual Review of Fluid Mechanics,38(1): 395—425

Michallet H,Barthélemy E,1998. Experimental study of interfacial solitary waves. Journal of Fluid Mechanics,366(1): 159—177

Osborne A R,Burch T L,1980. Internal solitons in the Andaman Sea. Science,208(4443): 451—460

Wei G,You Y X,Su X B,2009. Two-dimension numerical internal wave tank for navier-stokes equation model in the stratified fluid. In: Proceedings of the Fifth International Conference on Fluid Mechanics. Berlin Heidelberg: Springer,364—367

Zhang L,Wang L L,Yu Z Zet al,2012. Characteristics of non-linear internal waves in a three-dimensional numerical wave tank.Applied Mechanics and Materials,212—213: 1123—1130