0 引言

自然界中飞行及游动生物具有高效率、低噪音和高机动性等卓越运动能力，远远超过人类现有飞行器与水下航行器的性能。这激发着人类对自然界生物及其运动能力的不断探索，以期开发出更先进的飞行器与水下航行器。数值模拟具有可重复性强、测力简单、可以获得更多流场信息等优点，弥补了实验研究与理论分析的不足，逐渐成为当前仿生运动研究的热点。基于贴体网格求解NS方程首先在仿生流动领域得到应用^[1]。然而生物运动幅度大，贴体运动网格重构算法鲁棒性和运动网格质量难以保证，很容易导致计算效率降低、模拟精度下降甚至计算发散。而浸入边界法(Immersed boundary method，IBM)^[2]与格子玻尔兹曼方法^[3]均采用笛卡尔网格(Lattice Boltzmann method，LBM)，因此可以将格子玻尔兹曼方法与浸入边界法结合起来，使其具有格子玻尔兹曼方法的超高并行特性又具有浸入边界法灵活处理边界的能力。在此思想上发展起来的IB-LBM方法对于大变形边界与复杂形状几何边界具有良好的稳定性，因此被广泛用于仿生运动数值模拟^[4-7]。

目前大量仿生运动研究针对大展弦比模型，将研究焦点放在沉浮、俯仰拍动翼型的二维效应上^[8]。数值结果对于理解自然界中的信天翁(展弦比18)、盲蜘蛛(展弦比11)等大展弦比的生物具有重要意义。然而自然界中生物尺度千差万别，例如隆头鱼科展弦比介于1.5与3.5之间^[9]，银鮫展弦比为2.2^[10]。大量实验与计算表明三维仿生翼诱导产生的涡系结构远比二维翼型产生的涡系结构复杂。三维拍动翼尾缘中截面脱落的主涡，在向下游扩散的过程中，与翼型两端脱落的翼尖涡相互连接形成涡环结构。这种涡环结构在大量的研究中得到证明^[11-13]。因此对这些小展弦比生物进行仿生研究需考虑不同展弦比带来的三维效应，仅仅采用二维模型来进行模拟已经不足以反映其流动机理。

三维动边界非定常流动数值模拟耗时严重，目前对三维仿生运动翼的数值研究少于二维模型，而关于几何参数和流动参数对拍动翼三维流动效应影响的数值研究更为少见。Dong等人采用NS方程数值研究了椭球翼型在俯仰与沉浮运动下，不同展弦比与雷诺数等参数对流场涡系结构的影响^[14]。陈刚等人研究了拍动翼不同几何形状对推进性能与涡系结构的影响^[15]。本文在此基础上进一步研究不同展弦比对拍动翼推进性能与涡系结构的影响。针对三维复杂外形仿生运动翼高精度数值模拟耗时巨大的问题，首先发展一套面向超级计算机的三维IB-LBM并行求解器。然后在广州天河2号超级计算机上，对不同展弦比的NACA0012仿生拍动翼进行大规模数值模拟。通过分析不同展弦比仿生运动翼的非定常涡系结构及其演化规律对推进性能的影响，为仿生飞行器与潜水器时仿生运动翼参数设计提供借鉴。

1 数值方法与验证 1.1 IB-LBM方法

描述流体粒子运动带外力项的格子玻尔兹曼方程如下：

$ \begin{array}{l} {f_i}\left( {x + {e_i}\delta t, t + \delta t} \right) - {f_i}\left( {x, t} \right) = \\ - \frac{1}{\tau }\left[{{f_i}\left( {x, t} \right)-f_i^{{\rm{eq}}}\left( {x, t} \right)} \right] + \frac{2}{3}{\omega _i}f\cdot{e_i}\delta t \end{array} $

(1)

