基于CLE方法的有机材料对柔性摆动双翼性能影响研究

张 营 川， 马 骏， 王 涛＊， 郑 来 新， 张 春 晓， 吴 涛

（1.大连理工大学 船舶工程学院，辽宁 大连 116024；2.总装工程兵驻武汉地区军事代表室，湖北 武汉 430073）

0 引 言

当前，大部分船舶都采用螺旋桨作为推进器，但是作为传统的刚性推进方式，它不可避免地存在一些弊端，如机械传动复杂、推进效率受到限制、结构尺寸和重量大、环境噪音大等.涡动力推进技术在这样的背景下孕育而生，逐渐成为水下推进技术研究领域的热点之一，摆动单翼和摆动双翼就是其中的代表.它具有高效率、高机动性、完善的流体性能、低噪声、推进器和舵合二为一、驱动方式多样等优点，是深水机器人、鱼雷、潜艇等水下运载工具的理想推进器.然而，完全从理论上求解摆翼性能还存在相当的难度.因此，基于模型实验和数值模拟的方法，相关研究逐渐深入[1－3].但是以往的研究工作存在以下不足：学者大多都是基于单一的拉格朗日动网格技术及相应的网格修正方法开展流固耦合研究，如局部网格重划法等[4]；模型处理过程中，需要花费大量精力区分流体与固体之间边界[5]；在动网格技术中，流体网格扭曲、畸变影响计算精度等[6].

针对上述情况，本文在涡动力推进器研究领域采用相容拉格朗日－欧拉方法（compatible Lagrangian－Eulerian method，以下简称 CLE 方法）对柔性摆动双翼的流固耦合过程进行建模和求解；研究有机材料对柔性摆动双翼水动力性能和结构性能的影响.

1 流固耦合方法

本文采用的流固耦合方法为CLE方法[7]，这种方法来源于任意拉格朗日－欧拉方法（arbitrary Lagrangian－Eulerian method，以下简称 ALE 方法）.ALE方法中，流体用在空间任意变形和运动的网格描述，结构用拉格朗日网格描述；CLE方法中，流体被划分为空间位置不变的欧拉网格，结构仍然采用拉格朗日网格描述.

CLE方法有两个主要优点：一是建模过程中网格可以相互重叠；二是数值求解方面，流体网格不发生运动，通过流体网格的物质传递表达流体的流动，避免了动网格技术中流体网格扭曲的问题.

从数学表达式的角度讲，CLE方法与ALE方法唯一的不同点在于，流体网格速度u＝0，即表示流体域不发生运动，将之代入ALE方法的各个公式中，最后就会得到CLE方法的基本公式.

在ALE方法的描述中，引入Lagrange和Euler坐标之外的第3个任意参照坐标.与参照坐标相关的材料微商可以采用下式描述：

式中：f为ALE方法中与坐标和时间相关的任意变量；Xi为拉格朗日坐标；xi为欧拉坐标；wi为相对速度.相对速度通过下式表达：

式中：v为物质速度；u为网格速度.

从Navier－Stokes方程可以推导出ALE方法的控制方程，并以质量守恒、动量守恒和能量守恒方程的形式给出：

式中：ρ为流体密度；bi为单位质量的体积力分量；σij为Cauchy应力张量分量；E为内能.Cauchy应力张量分量σij能够以Stokes本构方程的形式表示，具体表达式如下：

式中：p为压力；δij为Kroneckerδ函数；μ为流体的动力黏性系数.

在流固耦合计算中，计算的第一个阶段执行拉格朗日过程，此时网格随物质运动.只考虑拉格朗日网格，由于没有物质流经单元边界，质量自动保持守恒.计算的第二个阶段，即对流相，对穿过单元边界的质量输运、内能和动量进行计算，这是将拉格朗日过程的移位网格重映射回其初始位置或任意位置.

根据CLE方法和LS－DYNA求出必需的流场参数和结构参数后，再根据下列公式就可以得到摆翼的水动力性能，而摆翼的结构动力响应可以直接在软件中显示.摆翼的沉浮和摆动运动方程表示如下（两种运动之间的相位角为π/2）：

式中：y0为沉浮运动幅值；θ0为摆动运动幅值；ω为沉浮和摆动运动的角频率.在摆翼研究中，Strouhal数通常定义如下[8]：

式中：U为来流速度.推力系数按照下式求出：

式中：n为周期个数；T为周期；L为升力；M为扭矩.

2 模型实验验证

本文从一系列经典的摆动单翼推进性能实验结果中选取部分实验结果[8]，作为应用CLE方法和LS－DYNA对涡动力推进器进行数值研究是否可靠的验证，数值仿真的输入数据基于选取的实验条件数据.

