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

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

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

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

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

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

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

0 引 言

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

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

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

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

對于主相，它的體積分數(shù)由（2）式約束計算得到。

體積分數(shù)方程通過顯式時間的離散來求解，采用當前時刻之前的時間步長帶入標準有限差分插值方案中進行體積分數(shù)值的計算。

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

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

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

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

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

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

圖1 計算模型示意圖Fig.1 Sketch of the calculation model

在計算中取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）。

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

其中：U為潛體運動速度，Δρ為兩層流體密度差，h′為上層流體厚度、h為下層流體厚度。

2 潛體勻速運動對內(nèi)波波高及波長的影響

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

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

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

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

圖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 潛體加/減速運動對內(nèi)波波幅的影響

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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

[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ù)研專項基金資助項目

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