周軍偉，梅 蕾，郭 斌，李明陽，潘小云

（哈爾濱工業(yè)大學(xué)（威海）海洋工程學(xué)院，山東威海 264209）

0 引 言

自然游泳者諸如魚類、水生哺乳動物等已經(jīng)進化得相當(dāng)完美，其神秘的游動機理帶來了低噪聲、高機動性和超高效率的游動[1]。近幾十年來，越來越多的科研工作者致力于揭開仿生推進的神秘“面紗”，不斷解析生物游動中暗藏的物理規(guī)律，并逐步應(yīng)用于工程實踐中。從產(chǎn)生推力的水動力特征來劃分，仿生推進大致可以分為三類，分別是噴射法，如烏賊等；擺動身體法，如鰻魚等；擺動翼法，如海龜?shù)?。其中，擺動翼法可同時執(zhí)行平移和旋轉(zhuǎn)運動，形成類似拍打的運動，在拍打過程中，渦流可以被卷起到尾流區(qū)域形成反向的von-Karman渦流道，以產(chǎn)生推力[2]。

由于擺動翼法運動形式自由度少，推進效率高，可達80%[3-6]，且便于工程實現(xiàn)，因而研究者們展開了大量的研究，包括自然觀察、理論分析、實驗研究和計算研究。早期的研究主要通過自然觀察進行的[7-9]，研究指出，游泳動物試圖維持某些參數(shù)以達到最佳效率。為了找出這些參數(shù)，已經(jīng)進行了許多有關(guān)擺翼推進的實驗和數(shù)值模擬研究。Lighthill[1]較早就發(fā)展了一種二維擺動翼的線性理論方法。隨后，在1995 年Streitlien 等[10]對該方法進行了改進，使其能夠計算大擺幅運動情況。而后，Anderson等[3]對比了Lighthill 的線性方法與Streitlien 的非線性方法在求解擺翼問題上的結(jié)果差異，并與自己的實驗測量結(jié)果進行了對比。而Schouveiler 等[4]采用與Anderson 同樣的方法進行了擺翼性能的實驗測量和數(shù)值計算，并對擺翼的工作機理進行了更加深入的研究。從以上研究者的數(shù)值與實驗的對比研究來看，基于勢流理論的數(shù)值方法與實驗結(jié)果存在一定的誤差。1996 年，Tuncer 等[11]引入粘流理論，采用求解二維粘流方程的方法得到了一組串列擺翼的性能。而后，Jones 等[12]采用與Tuncer 相同的粘流方程計算方法計算了某擺翼性能，并進行了實驗測量。但從結(jié)果對比來看，計算結(jié)果與實驗結(jié)果相差仍明顯。

為了更準(zhǔn)確地模擬擺翼的大擺幅和非定常運動特征，許多新的方法得以發(fā)展。Smith 等[13]發(fā)展了一種基于彈性網(wǎng)格的RANS求解器，使其能夠方便地求解三維擺動翼的水動力性能。Akhtar等[14]發(fā)展了一種浸入式邊界方法，以適應(yīng)擺翼的大幅度運動，并對一對串列的擺動平板的性能進行了研究。張曦等[6]以一個魚尾型的三維擺翼為研究對象，將面元法的計算結(jié)果與軟件Fluent的計算結(jié)果和實驗結(jié)果進行了對比，結(jié)果吻合較好。陳剛等[15]利用浸入邊界法和格子Boltzmann 方法來研究幾何參數(shù)對三維仿生運動擺翼推進性能的影響。Karbasiana 等[16]采用開源CFD 軟件OpenForm 的動網(wǎng)格功能實現(xiàn)了二維擺翼的水動力模擬。Mannam 等[17]研究了企鵝模式擺翼的推進性能，并將Fluent 計算結(jié)果與實驗結(jié)果進行了對比，吻合也較好。

出于繼續(xù)豐富擺翼研究的實驗數(shù)據(jù)和數(shù)值方法的目的，本文設(shè)計了一種簡易的兩自由度擺翼試驗平臺，實現(xiàn)了預(yù)定的兩自由度運動軌跡，并在循環(huán)水槽開展了擺動平板推進性能實驗研究；隨后基于動網(wǎng)格技術(shù)，利用流體數(shù)值仿真軟件對擺動平板推進性能進行了數(shù)值仿真分析，并與實驗結(jié)果進行了對比分析。最后，針對實驗中擺動平板效率較低的情況，數(shù)值分析了平板前后緣導(dǎo)圓、無限展長假設(shè)和大擺幅對效率的影響，以期待可以進一步地對試驗平臺進行改進。

