基于溫度異重流模型的連續(xù)分層流數(shù)值模擬方法研究

牛明昌，丁勇，馬衛(wèi)狀，韓盼盼

（1.哈爾濱工程大學船舶工程學院，哈爾濱150001；2.上海交通大學海洋工程國家重點實驗室，上海200240）

基于溫度異重流模型的連續(xù)分層流數(shù)值模擬方法研究

牛明昌1，丁勇1，馬衛(wèi)狀1，韓盼盼2

（1.哈爾濱工程大學船舶工程學院，哈爾濱150001；2.上海交通大學海洋工程國家重點實驗室，上海200240）

文章基于溫度異重流模型，通過指定流場中溫度的分布及密度與溫度的函數(shù)關(guān)系實現(xiàn)密度的連續(xù)分層，并假設(shè)流場中的熱傳導在數(shù)值過程中可以忽略不計以確保密度分層的穩(wěn)定。分別采用RANS方法和LES方法對低傅汝德數(shù)和高傅汝德數(shù)情況下的圓柱繞流進行數(shù)值模擬，其尾流場中的lee波、葫蘆狀流線以及高傅汝德數(shù)情況下的波包與實驗值相符，低傅汝德數(shù)情況下圓柱上下游速度剖面曲線的最小值與實驗值僅存在2%的誤差，且其分布規(guī)律相同。通過改變流體的導熱系數(shù)，探討了Peclet數(shù)對該方法收斂穩(wěn)定性的影響，得到了收斂的數(shù)值條件。

分層流；溫度異重流；尾流場；大渦模擬

0 引言

由于太陽的熱輻射及大洋鹽熱環(huán)流等因素，海洋通常是分層的，其密度、溫度和鹽度沿垂直方向有明顯的變化。當潛體在密度分層流體中運動時，流體質(zhì)點會受到體積效應下的浮力作用，產(chǎn)生異于均勻環(huán)境的水動力學現(xiàn)象，如可以明顯觀察到流場中源生內(nèi)波的存在。連續(xù)分層作為海洋中典型的密度分層形式受到學者們格外的關(guān)注。Chomaz,Bonneton和Hopfinger[1]通過對連續(xù)分層流中的球體進行拖曳實驗，研究了近場尾流在不同傅汝德數(shù)下的形態(tài)；王進[2]在連續(xù)分層流中對三種不同長徑比拖曳潛體激發(fā)的內(nèi)波進行了實驗，得到了內(nèi)波轉(zhuǎn)捩傅汝德數(shù)經(jīng)驗公式；Long[3]研究了連續(xù)分層流中不同尺度障礙物對定常繞流控制方程求解的影響；Miles[4]在Long模型的基礎(chǔ)上用漸近分析的方法研究了不同傅汝德數(shù)條件下半圓柱繞流形成的lee波；Stevenson[5-6]通過對連續(xù)分層流中的圓柱進行拖曳實驗，先后研究了尾流場中的定常內(nèi)波和流體粘性對lee波的影響，其結(jié)果與Lighthill[7]的彌散波模型相吻合；Boyer[8-9]在連續(xù)分層流水槽中系統(tǒng)地研究了圓柱繞流的流動現(xiàn)象；丁勇[10]基于多相流混合模型研究了線性分層流中的圓柱繞流現(xiàn)象。目前，國內(nèi)外對連續(xù)分層流的研究主要集中在理論和實驗方面，尤其缺少相關(guān)的數(shù)值模擬研究。

本文通過指定流場中溫度的分布及密度與溫度的函數(shù)關(guān)系實現(xiàn)密度的連續(xù)分層，并假設(shè)數(shù)值過程中通過熱傳導改變的溫度值很小，建立了密度連續(xù)分層流的數(shù)值模擬方法。采用RANS方法和LES方法對其中的圓柱繞流現(xiàn)象進行數(shù)值模擬，探討研究了該方法獲得收斂結(jié)果的數(shù)值條件、收斂結(jié)果的穩(wěn)定性以及尾流場時間歷程的演變，將流場特征與Boyer[8]的實驗值進行對比，證實了此種方法在數(shù)值模擬密度連續(xù)分層流體時的可行性。

