基于大渦模擬研究分形柵格激發(fā)湍流的傳熱特性

段成棟，魏超(中國礦業(yè)大學(xué)徐海學(xué)院，江蘇 徐州 221008)

0 引言

不同湍流狀態(tài)下的流動(dòng)對換熱效果影響顯著，例如：湍流強(qiáng)度、壓降、均勻各向同性等湍流參數(shù)對于換熱效果具有重要作用。湍流對圓柱體傳熱的影響在很多工程應(yīng)用中有重要意義[1]。在實(shí)際工程中，通常使用圓柱形幾何體在流體和壁之間進(jìn)行熱交換，例如管殼式熱交換器、壓水反應(yīng)堆、水-空氣散熱器等[2]。

1935年，G.I Taylor在風(fēng)洞中放置了多排柵格，當(dāng)均勻氣流垂直流過柵格時(shí)，會(huì)發(fā)生不規(guī)則擾動(dòng)。這種不規(guī)則擾動(dòng)向下游移動(dòng)過程中，由于沒有外部干擾，逐漸演化為均勻各向同性湍流。相似的均勻各向同性湍流已經(jīng)在風(fēng)洞或水洞中得到了廣泛的研究。2007年，Hurst等[3]使用各種多尺度(分形)柵格在風(fēng)洞中產(chǎn)生湍流，發(fā)現(xiàn)復(fù)雜的多尺度柵格結(jié)構(gòu)會(huì)極大地影響湍流行為，并通過實(shí)驗(yàn)方法探索了分形維數(shù)Df≤2、柵格有效尺度Meff和網(wǎng)格厚度比tr=tmax/tmin。這三個(gè)因素對分形柵格激發(fā)湍流的流動(dòng)機(jī)制、各向同性、均勻性和衰減特性有很大影響。

2016年，G.Melina等[4]研究了不同幾何形狀柵格激發(fā)的湍流場中心線上的流動(dòng)：分形柵格(FSG17)、單一方形柵格(SSG)、規(guī)則柵格(RG60)。結(jié)果表明，在湍流產(chǎn)生區(qū)，單一柵格結(jié)構(gòu)產(chǎn)生的湍流強(qiáng)度大于其他兩種情況，這是由于較強(qiáng)的渦旋脫落效應(yīng)，而分形柵格抑制了渦旋脫落。2015年Torrano等[5]將實(shí)驗(yàn)測量數(shù)據(jù)與基于RANS的柵格生成湍流衰減數(shù)值研究進(jìn)行了比較。結(jié)果表明，基于RANS的數(shù)值模擬捕捉湍流流場大尺度特征的能力在可接受的精度范圍內(nèi)。R.J.Hearst等[6]設(shè)計(jì)了一種特殊的分型柵格結(jié)構(gòu)，還增加了分形柵格實(shí)驗(yàn)的下游測量范圍。結(jié)果表明，特殊柵格產(chǎn)生的湍流行為與之前分形柵格附近處觀察到的類似，遠(yuǎn)下游湍流性能則與規(guī)則柵格產(chǎn)生的類似。

上述結(jié)論為研究湍流參數(shù)對傳熱效果的影響提供了理論依據(jù)和實(shí)驗(yàn)依據(jù)。Smith等[7]提出了一個(gè)理論模型，假設(shè)在前駐點(diǎn)附近，雷諾應(yīng)力與自由流中的湍流強(qiáng)度Tu和與壁面的距離成正比。他們的模型得到了實(shí)驗(yàn)結(jié)果的支持，實(shí)驗(yàn)結(jié)論表明前駐點(diǎn)處的弗羅斯林?jǐn)?shù)與湍流參數(shù)TuRe0.5成正比，其中Re是基于圓柱體直徑D和平均流向速度的雷諾數(shù)。Kestin等[8]以及Lowery等[9]分別測量了風(fēng)洞中柵格產(chǎn)生的湍流場中圓柱體的傳質(zhì)和傳熱。這兩項(xiàng)研究都將NuFsp/Re0.5的值與TuRe0.5的二次多項(xiàng)式函數(shù)相關(guān)聯(lián)。此外，幾項(xiàng)實(shí)驗(yàn)表明，流動(dòng)中的積分長度標(biāo)度Lu也非常重要[10]。Van Der Hegge-Zijnen[11]報(bào)道，在Re和Nu相同條件下，當(dāng)0＜LU/D＜1.6時(shí)，Nu隨LU/D增加；當(dāng)LU/D＞1.6時(shí)，Nu隨LU/D減小，其中“L”定義為湍流的長度尺度。最近，為了強(qiáng)化傳熱，人們廣泛研究了流場中的渦旋結(jié)構(gòu)，并通過在流場中添加渦流發(fā)生器實(shí)現(xiàn)了良好的傳熱效果。

