基于流固耦合的柔性翼型動態(tài)失速特性研究

戴 昱， 陳德龍， 信志強

（河海大學力學與材料學院，南京 211100）

隨著時代的進步，微型飛行器的機翼、飛機的機翼以及風機葉片都在飛速發(fā)展，這對能源的利用和信息的收集都有著巨大的幫助.因此對翼型進行深入研究和改良創(chuàng)新是十分重要的.

如何提高翼型的氣動性能是研究的核心問題.目前已有眾多研究采用主動控制的方法使翼型發(fā)生變形來提升翼型的氣動性能.張敬斌等［1］通過主動控制技術(shù)研究翼型尾緣部分的變形，發(fā)現(xiàn)翼型尾緣部分發(fā)生適當變形能夠有效地提高翼型的升力系數(shù)和升阻比.而Zhuang等［2］的研究發(fā)現(xiàn)，在一定范圍內(nèi)，隨著翼型尾緣部分主動變形的范圍增大，翼型的氣動性能也會隨之提高.對于俯仰振蕩翼型，張俊偉等［3］研究的結(jié)果表明，在翼型上仰過程中翼型主動變形量的增大可以抑制渦的分離.Kan等［4］則對翼型后緣的主動變形進行了進一步的研究，發(fā)現(xiàn)變形的角度如果發(fā)生周期性的變化能夠增大升力系數(shù)并且可以延遲失速.尾緣主動變形的翼型和普通的襟翼相比存在一些差異，在這之前，普通襟翼的研究已有一些，例如賈亞雷等［5-6］的研究，而由于無縫光滑的特征，翼型尾緣的變形具有更優(yōu)越的氣動性能［7］.通過這些研究可以看出，對翼型施加主動控制使翼型變形從而提高氣動性能的方法具有較為優(yōu)越的效果.雖然主動變形控制的方法對氣動性能的提升已經(jīng)呈現(xiàn)出較為優(yōu)越的效果，但它需要主動控制系統(tǒng)，這往往導致其結(jié)構(gòu)較為復雜，被動控制則相對簡單.被動控制的方式主要是對翼型采用柔性的結(jié)構(gòu)，使其可以產(chǎn)生自適應變形從而提升翼型的氣動性能.自適應變形的潛能是巨大的［8］.Macphee和Beyene［9］采用被動控制，利用流固耦合的方法對柔性翼型的氣動性能展開了研究.研究表明，在垂直軸風機中柔性翼型的功率系數(shù)要大于剛性翼型，即使是低彈性模量的柔性翼型，相比于剛性翼型依然可以提高翼型的氣動特性［10］.柔性翼型在低于風機的額定風速下所表現(xiàn)出來的性能要優(yōu)于剛性翼型［11］.前人對柔性翼型的研究中還未考慮到柔性翼型的動態(tài)失速特性，而動態(tài)失速現(xiàn)象在風力機和各種飛行器中又是普遍存在的［12-16］.因此對于柔性翼型在動態(tài)失速現(xiàn)象中所表現(xiàn)出來的氣動特性進行研究是十分有必要的.

本文的研究是在前人對柔性翼型研究的基礎(chǔ)上，進一步考慮柔性翼型在動態(tài)失速現(xiàn)象中是否還具有一定的價值.因此本文基于流固耦合的方法，研究柔性翼型在不同減縮頻率下的動態(tài)失速特性.通過對比剛?cè)嵝砸硇偷臍鈩犹匦浴⒘鲌鎏攸c以及分析柔性翼型的結(jié)構(gòu)響應，得到柔性翼型的動態(tài)失速特性及其優(yōu)勢.

1 問題描述

本文采用的柔性翼型的模型如圖1所示，以NACA0012為基礎(chǔ)翼型，弦長c=1 m，并在距離翼型前端0.25 m處設置了一個0.02 m×0.06 m的矩形孔洞作為旋轉(zhuǎn)軸，整個翼型除矩形孔洞外，其余部分均為柔性材料.

圖1 柔性翼型模型Fig.1 Flexible airfoil model

