丁 勇，段 菲，韓盼盼，牛明昌

（哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001）

兩層流中潛體運(yùn)動(dòng)與誘發(fā)內(nèi)波特征關(guān)系研究

丁 勇，段 菲，韓盼盼，牛明昌

（哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001）

在兩層流體中潛體運(yùn)動(dòng)會(huì)激發(fā)生成內(nèi)波。文章通過(guò)編程控制潛體運(yùn)動(dòng)，使用動(dòng)態(tài)分層法保證潛體運(yùn)動(dòng)時(shí)邊界層不受影響，基于RANS方程多相流模型，采用求解相的體積分?jǐn)?shù)方式捕捉內(nèi)界面波形。潛體定常運(yùn)動(dòng)計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果相符，通過(guò)此方法，對(duì)比不同潛深下潛體誘發(fā)內(nèi)波波高和波長(zhǎng)隨密度Froude數(shù)的變化特征，得到了誘發(fā)內(nèi)波波高最大處的臨界密度Froude數(shù)，重點(diǎn)研究了在臨界密度Froude數(shù)處潛體加速運(yùn)動(dòng)和減速運(yùn)動(dòng)對(duì)誘發(fā)內(nèi)波幅值及其波動(dòng)過(guò)程的影響。

兩層流體；潛體尾跡；非定常運(yùn)動(dòng)；內(nèi)波

0 引 言

水下航行體在分層流體中運(yùn)動(dòng)時(shí)，其周圍水質(zhì)點(diǎn)受到體積效應(yīng)下的浮力作用，在約化重力作用下，分層界面附近流場(chǎng)將產(chǎn)生波動(dòng)從而形成了內(nèi)波。在源生內(nèi)波的研究中，兩層流體模型受到人們的格外關(guān)注。在無(wú)粘假設(shè)下，Lu等[1-2]采用穩(wěn)定相方法給出了點(diǎn)源運(yùn)動(dòng)產(chǎn)生的內(nèi)界面波形解，Yeung[3]和Wei[4]等通過(guò)Green函數(shù)法探討了點(diǎn)源運(yùn)動(dòng)誘生的內(nèi)波模式與表面波模式對(duì)內(nèi)界面興波的影響，Xu等[5]運(yùn)用線性理論研究了含密越層的三層流體中卵形體產(chǎn)生的內(nèi)波速度梯度場(chǎng)，Zhu等[6]采用邊界元理論方法以求解球體在上層流體運(yùn)動(dòng)的內(nèi)波波形；Gourlay[7]以及Fritts[8]等采用直接模擬法研究了粘流場(chǎng)中漩渦對(duì)源生內(nèi)波的影響，Chernykh[9]采用RANS下的數(shù)值模擬方法模擬了分層流中有凈動(dòng)量的自航尾跡。目前絕大多數(shù)研究以潛體定常運(yùn)動(dòng)狀態(tài)為研究對(duì)象，以相對(duì)運(yùn)動(dòng)法來(lái)模擬潛體在分層流中的運(yùn)動(dòng)。

本文將通過(guò)編程控制潛體的絕對(duì)運(yùn)動(dòng)，運(yùn)用動(dòng)態(tài)分層法保證潛體運(yùn)動(dòng)時(shí)邊界層不受影響，確保此處有更穩(wěn)定的y+值以提高計(jì)算精度，基于RANS方程多相流模型，采用CFD技術(shù)，重點(diǎn)研究潛體加速和減速運(yùn)動(dòng)對(duì)內(nèi)波幅值及其波動(dòng)過(guò)程的影響。通過(guò)與實(shí)驗(yàn)值對(duì)比定常運(yùn)動(dòng)結(jié)果，確定本文模擬方法的可行性，對(duì)比分析不同潛深下潛體定常運(yùn)動(dòng)誘發(fā)內(nèi)波波高、波長(zhǎng)隨密度Froude數(shù)的變化特征，得到誘發(fā)內(nèi)波波高最大處的臨界密度Froude數(shù)；在臨界密度Froude數(shù)處的運(yùn)動(dòng)速度下，定義加速度無(wú)量綱數(shù)，研究潛體加速和減速運(yùn)動(dòng)對(duì)內(nèi)波幅值及其波動(dòng)過(guò)程的影響，進(jìn)一步定量討論加速度和減速度的大小與內(nèi)波幅值的關(guān)系。

1 數(shù)值計(jì)算方法與模型

