李 芃

(銅陵學(xué)院 電氣工程學(xué)院， 安徽 銅陵 244000)

外流是浸沒在流體中物體外圍的流動，流體流經(jīng)翼型表面會在垂直于表面的方向產(chǎn)生壓力，在平行于表面的方向產(chǎn)生剪切力。壓力與剪切力合力在流動方向的分量稱為阻力，在垂直于流動方向上的分量稱為升力，它們的研究對翼型設(shè)計具有指導(dǎo)意義。傳統(tǒng)方法可通過風(fēng)洞實驗測量出阻力和升力[1-2]，但風(fēng)洞實驗準(zhǔn)備周期較長，且耗資不菲。隨著計算機硬件的發(fā)展和數(shù)值模擬技術(shù)的進步，可采用計算流體力學(xué)軟件CFD求解流體運動控制方程獲得阻力和升力。該方法周期短、效率高、花費少，目前已成為研究翼型外流場的主要方法[3-4]。無量綱阻力系數(shù)Cd和升力系數(shù)Cl是用于表征物體阻力和升力特性的關(guān)鍵系數(shù)，本文以NACA0012翼型為研究對象，通過CFD軟件Ansys-fluent分析了翼型后緣形狀、湍流模型和遠場距離在不同攻角下對Cd和Cl數(shù)值仿真的影響，以期對翼型外流場的數(shù)值仿真研究有所借鑒。

1 仿真基礎(chǔ)

根據(jù)文獻[5-7]提供的自由來流條件和試驗數(shù)據(jù)分析不同參數(shù)對翼型外流場關(guān)鍵系數(shù)數(shù)值仿真的影響，分別采用Airfoil tools(UIUC Airfoil coordinates database)和NACA 4 digit airfoil generator設(shè)計了2種NACA0012二維翼型，翼型弦長L=1 m。根據(jù)Airfoil tools提供的132個數(shù)據(jù)點設(shè)計的翼型后緣是鈍型，采用3D curve導(dǎo)入Ansys DesignModeler時可看到后緣處沒有閉合，如果不手動將后緣兩點連接，而是直接生成二維平面，后緣處會形成類似橢圓形的鈍型，導(dǎo)致在mesh中劃分網(wǎng)格失敗。采用NACA 4 digit airfoil generator，選中close trailing edge設(shè)計包含200個數(shù)據(jù)點尖后緣翼型。2種模型如圖1所示，可以看到數(shù)據(jù)點位置存在明顯不同。試驗雷諾數(shù)Re=6 × 106，對外部流動，沿表面位置分布的雷諾數(shù)大于5 × 105可認(rèn)為流動狀態(tài)為湍流，因此要考慮湍流對翼型外流場的影響。

1.1 湍流模型

湍流是一種流體流動狀態(tài)，此時的流體做不規(guī)則運動。描述流體運動的NS方程是封閉方程組，但是直接數(shù)值模擬(DNS)需要非常大的計算資源，難以在實際中得到應(yīng)用。目前普遍采用雷諾平均(RANS)方法對湍流脈動進行時間項的平均化處理，通過簡化時間項降低對計算資源的要求。但同時增加了新的缺少控制方程的未知量導(dǎo)致NS方程不封閉，需要選擇合適的湍流模型來封閉描述湍流的方程組。選擇合適的模型描述湍流是數(shù)值仿真的關(guān)鍵，對大多數(shù)數(shù)值仿真應(yīng)用，Spalart-Allmaras、SST K-ω和K-ε是被廣泛接受的較為精確的RANS湍流模型。

(a)兩種不同的翼型模型示意圖 (b)數(shù)據(jù)點位置的近景圖

1.2 網(wǎng)格劃分

