旋轉(zhuǎn)圓球入水空泡特性與流場結(jié)構(gòu)的大渦模擬研究1)

嚴(yán)晨祎 陳 瑛

(上海交通大學(xué)船舶海洋與建筑工程學(xué)院,水動(dòng)力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室,上海 200240)

引言

物體穿越水面的入水問題與許多軍事應(yīng)用直接相關(guān),如:艦載發(fā)射航行體入水[1]、船舶砰擊[2]、飛機(jī)迫降[3]等.物體入水問題的研究已經(jīng)持續(xù)了近一個(gè)世紀(jì),這個(gè)經(jīng)典問題的影響因素有很多:物體幾何形狀[4-6]、沖擊速度和角度[5]、物體密度和表面材質(zhì)[6]等.當(dāng)引入變化的物體密度[7]和物體入水帶有旋轉(zhuǎn)時(shí)[8],問題變得更加復(fù)雜.研究圓球高速旋轉(zhuǎn)入水過程的空泡特性和精細(xì)流場結(jié)構(gòu),可為將來基于多物體入水[9]的新型航行體入水降載方式提供關(guān)鍵性的力學(xué)機(jī)理認(rèn)識(shí).

入水問題的關(guān)注點(diǎn)主要包括入水沖擊和入水空泡兩個(gè)方面.Von Karman[10]提出的附加質(zhì)量方法為入水沖擊研究提供了先驅(qū)性理論指導(dǎo),其后國際上的理論研究大多圍繞該理論展開.對(duì)于圓球入水問題,Watanabe[11]開展的相關(guān)試驗(yàn)為后人理論模型的發(fā)展提供了關(guān)鍵的驗(yàn)證數(shù)據(jù).入水空泡的研究最早始于19 世紀(jì)末Worthington 和Cole[12]的工作,其采用電子閃光攝像捕捉圓球的垂直入水過程,首次系統(tǒng)地闡述了水面飛濺、空泡閉合和射流等現(xiàn)象.Truscott等[8,13-14]對(duì)圓球高速旋轉(zhuǎn)入水作了系統(tǒng)性的試驗(yàn)探究,包括采用不同入水速度和入水轉(zhuǎn)速,球面采用不同接觸角的涂層以及采用不同密度的圓球進(jìn)行試驗(yàn).他們發(fā)現(xiàn)由于圓球表面的無滑移邊界,使得圓球旋轉(zhuǎn)會(huì)在空泡內(nèi)部形成橫向楔形射流,并導(dǎo)致圓球彈道會(huì)發(fā)生明顯的橫向偏轉(zhuǎn).

近年來,高性能計(jì)算機(jī)的發(fā)展使得數(shù)值模擬成為求解入水問題越來越重要的手段.數(shù)值模擬在湍流場細(xì)節(jié)和渦結(jié)構(gòu)時(shí)空演變的精細(xì)捕捉方面,相比試驗(yàn)技術(shù)有天然的優(yōu)勢(shì).而對(duì)于入水問題,數(shù)值模擬方法有雷諾平均(RANS)類方法,大渦模擬(LES)類方法,也有部分直接數(shù)值模擬方法(DNS)[15],還有一些基于拉格朗日體系的無網(wǎng)格方法[16-18]等.

RANS 類方法是最常用的湍流模擬方法,Yu等[19]采用RANS 方法模擬了圓球的入水過程,發(fā)現(xiàn)無量綱化之后的沖擊力與入水速度和圓球半徑都無關(guān),只和球的質(zhì)量有關(guān).周波等[20]使用SST 模型對(duì)不同旋轉(zhuǎn)方向強(qiáng)制運(yùn)動(dòng)圓球的傾斜入水過程進(jìn)行數(shù)值模擬,發(fā)現(xiàn)球體繞不同軸旋轉(zhuǎn)傾斜入水時(shí)的運(yùn)動(dòng)受力特性的差異顯著.

然而傳統(tǒng)RANS 方法采用的渦黏湍流模型對(duì)于空泡流這種包含大范圍流動(dòng)分離和渦旋流動(dòng)的情形難以獲得較好的計(jì)算結(jié)果.相對(duì)而言大渦模擬類的高精度數(shù)值模擬方法(LES,DES,LES/RANS 混合等) 在湍流空泡流領(lǐng)域的應(yīng)用在近年來取得了長足的進(jìn)步.Li等[21]建立了六自由度的圓盤撞擊水面空泡演化的數(shù)值計(jì)算模型,他們發(fā)現(xiàn)沖擊生成的空泡在圓盤前進(jìn)方向存在非對(duì)稱與旋轉(zhuǎn)的現(xiàn)象.Zhao等[22]對(duì)低弗勞德數(shù)下圓柱斜入水問題進(jìn)行了數(shù)值模擬和試驗(yàn)探究,他們利用LES 與VOF (volume of fluid)相結(jié)合的方法求解湍流場和氣液界面.并且成功捕捉了多種尺度下圓柱周圍的渦結(jié)構(gòu),同時(shí)對(duì)多種入水參數(shù)進(jìn)行了數(shù)值模擬.

