王亞龍，朱瀟瀟

（航天工程大學(xué)，北京 101416）

0 引言

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

1 算例的選用

1.1 研究狀況

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

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

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

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

1.2.1 幾何模型

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

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

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

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

1.2.2 數(shù)值方法

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

連續(xù)性方程：

動(dòng)量方程：

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

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

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

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

2 結(jié)果與討論

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

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

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

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

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

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

2.2 k-w模型和k-ε模型的對(duì)比

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

圖4 不同湍流模型下的流場(chǎng)

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

圖5 不同湍流模型下的渦量場(chǎng)

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

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

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

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

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

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

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

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

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

圖9 定常與非定常模k-w型渦量場(chǎng)

3 結(jié)論

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

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

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

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