1 數(shù)值模型

1.1 溫度場假設(shè)

為了與流體動力學控制方程中的對流、擴散項保持一致，本文中的熱對流是指由于對流引起的溫度物理量輸運，與熱力學中熱對流的意義相異，熱傳導是指由于擴散引起的的溫度物理量輸運。

流體的密度主要由溫度所決定，為了保證密度垂向的穩(wěn)定分層，需保證溫度在垂向的分布是穩(wěn)定的，而根據(jù)熱力學第二定律，經(jīng)過足夠長時間之后，溫度在垂向的分布必然趨于均勻。因此本文假設(shè)通過熱傳導改變的溫度值在潛體尾跡演化的時間跨度內(nèi)足夠小，即熱傳導的影響可忽略不計，溫度在垂向上的分布短時間內(nèi)是穩(wěn)定不變的，如此便可實現(xiàn)穩(wěn)定的密度分層形式。此時，溫度輸運方程（能量方程）為對流占優(yōu)擴散方程，方程如下所示。

其中：T為溫度，k為流體的導熱系數(shù)，c為流體的比熱容，Γ=k/c為相對于T的廣義擴散系數(shù)。方程左邊為溫度場的局部變化率，右邊兩項分別為熱對流及熱擴散項。熱傳導可以忽略的程度在方程中體現(xiàn)為對流項及擴散項對局部場變化率影響程度的相對大小。

對流與擴散的強度之比可以用Pe（Peclet數(shù)）來衡量，其表達式如下所示：

其中：u為特征速度，L為特征長度，α=Γ/ρ為特征擴散系數(shù)。Pe數(shù)越大，對流占優(yōu)程度越高，那么必然存在一個臨界Peclet數(shù)Pecr，當Pe＞Pecr時，假設(shè)成立，流場中可以形成穩(wěn)定的密度分層。

1.2 模擬方法

整個流場的溫度、密度分布及計算模型示意圖如圖1所示。規(guī)定該圖的原點在圓柱圓心位置，x為流向，z為垂向，H為計算域的高度，U為遠方來流速度，g為重力加速度，ρ0和T0分別為流體域中的平均密度和平均溫度。為了與實驗保持一致，計算域的高度設(shè)置為0.2 m，圓柱的直徑設(shè)置為0.024 m，入口及出口邊界足夠遠。

文中用的無量綱參數(shù)：內(nèi)傅汝德數(shù)：Fr=U/Nd；雷諾數(shù)：Re=Ud/ν；無量綱時間：t′= tU/d，t在數(shù)值模擬條件下為計算時間，無量綱時間用tn′表示，在實驗條件下為繞流時間，無量綱時間用te′表示；N為浮頻率，N=（g△ρ/ρ0H）1/2；ν為流體的運動粘性系數(shù)，其它參數(shù)的含義如圖1所示。

圖1 計算模型示意圖Fig.1 Sketch of the calculation model

溫度T的垂向分布形式如（3）式所示。當T=300 K時，ρ=1 018.58 kg/m3；當T= 300.2 K時，ρ=998 kg/m3。由此可得整個流域溫度的垂向分布函數(shù)如（4）式所示。

此時，浮頻率的值為1，與實驗相符。值得注意的是，流體的運動粘性系數(shù)ν雖為流體的固有屬性，但是受溫度的影響很大。根據(jù)第10屆ITTC關(guān)于淡水和鹽水粘性系數(shù)的規(guī)定，在T=300 K附近，如果溫度變化5 K，粘性系數(shù)相差約0.07×10-6m2·s-1，本文流場的最大溫差為0.2 K，粘性系數(shù)在該范圍內(nèi)的變化非常小，因此可以將粘性系數(shù)視作常數(shù)。