在兩層流體中求解非穩(wěn)態(tài)的連續(xù)方程和Navier-Stokes方程時(shí)，假設(shè)兩種液體互不相溶，以q流體為例，引入體積分?jǐn)?shù)變量αq，兩層流體交界面通過(guò)求解q流體相的體積分?jǐn)?shù)連續(xù)性方程（1）得到：

對(duì)于主相，它的體積分?jǐn)?shù)由（2）式約束計(jì)算得到。

體積分?jǐn)?shù)方程通過(guò)顯式時(shí)間的離散來(lái)求解，采用當(dāng)前時(shí)刻之前的時(shí)間步長(zhǎng)帶入標(biāo)準(zhǔn)有限差分插值方案中進(jìn)行體積分?jǐn)?shù)值的計(jì)算。

其中：n+1代表當(dāng)前的時(shí)間步長(zhǎng)；n代表前一間隔的時(shí)間步長(zhǎng)；αq，f是通過(guò)一階迎風(fēng)格式計(jì)算得到的q流體體積分?jǐn)?shù)在表面處的值；V為計(jì)算單元的體積；Uf是根據(jù)正常速度流體在表面處的體積通量。

由于表面力使用散度定理可以表示為體積力，而這個(gè)體積力正是添加給動(dòng)量方程的源項(xiàng)，此模型采用Brackbill等[10]提出的連續(xù)表面力模型。當(dāng)一個(gè)單元中只有兩相時(shí)，體積力為：

其中：σij為表面張力系數(shù)，為表面曲率，。

本文通過(guò)編程控制潛體運(yùn)動(dòng)規(guī)律，運(yùn)用動(dòng)態(tài)分層法保證潛體運(yùn)動(dòng)時(shí)邊界層不受影響，確保此處有更穩(wěn)定的y+值，通過(guò)理想高度和分割因子（合并因子）兩種參數(shù)來(lái)控制這一過(guò)程。（5）式表示當(dāng)網(wǎng)格被拉伸大于一定高度時(shí)，該層網(wǎng)格會(huì)分裂為兩層網(wǎng)格，（6）式表示當(dāng)網(wǎng)格被壓縮小于一定高度時(shí)，該層網(wǎng)格會(huì)與相鄰的另一層網(wǎng)格相融合變?yōu)橐粚泳W(wǎng)格。

其中：h是網(wǎng)格高度，hideal是理想高度，αs是分割因子，αc是合并因子。

為了使用這種方法，需要對(duì)計(jì)算域進(jìn)行變形。在潛體運(yùn)動(dòng)方向分割出一個(gè)條形內(nèi)域，內(nèi)域隨著潛體的水平運(yùn)動(dòng)進(jìn)行平動(dòng)，其他部分為外域，內(nèi)域與外域之間以interface相連，這樣內(nèi)域在平動(dòng)時(shí)會(huì)引起周圍網(wǎng)格的生成與消去，但不會(huì)對(duì)外域的網(wǎng)格造成影響。圖1為計(jì)算模型示意圖。

圖1 計(jì)算模型示意圖Fig.1 Sketch of the calculation model

在計(jì)算中取h′=0.15 m，h=0.25 m，潛體直徑D=0.05 m，ρ1= 998.2 kg/m3，μ1=1.003×10-3kg/（m·s ），ρ2=1 023.79 kg/m3，μ2= 0.981 1×10-3kg/（m·s）。

在密度分層的流體中潛體運(yùn)動(dòng)會(huì)誘發(fā)生成內(nèi)波，內(nèi)波特征取決于密度Froude數(shù)Frd。

其中：U為潛體運(yùn)動(dòng)速度，Δρ為兩層流體密度差，h′為上層流體厚度、h為下層流體厚度。

2 潛體勻速運(yùn)動(dòng)對(duì)內(nèi)波波高及波長(zhǎng)的影響

潛體勻速運(yùn)動(dòng)誘發(fā)內(nèi)波特征計(jì)算結(jié)果將與Arntsen[11]在分層水槽中的實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。圖2所示即為本文方法計(jì)算得到的潛體位于內(nèi)界面之上1.6D、Frd=0.6時(shí)的波形圖，內(nèi)波波高值取前兩個(gè)波中的最大值，波長(zhǎng)值取最大波高處對(duì)應(yīng)的波長(zhǎng)。本文重點(diǎn)考察潛體運(yùn)動(dòng)速度以及潛深的變化對(duì)內(nèi)波波形的影響，因此密度Froude數(shù)中的變量?jī)H改變運(yùn)動(dòng)速度，考察范圍為0.2＜Frd＜1.0，潛深的影響分別考察潛體位于內(nèi)界面之上2.3D、1.6D、1.0D以及位于內(nèi)界面之下0.9D、1.4D的情況。圖3為內(nèi)波波高隨Frd變化的對(duì)比結(jié)果，波長(zhǎng)變化結(jié)果見圖4。