黏性是流體內(nèi)部黏滯性的度量，所有流體流動均涉及某種程度的黏性效應(yīng)。黏性導(dǎo)致直接與翼型壁面接觸的流體被粘在表面無滑移，即無滑移條件。由于無滑移條件的作用，在臨近翼型壁面非常小的法向距離內(nèi)，流體速度會從較大的數(shù)值減少到與壁面的相同速度，導(dǎo)致壁面的法向速度具有非常大的梯度，形成速度剖面，黏性效應(yīng)產(chǎn)生的速度梯度顯著的區(qū)域稱為邊界層。邊界層是物體表面附近薄的粘性層，包含了所有的粘性效應(yīng)，而邊界層外面的流動是無粘的。相關(guān)系數(shù)的計算與粘性邊界層息息相關(guān)，為了使數(shù)值仿真的結(jié)果更加貼近真實情況，翼型臨近壁面的網(wǎng)格劃分十分關(guān)鍵，尤其是第一層網(wǎng)格高度yH。對翼型臨近壁面邊界層參數(shù)的求解一般有2種方法：(1)求解粘性子層，(2)利用壁面函數(shù)。Spalart-Allmaras(SA)、SST K-ω湍流模型采用NS方程求解黏性子層的流動物理量分布，不采用壁面函數(shù)，需要滿足y+≤1。K-ε湍流模型采用壁面函數(shù)，需要滿足30 ≤y+≤300。本文采用求解粘性子層的辦法，基于Spalart-Allmaras和SST K-ω湍流模型進行翼型外流場的網(wǎng)格劃分[8-9]，需要滿足y+≤1以確定第一層網(wǎng)格的高度。雷諾數(shù)Re計算公式如公式(1)所示，其中L是特征長度，μ是動力黏度，U是自由來流速度，ρ是空氣密度。

(1)

Cf是壁面摩擦系數(shù)，計算公式如下：

Cf=[2lg(Re)-0.65]-2.3

(2)

壁面切應(yīng)力τW的計算公式如下：

(3)

根據(jù)壁面切應(yīng)力可計算摩擦速度uτ：

(4)

y+的計算公式如下，其中yP是第一層網(wǎng)格中間位置到壁面的距離：

(5)

根據(jù)公式(1)～公式(5)可計算出yP，第一層網(wǎng)格高度yH=2yP。本文的試驗參數(shù)如下：Re=6 × 106，L=1 m，μ=1.81 × 10-5kg·m-1·s-1，ρ=1.225 kg·m-3，U=88.65 m·s-1，代入上述公式可計算得到y(tǒng)H=4.4 × 10-6m。

常用的網(wǎng)格劃分工具有mesh、ICEM CFD、pointwise等，本文采用Ansys內(nèi)置的ICEM CFD。官方建議翼型外流場數(shù)值仿真的遠場距離一般應(yīng)12～20倍于翼型的弦長L，分別選取了12 m、16 m和20 m繪制C型外流場并進行網(wǎng)格劃分，外流場輸入邊界和輸出邊界命名為inlet和outlet，翼型表面命名為airfoil。采用ICEM CFD對設(shè)計的2種翼型進行網(wǎng)格劃分，考慮到后緣具備鈍型和尖型2種形式，采用了類似的網(wǎng)格分塊形式：對鈍型后緣翼型，應(yīng)用類似“魚尾”的形式劃分網(wǎng)格，保留后緣處的塊。對尖后緣翼型，折疊后緣處的塊。由于邊界層第一層的網(wǎng)格高度對Cd等參數(shù)的計算有直接影響，所以在臨近翼型壁面處再進行一次分塊。鈍后緣和尖后緣處的網(wǎng)格劃分如圖2所示。

(a)鈍后緣 (b)尖后緣

