基于求解速度勢(shì)通量的指定壓力分布二維翼型設(shè)計(jì)方法

周 斌，唐登海

（中國(guó)船舶科學(xué)研究中心，江蘇無(wú)錫214082）

基于求解速度勢(shì)通量的指定壓力分布二維翼型設(shè)計(jì)方法

周 斌，唐登海

（中國(guó)船舶科學(xué)研究中心，江蘇無(wú)錫214082）

為了提升翼型的水動(dòng)力和空泡等性能，指定壓力分布的翼型剖面設(shè)計(jì)的方法多數(shù)集中在給定攻角下的翼型剖面的設(shè)計(jì)，該方法存在計(jì)算量較大，收斂性不理想，特別是推廣到三維問(wèn)題時(shí)，上述問(wèn)題尤為突出，限制了翼型設(shè)計(jì)的進(jìn)一步發(fā)展。文章以勢(shì)流理論面元方法為基礎(chǔ)，通過(guò)求解指定壓力分布條件下翼型表面的速度勢(shì)通量，獲得翼型表面形狀的修正量，并將修正量分解為攻角的變化以及剖面自身的變化兩部分，從而得到了翼型唯一的設(shè)計(jì)攻角和翼型剖面幾何（厚度分布、拱度分布）。文中采用上述方法對(duì)二維翼型問(wèn)題進(jìn)行了設(shè)計(jì)驗(yàn)證，表明該方法可以設(shè)計(jì)任意指定壓力分布的翼型剖面，理論上該方法可以用于全三維翼型的設(shè)計(jì)問(wèn)題。

翼型設(shè)計(jì)；指定壓力分布；勢(shì)流理論；新型翼剖面

0 引 言

翼型表面的壓力分布決定了翼型的流體動(dòng)力特征，對(duì)于水翼來(lái)說(shuō)，為了提升翼型的水動(dòng)力和空泡等性能，人們一直希望能夠設(shè)計(jì)出指定壓力分布的翼型幾何剖面。

指定任意壓力分布的翼型剖面設(shè)計(jì)方法已經(jīng)發(fā)展多年。Kennedy和Marsden（1978）[1]通過(guò)二維流函數(shù)求解法，在二維翼型表面采用渦層密度分布代替翼型剖面，獲得了指定攻角下，任意壓力分布的翼型剖面。Eppler和Shen（1979）[2]采用了保角變換的思想提出了在指定攻角下任意壓力分布的翼型剖面設(shè)計(jì)方法。

上述設(shè)計(jì)方法只適用于二維翼型剖面的設(shè)計(jì)。為了更精確地模擬物體表面，蘇玉民[3]采用面元法，在指定攻角條件下，通過(guò)求解控制點(diǎn)上幾何的細(xì)微擾動(dòng)對(duì)翼型剖面流動(dòng)的影響，獲取雅可比影響矩陣，然后對(duì)原始剖面進(jìn)行修改達(dá)到指定的壓力分布；李俊華[4]基于同樣的思想，采用B樣條來(lái)重新表達(dá)翼型剖面，通過(guò)改變控制點(diǎn)來(lái)調(diào)整翼型剖面的幾何，建立翼型剖面幾何與翼形表面壓力分布的關(guān)系，最終獲得滿足給定壓力分布的翼型剖面。上述關(guān)于翼型的設(shè)計(jì)方法都以面元法為基礎(chǔ)，可以在指定攻角條件下設(shè)計(jì)三維翼型的剖面。但是上述方法需要計(jì)算雅可比矩陣，計(jì)算量很大，假設(shè)三維翼型表面網(wǎng)格劃分N個(gè)，為獲得雅可比修正矩陣需要對(duì)N個(gè)面元分別進(jìn)行一次正問(wèn)題計(jì)算。如果需要迭代m次，則計(jì)算量正比于mO（N2），因此上述方法如果推廣至更復(fù)雜的全三維問(wèn)題，需要的計(jì)算量往往會(huì)很大。

Lee Chang-Sup（1994）[5]采用面元法求解指定攻角、指定壓力分布情況下，翼型表面的速度勢(shì)通量，然后修改原始翼型剖面獲得滿足設(shè)計(jì)要求的翼型剖面，并對(duì)二維問(wèn)題和三維的翼型設(shè)計(jì)問(wèn)題進(jìn)行計(jì)算，上述方法將翼型幾何的求解問(wèn)題轉(zhuǎn)換成求解速度勢(shì)通量，因此求解一次翼型設(shè)計(jì)問(wèn)題的工作量等于同樣求解一次正問(wèn)題，m次迭代的計(jì)算量為mO（N），較以往求解雅可比矩陣修正翼型剖面的方法減小了一個(gè)量級(jí)的計(jì)算量。

