變加速運動圓球附加質(zhì)量和阻力仿真分析

姚保太+康寧+鄭偉奇

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

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

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

人們已經(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控制方程組

水中運動物體滿足非定常不可壓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計算域、邊界條件和網(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: 117118.

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

［4］李寶元, 姜昱汐. 船舶在進(jìn)出船廂運動中的附加質(zhì)量的計算［J］. 計算力學(xué)學(xué)報, 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ù)值計算復(fù)雜外形物體附加質(zhì)量的新方法［J］.