綜上所述,目前入水問題的研究方法眾多,而旋轉(zhuǎn)圓球入水還有許多復(fù)雜的現(xiàn)象值得深入探究,在這一方面采用LES 類方法進(jìn)行模擬的相關(guān)研究較少.本文采用LES 結(jié)合動(dòng)網(wǎng)格與滑移網(wǎng)格技術(shù),基于壁面邊界層八叉樹各向同性網(wǎng)格,對(duì)低弗勞德數(shù)下旋轉(zhuǎn)圓球入水的自由運(yùn)動(dòng)過程開展數(shù)值模擬,并與試驗(yàn)結(jié)果[7-8]進(jìn)行對(duì)比,著重分析空泡形態(tài)演化和水中彈道隨入水旋轉(zhuǎn)速度的變化規(guī)律,揭示圓球高速旋轉(zhuǎn)入水過程的空泡流場精細(xì)結(jié)構(gòu),為多物體入水降載問題提供基礎(chǔ)力學(xué)機(jī)理的認(rèn)識(shí).

1 數(shù)學(xué)模型

1.1 基本控制方程

本文采用均質(zhì)多相流的VOF 模型,氣液混合介質(zhì)由各相體積分?jǐn)?shù)比例混合而成,并將這種多相混合物看作一種單一的變密度流體.各相介質(zhì)共享相同的速度、壓力、湍流量等各種流動(dòng)物理量,混合物的物性參數(shù)則通過各相加權(quán)平均獲得.

混合介質(zhì)的流場控制方程如下

1.2 大渦模擬濾波的控制方程

在大渦模擬(LES)中,物理量的求解是通過濾波核函數(shù)G(x,x′) 在空間上對(duì)上述方程進(jìn)行濾波來得到的.濾波過程能很好的拆分尺度大于濾波長度的渦與尺度小于濾波長度的渦.濾波之后的變量如下定義

濾波后小尺度變量為 ?′=?-?ˉ,而在有限體積法的數(shù)值求解過程中,網(wǎng)格尺度V是隱式的濾波尺度,即G(x,x′)=1/V通過對(duì)式(1)和式(2)進(jìn)行濾波運(yùn)算之后,能得到以下濾波之后的控制方程

其中τij為亞格子應(yīng)力,定義如下

1.3 亞格子尺度應(yīng)力模型(SGS)

對(duì)N-S 方程濾波之后的亞格子應(yīng)力項(xiàng)是未知的,所以需要對(duì)其進(jìn)行建模來求解.采用Boussinesq假定[23],亞格子應(yīng)力計(jì)算方法如下

其中μt為亞格子渦黏系數(shù).采用常規(guī)濾波方法,所以各向同性部分τkk可以與濾波之后的壓力項(xiàng)合并[24],為濾波之后的應(yīng)變率張量定義如下

而對(duì)于亞格子渦黏系數(shù)μt,采用Nicoud 與Ducros[25]提出的WALE (wall-adapting local eddyviscosity)模型.該模型優(yōu)勢(shì)在于剪切層流的μt=0,相比傳統(tǒng)的Smagorinsky-Lilly 模型[26]能較好地返回流體區(qū)域中層流的特性,并且在壁面擁有更好的對(duì)數(shù)律y3性質(zhì).該模型具體形式如下

其中馮卡門常數(shù)κ=0.418 7,d為單元中心到壁面距離,ΔV為單元體積.WALE 常數(shù)采用Cw=0.325,這個(gè)值對(duì)于大多數(shù)流動(dòng)擁有更好的適應(yīng)性.

1.4 表面張力模型

當(dāng)氣液固三相接觸時(shí),在三相交界處產(chǎn)生接觸線與接觸角.當(dāng)溫度和壓強(qiáng)給定,在平衡狀態(tài)下這個(gè)接觸角被稱為靜態(tài)接觸角.然而通常在試驗(yàn)中能觀察到移動(dòng)的接觸線,這個(gè)接觸線對(duì)應(yīng)的接觸角為動(dòng)態(tài)接觸角(如圖1 所示).而在入水過程中動(dòng)態(tài)接觸角對(duì)空泡的產(chǎn)生與形態(tài)有很重要的影響[27].

圖1 入水參數(shù)與浸潤性示意圖Fig.1 The schematic description for the impact and wettability parameters

