橢球體定常流動粘性流場和水動力計算方法研究

于向陽 姚凌虹 孟慶昌 劉巨斌 張志宏

（1海軍航空大學青島校區 青島 266000）（2海軍工程大學 武漢 430033）

1 引言

潛艇等水下航行體粘性流場和水動力對其操縱性預報至關重要。而潛艇等水下航行體作操縱運動時，其水動力主要集中在主體部分，同時復雜的三維流動（潛艇非線性水動力的主要來源）也主要發生在主體上，所以主體水動力預報的準確程度備受關注。6：1橢球體類似于潛艇等水下航行體的主體，其外形雖然簡單，但其帶攻角的繞流流動卻是一個復雜、典型的三維分離和渦流流動，而對其流動預報的結果可以間接檢驗數值計算模型對潛艇等水下航行體操縱性預報的能力［1～4］。

本文應用CFD軟件對6：1橢球體粘性流場和水動力進行數值模擬。采用基于非結構化網格的有限體積法離散控制方程和湍流模式，速度-壓力耦合采用SIMPLE（Semi-implicitMethodforPressure-LinkedEquations）方法處理，對流項采用二階迎風格式，采用不同湍流模型，包括k-ε湍流模型，SA一方程湍流模型，k-ω湍流模型，SST兩方程湍流模型，得到5°、10°、15°、20°、25°、30°這6個攻角的水動力數值計算結果，同時分析了不同計算域范圍，不同網格密度等對數值計算結果的影響，并與文獻［5］中的實驗結果進行了比較。

2 數值計算方法

2.1 數值離散方法

采用基于非結構化網格的有限體積法，數值求解定常的RANS方程，在SA模式中，生成項所采用的形式為

其中

2.2 計算對象與邊界條件設置

計算對象為長細比6：1的橢球體，采用直角坐標系，其原點位于橢球體的幾何中心，x軸與橢球體的對稱軸重合，方向與無攻角時來流方向一致，z軸豎直向上，y軸水平。考慮到對稱性，僅對半橢球體流動進行數值計算；計算域為半圓柱體區域，范圍：-a≤x/L≤b（a，b取值在各節中不同），-2≤z/L≤2，0≤ y/L≤2，其中L為橢球體長度，取值L=1.37，計算域如圖1所示。

圖1 模型的邊界條件設置

計算域的邊界條件：半圓柱面的下半部分AEFD（z≤0）和半圓柱體的左半圓面AEB（x/L=a處）定義為速度入口，給定速度Uo，方向沿x軸正向；半圓柱面的上半部分EBCF（z＞0）和半圓柱體的右半圓面CFD（x/L=b處）定義為壓力出口，考慮到在迭代過程中，有可能存在數值上的反向流動，設置出口回流方向的方式為垂直于出口面；物面為無滑移壁面，速度分量和湍動能為零；對稱面上，法向速度分量為零，平行于對稱面的速度分量的法向導數和所有標量的法向導數為零。

在所定義的坐標系下，規定來流的z方向速度分量為正時為正攻角，計算雷諾數Re=4.2×106，收斂準則為所有求解變量的無量綱殘差下降4個量級。

3 邊界層處理

3.1 邊界層理論

邊界層問題表現為近壁湍流問題，在工程實際中流體具有粘性，因而流體在物面上的速度相對于物面應為零，即無滑移條件。而在緊貼物面的薄層區域內，分子輸運現象起主導作用，致使一般湍流模型的很多基本假設不再成立；同時近壁處的速度、壓力等變量沿物面法向方向梯度很大。在數值計算中，為了獲得可靠的數值計算結果，不得不考慮上述問題。一種解決方法是在物面附近，使用壁面函數法來處理，但針對具有高雷諾數的流動問題；另一種解決方法是在物面附近，增加網格密度，以取得較密的節點，但這樣會增加計算量。

在數值計算中，引入無量綱值u+和 y+，進行近壁處網格的布置。

其中，τω是壁面剪切力，y是沿物面法向的距離。

如采用近壁模型法，要求第一層網格y+＜5，而且在法向雷諾數Re＜200的范圍內，至少布置10層網格，但這勢必要大大增加計算量。

圖2 近壁處理方法

3.2 網格劃分

計算域網格由邊界層網格和體網格兩部分組成。

CFD對網格選取有一定要求，一方面考慮到粘性效應近壁面采用較密的貼體網格，另一方面網格的疏密程度與流場參數（如壓力、速度等）的變化梯度趨于一致［6～11］。為捕捉壁面附近的強分離流動，使橢球體表面附近的壓力，速度等物理參數有較好的模擬效果，采用三角形網格劃分底層物面，沿垂直于壁面方向，邊界層進行加密處理，以求得高質量的貼體網格，如圖3所示；體網格采用基于四面體的非結構化網格形式，以適應橢球體表面的曲率變化［12～15］。估算 y+并確定第一層網格控制點離壁面的距離△y/L，控制在1×10-5。