Boyer的實驗表明，當Fr=0.4時，lee波已經(jīng)不是尾流場的主要特征，這時尾流場中出現(xiàn)明顯的渦結(jié)構(gòu)，所以規(guī)定當Fr＞0.4時為高傅汝德數(shù)，F(xiàn)r＜0.4時為低傅汝德數(shù)。對低傅汝德數(shù)下的情形，可通過RANS方法予以模擬實現(xiàn)，對高傅汝德數(shù)下的情形，則應通過LES方法實現(xiàn)。在高傅汝德數(shù)條件下，尾流場中湍流脈動成份的水動力作用會增強，而RANS方法抹去了瞬時脈動成份，LES方法則沒有這種弊端。理論研究[11-13]認為，當橫向距離大于πd時，利用大渦模擬進行計算，三維的影響就變得非常小，因此本文橫向設(shè)置為0.084 m（3.5d）。

采用O型網(wǎng)格對流體域進行空間離散，二維及三維情況下的網(wǎng)格數(shù)量分別為30萬和150萬。對于RANS方法，壓力離散格式采用“PRESTO！”，對流項離散格式采用“Second Order Upwind”，擴散項離散格式采用中心差分格式，時間離散格式采用“First Order Upwind”。對于LES方法，壓力離散格式采用“Standard”，對流項離散格式采用“Bounded Central Differencing”，擴散項采用中心差分格式，時間離散格式采用“Second Order Implicit”。同時兩種方法都是采用PISO算法進行求解，時間步長為0.005 s。

2 數(shù)值結(jié)果與討論

2.1 低傅汝德數(shù)條件下的數(shù)值結(jié)果

在低傅汝德數(shù)條件下，計算了三組工況，具體的參數(shù)如表1所示。

工況1條件下的流線圖如圖2所示，圖中（a）為數(shù)值結(jié)果，（b）為實驗結(jié)果。該圖表明低速條件下雖然在圓柱的下游形成一對與實驗值較為相符的lee波，但是下游的流線形狀并不能與實驗值完全相符。說明在低速條件下對流項的影響較弱，此時熱擴散并不能完全忽略不計。

工況2條件下結(jié)果穩(wěn)定時的流線圖及相同實驗條件下的流線圖如圖3所示，圖中（a）為數(shù)值結(jié)果，（b）為實驗結(jié)果。（a）圖中可以觀察到在圓柱的后方有兩對較明顯的lee波峰線，其特征與（b）圖的實驗結(jié)果相吻合。圖（a）、（b）在靠近x軸線處lee波的形狀與大小基本一致，不同的是數(shù)值結(jié)果中x軸線附近的流線上下兩側(cè)間距較大，實驗條件下x軸線附近上下兩側(cè)的流線在緊鄰圓柱及圓柱后方第一個波包后端幾近接觸。

表1 工況表（低傅汝德數(shù)）Tab.1 Conditions of calculation (low Froude number)

圖2 流線圖對比（Fr=0.018，Re=12）Fig.2 Comparison of streamline(Fr=0.018，Re=12)

圖3 流線圖對比（Fr=0.08，Re=60）Fig.3 Comparison of streamline(Fr=0.08，Re=60)

工況3條件下的流線圖及實驗條件下的流線圖如圖4所示，圖中（a）為數(shù)值結(jié)果，（b）為實驗結(jié)果。圖（a）表明緊鄰圓柱的尾流場流線呈葫蘆狀，與圖（b）中Boyer實驗得到的流線圖非常吻合，這是該階段尾流場的主要特征。另外，圖（a）中還存在斜向傳播的lee波，這種斜向傳播的lee波在實驗結(jié)果中也可以觀察到，只是lee波的位置略有差異，數(shù)值模擬過程中沒有觀察到轉(zhuǎn)子的出現(xiàn)。在該工況下，數(shù)值結(jié)果與實驗結(jié)果符合較好。

圖4 流線圖對比（Fr=0.17，Re=98.7）Fig.4 Comparison of streamline(Fr=0.17，Re=98.7)

2.2 高傅汝德數(shù)條件下的數(shù)值結(jié)果

在高傅汝德數(shù)條件下，計算了兩組工況，具體的參數(shù)如表2所示。

