李寧宇, 蘇玉民，劉葳興，劉鵬，曹建

哈爾濱工程大學(xué) 水下機(jī)器人技術(shù)重點(diǎn)實(shí)驗(yàn)室，黑龍江 哈爾濱 150001

?

基于浸入邊界法的拍動(dòng)翼尾渦結(jié)構(gòu)與水動(dòng)力分析

李寧宇, 蘇玉民，劉葳興，劉鵬，曹建

哈爾濱工程大學(xué) 水下機(jī)器人技術(shù)重點(diǎn)實(shí)驗(yàn)室，黑龍江 哈爾濱 150001

摘要：針對目前拍動(dòng)翼尾渦結(jié)構(gòu)與水動(dòng)力性能關(guān)系的研究有限，參數(shù)影響分析集中于水動(dòng)力性能方面，且拍動(dòng)翼等仿生流動(dòng)的數(shù)值模擬經(jīng)常涉及復(fù)雜移動(dòng)邊界問題的現(xiàn)狀，應(yīng)用自主開發(fā)的改進(jìn)浸入邊界法數(shù)值模擬拍動(dòng)翼的非定常運(yùn)動(dòng)，探討翼運(yùn)動(dòng)學(xué)、渦動(dòng)力學(xué)和力的產(chǎn)生之間的關(guān)系。結(jié)果表明，所提出的邊界條件重建算法兼顧計(jì)算效率和健壯性，開發(fā)的移動(dòng)邊界處理方法可有效解決浸入邊界法中移動(dòng)邊界相關(guān)問題。三維拍動(dòng)翼的尾渦由兩列形狀復(fù)雜的渦環(huán)組成，這兩列渦環(huán)與尾流中心線成一定傾角向下游對流。水動(dòng)力性能隨運(yùn)動(dòng)參數(shù)的變化取決于力產(chǎn)生表象之下的渦動(dòng)力學(xué)，包括渦的強(qiáng)度和方向、渦環(huán)的相互連接和粘性耗散。

關(guān)鍵詞：拍動(dòng)翼；浸入邊界法；尾渦結(jié)構(gòu)；水動(dòng)力性能；推進(jìn)機(jī)理；不可壓縮N-S方程

蘇玉民(1960-), 男，教授，博士生導(dǎo)師.

盡管水下機(jī)器人研究領(lǐng)域有很多給人印象深刻的創(chuàng)新，但是軍方和科研團(tuán)體都希望能夠研究出更靈活機(jī)動(dòng)的水下機(jī)器人獲益，拍動(dòng)翼推進(jìn)器有望成為可行的技術(shù)方案之一。目前，國內(nèi)外學(xué)者已對仿生拍動(dòng)翼進(jìn)行了大量研究。Mackowski等[1]實(shí)驗(yàn)測量了拍動(dòng)翼的推力和效率，并分析運(yùn)動(dòng)參數(shù)對推進(jìn)性能的影響。Politis等[2]采用邊界元方法計(jì)算了較大參數(shù)空間范圍內(nèi)無粘流場中拍動(dòng)翼的水動(dòng)力性能，并在此基礎(chǔ)上進(jìn)行仿生翼推進(jìn)器的設(shè)計(jì)。Wang等[3]利用FLUENT軟件研究運(yùn)動(dòng)參數(shù)對拍動(dòng)翼水動(dòng)力性能的影響，并對推進(jìn)機(jī)理進(jìn)行簡要探討。文敏華等[4]基于動(dòng)網(wǎng)格技術(shù)分析粘性力和壓差力對拍動(dòng)翼推力的貢獻(xiàn)。但現(xiàn)有研究多側(cè)重于水動(dòng)力性能方面，對粘性流場中三維拍動(dòng)翼的尾渦結(jié)構(gòu)特征及演化機(jī)理的探討還相對較少且不夠深入，而水動(dòng)力性能隨運(yùn)動(dòng)參數(shù)的變化正是由相應(yīng)的尾渦結(jié)構(gòu)差異造成的。因此，系統(tǒng)分析運(yùn)動(dòng)參數(shù)對拍動(dòng)翼尾渦結(jié)構(gòu)的影響，研究尾渦結(jié)構(gòu)與水動(dòng)力性能的關(guān)系具有重要意義，而現(xiàn)階段這方面的研究還是有限的。另一方面，拍動(dòng)翼等仿生流動(dòng)的數(shù)值模擬經(jīng)常涉及復(fù)雜移動(dòng)邊界問題，而這類問題的高效求解一直是計(jì)算流體力學(xué)的難點(diǎn)。浸入邊界法[5]為外形復(fù)雜結(jié)構(gòu)在粘性流場中運(yùn)動(dòng)的數(shù)值模擬提供了新的途徑，它利用拉格朗日方法跟蹤物體運(yùn)動(dòng)邊界，而流場控制方程在固定的笛卡爾網(wǎng)格上求解，不必按照物體形狀生成復(fù)雜的貼體網(wǎng)格，對于移動(dòng)邊界問題也不需要流場的動(dòng)網(wǎng)格更新。文中基于自主開發(fā)的改進(jìn)浸入邊界法數(shù)值求解非定常不可壓縮N-S方程，研究拍動(dòng)翼的尾渦結(jié)構(gòu)和推進(jìn)機(jī)理，并分析運(yùn)動(dòng)參數(shù)對水動(dòng)力性能的影響，探討翼運(yùn)動(dòng)學(xué)、渦動(dòng)力學(xué)和力的產(chǎn)生之間的關(guān)系，理解水動(dòng)力性能隨參數(shù)變化的內(nèi)在機(jī)制。

