基于浸入邊界法的翼型俯仰振擺動(dòng)態(tài)特性

闞 闞，李昊宇，呂 品，徐宇航，鄭 源，徐 輝

（1. 河海大學(xué)能源與電氣學(xué)院，南京 211100；2. 河海大學(xué)水利水電學(xué)院，南京 210098；3. 中國(guó)船舶集團(tuán)有限公司第七一〇研究所，宜昌 056002；4. 河海大學(xué)農(nóng)業(yè)科學(xué)與工程學(xué)院，南京 211100）

0 引 言

翼型俯仰振擺會(huì)引起流場(chǎng)顯著的非線性變化，在氣動(dòng)彈性和非定常流體動(dòng)力學(xué)領(lǐng)域受到越來越多的重視，被認(rèn)為是噴氣發(fā)動(dòng)機(jī)、風(fēng)力發(fā)電機(jī)和水力機(jī)械等動(dòng)力設(shè)備設(shè)計(jì)中的關(guān)鍵問題之一。在工程中翼型的振擺通常分為受迫振擺和主動(dòng)振擺兩種情況。其一，風(fēng)力機(jī)、水輪機(jī)等流體機(jī)械在運(yùn)行中，由于結(jié)構(gòu)約束不足或者變槳葉調(diào)節(jié)機(jī)構(gòu)失靈，葉片可能發(fā)生一定程度的振擺，從而導(dǎo)致流體機(jī)械的動(dòng)態(tài)性能發(fā)生劇烈波動(dòng)，甚至造成流體機(jī)械的破壞。其二，主動(dòng)振擺的襟翼可以用來調(diào)節(jié)風(fēng)力機(jī)、飛行器以及水下航行器等設(shè)備動(dòng)態(tài)性能，是近年來研究的工程熱點(diǎn)。因此對(duì)翼型俯仰振擺動(dòng)態(tài)特性的研究具有極高的理論價(jià)值和實(shí)際工程意義。

針對(duì)翼型作俯仰振擺運(yùn)動(dòng)的動(dòng)態(tài)特性，國(guó)內(nèi)外學(xué)者主要圍繞振擺頻率、振擺幅值、振擺平均攻角及俯仰軸位置等參數(shù)展開了研究。Amiralaei等對(duì)低雷諾數(shù)下翼型俯仰振擺的研究指出振擺頻率和振擺幅值對(duì)動(dòng)力學(xué)性能有著重要的影響。Tian等發(fā)現(xiàn)隨著翼型俯仰軸位置的前移，翼型產(chǎn)生的推力增加。薛臣等指出振動(dòng)軸位置主要改變翼型的有效迎角，而振擺頻率的增大會(huì)導(dǎo)致翼型周圍流場(chǎng)結(jié)構(gòu)產(chǎn)生變化。李曉華等通過數(shù)值模擬對(duì)俯仰振擺減縮頻率和平均攻角的研究表明減縮頻率對(duì)翼型氣動(dòng)特性的影響主要在翼型攻角的下行段。Kurtulus發(fā)現(xiàn)在小幅值的振擺中，振擺頻率會(huì)對(duì)升力波動(dòng)幅值產(chǎn)生影響。而對(duì)于上述參數(shù)對(duì)翼型俯仰振擺動(dòng)態(tài)特性的影響機(jī)理，前人也進(jìn)行了許多研究。Worasinchai等利用風(fēng)洞實(shí)驗(yàn)方法研究了垂直軸風(fēng)力機(jī)葉片翼型在高頻振擺運(yùn)動(dòng)中的表面非穩(wěn)態(tài)的壓力變化，研究指出了翼型的振擺運(yùn)動(dòng)導(dǎo)致的前緣渦的形成與脫落對(duì)于翼型表面的壓力分布有顯著影響。李紹斌等通過數(shù)值模擬研究了平均攻角和振幅對(duì)俯仰振蕩翼型氣動(dòng)特性的影響，指出翼型振擺過程中上仰階段前緣渦的產(chǎn)生和集中渦的穩(wěn)定附著是平均升力系數(shù)隨振幅增大呈階躍式提升的原因。

