張 營(yíng) 川， 馬 駿， 王 濤＊， 鄭 來(lái) 新， 張 春 曉， 吳 濤

（1.大連理工大學(xué) 船舶工程學(xué)院，遼寧 大連 116024；2.總裝工程兵駐武漢地區(qū)軍事代表室，湖北 武漢 430073）

0 引 言

當(dāng)前，大部分船舶都采用螺旋槳作為推進(jìn)器，但是作為傳統(tǒng)的剛性推進(jìn)方式，它不可避免地存在一些弊端，如機(jī)械傳動(dòng)復(fù)雜、推進(jìn)效率受到限制、結(jié)構(gòu)尺寸和重量大、環(huán)境噪音大等.渦動(dòng)力推進(jìn)技術(shù)在這樣的背景下孕育而生，逐漸成為水下推進(jìn)技術(shù)研究領(lǐng)域的熱點(diǎn)之一，擺動(dòng)單翼和擺動(dòng)雙翼就是其中的代表.它具有高效率、高機(jī)動(dòng)性、完善的流體性能、低噪聲、推進(jìn)器和舵合二為一、驅(qū)動(dòng)方式多樣等優(yōu)點(diǎn)，是深水機(jī)器人、魚雷、潛艇等水下運(yùn)載工具的理想推進(jìn)器.然而，完全從理論上求解擺翼性能還存在相當(dāng)?shù)碾y度.因此，基于模型實(shí)驗(yàn)和數(shù)值模擬的方法，相關(guān)研究逐漸深入[1－3].但是以往的研究工作存在以下不足：學(xué)者大多都是基于單一的拉格朗日動(dòng)網(wǎng)格技術(shù)及相應(yīng)的網(wǎng)格修正方法開(kāi)展流固耦合研究，如局部網(wǎng)格重劃法等[4]；模型處理過(guò)程中，需要花費(fèi)大量精力區(qū)分流體與固體之間邊界[5]；在動(dòng)網(wǎng)格技術(shù)中，流體網(wǎng)格扭曲、畸變影響計(jì)算精度等[6].

針對(duì)上述情況，本文在渦動(dòng)力推進(jìn)器研究領(lǐng)域采用相容拉格朗日－歐拉方法（compatible Lagrangian－Eulerian method，以下簡(jiǎn)稱 CLE 方法）對(duì)柔性擺動(dòng)雙翼的流固耦合過(guò)程進(jìn)行建模和求解；研究有機(jī)材料對(duì)柔性擺動(dòng)雙翼水動(dòng)力性能和結(jié)構(gòu)性能的影響.

1 流固耦合方法

本文采用的流固耦合方法為CLE方法[7]，這種方法來(lái)源于任意拉格朗日－歐拉方法（arbitrary Lagrangian－Eulerian method，以下簡(jiǎn)稱 ALE 方法）.ALE方法中，流體用在空間任意變形和運(yùn)動(dòng)的網(wǎng)格描述，結(jié)構(gòu)用拉格朗日網(wǎng)格描述；CLE方法中，流體被劃分為空間位置不變的歐拉網(wǎng)格，結(jié)構(gòu)仍然采用拉格朗日網(wǎng)格描述.

CLE方法有兩個(gè)主要優(yōu)點(diǎn)：一是建模過(guò)程中網(wǎng)格可以相互重疊；二是數(shù)值求解方面，流體網(wǎng)格不發(fā)生運(yùn)動(dòng)，通過(guò)流體網(wǎng)格的物質(zhì)傳遞表達(dá)流體的流動(dòng)，避免了動(dòng)網(wǎng)格技術(shù)中流體網(wǎng)格扭曲的問(wèn)題.

從數(shù)學(xué)表達(dá)式的角度講，CLE方法與ALE方法唯一的不同點(diǎn)在于，流體網(wǎng)格速度u＝0，即表示流體域不發(fā)生運(yùn)動(dòng)，將之代入ALE方法的各個(gè)公式中，最后就會(huì)得到CLE方法的基本公式.

在ALE方法的描述中，引入Lagrange和Euler坐標(biāo)之外的第3個(gè)任意參照坐標(biāo).與參照坐標(biāo)相關(guān)的材料微商可以采用下式描述：

式中：f為ALE方法中與坐標(biāo)和時(shí)間相關(guān)的任意變量；Xi為拉格朗日坐標(biāo)；xi為歐拉坐標(biāo)；wi為相對(duì)速度.相對(duì)速度通過(guò)下式表達(dá)：