1拍動(dòng)翼的計(jì)算模型

文中研究拍動(dòng)翼的升沉俯仰耦合運(yùn)動(dòng)，拍動(dòng)翼采用NACA0030翼型剖面，弦長為C，展長為S，置于均勻來流U中，如圖1所示。

圖1　拍動(dòng)翼的計(jì)算模型

拍動(dòng)翼在Y方向的升沉運(yùn)動(dòng)可以表示為

h(t)=h0cos(2πft)

式中：f為運(yùn)動(dòng)頻率，t為時(shí)間，h0為升沉幅度。俯仰運(yùn)動(dòng)是繞距離前緣點(diǎn)1/4弦長位置處的展向軸(Z方向)，可以表示為

式中：θ0為俯仰幅度，ψ為俯仰運(yùn)動(dòng)與升沉運(yùn)動(dòng)間的相位差，這里取為π/2，根據(jù)文獻(xiàn)[6]的研究，此時(shí)拍動(dòng)翼具有最佳的推進(jìn)性能。

在拍動(dòng)翼推進(jìn)中一個(gè)關(guān)鍵的參數(shù)是斯特勞哈爾數(shù)[7]，定義為

St=2h0f/U

雷諾數(shù)Re和翼的展弦比AR分別定義為

Re=UC/n

式中：ν為運(yùn)動(dòng)粘性系數(shù)。

AR=S/C

拍動(dòng)翼的推力系數(shù)為

Cx=Fx/0.5ρU2A

式中： Fx為推力，ρ為流體密度，A為翼的投影面積。一個(gè)運(yùn)動(dòng)周期內(nèi)平均推力系數(shù)可由下式計(jì)算:

式中：T為運(yùn)動(dòng)周期。拍動(dòng)翼的推進(jìn)效率可表示為

式中：Fy為升力，Mz為俯仰力矩。

2數(shù)值方法

2.1控制方程及其離散

流動(dòng)控制方程為三維非定常、不可壓縮N-S方程：

控制方程在固定的笛卡爾網(wǎng)格上離散，采用有限體積法求解，其中u表示單元中心的速度，p為壓力，ρ和μ分別表示流體的密度與動(dòng)力粘度，CV與CS分別表示控制體體積與控制體表面，n表示控制體表面的單位法向量。文中流場數(shù)值計(jì)算基于分裂步方法[8](fractional-step method)。在時(shí)間推進(jìn)方案和空間離散格式上分別采用二階隱式Crank-Nicolson格式和二階中心差分[8]。