1 擺翼實驗平臺設(shè)計

1.1 擺翼運動機理

簡單的擺翼運動可描述為勻速前進（X方向平移）、橫蕩（Y方向平移）和首搖（繞Z軸旋轉(zhuǎn)）三自由度耦合運動。擺翼運動使用絕對坐標(biāo)系xoy，其運動模式如圖1所示。圖中以某翼型代替本研究中的平板，VA為擺翼的恒定前進速度，T為擺動周期，B為擺動幅度，即掃掠寬度。其中，橫蕩運動和首搖均為正弦規(guī)律的運動，橫蕩位移Y和首搖角度θ分別表示為

圖1 擺翼運動和受力示意圖Fig.1 Schematic diagram of flapping foil motion and force

式中，Y0是橫蕩幅度（m），θ0是首搖幅度（rad），f是頻率（Hz），φ是橫蕩運動和首搖運動之間的相位差（rad），t是時間（s）。

1.2 性能參數(shù)定義

對比螺旋槳的進速系數(shù)定義，本文選用橫蕩掃掠寬度B作為特征長度，定義進程hp與橫蕩掃掠寬度B的比值為進速系數(shù)，以J表示：

由式（3）可知，在頻率f、掃掠寬度B一定的條件下，改變VA值即可得到不同的進速系數(shù)J，本文擺翼運動的不同工況即指不同進速系數(shù)的工況。許多文獻中，采用斯特勞哈爾數(shù)St代替進速系數(shù)來描述擺翼的不同工況，兩者關(guān)系表示為

參照螺旋槳性能推力系數(shù)的定義以及Floc’h等[5]的定義方法，以擺翼掃掠面積BH為參考面積，以fB為參考速度，其中H為擺翼的展長，重新定義推力系數(shù)KT為

按相對速度計算，擺翼推進器的工作Re數(shù)可表示為

式中，ν為水的運動粘度系數(shù)。本文各工況雷諾數(shù)Re約在1.6 × 105～5.0 × 105的范圍內(nèi)。

1.3 擺翼實驗平臺整體設(shè)計

由于擺翼的運動為三自由度耦合運動，一些研究者采用多個伺服電機來實現(xiàn)上述結(jié)構(gòu)[3-4,18], 但本文實驗中發(fā)現(xiàn)該結(jié)構(gòu)的脈動力明顯，為了減輕脈動，實現(xiàn)擺翼連續(xù)平滑的運動，并簡化試驗機構(gòu)，本設(shè)計參照劉江鹓等的工作[19]，采用曲柄連桿結(jié)構(gòu)和直流電機相結(jié)合的方式控制擺翼運動，如圖2 所示。本文采用曲柄連桿滑軌機構(gòu)來實現(xiàn)往復(fù)的橫蕩運動，采用第二套曲柄連桿滑軌來控制以實現(xiàn)一定的擺角θ。其中兩滑塊分別連接以同轉(zhuǎn)速異相位圓周運動的曲柄，以保證θ值呈近似正弦變化規(guī)律，兩曲柄間的相位差可控制θ的幅值大小，且可以保證橫蕩和首搖運動間的相位差φ為90°。水平擺幅的大小可通過調(diào)整曲柄長度實現(xiàn)。圖2中紅色鉸鏈繞X軸旋轉(zhuǎn)，藍(lán)色鉸鏈繞Z軸旋轉(zhuǎn)。

圖2 擺翼機構(gòu)運動原理示意圖Fig.2 Schematic diagram of motion principle of flapping foil mechanism

