柏勁松,王 濤,鄒立勇,李 平

(中國(guó)工程物理研究院流體物理研究所沖擊波物理與爆轟物理國(guó)防科技重點(diǎn)實(shí)驗(yàn)室,四川 綿陽(yáng)621900)

1 引 言

對(duì)于復(fù)雜湍流的研究,解析方法幾乎無能為力,實(shí)驗(yàn)測(cè)量和數(shù)值模擬是研究復(fù)雜湍流的必要手段。湍流的數(shù)值模擬主要包括直接數(shù)值模擬(direct numerical simulation,DNS)、雷諾平均數(shù)值模擬(Reynolds averaged Navier-Stokes,RA NS)和大渦數(shù)值模擬(large eddy simulation,LES)3 種。其中LES 是基于網(wǎng)格尺度封閉模型及對(duì)大尺度渦進(jìn)行直接求解NS 方程,可以模擬湍流發(fā)展過程。近年來,湍流大渦數(shù)值模擬普遍接受的觀念是:湍流統(tǒng)計(jì)特性能否與DNS 或?qū)嶒?yàn)結(jié)果符合良好是衡量大渦數(shù)值模擬準(zhǔn)確性的標(biāo)準(zhǔn),而不是追求它與DNS 結(jié)果的相關(guān)性[1-2]。大尺度脈動(dòng)到小尺度脈動(dòng)的能量傳遞是大渦數(shù)值模擬的關(guān)鍵量,只要亞格子模型正確模擬能量傳遞,大尺度脈動(dòng)的統(tǒng)計(jì)量就基本上能夠計(jì)算正確。

亞格子應(yīng)力模型是湍流大渦模擬的核心問題,Smagorisky 模型就是參照雷諾平均模式最早提出的唯象渦粘模型,至今仍在使用。渦粘模型的最大優(yōu)點(diǎn)是簡(jiǎn)單,調(diào)整模型系數(shù)能夠保證模型的亞格子耗散與實(shí)際亞格子耗散一致,但是Smagorisky 模型屬于耗散型,在流場(chǎng)中任意一點(diǎn)都是從可解尺度湍流向不可解尺度湍流輸送能量,即所謂正轉(zhuǎn),而不存在相反方向的能量傳遞,在實(shí)際復(fù)雜湍流中,已經(jīng)發(fā)現(xiàn)可能存在局部正轉(zhuǎn)。Smagorisky 模型總體上耗散過大,因此后來又出現(xiàn)了動(dòng)力模型[3]、理性亞格子模型[4]等。另一種渦粘模式亞格子應(yīng)力模型叫Vreman 模型,由A.W.Vreman[5]提出,該模型的耗散較S magorisky 模型小,在低雷諾數(shù)湍流混合層計(jì)算中的計(jì)算結(jié)果與動(dòng)力模型基本一致,更接近于DNS 的計(jì)算值。Vreman 模型保持了Smagorisky 模型的優(yōu)點(diǎn),在湍流混合層和槽道湍流大渦模擬中得到推廣使用[6]。

