多个柔性梁V 型集群行为节能机理研究1)

张磊 敖雷 裴志勇

(武汉理工大学绿色智能江海直达船舶与邮轮游艇研究中心,武汉 430063)

(绿色智能江海直达船舶湖北省工程研究中心,武汉 430063)

引言

集群是动物较为普遍的行为,如大多数鸟类和鱼类都会阶段性地集群,这种聚集行为产生的有利作用称为集群效应.生物学[1-4]从动物对资源需求和环境适应等方面阐述集群行为产生的原因,指出集群可以提高族群的防御力、捕食效率和繁殖率.然而,这些定性解释并不能满足人类对动物集群的认知需求,例如:迁徙过程中为什么大雁会阶段性地使用一字形排列的集群方式;为什么大雁等体型较大的鸟类常采用V 型集群方式,而小型鸟类不采用.亟需定量研究动物集群行为,提升对这一自然现象的认知水平.

近几十年国内外学者尝试从受力和能量的角度对动物集群行为和集群效应的内在机理进行定量研究,采用的方法包括理论分析、试验(活体试验和仿生学试验)研究和数值分析.首先,空气动力学和流体动力学的相关理论被用来研究V 型集群方式下鸟群个体之间和菱形排列鱼群个体之间的相对运动,指出集群中跟随运动的个体可以利用前面个体产生的涡街从而节省能耗[5-6].然后,这种节能观点也得到了试验验证.对V 型排列8 只鹈鹕的心率进行实时监测发现,群体飞行中个体心率较独自飞行的鹈鹕减少了11.4%～14.5%[7];采用全球定位系统(GPS)对14 只朱鹮的飞行轨迹进行记录发现,虽然个体位置不断调整,但群体呈现出比较完整的V 型排列[8];测量鲑鳟鱼在均匀流和圆柱尾流后的运动发现,鱼可以利用圆柱尾流涡降低自身肌肉的参与度[9].由于活体试验变量较难控制,且采集的数据(如心率和肌肉参与度等)只是集群节能观点的侧面证据,亟需能够表征受力和能耗的直接数据.学者们开展了大量的仿生学试验[10-17],根据仿生简化方式的不同将这些试验分为两类,一类是多个柔性体流致振动的研究,如交错排列的多个柔性体中领头柔性梁的受力较单独流致振动柔性梁的受力减少可以达到50%[13],本类研究没有考虑柔性梁自主摆动的推进作用;第二类是多刚性体集群行为的研究,如纵向排列2 个刚性翼采用转动或者升沉方式推进,根据不同的初始纵向间距会自然形成不同的集群模式,且合适的纵向间距下刚性翼的推进效率较单独推进的刚性翼高[15,17],本类研究忽略了动物柔性的影响.为了更加准确地描述动物集群行为,需要开展多个自推进柔性体自形成集群运动的仿生学试验,受制于目前的试验条件,实施起来难度较大.

相对于试验而言,数值方法求解动物集群这种多个自推进柔性体−流体介质相互作用问题更为经济,得到的流场细节有利于开展机理分析.文献[18-19]采用浸入边界法对均匀流中并排2 个柔性梁的流致振动进行了模拟,本文作者基于Fluent 二次开发提出了流固耦合数值算法,对均匀流中并排2 个柔性梁的流致振动进行了模拟[20],研究均发现不同的横向间距下柔性梁呈现同向和异向振动两种模式.同样,浸入边界法被用于均匀流中纵向排列2 个柔性梁的流致振动计算[21],后梁受力和幅值较前梁大.文献[14]采用面元法和振型叠加的联合仿真算法对均匀流中不同排列(并排、纵向或交错) 2 个柔性梁的流致振动问题进行讨论,得出了不同排列方式下柔性体受力的相关数据.上述研究中柔性体固定在均匀流中,受到流体力产生被动变形,而自然界中鸟类或鱼类不仅具有被动变形,它们还通过主动控制自身肌肉的形变(振翅或摆尾等)获得动力.文献[22]在流致振动前梁形成的尾流中引入自推进的柔性梁,采用浸入边界法进行计算发现,后梁在前梁尾流中的运动呈现绕涡模式和穿涡模式,分别对应低阻力和高阻力的状态.文献[23-24]采用浸入边界法对纵向多个柔性梁自推进的流固耦合问题进行求解,研究中柔性梁导边做同振幅、同频率的横向指定运动,最终形成了紧密排列集群模式,多个柔性梁群体较单独自推进梁具有更大的速度.同样,浸入边界法还用于并排2 个柔性体自推进的研究[25-26],根据横向间距的不同呈现交替跟随、交替领先和并驾齐驱的集群模式,其中交替跟随模式下柔性梁的推进速度和效率较高.本文作者采用基于Fluent 二次开发的流固耦合数值算法对并排排列多个(小于4)自推进柔性梁进行了仿真[27],研究中的并排集群方式并不能降低能耗.