隨著計(jì)算機(jī)技術(shù)的進(jìn)步，數(shù)值模擬已成為研究湍流的重要手段。常用的數(shù)值模擬方法有三種：直接數(shù)值模擬(DNS)、雷諾平均Navier-Stokes方程模擬(RANS)、大渦模擬(LES)。DNS可以獲得湍流統(tǒng)計(jì)的深度數(shù)據(jù)，包括壁面附近的高階量，而它只能用于求解非常有限類型的湍流。RANS能以最經(jīng)濟(jì)的方法正確預(yù)測高雷諾數(shù)的復(fù)雜流動(dòng)，但結(jié)果是時(shí)間平均的，不能反映流場的瞬時(shí)信息。這類模型的典型例子是k-ε或k-ω模型。最近，許多雷諾應(yīng)力模型都是通過修正湍流的壓力-應(yīng)變關(guān)系發(fā)展起來的[12-15]。LES直接預(yù)測大尺度渦旋運(yùn)動(dòng)，對較小的尺度(亞格子尺度)和更接近各向同性的渦旋[16]，采用亞格子尺度(SGS)模型。不僅可以獲得詳細(xì)的流場動(dòng)態(tài)信息，而且可以節(jié)省計(jì)算成本。因此，大渦模擬得到了越來越廣泛的應(yīng)用。

鑒于上述分析，文章采用大渦模擬的方法研究了三種不同截面柵格產(chǎn)生的湍流以及恒溫圓柱在流場中的傳熱。三種湍流發(fā)生裝置的橫截面形狀不同：三角形、倒三角形和正方形。通過分析流場中速度、湍流強(qiáng)度、壓降和努塞爾數(shù)等瞬時(shí)和統(tǒng)計(jì)量，揭示了不同擾動(dòng)下流體運(yùn)動(dòng)和傳熱機(jī)理，以及湍流特性與傳熱之間的關(guān)系。

2 數(shù)值模擬方法

2.1 計(jì)算域與網(wǎng)格劃分

基于LES研究流體通過分形柵格產(chǎn)生的不同流場對圓柱體產(chǎn)生不同的熱交換效應(yīng)。如圖1所示，具體的模擬方案如下：在長方體風(fēng)洞中距離入口0.1 m處放置三種不同形狀的分形柵格，第四種不放置任何擾動(dòng)作為參照。在所有工況中的同一位置各放置一個(gè)直徑為D=0.05 m、長度為0.2 m、溫度為323 K的恒溫圓柱體。為了消除尺寸因子的不可比性，計(jì)算域基于換熱圓柱直徑進(jìn)行無量綱化。無量綱后的計(jì)算域?yàn)長x×Ly×Lz=24D×4D×4D。

圖1 計(jì)算域

分形柵格的特征由以下參數(shù)定義：

—分形結(jié)構(gòu)的維數(shù)N，此研究中N=3；

—某一維度分形圖案的數(shù)量4j；

—柵格長度Lj=0.5jL0和橫向厚度tj=0.5jt0，L0=0.5Ly，t0=tmax=0.01 m；

—厚度比tr/tmax/tmin=4，即形成最大正方形的橫向厚度與最小正方形的橫向厚度之間的比率；

—三棱柱柵格的角度α=53°。

文章重點(diǎn)研究了三棱柱分形網(wǎng)格中湍流的演化規(guī)律。

考慮到三角柱相交處的結(jié)構(gòu)網(wǎng)格質(zhì)量難以保證，在工況1和工況2的整個(gè)計(jì)算域中使用非結(jié)構(gòu)網(wǎng)格。為了保證更高的計(jì)算精度，設(shè)置y+=Δy×u∞/ν=1，因此，第一層網(wǎng)格厚度只有0.005 6 mm，y+是從第一層網(wǎng)格的質(zhì)心到壁面的無量綱距離。在工況1和工況2中，使用MESH生成四面體非結(jié)構(gòu)化網(wǎng)格，選擇曲率尺寸函數(shù)來控制曲率較大的區(qū)域和網(wǎng)格生長，以及表面相鄰區(qū)域的分布。將最小尺寸、最大面尺寸、最大尺寸分別設(shè)置為2.0×10-5m、0.002 m和0.005 m?？偣?jié)點(diǎn)數(shù)為1 629 407，總網(wǎng)格數(shù)為8 917 056。