如圖3實物圖和裝配圖所示，擺翼機構(gòu)共分為四大部分：電機驅(qū)動部分、曲柄連桿部分、滑塊軌道部分和擺翼測量部分。電機驅(qū)動部分為整套裝置提供動力并通過聯(lián)軸器將電機與曲柄軸相連接，帶動曲柄做圓周運動；兩曲柄間也與聯(lián)軸器固接，以實現(xiàn)同轉(zhuǎn)速轉(zhuǎn)動；曲柄連桿帶動滑塊做水平運動；滑塊與擺翼通過鉸鏈鏈接，帶動擺翼作預(yù)定軌跡運動；測力裝置固連在擺翼上方，這里設(shè)計采用三分力天平測量擺翼受力。試驗中用平板代替擺翼，弦長c=0.1 m，豎直浸入水中，浸沒展長H=0.45 m。平板截面為矩形，厚度為5 mm，即弦長厚度比為20。平板實驗件與自制三分力天平如圖4所示。

圖3 曲柄連桿滑軌機構(gòu)的實物圖和示意圖Fig.3 Picture（left）and schematic diagram（right）of crank connecting rod slide rail mechanism

圖4 平板實驗件（左）與自制三分力測力計（右）Fig.4 Plate test specimen（left）and self-made three axis force sensor（right）

由于曲柄連桿方式不能得到標(biāo)準(zhǔn)的正弦曲線，因而必須通過各部分的幾何關(guān)系計算出擺翼運動中每一時刻的首搖角θ（t）和水平位置Y（t）。換算方法為：由擺翼曲柄的轉(zhuǎn)動方程，結(jié)合曲柄長度和連桿長度推出其中一曲柄控制的滑塊水平位置時變曲線，進而通過兩曲柄間的相位差可求出另一滑塊的水平位置時變曲線，兩水平位置中心便是擺翼的型心位置；通過兩個滑塊的位置、擺翼的弦長和鉸鏈的長度即可求得首搖角度θ（t）。具體的換算和推導(dǎo)過程詳見附錄A。設(shè)計電機轉(zhuǎn)速為1 r/s，實際實驗中測得電機轉(zhuǎn)速為1.011 r/s，即擺動頻率f=1.011 Hz，周期約為T=0.989 s。兩組曲柄設(shè)計長度均為80 mm，連桿長度都為200 mm，當(dāng)驅(qū)動電機勻速轉(zhuǎn)動，可得A點的橫蕩運動曲線y，如圖5（a）所示，橫蕩幅度為Y0=80 mm。圖2 中黑色L 型支架的長邊S1為148 mm，短邊S2為80 mm。定義兩個曲柄間相位差為Δβ，當(dāng)Δβ分別取0.4 rad 和0.6 rad 時，可得平板的首搖運動曲線θ（t），如圖5（b）所示，首搖擺角幅度分別約為0.23 rad和0.35 rad。從圖中可以明顯看出，兩個運動的曲線都不是簡諧的，即兩個運動的規(guī)律均未達到絕對的正弦曲線形式。由理論分析可知，增加連桿長度和S2長度可以改善這一現(xiàn)象，但為了減小裝置的整體尺寸，本文中并未嘗試。

圖5 擺動平板的橫蕩運動曲線和兩個不同相位差下的首搖運動曲線Fig.5 Heave motion and pitching angles under two different phase differences of flapping plate

1.4 測量和數(shù)據(jù)采集

實驗采用三分力天平作為測力計測量推力Fx、側(cè)向力Fy以及扭矩Mz，測力計固連在圖2中S1支架的中點位置，同時與平板固連，其示意圖如圖3（b）所示。當(dāng)平板運動時，實驗測量的力的矢量方向是隨時間變化的，因而實驗中實測的力為沿平板弦向的力Fs、垂直平板方向的力Fn，以及繞Z軸的轉(zhuǎn)矩Mz，坐標(biāo)系如圖2 所示，F(xiàn)s、Fn與Fx、Fy的關(guān)系如圖6 所示。已知測得的Fs和Fn，根據(jù)如下關(guān)系可得出Fx和Fy：

圖6 平板受力示意圖Fig.6 Sketch of forces on flapping foil

