基于DES模型圓柱彈丸水中繞流特性仿真

張弛宇，　郭　銳，　劉　敏，　劉榮忠，　周　昊，　陳　亮

（1.南京理工大學(xué) 機(jī)械工程學(xué)院， 江蘇 南京， 210094； 2. 陸軍西安軍事代表局， 陜西 西安， 710032）

基于DES模型圓柱彈丸水中繞流特性仿真

張弛宇1，郭銳1，劉敏2，劉榮忠1，周昊1，陳亮1

（1.南京理工大學(xué) 機(jī)械工程學(xué)院， 江蘇 南京， 210094； 2. 陸軍西安軍事代表局， 陜西 西安， 710032）

圓柱型裝藥是水下戰(zhàn)斗部的主要形式， 文中基于 Fluent軟件， 采用分離渦模擬（DES）方法對處于高雷諾數(shù)下圓柱彈丸進(jìn)行了二維水中繞流仿真， 重點(diǎn)研究了 25 m/s來流速度下圓形彈丸水中繞流特性， 得到了升力系數(shù)、阻力系數(shù)以及斯特勞哈爾數(shù)。通過仿真得到了周期性交替脫落的卡門渦街結(jié)構(gòu)， 分析了流場在不同截面時均速度規(guī)律， 以及圓柱周向時均穩(wěn)態(tài)壓力系數(shù)， 其圖像規(guī)律也與已知大湍渦模擬（LES）方法高度吻合。通過對不同來速下彈丸阻力系數(shù)和斯特勞哈爾數(shù)的分析， 結(jié)果表明， 隨著來流速度增大， 阻力系數(shù)減少， 圓柱彈丸在水中運(yùn)動所受的阻力減小， 而斯特勞哈爾數(shù)增大意味圓柱彈丸兩側(cè)渦脫落的頻率增大， 當(dāng)渦脫落頻率與物體固有頻率接近時， 會引起共振， 造成彈的損傷。文中的研究可為水下戰(zhàn)斗部彈丸設(shè)計提供參考。

水下戰(zhàn)斗部； 圓形彈丸繞流； 分離渦模擬（DES）湍流模型； 高雷諾數(shù)； 卡門渦街

0　引言

柱體繞流是經(jīng)典的流體力學(xué)課題， 于 20世紀(jì)60～90年代被廣泛關(guān)注。水下彈道也存在彈丸柱體繞流現(xiàn)象， 在特定水下彈道工況中柱體繞流會產(chǎn)生卡門渦街， 此時柱體上下邊界層交替脫落形成渦旋， 對彈丸作用周期性變化的脈動升阻力，導(dǎo)致柱體彈丸振動， 當(dāng)脈動升阻力頻率與柱體彈丸固有頻率接近時， 會誘發(fā)共振， 對彈丸造成結(jié)構(gòu)破壞， 并可能增大彈丸阻力， 產(chǎn)生噪音， 使彈丸偏離設(shè)計彈道。由于水下彈道通常處于高雷諾數(shù)的臨界區(qū)和過臨界區(qū)， 理論和實(shí)驗(yàn)研究存在諸多條件限制， 因而數(shù)值仿真已成為研究柱體繞流的有效手段。

目前， 國內(nèi)外學(xué)者針對柱體繞流的研究已有很多。史里希廷［1］給出了來流速度小于100m/s即雷諾數(shù)小于1×107的實(shí)驗(yàn)結(jié)果； E.Achenbach［2］對雷諾數(shù)處于 6×104＜Re＜2×106的光滑表面單圓柱繞流進(jìn)行了實(shí)驗(yàn)研究， 得到了圓柱周圍局部壓力和摩擦力的分布； Yokuda和 Ramaprian［3］對亞臨界區(qū)圓柱繞流進(jìn)行了實(shí)驗(yàn)研究， 但目前已知關(guān)于圓柱繞流的研究多采用大湍渦模擬 （large eddy simulation， LES）方法， 雷諾數(shù)一般在1×106以下， 且采用LES需要將計算網(wǎng)格劃分到慣性子尺度以內(nèi)， 計算時間和網(wǎng)格數(shù)量隨著雷諾數(shù)的增大呈幾何倍數(shù)關(guān)系［4］； 雷諾平均（Reynolds average navier Stockes， RANS）雖然計算量較小， 但文獻(xiàn)［5］的研究結(jié)果表明， 二維RANS模型難以合理地預(yù)測非定常柱體繞流； Spalart［6］在1997年提出了分離渦模擬（detached eddy simulation， DES）的方法，將 RANS與 LES相結(jié)合， 既能在附面層內(nèi)發(fā)揮RANS計算量小的優(yōu)點(diǎn)， 又可以在遠(yuǎn)離物面的區(qū)域?qū)Υ蟪叨确蛛x湍流流動進(jìn)行較好的模擬。

