姚保太+康寧+鄭偉奇

作者簡介： 姚保太（1973—），男，河南獲嘉人，高級工程師，研究方向?yàn)闄C(jī)械制造及仿真計(jì)算，(Email)yaobaotai@163.com

通信作者： 康寧（1963—），女，遼寧金州人，副教授，博士，研究領(lǐng)域方向?yàn)榱黧w仿真計(jì)算，(Email)kangning@buaa.edu.cn0引言

附加質(zhì)量和阻力計(jì)算是水下物體運(yùn)動分析的關(guān)鍵問題.［1］以往對附加質(zhì)量的計(jì)算主要是根據(jù)勢流理論［14］，如在無界無黏流體中球體等簡單形狀物體的附加質(zhì)量存在解析解,復(fù)雜形狀物體的附加質(zhì)量通過估算法或三維面元法求得.近年來，隨著計(jì)算流體力學(xué)的發(fā)展，開始用數(shù)值模擬的方法計(jì)算物體的附加質(zhì)量，發(fā)展基于動網(wǎng)格技術(shù)計(jì)算附加質(zhì)量的新方法.［5］數(shù)值方法有更廣泛的應(yīng)用空間，不僅可以計(jì)算理想流體中的附加質(zhì)量［56］，還可以計(jì)算黏性流體中的附加質(zhì)量［710］.

人們已經(jīng)熟知無界無黏流體中物體加速運(yùn)動時(shí)的附加質(zhì)量只與物體形狀及流體密度有關(guān)，與運(yùn)動加速度無關(guān)，目前還缺少系統(tǒng)地研究在有界區(qū)域中流體黏性和加速度對做變加速運(yùn)動物體附加質(zhì)量的影響.由于工程中實(shí)際問題的復(fù)雜性，比如水下導(dǎo)彈發(fā)射問題等，因此有必要研究在有界區(qū)域黏性流體中物體做變加速運(yùn)動的附加質(zhì)量問題.

本文采用數(shù)值模擬方法和動網(wǎng)格技術(shù)，計(jì)算管體域中黏性不可壓流體中圓球做變加速運(yùn)動時(shí)的附加質(zhì)量、阻力和阻力因數(shù)，深入分析和研究在管體域中圓球變加速運(yùn)動對附加質(zhì)量、阻力和阻力因數(shù)的影響.

1計(jì)算模型和計(jì)算方法

1.1控制方程組

水中運(yùn)動物體滿足非定常不可壓NS方程，

連續(xù)方程： vixi=0(1)動量方程：vit+vjvixj=Fi－1ρpxi+μρxjvixj(2)式中：vi和vj為速度；F為質(zhì)量力；p為壓力；ρ為流體密度；μ為流體動力黏性系數(shù).

1.2計(jì)算域、邊界條件和網(wǎng)格劃分

1.2.1管體域

管體域由管體、喇叭口和水域組成，圓球直徑為0.254 m，圓球在管體域中運(yùn)動前的計(jì)算域見圖1.通海水口邊界條件為壓強(qiáng)進(jìn)口，出口為壓強(qiáng)出口，壓強(qiáng)均為57 m水深處壓強(qiáng)574 826.8 Pa；根據(jù)圓球平移運(yùn)動特點(diǎn)，動網(wǎng)格模型選擇動態(tài)網(wǎng)格層變方法，圓球表面為無滑移運(yùn)動壁面條件；管體和喇叭口為無滑移靜止壁面條件，水域邊界為滑移靜止壁面條件.圓球在管體域中運(yùn)動前、后網(wǎng)格分布見圖2和3.采用六面體結(jié)構(gòu)網(wǎng)格，網(wǎng)格數(shù)約為58萬個(gè).

圖 1管體域中圓球運(yùn)動前的計(jì)算域

Fig.1Computation domain in tube domain before ball

movement

圖 2管體域中圓球運(yùn)動前的網(wǎng)格分布

Fig.2Mesh distribution of tube domain before ball movement

圖 3管體域中圓球運(yùn)動后的網(wǎng)格分布

