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

（1海軍航空大學(xué)青島校區(qū) 青島 266000）（2海軍工程大學(xué) 武漢 430033）

1 引言

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

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

2 數(shù)值計(jì)算方法

2.1 數(shù)值離散方法

采用基于非結(jié)構(gòu)化網(wǎng)格的有限體積法，數(shù)值求解定常的RANS方程，在SA模式中，生成項(xiàng)所采用的形式為

其中

2.2 計(jì)算對(duì)象與邊界條件設(shè)置

計(jì)算對(duì)象為長細(xì)比6：1的橢球體，采用直角坐標(biāo)系，其原點(diǎn)位于橢球體的幾何中心，x軸與橢球體的對(duì)稱軸重合，方向與無攻角時(shí)來流方向一致，z軸豎直向上，y軸水平?？紤]到對(duì)稱性，僅對(duì)半橢球體流動(dòng)進(jìn)行數(shù)值計(jì)算；計(jì)算域?yàn)榘雸A柱體區(qū)域，范圍：-a≤x/L≤b（a，b取值在各節(jié)中不同），-2≤z/L≤2，0≤ y/L≤2，其中L為橢球體長度，取值L=1.37，計(jì)算域如圖1所示。

圖1 模型的邊界條件設(shè)置

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

在所定義的坐標(biāo)系下，規(guī)定來流的z方向速度分量為正時(shí)為正攻角，計(jì)算雷諾數(shù)Re=4.2×106，收斂準(zhǔn)則為所有求解變量的無量綱殘差下降4個(gè)量級(jí)。

3 邊界層處理

3.1 邊界層理論

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

在數(shù)值計(jì)算中，引入無量綱值u+和 y+，進(jìn)行近壁處網(wǎng)格的布置。

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

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

圖2 近壁處理方法

3.2 網(wǎng)格劃分

計(jì)算域網(wǎng)格由邊界層網(wǎng)格和體網(wǎng)格兩部分組成。

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

圖3 邊界層網(wǎng)格

4 不同湍流模型下的粘性流場和水動(dòng)力數(shù)值計(jì)算結(jié)果比較

湍流模型的實(shí)質(zhì)是為RANS方程中出現(xiàn)的雷諾應(yīng)力尋求計(jì)算方法?？紤]到求解問題的復(fù)雜性及計(jì)算條件的局限性，必須選擇合適的湍流模型進(jìn)行數(shù)值計(jì)算。

4.1 湍流模型的選擇

由于湍流流動(dòng)問題的差異性，沒有一種湍流模型能夠?qū)λ泄こ虒?shí)際問題具有通用性，即使比較常用的湍流模型也只是對(duì)某一類或某一些特定問題適用，因此需要對(duì)所研究的特定問題選取不同湍流模型，以驗(yàn)證其對(duì)計(jì)算對(duì)象的適用性。對(duì)湍流模型的選擇通常注意以下幾個(gè)方面：1）所計(jì)算問題的特殊性，物理問題是否可以簡化；2）所擁有的計(jì)算資源，對(duì)計(jì)算時(shí)間的限制；3）對(duì)計(jì)算結(jié)果的要求，如計(jì)算精度的要求等［16～18］。

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

表1 采用的湍流模型

4.2 計(jì)算結(jié)果分析

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

表2 網(wǎng)格參數(shù)

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

圖4 升力系數(shù)與攻角的關(guān)系

圖5 阻力系數(shù)與攻角的關(guān)系

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

圖6 力矩系數(shù)與攻角的關(guān)系

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

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

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

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

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

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

5 不同計(jì)算域范圍對(duì)水動(dòng)力數(shù)值計(jì)算結(jié)果的影響

5.1 計(jì)算域范圍的選取