雖然對(duì)于翼型俯仰振擺動(dòng)態(tài)特性已經(jīng)開展了很多研究，但目前關(guān)于翼型在俯仰振擺運(yùn)動(dòng)過程中升力系數(shù)發(fā)生的波動(dòng)與流場(chǎng)中旋渦發(fā)展演化的相關(guān)性方面的研究仍然不夠深入。此外，前人對(duì)俯仰振擺翼型的研究所采用的數(shù)值方法大多基于貼體網(wǎng)格，其在處理快速移動(dòng)邊界時(shí)存在一定的局限性，可能因?yàn)榫W(wǎng)格拉伸過大、網(wǎng)格質(zhì)量下降等原因影響模擬結(jié)果的精確性。而浸入邊界法（Immersed Boundary Method，IBM），是處理此類動(dòng)邊界問題的一種有效方法，由Peskin于1972年首次提出。相比傳統(tǒng)的基于貼體網(wǎng)格的數(shù)值方法，因其無需生成貼體網(wǎng)格和進(jìn)行網(wǎng)格移動(dòng)與重生，從而具有較高的計(jì)算效率，可以相對(duì)高效地處理具有移動(dòng)邊界的流動(dòng)。王文全等通過自編程序驗(yàn)證了浸入邊界法可用于模擬流體剛體的耦合運(yùn)動(dòng)，但是采用浸入邊界法對(duì)作俯仰振擺運(yùn)動(dòng)的翼型的動(dòng)態(tài)特性進(jìn)行的研究還較少。

因此本文采用基于浸入邊界法的自編求解器，對(duì)雷諾數(shù)=1000條件下NACA0012翼型以弦長(zhǎng)1/2處為中心作不同初始攻角、不同振擺頻率、不同振擺幅值的俯仰振擺運(yùn)動(dòng)進(jìn)行二維流動(dòng)直接數(shù)值模擬（Direct Numerical Simulation，DNS），分析各工況下升力系數(shù)的波動(dòng)特性和升力系數(shù)波動(dòng)與流場(chǎng)演變的相關(guān)性，同時(shí)驗(yàn)證浸入邊界法對(duì)瞬態(tài)流動(dòng)捕捉的準(zhǔn)確性。研究結(jié)果擬為工程中風(fēng)力機(jī)等流體機(jī)械翼型顫振引起的流動(dòng)不穩(wěn)定性研究提供理論參考和工程依據(jù)。

1 數(shù)值方法及驗(yàn)證

1.1 控制方程

不可壓縮流體不考慮密度的變化，質(zhì)量守恒方程（連續(xù)性方程）為

動(dòng)量守恒方程（Navier-Stokes方程，N-S方程）為

式中為速度，m/s；為流體密度，kg/m；為流體運(yùn)動(dòng)黏度，m/s；為壓力，Pa；為對(duì)浸入邊界（Immersed Boundary，IB）點(diǎn)施加的外部源項(xiàng)體積力，m/s。

1.2 空間離散

采用直角坐標(biāo)系下的笛卡爾交錯(cuò)網(wǎng)格，基于二階中心差分格式的有限差分法進(jìn)行離散。、、分別為速度在軸方向的正交分量，定義在網(wǎng)格單元體面心；而其余物理量例如密度、壓力等定義在網(wǎng)格體心，如圖1所示。

圖1 空間離散網(wǎng)格示意圖 Fig.1 Schematic diagram of spatial discrete grid

1.3 時(shí)間步推進(jìn)與壓力泊松方程求解

采用投影法與二階龍格-庫塔法（Runge-Kutta Method，RKM）推進(jìn)計(jì)算時(shí)間步。二階龍格-庫塔法將求解分為兩個(gè)時(shí)間子步進(jìn)行迭代求解，在每一時(shí)間子步中利用投影法分步計(jì)算N-S方程各項(xiàng)對(duì)速度場(chǎng)的影響并疊加推進(jìn)計(jì)算進(jìn)行，同時(shí)保證無散度條件。以第一步龍格-庫塔法為例進(jìn)行說明（下文說明均以第一步龍格-庫塔法為例），其過程如公式（3）至（6）所示。

