引言

柔性拍翼推进策略广泛存在于昆虫和鸟类飞行, 这类飞行方式的优越气动表现和良好机动特性使其受到了广泛的关注^[1-9].对拍翼飞行时的流固耦合问题, 已经开展了大量的研究.例如, 研究表明, 昆虫飞行时其柔性翅膀将产生极大的变形, 这种大变形对于昆虫的良好气动表现具有十分重要的意义^[2-5].然而, 值得注意的是, 拍翼在取得优异的气动特性的同时产生的噪音也并不十分显著.因此, 对拍翼飞行的流固声耦合问题进行研究, 揭示其飞行时的声学特性, 对于应用这类飞行技术具有重要的指导意义^[10-16].

Sueur等采用实验对苍蝇飞行时的声学问题进行了研究, 指出其拍翼发声在前后两侧分别为拍翼拍动频率和两倍的拍动频率所主导^[13]. Bae等的二维数值模拟研究表明, 拍翼拍动频率和2倍的该频率分别主导水平和数值方向的声场信号^[11]. Inada等指出拍翼附近的声场与其表面的压力密切相关^[14]. Geng等运用声学分离数值方法对蝉在三维空间中前飞时的发声问题进行了数值模拟研究, 并指出刚性和柔性拍翼发声的方向不同, 与拍翼的动力学特性高度相关, 且拍翼的转动可以降低噪声^[10].值得指出的是, 此前大量的文献[1-9]针对柔性拍翼的气动特性进行了研究, 但鲜有对拍翼发声的流固声耦合问题的研究, 文献[10]中对拍翼发声的研究也是采用预给定运动的拍翼, 并非双向的流固声耦合.本文作者曾首次采用直接数值模拟方法对前飞的拍翼发声的流固声问题进行过较为系统的研究, 指出柔性拍翼在选择适当的质量比和刚度时可以产生一定的降噪效果^[12].在此基础上, 本文采用这种浸入边界法对悬停飞行时的弹性拍翼发声问题进行进一步研究.具体针对不同的拍翼-流体质量比和拍翼刚度对拍翼的气动和发声特性进行数值模拟研究, 并对结果进行讨论和分析.本文内容组织如下:第1节对物理问题进行描述; 第2节介绍数值模拟方法; 第3节对结果进行分析和讨论; 最后给出结论.

1 物理问题描述

本文在二维空间内对悬停拍翼的流固声耦合问题进行数值模拟研究. 图 1是该问题的示意图.如图所示, 拍翼的前缘顶点在平动的同时并转动, 具体由如下两式进行描述^{[4, 9]}

$ \begin{array}{*{20}{c}} {{X_0}(t) = 0.5{A_0}{\rm{cos}}(2\pi {f_0}t)[{\rm{cos}}\beta ,{\rm{sin}}\beta ]}\\ {\alpha (t) = 0.5{\alpha _{\rm{m}}}{\rm{sin}}(2\pi {f_0}t + \phi )} \end{array} $

其中, β为拍动面的倾角, f₀为拍动频率, A₀和α_m分别为平动和转动的幅值.本文计算中取A₀=1.25, β=0, 相位差ϕ=0.

该问题的无量纲控制参数包括Reynolds数、Mach数、质量比和拍动频率比, 具体定义如下

$ \begin{array}{*{20}{l}} {Re = {\rho _{\rm{f}}}UL/\mu ,\quad Ma = \pi U/c}\\ {m = {\rho _{\rm{s}}}/{\rho _{\rm{f}}}L,\omega = 2\pi {f_0}/{\omega _n}} \end{array} $

其中, L为拍翼的弦长, ρ_s和ρ_f分别为拍翼的线密度和流体的面密度, U=2f₀A₀为拍翼前沿顶点的平均平动速度, c为未扰动流体中的声速, ω_n为拍翼的1阶本征频率, $\omega_{n}=k_{n}^{2} / L^{2} \sqrt{K_{\mathrm{B}} / \rho_{\mathrm{s}}}$, 其中K_B为拍翼的抗弯刚度, k_n=1.875 1^[12].本文计算中, Re=100, Ma=0.1.升力系数, 动力系数和效率定义如下

$ {C_{\rm{L}}} = \frac{{2{F_{\rm{L}}}}}{{{\rho _{\rm{f}}}{U^2}L}},\quad {C_{\rm{P}}} = \frac{{ - 2\int_0^L \mathit{\boldsymbol{f}} \cdot \mathit{\boldsymbol{v}}{\rm{d}}l}}{{{\rho _{\rm{f}}}{U^3}L}},\eta = \frac{{{C_{{\rm{L,m}}}}}}{{{C_{{\rm{P,m}}}}}} $