本文中以可壓縮多介質(zhì)粘性流體界面不穩(wěn)定性誘發(fā)湍流混合研究為目的,在M FPPM[7](multi-fluid piecew ise parabolic method)和M VPPM[8]基礎(chǔ)上,利用Smagorinsky 和V reman 亞格子湍流模型,采用大渦數(shù)值模擬方法求解可壓縮流體NS 方程,通過算子分裂分步計(jì)算,以多介質(zhì)流體高精度(parabolic piecew ise method)方法計(jì)算無粘流場(chǎng),二階空間中心差方法和兩步Rung-Kutta 時(shí)間推進(jìn)相結(jié)合的方法計(jì)算分子動(dòng)力學(xué)粘性、湍流粘性以及熱流部分對(duì)流場(chǎng)的影響,研制了適用于可壓縮多介質(zhì)流體界面不穩(wěn)定性發(fā)展演化至湍流階段的計(jì)算方法和二維計(jì)算程序M VFT。文中在2 種亞格子湍流模型下采用大渦模擬方法計(jì)算了在LANL 激波管上進(jìn)行的單氣柱實(shí)驗(yàn)[9],給出了不同時(shí)刻氣柱的高度和寬度,流場(chǎng)中SF6上游邊界、下游邊界和渦邊界處的速度以及渦的特征。通過與LAN L 實(shí)驗(yàn)和計(jì)算結(jié)果的比較可見,在單氣柱激波管界面不穩(wěn)定性誘發(fā)至湍流混合大渦數(shù)值模擬研究中,Vreman 模型和Smagorinsky 模型計(jì)算的氣柱形狀基本一致,而在流場(chǎng)速度計(jì)算中Vreman 模型的計(jì)算結(jié)果與實(shí)驗(yàn)更接近一些,比較2 種模型計(jì)算渦的分布,V reman 模型耗散小于Smagorinsky 模型。通過比較和分析初步確認(rèn)MVFT 方法和計(jì)算程序可用于對(duì)界面不穩(wěn)定性誘發(fā)湍流混合流動(dòng)的大渦數(shù)值模擬。

2 控制方程和計(jì)算方法

2.1 基本控制方程

考慮熱傳導(dǎo)和粘性情況下的大渦模擬控制方程采用笛卡爾坐標(biāo)系下的NS 方程,經(jīng)過Faver 濾波得到,略去非線性項(xiàng)后的基本形式為

式中:i、j 分別代表x、y、z 等3 個(gè)方向,相同的i、j 表示求和分別為密度、速度和壓力,為單位質(zhì)量的總能量,N 為介質(zhì)的種類,φs 為第s 種介質(zhì)的體積分?jǐn)?shù),在含有N 種不同介質(zhì)的混合網(wǎng)格中,各體積分?jǐn)?shù)滿足。狀態(tài)方程采用p=(γ-1)ρe+γπ形式,對(duì)于氣體,γ為比熱比(γ＞1),π=0,對(duì)于液體和用沖擊Hugoniot 曲線描述的彈性密實(shí)介質(zhì),γ為材料的擬合常數(shù),π為具有粘性張量量綱的材料常數(shù)。為有效粘性應(yīng)力張量,為湍流有效粘性系數(shù)為流體動(dòng)力粘性系數(shù),為亞格子粘性系數(shù)為熱傳導(dǎo)在單位時(shí)間單位空間的能量流為流體溫度為有效導(dǎo)熱系數(shù),Prlam+μsgscp/Prsgs,λlam為流體導(dǎo)熱系數(shù),λsgs為亞格子導(dǎo)熱系數(shù),cp為比定壓熱容,Prlam和Prsgs為普朗特?cái)?shù)為熱導(dǎo)率。

2.2 2 種亞格子應(yīng)力模型

對(duì)于Smagorinsky 模型,其亞格子粘性系數(shù)

對(duì)于V reman 模型,其亞格子粘性系數(shù)

2.3 計(jì)算方法

采用算子分裂技術(shù),將方程組(1)描述的物理過程分解為3 個(gè)子過程進(jìn)行計(jì)算,即整個(gè)流量的計(jì)算分解為無粘流量、粘性流量和熱流量3 部分。那么方程組(1)就可以分裂為如下的2 個(gè)方程組