式中：v為物質(zhì)速度；u為網(wǎng)格速度.

從Navier－Stokes方程可以推導(dǎo)出ALE方法的控制方程，并以質(zhì)量守恒、動(dòng)量守恒和能量守恒方程的形式給出：

式中：ρ為流體密度；bi為單位質(zhì)量的體積力分量；σij為Cauchy應(yīng)力張量分量；E為內(nèi)能.Cauchy應(yīng)力張量分量σij能夠以Stokes本構(gòu)方程的形式表示，具體表達(dá)式如下：

式中：p為壓力；δij為Kroneckerδ函數(shù)；μ為流體的動(dòng)力黏性系數(shù).

在流固耦合計(jì)算中，計(jì)算的第一個(gè)階段執(zhí)行拉格朗日過(guò)程，此時(shí)網(wǎng)格隨物質(zhì)運(yùn)動(dòng).只考慮拉格朗日網(wǎng)格，由于沒(méi)有物質(zhì)流經(jīng)單元邊界，質(zhì)量自動(dòng)保持守恒.計(jì)算的第二個(gè)階段，即對(duì)流相，對(duì)穿過(guò)單元邊界的質(zhì)量輸運(yùn)、內(nèi)能和動(dòng)量進(jìn)行計(jì)算，這是將拉格朗日過(guò)程的移位網(wǎng)格重映射回其初始位置或任意位置.

根據(jù)CLE方法和LS－DYNA求出必需的流場(chǎng)參數(shù)和結(jié)構(gòu)參數(shù)后，再根據(jù)下列公式就可以得到擺翼的水動(dòng)力性能，而擺翼的結(jié)構(gòu)動(dòng)力響應(yīng)可以直接在軟件中顯示.擺翼的沉浮和擺動(dòng)運(yùn)動(dòng)方程表示如下（兩種運(yùn)動(dòng)之間的相位角為π/2）：

式中：y0為沉浮運(yùn)動(dòng)幅值；θ0為擺動(dòng)運(yùn)動(dòng)幅值；ω為沉浮和擺動(dòng)運(yùn)動(dòng)的角頻率.在擺翼研究中，Strouhal數(shù)通常定義如下[8]：

式中：U為來(lái)流速度.推力系數(shù)按照下式求出：

式中：n為周期個(gè)數(shù)；T為周期；L為升力；M為扭矩.

2 模型實(shí)驗(yàn)驗(yàn)證

本文從一系列經(jīng)典的擺動(dòng)單翼推進(jìn)性能實(shí)驗(yàn)結(jié)果中選取部分實(shí)驗(yàn)結(jié)果[8]，作為應(yīng)用CLE方法和LS－DYNA對(duì)渦動(dòng)力推進(jìn)器進(jìn)行數(shù)值研究是否可靠的驗(yàn)證，數(shù)值仿真的輸入數(shù)據(jù)基于選取的實(shí)驗(yàn)條件數(shù)據(jù).

