彭 聰

(成都大學(xué)生物產(chǎn)業(yè)學(xué)院，四川 成都 610106)

基于附面層的摩擦阻力數(shù)值算法研究

彭 聰

(成都大學(xué)生物產(chǎn)業(yè)學(xué)院，四川 成都 610106)

利用附面層的速度型求解壁面的空氣剪切應(yīng)力，進而求得摩擦阻力，由此計算了一維層流平板邊界層和二維層流NACA0012翼型的摩擦阻力.平板邊界層計算結(jié)果同布拉修斯理論解相比較，吻合性良好.翼型計算結(jié)果同實驗數(shù)據(jù)比較發(fā)現(xiàn)，小攻角氣流不產(chǎn)生分離的情況下，摩擦阻力值與實驗數(shù)據(jù)接近，隨著迎角增大，分離區(qū)的擴展，壓差阻力的比重增加，計算誤差明顯增加.

摩擦阻力;附面層;數(shù)值計算;空氣粘性

0 引言

從20世紀70年代末開始，在先進理論方法和成功工程實踐的基礎(chǔ)上，由于計算流體動力學(xué)(Computational Fluid Dynamics，CFD)計算方法的不斷突破和計算機技術(shù)的快速發(fā)展，CFD進入了蓬勃發(fā)展的新時期，開始越來越多地應(yīng)用到航空航天飛行器的研究和設(shè)計中［1－4］.目前，阻力計算一直是CFD中的一大難題和研究熱點，根據(jù)其產(chǎn)生機理，阻力主要分為壓差阻力和摩擦阻力2種.壓差阻力是飛行器各部分表面受到的壓力在速度方向上的合力，又分成波阻、渦阻、型阻等;而摩擦阻力是飛行器各部分表面受到的剪切力在速度方向的合力，它是由空氣的粘性而引起的.目前的CFD方法對于流場中物體表面壓力的分布能得到很好的計算結(jié)果，在壓差阻力的計算方面精度較高，但對于摩擦阻力的計算卻存在一定難度，這是因為摩擦阻力的計算對CFD有較高的要求，需要計算出飛行器合理的流場結(jié)構(gòu)，尤其需要正確捕捉壁面邊界層內(nèi)的速度型以準確計算壁面處的剪切力［5－6］.平板繞流問題是附面層流動最簡單也是最重要的情況，而對于曲面繞流，只要不發(fā)生顯著分離，摩擦阻力與平板情況相差不多.本研究采用層流平板邊界層和NACA0012翼型邊界層來驗證摩擦阻力計算的有效性及精度.

1 摩擦阻力計算方法

1.1 控制方程

本研究利用商用Fluent軟件求解N-S方程，不考慮體積力和外部熱源，其形式［6］為，

式中，Q為守恒變矢量;f，g，h分別為x，y，z坐標方向的通量.

1.2 壁面摩擦阻力

摩擦阻力由切向應(yīng)力引起，并跟附面層的速度型相關(guān).切向應(yīng)力為，

式中，μ為流體粘性系數(shù)，摩擦應(yīng)力τ與附面層間的切向速度u的變化有關(guān).

摩擦阻力為壁面摩擦應(yīng)力xF的積分，其形式為，

平板摩擦阻力系數(shù)CF根據(jù)定義為，

1.3 壁面應(yīng)力的精度

壁面應(yīng)力τw的計算采用二階精度的中心差分格式，其形式為:

2 算例及分析

2.1 層流平板附面層及摩擦阻力計算

層流平板附面層的計算模型如圖1所示.網(wǎng)格數(shù)為106×50(流向 ×法向)，遠場為4倍平板長.為準確地模擬粘性效應(yīng)，附面層內(nèi)布置了至少4個網(wǎng)格點.自由來流條件為:來流速度Ma∞=0.3，雷諾數(shù) ReL=1.0 ×105，攻角 α =0°，研究網(wǎng)格質(zhì)量對流場及摩擦阻力的影響.表1為平板摩擦阻力計算結(jié)果，并與布拉修斯解進行了比較，圖2為平板各占位的速度型.

圖1 層流平板附面層計算網(wǎng)格

表1 摩擦阻力計算結(jié)果

圖2 層流平板各站位的速度型

從圖2可以看出，平板附面層內(nèi)各站位的速度型與布拉休斯解曲線擬合得很好，由此計算得到的阻力系數(shù)值很接近.但由于阻力系數(shù)值本身為一個小量，網(wǎng)格的粗細和計算格式的粘性分辨率對計算精度都產(chǎn)生影響.從計算效率來說，對于本算例，附面層內(nèi)布置3～4個網(wǎng)格點就能很好地捕捉速度型.此對于大型網(wǎng)格計算具有指導(dǎo)意義，但對于復(fù)雜外形，網(wǎng)格數(shù)量更大的情形，還需要進一步驗證.

2.2 NACA0012翼型附面層及摩擦阻力計算

NACA0012翼型附面層的計算模型如圖3所示.網(wǎng)格數(shù)為497×100(C型網(wǎng)格)，遠場為10倍弦長.邊界層內(nèi)同樣布置了4個網(wǎng)格點.自由來流條件為:來流速度 V∞=21.1 m/s，雷諾數(shù) ReL=1.44×106，攻角分別為 0°、2°、4°、8°時，著重研究翼型繞流情況對摩擦阻力計算的影響.

圖3 NACA0012翼型計算網(wǎng)格

圖4、圖5、表2為模型在不同攻角下的計算與實驗結(jié)果.

圖4 不同攻角下的流線圖