2.2浸入邊界的處理

文中用ghost-cell方法[5]施加浸入邊界對流動(dòng)的影響，物體表面用適應(yīng)性較強(qiáng)的三角形網(wǎng)格劃分。

2.2.1網(wǎng)格單元屬性的判斷

首先利用射線法[9]進(jìn)行網(wǎng)格單元內(nèi)外狀態(tài)的判斷，然后可將網(wǎng)格單元分為以下3類(見圖2)：1)流體單元，即單元中心位于物體邊界之外的單元；2)ghost-cell，即位于物體邊界之內(nèi)且至少有一個(gè)相鄰單元位于流體內(nèi)的單元；3)固體單元，即單元中心位于物體邊界之內(nèi)除ghost-cell外的其余單元。

2.2.2ghost-cell公式

下面利用局部流動(dòng)插值建立ghost-cell上應(yīng)滿足的方程，以保證物面邊界條件的正確實(shí)施。在之前的浸入邊界法研究中，盡管各種邊界條件重建算法已經(jīng)被提出，然而對于這些方法的有效性和健壯性還沒有給予足夠關(guān)注。這里提出一個(gè)簡單、快速的重建方案，而且該算法非常健壯，足以應(yīng)對各種情況。

圖2　ghost-cell浸入邊界法示意

先從ghost-cell(G)向物面引垂線，如圖2所示，將垂線與物面的交點(diǎn)O記為邊界點(diǎn)，再將G與O的連線向流體內(nèi)部延伸，得到G關(guān)于浸入邊界在流體域內(nèi)的鏡像點(diǎn)I。在I點(diǎn)周圍局部區(qū)域的流動(dòng)變量φ可由如下的插值多項(xiàng)式表示:

φ=c1+c2x+c3y+c4z

式中：cj為多項(xiàng)式的系數(shù)，可用圍繞I點(diǎn)的立方體的8個(gè)點(diǎn)中的3個(gè)流體點(diǎn)，再加上邊界點(diǎn)O，如圖2所示，總共4個(gè)點(diǎn)作為插值數(shù)據(jù)來計(jì)算，即

其中

其中

多項(xiàng)式的系數(shù)一旦確定，在鏡像點(diǎn)處的φ值即可由下式給出:

(1)

式中：βj是由插值多項(xiàng)式確定的權(quán)系數(shù)。對于速度的Dirichlet邊界條件，根據(jù)中點(diǎn)公式有

(2)

將式(1)代入式(2)，可得:

(3)

對壓力的Neumann邊界條件，利用二階中心差分有

(4)

式中：Δl是法向線段的長度。將式(1)代入式(4)，可得:

(5)

最后ghost-cell上的流動(dòng)插值方程(3)和(5)與流體單元上的N-S方程聯(lián)立求解。在以上過程中，可能會(huì)遇到2種特殊情況。一是圍繞I點(diǎn)的8個(gè)點(diǎn)中有1個(gè)點(diǎn)是ghost-cell自身(如圖2左側(cè)那個(gè)ghost-cell所示)，這不會(huì)引起任何額外的問題，因?yàn)槠渌?個(gè)點(diǎn)中的3個(gè)流體點(diǎn)和邊界點(diǎn)O可作為插值點(diǎn)。二是圍繞I點(diǎn)的8個(gè)點(diǎn)中包含其他的ghost-cell(如圖2右側(cè)那個(gè)ghost-cell所示)，在這種情況下，如果圍繞I點(diǎn)的8個(gè)點(diǎn)中的流體點(diǎn)的數(shù)量N≥3，那么這N個(gè)流體點(diǎn)中的3個(gè)連同邊界點(diǎn)O可作為插值數(shù)據(jù)；若流體點(diǎn)的數(shù)量N<3，則插值模板由N個(gè)流體點(diǎn)、邊界點(diǎn)O和3-N個(gè)其他的ghost-cell構(gòu)成，盡管這確實(shí)暗示著一些ghost-cell上的φ值會(huì)彼此耦合，但并不會(huì)導(dǎo)致任何相容性問題，因?yàn)樵揼host-cell上的流動(dòng)插值方程與流體單元上的N-S方程連同其他ghost-cell上的方程是聯(lián)立在一起求解的。根據(jù)以上分析，文中所提出的邊界條件重建算法足夠健壯，可以處理各種情況。

