超空泡航行體流體動力延遲特性數值模擬研究

2023-11-21 12:56:58劉喜燕袁緒龍羅凱杜兆星

西北工業大學學報 2023年5期

劉喜燕, 袁緒龍, 羅凱, 杜兆星

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

由于空泡延遲效應[1]的存在,當超空泡航行體改變攻角時,新生成的空泡需經過一段時間才能到達航行體尾部。空泡狀態的更新相對于航行體姿態的變化存在一定的時間滯后,沾濕面積改變,從而引起流體動力的變化。空泡與航行體之間相對位置關系的實時變化,導致超空泡航行體機動航行過程中的流體動力展現出極強的復雜性和非定常特性,從而使超空泡航行體彈道預報和穩定控制充滿挑戰。因此,研究超空泡航行體流體動力的延遲特性,具有十分重要的理論與工程應用價值。

目前國內外學者圍繞超空泡航行體機動航行過程中的空泡形態和非定常流體動力特性進行了大量研究。李魁彬等[2]基于空泡截面獨立擴張原理、空泡鏡像模型和航行體動力學方程,推導了恒定推力作用和勻變速情況下超空泡形態的解析解,探索了記憶效應對空泡形態的影響,獲得了空化數、空化器尺寸和加速度對記憶效應的影響規律。周景軍等[3-4]基于CFD方法研究了超空泡航行體加速和操舵過程中的流體動力特性,獲得了航行體空化器、尾舵操舵、通氣量、重力效應以及航行體攻角等參數變化過程中的空泡形態及航行體瞬態流體動力特性變化規律。李杰等[5]通過建立超空泡航行體自由運動的非定常流場求解數值模擬方法,對尾拍流體動力的形成和演化規律進行了研究,發現尾拍力與空泡壁面形狀、航行速度及尾部沾濕面積等密切相關。陳誠[6]通過建立考慮空泡延遲效應的超空泡航行器動力學模型,對航行器的尾部位置力特性進行了仿真研究,得到了空化器舵角導致流體動力出現附加攻角的規律。王威、何廣華等[7-8]采用數值方法分別研究了垂直方向和水平方向的周期性來流擾動對通氣超空泡航行體流體動力特性的影響,分析了周期性來流的波長和波幅對空泡形態和流體動力的影響規律。羅凱等[9]通過建立定量求解通氣超空化的三維數值模型,研究了通氣量變化對空泡形態的影響,獲得了通氣量突變會使得空泡特征呈現不同程度時間遲滯特性的規律。

綜上所述,公開的文獻主要集中在考慮空泡延遲效應對動力學建模和流體動力特性的研究,缺乏對非定常流體動力延遲特性形成機理的研究。目前,課題組[10]已通過搭建連續變攻角測力試驗平臺,定量分析了航行速度、擺動頻率及空化器預置舵角等因素對超空泡航行體連續變攻角過程中延遲時間特性的影響。在此研究基礎上,本文通過建立連續變攻角運動數值模型,結合對空泡與流體動力變化的分析,進一步揭示不同影響因素下非定常流體動力延遲特性的形成機理,為超空泡航行體動力學建模和機動航行控制技術提供參考。

1 計算模型與方法

1.1 連續變攻角運動模型建立

圖1給出了連續變攻角運動模型坐標系示意圖。設置地面坐標系為Oxyz,初始時刻航行體質心位置為地面坐標系原點,且xOy平面與航行體的對稱面重合,航行體中軸線與x軸重合,空化器指向x軸負方向。設置航行體的轉軸中心Or位于航行體質心O后460 mm處(水洞試驗中尾支撐轉動位置),旋轉軸正方向垂直于對稱面向外,建模過程中不考慮重力的影響。

圖1 坐標系定義

航行體的擺動規律為

式中:θ為擺動角度;θmax為擺動的最大角度,本文中θmax=2°;f為擺動頻率;ωz為角速度。

1.2 控制方程