從圖4中的流線可以看出，0°和2°攻角的流場情況較為相似，氣流均貼著翼型呈明顯的層流狀態(tài);4°攻角時，翼型對來流產(chǎn)生擾動，上翼面1/2弦長位置出現(xiàn)局部氣泡，但氣流未出現(xiàn)明顯分離，至后緣氣流再附于翼型表面，基本保持層流狀態(tài);8°攻角時，擾動進一步增強，上翼面氣泡出現(xiàn)位置向前緣移動，且出現(xiàn)明顯的分離趨勢，至后緣與下翼面繞流形成分離渦結(jié)構(gòu).

圖5 NACA0012翼型計算結(jié)果同實驗數(shù)據(jù)的比較

表2 NACA0012翼型在不同攻角下阻力系數(shù)計算值與實驗值的比較

此外，從圖4中可看出，在對稱流動時，駐點位于翼的前緣，翼型繞流為附體狀態(tài)，阻力的主要貢獻為空氣粘性影響的翼型表面很薄的附面層，表現(xiàn)為摩擦阻力，表面壓強分布對阻力影響很小.層流狀態(tài)下，摩擦阻力也很小，即使在很小的攻角情況下，只要繞流不產(chǎn)生分離，摩擦阻力的增加不明顯.隨著迎角的增大，駐點向下翼面移動，并逐漸后移，上翼面最大速度點向前緣移動.對于圓頭較厚翼型，在中等迎角下，上翼面靠近后緣的附面層在逐漸增大的逆壓梯度作用下發(fā)生局部分離.隨著攻角增加，分離區(qū)向前擴展，當迎角達到某個臨界值后上翼面的附體流動被徹底破壞，翼型表面產(chǎn)生很大的壓強差，導(dǎo)致升力降低，阻力顯著增大.由此可見，層流到完全破壞的過程，壓差阻力的貢獻越來越大.

同時，從圖5與表2可以看出，在小攻角的層流狀態(tài)下，根據(jù)摩擦阻力計算的阻力系數(shù)與實驗值接近，誤差約為10%.攻角增大后，隨著分離氣流的出現(xiàn)，壓差阻力的增加，阻力系數(shù)計算誤差也顯著增加，4°攻角時為15%，8°攻角時達40%.這主要有2方面原因:一是分離氣流的速度方向不再沿著翼型流向，法向速度分量的壓力貢獻產(chǎn)生的合力增加;二是分離氣流本身的不穩(wěn)定性，計算模型在模擬分離流動的精度不夠而造成.采用式(2)的阻力計算方法建立在精確模擬流場特性的基礎(chǔ)上，對此，需尋求更好的計算模型和計算精度來進行模擬求解.

3 結(jié)論

本研究基于對一維平板邊界層和二維NACA0012翼型采用數(shù)值方法得到的摩擦阻力計算結(jié)果，可以看出:采用低階精度格式計算低速條件下的層流邊界層能夠得到精度很高的摩擦阻力值，與理論解和實驗數(shù)據(jù)相當吻合;采用壁面摩擦阻力計算方法在層流狀態(tài)下能夠得到準確的結(jié)果，但隨著流動分離的出現(xiàn)，壓差阻力的比重增加，采用本方法得到的結(jié)果誤差增大，適用性降低.基于上述2點，本研究下一步的工作思路是研究湍流狀態(tài)下的摩擦阻力計算方法，包括選擇合適的湍流模型，采用更高精度的計算格式，等等.

［1］Tinoco E N.An assessment of CFD prediction of drag and other longitudinal characteristics［C］//39th Aerospace Sciences Meeting and Exhibit.Reno，NV:AIAA，2011:1002.

［2］Peavey C.Drag prediction of military aircraft using CFD［C］//38th Aerospace Sciences Meeting and Exhibit.Reno，NV:AIAA，2000:383.

［3］葉建，林國華，鄒正平，等.低雷諾數(shù)下二維翼型繞流的流場特性分析［J］.航空動力學(xué)報，2003，18(2):38 －45.

［4］侯志勇，石磊，聶萬勝.對稱翼型低速繞流流場特性的數(shù)值分析研究［J］.科學(xué)技術(shù)與工程，2009，9(18):5610－5613.

［5］周丹，禹建軍，閆超.層流平板摩擦阻力的數(shù)值計算［J］.北京航空航天大學(xué)學(xué)報，2007，52(6):663－667.

［6］閆超.計算流體力學(xué)方法及應(yīng)用［M］.北京:航空航天大學(xué)出版社，2007.

［7］Rice M S.Hand book of airfoil sections for light aircraft［M］.Milwaukee，WI:Aviation Publications，1976.

Study on Numerical Algorithm for Skin Friction Based on Boundary Layers

PENG Cong
(School of Biological Industry，Chengdu University，Chengdu 610106，China)

In order to obtain skin friction of the boundary layers，we adopt the air velocity model to solve the viscid stress.The skin frictions of two cases，1-D laminar flat-plate boundary layer and 2-D laminar NACA0012 airfoil，are calculated by this method.The result of the flat-plate boundary layer matches Blasius solution quite well.Comparing the result of the airfoil calculation with the experimental data，we find that the calculation precision is acceptable at a relatively small angle of attack without airflow separation.As the angle of attack increases and the separation region expands and the proportion of pressure resistance increases，the calculation errors of this method will increase significantly.

skin friction;boundary layer;numerical calculation;air viscosity

V211.3

A

1004－5422(2014)01－0029－04

2013－09－05.

彭 聰(1980—)，女，從事計算機軟件編程與算法研究.