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

王亞龍，朱瀟瀟

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

0 引言

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

1 算例的選用

1.1 研究狀況

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

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

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

1.2 幾何模型及數值方法

1.2.1 幾何模型

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

圖1 幾何模型及計算域網格

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

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

1.2.2 數值方法

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

連續性方程：

動量方程：

升力系數Cd定義為：

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

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

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

2 結果與討論

2.1 網格無關性驗證

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

圖3 三種網格尺寸下的渦脫落頻率

根據公式（5）計算出St數，見表1。

表1 不同網格下的St數與理論誤差

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

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

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

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

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

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

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

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

2.3 與k-ε、k-w湍流模型對雷諾數的響應

為研究雷諾數對湍流模型的影響，分別計算了雷諾數為800和8000時的圓柱繞流案例。圖7和圖8分別是不同雷諾數下k-w模型與k-ε模型的渦量圖?？梢钥闯鲭S著雷諾數的增加渦量增加，但k-w模型在雷諾數為800時依然能夠清晰地捕捉到渦脫落，而k-ε模型在低雷諾數時渦量趨于穩定，沒有明顯的渦脫落行為，整個流場變成穩態。說明k-ε模型對低雷諾數流動不敏感，沒能成功捕捉到低雷諾數下的渦脫落行為。

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

圖7 k-w模型對雷諾數響應的渦量圖

圖8 k-ε模型對雷諾數響應的渦量圖

2.4 定常與非定常k-w模型的區別

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

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

3 結論

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

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

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

渦脫落為瞬態過程，對比穩態和瞬態的k-w模型，發現在穩態計算時流場趨于穩定，定常k-w模型沒能捕捉到渦脫落的過程。