2.3　CFD理論

最早的非定常不可壓縮流體Navier-Stokes方程是以原始變量的形式來處理流動問題。一方面，原始變量公式提供了一種獨特的視角認識質量守恒定律與邊界條件。類似的算法可以處理二維和三維的流動問題。在必要情況下可以直接使用直觀的湍流模型。另一方面，不可壓縮流動的原始變量公式缺乏一個物理意義明確的壓力方程。通常采用Poisson的壓力方程［30］，對動量方程進行發散。也可以采用另外一種解決辦法，添加一個與速度場直接相關的壓力修正項［31］，這個壓力修正項最后會出現在連續性方程中，左邊乘以一個線性算子使壓力方程能夠順利求解。算子的選擇是根據不同的壓力校正方法來確定的。只有當矩陣對角線數值被保留在算子上時才能得到一個人為的壓縮性公式［32］。如果采用不包含矩陣的共軛梯度法來處理原始變量方程組，就無法順利求解壓力方程［33］。

原始變量公式在經過反復離散化后存在速度-壓力解耦合的問題。交錯網格法［34］雖然可以進行補救，但大大增加了程序的復雜性。如果對速度和壓力進行不同程度的插值（尤其是對有限元分析［35］ ）或者在連續性方程中引入高階壓力耗散項［36］，非交錯網格的壓力振蕩則會被抑制。

在本文中，非定常不可壓縮Navier-Stokes方程的原始變量形式用壓力修正方法進行了處理。其中，壓力修正采用的是不包含矩陣的求解方法。處理層流問題的相應方程組已由Chaviaropoulos［37］提出，它適用于內部、外部二維流場的計算。曲線坐標系上，控制方程以非定常的形式在結構化網格中進行離散。速度分量在網格節點上計算，而壓力在網格中心計算。動量方程用一個基于MSIP算法［38］的不完全LU因式分解近似方法（An Incomplete LU Approximate Factorization Technique）進行隱式處理。通過建立壓力修正和速度修正的關聯，壓力修正項將會最終出現在連續性方程中，然后通過重啟線性GMRES （m）算法［39］求解壓力場。本文所采用方法的優勢是壓力場和速度場在交互運算的過程中具有良好的兼容性，因為壓力修正公式有“精確”的處理形式。由于這個原因，可以省去壓力修正算法中弛豫時間中的壓力參數，唯一需要用戶給定的參數就是局部（或全局）時間步長，而時間步長理論上可以取無限大的值。不包含矩陣的壓力求解方法和交錯網格離散化相結合，允許壓力場的迭代計算在沒有壓力邊界條件的情況下進行。值得注意的是，在過去的Navier-Stokes方程中GMRES算法已被用來加快標準SIMPLE壓力修正算法的運算速度［40］。

采用速度一壓力的交錯耦合本身不足以抑制壓力振蕩。因此，一個四階壓力耗散項添加到連續性方程中同樣被不包含矩陣的求解器進行隱式處理。這一項代表質量流量，不會改變連續性方程的二階精度。為了防止在高雷諾數的速度場中的奇偶解耦合可以采用乘方格式對動量通量進行離散［41］。

對于湍流問題的處理，本方法采用了三種湍流模型：Baldwin-Lomax［42］的代數模型、Jones和Launder［43］的k-ε模型和Wilcox［44］的k-ω模型。當遇到粗糙網格時，這三種模型都可適應“壁面函數”邊界條件。　目前還沒有考慮層流到湍流的過渡問題和低雷諾數的影響。然而，層流到湍流的過渡流動問題可以通過在固體壁面上設置過渡位置（強制過渡或自由過渡）來進行處理。

本方法適用于二維定常和非定常翼型繞流問題的研究。在這個例子中，非定常可能是由隨時間變化的來流條件或機翼變槳距運動引起的。當機翼變槳的時候，可以使用隨翼型運動的相對參考系來研究流動問題，此時，需要將相應的慣性項添加入動量方程中，同時也要對遠場的邊界條件進行重新定義。但是，為了簡化問題，本文將湍流視作準穩態問題來研究，沒有對湍流模型進行特別處理。由于這個原因，本方法對湍流的處理僅僅適用于定常的、不存在大量流動分離的情況。

2.3.1　控制方程和湍流模型

變槳翼型的雷諾平均Navier-Stokes方程可通過無量綱的矢量形式表達。相應地，可以將此方程運用在不同位移（比如軸向位移）的情況下。控制方程建立在隨翼型運動的相對坐標系上，翼型本身的形變不在考慮范圍內。考慮到翼型在絕對坐標系上進行變槳距運動，不可壓縮Navier-Stokes的連續性方程和動量方程為