Fig.3Mesh distribution of tube domain after ball movement

1.2.2小計(jì)算域

作為對比，對圓球在小計(jì)算域中的運(yùn)動情況進(jìn)行計(jì)算.小計(jì)算域?yàn)閳A柱形，取較大尺度，直徑為3.6 m，長度為14 m.縱向?qū)ΨQ面網(wǎng)格分布見圖4.六面體結(jié)構(gòu)網(wǎng)格，網(wǎng)格數(shù)約為139萬個(gè).計(jì)算域左側(cè)為壓強(qiáng)入口邊界條件，右側(cè)為壓強(qiáng)出口條件，相對壓強(qiáng)為0；其他邊界和動網(wǎng)格模型同管體域.

圖 4小計(jì)算域縱向?qū)ΨQ面的網(wǎng)格分布

Fig.4Mesh distribution in longitudinal symmetry plane of

small computation domain

1.2.3大計(jì)算域

為驗(yàn)證附加質(zhì)量仿真方法的可靠性，對圓球在大計(jì)算域中的運(yùn)動情況進(jìn)行計(jì)算.大計(jì)算域仍為圓柱形，直徑加大到14.4 m，長度加大到20.86 m.采用六面體結(jié)構(gòu)網(wǎng)格，網(wǎng)格數(shù)約為264萬個(gè).邊界條件和動網(wǎng)格模型同管體域.1.3圓球變速運(yùn)動參數(shù)

1.3.1先加速后減速

為模擬發(fā)射體在管體內(nèi)直至出管的運(yùn)動情況，對圓球先加速后減速運(yùn)動情況進(jìn)行計(jì)算.加速階段模擬發(fā)射時(shí)在主動力作用下的運(yùn)動情況，減速階段模擬在慣性力作用下直至出管的運(yùn)動情況，具體參數(shù)見表1.

表 1圓球先加速后減速的運(yùn)動參數(shù)

Tab.1Movement parameters when ball is accelerated motion first and then decelerated初速度/

(m/s)加速度/

(m/s2)時(shí)間/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)一直加速運(yùn)動工況進(jìn)行計(jì)算，加速度和初速度與先加速后減速的加速階段相同，位移與先加速后減速的總位移相近，見表2.

表 2圓球一直加速的運(yùn)動參數(shù)

Tab.2Movement parameters when ball is accelerated motion all the time初速度/

(m/s)加速度/

(m/s2)時(shí)間/s末速度/

(m/s)位移/m052.0000.58930.6209.020

1.4求解方法和計(jì)算參數(shù)

速度和壓強(qiáng)耦合方法采用SIMPLE方法，壓力修正方程離散格式采用Standard格式，對流項(xiàng)離散格式采用2階迎風(fēng)格式，湍流模型為RNG kε湍流模型，近壁區(qū)域處理采用非平衡壁面函數(shù).

海水密度為1 028 kg/m3，動力黏性系數(shù)為1.394×10-3 kg/(m?s)，參考壓強(qiáng)為1標(biāo)準(zhǔn)大氣壓101 325 Pa.計(jì)算初始速度為0，初始壓強(qiáng)為57 m水深處壓強(qiáng)574 826.8 Pa.

2仿真方法驗(yàn)證

2.1阻力

分2種情況對圓球阻力因數(shù)進(jìn)行驗(yàn)證，分別為水繞靜止圓球運(yùn)動和圓球在靜水中運(yùn)動.

圓球?qū)嶒?yàn)阻力因數(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ì)算.采用小計(jì)算域進(jìn)行計(jì)算.

2.1.1靜止圓球

計(jì)算域左側(cè)為速度入口邊界條件，速度大小分別為5.339 0和0.053 4 m/s；右側(cè)為壓強(qiáng)出口邊界條件，相對壓強(qiáng)為0.

圓球表面為無滑移靜止壁面邊界條件；除進(jìn)口、出口和圓球表面以外的邊界條件為滑移靜止壁面邊界條件.計(jì)算結(jié)果見表3.

表 3圓球靜止時(shí)的阻力因數(shù)

