姚保太+康寧+鄭偉奇
作者簡介: 姚保太(1973—),男,河南獲嘉人,高級工程師,研究方向為機械制造及仿真計算,(Email)yaobaotai@163.com
通信作者: 康寧(1963—),女,遼寧金州人,副教授,博士,研究領(lǐng)域方向為流體仿真計算,(Email)kangning@buaa.edu.cn0引言
附加質(zhì)量和阻力計算是水下物體運動分析的關(guān)鍵問題.[1]以往對附加質(zhì)量的計算主要是根據(jù)勢流理論[14],如在無界無黏流體中球體等簡單形狀物體的附加質(zhì)量存在解析解,復(fù)雜形狀物體的附加質(zhì)量通過估算法或三維面元法求得.近年來,隨著計算流體力學(xué)的發(fā)展,開始用數(shù)值模擬的方法計算物體的附加質(zhì)量,發(fā)展基于動網(wǎng)格技術(shù)計算附加質(zhì)量的新方法.[5]數(shù)值方法有更廣泛的應(yīng)用空間,不僅可以計算理想流體中的附加質(zhì)量[56],還可以計算黏性流體中的附加質(zhì)量[710].
人們已經(jīng)熟知無界無黏流體中物體加速運動時的附加質(zhì)量只與物體形狀及流體密度有關(guān),與運動加速度無關(guān),目前還缺少系統(tǒng)地研究在有界區(qū)域中流體黏性和加速度對做變加速運動物體附加質(zhì)量的影響.由于工程中實際問題的復(fù)雜性,比如水下導(dǎo)彈發(fā)射問題等,因此有必要研究在有界區(qū)域黏性流體中物體做變加速運動的附加質(zhì)量問題.
本文采用數(shù)值模擬方法和動網(wǎng)格技術(shù),計算管體域中黏性不可壓流體中圓球做變加速運動時的附加質(zhì)量、阻力和阻力因數(shù),深入分析和研究在管體域中圓球變加速運動對附加質(zhì)量、阻力和阻力因數(shù)的影響.
1計算模型和計算方法
1.1控制方程組
水中運動物體滿足非定常不可壓NS方程,
連續(xù)方程: vixi=0(1)動量方程:vit+vjvixj=Fi-1ρpxi+μρxjvixj(2)式中:vi和vj為速度;F為質(zhì)量力;p為壓力;ρ為流體密度;μ為流體動力黏性系數(shù).
1.2計算域、邊界條件和網(wǎng)格劃分
1.2.1管體域
管體域由管體、喇叭口和水域組成,圓球直徑為0.254 m,圓球在管體域中運動前的計算域見圖1.通海水口邊界條件為壓強進(jìn)口,出口為壓強出口,壓強均為57 m水深處壓強574 826.8 Pa;根據(jù)圓球平移運動特點,動網(wǎng)格模型選擇動態(tài)網(wǎng)格層變方法,圓球表面為無滑移運動壁面條件;管體和喇叭口為無滑移靜止壁面條件,水域邊界為滑移靜止壁面條件.圓球在管體域中運動前、后網(wǎng)格分布見圖2和3.采用六面體結(jié)構(gòu)網(wǎng)格,網(wǎng)格數(shù)約為58萬個.
圖 1管體域中圓球運動前的計算域
Fig.1Computation domain in tube domain before ball
movement
圖 2管體域中圓球運動前的網(wǎng)格分布
Fig.2Mesh distribution of tube domain before ball movement
圖 3管體域中圓球運動后的網(wǎng)格分布
Fig.3Mesh distribution of tube domain after ball movement
1.2.2小計算域
作為對比,對圓球在小計算域中的運動情況進(jìn)行計算.小計算域為圓柱形,取較大尺度,直徑為3.6 m,長度為14 m.縱向?qū)ΨQ面網(wǎng)格分布見圖4.六面體結(jié)構(gòu)網(wǎng)格,網(wǎng)格數(shù)約為139萬個.計算域左側(cè)為壓強入口邊界條件,右側(cè)為壓強出口條件,相對壓強為0;其他邊界和動網(wǎng)格模型同管體域.
圖 4小計算域縱向?qū)ΨQ面的網(wǎng)格分布
Fig.4Mesh distribution in longitudinal symmetry plane of
small computation domain
1.2.3大計算域
為驗證附加質(zhì)量仿真方法的可靠性,對圓球在大計算域中的運動情況進(jìn)行計算.大計算域仍為圓柱形,直徑加大到14.4 m,長度加大到20.86 m.采用六面體結(jié)構(gòu)網(wǎng)格,網(wǎng)格數(shù)約為264萬個.邊界條件和動網(wǎng)格模型同管體域.1.3圓球變速運動參數(shù)
1.3.1先加速后減速
為模擬發(fā)射體在管體內(nèi)直至出管的運動情況,對圓球先加速后減速運動情況進(jìn)行計算.加速階段模擬發(fā)射時在主動力作用下的運動情況,減速階段模擬在慣性力作用下直至出管的運動情況,具體參數(shù)見表1.
表 1圓球先加速后減速的運動參數(shù)
Tab.1Movement parameters when ball is accelerated motion first and then decelerated初速度/
(m/s)加速度/
(m/s2)時間/s末速度/
(m/s)位移/
m加速階段052.0000.32917.1002.814減速階段17.108-5.6690.39014.8976.241總和0.71931.8979.055
1.3.2一直加速
作為對比,對圓球在管體域內(nèi)一直加速運動工況進(jìn)行計算,加速度和初速度與先加速后減速的加速階段相同,位移與先加速后減速的總位移相近,見表2.
表 2圓球一直加速的運動參數(shù)
Tab.2Movement parameters when ball is accelerated motion all the time初速度/
(m/s)加速度/
(m/s2)時間/s末速度/
(m/s)位移/m052.0000.58930.6209.020
1.4求解方法和計算參數(shù)
速度和壓強耦合方法采用SIMPLE方法,壓力修正方程離散格式采用Standard格式,對流項離散格式采用2階迎風(fēng)格式,湍流模型為RNG kε湍流模型,近壁區(qū)域處理采用非平衡壁面函數(shù).
海水密度為1 028 kg/m3,動力黏性系數(shù)為1.394×10-3 kg/(m?s),參考壓強為1標(biāo)準(zhǔn)大氣壓101 325 Pa.計算初始速度為0,初始壓強為57 m水深處壓強574 826.8 Pa.
2仿真方法驗證
2.1阻力
分2種情況對圓球阻力因數(shù)進(jìn)行驗證,分別為水繞靜止圓球運動和圓球在靜水中運動.
圓球?qū)嶒炞枇σ驍?shù)的范圍[3]為0~106,根據(jù)定義可得圓球速度范圍為0~5.339 m/s,因此對最大阻力因數(shù)106和較小阻力因數(shù)104,其分別對應(yīng)最大速度5.339 m/s和較小速度0.053 4 m/s的情況進(jìn)行計算.采用小計算域進(jìn)行計算.
2.1.1靜止圓球
計算域左側(cè)為速度入口邊界條件,速度大小分別為5.339 0和0.053 4 m/s;右側(cè)為壓強出口邊界條件,相對壓強為0.
圓球表面為無滑移靜止壁面邊界條件;除進(jìn)口、出口和圓球表面以外的邊界條件為滑移靜止壁面邊界條件.計算結(jié)果見表3.
表 3圓球靜止時的阻力因數(shù)
Tab.3Drag factor of stationary ball速度/(m/s)實驗值計算值誤差/%0.053 40.400.4164.005.339 00.170.1652.94
由表3可知,在2種速度情況下計算的圓球阻力因數(shù)與實驗結(jié)果的誤差都在4%以內(nèi),說明圓球靜止時阻力計算結(jié)果可靠.
2.1.2運動圓球
計算域左側(cè)為壓強入口條件,右側(cè)為壓強出口條件,相對壓強為0.動網(wǎng)格模型及圓球表面邊界條件均同管體域中的運動圓球,速度大小分別為5.339 0和0.053 4 m/s;其余邊界為滑移靜止壁面條件;初始速度和壓強均為0.計算結(jié)果見表4.
表 4圓球運動時的阻力因數(shù)
Tab.4Drag facctor of moving ball速度/(m/s)實驗值計算值誤差/%0.053 40.400.3717.255.339 00.170.18710.00
由表4可知,在2種速度情況下計算的圓球阻力因數(shù)與實驗結(jié)果的誤差都在10%以內(nèi),說明圓球運動時阻力計算結(jié)果可靠.
2.2附加質(zhì)量
2.2.1附加質(zhì)量計算公式
根據(jù)文獻(xiàn)[2],附加質(zhì)量λ=ρv2iViv2o(3)式中:vi為第i控制體的流體速度;Vi為第i控制體的體積;vo為物體運動速度.
2.2.2圓球附加質(zhì)量的理論解
根據(jù)文獻(xiàn)[1],在無界無黏流體情況下,圓球直線加速運動附加質(zhì)量的理論解λ=23πρR3 (4)式中:R為圓球半徑.
圓球直徑、水密度和動力黏性因數(shù)同前,這時圓球附加質(zhì)量的理論解為4.41 kg.
2.2.3圓球附加質(zhì)量的計算結(jié)果
為驗證附加質(zhì)量仿真方法的可靠性,對小計算域中圓球在無黏流體中做不同加速度的運動進(jìn)行計算,圓球的初速度均為0,加速時間為1.1 s,加速度分別為0.01,0.10,1.00和5.00 m/s2,并與圓球附加質(zhì)量的理論解進(jìn)行對比.所采用的網(wǎng)格、邊界條件、初始條件、計算設(shè)置和動網(wǎng)格模型均與運動圓球阻力驗證時一樣.計算的附加質(zhì)量見圖5.
圖 5小計算域中加速度對圓球附加質(zhì)量的影響
Fig.5Effect of acceleration on added mass of ball in
small computation domain
由圖5可知,時間和加速度越大附加質(zhì)量越大.在加速度較小時(0.01,0.10和1.00 m/s2)計算的附加質(zhì)量與理論解較為接近;在加速度較大時(5.00 m/s2),運動初期計算的附加質(zhì)量與理論解較接近,隨著時間增大其值開始大于理論解.這是由于隨著加速度和時間的增大,運動速度增大,導(dǎo)致計算域邊界效應(yīng)的影響也隨之增大,附加質(zhì)量增加,因而與理論解的差別增大.
為進(jìn)一步研究計算域邊界效應(yīng)的影響,在大計算域中計算加速度為5.00 m/s2的工況,初速度和加速時間均同前文,無黏流體中的附加質(zhì)量的計算結(jié)果見圖6,并與之前小計算域結(jié)果進(jìn)行對比.由圖6可知,除運動初期,無黏流體中小計算域附加質(zhì)量均大于大計算域附加質(zhì)量,大計算域附加質(zhì)量均大于理論解,并且隨著時間增加速度加大,差別也越大.大計算域附加質(zhì)量與理論解較接近,小計算域附加質(zhì)量與理論解差別較大.由圖5和6的結(jié)果可知,加速度越小或計算域越大,邊界效應(yīng)越小,則計算的附加質(zhì)量越接近理論解.圖 6在a=5 m/s2時無黏流體中計算域尺度對圓球
附加質(zhì)量的影響
Fig.6Effect of computation domain scale on added mass of
ball at a=5 m/s2 in inviscid fluid
2.2.4橢球體附加質(zhì)量的實驗驗證
為驗證附加質(zhì)量的理論計算公式(3),將黏性流體中加速運動橢球體附加質(zhì)量的仿真結(jié)果與文獻(xiàn)[6]的實驗結(jié)果進(jìn)行對比.橢球長軸為0.4 m,短軸為0.2 m,與文獻(xiàn)[6]中的完全一樣.計算域為圓柱形,直徑為6 m,長度為14 m.計算時水溫度為5 ℃,密度為1 000 kg/m3,動力黏性系數(shù)為1.518×10-3 kg/(m?s).橢球體為六面體網(wǎng)格,網(wǎng)格數(shù)為106萬個,網(wǎng)格分布見圖7.
圖 7橢球體縱向?qū)ΨQ面的網(wǎng)格分布圖
Fig.7Mesh distribution in longitudinal symmetry
plane of ellipsoid
邊界條件、初始條件和其他設(shè)置均同阻力驗證時的運動圓球,橢球動網(wǎng)格模型同運動圓球.橢球體運動初速度為0,加速度為1 m/s2,加速時間分別為0,0.5和0.9 s.根據(jù)式(3)計算的附加質(zhì)量和與文獻(xiàn)[6]的對比見表5.由表5可知,計算的附加質(zhì)量與實驗結(jié)果比較接近,誤差在9%以內(nèi),說明附加質(zhì)量的計算結(jié)果基本正確.
表 5橢球體的附加質(zhì)量
Tab.5Added mass of ellipsoid加速時間/s實驗值/kg計算值/kg誤差/%0.11.701.742.350.51.701.774.120.91.701.858.82
綜合前文分析結(jié)果可知,附加質(zhì)量的仿真方法較為可靠.3計算結(jié)果及分析
對圓球在管體域黏性流體中先加速后減速和一直加速2種變速運動的附加質(zhì)量、阻力和阻力因數(shù)進(jìn)行仿真計算,并對小計算域黏性流體中圓球先加速后減速運動進(jìn)行仿真計算,同時給出圓球在無界無黏流體中附加質(zhì)量的理論解.
管體域和小計算域中圓球在黏性流體中附加質(zhì)量、阻力和阻力因數(shù)隨位移的變化情況見圖8~10,在管體域中4 m處圓球先加速后減速和一直加速的速度矢量圖和壓強云圖見圖11和12.
圖 8管體域和小計算域中圓球的附加質(zhì)量
Fig.8Added masses of ball in tube domain and small
computation domain
圖 9管體域和小計算域中圓球的阻力
Fig.9Drags of ball in tube domain and small computation
domain
圖 10管體域和小計算域中圓球的阻力因數(shù)
Fig.10Drag factors of ball in tube domain and small
computation domain(a)速度矢量圖,m/s
(b)壓強云圖,Pa
圖 11在管體域中4 m處先加速后減速圓球的計算結(jié)果
Fig.11Calculation results at 4 m in tube domain while ball is accelerated first and then decelerated
(a)速度矢量圖,m/s
(b)壓強云圖,Pa
圖 12在管體域中4 m處一直加速圓球的計算結(jié)果
Fig.12Calculation results at 4 m in tube domain while ball is accelerated all the time
由圖8可知,在運動開始時,小計算域中附加質(zhì)量仿真值與理論解4.41 kg吻合較好;隨著運動進(jìn)行仿真值開始大于理論解,最后達(dá)到20.8 kg.管體域中附加質(zhì)量均大于小計算域中附加質(zhì)量,并隨著運動進(jìn)行差別逐漸增大,在出管體(7 m)后達(dá)到最大值,然后開始變小.先加速后減速運動的最大值為58.6 kg,最終降為48.1 kg;一直加速運動的最大值為44.7 kg,最終降為33.9 kg.由于管體域尺寸小于小計算域,邊界效應(yīng)較大,因而附加質(zhì)量較大,出喇叭口后過流斷面變大,邊界效應(yīng)減小,導(dǎo)致附加質(zhì)量變小.管體域中圓球先加速后減速和一直加速的附加質(zhì)量在2.8 m之前完全一樣,因為該階段2種工況運動狀態(tài)完全一樣;在2.8 m之后,一直加速運動附加質(zhì)量仍保持連續(xù)增加態(tài)勢,而先加速后減速運動附加質(zhì)量在2.8 m處出現(xiàn)不連續(xù)現(xiàn)象,開始大于一直加速附加質(zhì)量.先加速后減速運動由加速狀態(tài)變?yōu)闇p速狀態(tài)時,運動狀態(tài)的改變遠(yuǎn)大于一直加速運動狀態(tài).由圖11和12可看出,在4 m處,先加速后減速圓球周圍速度雖然小于一直加速的速度,但尾跡流影響范圍明顯較大,在圓球尾部較遠(yuǎn)處還卷起較大尾渦,因而導(dǎo)致附加質(zhì)量非但沒有變小,反而增大.小計算域附加質(zhì)量在加速階段結(jié)束(2.8 m)時也出現(xiàn)不連續(xù)現(xiàn)象,只是由于域較大,現(xiàn)象沒有管體域的明顯.
由圖9可知,先加速后減速運動圓球在管體域中的阻力大于小計算域中的阻力.隨著運動進(jìn)行,阻力不斷增大,在加速階段結(jié)束(2.8 m)后達(dá)到最大,管體域中達(dá)5 359.5 N,小計算域中達(dá)1 976.3 N,然后出現(xiàn)突然減小而后又增大的震蕩現(xiàn)象.隨著速度減小阻力也變小,最終管體域中阻力減小到1 354.9 N,小計算域中減小到1 268.8 N.
由文獻(xiàn)[3]可知,無界黏性流體中加速運動的阻力F非取決于附加質(zhì)量、加速度和勻速運動時的阻力F定,即F非=λa+F定 (5)阻力中一部分為勻速運動阻力,而該部分阻力與速度平方成正比,因此速度越小,勻速運動阻力越小.在加速階段,速度和附加質(zhì)量都增加,因此阻力增加;在減速階段,速度減小,附加質(zhì)量增加,但增加得不太多,最終使得阻力變??;加速階段轉(zhuǎn)換為減速階段時出現(xiàn)震蕩的原因是由于圓球運動加速度發(fā)生突變,導(dǎo)致流體運動也發(fā)生劇變,因而阻力發(fā)生突變.管體域中圓球在出管體處(7 m)阻力減小程度加劇,這是由于圓球進(jìn)入喇叭口后,過流斷面開始變大,因此除速度變小的因素外,又增加一個降低阻力的因素.
在2.8 m之前,管體域中圓球一直加速運動與先加速后減速運動的阻力相同;在2.8 m之后,一直加速運動的速度繼續(xù)增大,因而阻力繼續(xù)增加,直到圓球在出管體處(7 m)達(dá)最大值11 238.7 N,之后阻力開始變小,最終減小到6 832.4 N.
由圖11和12可知,在4 m處先加速后減速圓球前后壓差小于一直加速的壓差,因此阻力較小.
由圖10可知,管體域中圓球阻力因數(shù)均隨著運動的進(jìn)行不斷減小,都大于小計算域中的阻力因數(shù).在運動初期速度較小,阻力因數(shù)與速度平方成反比,因此阻力因數(shù)較大.管體域中圓球先加速后減速阻力因數(shù)最終降為0.24,一直加速時最終降為0.28,小計算域中圓球先加速后減速時其最終降為0.22.管體域中圓球先加速后減速與一直加速運動的阻力因數(shù)在2.8 m之后也比較接近,原因是雖然一直加速的阻力大于先加速后減速,但速度也較大.
4結(jié)論
管體域中圓球的附加質(zhì)量均大于無限域中的附加質(zhì)量,并隨著運動的進(jìn)行差別增大,在出管體后達(dá)到最大值,然后開始變小.管體域中圓球先加速后減速與一直加速的附加質(zhì)量在加速階段完全相同,減速階段先加速后減速的附加質(zhì)量大于一直加速的附加質(zhì)量.
管體域中先加速后減速的圓球阻力大于無限域中的阻力.隨著運動的進(jìn)行阻力均不斷增大,并在加速階段結(jié)束后達(dá)到最大.隨著速度的減小阻力也變小.管體域中圓球在出管體處阻力減小的程度加劇.在加速階段,管體域中圓球一直加速運動與先加速后減速運動的阻力相同;在減速階段,一直加速運動的阻力繼續(xù)增加到圓球出管體處達(dá)最大,之后阻力開始變小.
管體域中圓球的阻力因數(shù)均隨著運動的進(jìn)行不斷減小,在運動初期阻力因數(shù)較大.管體域中圓球在出管體處阻力因數(shù)減小的程度加劇.管體域中的阻力因數(shù)比較接近,都大于無限域中的阻力因數(shù).
以上結(jié)論說明運動區(qū)域、加速度和流體黏性對附加質(zhì)量、阻力和阻力因數(shù)均有較大的影響.參考文獻(xiàn):
[1]王獻(xiàn)孚, 熊鰲魁. 高等流體力學(xué)[M]. 武漢: 華中科技大學(xué)出版社, 2003: 200.
[2]夏國澤. 船舶流體力學(xué)[M]. 武漢: 華中科技大學(xué)出版社, 2003: 117118.
[3]潘文全. 流體力學(xué)基礎(chǔ)[M]. 北京: 機械工業(yè)出版社, 1982: 288293.
[4]李寶元, 姜昱汐. 船舶在進(jìn)出船廂運動中的附加質(zhì)量的計算[J]. 計算力學(xué)學(xué)報, 2003, 20(3): 350354.
LI Baoyuan, JIANG Yuxi. Calculation of added mass for ship entering and leaving a shipbox[J]. Chin J Comput Mech, 2003, 20(3): 350354.
[5]馬燁, 單雪雄. 數(shù)值計算復(fù)雜外形物體附加質(zhì)量的新方法[J].