以上介紹的方法雖然采用不同的原理，但是都是在解決求解給定攻角、給定壓力分布的設(shè)計(jì)問(wèn)題。我們知道一個(gè)具有固定型值的翼型，在指定攻角下，其壓力分布是唯一的。反過(guò)來(lái)我們思考這樣的問(wèn)題，如果給定了一種壓力分布，那么相應(yīng)的翼型剖面和攻角是否是唯一的？

按照以往的給定攻角、給定壓力分布翼型設(shè)計(jì)方法，顯然認(rèn)為同樣的壓力分布可以有不同翼型和攻角的組合。本文以此為切入點(diǎn)，在Lee Chang-Sup（1994）[5]的基礎(chǔ)上，嘗試將翼型變化的修正量分解為翼型攻角的變化和翼型自身的修正量，這樣可以在翼型設(shè)計(jì)時(shí)獲得唯一的翼型的攻角和翼型剖面幾何（厚度分布、拱度分布）。數(shù)值計(jì)算表明本文方法較好地還原了給定壓力分布的翼型剖面和攻角，同時(shí)適用于任意指定壓力分布的翼型設(shè)計(jì)問(wèn)題，并且原理上適用于三維翼型的設(shè)計(jì)問(wèn)題，為面元法全面應(yīng)用于翼型設(shè)計(jì)問(wèn)題提供了新的思路。

1 設(shè)計(jì)計(jì)算方法

1.1 設(shè)計(jì)方法原理

由勢(shì)流理論可知，對(duì)于無(wú)界、無(wú)旋、定常、不可壓和無(wú)粘理想流體的有升力體繞流問(wèn)題，流域內(nèi)的擾動(dòng)速度勢(shì)可用下式表達(dá)：

其中：Sb表示物體表面，Sw表示有升力體的泄出渦表面即尾渦面，S∞表示無(wú)窮遠(yuǎn)處流體域界面，G表示格林函數(shù)，如圖1所示，無(wú)窮遠(yuǎn)處擾動(dòng)速度勢(shì)對(duì)物體表面的影響可視為一個(gè)常值速度勢(shì)φ0，于是（1）式可寫成

圖1 有升力體面元法流體計(jì)算域示意圖Fig.1 Definition sketch for lifting surface method

（2）式就是常見(jiàn)的求解正問(wèn)題的基于擾動(dòng)速度勢(shì)的面元法方程，由于擾動(dòng)勢(shì)可以表示為流體域總速度勢(shì)與來(lái)流速度勢(shì)的差，即有下式：

對(duì)于給定物體表面的壓力問(wèn)題，可以根據(jù)伯努利方程求解物體表面的速度分布，通過(guò)求解速度分布在物體表面的積分獲得物體表面的總的速度勢(shì)：

其中：vc表示物體表面沿主流方向（弦向）的速度分布。對(duì)于設(shè)計(jì)問(wèn)題，由于物體表面的壓力分布是指定的，因此其給定的Φsb也是指定的，將指定的Φsb代入（8）式時(shí)，可以求解方程中的未知項(xiàng)，其中的物理意義為指定壓力分布時(shí)，翼型表面的速度勢(shì)通量值。

對(duì)于翼型設(shè)計(jì)問(wèn)題，為了獲得調(diào)整后翼型剖面的攻角和翼型幾何參數(shù)，將翼型表面的調(diào)整量分為兩部分共同作用的結(jié)果。一部分是翼型攻角的作用；另一部是翼型幾何參數(shù)的變化。于是得到了如下翼型攻角α和翼型剖面自身的修正量Δt的計(jì)算公式：

其中：s_trail為翼型隨邊位置；s_lead為翼型導(dǎo)邊位置；s（x）為控制點(diǎn)弦向位置到隨邊的距離。

此外對(duì)于有升力翼型剖面的面元法求解問(wèn)題，需要在翼型隨邊尾緣給定庫(kù)塔條件，對(duì)于翼型剖面問(wèn)題可采用速度勢(shì)庫(kù)塔條件[6]，即：