實(shí)驗(yàn)是在長(zhǎng)30m、寬2.6m、深1.2m 的MIT Tow水池中進(jìn)行的.實(shí)驗(yàn)選取的剛性翼型為NACA 0012，翼展0.6m，弦長(zhǎng)0.1m.翼型在水平方向以速度U前進(jìn)，同時(shí)以角頻率ω?cái)[動(dòng)并且上下沉浮，如圖1（a）所示[8].翼型在拖曳機(jī)械的帶動(dòng)下水平前進(jìn)速度為0.4m/s，St為0.20、0.25、0.30、0.35、0.40，對(duì)應(yīng)的ω為3.35、4.19、5.03、5.86、6.07rad/s.

圖1 單翼運(yùn)動(dòng)示意圖以及有限元模型Fig.1 Sketch of single foil motion and finite element model

有限元模型如圖1（b）所示，由于計(jì)算能力所限，計(jì)算水域不能完全按照水池大小來(lái)建立，其在長(zhǎng)度方向上為6個(gè)翼型弦長(zhǎng)，寬度方向上為5個(gè)翼型弦長(zhǎng)，厚度方向上為1個(gè)有限單元的厚度；計(jì)算模型包含30 000個(gè)六面體單元，翼型運(yùn)動(dòng)區(qū)域周圍的流體網(wǎng)格全部采用了網(wǎng)格加密；計(jì)算模型右側(cè)邊界面施加了流入邊界條件，左側(cè)邊界面施加了流出邊界條件，其他邊界面施加的是無(wú)反射邊界條件，整個(gè)水域施加了初始流入速度，可以保證流場(chǎng)的連貫；計(jì)算時(shí)間為4s，最多包含4個(gè)翼型運(yùn)動(dòng)周期，時(shí)間步長(zhǎng)Δt受到Courant穩(wěn)定性條件的限制不能取得太大，Δt為1.67×10－6s，在計(jì)算過(guò)程中還采用了LS－DYNA中的質(zhì)量縮放方法來(lái)節(jié)省計(jì)算時(shí)間；流固耦合計(jì)算公式及具體迭代方法見(jiàn)第1章.

本文將數(shù)值模擬結(jié)果、實(shí)驗(yàn)數(shù)據(jù)處理結(jié)果、無(wú)黏理論計(jì)算結(jié)果[8]進(jìn)行了對(duì)比，在兩種最大攻角情況下得到了平均推力系數(shù)，如圖2所示，圖2（a）最大攻角為15°，圖2（b）最大攻角為30°.

從圖2中可以看出，在較低的St下，即St從0.2到0.4，推力系數(shù)隨St的增加而增大.對(duì)圖中的曲線進(jìn)行分析，無(wú)黏理論計(jì)算得到的結(jié)果明顯高估了推力系數(shù)，這主要是因?yàn)樵谠摕o(wú)黏理論中，忽略了翼型表面受到的黏性阻力，從而放大了翼型的推力，而且其計(jì)算流體在翼型尾部的分離時(shí)也存在一定的局限性，最終導(dǎo)致推力系數(shù)比實(shí)驗(yàn)值要大；同樣，本文數(shù)值模擬結(jié)果也存在一些偏差.原因在于建模過(guò)程中考慮到計(jì)算機(jī)的處理能力，沒(méi)有完全按照水池大小建立翼型與流體的三維模型，而是采用單層的立體網(wǎng)格建立模型，固體和流體網(wǎng)格也只能盡量細(xì)化，所以數(shù)值模擬結(jié)果并沒(méi)有與實(shí)驗(yàn)處理結(jié)果完全一致.盡管如此，總體上本文的數(shù)值模擬結(jié)果反映出了推力系數(shù)隨St的變化趨勢(shì)，而且求解ALE系列方程時(shí)考慮了流體的黏性，所以比無(wú)黏理論的計(jì)算結(jié)果更加接近實(shí)驗(yàn)數(shù)據(jù).

圖2 不同方法計(jì)算得到的推力系數(shù)曲線對(duì)比圖Fig.2 Comparison of thrust force coefficient curves obtained from different methods

3 材料類型對(duì)柔性擺動(dòng)雙翼性能影響的數(shù)值模擬

關(guān)于柔性擺動(dòng)雙翼的研究中，大部分文獻(xiàn)都是基于剛性翼型開(kāi)展的水動(dòng)力性能計(jì)算與分析，然而剛性翼型由于無(wú)法變形，也就無(wú)從研究雙翼結(jié)構(gòu)力學(xué)方面的性能，特別是不能提供雙翼內(nèi)部應(yīng)力分布情況，而這對(duì)于雙翼推進(jìn)器疲勞性能的研究是至關(guān)重要的.下面通過(guò)數(shù)值模擬方法，對(duì)上述問(wèn)題進(jìn)行計(jì)算分析.

通過(guò)反復(fù)在雙翼上施加上下的沉浮運(yùn)動(dòng)，并限制雙翼繞自身的轉(zhuǎn)動(dòng)，就可以依靠雙翼的旋渦發(fā)放所產(chǎn)生的渦激力推動(dòng)雙翼前進(jìn)，圖3給出了柔性擺動(dòng)雙翼推進(jìn)過(guò)程的示意圖（圖3（a））和流固耦合有限元模型（圖3（b））.

圖3 雙翼運(yùn)動(dòng)示意圖以及有限元模型Fig.3 Sketch of the dual－foils motion and finite element model

本文選用4種不同類型的材料分別組成4種柔性擺動(dòng)雙翼，分別是金屬材料、聚乙烯材料以及兩種橡膠材料，其中橡膠雙翼直接從LS－DYNA中選用非線性橡膠單元組成.可以看作剛性體的金屬材料，彈性模量為0.2×1012Pa，泊松比為0.3，密度為7 800kg/m3；有機(jī)材料聚乙烯，彈性模量為2.62×108Pa，泊松比為0.3，密度為920 kg/m3；有機(jī)材料硬度為60IRHD的橡膠，泊松比為0.499，密度為1 300kg/m3，彈性模量相關(guān)系數(shù)為0.70×106和0.35×106；有機(jī)材料硬度為40IRHD的橡膠，泊松比為0.499，密度為1 300 kg/m3，彈性模量相關(guān)系數(shù)為0.275×106和0.028×106.由于橡膠的種類繁多，資料中彈性模量通常不直接給出，而是給出相關(guān)系數(shù)，可以直接輸入LS－DYNA中進(jìn)行非線性單元計(jì)算.