圖3 邊界層網格

4 不同湍流模型下的粘性流場和水動力數值計算結果比較

湍流模型的實質是為RANS方程中出現的雷諾應力尋求計算方法。考慮到求解問題的復雜性及計算條件的局限性，必須選擇合適的湍流模型進行數值計算。

4.1 湍流模型的選擇

由于湍流流動問題的差異性，沒有一種湍流模型能夠對所有工程實際問題具有通用性，即使比較常用的湍流模型也只是對某一類或某一些特定問題適用，因此需要對所研究的特定問題選取不同湍流模型，以驗證其對計算對象的適用性。對湍流模型的選擇通常注意以下幾個方面：1）所計算問題的特殊性，物理問題是否可以簡化；2）所擁有的計算資源，對計算時間的限制；3）對計算結果的要求，如計算精度的要求等［16～18］。

本文采用的湍流模型如表1所示，SA模型適合應用于存在翼形、壁面邊界層等場合的湍流流動問題，不適合射流類等具有自由剪切流動的湍流問題，本文中的橢球體繞流問題屬于其適用范圍；標準k-ε模型應用較為廣泛，具有較好的穩定性，同時也具有較高的計算精度，適合高雷諾數湍流流動，但對旋流等各向異性較強的流動具有局限性；標準k-ω模型，包含低雷諾數影響和可壓縮性影響，適用于受壁面限制的流動附著邊界層湍流流動計算；剪切應力輸運（SST），綜合了k-ω模型在近壁區計算和k-ε模型在遠場計算的優越性，適用于帶逆壓梯度的流動計算。湍流模型的計算量與求解的方程數量有關，進行數值模擬時，需要考慮計算資源的合理利用；各湍流模型的計算速度依次為SA一方程模型，k-ε和k-ω雙方程模型。

表1 采用的湍流模型

4.2 計算結果分析

本節適度選取網格間距（網格參數如表2），對不同湍流模型的水動力數值計算結果進行比較，結果包括力矩系數、升力系數CL和阻力系數隨攻角的變化情況，并與國外實驗結果進行了比較，結果表明表1中列舉的幾種湍流模型都能夠反映流場的水動力和粘性分離特性。

表2 網格參數

圖4～圖6（圖中exp為實驗值，下文不再另行說明）分別給出了橢球體升力系數CL和阻力系數CD和力矩系數CM隨攻角的變化情況，其中力矩的參考點取在橢球體的中心。與實驗值比較，升力系數均小于實驗值，但SA（Cυ=20）和k-ω模型的計算值與實驗值較為接近，且較好地模擬了實驗值的趨勢。

圖4 升力系數與攻角的關系

圖5 阻力系數與攻角的關系

依據實驗本身曲線的總體趨勢來分析，隨著攻角的增加，非線性作用明顯；在力矩系數的比較中，SST和k-ε模型與實驗值符合均較好，SA（Cυ=20）模型存在偏差；從阻力系數的計算結果比較中，可以看到除k-εstand模型外，其余四種模型在低攻角時，變化較為平直，大攻角時曲線陡直，非線性作用增強。綜合比較5種湍流模型，對橢球體的水動力數值計算效果，SA和k-ω模型具有較好的計算效果。

圖6 力矩系數與攻角的關系

圖7 ～圖11分別給出了橢球體20°攻角時x/L=0.272截面速度向量圖，在背風面可以清晰地看到有分離渦出現，反映了不同湍流模型下流場的粘性分離特性。

圖7 SST湍流模型20°攻角時x/L=0.272截面速度向量圖

圖8 k-ωstand湍流模型20°攻角時x/L=0.272截面速度向量圖

圖9 k-εstand湍流模型20°攻角時x/L=0.272截面速度向量圖

圖10 SA（Cυ=20）湍流模型20°攻角時x/L=0.272截面速度向量圖

圖11 SA（vorticity）湍流模型20°攻角時x/L=0.272截面速度向量圖

5 不同計算域范圍對水動力數值計算結果的影響

5.1 計算域范圍的選取

采用數值方法求解控制方程時，需要在空間區域上進行離散，以得到離散方程組。空間上離散控制方程，需要進行網格劃分，在劃分網格時，需要選取適當的計算域，防止邊界條件產生不利影響。經過查找大量的相關領域的文獻，發現對計算域范圍并沒有定量的討論，本文在參考相關文獻的基礎上選擇不同的計算域范圍進行了數值計算，并對不同計算域的數值計算結果進行了分析和比較。

表3 不同計算域參數

5.2 計算結果分析

圖12～圖14分別給出了橢球體升力系數CL和阻力系數CD和力矩系數CM隨攻角的變化情況，其中力矩的參考點取在橢球體的中心。與實驗值比較，升力系數的數值計算結果均小于實驗值，但較好地模擬了實驗值的趨勢，不同的計算域對應的數值計算結果基本一致，但隨著攻角的增加，特別是大攻角處，存在差異性，計算域選取-2.5≤＜x/L≤4.5時，已收斂；在力矩系數的比較中，不同計算域對應的計算結果與實驗值符合均較好，計算域選取-2.5≤＜x/L≤4.5時，在大攻角處，顯示出較好的模擬效果；對阻力系數的數值計算中，大攻角時曲線陡直，不同計算域都捕捉到了其非線性作用。