此外（8）式中還存在未知量φ0，需要補(bǔ)充方程使得方程（8）獲得唯一解，根據(jù)調(diào)和函數(shù)性質(zhì)，補(bǔ)充方程的物理描述為在滿足給定壓力分布后，所有物體表面的速度勢(shì)通量與面元面積的乘積和為0，即有如下方程：

至此通過(guò)求解方程（8）、（10）、（11）、（12）式和（13）式，就可以獲得指定壓力分布下確定的翼型幾何參數(shù)和翼型攻角。

1.2 數(shù)值離散方程

（8）式是對(duì)一般問(wèn)題的描述，對(duì)于二維翼型剖面的設(shè)計(jì)問(wèn)題，可以沿翼型表面劃分Np個(gè)單元，為保證計(jì)算精度以及獲得光滑的翼型剖面，采用余弦劃分形式，以保證網(wǎng)格在導(dǎo)邊和隨邊附近加密。此時(shí)對(duì)于每一個(gè)網(wǎng)格單元，（8）式可以離散為

求解線性方程組（15）便可以求解每個(gè)翼型網(wǎng)格上的通量變化量，爾后根據(jù)（10）式、（11）式對(duì)翼型剖面進(jìn)行修改迭代，最終可得到滿足指定壓力分布的翼型剖面。

2 驗(yàn)證算例

2.1 已知翼型壓力分布的翼型參數(shù)復(fù)原

為了驗(yàn)證本方法求解的翼型剖面及攻角是唯一的，對(duì)已知翼型剖面幾何的二維翼型剖面設(shè)計(jì)問(wèn)題進(jìn)行了算例驗(yàn)證。

首先采用二維擾動(dòng)速度勢(shì)面元法對(duì)NACA0025翼型在2°攻角條件下的壓力分布進(jìn)行了計(jì)算，然后將該壓力分布賦值給NACA0010，0°攻角的翼型，開(kāi)始數(shù)值迭代設(shè)計(jì)。為了獲得精細(xì)的物體表面幾何特征，對(duì)翼型剖面弦向劃分了100個(gè)直線單元。

NACA0010，0°攻角的翼型剖面與迭代過(guò)程剖面的壓力分布見(jiàn)圖2，NACA0010，0°攻角的翼型幾何與NACA0025，2°攻角翼型幾何及迭代過(guò)程剖面幾何見(jiàn)圖3；從圖2、圖3可以看出，采用本方法經(jīng)過(guò)3次迭代后的翼型幾何與目標(biāo)翼型幾何便具有良好的重合度。數(shù)值測(cè)試共迭代了10次，獲得的翼型剖面攻角收斂在2.07°附近與設(shè)計(jì)目標(biāo)值NACA0025，2°的給定值非常接近。

圖2 NACA0010，0°攻角翼型至NACA0025，2°攻角翼型的壓力分布迭代過(guò)程Fig.2 Iterative process of pressure distribution from NACA0010 airfoil at 0 deg to NACA0025 at 2 deg

圖3 NACA0010，0°攻角翼型至NACA0025，2°攻角翼型的幾何迭代過(guò)程Fig.3 Iterative process of section geometry from NACA0010 airfoil at 0 deg to NACA0025 at 2 deg

通過(guò)上述算例的驗(yàn)算，較好地復(fù)原了已知翼型壓力分布情況下原翼型的幾何參數(shù)，也說(shuō)明了在指定壓力分布情況下，可以獲得唯一的翼型剖面和相應(yīng)的攻角。

2.2 指定壓力分布的翼型剖面設(shè)計(jì)

為了獲得翼型表面更好的流動(dòng)形態(tài)，在翼型設(shè)計(jì)中希望得到指定壓力分布形式的翼型參數(shù)。我們以NACA0010翼型在攻角3°時(shí)的翼型壓力分布為基礎(chǔ)，指定了新的翼型壓力分布形態(tài)，見(jiàn)圖4。其中“NACA0010，α=3°”表示NACA0010翼型在攻角3°時(shí)翼型表面壓力分布；“Target Cp”表示指定的目標(biāo)壓力分布形態(tài)（抑制翼剖面導(dǎo)邊的負(fù)壓峰值）。以此為目標(biāo)壓力，應(yīng)用本文介紹的方法進(jìn)行翼型的設(shè)計(jì)。

