基于圓柱繞流分析k-ε和k-w湍流模型的差異

王亞龍，朱瀟瀟

（航天工程大學，北京 101416）

0 引言

隨著計算機技術(shù)的普及和計算方法的發(fā)展。作為計算流體力學和計算傳熱學的新型研究方法，CFD技術(shù)得到了迅猛發(fā)展。因為相對于實驗研究，數(shù)值模擬有很多獨特的優(yōu)點。例如：成本低、周期短以及能獲得較為完整的數(shù)據(jù)。同時能夠模擬出實際運行過程中的各種測試數(shù)據(jù)。對于新產(chǎn)品的研發(fā)設(shè)計和改造有著重要的指導作用。所以，目前CFD技術(shù)不僅在電子，制冷等實際工程領(lǐng)域中的應(yīng)用越來越廣泛。在現(xiàn)代航空航天、核能工程等領(lǐng)域里的應(yīng)用也越發(fā)的廣泛與深入。首先，商業(yè)CFD軟件數(shù)值模擬功能是強大的，它是目前功能最全面適用性最廣使用最廣泛的流體軟件之一。是基于CFD軟件群思想，從用戶需求角度出發(fā)，針對各種復雜流動的物理現(xiàn)象所采用的數(shù)值解法。從而能夠在計算速度穩(wěn)定性和監(jiān)督等方面達到優(yōu)化組合[1]，高效的解決各個領(lǐng)域的復雜流動計算問題；模擬流動傳熱和化學反應(yīng)等諸多物理現(xiàn)象。它能夠提供非常靈活的網(wǎng)絡(luò)特性。讓使用者可以使用非結(jié)構(gòu)網(wǎng)格，包括三角形、四邊形、六面體、金字塔型網(wǎng)格以及混合型非結(jié)構(gòu)網(wǎng)格[2]，并允許用戶根據(jù)求解的具體情況對網(wǎng)格進行修改。包括對網(wǎng)格的細化和粗化，從而解決具有復雜幾何形狀的流動換熱問題。CFD已經(jīng)成為流體數(shù)值模擬過程的一種重要的手段，然而，湍流是著名的難題，湍流是相對于層流來說的，流體的層流被看作是比較平滑的流動，而湍流無論是在空間還是在時間域內(nèi)都是無規(guī)則的運動。這一定義是1883年由Reynolds實驗時發(fā)現(xiàn)并一直沿用至今，實驗研究在湍流研究中占有十分重要的地位[3]，從湍流的發(fā)現(xiàn)，層流到湍流的過渡，湍流擬序結(jié)構(gòu)的發(fā)現(xiàn)和研究都與實驗密切相關(guān)，同時湍流理論研究進展也積極推動著湍流實驗研究的深入。如20～30年代的各種湍流理論，如Prandtl動量輸運和Taylor的渦量輸運理論，以及后來的Karman相似理論為湍能的產(chǎn)生和耗散之間的平衡關(guān)系的實驗研究起了指導性作用。近年來對擬序結(jié)構(gòu)的研究進一步揭示了湍流內(nèi)在的一些重要機理，如湍流的擴散和發(fā)展不僅僅是小尺度渦旋隨機擴散的結(jié)果，更主要是由大尺度擬序結(jié)構(gòu)的相互作用和卷并所致[4]。由于流現(xiàn)象廣泛存在于自然界和工程技術(shù)的各個領(lǐng)域，因此湍流基礎(chǔ)理論研究取得的進展就可能為經(jīng)濟建設(shè)和國防建設(shè)的廣泛領(lǐng)域帶來難以估量的效益。盡管湍流研究相當困難，但是仍然有大量的國內(nèi)外學者致力于這一領(lǐng)域的工作。目前隨著計算機技術(shù)和測量技術(shù)的不斷發(fā)展，湍流的精細實驗正在進一步地展開，它對深入認識湍流的物理本質(zhì)至關(guān)重要[5]；相關(guān)學科的發(fā)展也推動了湍流的研究，經(jīng)過很多科研人員的研究，結(jié)果表明：選擇合適的湍流模型對于研究結(jié)果的合理性和正確性極為重要，本文基于CFD軟件分析了二維圓管流道內(nèi)的圓柱繞流問題。通過選取k-w和k-ε湍流模型，分析了兩種湍流模型在計算圓柱繞流時的差異。