針對(duì)這個(gè)問題,本文采用 Brackbill等[28]提出的CSS (continuum surface stress)表面張力模型描述球體表面的潤濕性,將表面張力作為VOF 動(dòng)量方程的源項(xiàng),具體模型如下

CSS 模型避免顯式計(jì)算界面曲率,將表面張力影響作為各向異性的毛細(xì)應(yīng)力來計(jì)算,從而能對(duì)界面上的尖角得到較好的模擬結(jié)果.而對(duì)于氣液界面與球體壁面的接觸線,也引入了壁面依附約束模型(wall adhesion) 來處理.將流體與壁面的接觸角用于計(jì)算距離壁面的1 個(gè)流體單元界面法向的條件,如下所示

其中θw為壁面依附角,為壁面法向單位矢量,為壁面切向單位矢量.

2 數(shù)值算法與驗(yàn)證

本文采用有限體積法對(duì)流動(dòng)控制方程進(jìn)行時(shí)空離散,同時(shí)采用基于壓力的分離式求解器(SIMPLE)對(duì)控制方程中耦合的非線性項(xiàng)進(jìn)行解耦迭代求解.LES 需要解析流場中大部分高波數(shù)的湍流脈動(dòng),從而需要高精度的數(shù)值格式以及高精度的網(wǎng)格,而偶階次的數(shù)值格式天生不具有耗散性,很好的規(guī)避了對(duì)SGS 模型在耗散上的影響.所以采用中心差分格式能很好的適應(yīng)LES 方法,但是中心差分格式會(huì)對(duì)流場引入數(shù)值震蕩,所以最終我們選擇有界中心差分格式[29]作為動(dòng)量的離散格式.而對(duì)于方程中變量的梯度項(xiàng),采用最小二乘格式來進(jìn)行計(jì)算.由于入水過程為低速所以忽略空化效應(yīng)的影響.本文采用了VOF 多相流模型對(duì)氣液界面進(jìn)行捕捉,而迎風(fēng)格式對(duì)耗散會(huì)過估計(jì),中心差分則有數(shù)值震蕩都無法很好的對(duì)界面進(jìn)行捕捉.為了克服這些問題,采用一種高分辨率界面捕捉算法(HRIC)進(jìn)行空泡面捕捉[30],相比QUICK 與其他二階格式擁有更高的計(jì)算精度,相比傳統(tǒng)的幾何重構(gòu)格式計(jì)算速度能更快.圓球的三自由度運(yùn)動(dòng)采用動(dòng)網(wǎng)格方法中的動(dòng)框架來實(shí)現(xiàn)計(jì)算域平移運(yùn)動(dòng)控制,圓球帶空泡在水中平移和旋轉(zhuǎn)運(yùn)動(dòng)采用相對(duì)于背景計(jì)算域的滑移網(wǎng)格來控制.

2.1 計(jì)算域及邊界條件

為了校驗(yàn)當(dāng)前數(shù)值模擬方案的可靠性,與Aristoff等[3]的垂直入水試驗(yàn)進(jìn)行了對(duì)比.球直徑為D=0.025 4 m,半徑為R,分別對(duì)比了鋼和聚丙烯兩種不同密度比m*=ρl/ρs=7.86,0.86 的情況.如圖2 所示,整個(gè)計(jì)算域采用球直徑R進(jìn)行標(biāo)準(zhǔn)化,垂直長度為100R,水平長度為40R,寬度為20R,球心距離頂部壓力邊界設(shè)置為80R,距離側(cè)面速度邊界為20R,在寬度方向選取了對(duì)稱面,同時(shí)在球面上設(shè)置為無滑移邊界.水面位置距離球?yàn)?R.坐標(biāo)軸原點(diǎn)同水面保持一致如圖3 所示,其中Vl代表圓球入水時(shí)左側(cè)壁面速度大小,Vr代表圓球入水時(shí)右側(cè)壁面速度大小,圓球各項(xiàng)初始化參數(shù)如表1 所示.

表1 Aristoff等[7]垂直入水參數(shù)Table 1 Aristoff et al.[7] experimental condition

圖2 計(jì)算區(qū)域示意圖Fig.2 Computational domain of sphere

圖3 對(duì)稱面坐標(biāo)軸示意圖Fig.3 Coordinate system in cut view

