2.1 Wing kinematics

A flapping wing with chord length $c$ and span width $s$ is depicted in Fig. 1. A global space-fixed coordinate system OXYZ and a local body-fixed coordinate system $o{x_0}{y_0}{z_0}$ are adopted for a more concise description of the motion. Assuming that the wing is in forward flight with air speed of $U$ , and that the incoming flow direction is along the X-axis, the propulsive performance of the wing is studied considering the flapping, as well as the chordwise deformation. The flapping is described by the heaving $h(t)$ in the Z-axis direction and the pitching $\theta (t)$ around the ${y_0}$ -axis.

Figure 1. Flapping wing model and its coordinates.

An efficient flapping motion^{(Reference Liu, Li, Zhao and Su23)} is adopted in the present study. The heaving $h(t)$ and pitching $\theta (t)$ are described as follows:

(1) \begin{align} h(t) & = {h_0}c\cos( {2\pi ft})\end{align}

(2) \begin{align} \theta (t) & = -{\theta _0}\sin( {2\pi ft})\end{align}

where ${h_0}$ and ${\theta _0}$ denote the amplitude of heaving and pitching, respectively. $f$ denotes the flapping frequency, and $t$ denotes the time.

The bending of the wing has a significant influence on its aerodynamic characteristics, so chordwise deformation that varies with time is introduced into the flapping. Assuming that the bending motion has the same frequency as the flapping, the bending of the chordwise deformable wing is then expressed as

(3) \begin{align}\Delta h(t) & = {{ - a(t)c} \over {{{({1 - p})}^2}}}\left[ {{{\left( {{{{x_0}} \over c}} \right)}^2} - 2p{{{x_0}} \over c} - {p^2}} \right]\left( {p \lt {{{x_0}} \over c} \lt 1} \right)\end{align}

(4) \begin{align} a(t) = {a_0}c\cos({2\pi ft + \psi })\qquad\qquad\qquad\qquad\qquad\end{align}

where $f$ is the bending frequency, $p$ denotes the dimensionless chord length where the wing starts to bend, ${a_0}$ denotes the dimensionless bending amplitude and $\psi $ denotes the phase lag between the flapping and bending.

2.2 Immersed-boundary lattice Boltzmann method

Instead of solving the Navier–Stokes equations, the lattice Boltzmann method (LBM) solves the discrete Boltzmann equation to simulate the airflow using a collision model. The fundamental quantity in the LBM is the discrete-velocity distribution function ${f_i}({\mathbf{x},t})$ , where ${f_i}$ is also known as the particle population. The LBM makes use of ${f_i}({\mathbf{x},t})$ to solve the lattice Boltzmann equation. The macroscopic flow density $\rho $ and momentum $\boldsymbol{\rho} \textbf{u}$ at position $\mathbf{x}$ and time $t$ are replaced with the following forms:

(5) \begin{align}\rho & = \sum\limits_i \,{{f_i}}\end{align}

(6) \begin{align}\boldsymbol{\rho} \textbf{u} & = \sum\limits_i {\textbf{c}_i{f_i}} \end{align}

where ${\mathbf{c}_i}$ is the discrete velocity of a mesoscopic particle along the i-th direction. For a given lattice space $\Delta x$ and time step $\Delta t$ , the discrete velocity vectors are depicted as follows:

(7) \begin{equation}{\mathbf{c}_i} = \left( {{c_{ix}}{{\Delta x} \over {\Delta t}},{c_{iy}}{{\Delta x} \over {\Delta t}},{c_{iz}}{{\Delta x} \over {\Delta t}}} \right)\end{equation}

In the present study, a model consisting of a three-dimensional lattice with 19 velocity vectors (known as the D3Q19 model) is adopted. The 19 velocity vectors in the D3Q19 model are depicted in Table 1. The lattice Boltzmann equation is then

(8) \begin{equation}{f_i}( {\mathbf{x} + {\mathbf{c}_i}\Delta t,t + \Delta t} ) - {f_i}({\mathbf{x},t}) = - {1 \over \tau }\left[ {{f_i}({\mathbf{x},t}) - f_i^{eq}({\mathbf{x},t})} \right]\end{equation}

where $f_i^{eq}$ denotes the local equilibrium distribution function and $\tau $ is the dimensionless relaxation time. The equilibrium distribution function is expressed as