2.2.3移動(dòng)邊界

在浸入邊界法中，流場控制方程在固定的笛卡爾網(wǎng)格上求解，因此在每一時(shí)間步的計(jì)算完成后，廣闊的外流場并不需要像貼體動(dòng)網(wǎng)格方法那樣從舊網(wǎng)格到新網(wǎng)格的插值，但有些上一時(shí)間步的固體內(nèi)部單元會(huì)由于邊界的移動(dòng)而成為流體單元，這些新出現(xiàn)的流體單元稱之為fresh-cell[10]，如圖3所示，由于它們沒有在流體中的有效時(shí)間歷史，因此將給時(shí)間推進(jìn)求解造成困難。

圖3　邊界移動(dòng)形成Fresh-Cell示意

當(dāng)前浸入邊界法多采用顯式時(shí)間推進(jìn)格式，或?qū)α黜?xiàng)顯式-擴(kuò)散項(xiàng)隱式時(shí)間推進(jìn)格式，為保證穩(wěn)定性，要受到CFL條件的限制，因此一個(gè)時(shí)間步內(nèi)邊界的移動(dòng)不會(huì)超過一層網(wǎng)格深度，Mittal等[10]對于fresh-cell不求解N-S方程，而是采用類似于ghost-cell的局部流動(dòng)插值來得到其流動(dòng)變量，Yang等[11]則提出了field-extension方法。在文中隱式時(shí)間推進(jìn)格式框架下，可選取較大的時(shí)間步以提高計(jì)算效率，而不必?fù)?dān)心穩(wěn)定性問題，另外考慮到物體可能作大幅度運(yùn)動(dòng)，這些情況可能會(huì)使一個(gè)時(shí)間步內(nèi)邊界的移動(dòng)超過一層網(wǎng)格深度，因此需開發(fā)適用于較大步長隱式時(shí)間推進(jìn)和物體作大幅度運(yùn)動(dòng)的fresh-cell數(shù)值計(jì)算方法。這里采用改進(jìn)Shepard插值[12]構(gòu)造解決Fresh-cell的數(shù)值方案，該插值方法獨(dú)立于網(wǎng)格拓?fù)浣Y(jié)構(gòu)，因而具有較強(qiáng)的適應(yīng)性。具體做法是：1)每一時(shí)間步的計(jì)算完成后，將所有的邊界點(diǎn)和ghost-cell作為插值數(shù)據(jù)，構(gòu)成改進(jìn)Shepard插值所需的散點(diǎn)插值集合；2)由改進(jìn)Shepard插值得到所有固體單元上的流動(dòng)變量。這樣在接下來的時(shí)間推進(jìn)求解中，所需的fresh-cell上的流動(dòng)變量及空間導(dǎo)數(shù)都得以解決。

3結(jié)果與討論

3.1數(shù)值方法的驗(yàn)證

為驗(yàn)證數(shù)值方法的正確性，將利用浸入邊界法程序計(jì)算得到的平均推力系數(shù)與Read[6]的實(shí)驗(yàn)結(jié)果及采用面元法得到的計(jì)算結(jié)果進(jìn)行比較，如圖4所示，關(guān)于面元法的詳細(xì)介紹可參考李寧宇等[13-14]之前的研究工作。

計(jì)算中拍動(dòng)翼的參數(shù)為AR=3，Re=1 000，St=0.2-0.5，h0=0.5C，θ0=5°?？梢钥闯?，浸入邊界法和面元法得出的Cxm變化趨勢與實(shí)驗(yàn)結(jié)果一致，Cxm隨St的增加而增大。由于考慮了流體粘性，相比于面元法，浸入邊界法可以更精確地計(jì)算拍動(dòng)翼在真實(shí)流動(dòng)環(huán)境中的水動(dòng)力。