本文通過數值模擬和水洞試驗驗證相結合的方法開展研究。受制于試驗條件,水洞中的水速最多只能達到12 m/s,形成自然空化能力較弱,因此采用主動通氣實現超空泡流動。借鑒目前通氣超空泡數值模擬研究經驗[11],采用流體體積函數(VOF)多相流模型進行水、汽兩相流動問題模擬研究。

1) 連續性方程

(3)

2) 動量方程

(4)

3) 湍流方程

本文選取SSTk-ω湍流模型,它在更廣泛的流動領域中具有更好的魯棒性和求解精度。

式中:k為湍流動能;ω為湍流耗散率;Γk和Γω分別為k和ω的有效擴散系數;Yk和Yω分別為湍流導致的k和ω的耗散;Dω為交叉擴散項;Gk,Gω為湍流產生項。

1.3 計算模型與邊界條件設置

航行體模型如圖2所示,主要由圓盤空化器、圓錐段、圓柱段以及尾噴管組成。

圖2 航行體模型示意圖

航行體空化器直徑D1=19 mm,航行體圓柱段直徑D2=54 mm,尾噴管直徑D3=32 mm,全長L=432 mm,錐段長L1=181 mm,航行體質心位置距離空化器前端L2=200 mm。

計算域及邊界條件設置如圖3所示,計算域為長方體,考慮流域邊界效應對空泡發展的影響,結合作者對壁面效應的研究經驗[12],計算域尺寸長、寬、高分別為34D2,10D2,20D2,航行體放置于對稱面上,空化器頭部距離計算域左側為1Lc(Lc為空泡最大長度),尾噴管端面距離計算域右側為2Lc。計算域左、上、下面設置為速度入口,前、后面以及航行體設置為壁面,右側面設置為壓力出口。

圖3 計算域及邊界條件

本文采用HEXPRESS軟件進行結構化網格劃分,重點對空化器、航行體肩部及尾流空泡區的網格進行加密處理。通過研究3套不同單元數量的網格(252萬,110萬,42萬)對計算精度和效率的影響,發現網格數量由42萬增加到110萬,最大空泡長度增加5.25%,最大空泡直徑增大7.67%,而當網格數量由110萬增加到252萬,最大空泡長度增加0.53%,最大空泡直徑增加0.62%,綜合分析,最終確定網格數量為110萬,網格劃分細節如圖4所示。

圖4 網格劃分細節

1.4 計算方法

本文中超空泡航行體做受迫連續變攻角運動,模型不僅具有相對于水速的向前速度,還同時有轉動速度。采用動計算域法[13]實現航行體的位置和姿態更新,用UDF控制超空泡航行體按照給定的正弦規律進行擺動,同時將計算獲得的航行體受力結果輸出至文件。

超空泡航行體流場與運動耦合仿真基于Fluent軟件實現,壓力-速度耦合算法選用PISO算法,壓力項離散格式選用PRESTO格式,動量、體積分數、湍流動能和湍流耗散率的離散選用一階迎風格式,計算時間步長保持在2.0×10-4s。

計算流程:0～1 s內,航行體保持靜止,使空泡充分發展;1～3 s內航行體連續往復擺動若干周期,然后航行體停止運動;最后0.5 s空泡恢復原狀態。

2 計算結果與分析

2.1 連續變攻角測力水洞試驗系統

依托西北工業大學高速水洞完成連續變攻角測力試驗平臺的搭建。試驗系統如圖5所示,由攻角調節與測量系統、通氣流量控制與測量系統、流體動力測量系統、壓力測量系統、圖像采集系統及模型主體組成。模型以尾支撐方式安裝于水洞工作段內的電動轉臺上,試驗開始前,使模型保持在0°攻角位置上。通過LABVIEW編寫綜合控制軟件,由單片機控制步進電機,實現模型按設定規律連續轉動。通氣管路、測壓管路和天平導線從導流罩內孔引出至洞體外,壓縮空氣經流量控制器接至模型通氣管。

圖5 試驗系統構成

2.2 數值模型驗證

1) 試驗結果重復性驗證

