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

姚保太+康寧+鄭偉奇

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

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

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

人們已經熟知無界無黏流體中物體加速運動時的附加質量只與物體形狀及流體密度有關，與運動加速度無關，目前還缺少系統地研究在有界區域中流體黏性和加速度對做變加速運動物體附加質量的影響.由于工程中實際問題的復雜性，比如水下導彈發射問題等，因此有必要研究在有界區域黏性流體中物體做變加速運動的附加質量問題.

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

1計算模型和計算方法

1.1控制方程組

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

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

1.2計算域、邊界條件和網格劃分

1.2.1管體域

管體域由管體、喇叭口和水域組成，圓球直徑為0.254 m，圓球在管體域中運動前的計算域見圖1.通海水口邊界條件為壓強進口，出口為壓強出口，壓強均為57 m水深處壓強574 826.8 Pa；根據圓球平移運動特點，動網格模型選擇動態網格層變方法，圓球表面為無滑移運動壁面條件；管體和喇叭口為無滑移靜止壁面條件，水域邊界為滑移靜止壁面條件.圓球在管體域中運動前、后網格分布見圖2和3.采用六面體結構網格，網格數約為58萬個.

圖 1管體域中圓球運動前的計算域

Fig.1Computation domain in tube domain before ball

movement

圖 2管體域中圓球運動前的網格分布

Fig.2Mesh distribution of tube domain before ball movement

圖 3管體域中圓球運動后的網格分布

Fig.3Mesh distribution of tube domain after ball movement

1.2.2小計算域

作為對比，對圓球在小計算域中的運動情況進行計算.小計算域為圓柱形，取較大尺度，直徑為3.6 m，長度為14 m.縱向對稱面網格分布見圖4.六面體結構網格，網格數約為139萬個.計算域左側為壓強入口邊界條件，右側為壓強出口條件，相對壓強為0；其他邊界和動網格模型同管體域.

圖 4小計算域縱向對稱面的網格分布

Fig.4Mesh distribution in longitudinal symmetry plane of

small computation domain

1.2.3大計算域

為驗證附加質量仿真方法的可靠性，對圓球在大計算域中的運動情況進行計算.大計算域仍為圓柱形，直徑加大到14.4 m，長度加大到20.86 m.采用六面體結構網格，網格數約為264萬個.邊界條件和動網格模型同管體域.1.3圓球變速運動參數

1.3.1先加速后減速

為模擬發射體在管體內直至出管的運動情況，對圓球先加速后減速運動情況進行計算.加速階段模擬發射時在主動力作用下的運動情況，減速階段模擬在慣性力作用下直至出管的運動情況，具體參數見表1.

表 1圓球先加速后減速的運動參數

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一直加速

作為對比，對圓球在管體域內一直加速運動工況進行計算，加速度和初速度與先加速后減速的加速階段相同，位移與先加速后減速的總位移相近，見表2.

表 2圓球一直加速的運動參數

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求解方法和計算參數

速度和壓強耦合方法采用SIMPLE方法，壓力修正方程離散格式采用Standard格式，對流項離散格式采用2階迎風格式，湍流模型為RNG kε湍流模型，近壁區域處理采用非平衡壁面函數.

海水密度為1 028 kg/m3，動力黏性系數為1.394×10-3 kg/(m?s)，參考壓強為1標準大氣壓101 325 Pa.計算初始速度為0，初始壓強為57 m水深處壓強574 826.8 Pa.

2仿真方法驗證

2.1阻力

分2種情況對圓球阻力因數進行驗證，分別為水繞靜止圓球運動和圓球在靜水中運動.

圓球實驗阻力因數的范圍［3］為0～106，根據定義可得圓球速度范圍為0～5.339 m/s，因此對最大阻力因數106和較小阻力因數104，其分別對應最大速度5.339 m/s和較小速度0.053 4 m/s的情況進行計算.采用小計算域進行計算.

2.1.1靜止圓球

計算域左側為速度入口邊界條件，速度大小分別為5.339 0和0.053 4 m/s；右側為壓強出口邊界條件，相對壓強為0.

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

表 3圓球靜止時的阻力因數

Tab.3Drag factor of stationary ball速度/(m/s)實驗值計算值誤差/%0.053 40.400.4164.005.339 00.170.1652.94

由表3可知,在2種速度情況下計算的圓球阻力因數與實驗結果的誤差都在4%以內,說明圓球靜止時阻力計算結果可靠.

2.1.2運動圓球

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

表 4圓球運動時的阻力因數

Tab.4Drag facctor of moving ball速度/(m/s)實驗值計算值誤差/%0.053 40.400.3717.255.339 00.170.18710.00

由表4可知,在2種速度情況下計算的圓球阻力因數與實驗結果的誤差都在10%以內,說明圓球運動時阻力計算結果可靠.

2.2附加質量

2.2.1附加質量計算公式

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

2.2.2圓球附加質量的理論解

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

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

2.2.3圓球附加質量的計算結果

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

圖 5小計算域中加速度對圓球附加質量的影響

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

small computation domain

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

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

附加質量的影響

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

ball at a=5 m/s2 in inviscid fluid

2.2.4橢球體附加質量的實驗驗證

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

圖 7橢球體縱向對稱面的網格分布圖

Fig.7Mesh distribution in longitudinal symmetry

plane of ellipsoid

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

表 5橢球體的附加質量

Tab.5Added mass of ellipsoid加速時間/s實驗值/kg計算值/kg誤差/%0.11.701.742.350.51.701.774.120.91.701.858.82

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

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

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

圖 8管體域和小計算域中圓球的附加質量

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管體域和小計算域中圓球的阻力因數

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

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

(b)壓強云圖,Pa

圖 11在管體域中4 m處先加速后減速圓球的計算結果

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處一直加速圓球的計算結果

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

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

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

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

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

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

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

4結論

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

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

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

以上結論說明運動區域、加速度和流體黏性對附加質量、阻力和阻力因數均有較大的影響.參考文獻：

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

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

［3］潘文全. 流體力學基礎［M］. 北京: 機械工業出版社, 1982: 288293.

［4］李寶元, 姜昱汐. 船舶在進出船廂運動中的附加質量的計算［J］. 計算力學學報, 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］馬燁, 單雪雄. 數值計算復雜外形物體附加質量的新方法［J］.