Tab.3Drag factor of stationary ball速度/(m/s)實(shí)驗(yàn)值計(jì)算值誤差/%0.053 40.400.4164.005.339 00.170.1652.94

由表3可知,在2種速度情況下計(jì)算的圓球阻力因數(shù)與實(shí)驗(yàn)結(jié)果的誤差都在4%以內(nèi),說明圓球靜止時(shí)阻力計(jì)算結(jié)果可靠.

2.1.2運(yùn)動圓球

計(jì)算域左側(cè)為壓強(qiáng)入口條件，右側(cè)為壓強(qiáng)出口條件，相對壓強(qiáng)為0.動網(wǎng)格模型及圓球表面邊界條件均同管體域中的運(yùn)動圓球，速度大小分別為5.339 0和0.053 4 m/s；其余邊界為滑移靜止壁面條件；初始速度和壓強(qiáng)均為0.計(jì)算結(jié)果見表4.

表 4圓球運(yùn)動時(shí)的阻力因數(shù)

Tab.4Drag facctor of moving ball速度/(m/s)實(shí)驗(yàn)值計(jì)算值誤差/%0.053 40.400.3717.255.339 00.170.18710.00

由表4可知,在2種速度情況下計(jì)算的圓球阻力因數(shù)與實(shí)驗(yàn)結(jié)果的誤差都在10%以內(nèi),說明圓球運(yùn)動時(shí)阻力計(jì)算結(jié)果可靠.

2.2附加質(zhì)量

2.2.1附加質(zhì)量計(jì)算公式

根據(jù)文獻(xiàn)［2］，附加質(zhì)量λ=ρv2iViv2o(3)式中：vi為第i控制體的流體速度；Vi為第i控制體的體積；vo為物體運(yùn)動速度.

2.2.2圓球附加質(zhì)量的理論解

根據(jù)文獻(xiàn)［1］，在無界無黏流體情況下，圓球直線加速運(yùn)動附加質(zhì)量的理論解λ=23πρR3 (4)式中：R為圓球半徑.

圓球直徑、水密度和動力黏性因數(shù)同前，這時(shí)圓球附加質(zhì)量的理論解為4.41 kg.

2.2.3圓球附加質(zhì)量的計(jì)算結(jié)果

為驗(yàn)證附加質(zhì)量仿真方法的可靠性，對小計(jì)算域中圓球在無黏流體中做不同加速度的運(yùn)動進(jìn)行計(jì)算，圓球的初速度均為0，加速時(shí)間為1.1 s，加速度分別為0.01，0.10，1.00和5.00 m/s2，并與圓球附加質(zhì)量的理論解進(jìn)行對比.所采用的網(wǎng)格、邊界條件、初始條件、計(jì)算設(shè)置和動網(wǎng)格模型均與運(yùn)動圓球阻力驗(yàn)證時(shí)一樣.計(jì)算的附加質(zhì)量見圖5.

圖 5小計(jì)算域中加速度對圓球附加質(zhì)量的影響

Fig.5Effect of acceleration on added mass of ball in

small computation domain

由圖5可知，時(shí)間和加速度越大附加質(zhì)量越大.在加速度較小時(shí)（0.01，0.10和1.00 m/s2）計(jì)算的附加質(zhì)量與理論解較為接近；在加速度較大時(shí)（5.00 m/s2），運(yùn)動初期計(jì)算的附加質(zhì)量與理論解較接近，隨著時(shí)間增大其值開始大于理論解.這是由于隨著加速度和時(shí)間的增大，運(yùn)動速度增大，導(dǎo)致計(jì)算域邊界效應(yīng)的影響也隨之增大，附加質(zhì)量增加，因而與理論解的差別增大.