式中各物理項(xiàng)的上標(biāo)為時(shí)間子步序號(hào)。u為計(jì)算第步的速度場(chǎng)，u為龍格庫塔法的中間速度，Δ為求解時(shí)間步長(zhǎng)。N-S方程右側(cè)的粘性項(xiàng)?與對(duì)流項(xiàng)-(??)之和記為RHS（Right Hand Side）。u表示時(shí)間子部中忽略固體存在性的速度場(chǎng)。對(duì)IB點(diǎn)進(jìn)行速度插值并通過公式（4）得到更新后的速度場(chǎng)~。求解壓力泊松方程時(shí)采用的求解器為 PETSc（Portable，ExtensibleToolkit for Scientific Computation），即科學(xué)計(jì)算可移植擴(kuò)展工具包，可以對(duì)偏微分方程組高效并行求解。為減少泊松方程迭代次數(shù)，根據(jù)公式（5），泊松方程求解結(jié)果為各時(shí)間子步計(jì)算前后的壓力差Δ，即相對(duì)壓力。將求解泊松方程得到的相對(duì)壓力Δp與p相加獲得+1/2時(shí)刻無量綱絕對(duì)壓力p，若求解節(jié)點(diǎn)為流體點(diǎn)則忽略體積力

1.4 基于水平集函數(shù)的浸入邊界法

為使流固交界處流體速度滿足無滑移邊界條件，在IB點(diǎn)施加一個(gè)外部源項(xiàng)的體積力，在求解N-S方程時(shí)需得到此力對(duì)流場(chǎng)的影響以更新速度場(chǎng)。

采用Sussman等的水平集方法（Level-Set）用于翼型界面的追蹤，將NACA0012翼型邊界定義為零水平集，見公式（7）。

首先需確定各節(jié)點(diǎn)的水平集函數(shù)值，將各速度分量所在節(jié)點(diǎn)分類為固體點(diǎn)、流體點(diǎn)與浸入邊界點(diǎn)。特別地，因使用交錯(cuò)網(wǎng)格，速度分量定義在網(wǎng)格面心，而水平集函數(shù)定義在網(wǎng)格體心，故在定義速度節(jié)點(diǎn)類型時(shí)需要通過相鄰的兩個(gè)網(wǎng)格的水平集函數(shù)值通過線性插值的方法得到網(wǎng)格體面的水平集函數(shù)值。水平集函數(shù)值定義如公式（8）所示。

式中為此節(jié)點(diǎn)沿水平集函數(shù)梯度方向相對(duì)固液邊界的坐標(biāo)。若一節(jié)點(diǎn)<0，則將其定義為固體點(diǎn)。在二維情況中，若一節(jié)點(diǎn)>0，且其4個(gè)相鄰節(jié)點(diǎn)均滿足>0，則將其定義為流體點(diǎn)，若一節(jié)點(diǎn)>0，且其4個(gè)相鄰中存在滿足<0的節(jié)點(diǎn)，則將其定義為IB點(diǎn)，見圖2。

圖2 浸入邊界法示意圖 Fig.2 Schematic diagram of immersed boundary method

以IB點(diǎn)4的計(jì)算為例，利用流固耦合邊界兩側(cè)節(jié)點(diǎn)水平集函數(shù)值求梯度，得到IB點(diǎn)4沿梯度向量在流固邊界上的投影點(diǎn)1，點(diǎn)2、3為對(duì)點(diǎn)4進(jìn)行速度插值的模板點(diǎn)。記IB點(diǎn)4在整個(gè)計(jì)算域的坐標(biāo)為(,)，結(jié)合水平集函數(shù)梯度以及網(wǎng)格長(zhǎng)度可計(jì)算得到點(diǎn)1、2、3的坐標(biāo)，分別為(,)、(,)、(,)。點(diǎn)1速度由固體本身的運(yùn)動(dòng)狀態(tài)決定，記為，點(diǎn)2、3速度分別記為、。通過求解線性方程組（9）得到插值系數(shù)、、，根據(jù)公式（10）對(duì)IB點(diǎn)速度基于無滑移邊界條件線性插值更新速度場(chǎng)。