对以上案例进行分析可知,首先,浸入边界法作为一种流固耦合边界处理方法,主要通过与特定的计算流体力学方法(如格子玻尔兹曼方法)结合用于一些基础研究,计算流体力学软件中尚未接入浸入边界法的相关模块,使用该方法解决复杂科学问题或者工程问题的门槛较高;然后,动物集群流固耦合问题的相关研究方兴未艾,集群模式也多集中在纵向或并排排列的2 个柔性体,缺少对多个柔性梁(大于2) V 型集群行为的研究.鉴于此,本文研究具有两个目的,一是采用基于Fluent 二次开发的流固耦合数值方法对多个自推进柔性梁自形成的V 型集群过程进行模拟;二是对V 型集群中柔性梁个体的受力和运动参数进行讨论,结合流场细节,揭示集群产生的原因和节能的内在机理,特别是对领头柔性梁的节能机理进行阐述.

1 数值模型和基本原理

1.1 V 型集群数值模型

目前对动物集群行为多进行二维多个柔性体−流体介质相互作用流固耦合问题的研究[23-27],动物个体近似表达为二维柔性梁.多个柔性梁V 型排列集群计算模型如图1 所示,x轴和y轴分别表示水平方向和垂直方向,多个二维柔性梁(即结构域Ωs)置于流体域Ωf中.本文柔性梁选取的最大数量为5,长为L,厚度为h.由于h≪L,柔性梁可看作欧拉伯努利梁.柔性梁左端(即导边Γl)固定并在y方向做指定的升沉运动y(t)=Acos(2πft),然后该升沉运动与流体的相互作用提供柔性梁水平方向自由运动的动力,最终形成稳定的V 型排列.多个柔性梁同相位升沉运动的幅值A和频率f保持一致.导边横向(即y方向上)间距设置为D,运动过程中保持不变.柔性梁i和柔性梁j导边纵向(即x方向上) 间距用Hij(t)表示,i,j=1,2,3,4,5.纵向间距初始值设置为Hij(t=0)=H0,在运动过程中纵向间距Hij(t)是随时间变化的,由多个柔性体与流体之间的相互作用决定.

图1 多个柔性体自推进V 型排列计算模型示意图Fig.1 Schematic of multiple self-propelled flexible beams in V configuration

1.2 流固耦合问题控制方程

不可压缩黏性流动力学控制方程包括质量连续性方程和动量方程,分别为

其中,uf是流体速度矢量,ρf是流体密度,pf表示压力,μ是动力黏性系数,gf是外部体积力.为了求解流体动力学方程,还需要补充相应的初始条件和边界条件,分别为

其中,式(2a)表示流体在初始时刻保持静止,式(2b)表示流固耦合面Γi=Ωf∩Ωs上的速度边界由柔性体的运动ui决定,式(2c)表示流体域边界Γw为无滑移边界条件.流体动力学方程采用有限体积法进行求解.

欧拉伯努利梁在x和y方向上的动力学控制方程分别为