圖2 Z0/D=1.6，F(xiàn)rd=0.6時(shí)內(nèi)波波形圖Fig.2 Internal wave at Z0/D=1.6,Frd=0.6

在本文考察的所有計(jì)算范圍內(nèi)，數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合良好，確定通過(guò)編程和動(dòng)態(tài)分層法模擬分層流中潛體運(yùn)動(dòng)的可行性。結(jié)果顯示，內(nèi)波波高與波長(zhǎng)隨密度Froude數(shù)的變化趨勢(shì)與潛深無(wú)關(guān)。內(nèi)波波高值隨運(yùn)動(dòng)速度的增加先增大后減小，在Frd≈0.6波高最大，在0.2＜Frd＜1.0的范圍內(nèi)，波高值減小幅度對(duì)比增大幅度較??；潛體潛深僅改變幅值大小，潛體運(yùn)動(dòng)距離內(nèi)界面越近，誘發(fā)內(nèi)波幅值越大，而當(dāng)Frd＞0.8時(shí)，潛體位于內(nèi)界面之下距離1.4D誘發(fā)的內(nèi)波幅值比0.9D的幅值大。

內(nèi)波波長(zhǎng)的度量以上層流體高度h′為參考，隨著潛體運(yùn)動(dòng)速度的增加波長(zhǎng)值呈指數(shù)形式增大，在0.2＜Frd＜1.0的的范圍內(nèi)，最大波長(zhǎng)值發(fā)生在Frd=1.0處，當(dāng)潛體在上層流體中運(yùn)動(dòng)時(shí)，潛體距離內(nèi)界面越近，相同速度下誘發(fā)生成的內(nèi)波波長(zhǎng)越小，而當(dāng)潛體在下層流體中運(yùn)動(dòng)時(shí)，結(jié)果相反；與實(shí)驗(yàn)結(jié)果相比，本文由于假設(shè)兩種液體沒(méi)有相互穿插，密越層厚度幾乎為0，而實(shí)驗(yàn)中由于淡水與鹽水的相容性，存在一定厚度的密越層，約為圓柱直徑的0.7～1.0倍，因此與實(shí)驗(yàn)結(jié)果存在一定偏差，但是本文計(jì)算方法有效地解決了由于水池長(zhǎng)度的限制Frd＞0.8時(shí)產(chǎn)生的內(nèi)波未達(dá)到有效的穩(wěn)定狀態(tài)，波長(zhǎng)無(wú)法獲得的缺陷，以及內(nèi)波波高太小難以測(cè)量等問(wèn)題。

圖3 不同潛深下內(nèi)波波高（H/D）隨Frd變化Fig.3 Variation in internal wave heights(H/D)with Frdat different depth

圖4 不同潛深下內(nèi)波波高（λ/h′）隨Frd變化Fig.4 Variation in internal wave heights(λ/h′)with Frdat different depth

3 潛體加/減速運(yùn)動(dòng)對(duì)內(nèi)波波幅的影響

本文重點(diǎn)研究潛體非定常運(yùn)動(dòng)時(shí)誘發(fā)內(nèi)波幅值與波動(dòng)過(guò)程特征，定義描述加速度無(wú)量綱數(shù)為，其中a為潛體運(yùn)動(dòng)加速度，υ=μ/ρ為流體運(yùn)動(dòng)粘性系數(shù)。潛體勻速運(yùn)動(dòng)結(jié)果顯示，在各個(gè)潛深處當(dāng)Frd≈0.6時(shí)內(nèi)波波高最大，因此當(dāng)潛體位于內(nèi)界面之上2.3D處以Frd=0.6勻速運(yùn)動(dòng)并達(dá)到穩(wěn)定狀態(tài)時(shí)，分別研究潛體以Ara=±0.908 5×103、Ara=±1.112 7×103、Ara=±1.573 6×103作加速運(yùn)動(dòng)和減速運(yùn)動(dòng)誘發(fā)內(nèi)波幅值與波動(dòng)過(guò)程。

定義在運(yùn)動(dòng)潛體正下方誘發(fā)的內(nèi)波第一波峰稱之為“First Hump”，內(nèi)波第一波谷稱之為“First Valley”也稱之為首波峰及首波谷，第二個(gè)波峰稱之為“Second Hump”，波谷稱之為“Second Valley”也稱之為次波峰及次波谷。

