ENERGY SAVING MECHANISM OF HYDRODYNAMIC COLLECTIVE BEHAVIOR OF MULTIPLE FLEXIBLE BEAMS IN V FORMATION
-
摘要: V型排列是自然界中常见的动物(如大雁等迁徙性鸟类)集群模式, 普遍推测该模式可以有效地降低能耗, 然而目前没有研究给出相关的直接数据证据. 开展其节能机理研究有助于提升对集群自然现象的认知水平, 为集群仿生应用打下基础. 本文采用基于Fluent二次开发的数值方法求解多个柔性体−流体介质相互作用的流固耦合问题, 其中流体动力学方程采用有限体积法进行求解, 柔性体动力学控制方程通过用户自定义模块(UDF)嵌入, 并采用模态叠加法和4阶龙格库塔法求解, 流固交界面形变使用动网格技术处理. 实现了多个自推进二维柔性梁自主形成V型集群运动过程的数值模拟, 并将得到的推进性能参数(平均速度, 输入功率和效率)与单独自推进柔性体的数据进行对比. 研究发现: 该V型集群运动中不仅后排柔性梁的速度和推进效率得到提升, 领头柔性梁性能也得到大幅提升, 增幅均超过14%. 此外, 对V型集群运动的流场细节(涡量和压力云图)开展分析, 揭示了多柔性梁V型集群行为产生的原因和节能的内在机理, 特别是对领头柔性梁的节能机理进行了阐述.Abstract: The phenomenon of aggregation of animals in V formation is ubiquitous in our daily life, such as bird flocks in migration. It is commonly recognized that this collective mode helps to save energy of the group. However, little direct evidence is given. Research on the energy saving mechanism of this collective behavior can not only help to improve the understanding of nature secret, but also lay a foundation for its bionic application. In this paper, a simulation method developed based on Fluent is adopted to solve this fluid-structure interaction problem of hydrodynamic collective behavior of multiple flexible beams in V formation. Specifically, finite volume method is used to simulate the flow field, governing equations of Euler-Bernoulli beam are complemented through user-defined function, and then solved by the mode superposition method and fourth-order Runge-Kutta method. Dynamic mesh technique is adopted to trace the coupling interface between flow field and structural field. The hydrodynamic aggregations of multiple (three or five) self-propelled 2D flexible beams in V configuration are simulated. Three propulsive properties (mean velocity, input power and efficiency) of beams in V formation are compared with the corresponding data of single self-propelled beam. It is found that not only the following beams in V formation possess the promotion of mean velocity and propulsive efficiency, the performance of leading beam also increases, and the growing rate surpasses 14%. Those data provide the direct evidence of energy saving in the collective behavior of V formation. In addition, in order to find out the mechanism of the formation of hydrodynamic aggregation behavior and the reason of the energy saving of beams (especially the leading beam) in V formation, the obtained flow details (vortices contour and pressure contour) are analyzed.
-
引 言
集群是动物较为普遍的行为, 如大多数鸟类和鱼类都会阶段性地集群, 这种聚集行为产生的有利作用称为集群效应. 生物学[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
$ \ll $ 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.2 流固耦合问题控制方程
不可压缩黏性流动力学控制方程包括质量连续性方程和动量方程, 分别为
$$ \nabla \cdot {{\boldsymbol{u}}^f} = 0 ,\quad {\rm{on}}\;{\varOmega ^f} \times (0,T) \tag{1a}$$ $$\begin{split} & {\rho ^f}\frac{{\partial {{\boldsymbol{u}}^f}}}{{\partial t}} + {\rho ^f}{{\boldsymbol{u}}^f} \cdot \nabla {{\boldsymbol{u}}^f} =- \nabla {p^f} +\\ & \qquad\mu {\nabla ^2}{{\boldsymbol{u}}^f} + {\rho ^f}{{\boldsymbol{g}}^f} , \quad {\rm{on}}\; {\varOmega ^f} \times (0,T) \end{split} \tag{1b}$$ 其中, uf是流体速度矢量, ρf是流体密度, pf表示压力, μ是动力黏性系数, gf是外部体积力. 为了求解流体动力学方程, 还需要补充相应的初始条件和边界条件, 分别为
$$ {\left. {{{\boldsymbol{u}}^f}} \right|_{t = 0}} = {\boldsymbol{0}} , \quad {\rm{on}}\; \varOmega _0^f \tag{2a}$$ $$ {{\boldsymbol{u}}^f} = {{\boldsymbol{u}}^i} , \quad {\rm{on}}\; {\varGamma ^i} \times (0,T) \tag{2b}$$ $$ {{\boldsymbol{u}}^f} = {\boldsymbol{0}} , \quad {\rm{on}}\; {\varGamma ^w} \times (0,T) \tag{2c}$$ 其中, 式(2a)表示流体在初始时刻保持静止, 式(2b)表示流固耦合面Γi=Ωf∩Ωs上的速度边界由柔性体的运动ui决定, 式(2c)表示流体域边界Γw为无滑移边界条件. 流体动力学方程采用有限体积法进行求解.
欧拉伯努利梁在x和y方向上的动力学控制方程分别为
$$ {\rho ^l}\frac{{{\partial ^2}X(s,t)}}{{\partial {t^2}}}{\text{ = }}{F_x}(s,t) ,\;\;\;{\rm{on}}\; {\varOmega ^{\text{s}}} \times (0,T) \tag{3a}$$ $$ $$ $$ {\rho ^l}\frac{{{\partial ^2}Y(s,t)}}{{\partial t^2}} + EI\frac{{{\partial ^4}Y(s,t)}}{{\partial {s^4}}} = {F_y}(s,t) , \quad {\rm{on}}\; {\varOmega ^s} \times (0,T) \tag{3b} $$ 其中, 二维柔性梁线密度ρl = ρsh, ρs为柔性体密度, 柔性梁位置矢量X(s, t) = (X(s, t), Y(s, t)), 流体外力矢量F(s, t) = (Fx(s, t), Fy(s, t)), s表示沿着柔性梁轴线的拉格朗日坐标, 柔性梁的杨氏模量为E, 转动惯量为I. 运动过程中认为柔性梁自身长度保持不变, 补充几何约束条件如下
$$ \frac{{\partial{\boldsymbol{ X}}}}{{\partial s}} \cdot \frac{{\partial {\boldsymbol{X}}}}{{\partial s}} = 1 , \quad {\rm{on}}\; {\varOmega ^s} \times (0,T) \tag{3c}$$ 同样地, 为了求解柔性梁的动力学方程还需要补充合适的初始条件和边界条件为
$${\boldsymbol{X}}(s,0) = (X(s,0),Y(s,0)) ,\quad {\rm{on}}\; \varOmega _{\text{0}}^s \tag{4a} $$ $$ {{\partial {\boldsymbol{X}}(s,0)} \mathord{\left/ {\vphantom {{\partial X(s,0)} {\partial t}}} \right. } {\partial t}} = ({{\partial X(s,0)} \mathord{\left/ {\vphantom {{\partial X(s,0)} {\partial t}}} \right. } {\partial t}},{{\partial Y(s,0)} \mathord{\left/ {\vphantom {{\partial Y(s,0)} {\partial t}}} \right. } {\partial t}}) , \quad {\rm{on}}\; \varOmega _{\text{0}}^s \tag{4b} $$ $$ Y(0,t) = y(t){\text{ + }}{y_i}{\text{ (}}i{\text{ = 1,\;2,\;3,\;4,\;5)}} ,\quad {\rm{on}}\; {\varGamma ^l} \times (0,T) \tag{4c}$$ $${\boldsymbol{F}}(s,t){\text{ = }}{p^f}(s,t) \cdot {{\boldsymbol{n}}^f} ,\quad {\rm{on}}\; {\varGamma ^i} \times (0,T) \tag{4d}$$ 其中, 式(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)可以看出, 流场产生的作用于结构物上的流体力将改变结构场的动力学条件, 反过来结构场运动状态(即速度)的改变提供影响流场的边界条件, 流固耦合计算中需要保持耦合面上两场速度和受力的连续性, 表达式分别为
$${{\boldsymbol{u}}^i}{\text{ = }}{{\partial {\boldsymbol{X}}(s,t)} \mathord{\left/ {\vphantom {{\partial X(s,t)} {\partial t}}} \right. } {\partial t}} ,\qquad {\rm{on}}\; {\varGamma ^i} \times (0,T) \tag{5a} $$ $${p^f}(s,t) \cdot {{\boldsymbol{n}}^f}{\text{ = }} - F(s,t) \cdot {{\boldsymbol{n}}^s} ,\qquad {\rm{on}}\; {\varGamma ^i} \times (0,T) \tag{5b}$$ 其中, ns是结构力单位矢量, 指向结构外. 根据如式(5a)耦合面上的运动采用动网格技术更新计算域.
1.3 松耦合方法
采用分区算法中的松耦合方法对多个柔性体−流体介质相互作用流固耦合问题进行求解. 在单个时间步内, 分别对流场和结构场进行一次计算, 然后两场之间数据也只交换一次, 具体计算流程如图2所示, ϕ为表示任意流体参量的通用变量. 从时间步tn到tn + 1的具体计算流程叙述如下:
(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无因次化参量, 则柔性体的线密度无因次量表示为
$ \bar \rho $ = ρl/(ρfL), 雷诺数为Re = ρfUrefL/μ, 弯曲刚度为K = EI/(ρf$U_{{ref}}^2$ ). 选择$ \bar \rho $ = 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为柔性梁指定升沉运动的周期.图4将计算得到的2个柔性梁速度和纵向间距的时历曲线与浸入边界法的结果进行对比可知, 两种方法计算得到的曲线吻合良好. 在t/T = 0至2区间内, 纵向排列2个柔性梁保持相同的速度水平移动, 接着后梁突然加速, 前后梁的纵向间距逐渐变小直到最后达到一个稳定的间距((H−L)/L = 5).
图 4 纵向排列柔性梁导边运动参数时历曲线对比: (a) 水平运动速度; (b) 纵向间距Figure 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型集群行为的模拟中采用相同的设置.
3. 数值模拟结果及分析
对多个柔性体V型排列进行仿真, 用来模拟自然界中动物V型集群行为. 文献[29-30]中指出自然界中动物的弯曲刚度系数近似为O(1), 这里仍然选用第三节验证案例中的柔性体, 参数为
$ \bar \rho $ = 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. 图6展示了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]的区间内作稳定的余弦变化.
图 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)Figure 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个柔性体的集群行为以及相应的涡量和压力云图. 如图7(a)~图7(d)所示, y方向上, 柔性梁导边从其升沉运动的平衡位置(如柔性梁1, Y1(s = 0, t = 0) = 0)向最大负幅值运动(如柔性梁1, Y1
$\Big(s=0,\; t= \dfrac{T}{4}\Big)$ = −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)给出相应典型时刻的流场压力云图, 远场参考压力设置为零, 特征压力pref = ρf$U_{ref}^2 $ 用来无因次化计算得到的压力值, 红色和蓝色分别表示最大正压力值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大.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型集群的动物活体试验结果得出的结论相符合.
图 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)Figure 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选取展示了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.
图11(a)~图11(l)给出了相应时刻的流场压力云图, 根据3.1节对压力场的分析方法, 边界层分离以及涡产生的地方形成低压区, 结合柔性体的形变, 可以对作用在柔性体上流体力进行说明, 这里不再详述. 图12给出了5个柔性体水平方向受力的时历曲线, 第三排柔性体4和5受到的推力(Fx < 0)幅值最大, 第二排柔性体2和3次之, 领头柔性体最小; 领头柔性体受到的阻力(Fx > 0)幅值最小, 后排柔性体较大.
3.3 V型集群性能分析
上述内容对多(3和5)个柔性体的V型集群运动进行了描述, 接下来引入三个主要参数对多柔性体自推进集群运动性能进行讨论[23-25]. 集群运动水平方向平均推进速度Ui, 平均输入功率Pi和推进效率ηi分别为
$$ {U_i} = - \frac{1}{{{T_c}}}\int_t^{t + {T_c}} {\frac{{\partial {X_i}(0,t)}}{{\partial t}}{\rm{d}}t} \tag{6a}$$ $$ {P_i} = \frac{1}{{{T_c}}}\int_t^{t + {T_c}} {\left[ {\int_0^L {{{\boldsymbol{F}}_i}(s,t) \cdot \frac{{\partial {{\boldsymbol{X}}_i}(s,t)}}{{\partial t}}} {\rm{d}}s} \right]{\rm{d}}t} \tag{6b} $$ $$ {\eta _i} = {{\frac{1}{2}{\rho ^l}LU_i^2} \mathord{\left/ {\vphantom {{\frac{1}{2}{\rho ^l}LU_i^2} {({P_i}{T_c})}}} \right. } {({P_i}{T_c})}} \tag{6c}$$ 其中, 下标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 beamsCollective behavior Velocity Input power Efficiency/% U P1 P2 P3 P4 P5 η1 η2 η3 η4 η5 single beam 1.74 1.22 7.90 3 beams 2.01 1.26 1.29 1.28 10.2 10.0 10.1 5 beams 2.02 1.24 1.44 1.44 1 27 1.28 10.5 9.05 9.05 10.3 10.2 接下来尝试对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$U_{ref}^2{L} $ 被用来无因次化计算得到的合力Figure 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$U_{ref}^2{L}$ 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型集群中不同位置的个体采用不同策略提升推进性能. 后排柔性梁节能机理与其他研究者给出的结论一样, 即其可利用前排个体的尾流涡进而节省能耗;本文还对集群中领头柔性梁的节能机理进行了阐述, 即在自身产生的尾涡与后排柔性体的相互作用下, 领头柔性梁自身形变变小和所受推力幅值变大.
为了使多个柔性体的集群运动更加接近自然界中丰富多样的动物集群现象, 研究具有不同长度和弯曲刚度的多个柔性梁采用不同指定运动(升沉、转动或两者相结合)方式自推进进而自发形成的集群运动行为将是下一步的工作方向.
-
图 4 纵向排列柔性梁导边运动参数时历曲线对比: (a) 水平运动速度; (b) 纵向间距
Figure 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]
图 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)
Figure 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
图 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)
Figure 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
图 13 领头柔性梁对比: (a) 水平受力时历曲线; (b) 横向形变时历曲线. 特征值Fref = (1/2)ρf
$U_{ref}^2{L} $ 被用来无因次化计算得到的合力Figure 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
$U_{ref}^2{L}$ 表 1 多个自推进柔性体性能参数比较
Table 1 Propulsive performance of multiple self-propelled flexible beams
Collective behavior Velocity Input power Efficiency/% U P1 P2 P3 P4 P5 η1 η2 η3 η4 η5 single beam 1.74 1.22 7.90 3 beams 2.01 1.26 1.29 1.28 10.2 10.0 10.1 5 beams 2.02 1.24 1.44 1.44 1 27 1.28 10.5 9.05 9.05 10.3 10.2 -
[1] Humston R, Ault JS, Lutcavage M, et al. Schooling and migration of large pelagic fishes relative to environmental cues. Fisheries Oceanography, 2010, 9(2): 136-146
[2] Olson RS, Hintze A, Dyer FC, et al. Predator confusion is sufficient to evolve swarming behavior. Journal of the Royal Society Interface, 2013, 10(85): 20130305 doi: 10.1098/rsif.2013.0305
[3] Peng H, Maldonado-Chaparro AA, Farine DR. The role of habitat configuration in shaping social structure: a gap in studies of animal social complexity. Behavioral Ecology and Sociobiology, 2019, 73(1): 1-14
[4] Shelton DS, Shelton SG, Daniel DK, et al. Collective behavior in wild zebrafish. Zebrafish, 2020, 17(4): 243-252 doi: 10.1089/zeb.2019.1851
[5] Lissaman PBS, Schollenberger CA. Formation flight of birds. Science, 1970, 168(3934): 1035
[6] Weihs D. Hydromechanics of fish schooling. Nature, 1973, 241: 290-291 doi: 10.1038/241290a0
[7] Weimerskirch H, Martin J, Clerquin Y, et al. Energy saving in flight formation. Nature, 2001, 413: 697-698 doi: 10.1038/35099670
[8] Portugal SJ, Hubel TY, Fritz J, et al. Upwash exploitation and downwash avoidance by flap phasing in ibis formation flight. Nature, 2014, 505(7483): 399 doi: 10.1038/nature12939
[9] Liao JC, Beal DN, Lauder GV, et al. Fish exploiting vortices decrease muscle activity. Science, 2003, 302: 1566-1569 doi: 10.1126/science.1088295
[10] Gopalkrishnan R, Triantafyllou MS, Triantafyllou GS, et al. Active vorticity control in a shear flow using a flapping foil. Journal of Fluid Mechanics, 1994, 274: 1-21 doi: 10.1017/S0022112094002016
[11] Beal DN, Hover FS, Triantafyllou MS, et al. Passive propulsion in vortex wakes. Journal of Fluid Mechanics, 2006, 549: 385-402 doi: 10.1017/S0022112005007925
[12] Lauder GV, Anderson EJ, Tangorra J, et al. Fish biorobotics: kinematics and hydrodynamics of self-propulsion. Journal of Experimental Biology, 2007, 210: 2767-2780 doi: 10.1242/jeb.000265
[13] Ristroph L, Zhang J. Anomalous hydrodynamic drafting of interacting flapping flags. Physical Review Letters, 2008, 101(19): 194502 doi: 10.1103/PhysRevLett.101.194502
[14] 王思滢. 柔性体与流体耦合运动的数值模拟和试验研究. [博士论文]. 合肥: 中国科学技术大学, 2010 Wang Siying. Numerical and experimental investigation on the interaction between moving fluid and flexible bodies. [PhD Thesis]. Hefei: University of Science and Technology of China, 2010 (in Chinese)
[15] Boschitshch BM, Dewey PA, Smith AJ. Propulsive performance of unsteady tandem hydrofoils in an in-line configuration. Physics of Fluids, 2014, 26: 051901 doi: 10.1063/1.4872308
[16] Becker AD, Masoud H, Newbolt JW, et al. Hydrodynamic schooling of flapping swimmers. Nature Communications, 2015, 6: 8514 doi: 10.1038/ncomms9514
[17] Ramananarivo S, Fang F, Oza A, et al. Flow interactions lead to orderly formations of flapping wings in forward flight. Physical Review Fluid, 2016, 1: 071201 doi: 10.1103/PhysRevFluids.1.071201
[18] Zhu LD, Peskin CS. Interaction of two flapping filaments in a flowing soap film. Physics of Fluids, 2003, 15(7): 1954-1960 doi: 10.1063/1.1582476
[19] Huang WX, Shin SJ, Sung HJ. Simulation of flexible filaments in a uniform flow by the immersed boundary method. Journal of Computational Physics, 2007, 226: 2206-2228 doi: 10.1016/j.jcp.2007.07.002
[20] Zhang L, Zou SY, Wang CZ, et al. A loosely-coupled scheme for the flow-induced flapping problem of two-dimensional flexible plate with strong added-mass effect. Ocean Engineering, 2020, 217: 107656 doi: 10.1016/j.oceaneng.2020.107656
[21] Zhu LD. Interaction of two tandem deformable bodies in a viscous incompressible flow. Journal of Fluid Mechanics, 2009, 635: 455-475 doi: 10.1017/S0022112009007903
[22] Uddin E, Huang WX, Sung HJ. Actively flapping tandem flexible flags in a viscous flow. Journal of Fluid Mechanics, 2015, 780: 120-142 doi: 10.1017/jfm.2015.460
[23] Zhu XJ, He GW, Zhang X. Flow-mediated interactions between two self-propelled flapping filaments in tandem configuration. Physical Review Letters, 2014, 113: 238105 doi: 10.1103/PhysRevLett.113.238105
[24] Peng ZR, Huang HB, Lu XY. Hydrodynamic schooling of multiple self-propelled flapping plates. Journal of Fluid Mechanics, 2018, 853: 587-600 doi: 10.1017/jfm.2018.634
[25] Peng ZR, Huang HB, Lu XY. Collective locomotion of two closely spaced self-propelled flapping plates. Journal of Fluid Mechanics, 2018, 849: 1068-1095 doi: 10.1017/jfm.2018.447
[26] 彭泽瑞. 仿生自主推进柔性板集群运动的流固耦合数值研究. [博士论文]. 合肥: 中国科学技术大学, 2018 Peng Zerui. Numerical investigation of hydrodynamic schooling of bio-inspired self-propulsive flexible plates. [PhD Thesis]. Hefei: University of Science and Technology of China, 2018 (in Chinese)
[27] Zhang L, Wang CZ, Sun JL, et al. Study on the hydrodynamic aggregation of parallel self-propelled flexible plates based on a loosely coupled partitioned algorithm. Ocean Engineering, 2021, 223: 108703 doi: 10.1016/j.oceaneng.2021.108703
[28] 克拉夫, 彭津. 结构动力学. 王光远译. 北京: 高等教育出版社, 2006 Clough R, Penzein J. Dynamics of Structures. Wang GY trans. Beijing: Higher Education Press, 2006 (in Chinese)
[29] Hua RN, Zhu L, Lu XY. Locomotion of a flapping flexible plate. Physics of Fluids, 2013, 25(12): 121901 doi: 10.1063/1.4832857
[30] Liu K, Huang H, Lu XY. Hydrodynamic benefits of intermittent locomotion of a self-propelled flapping plate. Physical Review E, 2020, 102: 053106
[31] Liu K, Huang H, Lu XY. Self-propelled plate in wakes behind tandem cylinders. Physical Review E, 2019, 100: 033114 doi: 10.1103/PhysRevE.100.033114
-
期刊类型引用(5)
1. 张磊,封少雄,谭昆,郭涛,宋成果,初秀民,苗洋. 基于船舶-流场运动耦合的内河航道设计方法. 上海交通大学学报. 2025(01): 89-98 . 百度学术
2. 王皓田,张磊,张智勇,谭昆,吕品,吴卫国. 基于船舶-流场运动耦合的内河典型碍航结构物绕航机理. 船舶工程. 2025(02): 44-52 . 百度学术
3. 谭昆,郭涛,宋成果,苗洋,初秀民,吴卫国,张磊. 基于船舶–流场耦合作用的内河航段船舶航行虚拟辅助驾引. 中国舰船研究. 2024(02): 71-80 . 百度学术
4. 张磊 ,陈华芳 ,刘斌 ,段健 ,裴志勇 . 均匀流下两柔性梁纵列行为变化研究. 力学学报. 2024(11): 3165-3177 . 本站查看
5. 刘小峰,陈果,刘宇,王希. 动物集群行为的机制和应用. 科学通报. 2023(23): 3063-3076 . 百度学术
其他类型引用(0)