具體計(jì)算網(wǎng)格采用結(jié)構(gòu)化與非結(jié)構(gòu)化相結(jié)合的方式,盡可能采用結(jié)構(gòu)化網(wǎng)格以減小數(shù)值誤差(如圖4 所示).為了實(shí)現(xiàn)精度較高的流場大渦模擬,近壁面邊界層附近采用盡量各向同性和均勻的網(wǎng)格,以減小大渦模擬低通濾波算法中濾波函數(shù)的交換誤差.計(jì)算網(wǎng)格采用盡量合理的網(wǎng)格疏密分布.在物體表面布置致密網(wǎng)格,網(wǎng)格的精細(xì)度達(dá)到能夠刻畫湍流邊界層結(jié)構(gòu)和典型湍流擬序結(jié)構(gòu)的程度.除了滿足數(shù)值模擬對(duì)網(wǎng)格質(zhì)量的一般性要求,即保證無量綱壁面垂向距離y+＜1 之外,還要盡量使壁面流向和展向的無量綱網(wǎng)格單元大小滿足Δx+≤100 和Δz+≤30.對(duì)于當(dāng)前模擬工況Δx+,Δz+列在表2 中,選取的時(shí)間為圓球剛剛沉入水面以下的時(shí)刻(t=10 ms 左右).

圖4 近壁面網(wǎng)格示意圖Fig.4 Illumination of near wall mesh

表2 流向與展向平均無量綱網(wǎng)格尺度表(t=10 ms)Table 2 Averaged nondimensional grid sizes of the finest mesh in the streamwise and spanwise (t=10 ms)

而另一種普遍采用的衡量LES 解析度的方式是Benard等[31]建議整體網(wǎng)格量能計(jì)算出流動(dòng)區(qū)域內(nèi)至少80%的湍動(dòng)能.將未解析湍動(dòng)能比率定義為

其中參數(shù)C選取為100.

雖然湍流脈動(dòng)速度和湍動(dòng)能ER在LES 中不是顯式求解獲得的物理量,但是對(duì)于那些宏觀時(shí)均流動(dòng)幾乎定常的流場,可以通過數(shù)據(jù)采樣統(tǒng)計(jì)的方法獲得和ER.雖然圓球入水空泡流動(dòng)并不是宏觀定常流場,但是可以采用同樣的網(wǎng)格對(duì)圓球在水中重力場中無空泡垂直降落過程這樣類似的時(shí)均定常流場進(jìn)行LES 模擬,進(jìn)行數(shù)據(jù)采樣統(tǒng)計(jì)從而獲得脈動(dòng)速度并計(jì)算湍動(dòng)能,以檢驗(yàn)所采用的網(wǎng)格分辨率的可靠性.所以模擬了圓球初始速度Vo=5.45 m/s 不帶空泡在水中下沉,其中t∈[0.07,0.18] (1100 時(shí)間步數(shù)據(jù)采樣),開啟數(shù)據(jù)采樣,得到均方根(RMS)的來計(jì)算ER與Esgs.

他逢到小小的失意也就哭，沒有一點(diǎn)隱藏。花生米翻落到地上了，自己嚼痛舌頭了，小貓不肯吃糕了，都是他哭的原因。有一天，他自己失手，把外婆送給他的泥人跌碎了，他哭得異常悲傷，和旁人死了親人沒有兩樣。

圖5 展示了對(duì)稱面上Qc與ER大小,可以發(fā)現(xiàn)在ER較高的圓球尾流區(qū)Qc遠(yuǎn)遠(yuǎn)小于0.2.所以可以認(rèn)為在流場中絕大多數(shù)區(qū)域Qc＜0.2 的要求都能被滿足.

圖5 對(duì)稱面上(a)大渦模擬解析湍動(dòng)能ER,(b)湍動(dòng)能未解析率Qc,采用數(shù)據(jù)采樣方法,采樣點(diǎn)1100 (Vo=5.45 m/s,ωo=0 rad/s)Fig.5 (a) Resolved turbulence kinetic energy ER and (b) its resolution Qc at the symmetry,obtained by the data sampling method using 1100 sampling points (Vo=5.45 m/s,ωo=0 rad/s)

根據(jù)經(jīng)典湍流理論,湍流場能量在慣性子區(qū)滿足-5/3 次方律(Kolmogorov 的K41 理論)[32]

而湍流壓力脈動(dòng)的能譜滿足文獻(xiàn)[33]的-7/3 次方律

這也是LES 研究論文中經(jīng)常采用的驗(yàn)證方法.為了進(jìn)一步驗(yàn)證所采用的網(wǎng)格分辨率的可靠性,實(shí)時(shí)監(jiān)測獲取穩(wěn)定之后流場中的壓力測點(diǎn)(在隨球運(yùn)動(dòng)的動(dòng)坐標(biāo)系內(nèi),測點(diǎn)坐標(biāo)x=0,y=3R,z=0)的壓力信息,采用快速傅里葉變換得到歸一化的功率譜密度估計(jì)PSD (見圖6).其中采樣頻率Fs=104,NFFT=1024.可以直接測量如下的雙對(duì)數(shù)坐標(biāo)PSD (power spectral density)圖中的峰值包絡(luò)線斜率,發(fā)現(xiàn)與上述經(jīng)典湍流理論中的-7/3 次方律非常接近.