實驗裝置在循環(huán)水槽中對擺動平板的受力進行測量，測量段長8 m，寬2 m，深1.2 m。擺動平板豎直浸入水中，通過控制來流速度VA實現(xiàn)不同工況的測量。采用東華DH5908 無線動態(tài)應(yīng)變采集器進行數(shù)據(jù)采集，采樣頻率為200 Hz，即采樣周期為0.005 s。為了去除測量結(jié)果中的雜波，采用FFT 低通濾波對測量數(shù)據(jù)進行處理。考慮測量信號中主要關(guān)注的是低頻分量，且實驗用循環(huán)水槽動力設(shè)備振動頻率為3～6 Hz，因而本文中截止頻率選用2.5 Hz，以消除設(shè)備振動對測量結(jié)果的影響。圖7 列出了兩個分別代表雜波較?。‵n的應(yīng)變測量數(shù)據(jù)）和雜波較大（Fs的應(yīng)變測量數(shù)據(jù)）的濾波情況，圖中工況都為Δβ=0.6 rad，來流速度為0.6 m/s，黑線為原始測量數(shù)據(jù)，紅線為濾波后數(shù)據(jù)。從濾波結(jié)果來看，圖7（b）中Fn應(yīng)變測量數(shù)據(jù)在濾波前后趨勢幾乎一致，僅濾掉了較小幅度的雜波；圖7（a）中Fs的應(yīng)變結(jié)果相差較大，部分幅度明顯的高頻被濾掉了，但考慮到高頻分量對擺動平板的時均性能影響較小，因而本文采用2.5 Hz的低通濾波截止頻率是可行的。

圖7 Δβ=0.6 rad，VA=0.6 m/s工況下弦向應(yīng)變與法向應(yīng)變的測量數(shù)據(jù)及濾波結(jié)果Fig.7 Measured data and filtering results of tangential strain and normal strain under Δβ=0.6 rad，VA=0.6 m/s working condition

2 擺翼推進性能數(shù)值分析

本文采用NUMECA 公司的CFD 軟件包FINE/Marine 求解擺翼平板周圍流場。該軟件采用動網(wǎng)格技術(shù)實現(xiàn)對運動物體周圍流場的模擬，采用高分辨率格式能夠很好地捕捉到自由液面，能夠滿足擺翼大幅度運動的流場的模擬要求。

2.1 控制方程

求解積分型不可壓粘性流體動力學(xué)方程，考慮網(wǎng)格單元的運動和重力的影響，連續(xù)方程和動量方程分別為

式中：Ω表示單元體積；Sj表示單元面積矢量S的各個分量；vi表示速度矢量v的各個分速度；τij為粘性應(yīng)力張量的各個分量；p為壓強；gi為單元的加速度分量，其僅在重力方向有分量，大小為重力加速度g。算例采用的湍流模型均為k-ω（SST）模型，湍流模型的具體參數(shù)參見Menter[20]的工作，該模型考慮了湍流切應(yīng)力傳遞的影響，在逆壓梯度和分離流中表現(xiàn)優(yōu)異。

2.2 網(wǎng)格劃分及無關(guān)性驗證

2.2.1 網(wǎng)格劃分

利用NUMECA 系列軟件中的HEXPRESS 模塊生成矩形體網(wǎng)格，逐次進行網(wǎng)格加密。首先在平板表面進行網(wǎng)格加密，使最大網(wǎng)格尺度小于1/64 弦長；在翼型尾流區(qū)域進行兩次網(wǎng)格加密，以保證尾渦的較精確模擬，其中在2c×10c×h的區(qū)域（加密區(qū)1），最大網(wǎng)格尺寸小于1/16弦長，在10c×20c×h的區(qū)域（加密區(qū)2），最大網(wǎng)格尺寸小于1/4 弦長，h為平板高度；在自由液面附近亦進行了網(wǎng)格加密。整個計算域大小約為25c×25c×50c，網(wǎng)格剖視如圖8 所示。計算中采用網(wǎng)格剛性平移+彈性扭曲的策略，以適應(yīng)平板的大幅度運動。在本文中，平板僅在Z軸方向發(fā)生一定角度的偏轉(zhuǎn)，隨平板偏轉(zhuǎn)而發(fā)生變形的網(wǎng)格示例如圖9所示。

圖8 網(wǎng)格剖示圖Fig.8 Schematic diagram of grid section

圖9 水下0.2米處平板在發(fā)生首搖時網(wǎng)格所發(fā)生的局部形變示例Fig.9 Example of grid local deformation in case of pitching motion of foil at 0.2 m under water

2.2.2 邊界條件