基于此， 文中采用S-A模型下的DES方法進(jìn)行數(shù)值仿真， 驗(yàn)證了DES模型在大雷諾數(shù)研究中的準(zhǔn)確性， 得出了圓柱彈丸的平均阻力系數(shù)和斯特勞哈爾數(shù)， 繼而研究卡門渦街現(xiàn)象以及流場中幾個截面時均速度分布以及圓柱周向穩(wěn)態(tài)壓力系數(shù)。最后分析了平均阻力系數(shù)和斯特勞哈爾數(shù)與來流速度的關(guān)系， 以期為后續(xù)工程研究提供依據(jù)。

1　仿真方法及驗(yàn)證

1.1模型建立與仿真實(shí)驗(yàn)方案

假設(shè): 實(shí)際水下彈道中彈丸是運(yùn)動的， 水流體介質(zhì)是靜止的， 文中假設(shè)柱體彈丸靜止不動，賦予流體水介質(zhì)相對速度即來流速度 u， 來研究彈丸柱體繞流， 并認(rèn)為二者是等效的。

圖1為計算流場區(qū)域模型圖， 參考文獻(xiàn)［4～10］并經(jīng)過多次仿真實(shí)驗(yàn)確定計算流場尺寸和柱體的相對位置。取坐標(biāo)原點(diǎn)于彈丸圓心（x=0， y=0）， 彈丸直徑為D=100 mm， 來流速度u的方向?yàn)閤軸正方向； 矩形計算流場為（40×16）D， 上游來流區(qū)域8D， 為觀察渦街脫落后的形態(tài)取下游尾流區(qū)域?yàn)?32D； 取坐標(biāo)原點(diǎn)到上下邊界的距離各為 8D，此時， 上下壁面邊界對計算結(jié)果的影響可以忽略不計。

圖1　計算區(qū)域示意圖Fig. 1　Schematic of calculation area

1.2網(wǎng)格劃分及湍流模型

柱體繞流是由柱體外圍邊界層復(fù)雜湍流流動引起的， 所以網(wǎng)格質(zhì)量直接影響數(shù)值計算結(jié)果。

由圖2所示， 將圓形計算流場劃分為9個區(qū)域， 在不影響數(shù)值計算結(jié)果的前提下， 為減少網(wǎng)格數(shù)目， 加快計算速度， 網(wǎng)格密度從中心的 5區(qū)向四周逐漸減小。為了準(zhǔn)確地觀察邊界層的復(fù)雜流動以及渦街脫落形態(tài)， 采用O型拓?fù)浣Y(jié)構(gòu)對柱體周圍加密， 比例因子為0.26， 如圖2所示。

圖2　流場網(wǎng)格及柱體外圍網(wǎng)格Fig. 2　Grids of flow filed and cylinder periphery

［5］中網(wǎng)格無關(guān)性檢驗(yàn)計算結(jié)果， 壁面第1層網(wǎng)格的厚度為0.000 01 m， 采用小于1.07的網(wǎng)格長度增長率， 保證柱體邊界層位置獲得高分辨率網(wǎng)格， 最終網(wǎng)格總數(shù)取5.6萬。

根據(jù)水介質(zhì)的壓縮性， 入口邊界采用速度入口邊界條件， 設(shè)置速度25 m/s； 出口采用outflow出流條件； 為避免壁面對流場速度的影響， 上下壁面采用滑移壁面， 滑移速度與入口速度保持一致；柱體表面為無滑移壁面， 采用SIMPLEC算法求解壓力速度耦合， 動量采用2階迎風(fēng)格式， 壓力方程為2階精度離散， 瞬態(tài)格式采用2階隱式。流體介質(zhì)取為水介質(zhì)， 密度為ρ=998.2 kg/m3， 動力粘度為ν=1.003×10-3kg/（m·s）， 雷諾數(shù)計算公式為Re=ρuD/ν。