实验是在长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，对应的ω为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个六面体单元，翼型运动区域周围的流体网格全部采用了网格加密；计算模型右侧边界面施加了流入边界条件，左侧边界面施加了流出边界条件，其他边界面施加的是无反射边界条件，整个水域施加了初始流入速度，可以保证流场的连贯；计算时间为4s，最多包含4个翼型运动周期，时间步长Δt受到Courant稳定性条件的限制不能取得太大，Δt为1.67×10－6s，在计算过程中还采用了LS－DYNA中的质量缩放方法来节省计算时间；流固耦合计算公式及具体迭代方法见第1章.

本文将数值模拟结果、实验数据处理结果、无黏理论计算结果[8]进行了对比，在两种最大攻角情况下得到了平均推力系数，如图2所示，图2（a）最大攻角为15°，图2（b）最大攻角为30°.

从图2中可以看出，在较低的St下，即St从0.2到0.4，推力系数随St的增加而增大.对图中的曲线进行分析，无黏理论计算得到的结果明显高估了推力系数，这主要是因为在该无黏理论中，忽略了翼型表面受到的黏性阻力，从而放大了翼型的推力，而且其计算流体在翼型尾部的分离时也存在一定的局限性，最终导致推力系数比实验值要大；同样，本文数值模拟结果也存在一些偏差.原因在于建模过程中考虑到计算机的处理能力，没有完全按照水池大小建立翼型与流体的三维模型，而是采用单层的立体网格建立模型，固体和流体网格也只能尽量细化，所以数值模拟结果并没有与实验处理结果完全一致.尽管如此，总体上本文的数值模拟结果反映出了推力系数随St的变化趋势，而且求解ALE系列方程时考虑了流体的黏性，所以比无黏理论的计算结果更加接近实验数据.

图2 不同方法计算得到的推力系数曲线对比图Fig.2 Comparison of thrust force coefficient curves obtained from different methods

3 材料类型对柔性摆动双翼性能影响的数值模拟

关于柔性摆动双翼的研究中，大部分文献都是基于刚性翼型开展的水动力性能计算与分析，然而刚性翼型由于无法变形，也就无从研究双翼结构力学方面的性能，特别是不能提供双翼内部应力分布情况，而这对于双翼推进器疲劳性能的研究是至关重要的.下面通过数值模拟方法，对上述问题进行计算分析.

通过反复在双翼上施加上下的沉浮运动，并限制双翼绕自身的转动，就可以依靠双翼的旋涡发放所产生的涡激力推动双翼前进，图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，弹性模量相关系数为0.70×106和0.35×106；有机材料硬度为40IRHD的橡胶，泊松比为0.499，密度为1 300 kg/m3，弹性模量相关系数为0.275×106和0.028×106.由于橡胶的种类繁多，资料中弹性模量通常不直接给出，而是给出相关系数，可以直接输入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种不同类型的材料组成的柔性摆动双翼推进过程进行数值仿真.

由于翼型结构内部的应力分布与结构约束以及动力施加位置有关，需要进行相关说明.本文中，摆翼的约束位置也是外部动力的施加位置，它在靠近翼型头部，翼型弦长的1/3处.改变约束位置和外部动力的施加位置，将会引起翼型内部应力分布的变化，详细的变化规律将在以后的论文中进行研究，本文基于上述固定的约束和外力位置进行讨论.

3.1 材料类型对双翼水动力性能影响的数值模拟计算与分析

由于双翼是关于对称轴对称的，为了更加清晰，只选取了双翼中的一个翼型给出数值模拟图像.受到计算软件当前版本的功能限制，对于含固体变形的流固耦合计算，暂时还不能直接输出涡量分布图，但可以提供流体的速度矢量云图作为替代，以便讨论材料类型对翼型后方发放的旋涡形状的影响，如图4所示.表1为4种柔性摆动双翼的水动力性能计算结果.

图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表明，由于有机材料的弹性模量比金属材料的弹性模量要小，有机材料组成的柔性摆动双翼在推进过程中与流体相互作用时比较容易发生变形，这种变形反过来又影响到旋涡运动的轨迹，因此从流体的速度云图上可以看出，金属双翼发放的旋涡轨迹明显有别于有机材料组成的柔性摆动双翼发放的旋涡轨迹；另外，聚乙烯双翼、60IRHD橡胶双翼和40IRHD橡胶双翼之间材料弹性模量的差别并不大，所以它们的尾部旋涡轨迹非常相似，导致它们水动力性能计算结果的差别也不大，但是都优于金属双翼的计算结果.

上述计算结果是在St为0.4时得到的，还不够全面.为了进一步研究材料对双翼水动力性能的影响，在St从0.2到0.4时，分别对由金属材料组成的双翼和由40IRHD橡胶材料组成的双翼进行数值计算，它们的推力系数曲线和推进效率曲线如图5所示.