圖4　平均推力系數(shù)計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對比

3.2尾渦結(jié)構(gòu)

本節(jié)研究拍動(dòng)翼的尾渦結(jié)構(gòu)特征及演化機(jī)理，數(shù)值模擬中三維空間渦結(jié)構(gòu)的識(shí)別采用Q判據(jù)[15]，Q為速度梯度張量的第二不變量，在Q>0的區(qū)域內(nèi)流體的旋轉(zhuǎn)比較應(yīng)變而言起主導(dǎo)作用。當(dāng)AR=1.5，Re=200，St=0.5，h0=0.5C，θ0=30°時(shí)，拍動(dòng)翼尾渦的透視圖、側(cè)視圖和俯視圖分別如圖5(a)～(c)所示，圖中根據(jù)計(jì)算出的渦矢量用箭頭標(biāo)出了渦的方向。結(jié)果表明尾渦由2列形狀復(fù)雜的渦環(huán)組成，這2列渦環(huán)與尾流中心線成一定傾角向下游對流，這是三維拍動(dòng)翼尾渦結(jié)構(gòu)的重要特征。在圖中標(biāo)出了處于上面一列渦環(huán)之中的環(huán)R1與位于下面一列渦環(huán)之中的環(huán)R2?？梢园l(fā)現(xiàn)，在向下游對流的過程中，尾渦在Z方向稍有變窄，在Y方向逐漸變寬，此現(xiàn)象的機(jī)理可通過研究圖5(a)～(c)中虛線所標(biāo)出的區(qū)域內(nèi)流向平面上的流動(dòng)加以解釋。圖5(d)為在此平面內(nèi)梢渦(TV1與TV2)和展向渦(TV2)的示意圖，其中用實(shí)心圓表示的梢渦垂直紙面向外，用空心圓表示的梢渦垂直紙面向里，其他3個(gè)梢渦在1個(gè)梢渦上的誘導(dǎo)速度方向用箭頭標(biāo)出。

從誘導(dǎo)速度方向可以看出，同處上面渦環(huán)中的2個(gè)梢渦TV1彼此互相靠近，同處下面渦環(huán)中的2個(gè)梢渦TV2也是如此，從而導(dǎo)致尾渦在Z方向稍有變窄；處于不同渦環(huán)中的梢渦TV1和TV2趨于在垂向遠(yuǎn)離，從而造成尾渦在Y方向逐漸變寬。梢渦誘導(dǎo)的渦V1向上運(yùn)動(dòng)是導(dǎo)致渦環(huán)傾斜的主要機(jī)理，而一旦渦環(huán)傾斜，它的自誘導(dǎo)速度趨于使渦環(huán)沿著它的軸線向下游對流，這是2列渦環(huán)與尾流中心線成一定傾角的原因。由于渦環(huán)相對于來流方向是傾斜的，因此渦環(huán)沿其軸向所誘導(dǎo)的流動(dòng)有一個(gè)流向成分和一個(gè)垂向成分，射流的流向成分與拍動(dòng)翼的推力產(chǎn)生直接相關(guān)。

圖5　尾渦結(jié)構(gòu)與誘導(dǎo)機(jī)理

3.3St對尾渦結(jié)構(gòu)和水動(dòng)力性能的影響