2　仿真結(jié)果與分析

根據(jù)彈丸的水下彈道工況設(shè)計仿真實(shí)驗(yàn)方案，以圓柱彈丸來流速度為25 m/s來研究彈丸水下流場特性、升力系數(shù)、阻力系數(shù)及斯特勞哈爾數(shù)的變化規(guī)律， 雷諾數(shù)為2.5×106。

2.1升力系數(shù)、阻力系數(shù)及斯特勞哈爾數(shù)

升力系數(shù)Cl和阻力系數(shù) Cd是描述繞流對柱體作用力的重要參數(shù)， 斯特勞哈爾數(shù)Str則是描述旋渦脫落非定常性的特征參數(shù)， 反映了繞流對柱體作用的非定常特征。三者定義分別為

式中: F1為圓柱受到的橫向力； Fd為圓柱受到的流動方向的力；f為旋渦脫落頻率； u0為均勻來流速度。

圖3是Re為2.5×106時升力系數(shù)和阻力系數(shù)在某段時間范圍內(nèi)的變化規(guī)律。由圖3可以看出，計算穩(wěn)定后的升力系數(shù)和阻力系數(shù)隨時間變化的曲線都近似于正弦曲線。

總體來說， 在圖 2所描述的時間范圍內(nèi)， Cd變化了15個周期， Cl變化了7.5個周期。也就是說升力變化頻率是阻力變化頻率的2倍， 與之前關(guān)于圓柱繞流的研究基本一致。這是由升力系數(shù)、阻力的方向和旋渦脫落特征決定的， 上、下渦的脫落各引起阻力改變一次， 而這2個渦共同影響升力變化一次。

圖3　升力系數(shù)阻力系數(shù)隨時間變化曲線Fig. 3　Curves of lift coefficient and drag coefficient versus time

圖4是升力自功率頻譜圖?？v坐標(biāo)A表示幅值。在文中雷諾數(shù)下， 自升力功率頻譜有一個明顯的峰值， 根據(jù)該值可計算對應(yīng)的Str。

來流速度為25 m/s時， 在該DES模擬方法下，得到Cl， Cd以及Str， 見表1。

圖4　升力自功率頻譜Fig. 4　Self-power spectrum of lift

文中和文獻(xiàn)［6］得到的平均阻力系數(shù)都大于H.史里希延［1］所提到的實(shí)驗(yàn)值， 文中的 Str數(shù)和文獻(xiàn)［1］數(shù)據(jù)結(jié)果更為接近（參見表 1）。文獻(xiàn)［9］中指出， 由于流體介質(zhì)為水， 存在負(fù)壓的卷吸作用，所以介質(zhì)為水時計算的阻力系數(shù)要大于實(shí)驗(yàn)值，這和文中得到仿真結(jié)果一致。

2.2卡門渦街現(xiàn)象及數(shù)值討論

在此雷諾數(shù)下， 尾流結(jié)構(gòu)出現(xiàn)了典型的卡門渦街， 如圖5所示。

圖6為一個周期內(nèi)的渦量圖， 由圖中可以看出， 邊界層分離區(qū)形成的渦旋周期性的交替從柱體上下側(cè)脫落， 在尾流區(qū)形成反向旋轉(zhuǎn)的渦對。

圖7為沿流場中軸線y=0上時均速度分布情況。從圖中可以看出， 回流區(qū)長度l/D=0.45（D為圓柱彈丸直徑）。當(dāng)l/D=0.35時回流速度達(dá)到最大，最大回流速度為umin/u= 0.16。圖中， 橫坐標(biāo)（x/D）為中軸線y=0區(qū)間段（0～800 mm）。

表1　計算結(jié)果與文獻(xiàn)對比Table 1　Comparison between calculation results and data from literature

圖5　卡門渦街現(xiàn)象圖Fig. 5　Diagram of Karman vortex street

圖6　一個周期內(nèi)的渦量圖Fig. 6　Vorticity distributions in one period

圖7　中軸線y=0上時均速度分布Fig. 7　Mean velocity distribution at y=0 on the medium axis

