劉　璐

(西安航空學(xué)院 飛行器學(xué)院，陜西 西安 710077)

?

基于空氣動力學(xué)N-S方程的正方體繞流問題研究

劉璐

(西安航空學(xué)院 飛行器學(xué)院，陜西 西安 710077)

繞流是一個經(jīng)典的空氣動力學(xué)問題，對于渦流的產(chǎn)生和形成一直是教學(xué)過程中的難點。通過對流體繞過正方體的研究，可以具體地展示出渦流隨時間產(chǎn)生和發(fā)展，幫助理解渦流的形成過程和規(guī)律。首先對動量方程進行分析推導(dǎo)，得到壓力泊松方程，使用交錯網(wǎng)格投影法將N-S方程離散化，并得到離散化的壓力泊松方程，通過使用SOR技術(shù)來求解。利用MATLAB建模，根據(jù)所得到的圖像來具體分析流體繞過正方體的特征。結(jié)果顯示流體從正方體的一角分離、卷起并在正方體后方立即形成一對反向旋渦，該反向漩渦隨即分解，后形成幾個分散的漩渦，形象地展示了渦流隨時間變化的過程。

N-S方程；正方體繞流；數(shù)值仿真

正方體是一種典型的鈍頭柱體，關(guān)于正方體繞流問題的研究一直是空氣動力學(xué)領(lǐng)域的熱點問題之一，比如在大型建筑物、冷卻塔、熱交換器等等方面[1]。這是因為空氣繞過固體表面物體產(chǎn)生和脫落的渦流會影響物體上被施加的載荷，如果分離頻率接近于物體的自激震蕩頻率，可能會使物體損壞[2]。

繞流現(xiàn)象是一類現(xiàn)實生活廣泛存在的空氣動力學(xué)問題。由于流體流過結(jié)構(gòu)時會產(chǎn)生分離現(xiàn)象，結(jié)構(gòu)會受到持續(xù)力的作用，因而產(chǎn)生周期性震動，可能會使結(jié)構(gòu)產(chǎn)生疲勞等現(xiàn)象。近年來，關(guān)于繞流現(xiàn)象的研究一直不斷。得益于計算機速度的提高和計算方法的發(fā)展，數(shù)值模擬成為研究繞流問題的重要方法。王亞玲等利用不可壓縮N-S方程，模擬了不可壓縮流體的圓柱體繞流問題[3]。馬金花等利用有限元方法求解雷諾平均N-S方程，對圓柱體繞流問題進行數(shù)值模擬[4]。但他們采用N-S方程的求解方法是流函數(shù)-渦量法，引入流函數(shù)，連續(xù)方程自動滿足；然后讓u動量方程和v動量方程分別對y和x求偏導(dǎo)，消掉壓力項，并利用渦量-速度的關(guān)系式合并相關(guān)速度項，就能得到渦量方程；同樣，再利用流函數(shù)-速度-渦量關(guān)系可得流函數(shù)的Poisson方程，這樣，原方程就簡化為了求流函數(shù)和渦量兩個變量的方程。這種方法相對于SOR法較簡單，不需要假設(shè)-迭代，但是邊界條件的處理尤其關(guān)鍵，當(dāng)邊界條件簡單的情況下，這種方法有效。但當(dāng)繞流固體為正方體時，流函數(shù)-渦量法計算較慢，還可能會出現(xiàn)不收斂的情況。

本文應(yīng)用MATLAB代碼定義正方體與流體接觸面，然后將邊界條件應(yīng)用于N-S方程，再使用交錯網(wǎng)格投影法，經(jīng)過迭代求解N-S方程，從而得到不同時間下的渦線圖。

1　N-S方程及求解方法

N-S方程，即納維-斯托克斯方程。本文中N-S方程是針對二維非定常不可壓縮流體繞過正方體的研究，方程具體形式如下所示：

(1)

(2)

(3)

其中，u和v是速度在兩個方向上的分量，p是壓力，而Re是雷諾數(shù)[5]。