為減小試驗中通氣流量和水洞水速脈動變化帶來的影響,確保結果準確性,試驗工況均做2～3組重復試驗。圖6給出水速8.5 m/s,頻率1.5 Hz,預置舵角0°工況時,模型擺動方向受力的數據處理前后對比情況。由圖6可以看出,2組數據的重合度較好,試驗的重復性良好,可以說明試驗結果具有較高的可靠性。

圖6 試驗重復性結果

2) 數值結果與試驗結果對比驗證

圖7給出了通氣流量Qg為80 L/min,空化數為0.05時的試驗與仿真計算結果對比。由圖7可以看出,參照試驗空泡結果,仿真空泡最大直徑相對誤差為5.3%,空泡最大長度相對誤差為7.3%,說明建立的數值模型基本合理。

圖7 水洞試驗與數值仿真結果對比

圖8給出了連續變攻角運動過程中擺動方向所受流體動力相對于運動相位之間的延遲效應示意圖。計算各個運動周期T內攻角曲線峰值與流體動力峰值之間的時間差Δt,取多個周期內時間差值Δti求均值作為該工況下的流體動力延遲時間t。為了方便對比,對延遲時間t進行無量綱化處理,定義為延遲系數Cyc,其表達式為

圖8 延遲時間示意圖

(7)

式中:Tck為無量綱參數;Tck=L/v,其中,L為航行體總長;v為參考速度。

圖9給出了來流速度v分別為7.0,7.5,8.0,8.5 m/s、擺動頻率f均為1.5 Hz、保持相同通氣量80 L/min、預置舵角0°等不同工況水洞試驗獲得的延遲系數Cyc與數值仿真的對比結果。由圖9可以看出,仿真結果與試驗結果變化規律一致,延遲系數隨著速度增加而增大。數值計算得到的延遲系數與試驗結果的最大相對誤差為10.4%,進一步證明建立的數值模型具有較好的準確性。

圖9 延遲系數Cyc隨速度v的變化規律

2.3 速度影響分析

圖10給出了水速8.5 m/s,空化數0.05,改變攻角方向,獲得攻角范圍為-2°～2°,每0.2°一個步長攻角下的航行體升力系數變化曲線。由圖10可以看出,在自由攻角以內,航行體所受升力系數始終為零,在攻角αfr=±0.8°,航行體圓柱段處于臨界沾濕狀態;隨著圓柱段沾濕,航行體所受升力系數隨攻角增大而增大,當到達臨界攻角αcr=±1.7°時,圓錐段開始沾濕并在肩部發展形成附體空泡,使得圓柱段沾濕面積減小,圓錐段沾濕面積開始增大,總沾濕面積繼續增大,但升力系數增大的斜率減小。總體而言,升力系數隨攻角增加而增大,沒有呈現遲滯效應。因此,定常流體動力特性不能用來預報航行體在非定常運動時受到的流體動力。

圖10 升力系數隨攻角變化曲線

圖11給出了水速為8.5 m/s,空化數0.05,擺動頻率1.5 Hz,擺動攻角范圍-2°～2°,3個擺動周期內航行體升力系數變化曲線。由圖11可以看出,隨著擺動攻角的周期性變化,航行體的升力系數呈現周期性變化規律。在一個擺動周期內,擺動攻角達到波峰或波谷(如圖中b時刻)時,對應此攻角下的升力系數未達到最大值,而在滯后于該時刻出現最大峰值(如圖中c時刻),說明航行體在連續變攻角運動過程中,流體動力具有延遲特性。

圖11 升力系數隨擺動攻角的變化規律

為了分析升力系數曲線的變化細節,提取了航行體連續變攻角運動下第一個周期內典型時刻的空泡特征。圖12給出了圖11中7個特征時刻下航行體對應的空泡形態。

圖12 一個周期內典型時刻空泡形態