圖4 指定的壓力分布與母翼型的壓力分布Fig.4 The specified flat rooftop pressure distribution and the original pressure distribution of NACA0010 at 3 deg

圖5 翼型參數(shù)隨著迭代過(guò)程的變化趨勢(shì)Fig.5 Iterative process of the airfoil parameters

首先以3°攻角NACA0010翼型為母翼型，對(duì)翼型剖面弦向劃分了100個(gè)直線單元，然后迭代計(jì)算獲取滿足指定壓力分布要求的翼型剖面。對(duì)迭代過(guò)程中獲取的翼型剖面的攻角、最大拱弧和最大厚度進(jìn)行提取，可以獲得翼型參數(shù)隨著迭代過(guò)程的變化趨勢(shì)，見(jiàn)圖5。其中α表示攻角，fmax表示翼型剖面的最大拱度，tmax表示翼型剖面的最大厚度。從翼型參數(shù)的變化可以看出，在進(jìn)行3次迭代后翼型參數(shù)便趨于收斂。迭代過(guò)程產(chǎn)生的翼型幾何見(jiàn)圖6，翼型導(dǎo)邊局部幾何見(jiàn)圖7，相應(yīng)的翼型表面壓力分布變化見(jiàn)圖8。經(jīng)過(guò)10次迭代后滿足指定壓力分布要求的翼型攻角為2.54°，最大拱度比為0.011，最大厚度比為0.116。

圖6 迭代過(guò)程產(chǎn)生的翼型幾何Fig.6 Iterative process of section geometry

圖7 迭代過(guò)程翼型導(dǎo)邊局部幾何Fig.7 Close view of the leading edge geometry

圖8 迭代過(guò)程中翼型表面壓力分布變化Fig.8 Iterative process of pressure distributions

圖9 面元法（BEM）和RANS方法對(duì)翼型壓力分布的計(jì)算結(jié)果Fig.9 Comparisons of pressure distributions obtained from BEM and RANS methods for the design airfoil

圖10 指定的“鋸齒形”的壓力分布Fig.10 The specified‘zigzag'pressure distribution compared with baseline NACA0010 airfoil at 3 deg

圖11 滿足“鋸齒形”的壓力分布的翼型幾何Fig.11 Designed section geometry with‘zigzag'pressure distribution compared with baseline NACA0010 airfoil

由于上述獲得的滿足指定壓力分布翼型幾何參數(shù)是全新的，除了采用面元法對(duì)其壓力分布進(jìn)行數(shù)值計(jì)算外，我們還采用RANS方法對(duì)母翼型和設(shè)計(jì)翼型的壓力分布形態(tài)分別進(jìn)行了數(shù)值計(jì)算。兩種方法的計(jì)算結(jié)果見(jiàn)圖9。從RANS和面元法的計(jì)算結(jié)果對(duì)比可以看出，在（0.0～0.9）弦長(zhǎng)區(qū)域內(nèi)RANS方法和面元法計(jì)算的壓力分布主要特征一致，重合較好，在隨邊（0.9～1.0）弦長(zhǎng)區(qū)域，壓力分布有一些區(qū)別，這個(gè)區(qū)別主要是因?yàn)槊嬖ㄊ腔趧?shì)流的計(jì)算方法，沒(méi)有考慮邊界層的影響。通過(guò)兩種方法的計(jì)算對(duì)比可以看出采用本方法可以較準(zhǔn)確地獲得滿足指定壓力分布的翼型。

在完成上述指定壓力分布翼型的剖面設(shè)計(jì)后，為了考察本方法的穩(wěn)定性和適用性，我們還對(duì)更加“特殊”的指定壓力分布進(jìn)行了翼型設(shè)計(jì)。“特殊”的指定壓力分布見(jiàn)圖10，從其壓力分布可以看出，給定的壓力分布在翼型上表面存在 “鋸齒形”的壓力分布區(qū)域。仍以NACA0010翼型為母翼型，迭代10次后獲得的翼型剖面見(jiàn)圖11。仍采用RANS方法和面元法分別對(duì)設(shè)計(jì)得到的翼型的壓力分布形態(tài)進(jìn)行了數(shù)值計(jì)算，兩種方法的計(jì)算結(jié)果見(jiàn)圖12。從兩種計(jì)算得到的壓力分布可以觀察到所計(jì)算的翼型上表面也存在“鋸齒形”的分布形態(tài)，說(shuō)明本方法具有較好的適應(yīng)性。

