1 引言
關於翼型結構設計和升力的理論研究十分久遠[1,2]，Abbot和Doenhoff[3]系統地分析了翼型的結構與攻角對機翼的升力和升力因數的影響。Critzos等[4]對NACA 0012翼型攻角從0變化到108度過程中的空氣動力學特性進行了研究。Lissaman[5]忽略可壓縮性的影響，對雷諾數處於104到106之間的二維翼型的流體力學效能進行了綜述，指出過渡、湍流、分離和回附是基礎的經典流體力學問題。Storms和Jang[6]採用古奈擾流板和渦發生器對翼型升力影響因素進行了試驗研究；實驗測量了NACA 4412 翼型的表面壓力分佈和尾流渦形狀；發現使用古奈擾流板可以提高翼型最大升力同時能夠減小尾翼脫流渦的強度。
Somers[7]設計了層流翼型S809的水平軸風力渦輪機裝置並理論和實驗驗證了此翼型的效能，解決了受限的升力和輪廓阻力問題。Li 等[8]對NACA0012翼型進行了表面壓力分佈和尾流分佈的實驗測量，從而確定不同配置下的升力、阻力和俯仰力矩係數；加入古奈擾流板使最大升力因數從1.37增加到1.74；還對邊界層剖面採用總壓探頭的前傾角吸力側90%弦位置進行了測量，發現有效古奈擾流板高度約為弦長的2%時，與無古奈擾流板的NACA0012翼型相比，提供了研究配置中最高的升阻比。
Timmer 和Van[9]基於代爾夫特理工大學(DUT)專用翼型的風力發電機的設計和風洞試驗結果；DUT型翼型的最初設計由XFOIL程式碼完成，還介紹了無古奈擾流板、後緣楔片、旋渦發生器(vg)和脫扣線對各種DUT型機翼特性的影響；最後給出了翼型前緣厚度與前緣分離的攻角的關係。Lanzafame和Messina[10]建立並改進了基於葉元動量理論的流體動力風力機設計數學模型，在設計工況下和非設計工況下，對整個風速範圍進行了模擬；重點分析了切向誘導因子和表示升力和阻力系數的模型。Gie等[11]設計了低壓渦輪以獲得效率和重量之間的最佳效果，繼而實現最低的發動機燃料損耗；分析了翼型升力水平對損失和效率的主要影響及機理，為確定最佳翼型升力選擇提供了依據。
Ormsbee和Chen[12]設計了不可壓縮流體在大雷諾數下的最大升力因數意義下的最優翼型，確定了最佳壓力分佈是雷諾數和尾緣速度的函式並透過一種直接迭代方法獲取能夠產生最優壓力分佈的翼型形狀。Bénard等[13]研究了等離子體作動器控制的軸對稱翼型的升力和阻力效能；在NACA 0015機翼模型的前緣安裝了一個介質阻擋放電(DBD)裝置，採用時間平均法研究了穩態和非穩態作用對升力和阻力系數的影響。Suresh等[14]考慮了動葉失速效應的辨識，採用複發性神經網路(RNN)方法，預測了高攻角時的升力因數(CZ)。
Somers[15]研製了翼型矩陣和最大升力對最大升力的靈敏度對前緣粗糙度的影響，在賓夕法尼亞州立大學的實驗中低速、低湍流風洞對翼型進行了設計、理論分析和驗證；研究發現粗糙度對最大升力的影響隨機翼厚度的增加而增大，隨最大升力的增加而略有下降。Whitcomb[16]基於基本的空氣動力學現象對美國宇航局超臨界翼型的大幅度減少上表面的中弦區曲率隨著增大而增大在尾緣附近彎曲現象進行了討論，結果發現NASA超臨界翼型的馬赫數隨阻力的增大而增大比同等的NACA 6系列機翼高0.1倍。Tetervin [17]採用二維低雷諾數的風洞對NACA 66,2-216,a=0.6翼型的邊界層效能進行了研究。Shur[18]提出了分離渦流模擬(DES)在三維意義上的第一個真正的應用, 進行了高攻角的翼型的模擬。Miao[19]研究了不同雷諾數組合減振頻率下弦向彎曲幅值對撲動翼型非定常氣動特性的影響，考慮八種不同的撓曲振幅來研究了柔性振幅對撲翼翼型氣動效能的影響。
2 單個NACA翼型的計算2.1 基於COMSOL軟體的模擬計算
2.1.1 翼型引數的設定氣體流動模型的設定
NACA翼型是一種特殊的形狀，可以由幾個特定引數進行描述。在本研究中，關注其翼長，最大拱高，最大拱高出現的位置以及翼的厚度。其定義如下圖所示。
圖2-1 翼型引數模型
2.1.2 氣體流動模型的設定
我們認為空氣為可壓縮氣體，有粘性。取空氣的Re數，即，當Re>4000時，流動為湍流，採用Spalart-Allmaras湍流模型進行計算，當Re<4000時，流動為層流，採用層流模型。在模擬時取的速度範圍在1m/s-9m/s之間，經計算均滿足湍流條件，因而著重介紹湍流模型的設定方法。
流動模型應用域為氣體流動域。定義初始值時氣體流動速度Ux(0)=u_in為可輸入的變數，並且同時定義Uy(0)=0。由於氣體粘性在翼型表面產生的渦量擴散能力可以用無衰減湍流運動學粘度來代替，其量綱為，可呼叫材料庫中的引數得到。模型邊界條件中，入口邊界條件為指定速度場，即入口速度保持與U(0)一致，出口邊界條件為壓力為0，使得氣體自由流出。
2.1.3 模擬引數的輸入及計算過程
首先利用軟體自帶的計算Re數的功能，透過輸入空氣流速、翼長、空氣密度以及空氣動力粘度進行Re數計算，如下圖所示。
圖2-2 Re數計算
根據Re數不同，來選取不同的氣體流動計算模型。而後調整翼型的最大拱高，最大拱高出現的位置以及翼的厚度。改變攻角的大小，進行迭代求解，如下圖所示。
圖2-3 模擬介面
迭代求解後即可得到速度分佈、壓力分佈以及升力大小和升力因數。設定翼長為1m。查表可知，20℃下，空氣密度為ρ=1.2kg/m³，動力粘度μ=Pa·s。設定空氣流速變化值為1、3、5、7、9m/s，攻角為0°，±2°，±4°，±6°，±8°。考慮了三種翼型，其引數如下表所示。
表2-1 三種翼型引數
在確定流場的外部條件之後，使用COMSOL軟體，輸入Case 1的引數條件和流場的流動條件，得到其升力因數隨流場流速和攻角變化的變化情況。計算結果如表2-2所示。將計算結果整理成三維圖如圖2-4所示。
表2-2 Case 1的模擬結果
圖2-4 升力因數隨攻角及速度的變化趨勢
透過表中資料，本文可以得到，在以上的輸入條件下，當攻角為正時，升力因數隨著流體流速的增大而上升，當攻角為負時，升力因數隨著流體流速的增大而下降。且升力因數的變化幅度小，以1m/s時的升力因數為基準，最大的升力因數相對變化量為12.4%。理論上，翼型的升力因數和流速無關，但是由於實際流體粘性、可壓縮性和渦的存在，此翼型的升力因數在流體流速上升的時候，升力因數也是處於上升的趨勢。在流體流速不變的情況下，升力因數隨著攻角的增大而增大，且相對攻角為0°有對稱關係。在攻角為負時，此翼型處於向下的運動趨勢，相應的翼型上端的平均流速小於下端流速，應用伯努利原理，上端壓力大於下端壓力，升力因數為負；在攻角為正時，相應的升力因數為正。同時，在0°時，升力因數很小，幾乎為0，這是因為在Case 1中，最大彎度為0，因此此翼型為對稱翼型，相應的流場也是對稱的，理論上在攻角為0°時，升力因數為0。
以流場流速U=3m/s，攻角為8°為例，分析翼型附近的壓力場和流場，在軟體中輸出翼型附近標準化壓力場，如圖2-5所示，輸出標準化等壓線分佈，如圖2-6所示，輸出單位化速度雲圖，如圖2-6所示。
圖2-5 翼型附近標準化壓力場（以U=3m/s，攻角為8°為例）
圖2-6 標準化等壓線分佈
圖2-7 單位化速度雲圖
在圖2-5中，翼型最前端有一個高壓區，這是因為流體接觸到翼型時，速度急劇下降，根據伯努利原理，壓力也急劇上升，在翼型的上沿，有較大一塊低壓區，壓力下降很大，這是翼型升力的來源，在翼型的下沿，亦有壓力下降，但是壓力下降較小。由於攻角的原因，翼型上沿很快便形成了脫體繞流，渦結構造成了很大一片低壓區。然而在翼型的下沿，固體邊界能夠和流體的流動良好貼合，壓力下降較小。翼型上下沿的壓力合成了翼型的升力。在圖2-6中，等壓線的分佈也和上文的分析吻合。在圖2-7中，翼型上沿的流體流速顯然大於來流的流體流速，在翼型下沿，流體的流速略大於來流流速，但是分佈顯然比上沿均勻。
再在COMSOL軟體中輸入Case 2的翼型引數，流場的引數輸入與Case 1相同。得到Case 2的模擬計算結果，如表2-3所示。再整理資料，得到升力因數隨攻角及速度的變化趨勢，如圖2-8所示。
表2-3 Case 2的模擬結果
圖2-8 升力因數隨攻角及速度的變化趨勢
透過表2-3中資料，在Case 2的輸入條件下，在流場流速不變的情況下，翼型的升力因數隨著攻角的增大而增大，而且有一個類似線性的增長過程。在-4°升力因數為負，在-2°時，升力因數為正，可以推測得到，在-4°~-2°之間某一個攻角，翼型的升力為0。考慮到Case 1中，升力因數為0時，攻角在0°附近，因此，Case 2的翼型即使在負攻角，依然可能具有上升的能力。Case 2的翼型相比於Case 1是非對稱的翼型，其最大彎度為0.05，能夠更能貼合流場。在攻角不變的情況下分析，當流速為3、5、7、9m/s時，其升力係數變化很小，可以忽略其變化。但是流速1m/s變化到3m/s時，升力係數有一個相對較大的變化，這是因為流體的雷諾數相對較小，粘性力的作用相對明顯導致的。
以流場流速U=3m/s，攻角為8°為例，分析翼型附近的壓力場和流場，在軟體中輸出翼型附近標準化壓力場，如圖2-9所示，輸出標準化等壓線分佈，如圖2-10所示，輸出單位化速度雲圖，如圖2-11所示。
圖2-9 翼型附近標準化壓力場（以U=3m/s，攻角為8°為例）
圖2-10 標準化等壓線分佈
圖2-11 單位化速度雲圖
綜合比較圖2-5、2-6、2-9、2-10，可以分析得到，Case 2翼型附近的流場壓力和速度分佈同Case 1類似，不同的是，Case 2翼型上沿的低壓區壓力更低，特別是在翼型上沿末端，翼型下沿壓力相對高，而翼型前段兩個駐點的位置相近，變化不大。這總體導致了Case 2翼型的升力因數更大。綜合比較圖2-11和圖2-7，可以得到Case 2翼型相對於Case 1翼型，其上沿的速度紅色區域更大更深，下沿的速度更接近來流速度，這說明了Case 2翼型能夠在翼型上沿產生更劇烈的渦運動和低壓區，而在下沿，Case 2翼型能夠良好更好的貼近流體運動，這兩部分導致了Case 2翼型的升力因數更大。
在COMSOL軟體中輸入Case 3的翼型引數，流場的引數輸入與Case 1相同。得到Case 3的模擬計算結果，如表2-4所示。再整理資料，得到升力因數隨攻角及速度的變化趨勢，如圖2-12所示。
表2-4 Case 3的模擬結果
圖2-12 升力因數隨攻角及速度的變化趨勢
透過表2-4中資料，在Case 3的輸入條件下，在流場流速不變的情況下，翼型的升力因數隨著攻角的增大而增大，而且有一個類似線性的增長過程，這和Case 1 和Case 2都是類似的。類似的，可以推測得到，在-4°~-6°之間某一個攻角，翼型的升力為0。考慮到Case 1中，升力因數為0時，攻角在0°附近，在Case 2中，升力因數為0時，攻角在-3°附近，因此，Case 3的翼型相當於Case 1和Case 2，在負攻角的情況下，具有更好的升力特性。Case 3的翼型相對於Case 2的翼型，其最大厚度減小，相對於Case 1，最大彎度增加，最大厚度減小。總體而言，Case 3相對Case 1和Case 2在每種列出的條件下，均有更大的升力因數，因此以升力因數為指標分析，Case 3翼型更佳。在攻角不變的情況下分析，類似於Case 2，當流速為3、5、7、9m/s時，其升力係數變化很小，可以忽略其變化。但是流速1m/s變化到3m/s時，升力係數有一個相對較大的變化，這是因為流體的雷諾數相對較小，粘性力的作用相對明顯導致的。
以流場流速U=3m/s，攻角為8°為例，分析翼型附近的壓力場和流場，在軟體中輸出翼型附近標準化壓力場，如圖2-13所示；輸出標準化等壓線分佈，如圖2-14所示；輸出單位化速度雲圖，如圖2-15所示。
圖2-13 翼型附近標準化壓力場（以U=3m/s，攻角為8°為例）
圖2-14 標準化等壓線分佈
圖2-15 單位化速度雲圖
綜合比較圖2-9、圖2-10、圖2-13、圖2-14，可以分析得到，Case 3翼型附近的流場壓力和速度分佈同Case 2類似，Case 3翼型上沿的壓力更低，Case 3翼型下沿的壓力更好，其壓力甚至略高於設定的氣壓，對於駐點，靠上的駐點位置幾乎不變，然而靠下的駐點則消失，這說明在翼型下沿，邊界層不產生脫體流動，流體良好的貼合翼型的邊緣流動，相應的壓力也較高。因此，Case 3翼型上下沿的壓差更大，升力因數也更高。綜合比較圖2-11和圖2-15，Case 3翼型上沿的流體運動更加劇烈。
在較大攻角、Re=66666.667的條件下的不同攻角對應的速度場和壓力場如圖2-16至圖2-20所示。在攻角較大的條件下，可以明顯地看到翼型後側出現了速度極小的區域，這一現象即為失速，此時翼型後會產生大量的渦結構，實際速度場與模擬結果是不匹配的，但整體的速度分佈情況是一致的，因而可以用模擬結果做速度參考。
攻角為15°的單位化速度雲圖
攻角為15°的標準化等壓線分佈
攻角為-15°的單位化速度雲圖
攻角為-15°的標準化等壓線分佈
攻角為30°的標準化單位化速度雲圖
3 結論由模擬結果可知，在改變鴨翼的攻角後，後方大翼型收到的升力發生了明顯的改變，因而增加鴨翼對提高升力因數是有幫助的，且鴨翼攻角越小，得到的升力因數越大。至於在攻角較大時升力因數的再次升高，是由於失速導致的，此時的升力因數計算已經不準確了。
參考文獻[1] Jacobs E N, Sherman A. Airfoil section characteristics as affected by variations of the Reynolds number[J]. 1937.
[2] Jameson A. Iterative solution of transonic flows over airfoils and wings, including flows at Mach 1[J]. Communications on pure and applied mathematics, 1974, 27(3): 283-309.
[3] Abbott I H, Von Doenhoff A E. Theory of wing sections, including a summary of airfoil data[M]. Courier Corporation, 1959.
[4] Critzos C C, Heyson H H, Boswinkle R W. Aerodynamic characteristics of NACA 0012 airfoil section at angles of attack from 0 to 180 degrees[M]. National Advisory Committee for Aeronautics, 1955.
[5] Lissaman P B S. Low-Reynolds-number airfoils[J]. Annual review of fluid mechanics, 1983, 15(1): 223-239.
[6] Storms B L, Jang C S. Lift enhancement of an airfoil using a Gurney flap and vortex generators[J]. Journal of Aircraft, 1994, 31(3): 542-547.
[7] Somers D M. Design and experimental results for the S809 airfoil[R]. National Renewable Energy Lab., Golden, CO (United States), 1997.
[8] Li Y, Wang J, Zhang P. Effects of Gurney flaps on a NACA0012 airfoil[J]. Flow, Turbulence and Combustion, 2002, 68(1): 27.
[9] Timmer W A, Van Rooij R. Summary of the Delft University wind turbine dedicated airfoils[J]. Journal of solar energy engineering, 2003, 125(4): 488-496.
[10] Lanzafame R, Messina M. Fluid dynamics wind turbine design: Critical analysis, optimization and application of BEM theory[J]. Renewable energy, 2007, 32(14): 2291-2305.
[11] Gier J, Franke M, Hübner N, et al. Designing LP turbines for optimized airfoil lift[C]//ASME Turbo Expo 2008: Power for Land, Sea, and Air. American Society of Mechanical Engineers, 2008: 1345-1358.
[12] Ormsbee A I, Chen A W. Multiple element airfoils optimized for maximum lift coefficient[J]. AIAA Journal, 1972, 10(12): 1620-1624.
[13] Bénard N, Jolibois J, Moreau E. Lift and drag performances of an axisymmetric airfoil controlled by plasma actuator[J]. Journal of Electrostatics, 2009, 67(2-3): 133-139.
[14] Suresh S, Omkar S N, Mani V, et al. Lift coefficient prediction at high angle of attack using recurrent neural network[J]. Aerospace Science and Technology, 2003, 7(8): 595-602.
[15] Somers D M. Effects of airfoil thickness and maximum lift coefficient on roughness sensitivity[J]. National Renewable Energy Laboratory, USA, 2005.
[16] Whitcomb R T. Review of NASA supercritical airfoils[J]. ICAS paper, 1974 (74-10): 25-30.
[17] Tetervin N. Investigation of the variation of lift coefficient with Reynolds number at a moderate angle of attack on a low-drag airfoil[J]. 1942.
[18] Shur M, Spalart P R, Strelets M, et al. Detached-eddy simulation of an airfoil at high angle of attack[M]//Engineering Turbulence Modelling and Experiments 4. 1999: 669-678.
[19] Miao J M, Ho M H. Effect of flexure on aerodynamic propulsive efficiency of flapping flexible airfoil[J]. Journal of Fluids and Structures, 2006, 22(3): 401-419.