為進(jìn)一步研究計(jì)算域邊界效應(yīng)的影響，在大計(jì)算域中計(jì)算加速度為5.00 m/s2的工況，初速度和加速時(shí)間均同前文，無黏流體中的附加質(zhì)量的計(jì)算結(jié)果見圖6，并與之前小計(jì)算域結(jié)果進(jìn)行對比.由圖6可知，除運(yùn)動初期，無黏流體中小計(jì)算域附加質(zhì)量均大于大計(jì)算域附加質(zhì)量，大計(jì)算域附加質(zhì)量均大于理論解，并且隨著時(shí)間增加速度加大，差別也越大.大計(jì)算域附加質(zhì)量與理論解較接近，小計(jì)算域附加質(zhì)量與理論解差別較大.由圖5和6的結(jié)果可知，加速度越小或計(jì)算域越大，邊界效應(yīng)越小，則計(jì)算的附加質(zhì)量越接近理論解.圖 6在a=5 m/s2時(shí)無黏流體中計(jì)算域尺度對圓球

附加質(zhì)量的影響

Fig.6Effect of computation domain scale on added mass of

ball at a=5 m/s2 in inviscid fluid

2.2.4橢球體附加質(zhì)量的實(shí)驗(yàn)驗(yàn)證

為驗(yàn)證附加質(zhì)量的理論計(jì)算公式（3），將黏性流體中加速運(yùn)動橢球體附加質(zhì)量的仿真結(jié)果與文獻(xiàn)［6］的實(shí)驗(yàn)結(jié)果進(jìn)行對比.橢球長軸為0.4 m，短軸為0.2 m，與文獻(xiàn)［6］中的完全一樣.計(jì)算域?yàn)閳A柱形，直徑為6 m，長度為14 m.計(jì)算時(shí)水溫度為5 ℃，密度為1 000 kg/m3，動力黏性系數(shù)為1.518×10-3 kg/(m?s).橢球體為六面體網(wǎng)格，網(wǎng)格數(shù)為106萬個(gè)，網(wǎng)格分布見圖7.

圖 7橢球體縱向?qū)ΨQ面的網(wǎng)格分布圖

Fig.7Mesh distribution in longitudinal symmetry

plane of ellipsoid

邊界條件、初始條件和其他設(shè)置均同阻力驗(yàn)證時(shí)的運(yùn)動圓球，橢球動網(wǎng)格模型同運(yùn)動圓球.橢球體運(yùn)動初速度為0，加速度為1 m/s2，加速時(shí)間分別為0，0.5和0.9 s.根據(jù)式（3）計(jì)算的附加質(zhì)量和與文獻(xiàn)［6］的對比見表5.由表5可知，計(jì)算的附加質(zhì)量與實(shí)驗(yàn)結(jié)果比較接近，誤差在9%以內(nèi)，說明附加質(zhì)量的計(jì)算結(jié)果基本正確.

表 5橢球體的附加質(zhì)量

Tab.5Added mass of ellipsoid加速時(shí)間/s實(shí)驗(yàn)值/kg計(jì)算值/kg誤差/%0.11.701.742.350.51.701.774.120.91.701.858.82

綜合前文分析結(jié)果可知，附加質(zhì)量的仿真方法較為可靠.3計(jì)算結(jié)果及分析

對圓球在管體域黏性流體中先加速后減速和一直加速2種變速運(yùn)動的附加質(zhì)量、阻力和阻力因數(shù)進(jìn)行仿真計(jì)算，并對小計(jì)算域黏性流體中圓球先加速后減速運(yùn)動進(jìn)行仿真計(jì)算，同時(shí)給出圓球在無界無黏流體中附加質(zhì)量的理論解.

管體域和小計(jì)算域中圓球在黏性流體中附加質(zhì)量、阻力和阻力因數(shù)隨位移的變化情況見圖8～10，在管體域中4 m處圓球先加速后減速和一直加速的速度矢量圖和壓強(qiáng)云圖見圖11和12.

圖 8管體域和小計(jì)算域中圓球的附加質(zhì)量

Fig.8Added masses of ball in tube domain and small

computation domain

圖 9管體域和小計(jì)算域中圓球的阻力

Fig.9Drags of ball in tube domain and small computation

domain

圖 10管體域和小計(jì)算域中圓球的阻力因數(shù)

Fig.10Drag factors of ball in tube domain and small