圖8為流場中沿x方向（順?biāo)鳎┠硯讉€豎直截面上時均流向速度。由圖中可以看出， 在圓柱近尾流區(qū)， x方向時均速度程V型， 且u/u0振幅很大， 隨著截面逐漸遠(yuǎn)離近尾流區(qū)， V型底處逐漸變緩， 當(dāng)x/D=3（x=300 mm）時， x方向時均速度基本為一條曲線， 說明距離圓柱越遠(yuǎn)， 流速分布受圓柱影響越小， 反之越大。

圖8　順?biāo)鞣较虿煌孛嫔蠒r均流向速度Fig. 8　Mean velocities at different sections along flow direction

圖9為圓柱后x=D方向上垂直于水流方向截面上時均流向速度。

圖9　x=100方向垂直于水流方向截面上時均流向速度Fig. 9　Mean velocity along the direction of x =100 perpendicular to flow

圖 10為表征圓柱繞流中的圓柱時均穩(wěn)態(tài)壓力系數(shù)Cp，且

時均壓力穩(wěn)態(tài)系數(shù)沿圓柱表面分布， 由圖10可以看出， 穩(wěn)態(tài)壓力沿彈丸周向?qū)ΨQ分布， 正對來流處的駐點(diǎn) Cp為 1， 隨著流速的恢復(fù)， 壓力系數(shù)值減少， 在θ為72°和288°時達(dá)到最小， 隨后壓力回升， 在彈丸背面回流區(qū)形成一個較為均勻的低壓區(qū)。

圖10　穩(wěn)態(tài)壓力系數(shù)周向分布Fig. 10　Circumferential distribution of steady pressure coefficient

2.3不同來流速度下的Cd和Str

以上文模型和設(shè)置研究不同來流速度對 Cd和Str的影響， 實(shí)驗(yàn)安排參見表2。

表2　實(shí)驗(yàn)安排表Table 2　Experiment schedule

數(shù)值仿真得到的圓柱彈丸Cd和Str隨來流速度變化規(guī)律見圖11和圖12。由圖中可知， Cd隨著來流速度的增大而減小， Str隨著來流速度的增大而增大。與文獻(xiàn)［1］中實(shí)驗(yàn)結(jié)果的變化趨勢圖相吻合。

圖11　不同來流速度下圓柱Cd示意圖Fig. 11　Schematic of cylinder Cdat different fllow rate

圖12　不同來流速度下圓柱Str示意圖Fig. 12　Schematic of cylinder Strat different flow rate

2.4柱體繞流渦街

圖13　不同來流速度下圓柱體流場渦量圖Fig. 13　Vortex flow field around cylinder at different flow rate

改變來流速度所得圓柱體繞流的流場渦量圖見圖13。x， y坐標(biāo)表示流場范圍， 單位為mm。由圖中可知， 當(dāng)來流速度處于25 m/s～200 m/s時，下游尾流區(qū)形成了規(guī)則的渦街脫落， 隨著來流速度的增大， 渦街保持規(guī)則的距離增大， 渦旋在y方向的尺度先增大后減小。

3　結(jié)束語

文中針對特定水下彈道工況下不同來速的彈丸柱體繞流進(jìn)行了 2D數(shù)值實(shí)驗(yàn)研究， 得到如下結(jié)論:

1） 采用 DES可以準(zhǔn)確科學(xué)的獲得圓形彈體水中流場特性， 通過與文獻(xiàn)中數(shù)據(jù)相比， DES所得到的結(jié)果與文獻(xiàn)中數(shù)據(jù)高度吻合；

2） 圓形彈丸繞流， 來流速度從25 m/s變化到200 m/s， Cd隨著來速的增加而減少， Str隨著來流速度的增大而增大。

研究結(jié)果表明， 在一定條件下， 增加彈丸在水中的速度， 可以減少彈丸在水中的阻力， 但此時隨速度的增加， 會提高彈丸水中渦脫落的頻率，需注意彈丸脫落不能接近其固有頻率， 以免引起共振。文中的研究結(jié)果對于彈丸設(shè)計有一定的參考價值。不足之處在于文中采用 2D模型替代圓柱彈丸， 具有一定局限性， 后期將開展 3D圓柱彈丸水中繞流研究。

［1］ C史里希廷. 邊界層理論［M］. 北京: 科學(xué)出版社， 1991.

［2］ Achenbach E. Distribution of Local Pressure and Skin Friction around a Circular Cylinder in Cross-Flow up to Re = 5×106［J］. Journal of Fluid Mechanics， 2000， 34（4）: 625-639.