通過(guò)上述算例的驗(yàn)算，可以看出本方法適用于任意指定壓力分布設(shè)計(jì)問(wèn)題的求解。

圖12 面元法（BEM）和RANS方法對(duì)設(shè)計(jì)翼型壓力分布形態(tài)的計(jì)算結(jié)果Fig.12 Comparisons of pressure distributions obtained from BEM and RANS methods for the design airfoil

3 結(jié)論與展望

本文以勢(shì)流理論、總速度勢(shì)面元方法為基礎(chǔ)，通過(guò)求解指定壓力分布設(shè)計(jì)問(wèn)題物體表面的速度勢(shì)通量，將翼型剖面的修正量分為因攻角變化引起的改變和翼型剖面自身改變兩個(gè)部分，由此在指定壓力分布的條件下，可以獲得唯一的翼型幾何參數(shù)和攻角，為面元法應(yīng)用于翼型設(shè)計(jì)問(wèn)題提供了新的思路。

數(shù)值驗(yàn)算結(jié)果表明，本方法可以較好地復(fù)原給定壓力分布的翼型剖面，同時(shí)本方法也適用于任意指定壓力分布翼型設(shè)計(jì)問(wèn)題求解。

從求解方法的計(jì)算量分析，本方法從基于擾動(dòng)速度勢(shì)的面元法方程推導(dǎo)得出，因此在翼型設(shè)計(jì)時(shí)，不需要計(jì)算壓力分布與幾何變化之間的雅可比矩陣，計(jì)算量與面元法求解正問(wèn)題處于一個(gè)量級(jí)。同時(shí)數(shù)值算例表明，本方法收斂性很好。

從求解方法的原理上分析，本方法如果將（8）式中格林函數(shù)項(xiàng)取為三維格林函數(shù)，可以應(yīng)用于三維翼型的設(shè)計(jì)問(wèn)題以及其他三維流動(dòng)問(wèn)題的設(shè)計(jì)，作者正在開(kāi)展這方面的設(shè)計(jì)與應(yīng)用研究。

[1]Kennedy J L,Marsden D J.A potential flow design method for multicomponent airfoil sections[J].Journal of Aircraft,1978, 15(1):47-52.

[2]Eppler R,Shen Y T.Wing sections for hydrofoils-part 1:Symmetrical profiles[J].Journal of Ship Research,1979,23(9): 209-217.

[3]Su Yumin,Ikehata M,Kai H.A numerical method for designing three-dimentional wing based on surface panel method[J]. Journal of the Society of Naval Architects of Japan,1997,182(1):39-46.

[4]Li Junhua,Tang Denghai,Dong Shitang.Propelelr design by prescribed pressure distribution[J].Journal of Ship Mechanics, 2010,14(2):10-19.

[5]Lee Chang-Sup,Kim Young-Gi,et al.A surface panel method for design of hydrofoils[J].Journal of Ship Research,1994, 38(9):175-181.

[6]Morino L.Subsonic potential aerodynamics for complex configurations-A general theroy[J].AIAA Journal,1974,12(2): 191-197.

A design method for 2-D hydrofoil by the specified pressure distribution based on velocity potential flux of the surface

ZHOU Bin,TANG Deng-hai
(China Ship Scientific Research Center,Wuxi 214082,China)

In order to improve the hydrodynamic and cavitation performance of the hydrofoil,some hydrofoil section design methods with specified pressure distribution are developed.However the section is designed at a given angle of attack in the most of such methods,which may result in large computation cost and difficulties in convergence especially when the methods are extended to 3-D problems.In the paper,a new 2-D hydrofoil section design method is proposed based on surface panel method of the potential theory.The velocity potential flux and geometric corrections of the blade surface can be obtained by solving total velocity potential equation with specified pressure distribution on the blade surface.The geometric correction could be decomposed into two parts.One is owing to the change of angle of attack,and the other is the change of profile itself.In this way,the unique design angle of attack and section geometry for specified pressure distribution can be attained.The numerical results show that the method is robust and valid.

hydrofoil profile design;specified pressure distribution;potential theory;new blade section

U661.1

：Adoi:10.3969/j.issn.1007-7294.2016.04.003

1007-7294（2016）04-0403-07

2015-07-07

周 斌（1986-），男，工程師，E-mail:htrmax@163.com；唐登海（1965-），男，研究員，博士生導(dǎo)師。