選取的翼型為NACA 0012，翼展為0.6m，弦長(zhǎng)為0.1m.兩翼型運(yùn)動(dòng)時(shí)最近距離為0.015 m，沉浮運(yùn)動(dòng)最大振幅為3/4弦長(zhǎng)，即0.075m.在來(lái)流速度為0.4m/s，St為0.4，翼型沉浮運(yùn)動(dòng)角頻率為6.7rad/s，沉浮運(yùn)動(dòng)周期為0.94s時(shí)，對(duì)4種不同類型的材料組成的柔性擺動(dòng)雙翼推進(jìn)過(guò)程進(jìn)行數(shù)值仿真.

由于翼型結(jié)構(gòu)內(nèi)部的應(yīng)力分布與結(jié)構(gòu)約束以及動(dòng)力施加位置有關(guān)，需要進(jìn)行相關(guān)說(shuō)明.本文中，擺翼的約束位置也是外部動(dòng)力的施加位置，它在靠近翼型頭部，翼型弦長(zhǎng)的1/3處.改變約束位置和外部動(dòng)力的施加位置，將會(huì)引起翼型內(nèi)部應(yīng)力分布的變化，詳細(xì)的變化規(guī)律將在以后的論文中進(jìn)行研究，本文基于上述固定的約束和外力位置進(jìn)行討論.

3.1 材料類型對(duì)雙翼水動(dòng)力性能影響的數(shù)值模擬計(jì)算與分析

由于雙翼是關(guān)于對(duì)稱軸對(duì)稱的，為了更加清晰，只選取了雙翼中的一個(gè)翼型給出數(shù)值模擬圖像.受到計(jì)算軟件當(dāng)前版本的功能限制，對(duì)于含固體變形的流固耦合計(jì)算，暫時(shí)還不能直接輸出渦量分布圖，但可以提供流體的速度矢量云圖作為替代，以便討論材料類型對(duì)翼型后方發(fā)放的旋渦形狀的影響，如圖4所示.表1為4種柔性擺動(dòng)雙翼的水動(dòng)力性能計(jì)算結(jié)果.

圖4 不同材料柔性擺動(dòng)雙翼后方流體的速度矢量云圖Fig.4 Contours of fluid resultant velocity vector of flexible flapping dual－foils with different materials

表1 不同材料柔性擺動(dòng)雙翼水動(dòng)力性能Tab.1 Hydrodynamic performance of flexible flapping dual－foils with different materials

圖4和表1表明，由于有機(jī)材料的彈性模量比金屬材料的彈性模量要小，有機(jī)材料組成的柔性擺動(dòng)雙翼在推進(jìn)過(guò)程中與流體相互作用時(shí)比較容易發(fā)生變形，這種變形反過(guò)來(lái)又影響到旋渦運(yùn)動(dòng)的軌跡，因此從流體的速度云圖上可以看出，金屬雙翼發(fā)放的旋渦軌跡明顯有別于有機(jī)材料組成的柔性擺動(dòng)雙翼發(fā)放的旋渦軌跡；另外，聚乙烯雙翼、60IRHD橡膠雙翼和40IRHD橡膠雙翼之間材料彈性模量的差別并不大，所以它們的尾部旋渦軌跡非常相似，導(dǎo)致它們水動(dòng)力性能計(jì)算結(jié)果的差別也不大，但是都優(yōu)于金屬雙翼的計(jì)算結(jié)果.

上述計(jì)算結(jié)果是在St為0.4時(shí)得到的，還不夠全面.為了進(jìn)一步研究材料對(duì)雙翼水動(dòng)力性能的影響，在St從0.2到0.4時(shí)，分別對(duì)由金屬材料組成的雙翼和由40IRHD橡膠材料組成的雙翼進(jìn)行數(shù)值計(jì)算，它們的推力系數(shù)曲線和推進(jìn)效率曲線如圖5所示.