其中,二维柔性梁线密度ρl=ρsh,ρs为柔性体密度,柔性梁位置矢量X(s,t)=(X(s,t),Y(s,t)),流体外力矢量F(s,t)=(Fx(s,t),Fy(s,t)),s表示沿着柔性梁轴线的拉格朗日坐标,柔性梁的杨氏模量为E,转动惯量为I.运动过程中认为柔性梁自身长度保持不变,补充几何约束条件如下

同样地,为了求解柔性梁的动力学方程还需要补充合适的初始条件和边界条件为

其中,式(4a)和(4b)分别表示初始时刻柔性梁的位置和运动状态,本文柔性梁初始时刻无形变且保持静止.式(4c)表示柔性梁导边Γl做指定的升沉运动,yi是升沉运动平横位置在oxy坐标系中的坐标值(见图1).式(4d)表示流体作用在流固耦合面上的力pf(s,t),nf是流体力单位矢量,指向流体外,该流体力通过求解流体动力学方程获得.柔性梁x轴方向上的动力学方程为二阶微分方程.柔性梁导边在y轴方向上做无相位差的升沉运动y(t)=Acos(2πft),即可将研究的柔性梁看作是导边具有相同升沉运动激励的悬臂梁,其y轴方向上的动力学方程可采用广义支座激励的振型叠加法进行处理.首先,将悬臂处支座激励看作是等效载荷的作用,等效载荷处理方法参见文献[28] 第17 章;然后,由于密度、尺寸、杨氏模量和转动惯量等参数都相同,则各悬臂梁振型模态都是一致的.柔性梁的微分方程采用4 阶龙格库塔法进行求解.

上述两场(即流场和结构场)通过耦合面Γi相互作用,从式(3)、式(2b)和式(4d)可以看出,流场产生的作用于结构物上的流体力将改变结构场的动力学条件,反过来结构场运动状态(即速度)的改变提供影响流场的边界条件,流固耦合计算中需要保持耦合面上两场速度和受力的连续性,表达式分别为

其中,ns是结构力单位矢量,指向结构外.根据如式(5a)耦合面上的运动采用动网格技术更新计算域.

1.3 松耦合方法

采用分区算法中的松耦合方法对多个柔性体−流体介质相互作用流固耦合问题进行求解.在单个时间步内,分别对流场和结构场进行一次计算,然后两场之间数据也只交换一次,具体计算流程如图2所示,ϕ为表示任意流体参量的通用变量.从时间步tn到tn+1的具体计算流程叙述如下:

图2 松耦合方法流程图Fig.2 Flow chart of the proposed loosely-coupled scheme

(1)时间步tn已知参量uf=uf,n,X(s,t)=X(s,tn),∂X(s,t)/∂t=∂X(s,tn)/∂t;

(2)时间步tn到tn+1:

(a)求解流体动力学方程积分获得作用在柔性梁上的流体力pf,n;

(b)求解结构动力学方程获得tn+1 时刻柔性体的位置矢量X(s,tn+1);

(c)根据结构物位置矢量采用动网格技术更新流体计算域.

2 验证实例

为了验证本文数值方法的有效性,对文献[23]中纵向2 个柔性体自推进的集群行为进行模拟.柔性体导边做同相位的指定升沉运动y(t)=Acos(2πft),选择特征长度L、速度Uref=2πAf和流体密度ρf无因次化参量,则柔性体的线密度无因次量表示为=ρl/(ρfL),雷诺数为Re=ρfUrefL/μ,弯曲刚度为K=选择=0.2,Re=200,K=0.8,A/L=0.5,H0/L=8 的工况进行计算.H0表示2 个柔性体之间的初始纵向间距.为了消除流体域无滑移壁面Γw阻塞效应的影响,选择足够大的流体计算域,x和y方向上尺寸分别为[−60L,30L] × [−20L,20L].为了减少网格量采用混合网格对流体域进行空间离散,如图3 所示.为了准确地获取柔性梁附近区域的详细流场信息,网格需要进行局部加密,图3 中红色方框内局部放大图显示了柔性梁周围非结构网格加密情况,其中最小网格边长为∆Lmin=0.02L.流场和结构场的计算时间步长选择统一值∆t=0.01 s(约为0.001T),T=1/f为柔性梁指定升沉运动的周期.

