水下發射航行體尾渦不穩定性分析1)

高山 施 瑤, 潘 光 權曉波 魯杰文

* (西北工業大學航海學院,西安 710072)

? (無人水下運載技術工業和信息化部重點實驗室,西安 710072)

** (中國運載火箭技術研究院,北京 100076)

引言

在實際工程中,尾流常常由復雜的三維鈍物體產生.即使是一些非常簡單的外形,其尾流渦旋結構演變就比圓柱、翼型以及球體等二維物體的尾流結構復雜.在水下垂直發射領域中,常常伴隨著復雜的三維湍流流動現象,具有強非線性、高湍流性以及混沌性等典型特征[1].尤其在近尾流湍流區域內,存在著大量尺度不一和強度各異的渦旋結構[2],這些渦旋結構在湍流生成和維持過程中發揮著極其重要作用,同時也是水下連續發射尾流干擾難題的關鍵影響因素[3].

目前,關于水下發射大多數成果集中在空化水動力與流固耦合特性[4-9].王一偉等[10]在水下發射非定常空化流動方面系統地總結了空化穩定性、潰滅載荷特性以及流動控制等重要物理機制.尤天慶等[11]基于勢流理論分析了穩態與非穩態尾空泡流場演變特性,探索了尾空泡瞬態收縮過程對上游航行體物面壓力的擾動機理.雖然空化流動現象在水下發射過程顯得至關重要,然而由于水下連續發射環境的特殊性,航行體流動干擾可能由其他因素所主導,其中尾流效應是其中最重要的一個影響因素.

長期以來,三維航行體的尾流演變機理都是學者們非常感興趣的研究課題.因此,學者們為了進一步認識軸對稱鈍物體的三維尾流渦旋結構,以圓柱、方柱等物體為例開展了大量的研究.李永光等[12]通過實驗結果和理論分析相結合方法,首次得出了當有穩定的氣液兩相渦街發生時,其渦街結構參數與單相流不一致.另外,王智慧等[13]采用PIV 測量手段研究了不同雷諾數下尾跡的速度矢量場和渦量場,發現尾渦的長度、寬度以及脫落頻率與雷諾數緊密相關.在此基礎上,學者們取得了大量的研究成果[14-18].隨著科學技術的發展和研究方法的進步,原本的很多限制逐步被打破,研究對象也從原本的二維結構轉向三維結構.高洋洋等[19]研究探索了三維圓柱在不同雷諾數和不同傾斜角度的尾渦流動特性.劉闖等[20]采用大渦模擬方法模擬了高雷諾數下三維圓柱的尾渦結構,發現了高雷諾數下圓柱繞流尾跡變化具有不穩定性.另外,大量研究表明三維圓柱繞流模擬明顯比二維更符合實際[21-22].Chen 等[23-24]在三維結構基礎之上,研究了湍流動能耗散率和溫度耗散率在尾跡中的空間分布,相平均結果表明湍動能耗散率和溫度耗散率均集中在卡門渦結構的內部.然而,實際中三維鈍物體的尾流常常以大尺度結構存在.Taneda 等[25]提供了油滴在靜止水中下落時的尾跡相圖,渦環從渦面卷起時立刻失去各自特征,相互生成、相互滲透.Rosenhead 等[26]基于實驗結果對此也有相關的解釋: 鈍物體的尾跡是由一系列不規則的渦環組成,他們的方位是完全隨機的,是由渦等脫落位置所確定.夏雪湔等[27]對70°斜切尾鈍頭旋成體的尾渦結構開展了相關實驗研究,發現這類鈍物體尾跡中的三維尾渦結構均呈現為多個發卡渦相互連接狀態.Shi 等[28-29]也針對水下發射尾流渦結構演變開展了細致的研究,即采用RANS 方法對水下連續發射尾流干擾進行了初步探究,發現尾流中沿軸向間隔排列組成的發卡渦包結構,并對次發航行體的表面壓力分布和運動穩定性均有較大的影響.