圖6 歸一化的測點(diǎn)壓力功率譜密度及其與經(jīng)典湍流理論-7/3 次方律的對(duì)比(Vo=5.45 m/s,ωo=0 rad/s)Fig.6 Normalized PSD of the pressure in comparison with the-7/3 power spectrum of Kolmogorov’s K41 theory

2.2 垂直入水與試驗(yàn)對(duì)比驗(yàn)證結(jié)果

圖7 給出了圓球垂直入水?dāng)?shù)值模擬結(jié)果與文獻(xiàn)[7]的試驗(yàn)結(jié)果的對(duì)比.從圖中可以看出,當(dāng)前模擬方案在總體上對(duì)圓球深度有些低估,但總體上模擬與試驗(yàn)結(jié)果是非常接近的.表3 選擇了幾個(gè)特定時(shí)間,對(duì)比了圓球模擬深度與試驗(yàn)深度,可以發(fā)現(xiàn)相對(duì)誤差不大于10%,總體來說誤差水平是可以接受的.

表3 垂直入水圓球深度試驗(yàn)[7]與模擬對(duì)比表(Vo=2.17 m/s)Table 3 Comparison between the experimental[7] and simulation vertical depths for the case at Vo=2.17 m/s

圖7 垂直入水圓球深度試驗(yàn)[7]與模擬對(duì)比(Vo=2.17 m/s)Fig.7 Comparison between the experimental[7] and simulation vertical depths for the case at Vo=2.17 m/s

圖8 展示了數(shù)值模擬和試驗(yàn)泡型的對(duì)比.此處初始時(shí)刻定義同文獻(xiàn)[7]維持一致,即將圓球球心通過水面位置作為初始時(shí)刻.雖然在圖8(a) 68.9 ms 與圖8 (b) 62.2 ms 接近空泡斷裂時(shí)上半部分空泡略有偏小.但總體對(duì)于定性的泡形,當(dāng)前的模擬方案也能和試驗(yàn)吻合較好.特別是在空泡斷裂之后(圖8 (a)75.9 ms 與圖8 (b) 68.7 ms),還能捕捉到空泡面的收縮與向上的“Worthington”射流.

圖8 圓球垂直入水試驗(yàn)結(jié)果(上)與數(shù)值模擬結(jié)果(下)的對(duì)比Fig.8 Experimental (up) and simulation (down) results of the evolution of cavity shapes from the side view

圖8 圓球垂直入水試驗(yàn)結(jié)果(上)與數(shù)值模擬結(jié)果(下)的對(duì)比(續(xù))Fig.8 Experimental (up) and simulation (down) results of the evolution of cavity shapes from the side view (continued)

3 不同轉(zhuǎn)速對(duì)圓球旋轉(zhuǎn)入水影響

本文主要考察相同入水速度Vo=5.45 m/s 下不同入水轉(zhuǎn)速ωo=50,100,150,199 rad/s 對(duì)于入水空泡演化、圓球彈道、流場各個(gè)物性參數(shù)的影響.所有初始時(shí)刻t=0 ms 選取在圓球剛剛接觸自由液面的時(shí)刻.

3.1 入水角速度對(duì)圓球運(yùn)動(dòng)學(xué)性質(zhì)影響

本節(jié)將詳細(xì)討論圓球高速旋轉(zhuǎn)入水的水動(dòng)力學(xué)特性,主要關(guān)注在圓球的彈道(圖9)、速度(圖10)、加速度(圖11) 以及圓球受力情況(圖12) 隨時(shí)間t的變化.

圖9 不同入水轉(zhuǎn)速圓球入水彈道對(duì)比(ωo=199,150,100,50 rad/s,Vo=5.45 m/s)Fig.9 Trajectories of all spinning cases (ωo=199,150,100,50 rad/s,Vo=5.45 m/s)

圖10 圓球質(zhì)心速度與角速度隨時(shí)間t 變化Vo=5.45 m/sFig.10 The velocity of sphere as a function of time with Vo=5.45 m/s

圖11 圓球加速度與角加速度隨時(shí)間t 變化 Vo=5.45 m/sFig.11 The acceleration of sphere as a function of time with Vo=5.45 m/s

圖12 圓球所受升阻力系數(shù)隨時(shí)間變化Vo=5.45 m/sFig.12 Time-history of the simulation coefficients with Vo=5.45 m/s