來流介質(zhì)為空氣，密度取為1.225 kg/m3，來流速度為4.5 m/s，動力黏度取1.837 5×10-5Pa·s，雷諾數(shù)為30 萬.參考前人對于柔性翼型材料參數(shù)的選擇，本文中的柔性材料密度為100 kg/m3，泊松比為0.37，這種材料可用聚氨酯泡沫實現(xiàn)，具有低密度和廣泛的彈性模量的特征［9］.在本次動態(tài)失速研究中，通過翼型繞矩形旋轉(zhuǎn)軸做俯仰運動來模擬攻角的非定常變化，攻角的變化運動規(guī)律如式（1）所示：

其中：α0為初始攻角；αt為翼型俯仰過程中最大攻角；T為俯仰運動周期；n為俯仰運動周期的個數(shù).前3 s 翼型固定不動以獲取穩(wěn)定流場，3 s后做角速度大小不變的俯仰運動.

減縮頻率是翼型動態(tài)失速的重要影響因素［17］，其表達式為k=，其中ω表示振蕩角速度，c表示弦長，u為來流速度.本文針對剛?cè)醿煞N翼型分別設置了三種減縮頻率，具體工況設置如表1所示.

表1 工況設置Tab.1 Parameters of various working conditions

2 數(shù)值方法及驗證

2.1 控制方程

2.1.1 流體控制方程 對于不可壓縮流體而言，流體密度為常數(shù)，質(zhì)量守恒方程可寫為

式中：i為下標，表示坐標軸的方向；u代表流體速度.

對于湍流情況，雷諾時均方程可表示為

式中：j與前者i同為下標；uˉ為時均流速；u′為脈動流速；t代表時間；f表示質(zhì)量力；ρ為流體密度；p代表壓強；μ為流體的動力黏度.而雷諾時均方程相比于N-S方程多出了一項雷諾應力，該項為對稱的二階張量，與速度梯度有關(guān).這導致了雷諾時均方程無法封閉，需要建立湍流模型.本文將采用SST K-omega湍流模型［9，18］，該模型結(jié)合了k-omega模型和k-epsilon模型，能夠精確地捕捉近壁面的流動分離.

2.1.2 固體控制方程 柔性翼型采用的是齊次和各向同性結(jié)構(gòu)，可使用以下方程［19］：

其中：u為位移場；ρs為固體的密度；P為皮奧拉-基爾霍夫第一應力張量；S表示皮奧拉-基爾霍夫第二應力張量；F代表變形梯度；E代表格林應變張量；I表示特征張量；μ和λ為拉梅參數(shù).

2.1.3 流固耦合 在耦合面上，流體模型和固體模型須滿足運動學條件和動力學條件進行耦合.運動學條件是指耦合面上的速度和位移必須連續(xù)，可用以下表達式進行描述：

式中：Vf、uf、Vs、us分別為液體的速度、位移，固體的速度、位移.

動力學條件指的是耦合面上的力必須處于平衡狀態(tài)，可用以下表達式描述：

式中：n為法線方向；σf、σs分別為液體、固體的應力.

流體側(cè)和固體側(cè)通過在耦合面上進行變量交換來完成耦合.在每一個時間步內(nèi)，都有流體部分提供壓力和黏性力給固體，固體部分提供位移增量和速度給流體［20］.一個時間步內(nèi)的耦合求解流程如圖2所示.

圖2 一個時間步內(nèi)耦合求解流程Fig.2 The process of fluid-structure interaction solution within one time step

2.2 數(shù)值方法驗證

為了驗證本文中流固耦合的流場及固體變形計算的準確性，采用Krawczyk等［11］的研究，定量對比不同來流速度下柔性翼型變形的大小和升阻力系數(shù).計算的模型、參數(shù)等均與參考文獻保持一致.即采用C型計算區(qū)域，柔性翼型以弦長1 m的NACA4412為基礎(chǔ)翼型.本文取文獻中彈性模量為2.5 MPa下的五種來流速度，分別為3.79、5.04、6.309、7.59、10.22 m/s，對比其變形大小和升阻力系數(shù).