目前關于尾流渦旋演變機理的相關研究對象主要集中在二維圓柱以及少數的三維鈍物體等,然而針對航行體水下發射三維尾流渦旋結構演變機理還處于起步研究階段.本文的研究工作擬利用改進型延遲分離渦方法對三維航行體水下發射尾流演變過程開展精細化模擬,通過渦識別方法對其尾流區渦旋結構進行識別,研究尾渦結構的演變機理和流場的脈動壓力規律,分析不同無量綱橫流強度和不同雷諾數下尾渦結構演變差異和脈動壓力規律,以期為解決水下連續發射尾流干擾問題提供分析方法和手段.

1 數值方法

1.1 控制方程

水下發射過程涉及水與非凝結氣體的混合相流動現象,各相之間存在著強烈的相互作用,導致水下發射過程尾流渦旋結構演變問題具有強烈的非線性和高度的耦合性.本文采用VOF 多相流模型對混合流場進行描述,基本控制方程形式如下

體積分數方程

其中,混合介質密度 ρm=αkρk.

動量方程

其中,μm,λm,F分別表示混合介質動力黏度、混合介質第二黏度以及體力項.

能量方程

其中,Ea和Ta分別代表混合介質平均能量和平均溫度;k和 Φ 分別是熱傳導率和黏性耗散項.

1.2 湍流模型

本文采用的延遲改進型分離渦IDDES (improved delayed detached eddy simulation)模型基于SSTk-ω模型進行構造,引入湍流長度尺度lIDDES,對模型中湍動能耗散項進行了修正,可以寫成

其中,km和 τ 分別是湍流動能和湍流剪切力,Sij是應變率張量,lIDDES的基本形式為

式中,lRANS和lLES分別是RANS 長度尺度和LES 濾波尺度,可表示為

式(10) 中,CDES為模型系數,通常取為0.65;Δmesh為網格尺度;dw為計算點到壁面的距離;hmax為網格最大邊長;hwn為垂直壁面方向的網格尺度.另外,本文研究航行體所處的復雜水下環境中最小空化數遠遠大于半球頭型航行體的初始空化數,加之抗空化頭型和試驗驗證.綜合分析下,本文未考慮空化影響.

1.3 邊界條件與網格細節

如圖1 所示,航行體直徑為D,長徑比L/D=6,其中L是航行體的軸向長度.發射筒軸向長度為 1.2L,直徑為D.即發射筒壁面與航行體壁面緊密接觸,主要防止筒內高壓氣體的泄漏.在整個水下發射模擬過程中,采用底部通入高壓氣體方式將航行體彈射出筒,其中航行體在筒內沿著Z軸正方向加速運動.當航行體尾端離開發射筒口之后,以多自由度方式在水中航行運動,從而完成整個水下發射過程.

圖1 計算域和邊界條件示意圖Fig.1 Schematic diagram of calculation domain and boundary conditions

本文定義橫向來流速度u∞沿著X軸正方向,重力方向沿著Z軸負方向.計算流域長度為 4L,寬度4L,高度 5L.另外,水域高度為 4L,空氣域高度為L.計算域左側為速度入口,底面是壁面,發射筒底端為滯止入口,其余邊界條件均為壓力出口.滯止入口的壓力大小隨時間變化,當航行體尾端即將出筒時下降到筒口附近靜水壓力.為了控制計算量和提高計算效率,沿縱向截取流場計算域,采用二分之一模型開展計算.圖2 顯示了網格劃分細節,背景域和重疊域網格均采用切割體網格,而邊界層網格為均勻拉伸的棱柱層網格.為了捕捉尾流區域渦旋結構演化細節,同時與改進型分離渦模型合理匹配,網格細化率取1.3.另外,第一層邊界層高度1.5 mm,背景加密域和重疊域網格大小均為0.025D,網格總數是2.8 ×106,并對航行體運動區域和水面附近區域進行網格加密.

圖2 網格劃分細節Fig.2 Meshing details

1.4 重疊網格技術

重疊網格方法,又稱為嵌套網格技術,其計算原理是: 計算域網格被分割為多塊具有重疊或嵌套部分的子網格,其中貼體網格跟隨運動體一起運動.數值模擬計算在各個分塊子網格上分別進行,而流場信息的傳遞在重疊邊界上通過插值的方式進行,有利于數值計算.因此重疊網格的網格生成難度大大降低,且對模擬多體耦合運動和大轉角運動有著很大優勢.計算流程如圖3 所示.

