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

彭 聰

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

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

彭 聰

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

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

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

0 引言

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

1 摩擦阻力計(jì)算方法

1.1 控制方程

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

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

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的計(jì)算采用二階精度的中心差分格式，其形式為:

2 算例及分析

2.1 層流平板附面層及摩擦阻力計(jì)算

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

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

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

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

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

2.2 NACA0012翼型附面層及摩擦阻力計(jì)算

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

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

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

圖4 不同攻角下的流線圖

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

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

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

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

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

3 結(jié)論

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

［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］葉建，林國(guó)華，鄒正平，等.低雷諾數(shù)下二維翼型繞流的流場(chǎng)特性分析［J］.航空動(dòng)力學(xué)報(bào)，2003，18(2):38 －45.

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

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

［6］閆超.計(jì)算流體力學(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—)，女，從事計(jì)算機(jī)軟件編程與算法研究.