2 計(jì)算設(shè)置與方法驗(yàn)證

2.1 計(jì)算設(shè)置

采用1 024×768×3的笛卡爾坐標(biāo)系下的交錯(cuò)網(wǎng)格，對(duì)NACA0012進(jìn)行二維流動(dòng)直接數(shù)值模擬。取典型低雷諾數(shù)流場(chǎng)，雷諾數(shù)如公式（11）計(jì)算所得。

式中U為來流速度，0.146 m/s；為特征長(zhǎng)度，此處即為翼型弦長(zhǎng)，0.01 m；為流體的運(yùn)動(dòng)黏度，此處取值1.46×10m/s，即雷諾數(shù)=1 000。為確保翼型振擺過程中的尾流場(chǎng)捕捉的準(zhǔn)確性，同時(shí)減少計(jì)算域邊界對(duì)翼型周圍流場(chǎng)的影響，計(jì)算域中方向的長(zhǎng)度為10，方向的長(zhǎng)度為7。翼型周圍采用加密網(wǎng)格計(jì)算，加密中心為(4,3.5)，方向加密長(zhǎng)度為4，方向加密長(zhǎng)度為3。計(jì)算域進(jìn)口為正方向的均勻入流，計(jì)算域右側(cè)出口邊界條件為對(duì)流出口，上下邊界均為自由滑移邊界條件。計(jì)算域示意圖見圖3。

圖3 計(jì)算域示意圖 Fig.3 Schematic diagram of computing domain

翼型振擺中心位于弦長(zhǎng)1/2處，攻角為，振擺由公式（12）控制，其角速度見公式（13）。

式中為初始攻角，(°)；為最大振擺幅度，(°)；為振擺頻率，Hz；為時(shí)間，s。

2.2 方法驗(yàn)證

本次研究首先通過對(duì)不同攻角下的NACA0012靜止翼型的升力系數(shù)、阻力系數(shù)、斯特勞哈爾數(shù)進(jìn)行計(jì)算并與已有文獻(xiàn)對(duì)比，以驗(yàn)證計(jì)算數(shù)值方法、計(jì)算網(wǎng)格的準(zhǔn)確性及可靠性。為進(jìn)行驗(yàn)證對(duì)比，結(jié)合翼型較為常見的運(yùn)行工況，分別選取靜止?fàn)顟B(tài)下攻角分別為=5°、10°、15°、20°的情況進(jìn)行流動(dòng)模擬。升力系數(shù)、阻力系數(shù)及斯特勞哈爾數(shù)的定義見公式（14）～（16）。

式中F、分別為翼型受到的升力和阻力；為翼型尾部渦脫落頻率。本文采用的求解器計(jì)算過程中各物理量均為無量綱，速度，長(zhǎng)度，時(shí)間的無量綱化過程見公式（17）。

式中為無量綱化后的速度；為無量綱化后的長(zhǎng)度；為無量綱化后的時(shí)間，同時(shí)無量綱壓力p與無量綱渦量定義如下：

式中為渦量，1/s。

本次計(jì)算結(jié)果與其他文獻(xiàn)中數(shù)值模擬的結(jié)果對(duì)比見表1，此處升力系數(shù)、阻力系數(shù)為時(shí)間平均系數(shù)。由表可見，本次計(jì)算結(jié)果均與各文獻(xiàn)結(jié)果吻合良好。

表1 不同攻角平均升力系數(shù)、阻力系數(shù)、斯特勞哈爾數(shù)對(duì)比 Table 1 Comparison of averaged lift coefficient, drag coefficient and Strouhal number at different attack angles (α)