圖3 重疊網格工作流程Fig.3 Workflow of overlapping grids

網格裝配主要實現流程為尋點、挖洞和建立插值關系等過程.其中尋點為尋找物理空間點在網格中的相對位置,通過識別包含該空間點的網格單元實現.挖洞則是將背景區域和重疊區域進行耦合,在這個過程中,完全取自重疊區域的網格單元被標記為背景區域中的非活動網格單元,并將非活動網格單元從計算域中刪除,從而實現挖洞過程.

1.5 數值計算方法驗證

本文設計并搭建了水下發射實驗測試裝置,并對航行體水下發射過程開展實驗與數值計算結果對比與分析.如圖4 所示為航行體水下發射試驗平臺實物圖,主要由發射系統、艇速控制系統以及數據采集系統組成.其中,數據采集系統由控制電腦、高速攝像機(Phantom 型)及相應線路等組成.調節白平衡、分辨率(640×1024)、幀率(4000 幀/秒)等參數,以拍到清晰的畫面.實驗中采用泡沫板來抵消航行體出水慣性載荷的沖擊,既保證了實驗人員的安全性,同時又防止航行體頭型發生撞擊而損壞.

圖4 實驗裝置Fig.4 Experimental device

對實驗圖像的后處理主要是對航行體所處的空間位置的計算和捕捉.根據航行體在整個運動過程中所處于圖像中的位置對圖像進行區域劃分,從而大大減小處理過程中的運算量,同時可以避免其余不必要的背景干擾,提高邊緣檢測精度和效率.為了獲得航行體質心在坐標系下的實際位置,在實驗圖像中,通過對航行體的標定,可測得圖像中每個像素點之間的距離與真實距離的比例關系,通過比例關系可以得到質心以及空泡邊緣在坐標系下的實際位置.然而,在實驗過程中因為高速攝像機拍攝的實驗圖像存在有因光線穿過水體和空氣間產生的折射誤差,需要對提取的邊緣散點進行折射校正,校正公式如下

其中,Hmeasure是運動體運動深度的測量值,Hreal是運動體運動深度真實值,d1是運動體距水箱前壁面的距離,d2是高速攝像機鏡頭距水箱前壁面的距離,n為折射率,如圖5 所示.

圖5 折射校正Fig.5 Refraction correction

基于上述介紹的水下發射實驗測試系統和數據處理方法,本文開展了直徑20 mm,長度120 mm 的航行體在出筒速度14.3 m/s、平臺移動速度0.5 m/s下的實驗和數值模擬對比驗證工作.圖6 和圖7 分別給出了數值模擬與實驗結果相圖對比和運動位移曲線對比.在水下航行階段,仿真結果和試驗吻合度較好.然而,試驗中尾空泡由于發射筒口不均勻氣團效應導致與模擬結果存在一定的誤差.在運動位移方面,數值模擬與實驗測試的最大誤差為7.52%,驗證了本文采用數值模擬方法的合理性.

圖6 實驗與數值計算相圖對比Fig.6 Comparison of phase diagram between experiment and numerical result

圖7 實驗與模擬計算運動位移對比Fig.7 Comparison of motion displacement between experiment and simulation result

2 結果分析

2.1 瞬態尾流場分析

航行體水下發射過程一般可分為三階段: 出筒階段、水中航行階段以及出水階段[30].如圖8 所示為水下發射過程氣液相體積分數演變云圖.在橫流效應下,航行體帶有一定的攻角離開發射筒口時,筒內的高壓氣體瞬間失去約束,從有界流域向無界流域迅速膨脹.流場中的氣團分為三部分: 與管出口相連的氣團、附著在航行體尾端的氣團(尾空泡),以及在流場內游離的筒口預置氣團.由于氣團的膨脹、收縮、運動和碰撞等狀態變化,將直接對尾流區湍流渦旋流場演變規律具有重要的影響.