本文采用的N-S方程求解方法為壓力泊松方程法。這種方法的核心還是迭代，先將u動量方程與v動量方程分別對x和y求偏導(dǎo)，然后相加，并利用連續(xù)方程合并相關(guān)項，這樣就得到了壓力泊松方程。先不需要考慮壓力，而是從動量方程得到一個中間速度，然后把此速度代入泊松方程，用泊松方程代替連續(xù)方程，這樣可以得到一個關(guān)于兩個方程的迭代式[6]。

2　泊松方程的求解方法

在本文中，使用逐次超松弛法來解泊松方程。逐次超松弛法(successive over relaxation method)，簡稱SOR方法[7]，是解線性方程組的常用迭代法之一，它由高斯-賽德爾迭代法經(jīng)線性加速處理而得到[8]。原理如下：

假設(shè)方程組Ax=b的系數(shù)矩陣AX=b

(4)

A是n×n的矩陣 元素從a11到ann

(5)

使用下面的迭代格式來建立解方程組Ax=b的逐次超松弛迭代法(SOR)。

(6)

這里引入ω，ω為松弛因子，其定義如下：

(7)

nx和ny是網(wǎng)格上x方向和y方向上格子的數(shù)目，這會比用傳統(tǒng)的高斯-賽德爾迭代法求解線性方程時擁有更快的趨同速率[9]。

3　計算步驟

使用MATLAM軟件進行求解：

(1)定義任意幾何形狀的圖形：創(chuàng)建主體形態(tài)，在定義域內(nèi)定義固體表面X、Y方向的數(shù)值，正方體的邊長為2。

(2)參數(shù)的定義和初始化變量：確定雷諾數(shù)和時間步長，選取雷諾數(shù)為100，時間步長為0.1。

(3)找到主體的指數(shù)：將從可變形狀中指定主體的位置轉(zhuǎn)換到u,v和p的儲存位置。

(4)檢查穩(wěn)定性標(biāo)準(zhǔn)：之前定義的時間步長將通過穩(wěn)定性的標(biāo)準(zhǔn)來驗證。

(5)設(shè)置域邊界的速度邊界條件。

(6)設(shè)置主體的速度邊界條件：首先將主體內(nèi)全部速度點設(shè)置為零，然后從先前定義過的參數(shù)信息中尋找所需要的邊界條件的位置。

(7)求解表面上的u和v：使用循環(huán)計算表面上u和v的臨時值。

(8)在n+1的時間階上求解壓力：用泊松算法求解離散泊松方程。

(9)求解u和v：在n+1的時間階上的速度分量是使用循環(huán)來計算出的。循環(huán)包含了除主體內(nèi)的所有速度點。

(10)設(shè)置壓力邊界條件：對計算域和主體邊界施加針對壓力的離散形式的諾依曼條件。

(11)更新t和呈現(xiàn)結(jié)果：時間自動跳過一個時間步長直到時間已經(jīng)超過開始時設(shè)置的最大時間，然后繪制出渦線圖。

4　計算結(jié)果及分析

設(shè)置正方體的邊長為2，氣流的雷諾數(shù)Re=100，在固體表面u=v=0，時間步長為0.1。然后采用MATLAB軟件計算t分別在10，50，100s時的正方體繞流瞬時渦線，如圖1～圖3所示。

圖1t=10s的瞬時渦線

圖2　t=50s的瞬時渦線

圖3　t=100s的瞬時渦線

圖中顯示了在t=10s，50s，100s時單位正方體繞流的瞬時渦線，這個正方體被放置在上下對稱的區(qū)域范圍內(nèi)，距離入口L=0.5。在入口處，指定了一個簡單均勻的速度，隨后流量在到達主體之前得到充分發(fā)展。當(dāng)雷諾數(shù)為100時，數(shù)值解可以由不均勻的細小網(wǎng)格求得。結(jié)果表明，下游流量是不穩(wěn)定的；流量從主體脫離后，在主體后形成旋渦。在第一個10s，交叉流從主體的一角分離、卷起并在主體后立即形成一對反向旋渦，如圖1所示。然后反向漩渦開始分解，隨后形成幾個分散的漩渦，如圖2和圖3所示。在t=100s后，形成非對稱穩(wěn)定渦旋。