同時(shí)，為進(jìn)一步驗(yàn)證求解器對(duì)流場(chǎng)信息捕捉的準(zhǔn)確性，選取=5°、=10°、=15°的渦量（無量綱）瞬時(shí)分布圖（圖4），升力系數(shù)隨時(shí)間變化圖（圖5a）和頻譜圖（圖5b）進(jìn)行觀察分析。由圖可見，5°攻角時(shí)，吸力面與壓力面邊界層穩(wěn)定，均未發(fā)生流動(dòng)分離，翼型尾部流跡穩(wěn)定，無脫落渦產(chǎn)生，計(jì)算收斂后流動(dòng)穩(wěn)定且升力系數(shù)不發(fā)生波動(dòng)。10°攻角時(shí)，吸力面產(chǎn)生了明顯的流動(dòng)分離，分離點(diǎn)由后緣前移并產(chǎn)生了周期性旋渦脫落，壓力場(chǎng)隨旋渦脫落發(fā)生變化使得升力系數(shù)曲線產(chǎn)生波動(dòng)，由升力系數(shù)頻譜可得旋渦脫落頻率為13.255 Hz。在翼型攻角增加至15°時(shí)，吸力面流動(dòng)分離點(diǎn)繼續(xù)向前緣移動(dòng)，脫落渦尺度和強(qiáng)度增大，渦形成位置更加靠近尾緣，翼型壓力面邊界層形成明顯的逆壓梯度流動(dòng)，頻譜圖中出現(xiàn)頻率為主頻10.309 Hz及其倍數(shù)的次頻。通過與已有文獻(xiàn)數(shù)值模擬結(jié)果升力系數(shù)、斯特勞哈爾數(shù)等參數(shù)以及渦量分布的對(duì)比結(jié)果表明本求解器可以準(zhǔn)確捕捉翼型周圍流體流動(dòng)過程中的升力波動(dòng)與旋渦的產(chǎn)生。

圖4 靜止?fàn)顟B(tài)不同攻角下瞬時(shí)渦量云圖 Fig.4 Instantaneous vorticity nephogram at different attack angles at static condition

圖5 不同攻角靜止?fàn)顟B(tài)升力系數(shù)變化曲線圖及升力系數(shù) 頻譜圖 Fig.5 Diagram of variation curve and spectrum of the lift coefficient at different attack angles at static condition

3 試驗(yàn)與結(jié)果

3.1 翼型振擺下的升力系數(shù)時(shí)頻特性

翼型振擺主要將初始攻角、振擺頻率、振擺幅值這三個(gè)主要因素作為變量進(jìn)行考慮。根據(jù)翼型在實(shí)際工程的主要運(yùn)行工況和相關(guān)文獻(xiàn)，初始攻角選擇=5°（小攻角）、=10°（中攻角）、=15°（大攻角），振擺頻率選擇1.46 Hz（低頻）、2.92 Hz（高頻），振擺幅值選擇=5°（小幅值）、=10°（大幅值）。各工況翼型振擺運(yùn)動(dòng)參數(shù)及平均升力系數(shù)見表2。

表2 不同工況翼型振擺運(yùn)動(dòng)參數(shù)及平均升力系數(shù) Table 2 Pitching motion parameters and averaged lift coefficient under various operating conditions

在整個(gè)振擺周期中平均升力系數(shù)主要受到振擺初始攻角影響，隨著的增加，翼型平均升力系數(shù)增大。在其余條件相同的情況下，隨著振擺頻率的增加，平均升力系數(shù)均增加，且振擺頻率對(duì)大初始攻角下的平均升力系數(shù)增加的影響更加顯著，這與白鵬等研究的結(jié)論一致。振擺頻率由1.46 Hz增加至2.92 Hz引起=5°情況升力系數(shù)增加7.9%，=10°情況升力系數(shù)增加5.2%，=15°情況升力系數(shù)增加17.2%，對(duì)于中、小攻角的影響均在10%以內(nèi)，但對(duì)于大攻角其影響大于10%。與靜止情況相比，除=5°與=15°的高頻小幅值振擺平均升力系數(shù)分別保持不變、增加4.3%外，其余情況均有所減小，減小程度均在13%以內(nèi)，因=5°低頻大幅值振擺中出現(xiàn)負(fù)攻角，故不在考慮范圍內(nèi)。