算例邊界條件如圖10 所示，擺翼壁面為無滑移壁面條件，擺翼前方為遠(yuǎn)場來流條件（來流速度為進速VA，壓強為當(dāng)?shù)仂o壓），左、右側(cè)和后方均為零壓梯度條件（第二類Neumann邊界條件，即?p/?n=0，速度為遠(yuǎn)場來流速度），上方和下方為描述壓力條件（第一類Dirichlet 邊界條件，壓強為給定當(dāng)?shù)仂o壓，速度為遠(yuǎn)場來流速度），V→mir×n→=0。

圖10 邊界條件示意圖Fig.10 Schematic diagram of boundary condition

2.2.3 網(wǎng)格無關(guān)性驗證根據(jù)課題組前期研究結(jié)果，擺翼尾流區(qū)域的網(wǎng)格大小對性能曲線影響較小，擺翼周圍網(wǎng)格大小對性能曲線的影響遠(yuǎn)低于壁面網(wǎng)格大小，可在資源充沛時適當(dāng)加密。因此本文中驗證網(wǎng)格的無關(guān)性時僅對比了平板壁面三種加密形式對性能的影響。網(wǎng)格細(xì)分方法如表1所示，三種網(wǎng)格計算得到的效率結(jié)果如圖11所示。驗證工作的首搖幅度均取為0.23 rad，湍流模型選用的是k-ω模型，時間步長為0.002 s。由圖11可知，網(wǎng)格1和網(wǎng)格2效率曲線差距較大，即棱邊加密對擺動平板效率的影響較明顯；而網(wǎng)格2與網(wǎng)格3 的效率曲線幾乎重合，即繼續(xù)細(xì)分平板表面網(wǎng)格對效率幾乎沒有影響。因此后文中將選用網(wǎng)格2來展開討論。

圖11 三種網(wǎng)格對效率的影響Fig.11 Performance curves of three different mesh sizes

表1 網(wǎng)格劃分情況Tab.1 Details of mesh refinement

2.3 數(shù)值方法驗證

為驗證文中數(shù)值方法在模擬擺翼推進方面的可靠性，本文與Anderson 等[3]、Read 等[21]和Schou?veiler 等[4]實驗中的工況進行了同工況模擬驗證，具體參數(shù)可參照相應(yīng)的文獻。效率和推力系數(shù)的結(jié)果對比如圖12所示，從對比結(jié)果來看，本文計算結(jié)果與三個文獻中的實驗結(jié)果均比較吻合。圖12（b）中的cT為Anderson等[3]定義的推力系數(shù)，詳見文獻[13]。同時針對Schouveiler等[4]觀測的NACA 0012翼型的擺翼進行了同工況數(shù)值模擬，擺翼為展弦比6、兩端相對壁面滑動的準(zhǔn)二維平直翼，上下擺幅Y0=0.075 m，擺動頻率f=1.2 Hz，St=0.45，最大攻角為20°，觀測位置為中間截面，其對比如圖13 所示。圖13（a）和圖13（b）中分別顯示了翼型中間截面的尾渦實驗觀測圖和模擬云圖，可以看出模擬結(jié)果和實驗結(jié)果的尾渦形狀非常接近。綜上，可以證明本文對擺翼性能的計算方法較為可靠。

圖12 文獻與本文仿真數(shù)據(jù)的對比驗證（紅色曲線為本文工作，●為Anderson實驗結(jié)果，▲為Read實驗結(jié)果，▽為Schouveiler實驗結(jié)果。）Fig.12 Comparison of simulation results with previous computational and experimental results

圖13 St=0.45，amax=20°工況下文獻試驗結(jié)果與本文數(shù)值結(jié)果的尾渦對比圖Fig.13 Comparison of vorticity patterns visualized in the foil wake for St=0.45,amax=20°between the experimental results of Schouveiler et al[4]and the numerical results in this paper

3 試驗與模擬結(jié)果對比分析

本文裝置在哈爾濱工業(yè)大學(xué)（威海）的循環(huán)水池完成了相關(guān)實驗，試驗和數(shù)值模擬均在擺動頻率f=1.011 Hz，相位差φ=90o，升沉幅度Y0=80 mm 的工況下進行。在曲柄相位差設(shè)定分別為Δβ=0.6 rad 和Δβ=0.4 rad，擺角幅度分別約為0.23 rad和0.35 rad時，在VA=0.2 m/s、0.4 m/s、0.6 m/s和0.8 m/s四個來流速度下進行測量，得到擺動平板實際Fx、Fy和Mz。