由于工況3和工況4中的幾何模型較為簡單，因此工況3和工況4中使用了在流向、法向方向和展向方向上網(wǎng)格數(shù)為624×200×160的結(jié)構(gòu)化網(wǎng)格。

2.2 控制方程

LES的基本思想是通過過濾將大尺度渦旋與小尺度渦旋分離。為了進(jìn)行湍流計(jì)算，連續(xù)性方程和Navier-Stokes方程過濾如下：

式(1)～(3)中：ρ為流體密度；ν為流體的動(dòng)力黏度；為速度的過濾值；為壓力的過濾值；τij為子網(wǎng)格應(yīng)力項(xiàng)(SGS)，反映了小尺度渦旋對大尺度渦旋的影響。

為了使過濾后的方程封閉，引入以下方程：

式(4)～(7)中：τkk為亞格子應(yīng)力的各向同性部分，在不可壓縮流中通常可以忽略，根據(jù)Kolmogorov’s-5/3定律，通常Cs取0.10～0.30之間；Δ≈0.86 mm為濾波寬度，用來區(qū)分大尺度和小尺度渦旋的特征長度尺度；μt為SGS黏性項(xiàng)；為平均應(yīng)變率張量。

2.3 邊界條件

LES是一種瞬態(tài)湍流模型，將時(shí)間項(xiàng)中引入其邊界條件。

(1)初始表壓值為0 Pa。

(2)所有情況下的入口速度都設(shè)置為u∞=5 m/s，雷諾數(shù)為Re=u∞l/ν=13 611，其中l(wèi)是長方體風(fēng)洞的水力直徑。工況1～3中，在速度入口處不增加擾動(dòng)；工況4中，在入口處設(shè)置初始擾動(dòng)。

(3)在風(fēng)洞出口邊界處使用出流邊界條件。

(4)風(fēng)洞壁面和格柵采用無滑移條件和恒溫壁面T=283 K。圓柱也被設(shè)置為無滑移壁面，溫度為323 K。

在本研究中，LES模擬使用有限體積求解器進(jìn)行，該求解器采用有限體積法(FVM)在瞬態(tài)假設(shè)下求解不可壓縮Navier-Stokes方程。壓力項(xiàng)采用二階離散格式，動(dòng)量方程采用有界中心差分離散格式。數(shù)值計(jì)算方法采用PISO算法，亞格子模型采用Smagorinsky-Lilly模型，且Cs=0.10[17]。

3 結(jié)果與討論

3.1 仿真驗(yàn)證

為了驗(yàn)證本文模擬結(jié)果的正確性，用相同的模擬方案對參考文獻(xiàn)[4]中研究的RG60流場進(jìn)行了模擬。在圖2中，比較了數(shù)值模擬和風(fēng)洞試驗(yàn)中沿風(fēng)洞中心線的流向平均速度和湍流強(qiáng)度，結(jié)果吻合良好，計(jì)算誤差小于25%。

圖2 沿風(fēng)洞中心線的平均速度和湍流強(qiáng)度

3.2 流動(dòng)換熱特性

當(dāng)流體繞圓柱流動(dòng)時(shí)，由于流體與壁面之間的摩擦阻力，圓柱周圍形成邊界層，邊界層區(qū)域內(nèi)的流動(dòng)為層流。各層之間沒有質(zhì)點(diǎn)摻混，因此主要通過熱傳導(dǎo)進(jìn)行傳熱。邊界層分離后，由于宏觀上的相互混合，熱量被渦旋流動(dòng)帶走，其中的熱交換主要是對流換熱。

3.2.1 速度分布分析

文章通過對瞬時(shí)流向速度云圖和數(shù)據(jù)的分析，得到了湍流場的特征參數(shù)和傳熱效果。由于計(jì)算域中的流動(dòng)在流動(dòng)開始時(shí)尚未完全展開，僅選擇0.35～0.55 s時(shí)間段的流場數(shù)據(jù)進(jìn)行分析。所有工況都使用相同的笛卡爾坐標(biāo)系，且原點(diǎn)相同，x方向?yàn)榱飨蚍较?，y方向?yàn)檎瓜蚍较?，z方向?yàn)榇瓜蚍较颉?/p>

定性分析方面，提供了x-y平面(y=0)上的等值面云圖以及y-z平面在x/D=2，x/D=10，x/D=12.5,的等值面云圖。而在定量分析上，展示了沿著風(fēng)洞中心線上和圓柱體周圍的各變量變化趨勢。