图3 计算域空间离散示意图Fig.3 Spatial discretization of the fluid domain

图4 将计算得到的2 个柔性梁速度和纵向间距的时历曲线与浸入边界法的结果进行对比可知,两种方法计算得到的曲线吻合良好.在t/T=0 至2 区间内,纵向排列2 个柔性梁保持相同的速度水平移动,接着后梁突然加速,前后梁的纵向间距逐渐变小直到最后达到一个稳定的间距((H−L)/L=5).

图4 纵向排列柔性梁导边运动参数时历曲线对比:(a) 水平运动速度;(b) 纵向间距Fig.4 Comparison of time history of (a) streamwise velocity and(b) longitudinal separation distance of two tandem beams between the present solutions and Ref.[23]

此外,对数值计算的网格和时间步长无关性进行说明.由于柔性梁动力学方程求解中使用了振型叠加的处理方式,还需要开展模态截断分析(即模态无关性分析).另外选择两套空间离散网格(最小网格边长分别为∆Lmin=0.01L和∆Lmin=0.04L)、两个时间计算步长(∆t=0.002 s 和∆t=0.05 s)和两个振型模态数(jmax=3 和jmax=10)进行计算.图5 对比了集群运动稳定后的柔性梁水平运动速度,仿真计算结果随着网格尺寸和时间步长的减小单调收敛,随着柔性梁模态数的增加单调收敛.同时,当网格尺寸选择∆Lmin=0.02L,时间步长∆t=0.01 s,模态数jmax=5 时,仿真结果满足计算的精度需求,在后面V 型集群行为的模拟中采用相同的设置.

图5 柔性梁导边水平运动速度时历曲线对比:(a) 网格无关性分析;(b) 时间步长无关性分析;(c) 模态截断分析Fig.5 Comparison of streamwise velocity curves of the leading edge in(a) grid independence study,(b) time-step independence study,and(c) mode truncation study

图5 柔性梁导边水平运动速度时历曲线对比:(a) 网格无关性分析;(b) 时间步长无关性分析;(c) 模态截断分析(续)Fig.5 Comparison of streamwise velocity curves of the leading edge in(a) grid independence study,(b) time-step independence study,and(c) mode truncation study (continued)

3 数值模拟结果及分析

对多个柔性体V 型排列进行仿真,用来模拟自然界中动物V 型集群行为.文献[29-30]中指出自然界中动物的弯曲刚度系数近似为O(1),这里仍然选用第三节验证案例中的柔性体,参数为=0.2,Re=200,K=0.8,A/L=0.5,柔性体的数量分别为3 个和5 个.

3.1 3 个柔性梁的稳定V 型排列

3 个柔性梁横向间距设置为D=0.475L,初始纵向间距设置为H12(t=0)=H13(t=0)=H0=2L.展示了V 型排列3 个自推进柔性梁之间水平速度差值Uij(t)=ui(t)−uj(t)和纵向间距Hij(t)=Xi(s=0,t)−Xj(s=0,t)的时历曲线,其中ui(t)和Xi(s=0,t)分别代表柔性梁导边的水平速度和位置.领头柔性梁1 与柔性梁2 和3 的水平速度差值在t/T=4 后逐渐趋于稳定,绕着零点振荡,幅值为0.5Uref.柔性梁2 和3 与领头柔性梁1 之间的纵向间距先增大然后快速减小,t/T=4 后基本保持一致且位于[−0.91L,−1.16L]的区间内作稳定的余弦变化.