3.1 性能和受力分析

圖14 給出了曲柄相位差分別為Δβ=0.6 rad 和Δβ=0.4 rad 兩個首搖運動曲線下效率和推力系數(shù)的實驗和數(shù)值模擬結(jié)果對比。從整體來看兩者趨勢一致，但實驗測得的效率和推力系數(shù)均略低于數(shù)值計算結(jié)果，推力系數(shù)的最左側(cè)點除外，曲柄相位差Δβ=0.6 rad 情況下的效率要比Δβ=0.4 rad 情況下更高一些。從圖14（a）可知，曲柄相位差的增加，即擺角的增加會在低進速系數(shù)下獲得較高的效率，但同時隨著進速系數(shù)的增加，效率下降較快。而小擺角擺翼卻能在較寬的進速系數(shù)下保持較高的效率，即小擺角擺翼的工作范圍較寬，這與課題組之前的研究結(jié)果相一致[22]。從圖14（b）可以發(fā)現(xiàn)，整體來看兩種不同曲柄相位差的推力系數(shù)KT差別不大，即首搖擺角幅度對推力系數(shù)影響不大，這是由于推力Fx在這兩組工況中變化很小，這可以從后面圖15的推力曲線時程圖對比得出；而在低進速系數(shù)，曲柄相位差Δβ=0.6 rad 的工況下推力系數(shù)KT較曲柄相位差Δβ=0.4 rad 的工況略大，但隨著進速系數(shù)的增加，曲柄相位差較大的工況推力系數(shù)下降較快，這與圖14（a）中大擺角擺翼效率下降較快相一致。

圖14 推進效率和推力系數(shù)的實驗結(jié)果與數(shù)值模擬結(jié)果的對比Fig.14 Comparison between experimental results and numerical simulation results of propulsion efficiency and thrust coefficient

為了詳細(xì)分析實驗結(jié)果與數(shù)值計算結(jié)果之間差異的原因，本節(jié)進一步對比了各工況下的實驗測試和數(shù)值模擬結(jié)果的推力時程曲線、側(cè)向力時程曲線和轉(zhuǎn)矩時程曲線，如圖15所示。首先，較普遍的現(xiàn)象是，與Fx和Mz相比，F(xiàn)y的實驗測量結(jié)果與數(shù)值計算結(jié)果更吻合一些。由公式（11）和公式（12）可知，F(xiàn)x的結(jié)果在很大程度上依賴Fs，而Fy的結(jié)果更依賴Fn，而從圖7中Fs和Fn的測量結(jié)果可以看出，弦向力Fs脈動十分劇烈且雜波較多，可靠性較低，這可能是測試系統(tǒng)受到循環(huán)水槽振動干擾引起的；而法向力Fn測試結(jié)果較為穩(wěn)定，且干擾較少，所以由其確定的Fy的實驗測量結(jié)果與數(shù)值計算結(jié)果更吻合一些。其次從橫向和縱向?qū)Ρ鹊慕Y(jié)果來看，當(dāng)曲柄相位差Δβ固定，即擺角固定時，隨著進速VA從0.2 m/s增加到0.8 m/s，試驗和數(shù)值模擬的推力Fx和側(cè)向力Fy的峰值變化很小，而轉(zhuǎn)矩Mz的峰值隨著進速的增大略有增加，約由0.3 N?m增加到0.4 N?m；而當(dāng)進速固定，曲柄相位差Δβ由0.4 rad增加到0.6 rad時可以發(fā)現(xiàn)，推力Fx峰值幾乎不變，而側(cè)向力Fy和轉(zhuǎn)矩Mz的峰值均有所下降。例如進速VA=0.2 m/s，曲柄相位差分別為0.4 rad 和0.6 rad 的圖15（a）和15（e）相比，側(cè)向力Fy從20 N 下降到15 N，而轉(zhuǎn)矩Mz從0.28 N?m下降到0.2 N?m，而推力Fx幾乎不變，約等于5 N。從效率的公式（8）可知，當(dāng)與輸出功率相關(guān)的推力Fx基本不變而與輸入功率相關(guān)的側(cè)向力Fy和轉(zhuǎn)矩Mz均有所下降時，效率必然增加，這正與前面圖14（a）中兩個工況效率曲線的對比結(jié)果相一致。最后，圖中還發(fā)現(xiàn)無論在試驗結(jié)果還是CFD模擬結(jié)果中，推力Fx峰值均呈現(xiàn)大小交替的周期變化，而且在試驗測試結(jié)果中更明顯，這是由于擺翼實驗裝置的首搖和橫蕩兩個運動的曲線都未達到絕對的正弦曲線形式所造成的。而側(cè)向力Fy和轉(zhuǎn)矩Mz的波形偏寬，也明顯呈不對稱狀態(tài)，是多個波形的組合，也是上述原因所造成。

