基于CLE方法的有機材料對柔性擺動雙翼性能影響研究

張 營 川， 馬 駿， 王 濤＊， 鄭 來 新， 張 春 曉， 吳 濤

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

0 引 言

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

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

1 流固耦合方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 模型實驗驗證

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

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

圖1 單翼運動示意圖以及有限元模型Fig.1 Sketch of single foil motion and finite element model

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

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

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

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

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

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

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

圖3 雙翼運動示意圖以及有限元模型Fig.3 Sketch of the dual－foils motion and finite element model

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

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

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

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

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

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

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

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

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

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

從圖中可以看出，第一，在較低的St下，隨著St的增加，無論是金屬雙翼還是橡膠雙翼的水動力性能曲線都呈上升的趨勢，也就是說雙翼的推力系數(shù)和推進效率都逐漸增大；第二，在相同的St下，橡膠雙翼的水動力性能明顯優(yōu)于金屬雙翼.從兩個方面來進行分析：在雙翼變形方面，圖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)有機材料雙翼的變形量比金屬雙翼更大；在雙翼尾流方面，圖4中金屬雙翼只有2個有標(biāo)記的旋渦，而有機材料雙翼有3個有標(biāo)記的旋渦，這說明有機雙翼的變形對雙翼后方的旋渦發(fā)放有較大的影響.因此，橡膠雙翼的水動力性能明顯優(yōu)于金屬雙翼的原因在于，橡膠雙翼在流固耦合過程中發(fā)生了明顯變形，從而改變了翼型的尾渦運動軌跡，也就改變了旋渦發(fā)放產(chǎn)生的渦激力的大小和方向，最終改進了雙翼的水動力性能.

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

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

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

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

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

4 結(jié) 論

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

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

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

[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] 周 岱，馬 駿，李 磊.流場和流固耦合問題網(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] 張曉慶，王志東，張振山.二維擺動水翼仿生推進水動力性能研究 [J].水動力學(xué)研究與進展，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] 李裕春，時黨勇，趙 遠.ANSYS 10.0/LS－DYNA基礎(chǔ)理論與工程實踐[M].北京：中國水利水電出版社，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.