周 斌，唐登海

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

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

周 斌，唐登海

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

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

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

0 引 言

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

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

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

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

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

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

1 設(shè)計計算方法

1.1 設(shè)計方法原理

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

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

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

（2）式就是常見的求解正問題的基于擾動速度勢的面元法方程，由于擾動勢可以表示為流體域總速度勢與來流速度勢的差，即有下式：

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

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

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

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

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

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

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

1.2 數(shù)值離散方程

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

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

2 驗證算例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖9 面元法（BEM）和RANS方法對翼型壓力分布的計算結(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ù)是全新的，除了采用面元法對其壓力分布進(jìn)行數(shù)值計算外，我們還采用RANS方法對母翼型和設(shè)計翼型的壓力分布形態(tài)分別進(jìn)行了數(shù)值計算。兩種方法的計算結(jié)果見圖9。從RANS和面元法的計算結(jié)果對比可以看出，在（0.0～0.9）弦長區(qū)域內(nèi)RANS方法和面元法計算的壓力分布主要特征一致，重合較好，在隨邊（0.9～1.0）弦長區(qū)域，壓力分布有一些區(qū)別，這個區(qū)別主要是因為面元法是基于勢流的計算方法，沒有考慮邊界層的影響。通過兩種方法的計算對比可以看出采用本方法可以較準(zhǔn)確地獲得滿足指定壓力分布的翼型。

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

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

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

3 結(jié)論與展望

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

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

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

從求解方法的原理上分析，本方法如果將（8）式中格林函數(shù)項取為三維格林函數(shù)，可以應(yīng)用于三維翼型的設(shè)計問題以及其他三維流動問題的設(shè)計，作者正在開展這方面的設(shè)計與應(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)師。