圖6與圖7分別為在不同初始攻角下翼型進(jìn)行低頻小幅值振擺（工況1～3），高頻小幅值振擺（工況4～6），高頻大幅值振擺（工況7～9）時(shí)的升力系數(shù)隨時(shí)間變化圖與頻譜圖。總體來說，受到翼型振擺的影響，其升力系數(shù)均呈現(xiàn)周期性波動(dòng)的情況，且波動(dòng)主頻均為翼型振擺運(yùn)動(dòng)頻率，同時(shí)頻譜圖中的出現(xiàn)的次頻均為翼型振擺頻率的整數(shù)倍，說明其升力系數(shù)的波動(dòng)現(xiàn)象的發(fā)生均受到翼型振擺頻率的影響。

對(duì)比圖6和圖7中不同初始攻角下的情況，隨著初始攻角的增加，頻譜圖中出現(xiàn)的高次次頻數(shù)量明顯增加，升力系數(shù)的波動(dòng)現(xiàn)象更加劇烈。在較大的初始攻角下，出現(xiàn)了更多與相關(guān)的倍頻，可能是因?yàn)檩^大攻角下翼型繞流中更易發(fā)生流動(dòng)分離形成小尺度旋渦，在振擺的影響下形成翼型表面壓力分布的波動(dòng)。

對(duì)比圖7a、7b，隨著振擺頻率增加，在=15°的情況下，升力系數(shù)波動(dòng)頻率主頻幅值和次頻2幅值均顯著增加，次頻數(shù)量明顯減少。但是對(duì)于較小的初始攻角=5°、=10°，振擺頻率的增加減少了次頻數(shù)量，而對(duì)主頻和次頻幅值的影響較小，說明較大的初始攻角下翼型升力系數(shù)波動(dòng)受高頻振擺的影響更為顯著。對(duì)比圖6a、圖6b，振擺頻率的變化明顯改變了中（=10°）、大（=15°）初始攻角時(shí)升力系數(shù)周期變化的整體趨勢(shì)。楊筱沛等指出翼型升力系數(shù)的變化規(guī)律與翼型尾渦脫落規(guī)律保持一致，因此振擺頻率的增加使較低頻波動(dòng)波峰波谷的數(shù)量明顯減少，有可能是因?yàn)殡S著振擺頻率的增加，當(dāng)繞流流體通過翼型尾緣處時(shí)，翼型邊界與流動(dòng)的相對(duì)速度增大，在對(duì)流影響下形成的旋渦在尾緣處尚未充分演化便迅速脫落，從而減弱了升力系數(shù)的低頻大幅值波動(dòng)。

對(duì)比圖6a與圖6c、圖7a與圖7c，可以得出振擺幅值對(duì)翼型升力系數(shù)的影響。在相同振擺頻率下，振擺幅值的增加使三種攻角下頻譜圖中高次次頻的數(shù)量增加，同時(shí)各頻率幅值也有所增加。雖然主頻幅值基本不受初始攻角影響，但隨著初始攻角增大，次頻中與相關(guān)的倍頻增多，這與低頻小幅值振擺的情況類似。這表現(xiàn)出在翼型低頻振擺情況下，振擺幅值占據(jù)主導(dǎo)作用。

圖6 不同工況升力系數(shù)變化曲線 Fig.6 Change curves of lift coefficient variation at different conditions

圖7 不同工況攻角升力系數(shù)頻譜圖 Fig.7 Spectrum of lift coefficient at different conditions

3.2 升力系數(shù)波動(dòng)與流場(chǎng)演變的相關(guān)性

在前文中已經(jīng)說明翼型升力系數(shù)存在周期性波動(dòng)，以初始攻角=15°，振擺頻率=1.46 Hz，振擺幅值=5°的情況為例進(jìn)行分析，其在一個(gè)振擺周期內(nèi)的速度矢量分布、渦量分布（無量綱）、壓力分布（無量綱）見圖8。翼型升力系數(shù)在單個(gè)振擺周期存在著較低頻的大幅值波動(dòng)，引起這種波動(dòng)的原因可能為壓力面尾緣形成的逆時(shí)針旋渦在翼型吸力面尾緣處的持續(xù)演化發(fā)展。

圖8 單個(gè)振擺周期內(nèi)速度矢量圖、渦量、壓力隨時(shí)間變化分布圖 Fig.8 Velocity vector diagram, vorticity and pressure distribution diagram with time in a pitching period