圖15 不同曲柄相位差和流速下，F(xiàn)x、Fy和Mz隨時間的變化曲線Fig.15 Curves of Fx,Fy and Mz with time under different crank phase differences and flow velocities

3.2 效率分析

從圖14 的實驗結(jié)果可知，實驗測得的最高效率僅有27%，遠(yuǎn)低于同類型文獻中的數(shù)據(jù)。為了探索效率較低的問題，在擺角軌跡Δβ=0.4 rad的工況下，本節(jié)對平板前后緣增加導(dǎo)圓、無限展長假設(shè)和大擺幅情況三方面進行簡單探討，設(shè)置了4種工況進行數(shù)值計算對比：

工況1：矩形截面平板，升沉幅度Y0=80 mm；

工況2：矩形截面平板前后增加倒圓，倒圓半徑為1/2 板厚，即2.5 mm，其他計算參數(shù)同工況1，網(wǎng)格細(xì)節(jié)見圖16；

圖16 工況1和工況2下的網(wǎng)格細(xì)節(jié)Fig.16 Grid generation under case 1 and case 2

工況3：采用無限展長假設(shè)，分析平板底端和自由液面對水動力性能的影響。導(dǎo)圓平板為研究對象，其他技術(shù)參數(shù)同工況1；

工況4：將平板的橫蕩和首搖運動都改為簡諧規(guī)律，其中擺角幅度不變，橫蕩幅度從80 mm增加至200 mm。采用工況3的計算模型。

4 種工況的數(shù)值模擬效率曲線對比如圖17 所示。首先可以發(fā)現(xiàn)，工況2 的最高效率點在50%左右，與工況1 相比提高約20%。這說明矩形截面平板前后緣的導(dǎo)圓處理能夠非常明顯地提高推進效率。從表2計算結(jié)果可知，相比于工況1，推力Fx在工況2中明顯提高，由1.558 N 增大到2.491 N，增大約60%；而與此同時，工況2 中的功率P略有下降，由5.692 W 降低到4.809 W。根據(jù)公式（7）可知，工況2 推進效率η必然會增大。從圖18 的對比云圖中也可以看出，平板擺翼前后緣的倒圓處理可以明顯減小分離，進而可以降低阻力。其次，由圖17的效率對比可以看出，無限長假設(shè)下的工況3與工況2相比，效率提高幅度約10%左右。從圖19 工況2 的平板擺翼自由液面處的波高云圖可以看出，自由液面處有明顯興波，造成了阻力的增加，而考慮擺動平板的下端面，有限展長的三維擺翼前進的過程中會在翼尖處形成一個拖拽的旋渦，升力下降的同時產(chǎn)生誘導(dǎo)阻力，進而影響效率[21-22]。這也與本課題組前期的研究結(jié)果相一致，即較大展弦比的擺翼推進效率更接近無限長假設(shè)的理想二維擺翼，而隨著展弦比的降低，效率會降低。最后，相比于工況3，工況4 也有約10%的效率提高，其最高效率達到了70.7%。這個效率的提高一方面是由于擺幅比的增大對效率的貢獻，另一方面也有對工況1 中擺翼非簡諧運動規(guī)律的修正對效率的影響，這與Read 等[21]和Hover 等[23]的相關(guān)研究結(jié)論是一致的。

圖17 4種工況下的效率對比Fig.17 Efficiency comparison under four working conditions

表2 工況1和工況2數(shù)值計算結(jié)果對比Tab.2 Numerical results of Case 1 and Case 2

圖18 平板擺翼倒圓前后的尾渦變化Fig.18 Variation of wake vortices before and after rounding of the flapping plate

圖19 自由液面處的尾流場Fig.19 Wake field at free surface