图7 展示了计算稳定后一个升沉运动周期内典型时刻3 个柔性体的集群行为以及相应的涡量和压力云图.如图7(a)～图7(d)所示,y方向上,柔性梁导边从其升沉运动的平衡位置(如柔性梁1,Y1(s=0,t=0)=0)向最大负幅值运动(如柔性梁1,=−A,然后向上运动再次经过平衡位置直到t/T=1 时达到最大正幅值处.x方向上,柔性梁2 和3 齐头并进紧随领头柔性梁1 后,这一现象与图6(b)红色方框内放大的纵向间距时历曲线一致.观察形成的涡量图,其中涡根据相应梁的编号来命名,涡量值用特征值ωref=Uref/L无因次化,正负号(即 +,-)分别表示涡绕逆时针和顺时针旋转.t/T=1/4 时刻3 个柔性梁上表面导边处边界层分离形成的负涡1−,2−和3−向x正向运动,随后在t/T=3/4 时脱离柔性梁在尾流后混合形成负的卡门涡群(1−+3−和2−).同理,t/T=3/4 时刻3 个柔性梁下表面导边处边界层分离形成的正涡1 +,2+和3+向x正向运动,随后在t/T=1+1/4 时脱离柔性梁在尾流后混合形成正的卡门涡群(1++2+和3 +).由此可见,一个完整升沉运动周期内3 个柔性梁形成的V 型集群模式呈现出一对完整的卡门涡街群.附录视频1 展示了3 个柔性体V 型集群一个周期内完整的涡量变化.图7(e)～图7(h)给出相应典型时刻的流场压力云图,远场参考压力设置为零,特征压力用来无因次化计算得到的压力值,红色和蓝色分别表示最大正压力值2 和最小负压力值−2.为了弄清3 个自推进柔性梁自形成V 型集群行为产生的原因,尝试结合涡量云图和流场压力云图给出解释.边界层分离以及涡产生的地方形成低压区,例如图7(a)中,t/T=0/4 时刻3 个柔性体上表面处边界层发生分离,同时柔性体1 下表面分离的边界层脱离其尾部,这四个区域如图7(e)所示均为负压区,而柔性体下表面具有较大的正压力,另外,该时刻柔性体从升沉运动平衡位置向下运动产生了向上变形,柔性体上下表面压力差在水平方向上形成了沿着x负方向的合力,与柔性体水平运动方向相同,该合力推动柔性体加速运动.t/T=1/4 时刻柔性体达到升沉运动的最大负振幅位置处开始向上运动,此时柔性体具有向下的变形(见图7(f)),结合上表面的负压和下表面的较大的压力形成了向着x正方向的合力,该合力与柔性体水平运动方向相反,该阻力使得柔性体减速运动.采用同样的方法也可以对t/T=1/2 和t/T=3/4 时刻的涡量云图和压力云图进行分析.由以上分析可知,自推进多个柔性梁之间通过流体介质相互作用(如涡的相互混合)产生不同的流场分布,不同的压力分布和柔性梁形变共同决定了柔性梁的水平受力,从而影响多个柔性体之间的相互运动.图8 给出了一个升沉运动周期内3 个柔性梁水平受力Fx的时历曲线,Fx<0 时水平力为推力,柔性梁加速运动;反之,柔性梁作减速运动.后排柔性梁2 和3 受到的推力和阻力幅值较柔性梁1 大.

图6 3 个柔性梁V 型排列运动参数时历曲线:(a) 速度差值Uij(t)=ui(t)−uj(t);(b) 纵向间距Hij(t)=Xi(s=0,t)−Xj(s=0,t) (i,j=1,2,3)Fig.6 Time history of (a) velocity difference Uij(t)=ui(t)−uj(t),and(b) longitudinal separation distance Hij(t)=Xi(s=0,t)−Xj(s=0,t) (i,j=1,2,3) of three beams in V configuration

图7 3 个柔性梁V 型排列流场细节图:(a～d) 涡量云图;(e～h) 压力云图.特征值ωref=Uref/L 和pref=ρf 分别被用来无因次化计算得到的涡量和压力Fig.7 Fluid details:(a～d) vortices contour and (e～h) pressure contour of three beams in V configuration.The vortices and pressures are normalized by ωref=Uref/L,pref=ρf ,respectively

图8 3 个柔性梁水平受力时历曲线,特征值Fref=(1/2)ρf 被用来无因次化计算得到的合力Fig.8 Hydrodynamic force experienced by three flexible beams in the horizontal direction during one heaving motion period.The forces are normalized by Fref=(1/2)ρf

3.2 5 个柔性梁的稳定V 型排列

5 个柔性梁横向间距为D=0.475L,初始纵向间距为H12(t=0)=H13(t=0)=H24(t=0)=H35(t=0)=H0=2L.图9 显示了V 型排列5 个自推进柔性梁之间水平速度差值和纵向间距的时历曲线.t/T=4 后,第二排的柔性体2 和3 与领头柔性体1 之间的速度差值以及第三排柔性体4 和5 与第二排柔性体2 和3 的速度差值逐渐趋于周期性变化,变化周期为3T.同时如图9(c)和9(d)所示,V 型排列三排柔性体之间的纵向间距时历曲线也呈周期为3T规律性波动,第二排与领头柔性体之间的纵向间距位于[−0.98L,−1.76L]区间内,第三排与第二排柔性体之间的纵向间距较小,位于[−0.45L,−1.03L]区间内.这说明5 个柔性体在自推进集群运动过程中虽然整体一直保持V 型排列,但个体之间的相互位置不断发生着调整,该结论与文献[8]中对14 只朱鹮V 型集群的动物活体试验结果得出的结论相符合.

图10 选取展示了t/T=5.5～8.5 内典型时刻5 个柔性体的集群行为以及相应的涡量云图.如图10(a)～图10(d)所示,y方向上,柔性梁导边从其升沉运动的最大负幅值(t/T=5.5)向上运动经过平衡位置(t/T=5.75)达到最大正幅值(t/T=6.0)处,然后向下运动经过平衡位置(t/T=6.25),图10(e)～图10(h)以及图10(i)～图10(l) 柔性体导边依次重复上述运动.x方向上,第二排柔性体与领头柔性体纵向间距先变大而后变小,而第三排柔性体与第二排柔性体纵向间距经历了先变小然后变大的过程,这一现象与图9(c)和9(d)中纵向间距时历曲线一致.观察5 个柔性体V 型集群运动产生的涡量图可以看出,一个升沉运动周期T内,柔性体上(下)表面边界层分离依次产生一对负(正)涡,涡在向尾流移动过程中受到多个柔性体之间通过流体介质产生的相互作用而挤压混合形成了复杂的涡量云图.对比3 个升沉周期相同时刻(如图10(a),图10(e)和图10(i))的涡量图发现,涡相互作用的模式基本相似,但由于多个柔性体之间相互位置的差别,造成涡受挤压或者混合位置以及大小都不尽相同,这些差别将会导致柔性体受力的不同.5 个柔性体V 型集群运动在三个指定升沉运动周期内完整的涡量变化详见附录视频2.

图9 5 个柔性梁V 型排列运动参数时历曲线:(a,b) 速度差值Uij(t)=ui(t)−uj(t);(c,d) 纵向间距Hij(t)=Xi(s=0,t)−Xj(s=0,t) (i,j=1,2,3,4,5)Fig.9 Time history of (a,b) velocity difference Uij(t)=ui(t)−uj(t),and(c,d) longitudinal separation distance Hij(t)=Xi(s=0,t)−Xj(s=0,t)(i,j=1,2,3,4,5) of five beams in V configuration

图10 5 个柔性梁V 型排列流场涡量云图.特征值ωref=Uref/L 被用来无因次化计算得到的涡量Fig.10 Vortices contour of five beams in V configuration.The vortices are normalized by ωref=Uref/L

图11(a)～图11(l)给出了相应时刻的流场压力云图,根据3.1 节对压力场的分析方法,边界层分离以及涡产生的地方形成低压区,结合柔性体的形变,可以对作用在柔性体上流体力进行说明,这里不再详述.图12 给出了5 个柔性体水平方向受力的时历曲线,第三排柔性体4 和5 受到的推力(Fx<0)幅值最大,第二排柔性体2 和3 次之,领头柔性体最小;领头柔性体受到的阻力(Fx>0)幅值最小,后排柔性体较大.

图11 5 个柔性梁V 型排列流场压力云图.特征值pref=ρf 被用来无因次化计算得到的压力Fig.11 Pressure contour of five beams in V configuration.The pressures are normalized by pref=ρf

图12 5 个柔性梁水平受力时历曲线,特征值Fref=(1/2)ρf L 被用来无因次化计算得到的合力Fig.12 Hydrodynamic force experienced by five flexible beams in the horizontal direction during three heaving motion periods.The forces are normalized by Fref=(1/2)ρf L

3.3 V 型集群性能分析

上述内容对多(3 和5)个柔性体的V 型集群运动进行了描述,接下来引入三个主要参数对多柔性体自推进集群运动性能进行讨论[23-25].集群运动水平方向平均推进速度Ui,平均输入功率Pi和推进效率ηi分别为

其中,下标i表示柔性体的编号,Tc表示多柔性体在升沉运动自推进作用下形成的稳定集群运动周期.观察图6 和图9 的运动参数时历曲线可知,3 个和5 个柔性体V 型集群水平速度差值和纵向间距曲线分别呈现T和3T的周期性变化,即3 个和5 个柔性体形成的V 型集群运动周期分别为T和3T.需要注意的是,式(6b)中平均输入功率包括两部分,一部分产生供应升沉运动需要的能量,一部分提供水平运动的能量.引入文献[23]中单柔性体自推进运动的数据进行对比,结果见表1.在多个柔性体形成稳定的集群运动后个体水平方向平均推进速度保持一致,即集群具有一个固定的平均推进速度(U=Ui).从表1 可以看出,在柔性体输入功率基本相同的情况下,V 型排列集群中个体自推进获得的水平方向速度较单柔性体提升比例超过15%,虽然5 个柔性体集群中第二排柔性体2 和3 需要更大的输入功率,但是其推进效率较单个柔性体提升比例超过14.5%.该对比包含两方面的内容:(1)在相同能量的输入下,集群中个体可以获得更高的推进速度;(2)即使集群中某些个体存在高能量输入的情况,但其推进的效率仍然较高.本文的计算结果给出了V 型排列集群行为具有降低能耗的直接数据证据.此外值得注意的是,本文研究的V 型集群中领头柔性梁的运动速度和效率均较单独柔性梁高,且是集群中效率最高的个体.这一结论打破了集群中领头“羊”位置最耗能的认知.

表1 多个自推进柔性体性能参数比较Table 1 Propulsive performance of multiple self-propelled flexible beams

接下来尝试对V 型集群行为节能机理进行阐述.文献[23-25,31]等对集群中后排柔性体节能的原因做出了解释,认为后排柔性体可以利用前排柔性体尾流中的涡,从本文3.1 节和3.2 节对集群流场(涡量云图和压力云图)的分析可以得出同样的结论,即前排柔性体尾涡在后排柔性体前端形成了低压区,造成柔性体前后具有更大的压差(如图10(f)和图10(j),图11(f)和图11(j)中柔性体3 周围的涡量图和压力场所示),从而为水平运动提供更大的推力.

这里需要重点对领头柔性梁节能的机理进行阐述,图13(a)对比了领头柔性梁的水平受力时历曲线,图13(b)给出了领头柔性体尾部和导边横向位移差值(即横向形变,Δy=Y(s=L,t)−Y(s=0,t))时历曲线的对比,可以看出,V 型集群领头柔性梁水平方向受力的幅值较单柔性体大,且横向形变幅值较小.图14分别给出了单柔性体和3 柔性体V 型集群的流场细节图.t/T=5.25 时刻,3 柔性体V 型集群运动领头柔性梁尾部脱落的涡受到第二排柔性梁2 的限制,沿着柔性梁2 的下表面运动,与单柔性体压力云图对比可知,其尾涡形成的低压区下移(如图14(d) 和14(e)中红色方框所示),且领头柔性梁的尾部上表面处压力较大,这导致t/T=5.25 时刻集群中领头柔性体的最大横向形变较单柔性体小.t/T=5.25 和t/T=5.3 时刻分别是单柔性体受到推力和3 个柔性体中领头柔性体受到推力最大的时刻,对比图14(a),图14(d)和图14(c),图14(f)可以看出,由于领头柔性体尾涡与第二排柔性梁2 和3 导边边界层分离涡的相互作用,集群领头柔性体上下表面压差增大,所以虽然t/T=5.3 时刻领头柔性体的横向形变较小,但是其水平方向的推力仍然比t/T=5.25 时刻单柔性体的推力略大.从上述分析得出,集群中领头柔性体与后排柔性体通过流体介质发生相互作用(如领头柔性梁尾涡与后排柔性梁边界层分离涡相互作用造成流场变化),一方面使得柔性体横向形变幅值变小,另一方面造成柔性体上下表面压力差变大.在这两方面因素(更小的迎水面和更大的推力)的共同作用下,领头柔性体获得较高的运动速度和推进效率.

图13 领头柔性梁对比:(a) 水平受力时历曲线;(b) 横向形变时历曲线.特征值Fref=(1/2)ρf L 被用来无因次化计算得到的合力Fig.13 The comparison of the leading beams:(a) hydrodynamic force in the horizontal direction and (b) lateral deformation during one heaving motion period.The forces are normalized by Fref=(1/2)ρf L

图14 单柔性体和3 个柔性梁流场细节图:(a～c) 涡量云图;(d～f) 压力云图.特征值ωref=Uref/L 和pref=ρf 分别被用来无因次化计算得到的涡量和压力Fig.14 Fluid details of single beam and 3 beams:(a-c) vortices contour and (d-f) pressure contour.The vortices and pressures are normalized by ωref=Uref/L,pref=ρf ,respectively

4 结论

本文采用基于Fluent 二次开发的数值方法求解多个柔性体−流体介质相互作用的V 型集群问题.首先,使用该算法对纵向2 个柔性体的二维集群问题进行计算,将结果与浸入边界法仿真结果进行对比说明松耦合方法的可靠性.然后,对自然界中常见的动物V 型集群行为进行模拟,分别计算了3 柔性体和5 柔性体升沉运动自推进的二维集群行为.在4 个升沉运动周期(即t/T=4)后,多个柔性体都形成了稳定的V 型排列,柔性体间横向间距保持一定,纵向间距在一定区域范围内呈现周期性变化,如5 柔性体V 型集群中第二排柔性体与领头柔性体纵向间距稳定在[−0.98L,−1.76L]区间内.该计算结果与文献[8]中朱鹮V 型集群的动物活体试验结论一致,即动物在自推进集群运动过程中虽然整体一直保持V 型排列,但个体之间的相互位置不断发生着调整.引入水平方向平均推进速度、输入功率和推进效率三个参数对集群中个体推进性能进行对比可知,不仅后排柔性梁的速度和推进效率得到提升,领头柔性梁也得到大幅提升,其中推进速度较单独柔性体自推进速度提升比例超过15%,推进效率提升比例也超过14.5%,该计算结果给出了V 型排列集群行为降低能耗的定量直接证据.值得注意的是,该V 型集群运动中领头柔性梁是集群中速度和效率最高的个体.这一结果打破了集群中领头“羊”位置最耗能的普遍认知.

为了弄清集群行为产生的原因和不同位置个体的节能机理,对数值仿真得到的二维集群流场细节进行了详细分析.多个柔性体通过流体介质发生相互作用(如边界层分离形成的涡相互挤压混合)造成压力场的变化,在此压力场和柔性体自身形变的共同作用下,多个柔性体逐渐形成了稳定的V 型集群排列.V 型集群中不同位置的个体采用不同策略提升推进性能.后排柔性梁节能机理与其他研究者给出的结论一样,即其可利用前排个体的尾流涡进而节省能耗;本文还对集群中领头柔性梁的节能机理进行了阐述,即在自身产生的尾涡与后排柔性体的相互作用下,领头柔性梁自身形变变小和所受推力幅值变大.

为了使多个柔性体的集群运动更加接近自然界中丰富多样的动物集群现象,研究具有不同长度和弯曲刚度的多个柔性梁采用不同指定运动(升沉、转动或两者相结合)方式自推进进而自发形成的集群运动行为将是下一步的工作方向.