圖9 展示了當(dāng)圓球入水速度為Vo=5.45 m/s 時(shí),不同入水旋轉(zhuǎn)角速度ωo=50,100,150,199 rad/s 對(duì)于圓球彈道的影響.其中為了和文獻(xiàn)[8]的試驗(yàn)結(jié)果進(jìn)行對(duì)比,使用密度比m*=1.74,半徑R=28.6 mm 的圓球進(jìn)行模擬.可以發(fā)現(xiàn)ω0=199 rad/s 的模擬結(jié)果與文獻(xiàn)[8]的試驗(yàn)結(jié)果吻合較好.圓球入水后在馬格努斯力的作用下彈道發(fā)生偏轉(zhuǎn),在入水的初始階段不同轉(zhuǎn)速ωo對(duì)彈道的影響并不是很大,但是隨著圓球進(jìn)入水面以下(t=10 ms 之后,見圖13),入水轉(zhuǎn)速越大的球,相應(yīng)的彈道偏轉(zhuǎn)也越大.當(dāng)空泡發(fā)生深閉合時(shí),空泡斷裂時(shí)圓球所在深度都非常近似.

為了進(jìn)一步探究彈道偏轉(zhuǎn)與入水轉(zhuǎn)速的關(guān)系,在圖10 中展示了圓球質(zhì)心速度與角速度隨時(shí)間變化關(guān)系圖,其中采用入水速度Vo與角速度ωo對(duì)速度和角速度進(jìn)行歸一化.通過比較歸一化之后的結(jié)果可以發(fā)現(xiàn),入水轉(zhuǎn)速的改變對(duì)水平方向速度vx/Vo的影響很大(見圖10 (a)).相應(yīng)的轉(zhuǎn)速越大,水平方向速度大小也越大,從而使得彈道偏轉(zhuǎn)也越大.對(duì)于ωo=50 rad/s,vx/Vo近似于線性下降,而對(duì)于ωo=100,150,199 rad/s 都隨時(shí)間t,先快速下降(t=10～20 ms),然后趨于常數(shù)(t＞ 60 ms).通過對(duì)比x方向加速度ax隨時(shí)間變化(見圖11(a))可以發(fā)現(xiàn),入水砰擊階段(t＜ 5 ms)ωo=50 rad/s 的加速度峰值不到ωo=100,150,199 rad/s 的30%,并且當(dāng)砰擊階段結(jié)束,ax也趨于平穩(wěn),這就導(dǎo)致了vx/Vo是近似的線性下降.而對(duì)于垂直方向速度vy/Vo,并不受到ωo改變的影響(見圖10 (b)),所有工況下都呈現(xiàn)在入水砰擊階段(t＜ 5 ms),下降較快,之后下降趨勢(shì)逐漸放緩的情況.而圓球轉(zhuǎn)速ω/ωo隨時(shí)間t變化總體不大(見圖10 (c),min (ω/ωo)=0.94).

在圖11 中展示了圓球加速度與角加速度隨時(shí)間t的變化,在入水砰擊階段 (見圖11(a)和圖11(b)t=0～4 ms),y方向加速度ay到達(dá)峰值之后快速下降,然后趨于不變并且不同ωo對(duì)于ay影響不大.但是對(duì)于ωo=199,150,100 rad/s 時(shí),ax繼續(xù)變大直到圓球大部分進(jìn)入自由液面以下(t=6 ms)而到達(dá)峰值.峰值的大小隨ωo變大而變大,而ωo=199,150 rad/s 兩個(gè)情況峰值大小非常接近,說明入水砰擊階段ax峰值大小受到Vo大小的限制.在t=20 ms 之后ax趨于穩(wěn)定直到空泡斷裂(t＞ 100 ms).圓球角加速度β(見圖11 (c)),在砰擊階段震蕩過后都基本維持不變,β的大小也是隨著ωo的變大而變大.結(jié)合圓球表面未沾濕面積之比Aair/Asphere隨時(shí)間t變化(見圖14),可以發(fā)現(xiàn)ωo越大,未沾濕比例越小,水依附在球面的面積越大,從而導(dǎo)致球面上受到黏性力也而越大,所以β隨著ωo的變大而變大.過了入水砰擊階段(t＞ 5 ms)Aair/Asphere基本維持不變,使得β維持不變.

圖14 圓球表面未沾濕面積Aair/Asphere 比例隨時(shí)間t 變化Vo=5.45 m/sFig.14 Time-history of non-wetted area ratio Aair/Asphere with Vo=5.45 m/s

基于上述運(yùn)動(dòng)學(xué)的討論,圓球質(zhì)心的受力可以做如下分析