3.3 功率分析

擺翼的推進效率為輸出功率與輸入功率之比。根據(jù)公式（7）和公式（8）可知，影響效率值的計算參數(shù)不僅包括機械參數(shù)Fx、Fy和Mz，還包括速度參數(shù)vy和ωz。當(dāng)擺翼工作時，這些參數(shù)隨時間而波動，參數(shù)之間的相位差對效率結(jié)果影響很大。

本文選取了3個代表性的工況點，在圖20中對比了其實驗與數(shù)值結(jié)果的Fy-vy和Mz-ωz相圖。3個工況分別為：實驗測得的最高效率點（Δβ=0.6 rad，VA=0.6 m/s）；最高效率點擺動規(guī)律下的較小進速工況（Δβ=0.6 rad，VA=0.2 m/s）；另一個擺動規(guī)律下的較大進速工況（Δβ=0.4 rad，VA=0.8 m/s）。從整體上看，無論是實驗結(jié)果還是數(shù)值計算結(jié)果，F(xiàn)y-vy和Mz-ωz都存在一定的相位差，F(xiàn)y和vy的相位差很小，而Mz和ωz的相位差較大，接近90°。同時，Mz與ωz峰值相乘的結(jié)果比Fy與vy峰值相乘的結(jié)果小得多，有數(shù)量級上的差別。綜合以上兩點說明y向力Fy在輸入功率中產(chǎn)生的功率更大，而繞Z軸的轉(zhuǎn)矩Mz幾乎沒有做功。

圖20 3個工況下的Fy-vy圖和Mz-ωz圖Fig.20 Phase diagrams of Fy-vy and Mz-ωz at 3 working points

4 結(jié) 論

本文設(shè)計了一種簡易的兩自由度擺翼試驗平臺，并在循環(huán)水槽開展了擺翼推進性能實驗研究，隨后對擺翼推進性能進行了數(shù)值仿真分析，并與實驗結(jié)果進行了對比分析。最后，針對實驗中擺翼效率較低的情況，數(shù)值分析了平板前后緣導(dǎo)圓、無限展長假設(shè)和大擺幅對效率的影響，主要結(jié)論如下：

（1）本文設(shè)計了一套利用曲柄連桿滑塊機構(gòu)實現(xiàn)擺翼推進的實驗裝置，并測量了擺動平板在不同來流速度和首搖幅度下的水動力性能。同時，對實驗工況進行了數(shù)值模擬，以便于對比研究；

（2）同時對實驗測量與數(shù)值分析的總體性能和各個分力進行了對比。從對比結(jié)果來看，實驗結(jié)果與數(shù)值結(jié)果之間略有偏差，實驗測得的效率和推力系數(shù)普遍要略低于數(shù)值計算結(jié)果。從力的對比來看，法向力Fn吻合較好，弦向力Fs和轉(zhuǎn)矩Mz稍差；

（3）為了分析擺動平板效率較低的原因，設(shè)計了另外3 個工況，通過數(shù)值模擬得到了與實驗工況下效率計算結(jié)果的對比。從對比結(jié)果來看，平板前后緣導(dǎo)圓、消除端壁和自由液面的影響和增大擺幅都能明顯地提高擺動平板的效率；

（4）對擺動平板的輸入功率分量進行了分析，從結(jié)果來看，在擺動平板的工作過程中，Mz幾乎不做功，所有輸入功率近乎由Fy提供。

附錄A：橫蕩運動與首搖運動的表達式推導(dǎo)過程

根據(jù)正文中圖2，通過給定A點和B點的位移曲線，計算平板中點的橫蕩運動位移曲線，及其首搖運動的擺角曲線。B點相對A點延后Δβ相位，并且相對位置右偏LBC長度。根據(jù)圖中坐標(biāo)，在曲柄連桿工作平面內(nèi)（YZ平面），可得兩點的運動曲線分別為

式中，y0為曲柄長度，L為連桿長度，其他變量見論文正文。

兩滑軌間距為S1=LAC，則在XY平面內(nèi)，可得C點的運動曲線。首先計算距離LAB

則擺翼中點坐標(biāo)為

文中取y在其平均值上下波動。此外，擺翼中點的X坐標(biāo)也會有微小的波動，本文中對此忽略不計。