1 算例的選用

1.1 研究狀況

一個世紀以來，圓柱繞流一直是眾多理論分析、實驗研究及數(shù)值模擬對象。但迄今對該流動現(xiàn)象物理本質(zhì)的理解仍是不完整的。圓柱繞流中，起決定作用的是雷諾數(shù)，但還受到許多因素，如阻塞比，來流湍流度，下游邊界條件等的影響[6]。隨著雷諾數(shù)的增加，粘性不可壓縮流體繞圓柱的流動會呈現(xiàn)各種不同的流動狀態(tài)，在小雷諾數(shù)時，流動是定常的，隨著雷諾數(shù)的增加，圓柱后會出現(xiàn)一對尾渦。當雷諾數(shù)較大時，尾流首先失穩(wěn)，出現(xiàn)周期性的振蕩。而后附著渦交替脫落，瀉入尾流形成Karman渦街，隨著雷諾數(shù)的增加，流動變得越來越復雜，最后發(fā)展為湍流[7]。

一般認為圓柱繞流有2種定常的流動圖案：雷諾數(shù)為較小時，圓柱后無尾渦；當雷諾數(shù)為較大時，圓柱后有一對對稱的尾渦。關(guān)于定常流失穩(wěn)以及出現(xiàn)湍流的臨界雷諾數(shù)主要是通過應(yīng)用流場顯示技術(shù)觀察流動形態(tài)得到的，所以不是準確值[8]。

圓柱周圍的強制對流傳熱在工程中有許多應(yīng)用，例如燃氣輪機冷卻，熱交換器，核燃料棒等。圓柱或其他形狀的圓柱的傳熱增強和錯流引起了研究人員的充分關(guān)注。Buyruk等對雷諾數(shù)為120和390的錯流圓柱流的層流和傳熱特性進行了數(shù)值和實驗研究[9]。對于390的雷諾數(shù)，與實驗相比，他們通過數(shù)值計算獲得了更好的預測結(jié)果。力由Hover等人測量。在跨度為60cm的剛性圓柱體的兩端，在進入的水流中在Re處執(zhí)行橫向振蕩。Catalano等研究了大渦模擬（LES）的可行性和精度，并考慮了在超臨界狀態(tài)下考慮圓柱體繞流的高雷諾數(shù)復雜湍流。但是，沒有捕獲雷諾數(shù)依賴性，并且隨著雷諾數(shù)的增加，解決方案的準確性降低[10]。Bouhairie和Chup使用二維模型研究了從圓柱表面到橫流的傳熱，其雷諾數(shù)范圍從Re=200到15550。結(jié)果表明，二維模型捕獲了不穩(wěn)定過程并產(chǎn)生了與可用實驗數(shù)據(jù)一致的結(jié)果[11]。它提供了相對正確的總體結(jié)果，前停滯和總傳熱速率。Patnana等研究在二維（2-D）非定常流動狀態(tài)下，將a浸沒在流動的冪律流體中的圓柱體的強制對流傳熱特性。根據(jù)他們的研究，無論流動行為指數(shù)如何，努塞爾數(shù)都隨著雷諾數(shù)或普朗特數(shù)的增加而增加[12]。Khan等使用Karman-Pohlhausen方法研究了無限圓柱周圍的流體流動和熱量傳遞。阻力系數(shù)和傳熱系數(shù)的結(jié)果與圓柱的實驗/數(shù)值數(shù)據(jù)非常吻合[13]。

1.2 幾何模型及數(shù)值方法

1.2.1 幾何模型