computation domain(a)速度矢量圖,m/s

(b)壓強(qiáng)云圖,Pa

圖 11在管體域中4 m處先加速后減速圓球的計(jì)算結(jié)果

Fig.11Calculation results at 4 m in tube domain while ball is accelerated first and then decelerated

(a)速度矢量圖,m/s

(b)壓強(qiáng)云圖,Pa

圖 12在管體域中4 m處一直加速圓球的計(jì)算結(jié)果

Fig.12Calculation results at 4 m in tube domain while ball is accelerated all the time

由圖8可知，在運(yùn)動開始時(shí)，小計(jì)算域中附加質(zhì)量仿真值與理論解4.41 kg吻合較好；隨著運(yùn)動進(jìn)行仿真值開始大于理論解，最后達(dá)到20.8 kg.管體域中附加質(zhì)量均大于小計(jì)算域中附加質(zhì)量，并隨著運(yùn)動進(jìn)行差別逐漸增大，在出管體（7 m）后達(dá)到最大值，然后開始變小.先加速后減速運(yùn)動的最大值為58.6 kg，最終降為48.1 kg；一直加速運(yùn)動的最大值為44.7 kg，最終降為33.9 kg.由于管體域尺寸小于小計(jì)算域，邊界效應(yīng)較大，因而附加質(zhì)量較大，出喇叭口后過流斷面變大，邊界效應(yīng)減小，導(dǎo)致附加質(zhì)量變小.管體域中圓球先加速后減速和一直加速的附加質(zhì)量在2.8 m之前完全一樣，因?yàn)樵撾A段2種工況運(yùn)動狀態(tài)完全一樣；在2.8 m之后，一直加速運(yùn)動附加質(zhì)量仍保持連續(xù)增加態(tài)勢，而先加速后減速運(yùn)動附加質(zhì)量在2.8 m處出現(xiàn)不連續(xù)現(xiàn)象，開始大于一直加速附加質(zhì)量.先加速后減速運(yùn)動由加速狀態(tài)變?yōu)闇p速狀態(tài)時(shí)，運(yùn)動狀態(tài)的改變遠(yuǎn)大于一直加速運(yùn)動狀態(tài).由圖11和12可看出，在4 m處，先加速后減速圓球周圍速度雖然小于一直加速的速度，但尾跡流影響范圍明顯較大，在圓球尾部較遠(yuǎn)處還卷起較大尾渦，因而導(dǎo)致附加質(zhì)量非但沒有變小，反而增大.小計(jì)算域附加質(zhì)量在加速階段結(jié)束（2.8 m）時(shí)也出現(xiàn)不連續(xù)現(xiàn)象，只是由于域較大，現(xiàn)象沒有管體域的明顯.

由圖9可知，先加速后減速運(yùn)動圓球在管體域中的阻力大于小計(jì)算域中的阻力.隨著運(yùn)動進(jìn)行，阻力不斷增大，在加速階段結(jié)束（2.8 m）后達(dá)到最大，管體域中達(dá)5 359.5 N，小計(jì)算域中達(dá)1 976.3 N，然后出現(xiàn)突然減小而后又增大的震蕩現(xiàn)象.隨著速度減小阻力也變小，最終管體域中阻力減小到1 354.9 N，小計(jì)算域中減小到1 268.8 N.

由文獻(xiàn)［3］可知，無界黏性流體中加速運(yùn)動的阻力F非取決于附加質(zhì)量、加速度和勻速運(yùn)動時(shí)的阻力F定，即F非=λa+F定 (5)阻力中一部分為勻速運(yùn)動阻力，而該部分阻力與速度平方成正比，因此速度越小，勻速運(yùn)動阻力越小.在加速階段，速度和附加質(zhì)量都增加，因此阻力增加；在減速階段，速度減小，附加質(zhì)量增加，但增加得不太多，最終使得阻力變??；加速階段轉(zhuǎn)換為減速階段時(shí)出現(xiàn)震蕩的原因是由于圓球運(yùn)動加速度發(fā)生突變，導(dǎo)致流體運(yùn)動也發(fā)生劇變，因而阻力發(fā)生突變.管體域中圓球在出管體處（7 m）阻力減小程度加劇，這是由于圓球進(jìn)入喇叭口后，過流斷面開始變大，因此除速度變小的因素外，又增加一個(gè)降低阻力的因素.