如圖8中的渦量分布云圖所示，翼型的一個(gè)振擺周期中，其壓力面均未發(fā)生流動(dòng)分離，壁面邊界層流動(dòng)穩(wěn)定。至+0.250區(qū)間內(nèi)翼型攻角增加并逐漸達(dá)到最大攻角，結(jié)合速度矢量與壓力分布可知在此期間翼型吸力面前緣處流體速度受到壁面的影響逐漸加劇，形成強(qiáng)烈的逆壓梯度，促發(fā)了流動(dòng)分離，同時(shí)由于沿壁面法向速度梯度的存在，吸力面前緣處流體受到主流的剪切力，在順壓梯度的驅(qū)動(dòng)下形成順時(shí)針的負(fù)值旋渦流動(dòng)。當(dāng)翼型攻角發(fā)生變化時(shí)，受到振擺翼型壁面轉(zhuǎn)動(dòng)的影響，流體壓力場(chǎng)發(fā)生變化，尾緣處壓力降低誘發(fā)壓力面近壁流動(dòng)卷起，形成逆時(shí)針正值旋渦流動(dòng)。至+0.250時(shí)刻，翼型尾緣處形成了同時(shí)覆蓋壓力面與吸力面的連續(xù)低壓區(qū)，出現(xiàn)正值旋渦回流至吸力面的情況，并與前緣產(chǎn)生的負(fù)值渦相互作用。+0.250至+0.750區(qū)間內(nèi)翼型攻角減小并達(dá)到最小攻角，翼型壓力面不斷有正值逆時(shí)針旋渦卷起至吸力面并與阻礙前緣負(fù)值旋渦向下游的發(fā)展，同時(shí)伴隨著尾緣處覆蓋壓力面與吸力面的連續(xù)負(fù)壓區(qū)域的出現(xiàn)與消失。單個(gè)振擺周期內(nèi)正負(fù)值渦在吸力面中的相互作用主要發(fā)生在此區(qū)間內(nèi)，同時(shí)由速度矢量圖可知在此期間翼型吸力面尾緣上方旋渦流動(dòng)區(qū)域相比于其他時(shí)間區(qū)間更大。+0.750至+1.000區(qū)間內(nèi)翼型逐漸回到起始攻角，此區(qū)間內(nèi)尾部仍存在交替的脫落渦，但此區(qū)間內(nèi)尾部產(chǎn)生的連續(xù)負(fù)壓區(qū)負(fù)壓值、正負(fù)值渦量范圍均較小。

圖9為初始攻角=15°的低頻小幅值振擺（工況3）單個(gè)振擺周期升力系數(shù)變化曲線圖，圖10為升力系數(shù)處于波峰、波谷時(shí)刻的渦量與壓力分布云圖。由圖可見，在升力系數(shù)處于峰值的時(shí)刻，尾緣處由壓力面產(chǎn)生并卷起至吸力面的正值渦均正處于從翼型脫落的過程中，且渦核位置已經(jīng)脫離尾緣，其所處位置的負(fù)壓區(qū)對(duì)翼型表面的壓力分布的影響減弱，而此刻吸力面形成的負(fù)值前緣分離渦受到卷起渦的阻礙并未脫離，其所在位置仍存在一個(gè)負(fù)壓區(qū)，從而導(dǎo)致了翼型兩面壓力差處于峰值。之后隨著壓力面下一個(gè)的負(fù)值渦的產(chǎn)生并卷起至吸力面，并與吸力面負(fù)值渦相互作用，誘導(dǎo)出兩個(gè)負(fù)值的分裂渦，單個(gè)渦的強(qiáng)度下降，升力系數(shù)逐漸下降。在靠近下游的分裂渦從翼型脫落的同時(shí)，升力系數(shù)達(dá)到谷值，由此時(shí)的壓力分布圖可知，此時(shí)翼型尾緣處形成了同時(shí)覆蓋壓力面與吸力面的的連續(xù)低壓區(qū)，引起了兩面壓力差的衰減。此時(shí)刻后隨著流動(dòng)的發(fā)展，卷起至吸力面的渦經(jīng)尾緣處脫離，翼型升力系數(shù)逐漸增加直至下一個(gè)波峰。