網(wǎng)格劃分的難點在于如何在較大的外流場范圍內(nèi)保證較好的網(wǎng)格質(zhì)量，且網(wǎng)格數(shù)量不能太大。評價參數(shù)aspect ratio和determinant是2個重要指標(biāo)，為了滿足y+≤1的條件，yH計算數(shù)值為4.4 × 10-6m，極易使aspect ratio的最大數(shù)值超過10 000，進而在fluent仿真環(huán)境下出現(xiàn)浮點溢出或計算發(fā)散等錯誤導(dǎo)致仿真失敗。對本文二維翼型和遠場條件而言，測試表明應(yīng)盡量使劃分網(wǎng)格aspect ratio的最大值在5 000以內(nèi)，determinant最小值大于0.8。為了保證數(shù)值仿真結(jié)果的可靠性進行了網(wǎng)格無關(guān)性測試，生成的網(wǎng)格數(shù)量為580 000。對鈍后緣，12 m、16 m和20 m的遠場距離下aspect ratio的最大值和determinant的最小值分別為(2 200,0.894)、(2 920,0.89)和(3 800,0.863)。對尖后緣，12 m、16 m和20 m的遠場距離下aspect ratio的最大值和determinant的最小值分別為(1 680,0.839)、(2 280,0.886)和(3 640,0.875)。

2 數(shù)值仿真

計算機仿真環(huán)境參數(shù)如下：CPU AMD Ryzen 7 5800X 8-core processor，RAM DDR4 64G，操作系統(tǒng)WIN10 64專業(yè)版，Ansys 2020 R2版本，F(xiàn)luent仿真環(huán)境采用雙精度、8核并行處理方式。啟動Ansys-fluent，采用press-based壓力基求解器進行穩(wěn)態(tài)計算，輸入邊界inlet類型為velocity-inlet，velocity specification method選擇magnitude and direction，攻角為0°時來流方向x分量為cos(0°)=1、y分量為sin(0°)=0，攻角為10°時來流方向x分量為cos(10°)=0.984 8、y分量為sin(10°)=0.173 6，攻角為15°時來流方向x分量為cos(15°)=0.965 9、y分量為sin(15°)=0.258 8。輸出邊界outlet類型為pressure-outlet，除turbulence外其他保持默認(rèn)參數(shù)。翼型表面airfoil類型為wall，選擇no slip無滑移條件。參考值的設(shè)置影響Cl和Cd的計算，因為是2D模型，depth默認(rèn)為1，area=depth*length，length是特征值長度。對二維翼型，length =L=1，所以area=1。流速和動力黏度取試驗數(shù)值，其他參數(shù)取默認(rèn)值。采用coupled耦合算法，flow選擇二階迎風(fēng)型。設(shè)置Cl和Cd的reports以便分析，force vector的取值與攻角直接相關(guān)：當(dāng)攻角為0°，Cl的x=-sin(0°)=0、y=cos(0°)=1，Cd的x=cos(0°)=1、y=sin(0°)=0；當(dāng)攻角為10°，Cl的x=-sin(10°)=-0.173 6、y=cos(10°)=0.984 8，Cd的x=cos(10°)=0.984 8、y=sin(10°)=0.173 6；當(dāng)攻角為15°，Cl的x=-sin(15°)=-0.258 8、y=cos(15°)=0.965 9，Cd的x=cos(15°)=0.965 9、y=sin(15°)=0.258 8。初始化選用基于拉普拉斯方程的hybrid非均勻方法，迭代5 000次。因NACA0012是對稱翼型，攻角為0°時升力非常小，升力系數(shù)近似為0。在進行數(shù)值仿真時計算得到的升力系數(shù)量級為10-6，可近似為0。

2.1 Spalart-Allmaras模型數(shù)值仿真

湍流模型選擇Spalart-Allmaras，輸入邊界inlet的turbulence specification method選擇turbulent viscosity ratio，turbulent viscosity ratio取值1，空間離散算法的相關(guān)參數(shù)都選擇二階迎風(fēng)型。翼型外流場遠場距離為12 m，攻角分別為0°、10°和15°時的數(shù)值仿真結(jié)果如下：(1)對鈍后緣，攻角為0°，Cd=0.008 155 016 4，Cl=0；攻角為10°，Cd=0.014 346 259，Cl=1.072 940 1；攻角為15°，Cd=0.024 874 452，Cl=1.523 896 3。(2)對尖后緣，攻角為0°，Cd=0.008 094 270 9，Cl=0；攻角為10°，Cd=0.014 131 926，Cl=1.067 112 5；攻角為15°，Cd=0.024 498 527，Cl=1.518 368。