圖5 兩種不同材料翼型的水動(dòng)力性能曲線Fig.5 Hydrodynamic performance curves of foils with two kinds of materials

從圖中可以看出，第一，在較低的St下，隨著St的增加，無(wú)論是金屬雙翼還是橡膠雙翼的水動(dòng)力性能曲線都呈上升的趨勢(shì)，也就是說(shuō)雙翼的推力系數(shù)和推進(jìn)效率都逐漸增大；第二，在相同的St下，橡膠雙翼的水動(dòng)力性能明顯優(yōu)于金屬雙翼.從兩個(gè)方面來(lái)進(jìn)行分析：在雙翼變形方面，圖4中金屬雙翼的最大有效應(yīng)變?yōu)?.00×10－7、聚乙烯雙翼的最大有效應(yīng)變?yōu)?.02×10－6、60IRHD橡膠雙翼的最大有效應(yīng)變?yōu)?.84×10－5、40IRHD橡膠雙翼的最大有效應(yīng)變?yōu)?.39×10－5，可以發(fā)現(xiàn)有機(jī)材料雙翼的變形量比金屬雙翼更大；在雙翼尾流方面，圖4中金屬雙翼只有2個(gè)有標(biāo)記的旋渦，而有機(jī)材料雙翼有3個(gè)有標(biāo)記的旋渦，這說(shuō)明有機(jī)雙翼的變形對(duì)雙翼后方的旋渦發(fā)放有較大的影響.因此，橡膠雙翼的水動(dòng)力性能明顯優(yōu)于金屬雙翼的原因在于，橡膠雙翼在流固耦合過(guò)程中發(fā)生了明顯變形，從而改變了翼型的尾渦運(yùn)動(dòng)軌跡，也就改變了旋渦發(fā)放產(chǎn)生的渦激力的大小和方向，最終改進(jìn)了雙翼的水動(dòng)力性能.

3.2 材料類型對(duì)雙翼結(jié)構(gòu)力學(xué)性能影響的數(shù)值模擬計(jì)算與分析

上述計(jì)算在輸出雙翼水動(dòng)力性能數(shù)據(jù)的同時(shí)，也輸出了對(duì)應(yīng)的結(jié)構(gòu)力學(xué)性能數(shù)據(jù).圖6為不同雙翼的內(nèi)部應(yīng)力云圖，主要展示的是雙翼變形程度和應(yīng)力集中的位置；表2為雙翼推進(jìn)過(guò)程中，提取出的內(nèi)部最大有效應(yīng)力數(shù)據(jù)，主要展示的是材料類型對(duì)結(jié)構(gòu)內(nèi)部應(yīng)力的改變程度.

圖6 不同材料雙翼有效應(yīng)力云圖Fig.6 Effective stress contours of dual－foils with different materials

表2 不同材料雙翼最大有效應(yīng)力Tab.2 Maximum effective stress of dual－foils with different materials

因?yàn)椴伙@示周圍的流體，并且將翼型進(jìn)行了放大，所以圖6比圖4更加清晰地反映了雙翼的變形情況，由于材料的彈性模量各不相同，金屬雙翼在整個(gè)推進(jìn)過(guò)程中幾乎不發(fā)生變形，而有機(jī)雙翼都有不同程度的形狀改變；另外，應(yīng)力云圖主要體現(xiàn)的是雙翼在工作時(shí)內(nèi)部的應(yīng)力分布情況，可以看出4種雙翼都出現(xiàn)了應(yīng)力集中的現(xiàn)象，而且應(yīng)力集中的區(qū)域基本上都在雙翼的中部.這是因?yàn)殡p翼的沉浮運(yùn)動(dòng)荷載加載在靠近雙翼頭部的1/3處，而雙翼的旋轉(zhuǎn)約束也在該處，所以雙翼運(yùn)動(dòng)時(shí)旋渦發(fā)放產(chǎn)生的渦激力最終在靠近約束的位置產(chǎn)生了應(yīng)力集中.表2定量計(jì)算了4種雙翼應(yīng)力集中的程度，其中，金屬雙翼的內(nèi)部最大有效應(yīng)力遠(yuǎn)遠(yuǎn)大于有機(jī)材料組成的柔性擺動(dòng)雙翼，而在有機(jī)雙翼中，40IRHD橡膠的內(nèi)部最大有效應(yīng)力最低.因此，在外部條件相同的情況下，材料的柔性是雙翼內(nèi)部最大有效應(yīng)力的主要影響因素.