本節(jié)進(jìn)行斯特勞哈爾數(shù)的影響分析。數(shù)值模擬中St的變化范圍是0.35～1.1，其他參數(shù)為AR=1.5，Re=200，h0=0.5C，θ0=30°。圖6給出了St=0.35和St=0.95時(shí)拍動(dòng)翼的尾渦結(jié)構(gòu)，可結(jié)合圖5中St=0.5時(shí)的情況來分析St對尾渦結(jié)構(gòu)的影響?？梢园l(fā)現(xiàn)，當(dāng)St=0.35時(shí)，一個(gè)明顯的特征是St=0.5時(shí)所看到的梢渦之間的某些連接消失了，該現(xiàn)象是由于在低St下所形成的梢渦有更低的強(qiáng)度，這減弱了3.2節(jié)所闡述的梢渦之間相互誘導(dǎo)的機(jī)理。當(dāng)St增加至0.5時(shí)，由拍動(dòng)翼左右兩側(cè)產(chǎn)生的梢渦向中間靠近，相互融合， 形成渦環(huán)，并且可以看到連接2個(gè)渦列的發(fā)卡渦出現(xiàn)。隨著St進(jìn)一步增加到0.95，渦環(huán)傾角的改變使渦環(huán)的軸向與流向更加接近，并且2列渦環(huán)之間發(fā)展出更復(fù)雜的相互連接，顯然隨著St的增加，梢渦的強(qiáng)度增大，梢渦之間更強(qiáng)的相互誘導(dǎo)作用導(dǎo)致尾流在垂向變得更寬。

圖6　尾渦結(jié)構(gòu)隨St的變化

平均推力系數(shù)和推進(jìn)效率隨St的變化如圖7所示，可以看到推力隨St的增加而增大，這是因?yàn)殡SSt的增加，梢渦的強(qiáng)度增大，并且渦環(huán)傾角的改變使渦環(huán)的軸向與流向更加接近，從而射流的流向動(dòng)量增加，根據(jù)動(dòng)量定理，拍動(dòng)翼產(chǎn)生的推力也相應(yīng)增大。而推進(jìn)效率先是隨St的增加而增大，這顯然是由推力的增大造成的，在St=0.8左右，效率達(dá)到峰值，約為19%，之后又逐漸減小，這是因?yàn)殡SSt的增加，兩列渦環(huán)之間產(chǎn)生更多的相互連接，下游的尾渦結(jié)構(gòu)變得更加復(fù)雜，導(dǎo)致相反符號(hào)渦量之間的粘性抵消作用加強(qiáng)。注意文中采用基于浸入邊界法的直接數(shù)值模擬探討低雷諾數(shù)、低展弦比拍動(dòng)翼的尾渦結(jié)構(gòu)和水動(dòng)力性能，故推進(jìn)效率低于李寧宇等[13-14]之前有關(guān)高雷諾數(shù)、高展弦比拍動(dòng)翼[15]的研究工作；根據(jù)Dong等[16]對低雷諾數(shù)、有限展弦比橢圓翼的研究，當(dāng)AR=1.27和AR=2.55時(shí)得到橢圓翼的峰值效率分別約為18%和21%。

圖7　水動(dòng)力性能隨斯特勞哈爾數(shù)的變化

3.4θ0對尾渦結(jié)構(gòu)和水動(dòng)力性能的影響

本節(jié)研究俯仰幅度的影響。計(jì)算參數(shù)為AR=1.5，Re=200，St=0.5，h0=0.5C，θ0=10°～40°。圖8給出了θ0=10°和θ0=40°時(shí)拍動(dòng)翼的尾渦結(jié)構(gòu)，可結(jié)合圖5中θ0=30°時(shí)的情況來研究θ0對尾渦結(jié)構(gòu)的影響?？梢钥闯觯Sθ0的增加，渦環(huán)傾角的改變使渦環(huán)的軸向與流向更加接近，這使得射流的流向動(dòng)量增加；另一方面，渦環(huán)耗散速度加快，以致θ0=40°時(shí)遠(yuǎn)下游的渦環(huán)已經(jīng)由于粘性耗散作用而消失。

圖8　尾渦結(jié)構(gòu)隨俯仰幅度的變化

圖9為平均推力系數(shù)和推進(jìn)效率隨θ0的變化，可以發(fā)現(xiàn)，推力先是隨θ0的增加而增大，這是由于射流的流向動(dòng)量增加，在θ0=20°附近推力達(dá)到最大值，而后因?yàn)槲擦髦袦u環(huán)的耗散速度增加，推力又開始減小。與推力的變化情況類似，推進(jìn)效率先是隨θ0的增加而上升，這顯然是因?yàn)橥屏Φ脑龃螅?dāng)θ0為20°～30°某一值時(shí)，效率取得最大值，之后因推力減小和渦環(huán)耗散加快，效率又逐漸下降。