由圖12可以看出,航行體在a時刻之前上表面處于臨界沾濕狀態,升力系數呈現出隨負向攻角增加而小幅度增大的趨勢。a時刻為升力系數轉折點,此時對應的航行攻角為-1.88°,結合圖10可知,在定攻角航行狀態下,航行體自由攻角αfr為±0.8°,說明空泡形態發展相對于攻角改變存在滯后性,這使得流體動力特性呈現出差異性;在b時刻,航行攻角達到最大幅值時(α=-2°),航行體肩部開始出現沾濕,此時升力系數未達到最大值,隨后航行體開始回擺,此時迎流面肩部的附體空泡開始形成,當法向力達到極值時,航行體回擺至攻角為α=-1.93°處,見c時刻,升力系數的峰值滯后于攻角的峰值,升力系數呈現遲滯特性,遲滯的相位差為0.043°。在d時刻,航行體升力系數出現拐點,分析認為,這是由于航行體圓柱段大部分被附體空泡包裹,整個航行體尾部沾濕面積的變化率降低;在e時刻,航行體沾濕面積達到最小,見該時刻空泡形態;在t=T/2時刻后,升力系數迅速增加,可知沾濕面積發生了快速的變化,這是因為形成附體空泡向錐段方向發展并與之融合,圓柱段尾端的空泡未能及時發展,空泡發生潰滅,從而使得沾濕面積增大,空泡融合發展如圖13所示;在f時刻,航行體圓柱段附體空泡與主空泡發展融合,在閉合過程中產生壓力峰,上下表面壓力的不對稱分布形成了作用于航行體的升力,隨后圓柱段上的空泡繼續發展,沾濕面逐漸消失,并在g時刻重新包裹整個航行體。

圖13 附體空泡融合、發展示意圖

表1統計了不同速度下航行體升力系數出現峰/谷值(主峰)和脈沖峰值(副峰)時與攻角的對應關系,其中,Δα為相位角差值,主峰和副峰如圖11中標志所示。由表1可以看出,隨著速度增大,升力系數主峰值相位均滯后于擺動攻角,且相位角差值增加,升力系數主峰值呈減小趨勢,延遲效應加劇;隨著速度增加,升力系數副峰值對應的相位角提前,升力系數主峰值與副峰值比值增大,在v=8.5 m/s時,高達2倍。

表1 升力系數與攻角對應關系

2.4 擺動頻率影響分析

圖14給出了擺動頻率分別為f1=1.0 Hz,f2=1.5 Hz,f3=2.0 Hz,速度v=8.5 m/s,空化數σ=0.05條件下,第一個T/2周期內航行體從α=-1°(↓)到向上α=-1°(↑)的運動過程中,相同攻角狀態時的空泡形態對比。

由圖14可以看出,從α=-1°(↓)擺動到α=-2°(→)過程中,3種擺動頻率下的空泡形態發展趨勢基本一致,擺動頻率的增加使得航行體自由攻角范圍增大,如圖14所示。α=-1.8°(↓)時的空泡形態,此時僅f1工況的航行體圓柱段出現沾濕;當攻角達到最大幅值攻角(α=-2°)時,沾濕面積有S(f1)>S(f2)>S(f3),且航行體錐柱結合位置開始出現沾濕。隨后回擺運動中,航行體圓柱段形成附體空泡,擺動頻率的增加加劇了附體空泡的發展,使得航行體沾濕面積減小。攻角由-1.8°(↑)向-1.0°(↑)連續變化過程中,空泡變化呈現明顯的遲滯發展特性,并隨著頻率的增大,遲滯作用減弱。具體體現為:在f1工況下,隨著攻角幅值減小,附體空泡在軸向擴張變化率較小,在徑向基本未發展變化,使得沾濕面積變化率較小,而其他2種工況下,隨著攻角幅值減小,沾濕面積呈現較明顯減小。