對(duì)于無粘流量部分包括對(duì)沖擊波、稀疏波及接觸間斷的計(jì)算,采用PPM 方法求解方程組(4),對(duì)于多物質(zhì)界面,由于不穩(wěn)定性后期界面混亂甚至破碎,因此無法進(jìn)行界面重構(gòu),采用VOF(volume of fluid)追蹤算法來捕捉,根據(jù)介質(zhì)的體積分?jǐn)?shù)區(qū)分界面附近的不同材料,采用壓力一致、速度相同基本假定確定混合網(wǎng)格中的計(jì)算參量,具體推導(dǎo)過程見文獻(xiàn)[8]。對(duì)于粘性流量和熱流量的計(jì)算,考慮湍流的亞格子應(yīng)力、牛頓流體粘性應(yīng)力和能量流的影響,所以在求解方程組(4)的基礎(chǔ)上,還需要求解方程組(5)。由于亞格子應(yīng)力、牛頓流體粘性應(yīng)力和能量流僅對(duì)流場(chǎng)的動(dòng)量、能量和界面擴(kuò)散產(chǎn)生影響,因此方程組(5)中第1 個(gè)方程描述的質(zhì)量保持守恒,離散求解時(shí)可以不予考慮,方程組(5)的第2 ～4 個(gè)方程采用二階空間中心差分方法計(jì)算,然后再利用兩步Rung-Kutta 方法進(jìn)行時(shí)間推進(jìn),求解過程見文獻(xiàn)[11]。

3 LANL 激波加載氣柱不穩(wěn)定性實(shí)驗(yàn)?zāi)Ｐ痛鬁u模擬

C.A.Zoldi[9]在LANL 激波管上進(jìn)行激波加載重氣柱實(shí)驗(yàn),實(shí)驗(yàn)測(cè)試段的布局如圖1,1.2Ma 的空氣沖擊波作用于SF6氣柱,由于在空氣和SF6氣柱界面處存在壓力和密度梯度,因此在激波作用下的RM 界面不穩(wěn)定性發(fā)展演化產(chǎn)生渦結(jié)構(gòu),最后形成空氣和SF 6 氣體的混合。渦通過速度的旋度進(jìn)行計(jì)算,即ω=▽×v。C.A.Zoldi[9]采用平面激光誘導(dǎo)熒光技術(shù)PLIF(planar laser induced fluorescence)給出了實(shí)驗(yàn)觀察到的6 個(gè)時(shí)刻SF6 氣柱形狀,如圖2(a)所示,運(yùn)用粒子圖像速度儀PIV(particle image velocimetry)技術(shù)給出了SF6氣柱3 個(gè)測(cè)試部位的速度,并應(yīng)用RAGE 計(jì)算程序進(jìn)行數(shù)值模擬,給出了對(duì)應(yīng)時(shí)刻的計(jì)算圖形,如圖2(b)所示。RAGE 是求解Euler 方程的計(jì)算程序,忽略了氣體的粘性和熱效應(yīng),未考慮湍流模型。本文中根據(jù)第2 節(jié)方法研制了用于求解NS 方程的多介質(zhì)粘性流體和湍流大渦模擬計(jì)算程序M VFT,既考慮了氣體的物理粘性和熱效應(yīng),還加入湍流模型模擬湍流混合過程,從物理上更適用于模擬激波加載重氣柱至混合實(shí)驗(yàn)?zāi)Ｐ汀?/p>

圖1 激波管測(cè)試段示意圖Fig.1 Test section of the shock tube

圖2 實(shí)驗(yàn)動(dòng)態(tài)圖像和RAGE 模擬結(jié)果Fig.2 Experimental dynamic images and RAGE simulations at different times