［3］ 李燕玲， 蘇中地. 高雷諾數(shù)下單圓柱繞流的 DES三維數(shù)值模擬［J］. 中國計量學(xué)院學(xué)報， 2013， 24（4）: 364-369. Li Yan-ling， Su Zhong-di. 3D Mumerical Simulation of Flow Over a Circular Cylinder at High Reynolds Numbers Using DES Method［J］. Journal of China University of Metrology， 2013， 24（4）: 364-369.

［4］ Kravchenko A G， Moin P. Numerical Studies of Flow Over a Circular Cylinder at Re D=3900［J］. Physics of Fluids，2000， 12（2）: 403-417.

［5］ 鄧小兵. 不可壓縮湍流大渦模擬研究［D］. 綿陽: 中國空氣動力研究與發(fā)展中心， 2008.

［6］ Spalart P R， Jou W H， Strelets M， et al. Comments on the Feasibility of LES for Wings， and on Hybrid RANS/LES approach［C］. ［s.l.］: Advances in DNS/LES， 1997.

［7］ Breuer M. A Challenging Test Case for Large Eddy Simulation: High Reynolds Number Circular Cylinder Flow［J］. International Journal of Heat & Fluid Flow， 2000，21（5）: 648-654.

［8］ 郝鵬， 李國棟， 楊蘭， 等. 圓柱繞流流場結(jié)構(gòu)的大渦模擬研究［J］. 應(yīng)用力學(xué)學(xué)報， 2012， 29（4）: 437-443.

［9］ 詹昊， 李萬平， 方秦漢， 等. 不同雷諾數(shù)下圓柱繞流仿真計算［J］. 武漢理工大學(xué)學(xué)報， 2008， 30（12）: 129-132. Zhan Hao， Li Wan-ping， Fang Qin-han. Numerical Simulation of the Flow around a Circular Cylinder at Varies Reynolds Number［J］. Journal of Wuhan University of Technology， 2008， 30（12）: 129-132.

［10］ 祝志文. 高Re數(shù)圓柱繞流二維RANS模擬適用性分析［J］. 振動與沖擊， 2013， 32（7）: 98-101. Zhu Zhi-wen. Feasibility Analysis of 2D RANS Simulations for of Circular Cylinders Aevodynamics at High Re Number［J］. Journal of Vibratrion and Shock， 2013， 32（7）: 98-101.

（責(zé)任編輯: 楊力軍）

Simulation on Flow around Cylindrical Projectile in Water Based on DES model

ZHANG Chi-yu1，GUO Rui1，LIU Min2，LIU Rong-zhong1，ZHOU Hao1，CHEN Liang1
（1. School of Mechanical Engineering， Nanjing University of Science and Technology， Nanjing 210094， China； 2. Army Military Representative Office， Xi′an 710032， China）

Based on the software Fluent， the flow around a cylindrical projectile in water with high Reynolds number is simulated by using the detached eddy simulation（DES） method. The flow characteristics of the cylindrical projectile are analyzed at the flow rate of 25 m/s， and the lift coefficient， the drag coefficient and the Strouhal number， as well as the periodic alternatively shed Karmen vortex structure， are obtained. The mean velocity of flow field in different cross section and the mean steady pressure coefficient at cylinder circumference are analyzed. The image law is well consistent with that from the known large eddy simulation（LES） method. Further， the projectile drag coefficient and the Strouhal number at different speeds are analyzed， and the results show that with the increase in flow rate， the drag coefficient and the drag against the cylindrical projectile moving in water decrease， while the Strouhal number increases，inferring higher frequency of vortex shedding at both sides of the cylindrical projectile. When the vortex shedding frequency is close to the natural frequency of the projectile， resonance will occur to harm the projectile. This research may provide reference for design of projectile in underwater weapon warhead.

underwater warhead； flow around cylindrical projectile； detached eddy simulation（DES）； high Reynolds number； Karman vortex street

TJ630； TJ410.33

A

1673-1948（2016）05-0362-06

10.11993/j.issn.1673-1948.2016.05.009

2016-07-26；

2016-08-25.

高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金（20133219110019）.

張弛宇（1991-）， 男， 在讀碩士， 研究方向?yàn)閺椝幘_化與智能化.