圖5 潛體減速運(yùn)動(dòng)誘發(fā)生成內(nèi)波的峰值變化Fig.5 Peak-curves of internal wave induced by submerged body deceleration motion

圖6 潛體做加速運(yùn)動(dòng)誘發(fā)生成內(nèi)波的峰值變化Fig.6 Peak-curves of internal wave induced by submerged body acceleration motion

圖5為不同減速度下潛體減速運(yùn)動(dòng)過(guò)程中內(nèi)波峰值變化曲線，每隔0.5 s取一次值，可以看出內(nèi)波峰值都是呈波動(dòng)性的變化。處于潛體正下方的首波峰峰值隨潛體速度減小H/D從0.41迅速衰減，當(dāng)潛體速度減為0時(shí)，以減速度Ara=-1.573 6×103運(yùn)動(dòng)導(dǎo)致的內(nèi)波首波峰值H/D減小到約0.12，顯然此時(shí)潛體雖靜止但首波峰還未完全衰減，而當(dāng)減速度為Ara=-0.908 5×103和Ara=-1.112 7×103時(shí)，速度減小接近0時(shí)首波峰幾乎為0；次波峰的變化伴隨著減速的開始逐漸變大，H/D從0.3變化到0.6，波動(dòng)幅度幾乎增大了一倍，以增大后的幅度波動(dòng)性地持續(xù)一段時(shí)間后迅速衰減，從圖5（b）可看出幾乎持續(xù)了三個(gè)波動(dòng)周期，對(duì)比不同的減速度，當(dāng)潛體速度減為0時(shí)，次波峰峰值隨減速度大小依次為0.55、0.35和0.20，減速度為Ara=-1.573 6×103時(shí)，次波峰持續(xù)增大，并未衰減。

圖6為不同加速度下潛體加速運(yùn)動(dòng)過(guò)程中內(nèi)波峰值變化曲線，當(dāng)Frd=1.0時(shí)完成加速運(yùn)動(dòng)，加速運(yùn)動(dòng)誘發(fā)生成的內(nèi)波峰谷值依舊呈波動(dòng)性變化，首波谷隨潛體運(yùn)動(dòng)速度的增加H/D最大從0.05增大到0.21，但是相應(yīng)的加速度最小為Ara=0.908 5×103，以最大加速度運(yùn)動(dòng)的潛體在完成加速運(yùn)動(dòng)時(shí)H/D僅從0.05增大到約0.10，首波峰與次波谷的變化趨勢(shì)與首波谷截然相反，均隨運(yùn)動(dòng)速度增加而減小，在加速度為Ara=1.573 6×103時(shí)潛體運(yùn)動(dòng)對(duì)首波峰影響不大，H/D僅從0.37減小到0.34，次波谷谷值從0.48減小到0.38，在最小加速度下，首波峰與次波谷衰減最明顯，加速完成時(shí)分別減小到0.23和0.33。

4 結(jié) 論

本文基于RANS方程多相流模型，采用CFD技術(shù)，通過(guò)編程和動(dòng)態(tài)分層法模擬潛體真實(shí)運(yùn)動(dòng)，確保在模擬過(guò)程中邊界層不受影響，通過(guò)與實(shí)驗(yàn)值對(duì)比定常運(yùn)動(dòng)計(jì)算結(jié)果，確定本文模擬方法的可行性。對(duì)比分析不同潛深下潛體定常運(yùn)動(dòng)誘發(fā)內(nèi)波波高、波長(zhǎng)隨密度Froude數(shù)的變化特征，進(jìn)一步討論了潛體加速和減速運(yùn)動(dòng)對(duì)內(nèi)波幅值及其波動(dòng)過(guò)程的影響，具體結(jié)論如下：

（1）本文計(jì)算方法可有效地模擬潛體運(yùn)動(dòng)誘發(fā)內(nèi)波實(shí)驗(yàn)，并且解決了由于水池長(zhǎng)度的限制當(dāng)潛體運(yùn)動(dòng)速度過(guò)大內(nèi)波波長(zhǎng)無(wú)法獲得的缺陷，以及內(nèi)波波高太小難以測(cè)量等問(wèn)題。