如圖1所示，計算域為二維等直圓管，中間為圓柱繞流體，圓柱直徑D為20mm。為有效捕捉渦脫落過程，在圓柱周圍和壁面都進行了網(wǎng)格加密，網(wǎng)格總數(shù)為9萬，最小網(wǎng)格尺寸為0.4mm，壁面邊界層y+＜1。管內(nèi)流體為自定義氣體，由于流速較低考慮氣體為不可壓流體，密度ρ為1kg/m3，動力粘度為μ=0.00025Pa·s。為保證流動為充分發(fā)展的湍流，如圖2所示，在入口根據(jù)式（1）利用UDF定義充分發(fā)展的湍流。當雷諾數(shù)Re為800時，圓管入口流速為1m/s；雷諾數(shù)Re為8000時，管的入口速度為10m/s。

圖1 幾何模型及計算域網(wǎng)格

式中：Umax—氣體最大流速；Dinlet—圓管入口直徑；y—縱向坐標。邊界條件如圖2所示，氣體為不可壓粘性流體，因此入口為速度入口，出口為壓力出口，流道內(nèi)圓柱和圓管壁面均為無滑移絕熱壁面。

圖2 入口速度分布及邊界條件

1.2.2 數(shù)值方法

本文中，不可壓粘性流體的控制方程包含連續(xù)性方程和動量方程

連續(xù)性方程：

動量方程：

升力系數(shù)Cd定義為：

斯特勞哈爾數(shù)St是振蕩流的無量綱度量，其定義如下：

式中：f—渦脫落的頻率；L—特征長度，本文中為流場中圓柱直徑，U是來流速度。理論指出在800＜Re＜20000時St數(shù)幾乎不變?yōu)?.2。通過計算St數(shù)并與理論值進行比對，比較不同湍流模型對渦脫的捕捉能力。

本文中采用PISO算法求解壓力速度耦合，壓力項采用PRESTO！離散。方程中的對流項采用二階迎風差分格式來離散。為提高計算精度動量方程的離散采用了QUICK格式；瞬態(tài)方程采用了二階隱式離散以提高渦脫落的捕捉能力[14]。

2 結(jié)果與討論

2.1 網(wǎng)格無關(guān)性驗證

基于k-w模型，在雷諾數(shù)為8000時進行了網(wǎng)格無關(guān)性驗證。通過監(jiān)測圓柱表面升力系數(shù)的變化來表示渦脫落過程。對升力系數(shù)進行快速傅里葉變換從而得到渦脫落頻率。在三種尺寸的網(wǎng)格上進行了計算。網(wǎng)格數(shù)量分別為2萬、9萬和30萬，見圖3。

圖3 三種網(wǎng)格尺寸下的渦脫落頻率

根據(jù)公式（5）計算出St數(shù)，見表1。

表1 不同網(wǎng)格下的St數(shù)與理論誤差

通過比較可以看出，網(wǎng)格在9萬和30萬時誤差較小，并且網(wǎng)格從9萬增加到30時頻率變化不足1%。考慮計算資源的限制，選取數(shù)量為9萬的網(wǎng)格進行計算。

2.2 k-w模型和k-ε模型的對比

當雷諾數(shù)為8000時，比較 了k-w和k-ε模型 在計算圓柱繞流時的差異。如圖4 所示為不同湍流模型下的流場，可以看出在使用k-ε模型計算時，圓柱后方存在較大區(qū)域的低速回流區(qū)，在向下游流動的過程中流場震蕩迅速衰減，流線較為平穩(wěn)。而k-w模型在捕捉雷諾數(shù)較高時捕捉到的渦的尺寸較小，但能夠捕捉到逐漸向下游衰減脫落的過程，流場震蕩幅值較大。

圖4 不同湍流模型下的流場

圖5 顯示了不同湍流模型下的渦量場，可以看出兩種模型在渦量的捕捉上存在巨大差異，k-w模型明顯更能精確地捕捉到渦量，而k-ε模型存在著嚴重的失真。k-ε模型計算的渦量嚴重偏小，說明此模型流場外圍渦的捕捉能力不強。