絕對坐標系中的速度矢量V和相對坐標系中的速度矢量W相對應，公式為

式中：Ω是旋轉速度矢里（變槳距速度）；r代表位置矢量；下標R代表相對坐標系。

速度被參考速度W∞歸一化為尤量綱參數，壓力（或者說壓力除以密度項ρ）被歸一化，VT是被V∞歸一化的湍流運動黏度，Re是雷諾數。長度被翼型弦長歸一化。

動量方程中包含有翼型變槳距所產生的三個慣性項，代表Coriolis力、角加速度、離心力。當Ω等于0時，W與V相等，公式將退化為標準Navier-Stokes方程組。

本算法中采用了三種湍流模型，分別是Baldwin-Lomax代數模型、k-ε湍流模型和k-ω湍流模型。Baldwin-Lomax代數模型嚴格遵循Baldwin和Lomax［42］介紹的使用規則，在此不再贅述。

Jones和Launder［43］提出的高雷諾數的標準k-ε湍流模型在此作簡要介紹。k-ε方程以如下無量綱形式進行求解

式中：W代表相對速度；Re代表雷諾數。

Pk表達式為

湍流運動黏度vT為

對于標準k-ε模型，相關常數為Prk=1，Pre=1.3，Ce1=1.44，Ce2=1.92，Cm=0.09。在沒有低雷諾數的影響下，對于k和ε方程可以通過經典的壁面函數（Wall-Function）確定。

除了k-ε模型，由Wilcox［44］提出來的k-ω模型也運用到程序中。k-ω的方程為

在k-ω模型中，相關常數為α=5/9，β=3/40，β*=9/100，σ=1/2，σ*=1/2。k-ω壁面條件通過壁面函數（wall function）來表示，當網格質量較好的時候也可直接給定壁面數值。當網格質量較好時，k和ω壁面條件為

式中：Δy是壁面到相鄰點的距離。

三種湍流模型都被內置于程序中，用戶可以在“wall function”邊界條件的選項中進行選取。對于壁面函數的運算方式，有以下步驟。

第一步：取平行于壁面的速度分量和壁面相鄰位置處的距離△y計算出翼型表面每個網格節點處的摩擦速度UT，摩擦速度UT的求解利用非線性系統的牛頓迭代法

如果則w/UT =κy+，否則w/UT=1/κln （Ey+ ）

式中：K=0.41，是馮卡門常數；的下標Is代表翼型表面的層流邊界層的邊緣。當使用Baldwin-Lomax模型時，摩擦速度是壁面渦量的函數。

第二步：當采用雙方程的k-ε模型或k-ω模型時，壁面相鄰處k，ε或ω的值通過下面的公式進行迭代計算。y+定義為

第三步：在動量方程中添加剪切應力。為簡化計算，剪切應力是在壁面中段節點位置的湍流運動黏度中人為地添加一個數值實現的。這個人為添加進去的壁面數值隨后通過線性外插值法進行迭代運算。

中段節點的數值如下

2.3.2　線性化與時間步長

對于絕對來流速度為V（假設Ω=0），下面將介紹層流非定常Navier-Stokes方程的線性化和時間步長算法。在此基礎上可以很容易地將此線性化和時間步長的算法擴展到湍流和變槳距情況。

考慮一個迭代過程，令n代表當前時間點，通過一階向后歐拉方法對時間導數項進行離散，同時將對流項在n+1附近線性化，由*表示，見圖2-1。

由此可得一組隱性時間離散的非定常Navier-Stokes方程為

在此引入的速度和壓力修正項δV和δp代表n+1和*時刻的差值。上式Navier-Stokes方程被寫成準定常的形式，還需要進行額外的迭代使n+1時刻收斂（δV，δp →0）。將速度修正和壓力修正引入動量方程，表達式為

式中：N代表*時刻的非定常（標量）Navier-Stokes算子；R代表*時刻的非定常殘差。

將上式引入連續性方程，可以得到“準確的”壓力修正公式

由于N-1的存在，壓力修正算子▽·｛N-1▽（）｝的解析表達式可以消去。然而在某些簡單情況下，例如線性Stokes問題，運用這個算子可以得到一些有用結論。例如，當公式中沒有對流項時