圖9　水動(dòng)力性能隨俯仰幅度的變化

4結(jié)論

文中應(yīng)用自主開發(fā)的改進(jìn)浸入邊界法數(shù)值模擬拍動(dòng)翼的運(yùn)動(dòng)，分析其尾渦結(jié)構(gòu)特征和演化機(jī)理，并探討運(yùn)動(dòng)參數(shù)對尾渦結(jié)構(gòu)和水動(dòng)力性能的影響，得到如下結(jié)論：

1) 所提出的邊界條件重建算法兼顧計(jì)算效率和健壯性，可保證物面邊界條件的正確實(shí)施。

2) 基于改進(jìn)Shepard插值構(gòu)造的移動(dòng)邊界處理方法可有效解決浸入邊界法中與移動(dòng)邊界有關(guān)的fresh-cell問題，并且該處理方法在選取較大時(shí)間步或物體作大幅度運(yùn)動(dòng)時(shí)也完全適用。

3) 實(shí)驗(yàn)結(jié)果和其他計(jì)算結(jié)果的比較表明，所開發(fā)的改進(jìn)浸入邊界法程序可很好地完成三維拍動(dòng)翼在粘性流場中非定常運(yùn)動(dòng)的數(shù)值模擬。

4) 三維拍動(dòng)翼的尾渦由2列形狀復(fù)雜的渦環(huán)組成，這2列渦環(huán)與尾流中心線成一定傾角向下游對流，且尾渦在展向稍有變窄，在垂向逐漸變寬。在渦環(huán)軸向誘導(dǎo)速度作用下，翼的下游形成2股傾斜射流，射流的流向動(dòng)量與推力的產(chǎn)生直接相關(guān)。

5) 隨St增加，梢渦強(qiáng)度增加，2列渦環(huán)之間產(chǎn)生更多的相互連接，且渦環(huán)傾角的改變使渦環(huán)的軸向與流向更加接近，尾流在垂向變寬。推力隨St增加而增大，存在一個(gè)最優(yōu)St使推進(jìn)效率最高。

6) 隨θ0的增加，渦環(huán)傾角的改變使渦環(huán)的軸向與流向更加接近，且渦環(huán)耗散速度加快。推力和推進(jìn)效率均隨θ0的增加而先增大后減小，存在一個(gè)最優(yōu)θ0使推進(jìn)性能最佳。

參考文獻(xiàn):

[1]MACKOWSKIAW,WILLIAMSONCHK.Directmeasurementofthrustandefficiencyofanairfoilundergoingpurepitching[J].JournalofFluidMechanics, 2015(765): 524-547.

[2]POLITISGK,TSARSITALIDISVT.Flappingwingpropulsordesign:anapproachbasedonsystematic3D-BEMsimulations[J].OceanEngineering, 2014, 84(7): 98-123.

[3]SuYM,WangZL,ZhangX,etal.Numericalsimulationforhydrodynamicnumericalcharacteristicsofabionicflappinghydrofoil[J].ChinaOceanEngineering, 2012, 26(2): 291-304.

[4]文敏華，胡文蓉，劉洪.翼沉浮運(yùn)動(dòng)推力來源的數(shù)值研究[J].水動(dòng)力學(xué)研究與進(jìn)展，2012, 27(2):154-161.

[5]MITTALR,IACCARINOG.Immersedboundarymethods[J]AnnualReviewofFluidMechanicS, 2005(37): 239-261.

[6]READMB.Performanceofbiologicallyinspiredflappingfoils[D].Cambridge:MassachusettsInstituteofTechnology, 2006:133-135.