式中：f_i为密度分布函数，f_i^eq为对应的平衡态分布函数，τ为无量纲松弛时间，δt为时间步长，f为添加体力密度，e_i是离散速度模型，ω_i是权系数，其中e_i与ω_i与离散速度模型有关。方程(1)通过Chapman-Enskog展开可还原为到带体积力项的不可压NS方程。LB方法采用笛卡尔直角网格离散流体区域，流体被看成流体微团，只能向其他网格点运动或者静止不动。本文采用三维D3Q19离散速度格子模型，其平衡态分布函数f_i^eq为：

$ f_i^{{\rm{eq}}} = \rho {\omega _i}[1 + 3{e_i}\cdot u + \frac{9}{2}{\left( {{e_i}\cdot u} \right)^2}-\frac{3}{2}u\cdot u] $

(2)

式(2)中对应的权系数为：

$ {{\omega _i} = \left\{ \begin{array}{l} 1/3\;\;\;\;\;\;\;\;\;i = 0\\ 1/18\;\;\;\;\;\;\;i = 1, 2, \cdots, 6\\ 1/36\;\;\;\;\;\;\;i = 7, 8, \cdots, 18 \end{array} \right.} $

(3)

对应本文采用的D3Q19模型，离散速度模型为：

$ {e_i} = \left\{ \begin{array}{l} \left( {0, 0, 0} \right)\;\;\;\;\;\;i = 0\\ \left( { \pm 1, 0, 0} \right), \left( {0, \pm 1, 0} \right), \left( {0, 0, \pm 1} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i = 1, 2, \cdots, 6\\ \left( { \pm 1, \pm 1, 0} \right), \left( {0, \pm 1, \pm 1} \right), \left( { \pm 1, 0, \pm 1} \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i = 7, 8, \cdots, 18 \end{array} \right. $

(4)

求解式(1)后，流体网格点上的速度压力密度等通过下式求得:

$ u = \frac{{{\sum _i}{e_i}{f_i} + 0.5f\delta t}}{\rho } $

(5)

$ \rho = \sum\limits_i {{f_i}} $

(6)

$ p = c_s^2\rho $

(7)

将浸入边界法与格子玻尔兹曼方法耦合的关键在于如何处理LB方程中与N-S方程对应的体力密度分布项f。浸入边界格子玻尔兹曼方法的基本思想是将浸入边界上集中力分散到周围流体欧拉网格点上，然后再求解带外力项的格子玻尔兹曼方程(1)。为了方便引入IB-LBM的概念，假设流体区域Ω用欧拉坐标x表示，浸入边界Γ用拉格朗日坐标s表示，将拉格朗日节点上构造的的集中作用力F(s, t)分散到作用范围内的欧拉网格点x上，即可得到欧拉网格点上的体力项f(x, t)：

$ f\left( {x, t} \right) = \int_\mathit{\Gamma } {F\left( {s, t} \right)\delta \left( {x - X\left( {s, t} \right)} \right){\rm{d}}\mathit{s}} $

(8)

式(6)中三维δ(r)函数如下式所示：

$ \delta \left( r \right) = \left\{ \begin{array}{l} \frac{1}{4}\left( {1 + {\rm{cos}}\left( {\frac{{\pi \left| r \right|}}{{2\Delta x}}} \right)} \right), \;\;\;\;\;\left| r \right| \le 2\Delta x;\\ 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left| r \right| \ge 2\Delta x; \end{array} \right. $

(9)

浸入边界法上某点的集中力即式(6)中可以由不同的方法来构造，可以分为虚拟边界模型^[16]，虚拟弹簧力模型^[17]和直接力模型^[18]。上述三个模型均可以模拟物面边界条件，但均有其适用范围。虚拟边界模型适合处理刚性物体及其主动运动，弹簧力模型和直接力模型适合处理流固耦合运动。本文研究的对象为刚性主动拍动翼数值模拟，因此选用虚拟边界模型。虚拟边界模型中集中力通过下式计算：

$ \begin{array}{l} F\left( {s, t} \right) = \alpha \int_0^t {\left[{u\left( {X\left( {s, t} \right), {t_1}} \right)-u} \right]{\rm{d}}{\mathit{t}_1} + } \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\beta \left[{u\left( {X\left( {s, t} \right), t} \right] -u} \right] \end{array} $

(10)

其中u为物体表面速度。对于指定运动物体其速度u(X(s, t), t)可事先通过附近欧拉点上的速度信息插值得到。α与β为反馈力系数分别对应速度时间积分与速度反馈。关于速度修正IB-LBM算法的详细介绍参考相关文献^[19]，这里不再赘述。

本文针对刚性仿生运动翼所建立的三维IB-LBM并行求解器算法流程如下：

① 设定初始参数求解带外力项的格子玻尔兹曼方程(1)；初始化流场，调用几何模型函数并初始化三维仿生翼位置。

② 求解运动翼的运动控制方程计算下一时间步模型所处位置；从公式(10)计算浸入边界上集中力，再由公式(8)将集中力分散到对应欧拉网格点。

③ 求解带外力项的格子玻尔兹曼方程(1)，获得欧拉网格点上速度压力等信息，并输出流场信息与结构位置。

④ 若未达到终止条件，回到步骤(2)进行下一时间步迭代；否则结束迭代，进行后处理。

1.2 求解器验证

Dong基于有限体积法结合浸入边界法求解不可压NS方程，对椭圆翼模型进行了系统的数值研究^[18]。椭球拍动翼模型如图 1所示，椭球翼在来流方向长度为1，展向方向的距离为2，最大厚度为0.12^[18]。其中X-Axis方向与来流方向重合，Y-Axis方向为椭球翼展向方向，Z-Axis为椭球翼厚度方向。本节通过与DONG的数值模拟结果对比来验证本本文三维IB-LBM并行求解器的有效性。

对自然界飞行和游动生物仿生运动进行数值模拟时，一般将其运动抽象为沉浮与扑动组合运动模型。根据对比文献，仿生翼中心在z方向按照式(9)沉浮运动，同时拍动翼绕其中心按照式(10)转动：

$ z\left( t \right) = {A_z}{\rm{sin}}(2\pi \mathit{f}t) $

(11)

$ \theta \left( t \right) = {A_\theta }{\rm{cos}}(2\pi \mathit{f}t) $

(12)

其中A_z为0.5，转动幅值A_θ为π/6，f为椭球翼的拍动频率。

仿生翼推力系数与升力系数定义为：

$ {C_T} = \frac{T}{{0.5\rho U_\infty ^2A}}, \ \ \ {C_L} = \frac{L}{{0.5\rho U_\infty ^2A}} $

(13)

式中：T为运动翼顺气流方向受力，L为运动翼受到法向升力，A为拍动翼面积。

经过网格收敛性检验后本算例所采用的计算网格数量约3300万。图 2给出了本文预测的推力系数和升力系数随时间变化与文献结果对比。二者吻合良好，为了验证本文求解器对流场信息的捕捉能力，进一步给出了椭球扑翼展向对称面上的涡量云图。如图 3所示，其中图 3(a)为本文结果，图 3(b)为文献结果。从图中可以看到两种方法捕捉到了十分相似的涡量云图，进而验证了本文求解器的可信度。关于本文求解器对平板、圆球、椭球翼和波动翼等复杂模型升推系数和涡系结构精细化捕捉能力的详细验证请参见相关文献^{[7, 15, 20]}，本文不再赘述。

2 结果与讨论 2.1 拍动翼数学模型

本文主要研究不同展弦比对NACA0012拍动翼的流场涡结构与动力学系数的影响。NACA0012仿生拍动翼的几何模型如图 4所示。仿生翼在X方向长度为1，Y方向长度分别为1，2，4，6，即拍动翼的展弦比分别1，2，4，6。拍动翼采用经典拍动推进运动模式，仿生翼中心在Z方向按照式(9)做沉浮运动，同时仿生翼围绕其中心按照式(10)做俯仰运动。为研究低雷诺数下仿生拍动翼的推动性能，并使其运动尽可能的接近生物的拍动特征，故将详细控制参数设置为：浮沉幅值A_z为0.5；俯仰幅值A_θ为π/6；拍动频率f为0.6；拍动翼雷诺数Re为200，其中雷诺数定义为Re=U_∞L/υ。Strouhal数St为0.6，其中Strouhal数定义为St=2A_zf/U_∞。流体的来流速度U_∞为1。当拍动翼处于运动轨迹中间位置时，其中心点坐标为(5.0，5.0，5.0)。流体计算域的尺寸为18×10×10。翼型在X方向30个单位长度，即dx=0.033，流体计算域的整体网格量为540×300×300，总计4860万，时间步长dt=0.001。采用广州超算中心天河2超级计算机运行本文发展的并行求解器。计算步数为2万步。每个任务采用192个CPU计算核心，一个算例大约需要4个小时。

2.2 平均速度场

四个不同展弦比仿生拍动翼顺气流方向的速度平均场如图 5所示。平均速度值可通过一个周期内的瞬时速度取算术平均得到。图 5(a)(b)与(c)和(d)所示的平均速度分布明显不同；图 5(a)所示展弦比为1的NACA0012拍动翼与图 5(b)展弦比为2的NACA0012拍动翼平均速度场分布形状非常相似，均呈分叉的射流，但是展弦比为2的拍动翼顺气流方向平均速度场强度明显高于展弦比为1的拍动翼。根据动量定理，预示着展弦比为2的拍动翼推力系数高于展弦比为1的NACA0012拍动翼。当展弦比为4时，尾流区的分叉射流结构消失，取而代之的是一条主射流，但是该射流的两侧的速度强度高于中心线附近速度强度。当展弦比为6时，尾流区仍然呈现为一条主射流，但是射流的角度小于展弦比为4时的射流的角度。四幅图具有相同的颜色空间，展弦比为6的拍动翼射流的角度比较小，但是其分布强度偏小，因此需要进一步分析动力学系数。

2.3 升推系数

为了方便后续讨论，将展弦比定义为AR。图 6给出了不同展弦比AR下拍动翼的动力学系数与频率和时间组合参数ft的关系曲线。该曲线取自翼型拍动的第11、12周期，此时流场已经达到稳定状态。由于NACA0012翼型关于xy平面运动对称，所以推力系数变化频率为升力系数变化频率的两倍。尽管升力系数随着组合参数ft在变化，其最大幅值达到在4左右，但四个不同展弦比运动翼的平均升力系数均为0。这与零迎角下对称翼型升力为零的结论吻合，也再次验证了本文求解器模拟结果的可用性。

而从图 6(a)所示的不同展弦比运动翼的推力系数变化情况来看，其随着展弦比变化的波动差异要远大于图 6(b)所示的升力系数曲线的波动差异情况。特别是在推力系数波峰和波谷，不同展弦比运动翼的推力曲线差异尤其明显。这表明在本文所选取状态下，展弦比对仿生拍动翼推力系数具有更显著的影响。AR=1，2，4，6的仿生翼的平均推力系数分别为0.339，0.453，0.456，0.304。图 7给出了平均推力系数随展弦比变化曲线。当AR从1到2变化时，平均推力系数增加，这与展弦比为2的仿生翼顺气流方向平均速度场强度较高一致。AR从2到4变化时，推力系数变化不明显。当AR从4变化到6时，推力系数下降。从对应AR=6的平均速度场图 5(d)来看，虽然翼型后部射流角度变小，但是流体速度同时也变小，从而导致推力系数下降。

2.4 涡系结构

为了更进一步揭示不同展弦比对仿生拍动翼推进性能的影响机理，本节对涡系结构进行分析。图 8所示为不同展弦比拍动翼展向对称面上的涡量图，四幅涡量图具有相同的颜色空间。对比涡量图发现，展弦比为1时拍动翼上侧的前沿涡较小，并且距离拍动翼较远；当展弦比增加时，拍动翼上侧的前沿涡大小形状基本不变与展弦比无关。对于拍动翼的尾涡，当展弦比为1时，尾涡比较细长并且距离拍动翼较远；当展弦比为2时，尾涡与展弦比为1时相比较粗且距离拍动翼较近，诱导出来的速度场在拍动翼尾部较强，但是仍呈现射流状，如图 5(a)与图 5(b)顺气流方向平均速度场所示。这是由于虽然展弦比增加，拍动翼上下表面诱导的涡强度随着展弦比的增加而增加，但是拍动翼两侧的翼尖涡的强度仍相对较大，呈现出相似的流场特性。当展弦比继续增加，拍动翼上下表面产生的涡的强度增加，拍动翼两端产生的翼尖涡相对强度降低，翼尖涡对主涡的作用进一步被削弱，从而涡量图呈现出与二维翼型相似的特性。甚至在图 8(d)中，翼型表面产生的涡与上周期产生的涡呈现出典型的反卡门涡街，涡诱导产生速度场夹角较小。

因为不同展弦比拍动翼推力系数在波峰和波谷位置差异最大，图 9给出了拍动翼处于波谷也即最低位置时由Q准则识别的涡系结构。Q准则被定义为$Q = \frac{1}{2}\left( {{\mathit{\Omega }_{\mathit{ij}}}{\mathit{\Omega }_{\mathit{ij}}}\mathit{ - }{\mathit{S}_{\mathit{ij}}}{\mathit{S}_{\mathit{ij}}}} \right)$，其中Ω_ij=(u_{i, j}-u_{j, i})/2与S_ij=(u_{i, j}+u_{j, i})/2分别为速度梯度∇u的对称与反对称张量。从图中可以发现拍动翼上下表面脱落涡的展向距离与拍动翼的展向长度相近；即随着展弦比的增加，涡环的展向长度增加。但是不同展弦比的拍动翼诱导的脱落涡在顺气流方向尺度基本一致。涡环在向下游演化的过程中，展向距离均缩小。但是四者仍有不同之处，当展弦比为1，2时，拍动翼上下表面脱落的主涡比较窄，即涡的强度相对较低，在翼尖涡的影响下，涡系结构呈现出涡环结构。当展弦比4，6时，拍动翼上下表面脱落的主涡展向距离较长，翼尖涡的强度相对于主涡比较弱，翼尖涡没有足够的能量将涡环封闭，从而呈现为互相交错的涡系结构。翼尖涡诱导速度分量产生，从而使拍动翼上下表面产生的涡互相远离，与图 9所示的展向对称面上大展弦比拍动翼尾涡的角度较小相一致。

3 结论

不同几何、来流和运动参数对仿生运动翼流动机理与升推性能的影响机理是高性能仿生飞行器和潜航器气动外形与运动控制设计的关键科学问题。三维复杂外形仿生运动翼高精度数值模拟存在耗时巨大的问题，而进行参数影响机理的精细化研究则需要海量大规模计算。本文采用基于速度修正的IB-LBM算法发展了三维IB-LBM仿生运动翼大规模并行求解器。在广州天河2号超级计算机上采用5000万量级网格，对低雷诺数下不同展弦比的NACA0012仿生拍动翼进行了大规模并行数值模拟，捕捉到了仿生拍动翼的精细化涡系结构并给出了其升推系数变化规律。

数值结果表明：随着展弦比的增加，NACA0012拍动翼的平均推力系数先增加后减小；尾迹区中的涡环结构的展长增加，而涡环顺气流方向长度基本不变。特别是随着展弦比的增加，展向对称面的涡量分布二维特性越明显，这表明小展弦比仿生运动翼其二维特性减弱，必须考虑其三维效应采用三维方法对其升推性能进行研究。在进行仿生航行器气动设计时，有必要细致考虑不同展弦比对仿生拍动翼推进性能与涡系结构的影响。后续将进一步研究拍幅度和排频等运动方式、Strouhal数和雷诺数等来流条件以及编队构型等对仿生拍动翼推进性能与涡系结构演化机制的影响。