圖8 不同時刻下流場氣-液相體積分數演化Fig.8 Evolution of gas-liquid volume fraction in flow field at different times

本文涉及的無量綱參數有無量綱時間T、無量綱時間橫流強度U以及雷諾數Re,分別定義如下

式中,ρ∞和 μ 分別是液態水的密度和黏性系數;u∞和ve分別代表是橫流速度和航行體出筒速度;t和μ分別是航行體運動時間和動力黏性系數.

航行體水中航行階段中,由于尾流高速流體核心區與低速自由流相互作用,呈現出Kelvin-Helmholtz不穩定現象的典型特征,這主要是由于兩種速度差混合流中剪切運動所引起的[31-32],如圖9 所示.在水中航行階段,航行體表面附近不斷積累渦量.在流動分離作用下,背流側積累的渦量較多,并在航行體尾端形成脫落渦結構環進入尾流區內.從圖10 可知,脫落渦環的核心區交替出現,且沿著X軸正方向運動并發生耗散現象.

圖9 不同時刻下流場速度幅值演化Fig.9 Evolution of velocity amplitude in flow field at different times

圖10 不同時刻下流場渦量幅值演化Fig.10 Evolution of vorticity amplitude in flow field at different times

圖11 所示為有無橫流條件下航行體頭部觸及自由液面時刻渦量幅值與等值面尾渦圖.由于Q方法難以識別衍生渦以及二次渦,而 λci方法在衍生渦以及二次渦識別過程表現較好[33].因此,本文采用λci準則對尾流區渦旋結構進行了識別.其中,λci準則是在 Δ 準則的基礎上進一步發展而來,其數學定義為

圖11 渦量與等值面尾渦演變Fig.11 FIig.11.Evolution of vorticity and isosurface wake vortex

當Δ＞0,其特征值為 λ1=λr,λ2,3=λcr±iλci其中

其中,P,Q,R是速度梯度張量 ?V的三個伽利略不變量. tr 代表矩陣的跡,det 代表矩陣的行列式.

無橫流條件下,尾渦結構近似沿著航行體軸線兩側隨機分布,而且脫落渦環極易發生破碎,并未形成結構完整的尾渦結構.然而,在橫流效應下,當航行體表面附著渦量到達其尾端時,由于其壁面發生旋轉,形成低壓渦旋,導致流動圍繞周圍渦核中心,并卷起形成渦環,與橫向來流相互作用,發生了準周期性脫落,形成多個渦環,也稱為渦環包[29].在此過程中,渦管沿著流向發生拉伸形成渦腿,渦環包與渦腿形成弧狀結構發卡渦[34].總體來看,尾渦結構由一系列“發卡渦”相互鎖定的“發卡渦包”組成[27],其形態具體表現為不規則狀.另外,隨著發卡渦包結構的形成,二次渦結構在尾流中也不斷衍生并發展.圖12顯示了不同時刻下等值面尾渦結構演變圖.從圖中可以發現,單個發卡渦一般不會在航行體尾流區單獨存在,而是由多個發卡渦沿軸向間隔排列,組成發卡渦包存在于航行體尾流中.

圖12 不同時刻下流場等值面尾渦結構演化Fig.12 Evolution of isosurface wake vortex at different times

為了研究尾流場脈動壓力的演變規律,在空間流場內等間距不同位置處設置檢測線,如圖13 所示,其中壓力系數Cp的具體定義為

圖13 流場監測線分布Fig.13 Distribution of flow field monitoring lines

式中,p∞分別是無窮遠處當地絕對靜壓值.

從圖14 可以看出,不同監測線出現的第一次峰值時刻均為航行體頭部穿過監測位置時刻.然而,進一步分析可以發現,不同監測線的脈動壓力曲線均出現了二次峰值,而P4監測線甚至出現了三次峰值,并逐漸減小.在脈動壓力的二次峰值出現時刻,此時航行體尾端已離開此位置,因此主要由尾流區渦旋結構演變所主導.

圖14 不同監測線脈動壓力演變Fig.14 Evolution of fluctuating pressure of different monitoring lines

2.2 不同橫流強度下尾流結構演變