計算的結(jié)果與文獻中的數(shù)據(jù)對比如圖3所示.圖中給出的變形大小及升阻力系數(shù)都是經(jīng)過歸一化處理的.以柔性翼型的尾緣部分在來流速度為7.59 m/s下的變形大小為0 m，其余工況下的變形大小為與其的差值，其中變形為正代表尾緣部分相對變形方向朝上，變形為負代表尾緣部分相對變形方向朝下；而升阻力系數(shù)則為與剛性翼型在來流速度為7.59 m/s工況下所得的升阻力系數(shù)的比值.可以看出本文的計算結(jié)果與文獻相比整體吻合較好，可以證明本文對柔性翼型的流固耦合計算是準確的.

圖3 驗證算例結(jié)果對比Fig.3 Comparison of calculation results of verification cases

2.3 計算域網(wǎng)格劃分及邊界條件

計算域的左半邊為半圓形，其半徑為3.5 m，右半邊為長8.5 m、寬7 m 的矩形.翼型的尾部與半圓形的圓心重合，整個計算域的C形網(wǎng)格劃分如圖4所示.左端設定為速度入流邊界條件：u=4.5 m/s，v=0 m/s；上下設定為對稱邊界條件：?u/?y=v=0；右端出口設定為紐曼邊界條件：?u/?x=?v/?x=0；翼型則設定為流固耦合邊界條件.為了使邊界層區(qū)域內(nèi)的流場求解更加精確，壁面第一層網(wǎng)格y+設置為0.7，增長率取1.15，同時，為了驗證網(wǎng)格的無關(guān)性，對計算域劃分了三種不同數(shù)量的網(wǎng)格，網(wǎng)格數(shù)量分別為15 700、27 890和43 560并進行了計算.結(jié)果表明，后兩者的升阻力系數(shù)最大僅相差0.5%，固體最大變形相差0.7%，而前兩者升阻力系數(shù)相差23%，固體最大變形相差18%，考慮到計算資源和時間成本，本文采用網(wǎng)格數(shù)量為27 890的網(wǎng)格進行計算.對于時間步長，取庫朗數(shù)為0.1時計算得出的時間步長為1.5×10-5s，滿足計算所需的精度要求.

圖4 計算域網(wǎng)格Fig.4 Grid of the computational domain

3 計算結(jié)果與分析

3.1 不同彈性模量下的固定攻角柔性翼型氣動特性對比

為了便于研究減縮頻率對翼型動態(tài)失速特性的影響，首先對四組不同彈性模量下的固定攻角柔性翼型進行模擬，分析不同彈性模量下的升阻比，選取使翼型具有最佳氣動性能的彈性模量作為動態(tài)失速時柔性翼型的彈性模量.

如圖5 所示為四種彈性模量下的升阻比.可以看出，在2.5°～17.5°攻角范圍內(nèi)，四種彈性模量的柔性翼型的升阻比均高于剛性翼型，這表明柔性翼型相對于剛性翼型普遍具有較優(yōu)的氣動性能.在本次模擬的四種彈性模量中，彈性模量為3×107Pa 的柔性翼型在5°～17.5°攻角范圍內(nèi)所具有的氣動特性都是最優(yōu)的，因此，對于俯仰振蕩中的柔性翼型，彈性模量均選用3×107Pa.

圖5 柔性翼型在不同彈性模量下的升阻比Fig.5 Lift-to-drag ratios of flexible airfoil under different elastic moduli

3.2 不同減縮頻率下剛?cè)嵋硇偷臍鈩犹匦苑治?/h3>

柔性翼型在不同減縮頻率下的氣動特性具有明顯的區(qū)別，且同一減縮頻率下，剛性翼型和柔性翼型也存在區(qū)別.