進(jìn)一步的計(jì)算給出了40IRHD橡膠材料組成的雙翼，在較低的St下，即St從0.2到0.4時(shí)，計(jì)算得到的最大有效應(yīng)力（σe，max）曲線和最大剪切應(yīng)力（τmax）曲線（如圖7所示）.隨著St的增加，橡膠雙翼的最大有效應(yīng)力和最大剪切應(yīng)力都會(huì)變大，它們代表的是整個(gè)翼型各個(gè)截面當(dāng)中應(yīng)力的最大值.

圖7 40IRHD橡膠雙翼最大應(yīng)力曲線Fig.7 Maximum stress curves of rubber dual－foils with 40IRHD

4 結(jié) 論

（1）有機(jī)材料能改善雙翼水動(dòng)力性能.與金屬雙翼相比，有機(jī)雙翼的彈性模量較低，翼型結(jié)構(gòu)更柔軟，更加容易受到旋渦運(yùn)動(dòng)的影響產(chǎn)生較大的變形，這種變形反過(guò)來(lái)又改變了雙翼尾流區(qū)域旋渦的形狀，調(diào)整了渦激力的大小和方向，從而提高了雙翼渦動(dòng)力推進(jìn)器的推進(jìn)性能.

（2）有機(jī)材料能提高雙翼抗疲勞破壞能力.在雙翼沉浮運(yùn)動(dòng)過(guò)程中，雙翼在周期性流體力的反復(fù)作用下會(huì)出現(xiàn)應(yīng)力集中現(xiàn)象，翼型結(jié)構(gòu)交變應(yīng)力集中區(qū)域位于雙翼的中部.有機(jī)雙翼結(jié)構(gòu)比較柔軟而且容易變形，正是這種變形緩解了應(yīng)力集中，最后降低了結(jié)構(gòu)受到的最大交變應(yīng)力值，這能夠優(yōu)化柔性擺動(dòng)雙翼的抗疲勞性能.

綜上所述，有機(jī)材料對(duì)雙翼性能的影響是有利的，選擇合適的有機(jī)材料制造翼型結(jié)構(gòu)，不僅可以提高推進(jìn)性能，也能夠延長(zhǎng)其使用壽命.

[1] Hover F S，Haugsdal ，Triantafyllou M S.Effect of angle of attack profiles in flapping foil propulsion [J].Journal of Fluids and Structures，2004，19（1）：37－47.

[2] von Ellenrieder K D，Parker K，Soria J.Fluid mechanics of flapping wings [J].Experimental Thermal and Fluid Science，2008，32（8）：1578－1589.

[3] Jones K D，Plazter M F.An experimental and numerical investigation of flapping－wing propulsion[C]//37th Aerospace Sciences Meeting ＆ Exhibit.Reston：AIAA，1999.

[4] Miao Jr－ming，Sun Wei－h(huán)sin.Numerical analysis on aerodynamic force generation of biplane counterflapping flexible airfoils [J].Journal of Aircraft，2009，46（5）：1785－1794.

[5] 周 岱，馬 駿，李 磊.流場(chǎng)和流固耦合問(wèn)題網(wǎng)格剖分與更新的新方法[J].工程力學(xué)，2009，26（11）：10－15.ZHOU Dai，MA Jun，LI Lei.A novel technique for mesh generation and update in the analysis of flow field and fluid－structure interaction [J].Engineering Mechanics，2009，26（11）：10－15.（in Chinese）

[6] 張曉慶，王志東，張振山.二維擺動(dòng)水翼仿生推進(jìn)水動(dòng)力性能研究 [J].水動(dòng)力學(xué)研究與進(jìn)展，2006，21（5）：632－639.ZHANG Xiao－qing， WANG Zhi－dong， ZHANG Zhen－shan.Hydrodynamic study of bionic propulsion for 2－D flapping foil[J].Journal of Hydrodynamics，2006，21（5）：632－639.（in Chinese）

[7] 李裕春，時(shí)黨勇，趙 遠(yuǎn).ANSYS 10.0/LS－DYNA基礎(chǔ)理論與工程實(shí)踐[M].北京：中國(guó)水利水電出版社，2008.LI Yu－chun，SHI Dang－yong，ZHAO Yuan.Basic Theory and Engineering Practice of ANSYS 10.0/LSDYNA[M].Beijing：China Waterpower Press，2008.（in Chinese）

[8] Schouveiler L， Hover F S，Triantafyllou M S.Performance of flapping foil propulsion[J].Journal of Fluids and Structures，2005，20（7）：949－959.