图5 两种不同材料翼型的水动力性能曲线Fig.5 Hydrodynamic performance curves of foils with two kinds of materials

从图中可以看出，第一，在较低的St下，随着St的增加，无论是金属双翼还是橡胶双翼的水动力性能曲线都呈上升的趋势，也就是说双翼的推力系数和推进效率都逐渐增大；第二，在相同的St下，橡胶双翼的水动力性能明显优于金属双翼.从两个方面来进行分析：在双翼变形方面，图4中金属双翼的最大有效应变为1.00×10－7、聚乙烯双翼的最大有效应变为1.02×10－6、60IRHD橡胶双翼的最大有效应变为1.84×10－5、40IRHD橡胶双翼的最大有效应变为2.39×10－5，可以发现有机材料双翼的变形量比金属双翼更大；在双翼尾流方面，图4中金属双翼只有2个有标记的旋涡，而有机材料双翼有3个有标记的旋涡，这说明有机双翼的变形对双翼后方的旋涡发放有较大的影响.因此，橡胶双翼的水动力性能明显优于金属双翼的原因在于，橡胶双翼在流固耦合过程中发生了明显变形，从而改变了翼型的尾涡运动轨迹，也就改变了旋涡发放产生的涡激力的大小和方向，最终改进了双翼的水动力性能.

3.2 材料类型对双翼结构力学性能影响的数值模拟计算与分析

上述计算在输出双翼水动力性能数据的同时，也输出了对应的结构力学性能数据.图6为不同双翼的内部应力云图，主要展示的是双翼变形程度和应力集中的位置；表2为双翼推进过程中，提取出的内部最大有效应力数据，主要展示的是材料类型对结构内部应力的改变程度.

图6 不同材料双翼有效应力云图Fig.6 Effective stress contours of dual－foils with different materials

表2 不同材料双翼最大有效应力Tab.2 Maximum effective stress of dual－foils with different materials

因为不显示周围的流体，并且将翼型进行了放大，所以图6比图4更加清晰地反映了双翼的变形情况，由于材料的弹性模量各不相同，金属双翼在整个推进过程中几乎不发生变形，而有机双翼都有不同程度的形状改变；另外，应力云图主要体现的是双翼在工作时内部的应力分布情况，可以看出4种双翼都出现了应力集中的现象，而且应力集中的区域基本上都在双翼的中部.这是因为双翼的沉浮运动荷载加载在靠近双翼头部的1/3处，而双翼的旋转约束也在该处，所以双翼运动时旋涡发放产生的涡激力最终在靠近约束的位置产生了应力集中.表2定量计算了4种双翼应力集中的程度，其中，金属双翼的内部最大有效应力远远大于有机材料组成的柔性摆动双翼，而在有机双翼中，40IRHD橡胶的内部最大有效应力最低.因此，在外部条件相同的情况下，材料的柔性是双翼内部最大有效应力的主要影响因素.

进一步的计算给出了40IRHD橡胶材料组成的双翼，在较低的St下，即St从0.2到0.4时，计算得到的最大有效应力（σe，max）曲线和最大剪切应力（τmax）曲线（如图7所示）.随着St的增加，橡胶双翼的最大有效应力和最大剪切应力都会变大，它们代表的是整个翼型各个截面当中应力的最大值.

图7 40IRHD橡胶双翼最大应力曲线Fig.7 Maximum stress curves of rubber dual－foils with 40IRHD

4 结 论

（1）有机材料能改善双翼水动力性能.与金属双翼相比，有机双翼的弹性模量较低，翼型结构更柔软，更加容易受到旋涡运动的影响产生较大的变形，这种变形反过来又改变了双翼尾流区域旋涡的形状，调整了涡激力的大小和方向，从而提高了双翼涡动力推进器的推进性能.

（2）有机材料能提高双翼抗疲劳破坏能力.在双翼沉浮运动过程中，双翼在周期性流体力的反复作用下会出现应力集中现象，翼型结构交变应力集中区域位于双翼的中部.有机双翼结构比较柔软而且容易变形，正是这种变形缓解了应力集中，最后降低了结构受到的最大交变应力值，这能够优化柔性摆动双翼的抗疲劳性能.

综上所述，有机材料对双翼性能的影响是有利的，选择合适的有机材料制造翼型结构，不仅可以提高推进性能，也能够延长其使用寿命.

[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－hsin.Numerical analysis on aerodynamic force generation of biplane counterflapping flexible airfoils [J].Journal of Aircraft，2009，46（5）：1785－1794.

[5] 周 岱，马 骏，李 磊.流场和流固耦合问题网格剖分与更新的新方法[J].工程力学，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].水动力学研究与进展，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基础理论与工程实践[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.