為了研究橫流強度對尾流場渦旋結構演變機制的影響規律.本文模擬了不同無量綱橫流強度U下精細化尾流場演變過程.基于上述無量綱橫流強度定義,本文開展了U=0.00,U=0.03,U=0.06,U=0.09 以及U=0.12 下尾流場渦旋結構演變過程研究.如圖15 和圖16 所示是航行體頭部觸及自由液面時不同橫流強度下尾流場速度與渦量云圖.當無量綱橫流強度U=0.00 時,即沒有橫流作用下,尾流高速度核心區與渦核極易發生隨機破碎現象,并以較小形態游離在尾流場中,近似沿著航行體軸線兩側分布.隨著橫流強度的增大,尾流場高速度區破碎現象逐漸減小.當無量綱橫流強度U=0.03 時,尾流區速度核心區和渦核雖然結構上發生分離現象,但有緊密靠近趨勢.隨著橫流強度進一步增大,航行體背流側流動分離強度增加,尾流速度核心區和渦核強度均不斷增加,其形態上并未發生明顯的分離現象.同時,Kelvin-Helmholtz 不穩定現象的特征愈發明顯.另外,當橫流強度較小時(U=0.03),渦核近似以小尺度球狀形態存在;在橫流強度較大時(U=0.12),渦核近似以大尺度片狀形態存在于尾流中.

圖15 不同橫流強度下尾流場速度分布Fig.15 Velocity distribution of wake field under different crossflow intensity

圖16 不同橫流強度下尾流場渦量分布Fig.16 Vorticity distribution of wake field under different crossflow intensity

圖17 給出了不同橫流強度下等值面尾渦云圖.當無橫流強度下,尾渦結構主要以少數脫落渦環以及大量游離的小尺度渦結構組成,并未形成明顯的發卡渦結構.當橫流強度為U=0.03 時,渦管沿著流向拉伸并向上翹起形成渦腿,并與脫落的渦環形成發卡渦結構.進一步可以發現,此時尾流區主要存在兩級渦結構,主渦較為明顯的發卡渦包結構,而二級渦結構從尾流中脫落并在流場中形成.在主渦影響下,尾渦結構主要以大尺度不規則渦結構存在于尾流中.隨著橫流強度不斷增大,主渦中發卡渦包尺度不斷增強,形成多級準周期的發卡渦結構.當橫流強度為U=0.12 時,航行體尾部渦環的脫落頻率達到峰值,此時渦管形成渦腿的頻率不足以連接已發生脫落渦環的,從而導致發卡渦包形態變得極其不規則.隨著橫流強度的增加,航行體出筒時刻初始攻角越大,繞流現象愈發明顯,導致航行體尾端脫落的渦環頻率增加、尺度增大,從而導致尾流中發卡渦包的穩定性減小.

圖17 不同橫流強度下等值面尾渦分布Fig.17 Wake vortex distribution of isosurface under different crossflow intensity

圖18～圖20 分別是不同橫流強度下P2,P3和P4監測線處脈動壓力演變曲線.不同橫流強度下,由于尾流引起的脈動壓力二次峰值均出現.隨著橫流強度的增大,P2監測線二次脈動壓力峰值先增大后減小,P3監測線二次脈動壓力峰值不斷增大,而P4監測線出現了三次壓力峰值現象,且峰值隨著橫流強度的增加而增加.因此,可以發現,隨著橫流強度的增加,尾渦結構演變可近似總結為離散的小尺度渦結構-多級準周期的發卡渦包-大尺度不規則多級渦環碰撞,導致壓力出現多次峰值.

圖18 不同橫流強度下P2 脈動壓力演變Fig.18 Evolution of P2 fluctuating pressure under different crossflow intensities

圖19 不同橫流強度下P3 脈動壓力演變Fig.19 Evolution of P3 fluctuating pressure under different crossflow intensities

圖20 不同橫流強度下P4 脈動壓力演變Fig.20 Evolution of P4 fluctuating pressure under different crossflow intensities

2.3 不同雷諾數下尾流結構演變