圖12 升力系數與攻角的關系

圖13 阻力系數與攻角的關系

圖14 力矩系數與攻角的關系

6 不同網格密度對水動力數值計算結果的影響

隨著CFD技術應用的不斷深入，許多復雜的工程問題得到解決，特別是幾何外形復雜的實際問題更多的提到日程上來，在這里我們不得不提到CFD前處理軟件，為生成高效率的滿足計算精度的網格所做的貢獻。

6.1 網格無關解

用有限體積法離散得到的代數方程與原偏微分方程比較，存在截斷誤差。截斷誤差的大小與網格的間距（或網格數目）有關，網格間距越小（網格數目越多）截斷誤差就越小。因此就減小截斷誤差而言，網格的數目越多越好。

事實上，由于計算機內存和計算時間的限制，我們不可能無止境地減小網格間距，增加網格數目，同時，由于數值計算中總會有舍入誤差的出現，沒有必要片面地對截斷誤差作太苛刻的要求。

工程實際中，為了獲得網格無關解，劃分網格時，需要對同一個計算問題劃分多套網格，通過逐漸增加網格數目，來考察網格數目對計算結果的影響。

如果在兩個不同數目的網格上獲得的解基本相同，說明此時進一步加大網格數目，增加網格密度，已無必要，或者說我們已經得到了與網格無關的解。

6.2 網格密度的選擇

為了數值模擬得到與網格無關的收斂解，應該選取密度不同的網格進行數值計算，作比較以檢驗網格無關性。網格密度的選取一般要考慮所使用計算機的運行速度和內存容量，網格的生成往往占用整個CFD任務70%以上的工作量，所以在不影響計算精度的前提下，盡量節省計算資源，選取適度的網格密度。通過參考國內外文獻，確定以下三種網格進行數值計算（如表4所示）。

表4 不同網格密度

6.3 計算結果分析

圖15～圖17分別給出了橢球體升力系數CL和阻力系數CD和力矩系數CM隨攻角的變化情況，其中力矩的參考點取在橢球體的中心，即坐標系原點處。數值計算方法如前所述，湍流模式均采用SST模式。與實驗值比較，升力系數與力矩系數的計算值均小于實驗值，在小攻角下與實驗值符合較好，只是在30°時，不同的網格密度在數值模擬結果存在偏差，不同網格密度計算結果

圖15 升力系數與攻角的關系

總體趨勢保持一致，網格密度較小（Δx/L=0.05）的計算值在本例中并沒有表現出眾的模擬效果，網格密度取Δx/L=0.07時，已經取得較好的計算值，與網格密度為Δx/L=0.05時相比，相對偏差在1.5%左右，表明此時計算結果已收斂，網格密度選取已合理；阻力系數的比較中，不同網格密度的結果總體趨勢保持一致。

圖16 阻力系數與攻角的關系

圖17 力矩系數與攻角的關系

綜合考慮計算資源與數值計算精度的要求，網格密度取Δx/L=0.07是可行的，以相對偏差1.5%的結果換取百萬網格的代價，解決現有問題是可以考慮的。

7 結語

6：1橢球體類似于潛艇等水下航行體的主體，對其流動預報的結果可以間接檢驗數值計算模型對潛艇等水下航行體操縱性預報的能力。本文應用CFD軟件對長細比為6：1橢球體水動力進行數值模擬。采用基于非結構化網格的有限體積法離散控制方程和湍流模式，對流項采用二階迎風格式，速度-壓力耦合采用SIMPLE方法處理，得到5°、10°、15°、20°、25°、30°這6個攻角的水動力數值計算結果，并與已知的實驗結果進行了比較。

分析了不同湍流模式對數值計算結果的影響，綜合比較5種湍流模型對橢球體的水動力數值計算結果，認為SA和SST模型具有較好的計算效果。

分析了不同計算域，不同網格密度等對數值計算結果的影響。不同的計算域對應的數值計算結果基本一致，但隨著攻角的增加，特別是大攻角處，存在差異性，計算域選取-2.5≤＜x/L≤4.5時，略顯優勢；不同的網格密度在數值模擬上存在微小偏差，較小網格密度（Δx/L=0.05）的計算結果在本例中并沒有表現特別出眾的模擬效果，網格密度取Δx/L=0.07時，已經取得較好的計算結果，與網格密度為Δx/L=0.05的計算結果相比，相對偏差在1.5%左右，本文認為以相對偏差1.5%的結果換取百萬網格的代價，解決現有問題是可以考慮的。本文的計算結果和結論為后續湍流模式、計算域范圍和網格密度的研究提供了參考依據。