圖5為Fr=0.88，Re=480條件下流線圖的實驗結(jié)果及數(shù)值結(jié)果，圖（a）為實驗結(jié)果，圖（b）為數(shù)值結(jié)果。該圖表明數(shù)值和實驗條件下的尾流場流線形狀大致相同，在緊鄰圓柱的后方都有一對渦出現(xiàn)，在距離圓柱6倍直徑和16倍直徑的尾流場中存在著兩個波包，第一個波包下包絡(luò)著一對渦，圖（b）中第二個波包所在的位置比圖（a）中第二個波包出現(xiàn)的位置距離圓柱近，這是由于波包是一個不穩(wěn)定的尾流場特征，它在水平及垂向位置上不斷震蕩，震蕩幅度可以達到5倍的圓柱直徑。該數(shù)值方法基本模擬出了實驗工況下的結(jié)果。

表2 工況表（高傅汝德數(shù)）Tab.2 Conditions of calculation (high Froude number)

圖5 流線圖對比（Fr=0.88，Re=480）Fig.5 Comparison of streamline(Fr=0.88，Re=480)

圖6 渦量圖對比（Fr=1.77，Re=960）Fig.6 Comparison of vortex(Fr=1.77，Re=960)

表3 不同導熱系數(shù)下的工況表Tab.3 Conditions of calculation at different heat conductivity

圖7為Fr=0.018，Re=12時不同Pe數(shù)條件下圓柱繞流流線圖，圖（a）為Pe=7 212的情形，圖（b）為Pe=72 120的情形，兩幅圖中的流線明顯不同于Pe=72.12時的情形（圖2（a））。當Pe=7 212及Pe=72 120時，尾流場流線形狀與實驗結(jié)果更加相符（圖2（b））。圖8為Fr=0.08，Re=60時不同Pe數(shù)條件下圓柱繞流流線圖，圖（a）為Pe=32 054的情形，圖（b）為Pe=320 540的情形，兩幅圖也不同于Pe= 320.54的情形（圖3（a）），當Pe＞32 054時，尾流場流線形狀與實驗結(jié)果相符（圖3（b））。圖9為Fr= 0.17，Re=98.7時不同Pe數(shù)條件下圓柱繞流流線圖，圖（a）為Pe=68 114的情形，圖（b）為Pe=681 140的情形，結(jié)果與圖7、圖8的規(guī)律相符。但是當傅汝德數(shù)增大到一定程度時，Pe數(shù)的值已經(jīng)很大了，改變流體導熱系數(shù)對流場的影響變得非常有限。由于連續(xù)分層流中的圓柱繞流是一個準穩(wěn)態(tài)過程，流線隨時間變化，圖7～9各圖中（a）與（b）數(shù)值時刻均相同。

圖7 流線圖（Fr=0.018，Re=12）Fig.7 Diagram of streamline(Fr=0.018，Re=12)

圖8 流線圖（Fr=0.08，Re=60）Fig.8 Diagram of streamline(Fr=0.08，Re=60)

圖9 流線圖（Fr=0.17，Re=98.7）Fig.9 Diagram of streamline(Fr=0.17，Re=98.7)

分析表明，在同樣導熱系數(shù)條件下，隨著來流速度的增加，Pe數(shù)越來越大，對流占優(yōu)程度越高。雖然低傅汝德數(shù)工況下結(jié)果與實驗值不完全相符，但在高傅汝德數(shù)工況下能夠取得與實驗一致的結(jié)果。因此如果在低速條件下的數(shù)值結(jié)果能夠與實驗值吻合，那么高速條件下該方法理論上也能夠適用。

2.4 尾流場時間歷程的研究