[7]PAN Y L, DONG X X, ZHU Q, et al. Boundary-element method for the prediction of performance of flapping foils with leading-edge separation[J]. Journal of Fluid Mechanics, 2012(698) : 446-467.

[8]FERZIGER J H，PERIC M. Computational methods for fluid dynamics[M]. New York: Springer-Verlag, 2002: 56-67.

[9]BORAZJANI I. Numerical simulations of fluid-structure interaction problems in biological flows[D]. Minneapolis, The Faculty of the Graduate School of the University of Minnesota,2008:1-273.

[10]MITTAL R, DONG H, BOZKURTTAS M, et al. A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[J]. Journal of Computational Physics, 2008, 227(10): 4825-4852.

[11]YANG J M, BALARAS E. An embedded-boundary formulation for large-eddy simulation of turbulent flows interacting with moving boundaries[J]. Journal of Computational Physics, 2006, 215(1):12-40.

[12]THACKER W I, ZHANG J, WATSON L T, et al. Algorithm 905: SHEPPACK: modified Shepard algorithm for interpolation of scattered multivariate data[J]. ACM Transactions on Mathematical Software, 2010(37):1-20.

[13]李寧宇，蘇玉民，王兆立，等，三維拍動(dòng)翼推進(jìn)效率的計(jì)算方法[J].上海交通大學(xué)學(xué)報(bào)，2012, 46(2):323-328.

[14]李寧宇.基于面元法的拍動(dòng)翼水動(dòng)力性能計(jì)算[D].哈爾濱:哈爾濱工程大學(xué)，2012:1-71.

[15]HUNT J C R, WRAY A A. Moin P, Eddies, stream, and convergence zones in turbulent flows, Center for Turbulence Research Report No. CTR-S88[R], Center for Turbulence Research, Stanford, USA, 1988:193-208.

[16]DONG H, MITTAL R, BOZKURTTAS M. Wake structure and performance of finite aspect-ratio flapping foils[C]. //43rd AIAA Aerospace Sciences Meeting, Reno, Nevada, 2005:1-9.

網(wǎng)絡(luò)出版地址：http://www.cnki.net/kcms/detail/23.1191.U.20151205.2117.006.html

Hydrodynamic analysis and trailing vortex structure of flapping

wings based on the immersed boundary method

LI Ningyu，SU Yumin，LIU Weixing，LIU Peng，CAO Jian

Science and Technology on Underwater Vehicle Laboratory, Harbin Engineering University, Harbin 150001, China

Abstract：At present, the number of studies that have examined the relationship between flapping wing trailing vortex structure and hydrodynamic performance is limited, and the analysis of the effect of parameters is focused on hydrodynamic performance. Parameter simulation for flapping wing and other bionic flows always involves in issues of the complex moving boundary. Self-developed improved immersed boundary method is applied to numerically simulate the non-stationary motion of the flapping wing. Relationship between wing kinematics, vortex dynamics and generation of force is discussed. Results show that the reconstruction algorithm of boundary conditions gives consideration to both calculation efficiency and robustness. The treatment method of the moving boundary developed in this paper can effectively solve the issues related to the moving boundary of the immersed boundary method. Trailing vortex of the 3D flapping wing is constituted by two lines of complex vortex rings. The two lines of complex vortex rings and center line of trailing flow form a dip angle and flow to downstream. Hydrodynamic performance varying with kinematic parameters is decided by the vortex dynamics under the appearance of force production, including the strength and direction of vortex, interconnection of vortex rings and viscous dissipation.

Keywords：flapping wing; immersed boundary method; trailing vortex structure; hydrodynamic performance; propulsion mechanism; incompressible N-S equation

通信作者：李寧宇，E-mail：lny123d@163.com.

作者簡介：李寧宇(1986-), 男，博士研究生；

基金項(xiàng)目：國家自然科學(xué)基金資助項(xiàng)目(51479039).

收稿日期：2015-08-10.網(wǎng)絡(luò)出版日期：2015-12-05.

中圖分類號(hào)：O357.1

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1009-671X(2015)06-026-06

doi:10.11991/yykj.201508009