(9) \begin{equation}f_i^{eq} = {w_i}\rho \left[ {1 + {{{\mathbf{c}_i} \cdot \mathbf{u}} \over {c_s^2}} + {{{{( {{\mathbf{c}_i} \cdot \mathbf{u}} )}^2}} \over {2c_s^4}} - {{{\mathbf{u}^2}} \over {2{c_s}}}} \right]\end{equation}

where $\mathbf{u}$ is the fluid velocity, ${w_i}$ is the weight factor and ${c_s}$ is the speed of sound. The speed of sound ${c_s}$ in LBM^{(Reference Chen and Doolen26)} is expressed as

(10) \begin{equation}{c_s} = {1 \over {\sqrt 3 }}{{\Delta x} \over {\Delta t}}\end{equation}

The weight factors in the D3Q19 model are also presented in Table 1.

The dimensionless relaxation time is calculated from the kinematic viscosity $\nu $ of the incompressible viscous fluid^{(Reference KrÜger, Kusumaatmaja and Kuzmin27)}:

(11) \begin{equation}\tau = {\nu \over {c_s^2\Delta t}} + {1 \over 2}\end{equation}

When considering the force applied on the wing, the external body force is applied to Equation (8). The relation between the external body force $\mathbf{g}(\mathbf{x},t)$ and the discrete-velocity distribution ${f_i}$ in LBM is depicted as follows:

(12) \begin{align}f_i^*( {\mathbf{x} + {\mathbf{c}_i}\Delta t,t + \Delta t} ) & = {f_i}({\mathbf{x},t}) - {1 \over \tau }[ {{f_i}({\mathbf{x},t}) - f_i^{eq}({\mathbf{x},t})} ]\qquad\quad\end{align}

(13) \begin{align}{f_i}( {\mathbf{x},t + \Delta t} ) & = f_{i}^*( {\mathbf{x},t + \Delta t} ) + 3\Delta x{\mathbf{c}_i}{w_i} \cdot \mathbf{g}( {\mathbf{x},t + \Delta t} )\end{align}

The LBM is coupled with the immersed boundary to make full use of the Eulerian coordinate and Lagrangian coordinate, which are suitable for the moving surface. In the present study, the LBM with a multi-direct forcing immersed boundary^{(Reference Wang, Fan and Luo28,Reference Suzuki and Inamuro29)} is adopted. The position is then depicted by $\mathbf{x}$ in Eulerian coordinates, and by $\mathbf{X}$ in Lagrangian coordinates. The fluid velocity at the k-th Lagrangian points ${\mathbf{X}_k}$ is then determined by

(14) \begin{equation}{\mathbf{u}^*}( {{\mathbf{X}_k},t + \Delta t} ) = \sum\limits_\mathbf{x} {{\mathbf{u}^*}( {\mathbf{x},t + \Delta t} )\delta ( {\mathbf{x} - {\mathbf{X}_k}} ){{( {\Delta x} )}^3}} \end{equation}

where ${\mathbf{u}^*}({\mathbf{x},t + \Delta t})$ denotes the fluid velocity at Eulerian points x, which is determined by Equations (6) and (12).

Table 1 Parameters of the D3Q19 model

The kernel interpolation function $\delta ({x,y,z})$ between the Eulerian and Lagrangian coordinates is expressed as^{(Reference Peskin30)}

(15) \begin{equation}\delta ({x,y,z}) = \phi (x) \cdot \phi (y) \cdot \phi (z) \cdot {1 \over {{{({\Delta x})}^3}}}\end{equation}

where

(16) \begin{equation}\phi (r) = \begin{cases} \dfrac{1}{8}\left( {3 - 2\left| r \right| + \sqrt {1 + 4\left| r \right| - 4{r^2}} } \right) & {\left| r \right| \le \Delta x}\\[12pt] \dfrac{1}{8}\left( {5 - 2\left| r \right| - \sqrt { - 7 + 12\left| r \right| - 4{r^2}} } \right)& {\Delta x \le \left| r \right| \le 2\Delta x} \\[12pt] 0 & otherwise \end{cases}\end{equation}

The body force $\mathbf{g}({{\mathbf{X}_k},t + \Delta t})$ at the k-th Lagrangian points ${\mathbf{X}_k}$ is then determined by