在模擬過程中，數(shù)值時間并不能與實驗時間完全相對應，但在某數(shù)值時刻，若上游速度剖面曲線與某一實驗時刻結(jié)果相對應，同時刻下游速度剖面數(shù)值結(jié)果也應與同實驗時刻的結(jié)果相對應。圖10～13為Fr=0.018，Re=12時數(shù)值及實驗條件下上下游（x=±7.5d）速度剖面對比圖，圖10、圖11為Pe= 7 212時的情形，圖12、圖13為Pe=72 120時的情形。速度曲線的極值代表該剖面處速度的最大值或最小值，在上游速度剖面曲線吻合的時刻，下游速度剖面基本吻合，四組結(jié)果中波谷值的相對誤差最大為2%。另外，在Pe=7 212及Pe=72 120兩種情形下，實驗時刻te′=62對應的數(shù)值時刻均是實驗時刻te′=21對應數(shù)值時刻的3倍左右（實驗時刻對應的倍數(shù)也為3），在時間演化上，數(shù)值結(jié)果與實驗結(jié)果相符。

圖10 上下游速度剖面（Pe=7 212，tn′=7.88，te′=21)Fig.10 Velocity profiles of upstream and downstream(Pe=7 212,tn′=7.88,te′=21)

圖11 上下游速度剖面（Pe=7 212，tn′=24.06，te′=62)Fig.11 Velocity profiles of upstream and downstream(Pe=7 212,tn′=24.06,te′=62)

圖12 上下游速度剖面（Pe=72 120，tn′=6.64，te′=21)Fig.12 Velocity profiles of upstream and downstream(Pe=72 120,tn′=6.64,te′=21)

圖13 上下游速度剖面（Pe=72 120，tn′=18.9，te′=62)Fig.13 Velocity profiles of upstream and downstream(Pe=72 120,tn′=18.9,te′=62)

3 結(jié)論

基于溫度異重流模型，通過指定流場中溫度的分布及密度與溫度的函數(shù)關(guān)系實現(xiàn)了密度的連續(xù)分層，并假設(shè)流場中的熱傳導在數(shù)值過程中忽略不計以確保連續(xù)分層的穩(wěn)定性。分別采用RANS方法和LES方法數(shù)值模擬了低傅汝德數(shù)和高傅汝德數(shù)工況下的圓柱繞流，其尾流場特征與實驗相符合，證實了這種方法在數(shù)值模擬密度連續(xù)分層流時的可行性。具體結(jié)論如下：

（1）隨著傅汝德數(shù)的增大，Peclet數(shù)逐漸變大，能量方程變?yōu)閷α髡純?yōu)擴散方程，流場中的熱傳導可以忽略不計，密度分層能在較長時間內(nèi)保持穩(wěn)定，因此這種方法在數(shù)值模擬連續(xù)分層流時的穩(wěn)定性和精度都比較高。

（2）采用RANS方法成功模擬出了低傅汝德數(shù)條件下尾流場中的lee波，緊鄰圓柱的尾流場流線呈現(xiàn)葫蘆狀，其特征與實驗值相符。圓柱上下游（x=±7.5d）速度剖面曲線出現(xiàn)了一個極小值和兩個相等的極大值，其變化規(guī)律與實驗定量一致，峰值的誤差在2%以內(nèi)，且在時間歷程演變方面與實驗相符合。

（3）采用LES方法對連續(xù)分層流中高傅汝德數(shù)條件下的圓柱繞流進行了數(shù)值模擬。當Fr=0.88時，距離圓柱6倍直徑和16倍直徑的尾流場中分別出現(xiàn)了一個波包，且其在水平及垂向位置上不斷發(fā)生震蕩。當Fr=1.77時，尾流場進入了完全湍流狀態(tài)，圓柱后方出現(xiàn)了很強的渦結(jié)構(gòu)，由于分層的抑制作用，渦街在垂向的寬度保持不變。

（4）RANS方法能夠捕捉到流場中流線的變化，但對渦量的描述不夠細致，LES方法則彌補了這一弊端，能夠精確模擬出流場中的各種成分。因此，RANS方法適合于低傅汝德數(shù)下的流動，LES方法則同時適合低傅汝德數(shù)和高傅汝德數(shù)下的流動。

[1]Chomaz J M,Bonneton P,Hopfinger E J.The structure of the near wake of a sphere moving horizontally in a stratified fluid[J].Journal of Fluid Mechanics,1993,254:1-21.