采用數(shù)值方法求解控制方程時(shí)，需要在空間區(qū)域上進(jìn)行離散，以得到離散方程組?？臻g上離散控制方程，需要進(jìn)行網(wǎng)格劃分，在劃分網(wǎng)格時(shí)，需要選取適當(dāng)?shù)挠?jì)算域，防止邊界條件產(chǎn)生不利影響。經(jīng)過查找大量的相關(guān)領(lǐng)域的文獻(xiàn)，發(fā)現(xiàn)對(duì)計(jì)算域范圍并沒有定量的討論，本文在參考相關(guān)文獻(xiàn)的基礎(chǔ)上選擇不同的計(jì)算域范圍進(jìn)行了數(shù)值計(jì)算，并對(duì)不同計(jì)算域的數(shù)值計(jì)算結(jié)果進(jìn)行了分析和比較。

表3 不同計(jì)算域參數(shù)

5.2 計(jì)算結(jié)果分析

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

圖12 升力系數(shù)與攻角的關(guān)系

圖13 阻力系數(shù)與攻角的關(guān)系

圖14 力矩系數(shù)與攻角的關(guān)系

6 不同網(wǎng)格密度對(duì)水動(dòng)力數(shù)值計(jì)算結(jié)果的影響

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

6.1 網(wǎng)格無關(guān)解

用有限體積法離散得到的代數(shù)方程與原偏微分方程比較，存在截?cái)嗾`差。截?cái)嗾`差的大小與網(wǎng)格的間距（或網(wǎng)格數(shù)目）有關(guān)，網(wǎng)格間距越?。ňW(wǎng)格數(shù)目越多）截?cái)嗾`差就越小。因此就減小截?cái)嗾`差而言，網(wǎng)格的數(shù)目越多越好。

事實(shí)上，由于計(jì)算機(jī)內(nèi)存和計(jì)算時(shí)間的限制，我們不可能無止境地減小網(wǎng)格間距，增加網(wǎng)格數(shù)目，同時(shí)，由于數(shù)值計(jì)算中總會(huì)有舍入誤差的出現(xiàn)，沒有必要片面地對(duì)截?cái)嗾`差作太苛刻的要求。

工程實(shí)際中，為了獲得網(wǎng)格無關(guān)解，劃分網(wǎng)格時(shí)，需要對(duì)同一個(gè)計(jì)算問題劃分多套網(wǎng)格，通過逐漸增加網(wǎng)格數(shù)目，來考察網(wǎng)格數(shù)目對(duì)計(jì)算結(jié)果的影響。

如果在兩個(gè)不同數(shù)目的網(wǎng)格上獲得的解基本相同，說明此時(shí)進(jìn)一步加大網(wǎng)格數(shù)目，增加網(wǎng)格密度，已無必要，或者說我們已經(jīng)得到了與網(wǎng)格無關(guān)的解。

6.2 網(wǎng)格密度的選擇

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

表4 不同網(wǎng)格密度

6.3 計(jì)算結(jié)果分析

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

圖15 升力系數(shù)與攻角的關(guān)系

總體趨勢保持一致，網(wǎng)格密度較?。é/L=0.05）的計(jì)算值在本例中并沒有表現(xiàn)出眾的模擬效果，網(wǎng)格密度取Δx/L=0.07時(shí)，已經(jīng)取得較好的計(jì)算值，與網(wǎng)格密度為Δx/L=0.05時(shí)相比，相對(duì)偏差在1.5%左右，表明此時(shí)計(jì)算結(jié)果已收斂，網(wǎng)格密度選取已合理；阻力系數(shù)的比較中，不同網(wǎng)格密度的結(jié)果總體趨勢保持一致。

圖16 阻力系數(shù)與攻角的關(guān)系

圖17 力矩系數(shù)與攻角的關(guān)系

綜合考慮計(jì)算資源與數(shù)值計(jì)算精度的要求，網(wǎng)格密度取Δx/L=0.07是可行的，以相對(duì)偏差1.5%的結(jié)果換取百萬網(wǎng)格的代價(jià)，解決現(xiàn)有問題是可以考慮的。

7 結(jié)語

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

分析了不同湍流模式對(duì)數(shù)值計(jì)算結(jié)果的影響，綜合比較5種湍流模型對(duì)橢球體的水動(dòng)力數(shù)值計(jì)算結(jié)果，認(rèn)為SA和SST模型具有較好的計(jì)算效果。

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