(17) \begin{equation}\mathbf{g}( {{\mathbf{X}_k},t+\Delta t} )={{{\mathbf{U}_{\textit{\textbf{k}}}}( {{\mathbf{X}_{\textit{\textbf{k}}}},t+\Delta t}) - {\mathbf{u}^*}( {{\mathbf{X}_k},t+\Delta t})} \over {\Delta t}}\end{equation}

where ${\mathbf{U}_k}({{\mathbf{X}_k},t + \Delta t})$ denotes the fluid velocity at Lagrangian points ${\mathbf{X}_k}$ .

The relation between the body force $\mathbf{g}({\mathbf{x},t + \Delta t})$ in Eulerian coordinates and the body force $\mathbf{g}({{\mathbf{X}_k},t + \Delta t})$ in Lagrangian coordinates is

(18) \begin{equation}\mathbf{g}({\mathbf{x},t + \Delta t}) = \sum\limits_{k = 1}^N {\mathbf{g}({{\mathbf{X}_k},t + \Delta t})\delta ({\mathbf{x} - {\mathbf{X}_k}}){S_k}} \end{equation}

where ${S_k}$ is the correlative surface area at Lagrangian points ${\mathbf{X}_k}$ .

The body force is calculated by an iterative approach^{(Reference Suzuki and Inamuro29,Reference Kotsalos, Latt and Chopard31)} . With the body force, the total force on the wing is then obtained as

(19) \begin{equation}\mathbf{F}(t) = - \sum\limits_\mathbf{x} {\mathbf{g}({\mathbf{x},t}){{({\Delta x})}^3}} \end{equation}

The total torque on the wing is

(20) \begin{equation}\mathbf{T}(t) = - \sum\limits_\mathbf{x} {( {\mathbf{x} - {\mathbf{X}_C}(t)})\mathbf{g}({\mathbf{x},t}){{({\Delta x})}^3}} \end{equation}

where ${\mathbf{X}_C}(t)$ is the centre of rotation.

For the whole flapping process, the velocities, total forces and total torques at every time step are obtained by the IB-LBM simulation. In the present study, we use the open-source LBM-solver Palabos^{(Reference Latt, Malaspinas, Kontaxakis, Parmigiani, Lagrava, Brogi, Belgacem, Thorimbert, Leclaire, Li, Marson, Lemus, Kotsalos, Conradin, Coreixas, Petkantchin, Raynaud, Beny and Chopard32)} for the simulation.

2.3 Propulsive performance

The propulsive performance of a flapping wing includes the thrust, energy dissipation and efficiency. The thrust coefficient is defined as

(21) \begin{equation}{C_T} = {{ - {F_x}(t)} \over {0.5\rho {U^2}{A_w}}}\end{equation}

where ${A_w}$ denotes the projected area of the wing and ${A_w} = c \cdot s$ .

The average thrust coefficient is then determined as

(22) \begin{equation}{\bar C_T} = {1 \over T}\int_{t + T}^t {{C_T}(t)dt} \end{equation}

where $T$ is the period time of flapping.

The energy dissipation for the deformable flapping wing is defined as

(23) \begin{equation}P(t) = - \int_{{S_w}} {\boldsymbol{\sigma} (t)} \cdot {\mathbf{u}_{local}}(t)dS\end{equation}

where $\boldsymbol{\sigma} (t)$ is the stress on the local surface at time $t$ , ${\mathbf{u}_{local}}(t)$ is the local velocity of wing at time $t$ and ${S_w}$ is the total surface of the wing.

The energy dissipation coefficient for the deformable flapping wing is defined as

(24) \begin{equation}{C_P} = {{P(t)} \over {0.5\rho {U^3}{A_w}}}\end{equation}

The average energy dissipation coefficient is

(25) \begin{equation}{\bar C_P} = {1 \over T}\int_{t + T}^t {{C_P}(t)dt} \end{equation}

With the average thrust coefficient and the average energy dissipation coefficient, the propulsive efficiency can be defined as

(26) \begin{equation}{\eta _T} = {{{{\bar C}_T}} \over {{{\bar C}_P}}}\end{equation}

Besides the propulsive performance, there are some other parameters that are important for the flapping wing. These include the lift coefficient, average velocity, vorticity and Q-criteria. The lift coefficient ${C_L}$ , average velocity $\bar{\mathbf{U}}$ and vorticity ${\boldsymbol{\omega}^{\rm{*}}}$ are defined as follows:

(27) \begin{align}{C_L} & = {{{F_z}(t)} \over {0.5\rho {U^2}{A_w}}}\end{align}