其中, F_L为流体作用在拍翼上的升力(竖直方向), f为作用在拍翼各节点上的力, v为拍翼上各节点的速度, C_{L, m}和C_{P, m}分别为平均升力和动力系数.声压定义如下

$ \Delta \tilde p = p - \bar p $

其中, p为瞬时压力, p为时间平均压力.另外, 瞬时波动压力Δp=p-p₀, p₀为未扰动流场中的初始压力^{[17, 26]}.

2 数值方法

为了模拟弹性拍翼发声的流固声耦合问题, 本文采用一种浸入边界法在二维计算域进行直接数值模拟研究.该方法包括3个部分:流体求解器、固体求解器和流固耦合算法.

本文考虑可压缩流体, 二维空间内的控制方程为^{[12, 18]}

$ \frac{{\partial \mathit{\boldsymbol{Q}}}}{{\partial t}} + \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial x}} + \frac{{\partial \mathit{\boldsymbol{G}}}}{{\partial y}} - \frac{1}{{Re}}\left( {\frac{{\partial {\mathit{\boldsymbol{F}}_u}}}{{\partial x}} + \frac{{\partial {\mathit{\boldsymbol{G}}_v}}}{{\partial y}}} \right) = \mathit{\boldsymbol{S}} $

其中

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}} = {{[{\rho _{\rm{f}}},{\rho _{\rm{f}}}u,{\rho _{\rm{f}}}v,E]}^{\rm{T}}}}\\ {\mathit{\boldsymbol{F}} = {{[{\rho _{\rm{f}}}u,{\rho _{\rm{f}}}{u^2} + p,{\rho _{\rm{f}}}uv,(E + P)u]}^{\rm{T}}}}\\ {\mathit{\boldsymbol{G}} = {{[{\rho _{\rm{f}}}v,{\rho _{\rm{f}}}uv,{\rho _{\rm{f}}}{v^2} + p,(E + P)v]}^{\rm{T}}}}\\ {{\mathit{\boldsymbol{F}}_u} = {{[0,{\tau _{xx}},{\tau _{xy}},{b_x}]}^{\rm{T}}}}\\ {{\mathit{\boldsymbol{G}}_v} = {{[0,{\tau _{xy}},{\tau _{yy}},{b_y}]}^{\rm{T}}},}\\ {{b_x} = u{\tau _{xx}} + v{\tau _{xy}},{b_y} = u{\tau _{xy}} + v{\tau _{yy}}} \end{array} $

τ_ij为剪切应力, S为源项.为了有效捕捉微弱的声压信号, 流体求解器采用5阶精度的加权本质无振荡(weighted essentially non-oscillatory，WENO)算法, 具体算法实现可参考文献[18-19].

二维空间内弹性拍翼动力学的控制方程有如下形式^{[12, 16-18]}

$ {\rho _{\rm{s}}}\frac{{\partial {\mathit{\boldsymbol{X}}^2}}}{{\partial {t^2}}} + \frac{\partial }{{\partial s}}\left[ {{K_{\rm{S}}}\left( {\left| {\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial s}}} \right| - 1} \right)\frac{{\partial \mathit{\boldsymbol{X}}}}{{\partial s}}} \right] + {K_{\rm{B}}}\frac{{\partial {\mathit{\boldsymbol{X}}^4}}}{{\partial {s^4}}} = {\mathit{\boldsymbol{F}}_{\rm{f}}} $

其中, s为弧长坐标, X为固体节点坐标, K_S为拍翼的抗拉刚度.本文采用基于绝对节点坐标法的有限元法进行求解, 具体可参考文献[16, 18, 20].

在本文采用的数值方法中, 流体和固体分别独立求解, 然后通过一种罚浸入边界法(immersed boundary method)来实现流固界面的无滑移边界条件.由于其较高的计算效率, 浸入边界法在流固耦合中已经得到广泛应用^[21-28].在罚浸入边界法中, 流体作用在固体上的力可通过下式进行显示计算^{[12, 18, 29]}

$ {\mathit{\boldsymbol{F}}_{\rm{f}}} = \alpha \int_0^t {({\mathit{\boldsymbol{U}}_{{\rm{ib}}}} - \mathit{\boldsymbol{U}})} {\rm{d}}t + \beta ({\mathit{\boldsymbol{U}}_{{\rm{ib}}}} - \mathit{\boldsymbol{U}}) $