計(jì)算采用的空氣和SF6氣體參數(shù)如表1 所示,普朗特?cái)?shù)取為0.67 和0.7。在數(shù)值模擬中考慮到SF 6 氣柱在空氣中的耗散,在SF 6 和空氣之間設(shè)置了一個(gè)高斯函數(shù)型的耗散界面過渡層ITL(interface transition layer),如圖3 所示,過渡層內(nèi)的密度運(yùn)用高斯函數(shù)計(jì)算,其中Rc=0.12 cm,δ=0.275 cm,計(jì)算區(qū)域網(wǎng)格為600×100。計(jì)算采用Smagorinsky 和Vreman 亞格子湍流模型,分別給出了模型參數(shù)CS=0.17 和CV=0.07 時(shí)SF6 氣柱各對(duì)應(yīng)時(shí)刻的密度、渦等位圖,如圖4、5 所示。定性比較2 模型給出的密度和渦圖形,在早期區(qū)別較小,只是在后期圖形出現(xiàn)明顯差異,在750 μs時(shí)刻V reman 模型給出的SF 6 氣柱形狀與相同時(shí)刻實(shí)驗(yàn)圖形更接近一些,為了定量比較2 種亞格子應(yīng)力模型計(jì)算結(jié)果以及與實(shí)驗(yàn)和RAGE 計(jì)算結(jié)果,圖6 中給出了SF 6 氣柱在幾個(gè)時(shí)刻的寬度l 和高度h的比較,從圖6 可見,RAGE 的計(jì)算結(jié)果偏離實(shí)驗(yàn)最大,M VFT 所用的2 個(gè)亞格子模型計(jì)算結(jié)果基本一致,均較RAGE 接近實(shí)驗(yàn)值,但實(shí)驗(yàn)值均大于2 程序的計(jì)算結(jié)果;對(duì)于流場(chǎng)中SF6 氣柱運(yùn)動(dòng)速度的比較,給出了圖7 中定義的U E(upstream edge)、DE(dow nstream edge)和VE(vortex edge)3 個(gè)位置的位移歷史,圖7(a)、(b)分別為Smagorinsky 模型和Vreman 模型的計(jì)算結(jié)果,表2 列出了3 個(gè)位置平均運(yùn)動(dòng)速度vU E、vDE、vV E的M VFT 計(jì)算值、實(shí)驗(yàn)值和RAGE 程序計(jì)算值。通過比較可見,M VFT 計(jì)算中Vreman 模型(VM)給出的速度值大于Smagorinsky 模型(SM),但均小于RAGE 程序的計(jì)算值,MVFT 和RAGE 的計(jì)算值均遠(yuǎn)小于實(shí)驗(yàn)值。

表1 空氣和SF6 計(jì)算參數(shù)[9]Table 1 Computational parameters for air and SF6[9]

對(duì)于計(jì)算和實(shí)驗(yàn)差異,C.A.Zoldi[9]認(rèn)為實(shí)驗(yàn)中空氣沖擊波馬赫數(shù)略大于1.2 是造成實(shí)驗(yàn)值偏大的原因之一,另外如果實(shí)驗(yàn)所用的SF6 氣柱密度小于計(jì)算所取的SF6 氣柱密度也會(huì)導(dǎo)致實(shí)驗(yàn)值偏大;由于MVFT 考慮了氣體的物理粘性、熱效應(yīng)以及湍流模型,因此其計(jì)算值小于RAG E 的計(jì)算值是合理的;對(duì)于MV FT 所用的2 個(gè)亞格子模型計(jì)算的差異,則是由于V reman 模型的耗散小于Smagorinsky模型所引起,通過圖8 比較470 μs 時(shí)刻經(jīng)過最大渦中心沿x 方向的分布圖可以進(jìn)一步證實(shí),前者的渦峰大于后者,最左邊的負(fù)峰值和正峰值則分別對(duì)應(yīng)界面過渡層中上游邊界和下游邊界位置。

表2 運(yùn)動(dòng)氣柱速度Table 2 Velocities of the evolving cylinder

圖3 高斯函數(shù)型耗散界面過渡層Fig.3 Gaussian interfacial transition layer

圖4 不同時(shí)刻MVFT 計(jì)算的密度等位圖Fig.4 Density contour images by M VFT at different times

圖5 不同時(shí)刻M VFT 計(jì)算的渦等位圖Fig.5 Vortex contour images by M VFT at different times

圖6 SF 6 氣柱寬度和高度的實(shí)驗(yàn)和計(jì)算結(jié)果Fig.6 Width and height measurements and calculations of the SF 6 evloving cylinder

圖7 SF6 氣柱上游、下游和渦邊界的位置Fig.7 Location and displacement of the upstream edge(UE),the downstream edge(DE),and the vortex edge(VE)on the SF6 evolving cylinder