圖6 所示為剛性翼型和柔性翼型在不同減縮頻率下的升阻力系數(shù)遲滯回線圖.柔性翼型的升阻力系數(shù)遲滯回線除了局部具有波動外，整體變化趨勢與剛性翼型保持一致，均有以下規(guī)律：當k=0.08 時，升力系數(shù)在上仰過程中一直保持較低的狀態(tài)，在下俯過程中一直保持較高的狀態(tài)，阻力系數(shù)在整個過程中沒有出現(xiàn)大幅度的改變；當k=0.024時，在攻角大于7°后升力系數(shù)一直大于k=0.08時的值，但在攻角變小的過程中卻突然驟降，直到18.5°時才回升，阻力系數(shù)在20°之后明顯提升；當k=0.009 7時，上仰過程中升力系數(shù)在22°時發(fā)生了大幅度下降，在24°時出現(xiàn)回升，下俯過程中由25°至24°時有短暫的上升，隨后一直處于下降趨勢，阻力系數(shù)的變化趨勢則與升力系數(shù)相同.

圖6 剛性翼型和柔性翼型升阻力系數(shù)遲滯回線Fig.6 Hysteresis loops of the lift-drag coefficients for rigid airfoil and flexible airfoil

再通過比較同一減縮頻率下，剛性翼型和柔性翼型的升阻力系數(shù)遲滯回線可以發(fā)現(xiàn)，在k較小時，二者無明顯區(qū)別；而k較大時，即k=0.08時，柔性翼型的升力系數(shù)整體略高于剛性翼型，阻力系數(shù)整體小于剛性翼型，且阻力系數(shù)的最大值明顯小于剛性翼型.由此可見，在高減縮頻率下，柔性翼型的氣動特性要顯著優(yōu)于剛性翼型，具有高升阻比和低阻力的特點.

3.3 柔性翼型的結(jié)構(gòu)響應分析

在柔性翼型的俯仰過程中，柔性翼型會受到流體的作用力而發(fā)生自適應變形，柔性翼型的變形具有一定的規(guī)律.

圖7 展示了柔性翼型在k=0.08 時整個俯仰過程中尾緣變形量的變化，其中上仰至25°攻角時的變形及應力在圖8 中進行展示.通過圖8 中展示的應力可以看出高應力主要集中在矩形固定孔洞附近及離孔洞最近的上下表面處，從孔洞至尾緣應力逐漸減小.觀察圖8中的變形分布規(guī)律可以明顯看到，在整個過程中柔性翼型的變形主要集中在尾緣部分，其余部位的變形量接近零.變形量從翼型的中間至尾緣逐漸增大，在尾部達到最大值.從定性的角度來看，變形的分布規(guī)律與Macphee和Beyene［9］的研究結(jié)果不謀而合，但從定量的角度來分析，動態(tài)失速中的柔性翼型的變形在數(shù)值及變化規(guī)律上又有其獨有的特征.通過圖7可以看出在整個俯仰過程中，變形的大小隨著攻角的增大（減小）快速地上下波動，這表明柔性翼型的變形十分靈活，對氣動荷載較為敏感.在上仰過程中，隨著攻角的增大，變形的數(shù)值整體趨勢是略有降低的，而在下俯過程中卻相反，隨著攻角的減小，其變形的整體趨勢也隨之減小.上仰階段和下俯階段的變形差距存在一定的規(guī)律，在小攻角范圍內(nèi)（0°～12.5°），上仰階段和下俯階段變形大小的差距并不明顯，而在大攻角范圍內(nèi)（12.5°～25°），下俯階段的變形在數(shù)值上大約是上仰階段的1.5倍至2倍，差距十分明顯.

圖7 柔性翼型俯仰過程中尾部變形量的變化Fig.7 Deformation variations of flexible airfoil trailing edges during pitching

圖8 柔性翼型上仰至25°攻角的變形及應力Fig.8 Deformations and stresses of flexible airfoil pitching up to 25°angle of attack

3.4 不同減縮頻率下剛?cè)嵋硇偷牧鲌龇治?/h3>

本小節(jié)將通過分析流場來對柔性翼型動態(tài)失速的機理以及剛?cè)嵋硇椭g氣動特性存在差異的原因進行研究.