氣流在正方體前部發(fā)生分離，一部分是因為物體的阻礙作用，另一部分是因為氣流與固體表面有粘性力作用，消耗了流體的動能，使得氣流不能沿固體表面流動，從而發(fā)生了氣流分離。而在正方體的后方形成了一定的低壓區(qū)內(nèi)，不穩(wěn)定氣流形成了渦旋，隨著時間的增長，渦旋也不斷地發(fā)展、分解，最后形成分散的非對稱渦旋。

5　結(jié)語

本文應(yīng)用投影網(wǎng)格交錯法，求解了在二維定常流動下，雷諾數(shù)為100時，流體繞過正方體后，在t=10s，50s，100s時，渦流形成的情況，與實驗結(jié)果較為吻合。今后可以將MATLAB應(yīng)用于求解流體繞過任意形狀物體的過程，以及任意氣流條件下的繞流，可以更加深入了解渦流形成的過程，以減少渦流的形成和影響，可以發(fā)現(xiàn)更加方便計算和收斂N-S方程的數(shù)值方法。

[1] 劉宇，蘇中地.不同雷諾數(shù)下方柱繞流的數(shù)值模擬[J].中國計量學(xué)院學(xué)報，2006，17(1)：40-49

[2] 謝志剛，徐春曉，崔桂香，等.方柱繞流大渦模擬[J].計算物理，2007,24(2)：171-180.

[3] 王亞玲，劉應(yīng)中，繆國平.圓柱繞流的三維數(shù)值模擬[J].上海交通大學(xué)學(xué)報，2001，35(10)：1464-1469.

[4] 馬金花，金生，賀德馨，等.圓柱體繞流的數(shù)值模擬[J].山東建筑工程學(xué)院學(xué)報，2001，16(2):45-47.

[5] 王建，金志浩，張玉嬌.方柱繞流的數(shù)值模擬[J].沈陽化工學(xué)院學(xué)報，2007，21(2)：117-120.

[6] 景思睿，張鳴遠.流體力學(xué)[M].西安：西安交通大學(xué)出版社，2001:210-211.

[7] Woodward P,Colella P.The Numerical Simulation of Two-dimensional Fluid Flow with Strong Shocks[J].Journal of Computational Physics,1984,54(1):115-173.

[8] Xu X P,Needleman A.Numerical Simulations of Fast Crack Growth in Brittle Solids[J].Journal of the Mechanics & Physics of Solids,1994,42(42):1397-1434.

[9] Kopera M A,Cantwell C,Kerr R M,et al.Direct Numerical Simulation of Turbulent Flow over A Backward-facing Step[C]//Meeting of the APS Division of Fluid Dynamics.2011,330(1):349-374.

[責(zé)任編輯、校對：李琳]Analysis of Fluid Flow around a Cube Based on Aerodynamic N-S Equations

LIULu

(School of Aircraft,Xi'an Aeronautical University,Xi'an 710077,China)

Flow around is a typical aerodynamic issue,and the generation and formation of vortexes is usually understood with difficulty.This paper aims to study the fluid flow around a cube to demonstrate the generation and evolution of a vortex.Based on the Poisson pressure equation,a discretized Poisson equation was proposed by discrete N-S equation using the collocated mesh projection method,and the results were obtained through SOR.The modeling of this flow was conducted with MATLAB numerical calculation and image processing.The characteristics of fluid flow around a square were analyzed according to the images generated.The results show that the fluid was separated from the tip of the cube,curled up afterward,and generated a reverse vortex,which was immediately decomposed into several discrete vortexes thereafter.Analyzing through such numerical model illustrates the detailed changing process of vortex.

Navier-Stokes equations;flow around cube;numerical simulation

2016-09-01

劉璐(1992-)，女，陜西延安人，博士，主要從事空氣動力學(xué)以及計算流體力學(xué)研究。

O357.5

A

1008-9233(2016)05-0026-03