莫慧黠,黨建軍,羅 凱,黃 闖,黃 標(biāo)
(西北工業(yè)大學(xué) 航海學(xué)院,陜西 西安,710072)
水下航行器附加質(zhì)量數(shù)值計算方法
莫慧黠,黨建軍,羅 凱,黃 闖,黃 標(biāo)
(西北工業(yè)大學(xué) 航海學(xué)院,陜西 西安,710072)
附加質(zhì)量力是水下航行器在非定常運動中所受到的流體慣性力,其實驗或數(shù)值計算過程均涉及非定常運動,結(jié)果獲取難度大。為了進一步提高水下航行器附加質(zhì)量參數(shù)的計算精度和效率,基于動參考系的計算思想,耦合3D N-S方程和k-epsilon 湍流模型,建立了水下航行器非定常運動的流場計算模型; 對橢球、圓柱標(biāo)準(zhǔn)模型的慣性力特性以及水下航行器的全部附加質(zhì)量特性進行了數(shù)值仿真。結(jié)果表明,橢球、圓柱附加質(zhì)量的數(shù)值計算結(jié)果精確,且水下航行器附加質(zhì)量的計算誤差不超過10%。文中所提方法將有助于水下航行器總體設(shè)計和水動力計算。
水下航行器; 附加質(zhì)量; 動參考系; 數(shù)值仿真
附加質(zhì)量力是水下航行器在非定常運動中所受到的流體慣性力,水下航行器非定常運動相對于定常運動,周圍流體作用于水下航行器的力會產(chǎn)生一個增量。這一部分增量在人為設(shè)定中可以擁有慣性度量的量綱,被稱為附加質(zhì)量力。附加質(zhì)量對于水下航行器的運動穩(wěn)定控制有著重要影響。準(zhǔn)確求解附加質(zhì)量對于水下航行器運動精確控制具有重要意義。
中船重工 705所周景軍等[1]基于相對運動,采用添加動量源項的方法,求解了附加質(zhì)量; 上海交大的姚保太[2]、傅慧萍[3]等人通過求解 3D N-S方程,采用網(wǎng)格重構(gòu),進行了簡單無附體附加質(zhì)量求解; 西北工業(yè)大學(xué)王鵬等[4-5]對復(fù)雜航行器的位置力和阻尼力進行了計算,并未對慣性力進行研究; 而簡單的球、橢球、圓柱等形狀物體可以用基于勢流理論的理想流體附加質(zhì)量數(shù)值計算,獲得理論精確解[6-9]; 此外基于物體外形,采用細(xì)長體假設(shè)及切片理論,使用修正系數(shù)考慮3D效應(yīng)的工程計算方法,也可以計算得到附加質(zhì)量[10-12]; 試驗通過拖曳旋轉(zhuǎn)水池試驗的方法來獲取附加質(zhì)量,但除平動附加慣性力比較準(zhǔn)確之外,附加慣性矩和附加靜矩誤差較大[13]。
因此,為了進一步提高水下航行器附加質(zhì)量參數(shù)的計算精度和效率,基于動參考系的思想,耦合3D N-S方程和k-epsilon湍流模型,建立了水下航行器非定常運動的流場計算模型,通過讓外壁和雷體一起在流域中運動,模擬幾何體的平動與旋轉(zhuǎn),得到了一種新的附加質(zhì)量計算流體力學(xué)(computational fluid dynamics,CFD)數(shù)值解法。無需添加動量源項,不采用動網(wǎng)格以及網(wǎng)格更新技術(shù),僅僅使用一套網(wǎng)格,就能求出復(fù)雜外形航行器的全部附加質(zhì)量,并且經(jīng)過橢球、圓柱以及一個標(biāo)準(zhǔn)的復(fù)雜外形水下航行器進行了檢驗,結(jié)果準(zhǔn)確,計算效率高。
附加質(zhì)量一般取決于流場中物體形狀及運動方向,與運動速度無關(guān)[1]。附加慣性力的影響只有通過非定常運動才能有所表現(xiàn),一般求解思路如下。
水下航行器加速直線運動、變角速度轉(zhuǎn)動等非定常運動過程的數(shù)值仿真最直接的方法是采用動網(wǎng)格技術(shù),但是對于帶有泵噴推進器、對轉(zhuǎn)槳的水下航行器,采用動網(wǎng)格技術(shù)由于涉及到網(wǎng)格重生,網(wǎng)格數(shù)量大,同時質(zhì)量難以保證,從而嚴(yán)重影響求解精度[1],或者通過在動量方程中添加動量源項的方法保證了整個流場壓力分布的真實性,但是在計算附加慣性矩時,網(wǎng)格需要重新劃分,并且旋轉(zhuǎn)運動時,需要編寫復(fù)雜的UDF動量源項,在計算旋轉(zhuǎn)過程中,航行器的受力復(fù)雜,受網(wǎng)格精度的影響很大。
而采用動參考系法,將邊界與研究對象固定為體坐標(biāo)系,流域與地面坐標(biāo)系綁定,邊界與研究對象共同運動,相對于研究對象速度不為零,邊界條件采用速度進口與壓力出口,需要編寫的UDF進行二次開發(fā)。地面系不賦予任何運動特征,設(shè)置默認(rèn),速度為零,直接賦予研究對象所在體系相應(yīng)的運動初始速度與加速度即可。
關(guān)于其方法的基本內(nèi)涵方程如下
式中:vr是流域相對地面坐標(biāo)的運動速度;ur是研究對象相對流域的運動速度;v是研究對象相對于地面坐標(biāo)的運動速度;ur1是研究對象相對邊界的運動速度;ur2是邊界相對于流域的速度;ω1是研究對象相對流域的轉(zhuǎn)速;ω2是流域相對于地面坐標(biāo)系的轉(zhuǎn)速,ω是研究對象相對于地面坐標(biāo)系的轉(zhuǎn)速。
該方法的本質(zhì)是運動的相對性,與添加動量源項的方法不同的是,該方法是直接將速度特性賦予研究對象與邊界。無需進行第2套網(wǎng)格劃分,賦予 UDF二次開發(fā)的平動與轉(zhuǎn)動,將不同網(wǎng)格的影響徹底排除。
文中采用雷諾時均 N-S方程,湍流模型采用Realizablek-ε模型。不可壓縮流動控制方程主要包括連續(xù)性方程和N-S方程。
1) 連續(xù)性方程
2) 運動方程(N-S方程)
式中:U為速度矢量;ρ,p,g,μ分別為密度、壓強、重力加速度和動力粘度。
劉丹[13]、劉成剛[14]、馬燁[15]等人的論文中提到關(guān)于Realizablek-ε是為了應(yīng)對標(biāo)準(zhǔn)k-ε模型對時均應(yīng)變率特別大時,有可能導(dǎo)致負(fù)的正應(yīng)力出現(xiàn),為使流動符合湍流的物理定律,對正應(yīng)力進行某種數(shù)學(xué)約束的一個改進模型。在 Realizablek-ε模型中,關(guān)于k和ε的輸運方程如下。
關(guān)于湍流強度k方程
關(guān)于湍流耗散率ε方程
式中:μt為湍動粘度;μ為流體的時均速度;kσ,εσ分別為k-ε方程的湍流能量普朗特數(shù);C1,C2為經(jīng)驗常數(shù);E為時均應(yīng)變率;ν為運動粘度;xi,xj為各方向距離。
建立與控制方程相應(yīng)的離散方程,文中采用工程上應(yīng)用廣泛的壓力耦合方程組的半隱式方法——SIMPLE算法,通過構(gòu)造壓力修正方程以及速度修正方程,直至求得收斂解為止。
網(wǎng)格劃分是數(shù)值方法的基礎(chǔ),以某典型外形水下航行器為例,對其進行網(wǎng)格劃分。該航行器由于附體形狀比較復(fù)雜,因此網(wǎng)格劃分是工作重點。此次網(wǎng)格劃分采用 ICEM軟件,流域為圓柱體,流域長度為8l,l為產(chǎn)品的長度,直徑為20d,d為產(chǎn)品圓柱段的橫截面直徑,以確保其流域足夠大。航行器附近流域為加密區(qū)域,加密區(qū)邊界層起始厚度為0.1 mm,比例為1.2倍,共有5層,最終得到了一套完整的高質(zhì)量網(wǎng)格,網(wǎng)格質(zhì)量均在 0.5以上,可以滿足計算要求。網(wǎng)格如圖 1和圖2所示。
圖1 流域網(wǎng)格Fig. 1 Mesh of flow field
圖2 航行器網(wǎng)格Fig. 2 Mesh of an undersea vehicle
λ11,λ22,λ33,λ62計算邊界條件設(shè)置: 3 個平移方向的附加質(zhì)量通過3個方向加速運動與定常勻速運動的慣性力作差獲得,而λ62是由于平動產(chǎn)生的矩,同樣通過加速運動與定常勻速運動的慣性靜矩作差可以得到。在動參考系方法中,無需添加動量源項,如圖1,速度進口不給予設(shè)置,速度設(shè)為零,其余默認(rèn),尾部的邊界條件設(shè)置為壓力出口,參考壓力設(shè)置為零即可,直接在動坐標(biāo)系域中賦予研究對象速度特性與加速度特性,導(dǎo)入UDF,便能獲得真實的平動運動壓力場分布。
λ44,λ55,λ66,λ35計算邊界條件設(shè)置: 3 個旋轉(zhuǎn)方向的附加質(zhì)量通過3個方向繞質(zhì)心的加速旋轉(zhuǎn)運動與定常勻速旋轉(zhuǎn)的附加慣性矩作差獲得,λ35是旋轉(zhuǎn)產(chǎn)生的力,同樣通過加速運動與定常勻速運動的慣性力作差可以得到。在動參考系中,無需添加動量源項,設(shè)置見圖 1,直接在動坐標(biāo)系域中賦予研究對象旋轉(zhuǎn)角速度與角加速度特性即可,便能獲得真實的旋轉(zhuǎn)運動壓力場分布。
若用相對運動法,在計算非定常運動時,添加動量源項才能保證流場的真實性,并且在計算λ62和λ35時還要進行特殊處理,比如為了避免流動發(fā)生分離,采用理想流體計算; 而動網(wǎng)格技術(shù)需要網(wǎng)格重生,計算時間很長。該方法避開了 2種方法的不足,較以往方法有所改進。
文中全部采用了結(jié)構(gòu)化網(wǎng)格,一是對網(wǎng)格數(shù)量進行了有效控制,二是對網(wǎng)格質(zhì)量進行了確保,從而在效率和結(jié)果上都確保了可行性。使用的數(shù)值方法為了證明其可用性與準(zhǔn)確性,計算模型首先采用了有精確解的橢球。橢球長半軸長500 mm,短半軸長250 mm,網(wǎng)格數(shù)量分別選取20萬、40萬以及60萬,均采用有粘模型進行數(shù)值計算。
分別對 3種數(shù)量的網(wǎng)格進行計算,得到了橢球的阻力系數(shù),對其進行了網(wǎng)格無關(guān)性驗證。最后選取數(shù)量為20萬的橢球網(wǎng)格用于計算附加質(zhì)量。
給出了對應(yīng)尺寸橢球附加質(zhì)量理論解以及數(shù)值計算解的對比見表 1。對比計算結(jié)果,誤差在0.5%以內(nèi),精確度極高。為了驗證該方法的適用性,又使用圓柱進行了驗證。圓柱與橢球為 2個擁有精確理論的簡單模型。
表1 橢球附加質(zhì)量計算結(jié)果Table 1 Calculation results of ellipsoid′s additional mass
由于圓柱的理論計算公式要求圓柱的長細(xì)比較大,選取長細(xì)比L/D=200。將圓柱的旋轉(zhuǎn)中心定在圓柱的形心,此時,圓柱的附加質(zhì)量中不存在λ26,λ35。而λ11,λ44則基本可以忽略。表 2 給出了對應(yīng)尺寸的圓柱附加質(zhì)量理論解以及數(shù)值計算解的對比結(jié)果,可以發(fā)現(xiàn)其理論解與數(shù)值計算結(jié)果相對誤差在2%以內(nèi),
表2 圓柱附加質(zhì)量計算結(jié)果Table 2 Calculation results of circular cylinder′s additional mass
用相同的模型,將其旋轉(zhuǎn)中心定在形心偏左1 m的地方,這樣λ26和λ35就會存在。再次將對應(yīng)尺寸的理論參考值進行求解,再用模型進行驗證,得出結(jié)果如表3所示。
表3 形心偏左1 m時圓柱附加質(zhì)量計算結(jié)果Table 3 Calculation results of circular cylinder′s additional mass when the centroid is 1 m left
由橢球和圓柱計算結(jié)果對比,可以初步驗證該方法的計算結(jié)果是準(zhǔn)確的,進一步驗證了該方法求解附加質(zhì)量的可行性。
對所劃分的計算域網(wǎng)格進行無關(guān)性驗證后,開始流場計算,邊界條件見 1.2.3節(jié)求解λ11時所得到的對應(yīng)航行器對稱面壓力云圖如圖3所示。該圖的本質(zhì)是航行器做勻減速運動,參考表壓力為零,由靜壓曲線圖 4表明,航行器運動使流體有了速度動壓,所以靜壓曲線會出現(xiàn)負(fù)值。從而可知頭部高壓至尾部的壓力恢復(fù),頭部為滯止點,壓力高,計算流體有粘性及粘性損耗,尾部壓力無法恢復(fù)至與頭部相同的壓力,形成壓差阻力。圖3表明,求解λ11的原理與源項法及動網(wǎng)格法相近,均為求解非定常直線運動與定常直線運動受力,作差求解λ11。
圖3 航行器水平運動壓力云圖Fig. 3 Pressure contour of an undersea vehicle in horizontal motion
圖4 航行器水平運動軸向壓力特性曲線Fig. 4 Characteristic curve of axial pressure of an undersea vehicle in horizontal motion
而根據(jù)數(shù)值計算得到的仿真結(jié)果得出,不同的加速度計算附加質(zhì)量有 2%~5%的誤差,為了消除系統(tǒng)誤差,保證結(jié)果的準(zhǔn)確性,選取不同的初速度和加速度,即對多個不同運動工況進行求解,得到多個附加質(zhì)量參數(shù)λ11,文中提到附加質(zhì)量參數(shù)只與形狀運動方向有關(guān),與速度大小無關(guān),附加質(zhì)量參數(shù)λ11進行算術(shù)平均,如表4所示。
求解λ22,λ33以及λ26相關(guān)的壓力云圖如圖5所示,邊界條件見 1.2.3,向z側(cè)移動,賦予航行器z向平移速度,航行器殼體的表面靜壓壓力特性曲線如圖6所示。壓力特性曲線中出現(xiàn)負(fù)值,是因為航行器運動過程中使水有了速度(0.4~0.6 m/s),而參考表壓力設(shè)置為零,水有了動壓,因此在參考表壓力為零的情況下,靜壓出現(xiàn)了負(fù)值。一側(cè)高壓,一側(cè)低壓,流體依舊有粘性及粘性損失,與求解λ11時相同,對多組工況進行求解,得到附加質(zhì)量參數(shù)λ22,λ33以及λ26的算術(shù)平均值(見表 4)。
表4 航行器附加質(zhì)量求解結(jié)果Table 4 Calculation results of undersea vehicle′s additional mass
圖5 航行器側(cè)向運動壓力云圖Fig. 5 Pressure contour of an undersea vehicle in lateral motion
圖6 航行器側(cè)向運動軸向壓力特性曲線Fig. 6 Characteristic curve of axial pressure of an undersea vehicle in lateral motion
求解λ55,λ66以及λ35得到的壓力云圖如圖 7所示,邊界條件見1.2.3節(jié),直接賦予航行器角速度與角加速度即可,航行器殼體表面的靜壓壓力特性曲線如圖8所示。運動為繞航行器浮心位置轉(zhuǎn)動殼體一側(cè)的壓力不相等,航行器頭部受到流體作用的一個低頭力矩,航行器尾部受到一個抬尾力矩,殼體需要克服逆時針的旋轉(zhuǎn)力矩,同求解λ11過程,需要選取不同角速度與角加速度的工況,經(jīng)過求解,得到多個附加質(zhì)量參數(shù)λ55,λ66以及λ35的仿真結(jié)果,為了消除系統(tǒng)的誤差,最終得到各個附加質(zhì)量參數(shù)對應(yīng)的算術(shù)平均值(見表4)。
圖7 航行器繞浮心轉(zhuǎn)動壓力云圖Fig. 7 Pressure contour of an undersea vehicle rotating around center of buoyancy
圖8 航行器繞浮心轉(zhuǎn)動軸向壓力曲線Fig. 8 Characteristic curve of axial pressure of an undersea vehicle rotating around center of buoyancy
對該復(fù)雜外形航行器的附加質(zhì)量用文中方法進行解算得到的結(jié)果中,可得結(jié)論: 在附加質(zhì)量主要的幾個參數(shù),λ22,λ33,λ26,λ55,λ66,λ35結(jié)果都是滿意的,再次驗證了該方法的求解附加質(zhì)量精度相當(dāng)高,并且相比較前人的相對運動法和動網(wǎng)格方法,求解效率有了質(zhì)的變化。
文中基于運動參考系的思想,將邊界與研究對象固定為體坐標(biāo)系,流域與地面坐標(biāo)系綁定,邊界與研究對象共同運動,邊界控制選取與速度無關(guān)的方程,湍流模型離散成與邊界控制相對的方程,得到了一套完整的水下航行器附加質(zhì)量計算方法。
文中基于運動參考系和 k-epsilon湍流模型建立了一種適應(yīng)性強、精度高的水下運動體全部附加質(zhì)量參數(shù)解算方法。采用該方法,計算橢球和圓柱2種標(biāo)準(zhǔn)模型的附加質(zhì)量具有較高的精度,相對誤差不超過 2%; 對于有附體水下航行器等比例模型,文中方法的計算精度與標(biāo)稱值[16]誤差不超過10%。
文中避免了采用相對運動法需要添加動量源項以及不同附加質(zhì)量求解需要重新劃分網(wǎng)格的問題,也避免了動網(wǎng)格需要網(wǎng)格更新而帶來的計算緩慢,網(wǎng)格質(zhì)量難以保證的問題,解決了傳統(tǒng)水下航行器附加質(zhì)量求解難度大、計算效率低的問題。
[1] 周景軍,李育英. 一種水下航行體附加質(zhì)量數(shù)值計算方法[J]. 魚雷技術(shù),2013,21(4): 246-249.Zhou Jing-jun,Li Yu-ying. A Numerical Computation Method of Additional Mass for Underwater Vehicle[J].Torpedo Technology,2013,21(4): 246-249.
[2] 姚保太,康寧,鄭偉奇. 變加速運動圓球附加質(zhì)量和阻力仿真分析[J]. 計算機輔助工程,2014,23(3): 82-87.Yao Bao-tai,Kang Ning,Zheng Wei-qi. Simulation and Analysis on Added Mass and Drag of Varying Accelerated Motion Ball[J]. Computer Aided Engineering,2014,23(3): 82-87.
[3] 傅慧萍,李杰. 附加質(zhì)量CFD計算方法研究[J]. 哈爾濱工程大學(xué)學(xué)報,2011,32(2): 148-152.Fu Hui-ping,Li Jie. Numerical Studies of Added Mass Based on the CFD Method[J]. Journal of Harbin Engineering University,2011,32(2): 148-152.
[4] 王鵬,寧騰飛,杜曉旭,等. 帶復(fù)雜外形附體的AUV流體動力數(shù)值計算[J]. 兵工學(xué)報,2013,34(2): 223-228.Wang Peng,Ning Teng-fei,Du Xiao-xu,et al. Numerical Calculation of Hydrodynamics of AUV with the Complex Shape Appendage[J]. Acta Armamentarii,2013,34(2):223-228.
[5] 王鵬,翟繼瑩,寧騰飛. 帶復(fù)雜外形附體的 AUV 旋臂水池數(shù)值計算[J]. 西北工業(yè)大學(xué)學(xué)報,2013,31(5):764-769.Wang Peng,Zhai Ji-ying,Ning Teng-fei. Calculating Rotary Derivatives of AUV with Complex Shape Appendages[J]. Journal of Northwestern Polytechnical University,2013,31(5): 764-769.
[6] Causin P,Gerbeau J F,Nobile F. Added Mass Effects in the Design of Partitioned Algorithms for Fluid-Structure Problems with Application to Blood Flows[J]. Computer Methods in Applied Mechanics & Engineering,2005,194(42-44): 4506-4527.
[7] Lu L,Yang Y R,Li P,et al. Added Mass,Added Stiffness and Added Damping Coefficients for a Parallel Plate-Type Structure[J]. Applied Mechanics & Materials,2011,66-68: 1738-1742,2011.
[8] Zhu J,Lin Z,Liu Q,et al. Calculation of the Added Mass of a Liquid Tank′s Bulkheads[J]. Journal of Marine Science and Application,2014,13(1): 41-48.
[9] He C,Duan Z,Ou J,et al. Coefficients of Added Mass/damping of a Low-mass-damping-ratio Cylinder Subjected to Vortex-induced Vibration[C]//第十四屆中國海洋(岸)工程學(xué)術(shù)討論會論文集(上冊). 呼和浩特: 中國海洋,2009.
[10] Wang K,Zhang X,Zhang Z Q,et al. Numerical Analysis of Added Mass and Damping of Floating Production,Storage and Offloading System[J]. Acta Mechanica Sinica,2012,28(3): 870-876.
[11] Watts P B,Jensen R L,David M,et al. Vertical Hand Force and Forearm Emg during a High-step Rock-on Climbing Move With and without Added Mass[J].Medicine & Science in Sports & Exercise,2005,37(S5):S122.
[12] 馮雙雙,魏曉娟,孫磊,等. 流場中沿軸向運動圓柱的附加質(zhì)量計算[J]. 計算機輔助工程,2015,24(2): 42-46.Feng Shuang-shuang,Wei Xiao-juan,Sun Lei,et al.Calculation on Added Mass of Cylinder Moving along Its Axis in Flow Field[J]. Computer Aided Engineering,2015,24(2): 42-46.
[13] 劉丹,王曉亮,單雪雄. 平流層飛艇的附加質(zhì)量及其對飛艇運動的影響[J]. 計算機仿真,2006,23(6): 52-56.Liu Dan,Wang Xiao-liang,Shan Xue-xiong. Added Mass to Stratospheric Airship and Its Effect on Motion[J].Computer Simulation,2006,23(6): 52-56.
[14] 劉成剛. 潛艇附加質(zhì)量計算及其水中振動特性研究[D].哈爾濱: 哈爾濱工程大學(xué),2011.
[15] 馬燁,單雪雄. 數(shù)值計算復(fù)雜外形物體附加質(zhì)量的新方法[J]. 計算機仿真,2007,24(5): 75-78.Ma Ye,Shan Xue-xiong. A New Numerical Computation Method for Added Masses of Complicated Object[J].Computer Simulation,2007,24(5):75-78.
[16] 嚴(yán)衛(wèi)生. 魚雷航行力學(xué)[M]. 西安: 西北工業(yè)大學(xué)出版社,2005: 238.
Numerical Calculation of Additional Mass for Undersea Vehicle
MO Hui-xia,DANG Jian-jun,LUO Kai,HUANG Chuang,HUANG Biao
(School of Marine Science and Technology,Northwestern Polytechnical University,Xi′an 710072,China)
Additional mass force on undersea vehicle is a fluid inertia force due to unsteady motion,which is difficult to obtain from experiment or numerical calculation. In this paper,an unsteady flow field calculation model of an undersea vehicle is presented in the moving coordinates system by making use of the 3D N-S equation and the k-epsilon turbulence model. The proposed model is verified by inertia force solutions of ellipsoid and circular cylinder models. Then,based on the present numerical method,the additional mass of an undersea vehicle is calculated. Simulation results show that the relative errors of the obtained additional mass are less than 10%. This numerical calculation method of additional mass may be helpful to the overall design and hydrodynamic calculation of an undersea vehicle.
undersea vehicle; additional mass; moving coordinates system; numerical simulation
TJ630.1; O351.2
A
2096-3920(2017)03-0250-06
莫慧黠,黨建軍,羅凱,等. 水下航行器附加質(zhì)量數(shù)值計算方法[J]. 水下無人系統(tǒng)學(xué)報,2017,25(3): 250-255.
10.11993/j.issn.2096-3920.2017.03.006
2017-03-31;
2017-05-25.
國家自然科學(xué)基金項目(51579209、51409215、51679202).
莫慧黠(1992-),男,在讀碩士,主要研究方向為水下航行器流體動力仿真與動力機械仿真.
(責(zé)任編輯: 許 妍)