圖9為柔性翼型在k=0.08時上仰過程的渦量和流線.圖10為柔性翼型在k=0.009 7時上仰過程的渦量和對應的流線.通過對比圖9和圖10中的流場可以發(fā)現(xiàn)，上仰至15°攻角時，二者尾流均出現(xiàn)了波動，尾緣上表面出現(xiàn)回流現(xiàn)象，且k=0.08的流場尾流部分已經(jīng)出現(xiàn)渦街；上仰至25°時，二者出現(xiàn)了明顯區(qū)別，k=0.08的流場尾流部分的渦街增強，翼型尾緣上表面出現(xiàn)流動分離現(xiàn)象，而k=0.009 7的流場已經(jīng)出現(xiàn)渦旋從翼型表面大規(guī)模脫落的情況，翼型已經(jīng)處于失速狀態(tài)，這導致升力系數(shù)在上仰過程中突然下降.低減縮頻率下的翼型在上仰過程出現(xiàn)失速現(xiàn)象，表明了高減縮頻率下的翼型相比低減縮頻率具有失速延遲的特性.

圖9 柔性翼型上仰過程渦量、流線、變形及應力，k=0.08Fig.9 Vorticity contours，streamlines，deformations and stresses of the flexible airfoil during pitching up，k=0.08

圖10 柔性翼型上仰過程渦量及流線，k=0.009 7Fig.10 Vorticity contours and streamlines of the flexible airfoil during pitching up，k=0.009 7

圖11所示為剛性翼型在k=0.08時上仰至25°攻角的渦量，與圖9（b）中相同條件下的柔性翼型相比，二者的區(qū)別十分明顯，剛性翼型此時已經(jīng)在上表面產(chǎn)生了集中強渦旋，而柔性翼型上表面僅僅是尾緣部分出現(xiàn)了流動分離，并未出現(xiàn)渦旋.圖8展示了此時柔性翼型因為受到流體的作用力所產(chǎn)生的自適應變形，該變形會反作用于流場而改變流場的特性，此時柔性翼型的變形抑制了翼型表面渦旋的生成，對流場起到了一定的調(diào)控作用，使柔性翼型表面僅出現(xiàn)了流動分離.通過對流場差異性的分析可以得出，剛?cè)嵋硇偷臍鈩犹匦圆町愂怯勺冃螌α鲌龅恼{(diào)控導致的.再通過對應圖6（a）和（b）中上仰至25°時的升阻力系數(shù)發(fā)現(xiàn)柔性翼型的升力系數(shù)僅有微小提升，而阻力系數(shù)出現(xiàn)了大幅度降低.

圖11 剛性翼型上仰至25°攻角渦量及流線，k=0.08Fig.11 Vorticity contours and streamlines of the rigid airfoil pitching up to 25°of attack，k=0.08

4 結(jié)論

本文采用流固耦合數(shù)值模擬方法，研究了柔性翼型在不同減縮頻率下其升阻力系數(shù)和流場的變化規(guī)律，通過結(jié)構(gòu)響應分析揭示了剛?cè)嵋硇蜌鈩犹匦跃哂忻黠@差異的物理機理.研究發(fā)現(xiàn)即使在動態(tài)失速中柔性翼型依然具有較優(yōu)的氣動性能，這對未來翼型的發(fā)展具有一定的指導意義.具體的研究結(jié)論如下：

1）柔性翼型相比剛性翼型在氣動特性上所展現(xiàn)的優(yōu)勢是大幅度減阻和小幅度增升.2）柔性翼型的被動變形主要產(chǎn)生在尾緣部分，變形在動態(tài)失速下的分布規(guī)律與前人研究的靜態(tài)中的變形分布規(guī)律相同，但區(qū)別在于動態(tài)失速中變形的大小會隨攻角的變化上下波動.3）柔性翼型良好的氣動特性取決于柔性翼型變形對流場調(diào)控的結(jié)果.相比于剛性翼型，該變形能夠抑制表面渦旋的產(chǎn)生.