圖9 工況3單個(gè)周期內(nèi)升力系數(shù)變化曲線 Fig.9 Variation curves of lift coefficient in a single pitching period of condition 3

圖10 升力系數(shù)峰值、谷值時(shí)刻渦量及壓力分布圖 Fig.10 Distribution diagram of vorticity and pressure at peak and valley of lift coefficient

根據(jù)湯姆孫定理，正壓有勢(shì)的理想流體繞翼型速度環(huán)量不發(fā)生變化，而在非理想流體中黏性存在導(dǎo)致的旋渦流動(dòng)改變了速度環(huán)量，升力由繞翼型的逆時(shí)針?biāo)俣拳h(huán)量提供。圖11為升力峰值谷值時(shí)刻翼型周圍流體速度矢量圖，由圖分析可知，當(dāng)升力系數(shù)處于峰值時(shí)刻，尾緣處存在一個(gè)順時(shí)針的旋渦流動(dòng)，在壁面處其速度方向與來流方向相反，從而導(dǎo)致了繞翼型逆時(shí)針?biāo)俣拳h(huán)量增加，進(jìn)而提高了翼型的升力；當(dāng)升力系數(shù)處于谷值時(shí)刻，翼型尾緣處僅存在一個(gè)逆時(shí)針的旋渦流動(dòng)，此流動(dòng)的存在減弱了速度環(huán)量，導(dǎo)致了翼型升力的衰減。

圖11 升力系數(shù)峰值、谷值時(shí)刻速度矢量分布圖 Fig.11 Distribution diagram of velocity vector at peak and valley of lift coefficient

因此，翼型在振擺過程中出現(xiàn)的升力大幅值波動(dòng)的原因?yàn)閴毫γ婺鏁r(shí)針旋渦由翼型尾緣分離后卷起至壓力面發(fā)展成逆壓梯度流動(dòng)，形成正值渦量，與翼型吸力面流動(dòng)分離脫落的負(fù)值渦旋渦發(fā)生強(qiáng)烈相互作用，在此過程中導(dǎo)致翼型壁面壓力分布劇烈變化引起了升力的大幅值波動(dòng)。分析結(jié)果可為翼型振擺過程中的升力不穩(wěn)定性問題的解決提供一定思路。

4 結(jié) 論

本文基于浸入邊界法對(duì)二維NACA0012翼型靜止和俯仰振擺下的低雷諾數(shù)瞬態(tài)流場(chǎng)進(jìn)行了數(shù)值模擬，驗(yàn)證了浸入邊界法對(duì)于瞬態(tài)流動(dòng)捕捉的準(zhǔn)確性。主要結(jié)論如下：

1）升力系數(shù)的波動(dòng)現(xiàn)象主要受翼型俯仰振擺頻率的影響。俯仰振擺頻率由1.46 Hz增加至2.92 Hz導(dǎo)致的翼型的平均升力系數(shù)增加范圍在5.2%～17.2%，并減弱了升力系數(shù)的次頻幅值，且影響程度隨著初始攻角的變大而增加。

2）在振擺頻率，振擺幅值相同時(shí)，隨著初始攻角的增加，翼型升力系數(shù)波動(dòng)現(xiàn)象更加劇烈，次頻增多且波動(dòng)幅值增大。在翼型低頻振擺的情況下，振擺幅值對(duì)升力系數(shù)波動(dòng)現(xiàn)象占據(jù)主導(dǎo)作用。

3）升力系數(shù)的周期性波動(dòng)是因?yàn)樾郎u演化導(dǎo)致的翼型邊界壓力分布的周期性變化。在翼型單個(gè)振擺周期中，中（=10°）、大（=15°）初始攻角下升力系數(shù)出現(xiàn)的大幅值次頻波動(dòng)，主要是由于壓力面的流動(dòng)分離及尾緣處負(fù)壓梯度回流引起尾緣旋渦卷起后脫落，進(jìn)而引起了翼型壓力面壓力分布波動(dòng)。