前一種情況適用于Re＞＞1的情況，壓力修正算子退化成Laplace算子，這也解釋了為什么Poisson方程經常被用在壓力修正算法中。這個結論只在忽略對流項影響的前提下有效。其他的算法中一般通過N （）矩陣的對角線數值進行近似

上述壓力修正算子的任一近似方法都對壓力修正項和速度修正項的兼容性有影響，同時也損害迭代過程的收斂性。本論文的想法是，沒有必要對壓力修正算子進行近似，也無需獲得解析表達式，因為式（2-31）可以通過不包含矩陣的共軛梯度方法進行求解。由于壓力修正算子是非對稱的，在此處選用重啟GMRES （m）共軛梯度法，因為這種方法具有較好的收斂性。值得注意的是，這種方法必須進行矩陣N求逆來保持速度場和壓力場的交互迭代。為了簡化數值運算，矩陣N由因式分解后的近似矩陣M代替。MSIP算法用來實現矩陣N的不完全LU因式分解。

下列過程實現不可壓縮Navier-Stokes方程的數值運算過程。

第一步：通過N的值初始化*的值。

第二步：運算矩陣N且對其因式分解，得到矩陣M。

第三步：通過下式計算中間速度場。

第四步：通過重啟型GMRES （m）算法對壓力修正式（2-31）進行求解，得到δp，更新壓力場。

第五步：通過下式更新速度場。

第六步：將*的值替代n+1的值，返回第二步迭代直到非定常動量和連續性殘差滿足收斂條件。

第七步：如果完成收斂，繼續進入下一個時間點，從步驟一開始。

對于定常流動的運算，控制方程在每個時間點并不嚴格收斂。在每個時間步長內，一般需要GMRES （m）的7到20次內部迭代才能夠實現最終收斂，其中m代表選定的Krylov子空間的維度。顯然，對于非定常流動問題的運算需要使用內部迭代。根據選定時間步長的不同，內部迭代的次數在5到20次左右。需要注意的是，算法中耗時的部分是速度迭代以及矩陣N的因式分解，因為在每一個內部迭代中都要重復一次。當矩陣M保存之后，壓力修正公式的計算都能快速完成。由于壓力和速度的更新是兼容的，為了實現結果的收斂一般不需要外部迭代。

在雙方程的湍流運算中，湍流模型k-ε或k-ω一般使用類似的方式進行線性化和離散化。

2.3.3　坐標變換

控制方程在計算域的空間中進行離散，計算域是由貼體空間坐標變換而來。（ξ，η）作為二維歐幾里得空間（x，y）的貼體參數進行坐標變換x=x （ξ，η），y=y （ξ，η），而（x，y）表示位置矢量r的笛卡爾坐標。共變量為gi，逆變量為gj，i，j=1，2是（ξ，η）的基，其定義如下，下標代表偏導數。

這些基是正交的，滿足方程，其中是Kronecker符號。j與關系為

由于正交性條件，有（gij）=1/（gij）。Jacobian行列式為

速度V的散度和對流—擴散算子N（）在（ξ，η）坐標系下表示為

式中：U和V是速度的逆變量。表示為

（u，v）是笛卡爾速度分量。壓力梯度項可表述為

2.3.4　離散化

運用有限體積法，控制方程可以在（ξ，η）坐標系的均勻網格上進行離散，此時Δξ=Δη=1。令i和j分別作為ξ和η方向上網格節點的下標，笛卡爾速度分量被存儲在網格節點（i，j）上，壓力在網格中心位置（i+1/2，j+1/2）進行運算。動量方程的控制體積的中心與網格節點重合，而連續性方程的控制體積的中心與網格細胞中心位置重合，見圖2-2。

由于通量（JU）和（JV）在網格細胞中線性分布，速度V的散度如下

動量方程中的對流—擴散算子可以通過Patankar提出的乘方格式進行離散。令

根據乘方格式，可得

其中，Peclet數可以定義為Pe=C/D。對流項和擴散項為

混合的導數系數如ai+1，j+1與1/Re （Jg12）i+1/2j相關。矩陣N （）的對角線值，可以由下式給出

值得注意的是，動量方程的殘差也可以運用類似的方法進行偏導離散化。

2.3.5　速度-壓力耦合

通過計算離散div.grad （）算子可以發現，速度—壓力交錯方法本身不足以抑制壓力震蕩。需要在連續性方程中增加一個四階線性壓力耗散項，式（2-31）由下式代替

式中：L4 （）代表ξ和η方向上四階導數的和；ε是用戶自定義的一階參數。