翼型外流場遠場距離為16 m時，攻角分別為0°、10°和15°時的數(shù)值仿真結(jié)果如下：(1)對鈍后緣，攻角為0°，Cd=0.008 195 890 1，Cl=0；攻角為10°，Cd=0.013 932 938，Cl=1.077 611 6；攻角為15°，Cd=0.024 260 434，Cl=1.525 901。(2)對尖后緣，攻角為0°，Cd=0.008 117 424 6，Cl=0；攻角為10°，Cd=0.013 799 076，Cl=1.071 497；攻角為15°，Cd=0.024 209 021，Cl=1.516 862。

翼型外流場遠場距離為20 m時，攻角分別為0°、10°和15°時的數(shù)值仿真結(jié)果如下：(1)對鈍后緣，攻角為0°，Cd=0.008 158 128 6，Cl=0；攻角為10°，Cd=0.013 608 489，Cl=1.078 626 4；攻角為15°，Cd=0.023 542 659，Cl=1.529 061 9。(2)對尖后緣，攻角為0°，Cd=0.008 094 481 7，Cl=0；攻角為10°，Cd=0.013 821，Cl=1.067 038；攻角為15°，Cd=0.025 963 382，Cl=1.493 103 7。

2.2 SST K-ω模型數(shù)值仿真

湍流模型選擇SST K-ω，輸入邊界inlet的turbulence specification method選擇intensity and viscosity ratio，turbulence intensity取值0.052%，turbulent viscosity ratio取值0.009，空間離散算法的相關(guān)參數(shù)都選擇二階迎風(fēng)型。翼型外流場遠場距離為12 m，攻角分別為0°、10°和15°時的數(shù)值仿真結(jié)果如下：(1)對鈍后緣，攻角為0°，Cd=0.008 060 052，Cl=0；攻角為10°，Cd=0.014 627 732，Cl=1.062 077 3；攻角為15°，Cd=0.026 229 661，Cl=1.485 931 8。(2)對尖后緣，攻角為0°，Cd=0.007 993 417，Cl= 0；攻角為10°，Cd=0.014 367 92，Cl=1.054 641 4；攻角為15°，Cd=0.025 656 593，Cl=1.478 501 9。

翼型外流場遠場距離為16 m，攻角分別為0°、10°和15°時的數(shù)值仿真結(jié)果如下：(1)對鈍后緣，攻角為0°，Cd=0.008 072 766，Cl=0；攻角為10°，Cd=0.014 215 899，Cl=1.067 466 4；攻角為15°，Cd=0.025 651 792，Cl=1.488 497 1。(2)對尖后緣，攻角為0°，Cd=0.008 016 159 6，Cl=0；攻角為10°，Cd=0.014 027 672，Cl=1.057 444 5；攻角為15°，Cd=0.025 418 779，Cl=1.474 850 7。

翼型外流場遠場距離為20 m，攻角分別為0°、10°和15°時的數(shù)值仿真結(jié)果如下：(1)對鈍后緣，攻角為0°，Cd= 0.008 195 486，Cl=0；攻角為10°，Cd= 0.013 888 95，Cl=1.063 285 3；攻角為15°，Cd=0.024 878 235，Cl=1.487 204 1。(2)對尖后緣，攻角為0°，Cd=0.007 992 934，Cl=0；攻角為10°，Cd=0.013 878 605，Cl=1.056 018 2；攻角為15°，Cd=0.026 274 862，Cl=1.457 931。

3 仿真結(jié)果分析

3.1 湍流模型