圖3是t=0.35 s時(shí)的速度云圖，U是瞬時(shí)流向速度。圖中顯示，當(dāng)流體通過格柵時(shí)，格柵后部會(huì)產(chǎn)生一個(gè)局部低速區(qū)，渦旋逐漸在低速區(qū)產(chǎn)生，并沿流向方向交替脫落。對比3種工況，工況2產(chǎn)生的局部低速區(qū)域明顯大于其他兩種情況。從x/D=2平面上看，在柵格下游附近流動(dòng)速度云圖具有非常明顯的幾何印記，其中工況3最為明顯。在工況3和工況4中，流體繞流圓柱時(shí)，圓柱后部會(huì)出現(xiàn)較為明顯的卡門渦街。原因是：在這兩種情況下，圓柱體周圍環(huán)境的流動(dòng)湍流強(qiáng)度較小，一旦受到擾動(dòng)就會(huì)發(fā)生強(qiáng)烈地速度脈動(dòng)。

圖3 在0.5s時(shí)，x-z平面(y=0)和不同橫流平面(x/D=2,10,12.5)速度云圖

圖4 顯示了沿流場中心線的流向平均速度分布。對工況1～3，由于柵格阻塞的原因，0≤x/D≤2范圍內(nèi)流向平均速度明顯增加。對工況2，平均速度最高可達(dá)8.3 m/s，并且所有工況中速度衰減最小。此外，對工況3和工況4，在圓柱后部產(chǎn)生渦旋，并沿流向交替脫落。圓柱后部附近的流向平均速度急劇增加，并伴隨著震蕩，直到速度恢復(fù)到與進(jìn)口流速相當(dāng)。

圖4 沿中心線的平均速度剖面

3.2.2 湍流強(qiáng)度

圖5顯示了四種情況下沿中心線的湍流強(qiáng)度(tu)分布。湍流強(qiáng)度由公式(8)計(jì)算：

式中：u′為瞬時(shí)速度的均方根為時(shí)間平均流向速度。

從公式(8)中可以看出，湍流強(qiáng)度與隨機(jī)湍流波動(dòng)直接相關(guān)。較高的湍流強(qiáng)度意味著更大的隨機(jī)湍流波動(dòng)。

如圖5所示，工況1～3中圓柱前面區(qū)域湍流強(qiáng)度均有不同程度的遞增趨勢。在第三種工況下，流場產(chǎn)生的最高的湍流強(qiáng)度可達(dá)20%。但是，在沒有柵格的情況下湍流強(qiáng)度是恒定的。在圓柱后部區(qū)域，由于流場受到圓柱的擾動(dòng)，湍流效果會(huì)進(jìn)一步加強(qiáng)。其中，工況4產(chǎn)生的湍流強(qiáng)度最大約為56%，在工況1中約為48%。而且，在工況3和工況4的圓柱后部附近，湍流強(qiáng)度呈波浪式減小趨勢。

圖5 沿中心線的湍流強(qiáng)度

3.2.3 壓降分析

為了表示流場中圓柱周圍的壓力分布，文章引入了無量綱壓降公式，定義如下：

式中：P為圓柱周圍的總壓力；Pinlet為入口的總壓力；ρ、u分別為入口處的平均速度和密度。

圖6顯示了氣流在圓柱體周圍流動(dòng)時(shí)的邊界層分離現(xiàn)象。圖7顯示了圓柱體周圍的速度分布。在所有工況下，平均速度沿拋物線分布，先增大后減小，最后趨于穩(wěn)定。當(dāng)內(nèi)壓梯度小于零時(shí)，邊界層內(nèi)流動(dòng)加速，流動(dòng)動(dòng)能增大。當(dāng)內(nèi)壓梯度大于零時(shí)，邊界層中的流動(dòng)減速，動(dòng)能迅速耗盡。當(dāng)邊界層中的剪切應(yīng)力大于上游流體中的剪切應(yīng)力，就會(huì)出現(xiàn)與初始流動(dòng)方向相反的回流，從而形成了邊界層分離。對工況1和工況2，邊界層分離發(fā)生在φ≈100°；對工況3和工況4，邊界層分離發(fā)生在φ≈160°。

圖6 邊界層形成示意圖

圖7 圓柱周圍的速度

在圖 8中[18]，展示了圓柱周圍 0°≤φ≤180°的壓降系數(shù)分布。其中工況1和工況2在圓柱周圍呈線性下降趨勢，對于情況1，壓降系數(shù)的變化非常小。對情況3和情況4，壓降系數(shù)的趨勢幾乎相同。壓降系數(shù)在φ≈0°時(shí)最大，在 0°≤φ≤180°時(shí)呈指數(shù)下降。φ≈80°壓降系數(shù)降至最低點(diǎn)，并開始出現(xiàn)回升。在80°≤φ≤100°階段壓降系數(shù)略微上升，然后變得穩(wěn)定。同時(shí)，邊界層動(dòng)能耗盡開始出現(xiàn)分離。最終，工況3的最小壓降系數(shù)維持在-4.5左右，低于工況4的-3.5。