其中ms為小球質(zhì)量,Fhx,Fhy為球面上所受黏性力與壓力的合力,球面上還應(yīng)有表面張力Fst對(duì)球質(zhì)心的作用,根據(jù)式(15)表面張力大小的最大值不會(huì)超過Fst=πDσ,其大小不到重力1%:6σ/(D2ρsg) ≈ 10-3,ρs為小球密度.所以為了簡化分析,忽略了表面張力對(duì)球受力的影響[15,34].

圓球所受升力Fl,阻力Fd的方向定義如圖15所示,故根據(jù)式(21)和式(22)可以將升阻力化簡為如下形式

圖15 圓球所受升阻力方向簡圖Fig.15 Schematic illumination of direction of drag and lift forces

其中彈道切線角度γ=arctan(-vx/vy),從而得到相應(yīng)的升阻力系數(shù)

其中A=0.5πR2為圓球截面面積.

根據(jù)上述定義在圖12 中展示了圓球所受升力系數(shù)Cl,阻力系數(shù)Cd隨時(shí)間t變化.可以發(fā)現(xiàn)Cl與Cd的變化趨勢(shì)分別與ax,ay比較接近.Cd在度過入水砰擊階段之后(t＞5 ms)從Cd=0.2 開始隨時(shí)間線性下降直到空泡斷裂,并且ωo對(duì)其影響并不是很明顯.但是對(duì)于升力系數(shù)Cl,ωo越大入水初始階段(t=0～10 ms)的峰值越大.

3.2 入水角速度對(duì)空泡形態(tài)的影響

本節(jié)展示了ωo=199,150,100,50 rad/s 的空泡泡形(圖13),對(duì)稱面上壓力、速度云圖,以及Q準(zhǔn)則等值面,其中Q準(zhǔn)則為速度梯度張量的第二不變量定義如下[35].

圖13 圓球入水空泡演化對(duì)比圖(Vo=5.45 m/s)Fig.13 Time evolution of cavity during water entry of the spinning spheres from the longitudinal view (Vo=5.45 m/s)

圖13 展示了空泡隨時(shí)間的演化序列,其中球面上為水氣混合密度的云圖,即不同顏色的接觸線代表了圓球的水氣接觸線,在圖13(a) 中定性比較了文獻(xiàn)[8]的試驗(yàn)與當(dāng)前的模擬結(jié)果.可以發(fā)現(xiàn)模擬泡形相比試驗(yàn)結(jié)果略小一些,這在空泡斷裂之后較為明顯(圖13 (a)實(shí)驗(yàn)111 ms 與模擬110.8 ms).這可能是由于隨著圓球的下沉網(wǎng)格逐漸稀疏導(dǎo)致的.同時(shí)在試驗(yàn)中飛濺冠在液面以上先行閉合,而在當(dāng)前模擬中無法很好的模擬水面以上的飛濺冠使得這一現(xiàn)象無法觀察到.在圖13 和圖16 可以看到有橫向楔形劈尖在入水的初始階段(t=10 ms 附近)開始生成侵入空泡內(nèi)部,結(jié)合圖11中ax隨t變化可以發(fā)現(xiàn),橫向劈尖在球面上的攀升是ax變大的主要原因.當(dāng)t＞20 ms 時(shí),橫向劈尖在球面上不再繼續(xù)延展 (見圖13(a)～圖13(c))ax也趨于穩(wěn)定.隨著圓球的下沉,橫向劈尖逐漸發(fā)展直到擊穿泡面(ωo=199,150 rad/s).而相同時(shí)刻,ωo越大橫向劈尖也越大(見圖16).而對(duì)于ωo=100 rad/s,上半部分楔形劈尖還未到達(dá)另一邊空泡面空泡已經(jīng)發(fā)生斷裂.

圖16 圓球入水空泡演化俯視圖(Vo=5.45 m/s)Fig.16 Time evolution of cavity during water entry of the spinning spheres from the top view ( Vo=5.45 m/s)

對(duì)于ωo=50 rad/s 的情況(圖13,圖16 (d)),能觀察到空泡在自由液面上發(fā)生表面閉合,也沒有明顯的橫向劈尖產(chǎn)生(圖16 (d)),這是由于ωo=50 rad/s 對(duì)應(yīng)的無量綱旋轉(zhuǎn)參數(shù)S=0.262＜1,對(duì)應(yīng)Vr＜0,也就是接觸線大部分都在圓球赤道以下,所以沒有橫向劈尖的產(chǎn)生,與文獻(xiàn)[8]的結(jié)論保持一致.