采用Spalart-Allmaras湍流模型如圖3(a)所示，12 m/16 m/20 m的遠場距離+鈍后緣/尖后緣的翼型配置在不同攻角條件下對Cl的仿真效果均較好，誤差率最大不超過2%。20 m遠場距離+鈍后緣翼型配置對Cd的仿真效果最好，誤差率最大不超過24%。比較不同參數(shù)配置下Cd和Cl誤差率的均值和方差，鈍后緣翼型+20 m遠場距離的Cd誤差率數(shù)值最佳，是(0.126 868 385,0.013 330 791)，尖后緣翼型+16 m遠場距離的Cl誤差率數(shù)值最佳，是(0.001 681 592,1.756 76 × 10-6)?？紤]到Cd誤差率量級大于Cl誤差率，鈍后緣翼型+20 m遠場距離的仿真效果最接近真實情況, 在0°、10°和15°攻角下的系數(shù)誤差率分別是(0.84%,0.00%)、(13.31%,0.74%)、(23.91%,1.07%)，仿真效果最接近真實情況。采用SST K-ω湍流模型如圖3(b)所示，12 m/16 m/20 m的遠場距離+鈍后緣/尖后緣的翼型配置在不同攻角條件下對Cl的仿真效果均較好，誤差率最大不超過4%。20 m遠場距離+鈍后緣翼型配置對Cd的仿真效果最好，誤差率最大不超過31%。比較不同參數(shù)配置下Cd和Cl誤差率的均值和方差，結(jié)合Cd誤差率量級大于Cl誤差率，鈍后緣+20 m遠場距離的誤差率數(shù)值相對最佳，在0°、10°和15°攻角下的系數(shù)誤差率分別是(1.30%,0.00%)、(15.65%,0.69%)、(30.94%,1.70%)，仿真效果最接近真實情況。

(a)Spalart-Allmaras模型 (b)SST K-ω模型

3.2 翼型后緣

采用鈍后緣翼型如圖4(a)所示，12 m/16 m/20 m的遠場距離+SA/SST K-ω湍流模型配置在不同攻角條件下對Cl的仿真效果均較好，誤差率最大不超過2%。20 m遠場距離+SA湍流模型配置對Cd的仿真效果最好，誤差率最大不超過24%。比較不同參數(shù)配置下Cd和Cl誤差率的均值和方差，結(jié)合Cd誤差率量級大于Cl誤差率，SA湍流模型+20 m遠場距離的誤差率數(shù)值相對最佳，在0°、10°和15°攻角下的誤差率分別是(0.84%,0.00%)、(13.31%,0.74%)、(23.91%,1.07%)，仿真效果最接近真實情況。采用尖后緣翼型如圖4(b)所示，12 m/16 m/20 m的遠場距離+SA/SST K-ω湍流模型配置在不同攻角條件下對Cl的仿真效果均較好，誤差率最大不超過4%。16 m遠場距離+SA模型配置對Cd的仿真效果最好，誤差率最大不超過28%。比較不同參數(shù)配置下Cd和Cl誤差率的均值和方差，結(jié)合Cd誤差率量級大于Cl誤差率，SA湍流模型+16 m遠場距離的誤差率數(shù)值相對最佳，在0°、10°和15°攻角下的誤差率分別是(0.34%,0.00%)、(14.90%,0.07%)、(27.42%,0.26%)，仿真效果最接近真實情況。

(a)鈍后緣 (b)尖后緣

3.3 遠場距離