(28) \begin{align}\bar{\textbf{U}} & = {1 \over T}\int_{t + T}^t {\mathbf{u}dt}\end{align}

(29) \begin{align}{\boldsymbol{\omega}^*} & = \nabla \times \mathbf{u}\qquad\end{align}

The Q-criteria^{(Reference Chakraborty, Balachandar and Adrian33)} is used to identify the three-dimensional vortex structure^{(Reference Guerrero34,Reference Cheng, Roll, Liu, Troolin and Deng35)} , where Q is defined in incompressible flows as

(30) \begin{equation}Q = {1 \over 2}\left( {\left\| {{\boldsymbol{\Omega} ^2}} \right\| - \left\| {{\mathbf{S}_r}^2} \right\|} \right)\end{equation}

where ${\mathbf{S}_r} = {1 \over 2}\left[ {\nabla \mathbf{u} + {{( {\nabla \mathbf{u}} )}^T}} \right]$ is the strain rate tensor and $\boldsymbol{\Omega} = {1 \over 2}\left[ {\nabla \mathbf{u} - {{( {\nabla \mathbf{u}} )}^T}} \right]$ denotes the rotation of the fluid.

3.1 Three-dimensional numerical model

A three-dimensional flapping wing with NACA0012 aerofoil is studied. The chord and aspect ratio of the wing are chosen as $c = 1\text{m}$ and $s/c = 2$ . The wing is in forward flight with speed of $U = 1{\rm{m/s}}$ . The Reynolds number is $\text{Re} = {{Uc} \over \nu } = 200$ . The Strouhal number is $St = {{2{h_0}cf} \over U} = 0.5~1.0$ when the flapping frequency $f$ ranges from 0.5Hz to 1.0Hz. The length scale of the flow field in the simulation is set to be $13.5c \times 9c \times 9c$ to avoid wall effects. The inlet and outlet of the flow field are assumed to have a constant flow velocity $U$ . Periodic boundary conditions are adopted for the flow field. A heaving amplitude of 0.5 and a pitching amplitude of 30° are chosen in the simulation. The flapping wing model and the simulated flow field are shown in Fig. 2. The simulation is carried out on a 144-core parallel high-performance computer provided by Xi’an Jiaotong University.

3.2 Mesh convergence

Mesh convergence is studied by comparing the propulsive results for the rigid flapping wing model with different mesh densities. Parameters such as the flapping frequency $f$ , heaving amplitude ${h_0}$ , pitching amplitude ${\theta _0}$ , Reynolds numbers, and bending amplitude of the model are presented in Table 2.

The flow field is meshed with uniform grid distance. The number of grids in a chord length can be used to quantify the mesh density. Eight flapping wing models with grid numbers of n = 20, 25, 30, 35, 40, 45, 50 and 55 are established to study the mesh convergence. The results of these simulation are shown in Fig. 3. The convergence of the average thrust coefficient and the propulsive efficiency is obvious. Considering that the errors in both the average thrust coefficient and propulsive efficiency for the models with n = 45, n = 50 and n = 55 are less than 2%, the mesh density with n = 45 is adopted for further simulations.

3.3 Verification

The present IB-LBM simulation is further verified by comparing the simulation results with data published by Dong et al.^{(Reference Dong, Mittal and Najjar22)}. The aspect ratio of the rigid flapping wing is $s/c = 1$ . The thickness of the wing is 0.12 $c$ . Other parameters include ${\theta _0} = {30^ \circ }$ , $St = 0.6$ , $\text{Re} = 200$ and ${h_0} = 0.5$ .

Table 2 The parameters of the rigid flapping wing

Figure 4 shows the variation of the lift and thrust coefficients with time. The results from Dong et al.^{(Reference Dong, Mittal and Najjar22)} are also depicted in Fig. 4. It is seen that both the lift and thrust coefficients from the present simulation agree with Dong’s data very well. The difference between the present results and Dong’s data is less than 3.7% for the lift coefficient and less than 4% for the thrust coefficient.

Figure 2. CFD model for the chordwise-deformable flapping wing.

Figure 3. The result convergence with grid numbers: (a) average thrust coefficient and (b) propulsive efficiency.

Table 3 Chordwise deformation parameters and simulating cases

Figure 4. The numerical results for a rigid flapping wing: (a) lift coefficient and (b) thrust coefficient.