入水初始階段(t=10 ms 附近)與空泡斷裂后(t=110 ms 附近)不同ωo的壓力云圖展示在圖17 中,其中實(shí)線代表混合密度ρm=400 kg/m3自由液面,可以看出由于馬格努斯效應(yīng)圓球一側(cè)形成低壓區(qū)并且隨著ωo變大,低壓區(qū)也在變大,這也使得彈道的橫向位移也越大.同時(shí)ωo的變化對(duì)于空泡斷裂后的壓力場有顯著的影響.對(duì)于ωo=100 rad/s (見圖17(c)),在空泡的斷裂位置有明顯的高壓區(qū),但是對(duì)于ωo=199 rad/s (見圖17(a)),并不能觀察到明顯的高壓區(qū).這可以如下解釋:橫向劈尖侵入泡內(nèi)從而影響空泡深閉合的斷裂,而ωo越大橫向楔形越大,從而使得水進(jìn)入泡內(nèi)更快.最終結(jié)果就是空泡因靜水壓夾斷時(shí),產(chǎn)生的高壓區(qū)更小乃至幾乎沒有.

圖17 圓球?qū)ΨQ面壓力云圖(Vo=5.45 m/s)Fig.17 Pressure contours of the spinning spheres from symmetric view(Vo=5.45 m/s)

圖18 展示了ωo=199,150,100 rad/s 在100 ms 附近時(shí)的對(duì)稱面速度云圖與矢量場,對(duì)比圖13可以發(fā)現(xiàn)雖然此時(shí)空泡還未斷裂,但在對(duì)稱面上由于橫向楔形劈尖的存在,已經(jīng)接近斷裂.觀察速度云圖可以發(fā)現(xiàn)ωo=100 rad/s 在自由液面上沒有出現(xiàn)明顯的橫向氣流射出與ωo=199,150 rad/s 有所不同,這可以被解釋為ωo較小,橫向劈尖速度小,使得自由液面變形速度較慢從而抑制了橫向氣流的生成.

圖18 圓球?qū)ΨQ面速度大小與矢量場(Vo=5.45 m/s)Fig.18 Velocity and vector field of the spinning spheres from the symmetric view ( Vo=5.45 m/s)

為了探究渦結(jié)構(gòu)對(duì)空泡泡形的影響,圖19 展示了Q=100 000 s-2的等值面,可以發(fā)現(xiàn)渦結(jié)構(gòu)生成的時(shí)間點(diǎn)有兩個(gè):入水初始階段(t=11 ms 附近)與空泡斷裂階段(t=90～111 ms).這一點(diǎn)與垂直入水的情況保持一致,但橫向劈尖的存在使得空泡上半部分始終能觀察到相應(yīng)的渦結(jié)構(gòu).ωo的增加使入水初始階段的渦結(jié)構(gòu)更加趨向于非對(duì)稱,但是到了空泡斷裂階段,ωo越大渦結(jié)構(gòu)越少.這可以被解釋為小球橫向位移的增加,使得泡內(nèi)空氣與外界交換受到阻礙,抑制了渦結(jié)構(gòu)的生成.

圖19 旋轉(zhuǎn)圓球入水Q=100 000 s-2 等值面 (Vo=5.45 m/s)Fig.19 Vortex Q iso-surface of the spinning spheres from longitudinalview (Vo=5.45 m/s)

4 結(jié)論

本文采用大渦模擬對(duì)圓球高速旋轉(zhuǎn)入水建立了數(shù)值模型.研究了圓球在相同入水速度Vo=5.45 m/s 下不同入水轉(zhuǎn)速ωo=50,100,150,199 rad/s的空泡演化與水動(dòng)力學(xué)特性,主要結(jié)論如下.

(1)數(shù)值模擬結(jié)果通過與垂直入水,高速旋轉(zhuǎn)入水試驗(yàn)對(duì)比,在圓球彈道,空泡泡形上與試驗(yàn)結(jié)果有較好的一致性,驗(yàn)證了數(shù)值模型的正確性與可靠性.

(2)入水轉(zhuǎn)速的改變對(duì)圓球的水動(dòng)力學(xué)特性有顯著的影響.隨著入水轉(zhuǎn)速的增加,圓球彈道水平位移更大,圓球水平方向獲得的速度更大,在入水砰擊階段圓球所獲得的升力也更大,但升力的峰值受到入水速度的限制.而入水轉(zhuǎn)速對(duì)于圓球垂直方向特性影響不大,包括垂直方向速度,加速度,空泡斷裂時(shí)圓球所在深度.

(3)入水轉(zhuǎn)速對(duì)于空泡形態(tài)與流場特性得到詳細(xì)的探究.在模擬中觀察到ωo=50 rad/s 的表面閉合ωo=100,150,199 rad/s 的深閉合.對(duì)于深閉合,入水轉(zhuǎn)速的增加帶來了更強(qiáng)的橫向楔形射流并且抑制了空泡斷裂帶來的高壓以及相應(yīng)渦結(jié)構(gòu)的生成.