（2）潛體勻速運(yùn)動(dòng)誘發(fā)內(nèi)波結(jié)果顯示，在0.2＜Frd＜1.0的范圍內(nèi)，內(nèi)波波高值先增大后減小，在Frd≈0.6時(shí)波高最大，不同潛深下潛體距離內(nèi)界面越近，內(nèi)波波高越大；潛體勻速運(yùn)動(dòng)誘發(fā)內(nèi)波波長(zhǎng)呈指數(shù)形式增大，潛體在上層流體中運(yùn)動(dòng)時(shí)，距離內(nèi)界面越近，誘發(fā)生成的內(nèi)波波長(zhǎng)越小，而當(dāng)潛體在下層流體中運(yùn)動(dòng)時(shí)，結(jié)果相反。

（3）潛體加速運(yùn)動(dòng)和減速運(yùn)動(dòng)過(guò)程中內(nèi)波呈現(xiàn)波動(dòng)性變化，潛體正下方的內(nèi)波波幅變化趨勢(shì)與速度變化相同，即加速運(yùn)動(dòng)時(shí)首波峰（谷）值增大，減速運(yùn)動(dòng)時(shí)首波峰（谷）值減小，此外，潛體加/減速度越小，運(yùn)動(dòng)完成時(shí)峰谷幅值變化越大。

（4）當(dāng)潛體從誘發(fā)內(nèi)波最大的速度處（即Frd≈0.6）開始減速時(shí)，次波峰先增大后減小，加速時(shí)，首波峰以及次波谷均呈波動(dòng)性減小。

[1]Lu D Q,Chen T.Surface and interfacial gravity waves induced by an impulsive disturbance in a two-layer inviscid fluid [J].Journal of Hydrodynamics,Ser.B,2009,21(1)：26-33.

[2]Lu D Q,Chwang A T.Interfacial waves due to a singularity in a system of two semi-infinite fluids[J].Physics of Fluids (1994-present),2005,17(10)：102-107.

[3]Yeung R W N.Waves generated by a moving source in a two-layer ocean of finite depth[J].Journal of Engineering Mathematics,1999,35(1)：85-107.

[4]Gang Wei,Jiachun L,Shiqiang D.Surface effects of internal wave generated by a moving source in a two-layer fluid of finite depth[J].Applied Mathematics and Mechanics,2003,24(9)：1025-1040.

[5]Xu Z,Chen X,Sturova I V.Surface characters of internal waves generated by Rankine ovoid[J].Acta Mechanica Sinica, 2006,22(5)：417-423.

[6]Zhu W,You Y X,Miao G,et al.Waves generated by a 3 D moving body in a two-layer fluid of finite depth[J].Journal of Hydrodynamics,Series B,2005,17(1)：92-101.

[7]Gourlay M J W J F D.Direct numerical simulations of initially turbulent wakes[Z].Biloxi,MS：200118.

[8]Fritts D,Gourlay M,Orlando W,et al.Numerical simulation of late wakes in stratified and sheared flows[Z].2003：206-210.

[9]Chernykh G G,Moshkin N P,Fomina A V.Dynamics of turbulent wake with small excess momentum in stratified media [J].Communications in Nonlinear Science and Numerical Simulation,2009,14(4)：1307-1323.

[10]Brackbill J U,D.Kothe B,Zemach C.A continuum method for modeling surface tension[J].J Comput.Phys.,1992,100：335-354.

[11]Arntsen ? A.Disturbances,lift and drag forces due to the translation of a horizontal circular cylinder in stratified water [J].Experiments in Fluids,1996,21(5)：387-400.

Research on the relationship between moving patterns of submerged body and the features of induced internal waves in two layer fluid

DING Yong,DUAN Fei,HAN Pan-pan,NIU Ming-chang
(College of Shipbuilding Engineering,Harbin Engineering University,Harbin 150001,China)

The motion of submerged body in stratified fluids will excite the generation of internal waves.In this paper,by control the motion of submerged body using user defined files,using layering to ensure the boundary layer will not be influenced by the motion of submerged body.The waveform of interface was captured by means of solving the volume fraction of each fluids based on multiphase flow model of RANS equation.The calculation results of uniform motion were in good agreement with experimental results.Comparing the the characteristics of variation in internal wave heights and lengths with densimetric Froude number at different depth the critical densimetric Froude number induced the maximum internal wave heights was obtained by this method.The interaction of submerged body acceleration and deceleration motion with induced internal wave at critical densimetric Froude number was mainly studied.

two layer fluid;submerged body wake;unsteady motion;internal wave

U661.1

：Adoi：10.3969/j.issn.1007-7294.2016.05.002

1007-7294（2016）05-0523-07

2016-01-03

裝備預(yù)研專項(xiàng)基金資助項(xiàng)目

丁 勇（1959-），男，博士，教授，E-mail：dingyong@hrbeu.edu.cn；段 菲（1990-），女，碩士研究生。