其中, U_ib为由流体节点插值得到的固体节点上的虚拟速度, U为固体节点的真实速度. α和β为罚系数.固体作用在流体上的力为-F_f.虚拟速度和流体节点上的分布力可由以下两式进行计算

$ \begin{array}{l} {\mathit{\boldsymbol{U}}_{{\rm{ib}}}}(s,t) = \int_V \mathit{\boldsymbol{u}} (x,t){\delta _h}(\mathit{\boldsymbol{X}}(s,t) - x){\rm{d}}x\\ \mathit{\boldsymbol{f}}(x,t) = - \int_\varGamma {{F_{\rm{f}}}} (s,t){\delta _h}(\mathit{\boldsymbol{X}}(s,t) - x){\rm{d}}s \end{array} $

其中, u和x分别为流体节点的速度和坐标, δ_h为Dirac函数

$ {\delta _h}(x,y) = \frac{1}{{{h^2}}}\lambda \left( {\frac{x}{h}} \right)\lambda \left( {\frac{y}{h}} \right) $

其中, λ为插值函数, 本文采用Peskin等^[24]提出的4点插值函数, 其具体形式见文献[18, 30].

为提高计算效率, 流场采用拉伸网格, 如图 2所示.计算尺寸为125L×125L, 包含1个均匀网格区(总尺寸为4L×3L, 网格尺寸为L/40)用于流固耦合计算.

该算法已经过大量的验证, 包括圆柱扰流问题, 丝线流固耦合问题, 激波作用下的弹性板变形问题以及微型无人飞行器的流固声耦合问题.具体内容可参考文献[12, 18, 26], 这里不再对该算法进行验证.

3 结果及讨论 3.1 刚性拍翼的气动和声学特性

首先, 考虑刚性拍翼在悬停时的发声问题.计算中, 保持平动振幅A₀=1.25, 改变转动的幅值来研究其对拍翼气动和发声特性的影响.共考虑α_m=π/6, π/4, π /3和π/2这4个转动幅值.

图 3所示的是升力系数和动力系数在两个周期内的时间历程曲线.可以看出, 当α_m取π/6, π/4, π/3和π/2时, C_L和C_P具有很好的相似性.随着α_m的增加, 升力和动力系数峰值逐渐增加, 尤其是当α_m= π/2时, 升力和动力系数的峰值显著增加, 但同时在t=2T/4~3T/4时也产生了明显的反向升力, 因此采用平均升力能够更直接地反映出拍翼的气动表现, 后面将做进一步讨论.

图 4是在距离原点40L处测得的波动压力平方根的极坐标图.结果显示, 总体来看在α_m=π/6时声压最小.其中, 当α_m从π/6增加到π/3时, 波动压力的变化不大; 当α_m进一步增加到π/2时, 波动压力显著增大.另外, 图 4表明声压波动在左(90°-180°-270°)略小于右(90°-0°-270°)半侧, 且随着α_m的增大, 可以观测到波动压力逆时针方向略有转动.这些现象也可以直接从图 5所示的瞬时波动压力云图看出.值得注意的是这一现象在文献中研究前飞拍翼时并未观察到, 该现象很可能是由于拍翼在悬停时, 局部的涡不会如前飞时那样在来流作用下迅速向后消散, 如图 6所示.这些局部的涡结构具有很强的随机分布特性, 因此而引起的涡声也使得近场的声场受到了明显的影响.另外, 图 6中的涡量(采用拍翼前缘平均速度和翼长无量纲化)显示, 当转动角度增大时, 局部的涡量也明显增加, 这也很好地解释了在大转动角度时近场声场分布的不对称性更加显著.值得注意的是昆虫在飞行时也会利用其中的涡结构, 如前缘涡, 来加强其飞行的动力学表现.

为了分析拍翼拍动产生的声压的频率特性, 进一步采用快速Fourier变换对压力周期性波动后的10个周期数据进行分析, 得到波动压力不同频率的极坐标图, 如图 7所示.可以看出, 当α_m=π/6, π/4, π/3时, 声场由频率f₀主导, 当α_m增加到π/2后, 高阶频率(如2f₀和3f₀)的声压也明显增大.

表 1给出了拍翼在不同转动幅值时的平均升力系数, 动力系数和效率.结果显示, 随着α_m的增大, 平均升力系数逐渐增大, 平均动力系数先减小后增大, 因此效率先增大后减小, 在α_m= π/6时取得最大值, 这与文献[1]中前飞拍翼的结果一致.综合图 4中的声压结果可以看出, 当α_m= π/3时拍翼效率最高且产生的噪音最小.