不同遠場距離下數(shù)值仿真誤差率如圖5所示。在12 m的遠場距離下，鈍后緣/尖后緣翼型+SA/SST K-ω湍流模型配置在不同攻角條件下對Cl的仿真效果均較好，誤差率最大不超過2%。尖后緣翼型+SA湍流模型配置對Cd的仿真效果最好，誤差率最大不超過29%。比較不同參數(shù)配置下Cd和Cl誤差率的均值和方差，尖后緣翼型+SA湍流模型的誤差率數(shù)值相對最佳，在0°、10°和15°攻角下的數(shù)值分別是(0.05%,0.00%)、(17.67%,0.34%)、(28.94%,0.36%)，仿真效果最接近試驗數(shù)據(jù)。在16 m的遠場距離下，鈍后緣/尖后緣翼型+SA/SST K-ω湍流模型配置在不同攻角條件下對Cl的仿真效果均較好，誤差率最大不超過3%。尖后緣翼型+SA湍流模型配置對Cd的仿真效果最好，誤差率最大不超過28%。比較不同參數(shù)配置下Cd和Cl誤差率的均值和方差，鈍后緣翼型+SA與尖后緣翼型+SA的Cd誤差率均值方差接近，分別是(0.150 021 789,0.017 470 786)和(0.142 171 489,0.018 363 59)。尖后緣翼型+SA的Cl誤差率均值方差最小，為(0.001 681 592,1.756 76 × 10-6)，所以尖后緣翼型+SA的仿真效果最接近試驗數(shù)據(jù)，在0°、10°和15°攻角下的誤差率分別是(0.34%,0.00%)、(14.90%,0.07%)、(27.42%,0.26%)。如圖5(c)所示，在20 m的遠場距離下，鈍后緣/尖后緣翼型+SA/SST K-ω湍流模型配置在不同攻角條件下對Cl的仿真效果均較好，誤差率最大不超過4%。鈍后緣翼型+SA湍流模型配置對Cd的仿真效果最好，誤差率最大不超過24%。比較不同參數(shù)配置下Cd和Cl誤差率的均值和方差，鈍后緣翼型+SA湍流模型的數(shù)值最小，在0°、10°和15°攻角下的誤差率分別是(0.84%,0.00%)、(13.31%,0.74%)、(23.91%,1.07%)，仿真效果最接近試驗數(shù)據(jù)。

(a)12 m (b)16 m

(c)20 m

綜上所述，從湍流模型的角度出發(fā)Spalart-Allmaras的表現(xiàn)優(yōu)于SST K-ω。從不同翼型后緣形狀的選擇出發(fā)，翼型形狀對仿真數(shù)值有明顯影響。鈍后緣翼型的數(shù)值仿真效果隨著遠場距離的增大而變好，尖后緣翼型的數(shù)值仿真效果先隨著遠場距離的增大變好，達到一個峰值(本文是16 m)后再隨著距離的增大變差。從不同遠場距離的選擇出發(fā)，在距離相對較小時，仿真誤差率相對較大，尖后緣的仿真效果優(yōu)于鈍后緣。隨著距離的增大，誤差率降低，尖后緣和鈍后緣的仿真效果逐漸接近，達到某個數(shù)值(本文是20 m)后鈍后緣的仿真效果優(yōu)于尖后緣。

4 結(jié) 語

對NACA0012翼型，分析了在攻角為0°、10°和15°條件下不同后緣形狀、遠場距離和湍流模型對翼型外流場Cd和Cl的數(shù)值仿真結(jié)果的影響，有如下結(jié)論：

(1)優(yōu)先考慮使用Spalart-Allmaras湍流模型。

(2)翼型形狀對仿真效果的影響隨著遠場距離的變化而變化，遠場距離較小時尖后緣仿真效果優(yōu)于鈍后緣。所以在繪制翼型外流場時，如果遠場距離較小，那么要保證后緣形狀盡量準(zhǔn)確。但是如果遠場距離較大，對后緣形狀準(zhǔn)確度的要求可以降低。遠場距離的增大雖然會導(dǎo)致網(wǎng)格質(zhì)量指標(biāo)aspect ratio增大，但只要能控制在合理的范圍內(nèi)(本文是5 000)，隨著遠場距離的增大，誤差率可進一步降低。綜合來看，在20倍的遠場距離下采用Spalart-Allmaras湍流模型+鈍后緣翼型可獲得相對最佳的數(shù)值仿真結(jié)果，相應(yīng)的誤差率分別是(0.84%,0.00%)、(13.31%,0.74%)和(23.91%,1.07%)。