圖8 圓柱周圍的壓降系數(shù)分布

3.2.4 渦結(jié)構(gòu)

圖9顯示了四種情況下的渦結(jié)構(gòu)。在工況1和工況2中，流場中柵格附近和圓柱后部產(chǎn)生了少許的渦結(jié)構(gòu)，而工況3無論在柵格附近還是在圓柱周圍都產(chǎn)生了較多的渦旋結(jié)構(gòu)，這是因?yàn)楣r3的柵格結(jié)構(gòu)前后壓強(qiáng)差較大，更容易在流場中誘導(dǎo)渦旋產(chǎn)生。

圖9 介于-51 000和-50 000之間渦結(jié)構(gòu)

3.3 換熱特性

圖10 顯示了加熱圓柱表面的Nu分布的俯視圖、前視圖(面向入流方向)、側(cè)視圖和后視圖。圖中可以看出前三種工況在邊界層處均有對流換熱和熱傳導(dǎo)現(xiàn)象，而且邊界層中的換熱效果要強(qiáng)于尾流區(qū)域中渦流脫落產(chǎn)生的熱對流。這是由于柵格結(jié)構(gòu)在上游流場產(chǎn)生的一定數(shù)量的渦結(jié)構(gòu)導(dǎo)致。θ是邊界層分離角，即x方向與邊界層分離點(diǎn)和圓心之間的線的角度。直觀上看，工況1和工況2分離點(diǎn)均發(fā)生在，遠(yuǎn)大于工況2和工況3的θ≈85°。此外，工況1～3中Nu的分布在邊界層區(qū)域是非常不均勻的，特別是在案例3中，存在相當(dāng)一部分Nu≥6 000的區(qū)域，這些良好的換熱效果是由于流場中大量的渦結(jié)構(gòu)存在。

圖10 圓柱體的表面努塞爾數(shù)

圖11 是平面的溫度云圖。直觀上看，圖11中得出的結(jié)論與圖8、圖9和圖10有很好的匹配。圖12顯示了圓柱表面0°≤φ≤180°努塞爾數(shù)的分布情況。在 0°≤φ≤78°，工況 3的最大Nu約為 6 200，比工況4高約30%，比工況1和工況2高約63%。然而，它在105°≤φ≤180°之間急劇下降，達(dá)到最低值1 000后逐漸恢復(fù)，最后達(dá)到約2 600。這表明，工況3的熱交換效果是不均勻的。從面平均努塞爾數(shù)上來看，工況3的換熱效果是最好的，平均能達(dá)到3 870。

圖11 t=0.5 s時(shí)平面z=0上的溫度等值線圖

圖12 圓柱體周圍的Nu

工況1和工況2中兩個(gè)分布的趨勢和數(shù)值大致相同，沒有急劇上升和下降。最大值和最小值之間的差距約為2 300，比工況3的5 000和工況4的4 000小得多。這表明，在工況1和工況2中，流場可以產(chǎn)生更均勻的傳熱效果，面平均努塞爾數(shù)分別為2 493、2 666。

表1 面平均努塞爾數(shù)

4 結(jié)語

基于LES模擬了四個(gè)流場(①三角形柵格；②倒三角形柵格；③方形柵格；④無柵格)中圓柱體周圍的流動(dòng)演變和傳熱。分析了這四個(gè)流場的流動(dòng)特性和傳熱特性，并得出以下結(jié)論：

(1)文章比較了三角形柵格、倒三角形柵格、方形柵格和無柵格擾動(dòng)下的流場。結(jié)果發(fā)現(xiàn)，倒三角柵格產(chǎn)生的局部低速區(qū)域明顯大于其他兩種情況，因此產(chǎn)生的湍流強(qiáng)度較高。

(2)工況1和工況2中，圓柱周圍的壓降系數(shù)呈現(xiàn)線性下降趨勢，在工況1中，壓降系數(shù)的變化非常小，幾乎沒有變化。

(3)在工況1和工況2中，流場能產(chǎn)生比較均勻的傳熱效果，但整體傳熱效果不如工況4。

(4)在工況3中，邊界層提前分離，更多的渦旋被激發(fā)，增強(qiáng)了邊界層區(qū)域的傳熱效果。此外，從面平均努塞爾數(shù)高達(dá)3 870來看，傳熱效果最好。