圖8 470 μs 時(shí)刻過最大渦中心沿x 方向分布圖Fig.8 Vorticity profile alone the x axis through the maximum vorticity at t=470 μs

4 結(jié) 論

(1)將Smagorinsky 和V reman 亞格子模型引入到可壓縮多介質(zhì)粘性流體動(dòng)力學(xué)高精度計(jì)算程序M VPPM 中,實(shí)現(xiàn)了RM 界面不穩(wěn)定性后期湍流混合的大渦數(shù)值模擬,研制了適用于可壓縮多介質(zhì)流體界面不穩(wěn)定性發(fā)展演化至湍流階段的計(jì)算方法和二維計(jì)算程序M VFT;(2)分析了Smagorinsky 和V reman 渦粘型亞格子模型對(duì)能量的耗散,通過對(duì)LA NL 激波管單氣柱RM不穩(wěn)定性實(shí)驗(yàn)中氣柱上下游邊界、渦邊界速度的定量計(jì)算及與實(shí)驗(yàn)和RAGE 程序計(jì)算結(jié)果比較可見,V reman亞格子模型對(duì)能量的耗散小于Smagorinsky 模型;(3)RM 界面不穩(wěn)定性后期湍流混合大渦數(shù)值模擬還存在一些不確定因素,Smagorinsky 模型和Vreman 模型不包含轉(zhuǎn)捩機(jī)制,模型選取和模型參數(shù)選擇具有不唯一性等。

[1] Jimenez J,Moser R D.Large-eddy simulation:Where are w e and w hat can we expect?[J].AIAA Journal,2000,38(4):605-612.

[2] Pope S.Ten questions concerning the large eddy simulation of turbulent flow s[J].New Journal of Physics,2004,6(1):1-24.

[3] 張兆順,崔桂香,許春曉.湍流理論與模擬[M].北京:清華大學(xué)出版社,2005.

[4] Cui G X, Zhou H B,Zhang Z S.A new dynamic subgrid eddy viscosity model w ith application to tubulent channel flow[J].Physics of Fluids,2004,16(8):2835-2942.

[5] Vreman A W.An eddy-viscosity subgrid-scale model for turbulent shear flow:Algebraic theory and applications[J].Physics of Fluids, 2004,16(10):3670-3681.

[6] Park N, Lee S, Lee J, et al.A dynamic subgrid-scale eddy viscosity model w ith a global model coefficient[J].Physics of Fluids,2006,18(2):125109.

[7] Bai J S, Li P,Zhang Z J,et al.Application of the high-resolution Godunov method to the multifluid flow calculations[J].Chinese Physics, 2004,13(12):1992-1997.

[8] 柏勁松,李平,王濤,等.可壓縮多介質(zhì)粘性流體的數(shù)值計(jì)算[J].爆炸與沖擊,2007,27(6):515-521.BAI Jing-song, LI Ping,WANG Tao,et al.Computation of compressible multi-viscosity-fluid flow s[J].Explosion and Shock Waves, 2007,27(6):515-521.

[9] Zoldi C A.A numerical and experimental study of a shock accelerated heavy gas cylinder[D].New York:Applied Mathematics and Statistics, State University of New York at Stony Brook,2002.

[10] Lilly D K.The representation of small-scale turbulence in numerical simulation experiments[C]∥Goldstine H H.Proceedings of IBM Scientific Com puting Symposium on Environmental Sciences.Yorktown Heights,New York,1967:195-210.

[11] 柏勁松,李平,鄒立勇,等.界面不穩(wěn)定性引起混合過程的二維數(shù)值計(jì)算[J].力學(xué)學(xué)報(bào),2008,40(4):464-471.BAI Jing-song, LI Ping,ZOU Li-yong,et al.Two-dimensional simulation of the mixing induced by interface instabilities[J].Chinese Journal of Theoretical and Applied Mechanics,2008,40(4):464-471.