在2.8 m之前,管體域中圓球一直加速運(yùn)動與先加速后減速運(yùn)動的阻力相同；在2.8 m之后,一直加速運(yùn)動的速度繼續(xù)增大，因而阻力繼續(xù)增加,直到圓球在出管體處（7 m）達(dá)最大值11 238.7 N，之后阻力開始變小，最終減小到6 832.4 N.

由圖11和12可知，在4 m處先加速后減速圓球前后壓差小于一直加速的壓差，因此阻力較小.

由圖10可知，管體域中圓球阻力因數(shù)均隨著運(yùn)動的進(jìn)行不斷減小，都大于小計(jì)算域中的阻力因數(shù).在運(yùn)動初期速度較小，阻力因數(shù)與速度平方成反比，因此阻力因數(shù)較大.管體域中圓球先加速后減速阻力因數(shù)最終降為0.24，一直加速時(shí)最終降為0.28，小計(jì)算域中圓球先加速后減速時(shí)其最終降為0.22.管體域中圓球先加速后減速與一直加速運(yùn)動的阻力因數(shù)在2.8 m之后也比較接近，原因是雖然一直加速的阻力大于先加速后減速,但速度也較大.

4結(jié)論

管體域中圓球的附加質(zhì)量均大于無限域中的附加質(zhì)量，并隨著運(yùn)動的進(jìn)行差別增大，在出管體后達(dá)到最大值，然后開始變小.管體域中圓球先加速后減速與一直加速的附加質(zhì)量在加速階段完全相同，減速階段先加速后減速的附加質(zhì)量大于一直加速的附加質(zhì)量.

管體域中先加速后減速的圓球阻力大于無限域中的阻力.隨著運(yùn)動的進(jìn)行阻力均不斷增大，并在加速階段結(jié)束后達(dá)到最大.隨著速度的減小阻力也變小.管體域中圓球在出管體處阻力減小的程度加劇.在加速階段,管體域中圓球一直加速運(yùn)動與先加速后減速運(yùn)動的阻力相同；在減速階段,一直加速運(yùn)動的阻力繼續(xù)增加到圓球出管體處達(dá)最大，之后阻力開始變小.

管體域中圓球的阻力因數(shù)均隨著運(yùn)動的進(jìn)行不斷減小，在運(yùn)動初期阻力因數(shù)較大.管體域中圓球在出管體處阻力因數(shù)減小的程度加劇.管體域中的阻力因數(shù)比較接近，都大于無限域中的阻力因數(shù).

以上結(jié)論說明運(yùn)動區(qū)域、加速度和流體黏性對附加質(zhì)量、阻力和阻力因數(shù)均有較大的影響.參考文獻(xiàn)：

［1］王獻(xiàn)孚, 熊鰲魁. 高等流體力學(xué)［M］. 武漢: 華中科技大學(xué)出版社, 2003: 200.

［2］夏國澤. 船舶流體力學(xué)［M］. 武漢: 華中科技大學(xué)出版社, 2003: 117118.

［3］潘文全. 流體力學(xué)基礎(chǔ)［M］. 北京: 機(jī)械工業(yè)出版社, 1982: 288293.

［4］李寶元, 姜昱汐. 船舶在進(jìn)出船廂運(yùn)動中的附加質(zhì)量的計(jì)算［J］. 計(jì)算力學(xué)學(xué)報(bào), 2003, 20(3): 350354.

LI Baoyuan, JIANG Yuxi. Calculation of added mass for ship entering and leaving a shipbox［J］. Chin J Comput Mech, 2003, 20(3): 350354.

［5］馬燁, 單雪雄. 數(shù)值計(jì)算復(fù)雜外形物體附加質(zhì)量的新方法［J］.