[2]王進,尤云祥,胡天群,等.密度分層流體中不同長徑比拖曳潛體激發(fā)內(nèi)波特性試驗[J].科學通報,2012,08(57): 606-617. Wang Jin,You Yunxiang,Hu Tianqun,et al.The characteristics of internal waves excited by towed bodies with different aspect ratios in a stratified fluid[J].Chin Sci Bull,2012,08(57):606-617.

[3]Long R R.Some aspects of the flow of stratified fluids[J].Tellus,1955,7(3):341-357.

[4]Miles J W H.Lee waves in a stratified flow.Part 2.Semi-circular obstacle[J].Journal of Fluid Mechanics,1968,33(4): 803-814.

[5]Stevenson T N.Some two-dimensional internal waves in a stratified fluid[J].Journal of Fluid Mechanics,1968,33(04): 715-720.

[6]Stevenson T N,Chang W L,Laws P.Viscous effects in lee waves[J].Geophysical&Astrophysical Fluid Dynamics,1979, 13(1):141-151.

[7]Lighthill M J.On waves generated in dispersive systems to travelling forcing effects,with applications to the dynamics of rotating fluids[M].Hyperbolic Equations and Waves,Springer,1970:124-152.

[8]Boyer D L,Davies P A,Fernando H,et al.Linearly stratified flow past a horizontal circular cylinder[J].Philosophical Transactions of the Royal Society of London A:Mathematical,Physical and Engineering Sciences,1989,328(1601):501-528.

[9]Xu Y,Fernando H J,Boyer D L.Turbulent wakes of stratified flow past a cylinder[J].Physics of Fluids(1994-present), 1995,7(9):2243-2255.

[10]丁勇,韓盼盼,段菲,馬衛(wèi)狀.線性分層流中圓柱繞流數(shù)值模擬方法研究[J].哈爾濱工程大學學報,2016,9:1-5. Ding Yong,Han Panpan,Duan Fei,Ma Weizhuang.Numerical study of linearly stratified flow past a cylinder based on multiphase mixture model[J].Journal of Harbin Engineering University,2016,9:1-5.

[11]M B.Large eddy simulation of the subcritical flow past a circular cylinder:numerical and modeling aspects[J].International Journal for Numerical Methods in Fluids,1998,28(9):1281-1302.

[12]M B.Numerical and modeling influences on large eddy simulations for the flow past a circular cylinder[J].International Journal of Heat And Fluid Flow,1998,19(5):512-521.

[13]Kravchenko A G M P.Numerical studies of flow over a circular cylinder at ReD=3 900[J].Physics of Fluids,2000,12(2): 403-417.

Research on the numerical simulation methods of continuously stratified flows based on thermal density current model

NIU Ming-chang1,DING Yong1,Ma Wei-zhuang1,Han Pan-pan2
(1.College of Shipbuilding Engineering,Harbin Engineering University,Harbin 150001,China;2.State Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

Based on thermal density current model,the continuous stratification of density is achieved by specifying the distribution of temperature and regarding density as a function of temperature.It is assumed that the heat conduction is negligible to ensure the stability of density stratification.The RANS method and LES method are used respectively to simulate the flows past a cylinder at low Froude numbers and high Froude numbers.The lee waves,gourd-shaped streamlines and wave packets are consistent with the experimental results.For the minimum values of velocity profile curves,there are only 2%errors between the calculated and experimental values,and the curve shapes are the same.The influences of Peclet numbers on the stability of this method are discussed by changing the heat conductivity coefficient of fluid.Finally,the numerical conditions of convergence are obtained.

stratified fluid;thermal density current;wake field;large eddy simulation

U661.1

A

10.3969/j.issn.1007-7294.2017.08.002

1007-7294（2017）08-0941-09

2017-03-26

××減震降噪工程專項計劃

牛明昌（1993-），男，碩士研究生；丁勇（1959-），男，博士，教授，通訊作者，E-mail:dingyong@hrbeu.edu.cn。