3.2 柔性对拍翼声学特性的影响

昆虫飞行时, 柔性拍翼会产生较大的变形, 大量文献研究表明这种大变形有助于产生更大的升力和推力, 并提高拍翼的效率^[2-4].调研发现, 前人对拍翼的气动特性表现研究较多, 但是对其发声问题研究很少^[11-12].本节在前人研究的基础上, 同时对柔性拍翼的气动特性和声学特性进行了研究, 首先考虑质量比m=5.0的拍翼.

图 8所示的是升力系数和动力系数的时间历程曲线, 结果显示升力和能量消耗的峰值均随着拍翼弹性的增加而增加.

为量化其气动表现, 表 2给出了不同弹性参数时的平均升力系数, 平均动力系数和效率.结果表明, 引入弹性后的拍翼其气动表现(升力系数和效率)均可以显著增加, 例如, 当ω=0.4时, 升力系数和效率相对于刚性拍翼(ω=0)分别增加了约62.9%和12.0%.另外, 对比不同弹性参数时的数据可以看出, 当ω=0.3~0.4时, 拍翼的效率最高.另外, 表 2中的数据表明当ω > 0.4时, 动力系数的增加大于升力系数的增加, 因此拍翼效率在弹性较大时反而会降低, 这与文献[4]中的观测结果也是一致的.

为了量化说明柔性拍翼的声学表现, 图 9给出了不同弹性参数时的波动压力平方根的周向分布.结果显示, 声压总体上随着拍翼弹性的增加而增加, 例如, 当ω从0增加到0.5时, 拍翼产生的最大声压从2.42×10^-4增加到4.02×10^-4, 增加约66.1%.综合拍翼的气动表现, 可以看出适当的弹性可以提高拍翼的气动表现同时降低噪音的产生.

3.3 质量比对拍翼声学特性的影响

为研究拍翼-流体质量比对拍翼动力学和声学的影响, 对另一组m=1.0 (作用在拍翼上的气动力和拍翼的惯性相当)时不同刚度的拍翼进行了计算. 图 10所示的是不同弹性参数时的波动压力平方根的周向分布.结果显示, 质量比m=1.0时, 拍翼产生的声压随着ω的增加先增后减, 但总体的变化并不明显.另外, 为了量化说明柔性拍翼的声学表现, 定义声能$I=\int_{\mathit{\Gamma}} \Delta \tilde{p}^{2} \mathrm{d} s$, 这里Γ表示的是距离原点40L处的圆环.

为了比较不同质量比时, 弹性对拍翼气动表现和声学特性得影响, 图 11给出了平均升力, 平均动力系数, 效率和声能随ω变化的曲线.当气动力和拍翼惯性力相当(m=1.0)时, 升力系数和动力系数随着ω的增加先增大后减小, 与惯性力主导时的单调增加不同.效率是升力系数和动力系数综合作用的结果, 图 11表明m=1.0时, 效率同样先增大后减小, 在ω=0.4附近取得最大值.另外, 声能在m=1.0时随着ω的变化不大; 当m=5.0时, 声能随ω单调增加.综合考虑气动和声学表现, 可以看出当ω=0.3~0.4时, 低质量拍翼(即m=1.0)产生的声音较小, 同时又具备较高的效率.

4 结论

本文采用一种浸入边界法对拍翼悬停时的流固声耦合问题进行了数值模拟研究.具体针对刚性拍翼和不同刚度、质量比的柔性拍翼进行了数值模拟, 分析了拍翼刚度和质量比对拍翼悬停时的升力和声学特性的影响.主要结论如下:

(1) 拍翼的转动能有效增加升力, 提高效率并降低拍翼运动产生的声音;

(2) 引入适当的弹性可有效提高拍翼在悬停时的气动表现, 包括提高升力系数和效率;

(3) 拍翼悬停时的近场声受涡量影响十分显著, 增大转动角度时会使得近场声不对称性加剧, 同时加强高频声的产生(如2倍和3倍的拍动频率)；

(4) 综合考虑气动和声学表现, 可以看出当ω=0.3~0.4时, 低质量的拍翼(即m=1.0)产生的声音较小, 同时又具备较高的效率.

致谢 作者感谢Australian Research Council Discovery Early Career Researcher Award (编号: DE160101098)对本文工作的资助.