為了研究雷諾數對尾流場渦旋結構演變機制的影響規律,本文模擬了不同雷諾數下精細化尾流場演變過程.基于上述雷諾數定義,本文開展了Re=1.68×105,Re=2.77×105,Re=3.16×105,Re=3.56 ×105下尾流場渦旋結構演變過程研究.

圖21 和圖22 所示為當航行體頭部觸及自由液面時不同雷諾數下尾流場速度與渦量云圖.隨著雷諾數的增大,渦環脫落頻率逐漸減小,但尺度逐漸增大.即隨著雷諾數的增大,尾渦結構演變可近似總結為: 小尺度-中尺度-大尺度,隨機性也不斷增強.圖23 為不同雷諾數下等值面尾渦云圖.另外,隨著雷諾數的增大,發卡渦包中渦頭的數量明顯減小,不規則性增強.值得注意的是: 在高雷諾數下,衍生渦結構逐漸明顯,其中衍生渦結構由圓柱形渦和U 型渦組成.當雷諾數較大時(Re=3.56×105),尾流區渦結構演變較為激烈,發卡渦之間相互作用,隨機性加強,其基本形態特征消失,逐漸演變為大尺度不規則渦結構.另外,衍生渦結構極為明顯,并與主渦結構相互作用,相互影響.圖24～圖26 分別是不同橫流強度下P2,P3和P4監測線處脈動壓力演變曲線.可以發現不同監測處脈動壓力二次峰值不盡相同,但其渦結構演變對尾流場均產生了一定的擾動現象.

圖21 不同雷諾數下尾流場速度分布Fig.21 Velocity distribution of wake field under different Reynolds number

圖22 不同雷諾數下尾流場渦量分布Fig.22 Vorticity distribution of wake field under different Reynolds number

圖23 不同雷諾數下等值面尾渦云圖Fig.23 Wake vortex distribution of isosurface under different Reynolds number

圖24 不同雷諾數下P2 脈動壓力演變Fig.24 Evolution of P2 fluctuating pressure under different Reynolds number

圖25 不同雷諾數下P3 脈動壓力演變Fig.25 Evolution of P3 fluctuating pressure under different Reynolds number

圖26 不同雷諾數下P4 脈動壓力演變Fig.26 Evolution of P4 fluctuating pressure under different Reynolds number

3 結論

本文基于改進型分離渦模型、VOF 多相流模型以及嵌套網格零間隙技術,建立了橫流效應下航行體水下發射數值計算模型,開展了水下發射精細化尾流場數值模擬研究,分析了瞬態尾流場中氣液相體積分數、速度、渦量以及渦結構演變,并且進一步討論了無量綱橫流強度(U=0～0.12)和雷諾數(Re=1.68×105～3.56×105)對尾流場中渦旋結構和脈動壓力分布特性的影響,具體的結論如下.

(1)水下發射過程尾流區高速流體核心區與低速自由流相互作用,在兩種速度差混合流中剪切效應引起尾流區呈現明顯的Kelvin-Helmholtz 不穩定現象特征.在橫流條件下,渦管沿著流向發生拉伸形成渦腿,脫落的渦環包與渦腿形成發卡形的弧狀結構發卡渦.總的來看,發卡渦一般不會單獨存在于水下發射尾流中,其結構形態近似呈現由一系列“發卡型”渦相互鎖定的不規則的“發卡渦包”.

(2)隨著橫流強度的增大,尾流速度核心區和渦核強度增加,緊密貼合程度增強.另外,隨著橫流強度增大,主渦中發卡渦包尺度增強,形成多級準周期的發卡渦包結構.當渦環的脫落頻率較大時,渦管形成渦腿不足以連接多個渦環,從而導致發卡渦包形態特征變得不規則,且隨機分布性增強.不同橫流強度下,導致脈動壓力二次峰值均出現的主要原因是尾流渦旋流場演變引起的.

(3)隨著雷諾數的增大,衍生渦結構逐漸顯現,主要由圓柱形渦和U型渦組成.當雷諾數較大時,發卡渦之間相互作用的強度明顯增強,導致其基本形態特征消失,逐漸演變為大尺度不規則渦結構,不穩定性增強.