基于浸入邊界法的拍動翼尾渦結構與水動力分析

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

哈爾濱工程大學 水下機器人技術重點實驗室，黑龍江 哈爾濱 150001

?

基于浸入邊界法的拍動翼尾渦結構與水動力分析

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

哈爾濱工程大學 水下機器人技術重點實驗室，黑龍江 哈爾濱 150001

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

關鍵詞：拍動翼；浸入邊界法；尾渦結構；水動力性能；推進機理；不可壓縮N-S方程

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

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

1拍動翼的計算模型

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

圖1　拍動翼的計算模型

拍動翼在Y方向的升沉運動可以表示為

h(t)=h0cos(2πft)

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

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

在拍動翼推進中一個關鍵的參數(shù)是斯特勞哈爾數(shù)[7]，定義為

St=2h0f/U

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

Re=UC/n

式中：ν為運動粘性系數(shù)。

AR=S/C

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

Cx=Fx/0.5ρU2A

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

式中：T為運動周期。拍動翼的推進效率可表示為

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

2數(shù)值方法

2.1控制方程及其離散

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

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

2.2浸入邊界的處理

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

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

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

2.2.2ghost-cell公式

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

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

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

φ=c1+c2x+c3y+c4z

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

其中

其中

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

(1)

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

(2)

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

(3)

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

(4)

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

(5)

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

2.2.3移動邊界

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

圖3　邊界移動形成Fresh-Cell示意

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

3結果與討論

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

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

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

圖4　平均推力系數(shù)計算結果與實驗數(shù)據(jù)對比

3.2尾渦結構

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

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

圖5　尾渦結構與誘導機理

3.3St對尾渦結構和水動力性能的影響

本節(jié)進行斯特勞哈爾數(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時拍動翼的尾渦結構，可結合圖5中St=0.5時的情況來分析St對尾渦結構的影響。可以發(fā)現(xiàn)，當St=0.35時，一個明顯的特征是St=0.5時所看到的梢渦之間的某些連接消失了，該現(xiàn)象是由于在低St下所形成的梢渦有更低的強度，這減弱了3.2節(jié)所闡述的梢渦之間相互誘導的機理。當St增加至0.5時，由拍動翼左右兩側產(chǎn)生的梢渦向中間靠近，相互融合， 形成渦環(huán)，并且可以看到連接2個渦列的發(fā)卡渦出現(xiàn)。隨著St進一步增加到0.95，渦環(huán)傾角的改變使渦環(huán)的軸向與流向更加接近，并且2列渦環(huán)之間發(fā)展出更復雜的相互連接，顯然隨著St的增加，梢渦的強度增大，梢渦之間更強的相互誘導作用導致尾流在垂向變得更寬。

圖6　尾渦結構隨St的變化

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

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

3.4θ0對尾渦結構和水動力性能的影響

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

圖8　尾渦結構隨俯仰幅度的變化

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

圖9　水動力性能隨俯仰幅度的變化

4結論

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

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

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

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

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

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

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

參考文獻:

[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]文敏華，胡文蓉，劉洪.翼沉浮運動推力來源的數(shù)值研究[J].水動力學研究與進展，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]李寧宇，蘇玉民，王兆立，等，三維拍動翼推進效率的計算方法[J].上海交通大學學報，2012, 46(2):323-328.

[14]李寧宇.基于面元法的拍動翼水動力性能計算[D].哈爾濱:哈爾濱工程大學，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)絡出版地址：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-), 男，博士研究生；

基金項目：國家自然科學基金資助項目(51479039).

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

中圖分類號：O357.1

文獻標志碼：A

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

doi:10.11991/yykj.201508009