表2統計了第一個擺動周期內不同擺動頻率下升力系數與攻角的對應關系。圖15給出了3種擺動頻率下升力特性隨擺動攻角的變化曲線。結合表2和圖15可以看出,在同一擺動頻率周期內,航行體升力系數呈周期性變化,隨著擺動頻率增大,與擺動攻角的相位差Δα減小,升力系數的遲滯特性減小,升力系數的主峰峰值增大。隨著擺動頻率增加,升力系數副峰在回擺過程中逐步向后移動,其副峰峰值減小。從第一個周期來看,在f1工況下,升力系數副峰峰值出現在航行體負向恢復至零攻角之前,所處攻角為α=-0.53°;在f2工況下,升力系數副峰峰值出現在正向擺動至最大攻角過程中,所處攻角為α=0.19°;在f3工況下,未出現明顯的升力系數脈動,但在航行體在正向擺動至最大攻角過程中仍存在一個較小的脈動峰值,所處攻角為α=0.90°。分析認為,擺動頻率增加加劇了對空泡流場的擾動[11],使得空泡內外壓差發生變化,圓柱段附體空泡發展過程不同,f1和f2工況下圓柱段附體空泡形成了孤立氣泡,最終發生潰滅而產生較大的沖擊載荷[14],而f3工況的附體空泡發展融合于主空泡,擺動頻率的增大促進了空泡的融合。

表2 第一個周期內升力系數與攻角對應關系

圖15 升力系數隨擺動攻角變化規律

2.5 空化器預置舵角影響分析

圖16給出了空化器預置舵角β分別為10°和20°時,航行體0°攻角、空化數為0.05工況下的空泡形態。

圖16 航行體0°攻角狀態下的空泡形態

由圖16可以看出,空化器預置舵角為10°時,空泡軸線向上傾斜,航行體圓柱段仍包裹于空泡內。空化器預置舵角為20°時,空泡軸線上傾狀態更加明顯,航行體圓柱段出現沾濕(見圖中紅色區域),這是由于空化器預置舵角的存在,導致航行體周圍的空化流場呈上下不對稱分布,同時空化器法線與來流速度的夾角Δδ,相當于產生了附加攻角,空化器產生的升力會造成空泡軸線偏斜[15]。

圖17給出了空化器預置舵角β為10°和20°,擺動頻率f=1.5 Hz,航行體連續擺動過程中升力系數隨擺動攻角的變化曲線。表3統計了擺動周期內的升力系數峰值與攻角對應關系。可以看出,首次負攻角擺動過程中,航行體受到正向的法向力沖擊,量值與舵角呈正相關,隨著預置舵角增大,升力沖擊峰值增大,且出現的時間提前;隨后的往復回擺運動過程中,升力系數為恒定值時,航行體處于空泡完全包裹狀態,此階段航行體受力與航行攻角的方向無關,表現為穩定的低頭升力或抬頭升力,這與袁緒龍等[16]獲得的預置舵角能夠產生穩定的抬頭升力從而實現彈道轉平的結論一致。在航行體往復回擺過程中,空化器預置舵角的存在,使得流體動力延遲特性在經歷不同正/負攻角時呈現不同特性。正攻角運動階段中,隨著預置舵角增大,延遲效應減弱,甚至出現相位提前,升力峰值減小;負攻角運動階段中,預置舵角導致延遲效應加劇,升力峰值增大。

表3 升力系數峰值與攻角對應關系

圖17 升力系數隨擺動攻角的變化曲線

為了進一步分析預置舵角對航行體流體動力延遲特性影響的形成原因,圖18給出了圖17中特征點1,2前后的空泡特征。由圖18分析特征點1的形成原因。在t=1.030 s時刻,航行體只有空化器沾濕,受到由預置舵角形成的恒定負向法向力,并隨著航行攻角的減小而減小。隨著航行體下擺運動,發生空泡上飄現象,加之空泡的改變相對于航行體運動的滯后性,導致了航行體尾部背流面開始沾濕,航行體受到了較大的正向升力,見t=1.045 s時刻。隨著航行體負向攻角的繼續增大,空泡上飄現象逐漸被抑制,空泡軸線在流體作用下開始下偏,見t=1.060 s時刻。特征點2的形成原因分析:當航行體進入下一個運動周期時,空泡延遲效應產生的空泡輪廓彎曲導致圓柱段下方出現沾濕,且沾濕面積逐漸擴大,見t=1.701 s時刻。在t=1.705 s時,航行體圓柱段下側沾濕面積達到最大后迅速降低造成了法向力曲線上的尖點,擺動到t=1.730 s時航行體被空泡完全包裹。