圖5 不同湍流模型下的渦量場

圖6 顯示了渦脫落的幅頻圖，從中可以看出用k-ε模型計算的渦脫落頻率為81Hz，而k-w模型計算的為99Hz，對應(yīng)的St數(shù)分別為0.162和0.198。與理論值0.2相比k-ε模型誤差較大，說明k-ε模型不適應(yīng)于圓柱繞流的仿真。

圖6 不同湍流模型下的頻譜圖

2.3 與k-ε、k-w湍流模型對雷諾數(shù)的響應(yīng)

為研究雷諾數(shù)對湍流模型的影響，分別計算了雷諾數(shù)為800和8000時的圓柱繞流案例。圖7和圖8分別是不同雷諾數(shù)下k-w模型與k-ε模型的渦量圖。可以看出隨著雷諾數(shù)的增加渦量增加，但k-w模型在雷諾數(shù)為800時依然能夠清晰地捕捉到渦脫落，而k-ε模型在低雷諾數(shù)時渦量趨于穩(wěn)定，沒有明顯的渦脫落行為，整個流場變成穩(wěn)態(tài)。說明k-ε模型對低雷諾數(shù)流動不敏感，沒能成功捕捉到低雷諾數(shù)下的渦脫落行為。

比較圖7（a）和圖8（b）可以看出，在雷諾數(shù)為800時利用k-w模型捕捉的渦量場與雷諾數(shù)為8000時用k-ε捕捉的渦量場十分相似，渦量的強度也相差不大，這說明k-ε模型只有在雷諾數(shù)達到一定程度后才能捕捉到渦脫落，且渦量強度誤差較大。由此可以推測k-ε模型可能更適應(yīng)于高雷諾數(shù)流動。

圖7 k-w模型對雷諾數(shù)響應(yīng)的渦量圖

圖8 k-ε模型對雷諾數(shù)響應(yīng)的渦量圖

2.4 定常與非定常k-w模型的區(qū)別

為研究定常與非定常k-w模型在計算圓柱繞流時的差異，在雷諾數(shù)為800時，分別使用穩(wěn)態(tài)和瞬態(tài)的k-w模型進行了圓柱繞流的仿真。如圖9所示顯示了雷諾數(shù)為800時定常與非定常k-w模型的渦量場，可以看出即使流動為層流，流場中依然存在非定常的渦脫落行為。在用定常計算時流場迅速收斂達到穩(wěn)態(tài)，無法捕捉到渦脫落這種瞬態(tài)行為，而非定常k-w模型即使在雷諾數(shù)較低的情況下也能捕捉到渦脫落行為。

圖9 定常與非定常模k-w型渦量場

3 結(jié)論

本文基于CFD軟件，分別選擇在k-ε和k-w模型下對圓柱繞流問題進行了數(shù)值模擬。首先驗證了網(wǎng)格無關(guān)性，然后模擬了兩種不同模型下的渦脫落頻率和流場云圖。比較了兩種湍流模型對不同雷諾數(shù)的響應(yīng)。得出結(jié)論如下：

通過監(jiān)測圓柱表面的升力系數(shù)，得到了渦脫落頻率，發(fā)現(xiàn)Re=8000時k-w模型能夠有效捕捉渦脫落，對應(yīng)的St數(shù)為0.198與理論值吻合良好。采用k-ε模型計算圓柱繞流誤差較大。

對比了不同雷諾數(shù)下k-ε和k-w模型渦脫落的捕捉能力，發(fā)現(xiàn)在雷諾數(shù)低于一定值后k-ε沒能捕捉到渦脫落行為，說明k-ε模型不適合低雷諾數(shù)流動。

渦脫落為瞬態(tài)過程，對比穩(wěn)態(tài)和瞬態(tài)的k-w模型，發(fā)現(xiàn)在穩(wěn)態(tài)計算時流場趨于穩(wěn)定，定常k-w模型沒能捕捉到渦脫落的過程。