# Computational aerodynamics of insect flight using volume penalization

Thomas Engels, Hung Truong, Marie Farge, Dmitry Kolomenskiy, Kai Schneider

► To cite this version:

Thomas Engels, Hung Truong, Marie Farge, Dmitry Kolomenskiy, Kai Schneider. Computational aerodynamics of insect flight using volume penalization. Comptes Rendus. Mécanique, 2022, 350 (S1), pp.1-20. 10.5802/crmeca.129. hal-04510332

HAL Id: hal-04510332

https://hal.science/hal-04510332v1

Submitted on 22 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Mécanique

Thomas Engels, Hung Truong, Marie Farge, Dmitry Kolomenskiy and Kai Schneider

**Computational aerodynamics of insect flight using volume penalization**

Published online: 28 November 2022

https://doi.org/10.5802/crmeca.129

**Part of Special Issue:** More than a half century of Computational Fluid Dynamics

**Guest editor:** Mohammed El Ganaoui (Professeur des Universités, Spécialiste : Mécanique des Fluides et Transferts de Chaleur et de Masse, Université de Lorraine)

<figure><img src="image_2.png" /></figure>

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/

<figure><img src="image_3.png" /></figure>

Les Comptes Rendus. Mécanique sont membres du
Centre Mersenne pour l'édition scientifique ouverte

www.centre-mersenne.org

e-ISSN : 1873-7234
<figure><img src="image_4.png" /></figure>

More than a half century of Computational Fluid Dynamics / Plus d'un
demi-siècle de mécanique des fluides numérique

Computational aerodynamics of insect flight
using volume penalization

Thomas Engels<sup>a, ∗</sup>, Hung Truong<sup>b</sup>, Marie Farge<sup>c</sup>,
Dmitry Kolomenskiy<sup>d</sup> and Kai Schneider<sup>∗, b</sup>

<sup>a</sup> ISTA, TU Berlin, Berlin, Germany

<sup>b</sup> I2M-CNRS, Aix-Marseille Université, Marseille, France

<sup>c</sup> LMD-CNRS, Ecole Normale Supérieure, Paris, France

<sup>d</sup> CMT, Skolkovo Institute of Science and Technology, Moscow, Russia

*Current address:* ICube, Université de Strasbourg, Strasbourg, France (H. Truong)

E-mails: thomas.engels@ens.fr (T. Engels), dh.truong@unistra.fr (H. Truong),
marie.farge@ens.fr (M. Farge), d.kolomenskiy@skoltech.ru (D. Kolomenskiy),
kai.schneider@univ-amu.fr (K. Schneider)

**Abstract.** The state-of-the-art of insect flight research using advanced computational fluid dynamics techniques on supercomputers is reviewed, focusing mostly on the work of the present authors. We present a brief historical overview, discuss numerical challenges and introduce the governing model equations. Two open source codes, one based on Fourier, the other based on wavelet representation, are succinctly presented and a mass-spring flexible wing model is described. Various illustrations of numerical simulations of flapping insects at low, intermediate and high Reynolds numbers are presented. The role of flexible wings, data-driven modeling and fluid-structure interaction issues are likewise discussed.

**Keywords.** Numerical simulation, Insect flight, Fluid-structure interaction, Turbulence, Flexible wings.

Published online: 28 November 2022

**1. Introduction**

Insects were the first animals to learn how to fly and they remain the only invertebrates capable of
doing so. Giant dragonflies discovered flapping wing flight 350 million years ago and this remains
the most efficient way for quick maneuvering. Today insects flap their wings at a frequency which
varies from 5 to 1000 Hz depending on the species. Since the beginning of mankind humans quite
probably have wondered why birds and insects are able to fly, but not themselves. Some have
tried to fly, but often at the cost of their lives, as told in the Greek myth of Icarus. In the 16th
century Leonardo da Vinci, by studying birds, determined that flight requires three ingredients:
a light and powerful engine, wings capable of generating sufficient aerodynamic forces, and a
control system to keep the body in the air. The first recognized successes date from the end of the

<footnote>* Corresponding author.</footnote>
19th century, with Otto Lilienthal in 1891 for gliding, Clément Ader in 1897 for the powered flight without trajectory control, and the Wright brothers in 1903 for controlled flight.

In the 1860s, Etienne-Jules Marey, a physician and biomechanist, analyzed the flapping flight of insects and showed that, with each flap, their wings describe a double ellipse similar to a figure eight [1]. To understand how insects carry their weight and move, he built an artificial insect, whose body contains compressed air, which flaps its wings in the shape of the figure eight, and can fly up and down, which he presented at the "Académie des Sciences de Paris" on 15 March 1869. In 1882, he invented chronophotography (the technical basis of cinema) to decompose complex movements, which enabled him to prove that the flight of birds is very similar to that of insects. In 1899, he created a machine that produced 20 smoke trails to visualise the flow of air around a moving body. Through all his inventions, Etienne-Jules Marey played an important role in the birth of aeronautics; in 1874, he was elected vice-president of the Société de Navigation Aérienne, and later president in 1884.

In the 1960s the entomologist Weis-Fogh discovered a new mechanism for insect flight, thanks to the advent of fast cameras, which he called “clap-fling-sweep” [2]. The novelty of this mechanism lies in the change of topology that occurs during downstroke, when their wings separate at their hinge as “fling” gives way to “sweep”, thus producing opposite circulations around the wings during upstroke, until their wings clap together before the next downstroke. He was able to show that various tiny insects, such as the wasp *Encarsia formosa* that he studied, improve their performance and are able to hover using this mechanism, which some larger insects and birds, such as pigeons, also use when they need to escape. In the same year, the mathematician Sir James Lighthill explained the generation of lift produced by this mechanism using only 2D inviscid fluid dynamics [3]. In the 1970s Maxworthy [4] showed experimentally that such a hovering mechanism generates leading edge vortices that are much stronger than those of normal hovering flight. In 2011, Kolomenskiy *et al.* [5] explained the instantaneous generation of lift at wing separation by a viscous effect producing a singularity of the pressure gradient at that moment.

Further significant advances in the field of biofluidynamics were made about fifty years ago with the publication of Sir James Lighthill's seminal monograph on mathematical biofluidynamics [6], and the groundbreaking immersed boundary method-based simulations of flow in a beating heart presented in 1972 by Charles Peskin in his PhD thesis [7, 8]. Since then, the field of computational biofluidynamics has grown considerably, with applications extending to virtually all areas of physiological fluid dynamics, as well as the fluid physics of biolocomotion [9–12]. This expansion has been further boosted by physical discoveries and significant improvements in physics-based, as well as data-driven methodologies, including models that couple flow with structural dynamics, acoustics, and other physical domains. In the field of robotics remarkable progress has been reported, highlighting the developed flapping Nano Air Vehicle (NAV), so called nano hummingbirds [13]. This review aims to highlight the authors' recent advances in computational as well as data-enabled techniques in the field of biofluid dynamics, and to present the state-of-the-art in the application of these methods to study the flows generated by the flapping of insect wings and their biolocomotion.

In the following we give a succinct overview of the state-of-the-art, including some recent advances, in the field of computational fluid dynamics of insect flight. The potential of such numerical simulations lies in their inherent ability to provide the entire instantaneous flow field around the flapping insect, which is difficult to realize experimentally. This difficulty stems from the relatively high flapping frequency, which is ≈200 Hz for a fruit fly [14], ≈150 Hz for a bumblebee [15] and still ≈26 Hz for a hawkmoth [16]. In addition, as the insects are rather small, the flapping wings hinder the observation by the camera and the animal must be motivated to fly in the test section. These difficulties do not appear in numerical simulations, although the complexity of the problem is high. First of all, the geometry of insects is complicated
and therefore appropriate methods must be used that allow to enforce the no-slip boundary conditions on these complicated boundaries. Furthermore, as the flapping motion is dynamic, the geometry is also time-dependent. Additional challenges for numerical simulations are the inherent unsteady and 3D nature of the flow.

The first and most frequent approximation is to assume the wings to be perfectly rigid. Earlier work focused on 2D configurations, see e.g., [11]. Liu and Kawachi [17] were among the first to present 3D computations. Their approach is based on the finite volume method and moving, body-fitted grids. The approach is still used [14, 16], and its advantage is a typically lower memory requirement, at least for small enough Reynolds number *Re*, compared to approaches relying on non-body fitted grids. This advantage, however, comes at a cost of limited CPU parallelism and large numerical dissipation of the small scale turbulent structures in the wake. Ramamurti and Sandberg [18, 19] used the incompressible Navier-Stokes equations in Arbitrary Lagrangian Eulerian (ALE) formulation to simulate hovering fruit flies and a numerical investigation of the timing of stroke reversal described by Dickinson *et al.* [20]. These numerical methods involve significant computational and implementation overhead, whose reduction has motivated the development of methods that allow for a geometry-independent discretization, such as the immersed boundary method (IBM) and related techniques, including the volume penalization method (VPM). The former was established by Peskin [7, 8, 21], while the latter is based on the work of Angot and co-workers [22]. Reviews on different immersed boundary techniques can be found for instance in [21, 23, 24].

The IBM/VPM methods can achieve higher resolutions and performance, and are simpler to implement. However, the boundary conditions are satisfied only approximately. For classical discretization methods, grid-based methods can be used including finite differences, finite volumes or finite elements. Likewise discrete lattice-based methods, i.e., lattice-gas cellular automata leading to lattice-Boltzmann methods are employed, e.g., for simulating butterfly-like flapping wings [25]. Yokoyama *et al.* [26] used the immersed boundary method to simulate a rigid-winged butterfly in free forward flight with a state-of-the-art resolution of 1 billion grid points. They showed that the abdominal motion plays an important role in flight stabilization.

Deformable wings introduce additional complexity to this problem. Not only is the geometry even more complicated, but computations also require suitable solid models and stable fluid-structure interaction (FSI) coupling mechanisms. Here, IBM/VPM methods enjoy an even more significant reduction in implementation and computation time than for rigid wings. Computational studies with flexible wings can be divided into “active” and “passive” flexibility, depending on whether the deformation is actively computed from fluid forces or measured externally and then imposed in the simulation [27, 28]. Active FSI were mostly done in simplified 2D frameworks [29, 30], considering chord-wise flexibility. More recently, Nakata *et al.* presented 3D computations with flexible wings [31, 32], and likewise our team in [33].

The remainder of the article is organized as follows: numerical challenges and different approaches are discussed in Section 2. Fourier and wavelet-based fluid solvers are described and the flexible wing model is introduced. Then we present in Section 3 some numerical results for tiny and also larger insects. Some results for insects with flexible wings are given and high resolution computations of insects flying in the wake of a cylinder and of a fractal tree are shown. Finally some conclusions and future perspectives are presented in Section 4.

## 2. Numerical challenges and different approaches

### 2.1. *Challenges*

There are multiple numerical challenges for highly faithful simulations of insect flight. The spatial scales of range of motion have a broad spectrum and range from the size of the domain
(the numerical wind tunnel, or obstacles like plants generating turbulence), via the insect size, its body and wing dimensions, down to the Kolmogorov scale in the case of turbulent flow. Likewise the temporal scales call for a resolution of the flapping frequency, typically of the order of 10–200 Hz, as well as the time scale of turbulent motion. Moreover, there are FSI issues, the coupling of the air flow and the wing deformation, to be considered. Models for the flexible wings and their motion, i.e., the wing kinematics, are required and material properties must be provided. Since we study living animals, knowledge of biological behavioral issues, including the control mechanisms of insects are necessary, which illustrate the complexity and call for interdisciplinary expert knowledge and high demands on mathematical, computational, biological and engineering scientists.

## 2.2. Models using volume penalization

Modelling flapping insect flight implies that complicated, time varying geometries with FSI need to be considered when solving the aerodynamic equations. For solving the governing partial differential equations (PDEs) in such geometries efficiently, typically immersed boundary techniques, including direct forcing, fictitious domains or penalization approaches, are used. Thus the PDEs or their discretized versions are modified in an adequate way and the boundary conditions, typically no-slip for viscous fluids, are appended or imposed by an additional term. The resulting equations can then be solved in a simple geometry using e.g., Cartesian grids and applying the favorite numerical scheme. Thus all geometrical stiffness is removed, or better hidden in the penalty terms, which have been added to the equations. Here we mostly focus on the volume penalization, originally introduced for fixed geometries in [22, 34], extended for moving obstacles in [35] and generalized for deforming objects in [36]. The main idea of volume penalization (likewise called Brinkman or Darcy penalization) is to use spatially variable permeability to model the solid obstacle (in our case the insect) or solid walls as a porous medium. All geometrical information, which may change in time, is encoded in a mask (or indicator) function $\chi$, illustrated in Figure 1. The penalized versions of the non-dimensionalized incompressible Navier–Stokes equations, which are solved in the larger domain $\Omega = \Omega_f \cup \Omega_s$, where $\Omega_s$ is the solid domain (see Figure 1, left), read

$$ \partial_t \underline{u} + \underline{u} \cdot \nabla \underline{u} + \nabla p - \nu \nabla^2 \underline{u} = - \frac{\chi}{C_{\eta}} (\underline{u} - \underline{u}_s) \quad (1) $$

$$ \nabla \cdot \underline{u} = 0. \quad (2) $$

Here $\chi$ is the indicator function, $C_{\eta} \ll 1$ is the penalization parameter (permeability) and $\underline{u}_s$ the solid body velocity field. The indicator function is $\chi = 1$ inside $\Omega_s$ and $\chi = 0$ in $\Omega_f$, and in the shown computations we use a thin smoothing layer with thickness proportional to the grid spacing in the case of moving [35] or deforming obstacles [36]. The linear damping term in the momentum equation is thus mostly present in the solid domain. In the limit of vanishing permeability, the penalization parameter $C_{\eta}$ tends to zero, and the porous region tends towards a solid wall. Correspondingly the penalized Navier–Stokes equations converge to their non-penalized counterpart imposing no-slip boundary conditions, as shown rigorously in [22, 37]. The observed convergence has a slow rate $O(C_{\eta}^{1/2})$ and thus necessitates sufficiently small, but not too small, permeability values [38]. This approach suffers from several further limitations. A penalization boundary layer is generated at the interface of the solid–fluid region and the penalized equations become stiff not only in space but also in time when $C_{\eta}$ is small. This requires an adequate choice of the parameters, coupling thus the numerical discretization and the penalization parameters, see e.g., [39]. The error can be further reduced using smooth mask functions [35, 38]. Finally, the resulting numerical scheme is low order and yields typically
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Wingbeat of a model bumblebee in forward flight. The wingbeat cycle is divided into two parts, the down- and upstroke. The computational domain $\Omega$ consists of the fluid domain $\Omega_f$ and the solid domain $\Omega_s$ (the insect).</figcaption></figure>

a first, at most a second-order spatial convergence only, at least close to the boundary. An interesting direction for upgrading VPM using smoothing prescription for the damping has been proposed in [40] using multiple-scales matched-asymptotics and applying Richardson extrapolation. Therewith one gets second-order convergence in $C_\eta$ and higher-order spatial convergence can be obtained. An essential aspect of this technique is that the position of the mask functions needs to be known a priori which is for example not the case in FSI. Wake removal techniques for imposing in and outflow conditions based on VPM have been proposed using sponges [38, 41].

Let us mention that there are some generalizations of VPM for imposing Neumann boundary conditions, either homogeneous [42] or inhomogeneous [43] and Robin conditions by linear combination. This is important for passive scalar or heat transport in complex geometries, an example relevant for insects are odour sources (pheromones).

While the penalization method (similarly to the immersed boundary method) thus has some disadvantages, its benefits outweigh them. Besides the relative ease of implementation, those methods are particularly suited for obtaining quick results. While we have methods with higher accuracy available, their computational cost and implementation overhead is larger, as well as the time to set up a simulation. In the field of biofluidynamics, we are often interested in performing many simulations with medium precision rather than a single, very precise one. This stems from the biological variability: no insect looks like the other and their wingbeat is not a strict Fourier series but rather variable. Therefore, VPM is an excellent choice for numerical simulations of insects.

To give an example, Figure 1 illustrates the mask function $\chi$ of a bumblebee model. We show the $\chi = 0.5$ isosurface. In addition, we illustrate the wingbeat during the up- and downstroke. This model insect is in forward flight, thus the stroke plane in which the wings flap is inclined with respect to the horizontal plane. Note that like a helicopter, hovering insects move their wings mainly in a plane parallel to the ground.

## 2.3. *Fourier-based Fluid-Structure Interaction solver: FluSI*

The developed FluSI (Fourier-based Fluid-Structure Interaction solver) open source software package is a fully parallel Fourier pseudospectral method for high resolution computation of 3D flapping flight in viscous flow [41] on massively parallel computers with distributed memory architectures. The penalized incompressible Navier-Stokes equations (2) are discretized with
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Numerical simulation of a bumblebee model with the spectral code FluSI. This code relies on an equidistant grid with constant, fine resolution everywhere which is shown on the left. Only every 128th grid point is shown for visibility. Right part (from [41]) shows a strong scaling test on an IBM BlueGene/Q supercomputer.</figcaption></figure>

a classical Fourier pseudospectral method with adaptive time stepping [44, 45]. Solving a Poisson equation in Fourier space thus becomes trivial and the incompressibility constraint can be easily satisfied, since the Riesz projector is diagonal in Fourier space. Artificial numerical diffusion and dispersion are absent, but the no-slip (and outflow) boundary conditions are imposed with the volume penalization method and modeling error arises. The computational complexity is almost entirely determined by the FFT, for which excellent scaling can be obtained on supercomputers using highly optimized libraries, e.g., P3DFFT. Nevertheless Cartesian equidistant grids limit the number of grid points and the maximum domain size by the available computational resources. The inherent periodicity can be alleviated using penalization-based sponge techniques for wake removal. Moreover the numerics (and the physics) implies small time steps, $\Delta t < C_\eta$, even though Runge-Kutta-Chebychev (RKC) schemes could be included, particularly helpful for low Reynolds number simulations where the penalization method imposes prohibitively small time step sizes.

Different modules are available in the FluSI software package. Flapping flight with rigid wings, either tethered, but a free flight option is likewise available, however without active control so far. Semi-implicitly coupled FSI is also implemented and for handling the solid part of the deformable wings the mass spring-model described in Section 2.5 is used. The modular code structure allows easily to take into account different complex geometries, and thus insects with wings of different shape can be simulated easily. The wing kinematics and flight mechanisms are defined in further modules.

Extensive details on the numerical method including thorough validation and benchmarking can be found in [41]. The FluSI code has been used in particular for computing bumblebees in turbulence [46], for free flight in turbulence [47], FSI of Calliphora wings [48]. In Figure 2 we show a high resolution computation of a tethered bumblebee with FluSI and the corresponding strong scaling test on an IBM Blue Gene with up to 16384 CPUs. Here, a rectangular domain is used, and the reference run is using 1024 CPU cores. Good scaling is observed up to 8192 CPU cores. The software is available on Github at the following address, https://github.com/pseudospectators/FLUSI.
## 2.4. Wavelet-Adaptive Block-Based solver: WABBIT

Motivated by the limitations of the FluSI package, the Fourier discretization requires equidistant grids, which make large scale computations extremely expensive, we developed a second fully parallel code called WABBIT, which stands for Wavelet-Adaptive Block-Based solver for Interactions with Turbulence [49]. It has been designed for computing 3D multiscale flows in complex, time-varying geometries using a similar modular structure as in FluSI. The key feature is adaptivity with cubic blocks and locally equidistant grids endowed with an efficient implementation on last generation massively parallel architectures using tree-like data-structures. The major difference relative to FluSI is that we relax the incompressibility condition and we use instead artificial compressibility. Introducing a finite speed of sound avoids the solution of an elliptic problem but yields another modeling error. The boundary conditions, no-slip and in/outflow are modeled similarly to FluSI with the volume penalization. The governing equations read,

$$ \partial_t \underline{u} + \underline{u} \cdot \nabla \underline{u} + \nabla p - \nu \nabla^2 \underline{u} + \frac{\chi}{C_{\eta}} (\underline{u} - \underline{u}_s) + \frac{\chi_{sp}}{C_{sp}} (\underline{u} - \underline{u}_{\infty}) = 0 \quad (3) $$

$$ \partial_t p + C_0^2 \nabla \cdot \underline{u} + C_{\gamma} p + \frac{\chi_{sp}}{C_{sp}} (p - p_{\infty}) = 0. \quad (4) $$

In the limit when $C_0$ tends to infinity, we recover the incompressible equations (2). The spatial grid is block-structured, cubic-blocks of fixed size with equidistant Cartesian grid, but different step size. For spatial discretization, centered finite differences are used and for time integration, a Runge-Kutta scheme, both of fourth order. Each block of the computational domain is decomposed into biorthogonal wavelets of Cohen-Daubechies-Feauveau type [50]. The wavelet coefficients of the blocks are then used as refinement and prediction operators. Therewith a dynamically evolving grid is generated and an adaption strategy allows to track the solution in space and scale. Numerical analysis justifies the choice of the different modeling parameters (the artificial speed of sound, the porosity for the penalization) and the numerical parameters (thresholding parameter of the wavelet coefficients, space and time steps) to balance the different error contributions and to obtain globally a fourth-order convergent method in space and time. The explicit time marching approach together with the locality of the discretization in space allows to distribute the blocks among different MPI processes. Tree-like data structures and load-balancing with space-filling Hilbert curves then distributes contiguous chunks of blocks to the processes. For further details on the numerical method, its numerical analysis including thorough validation and benchmarking of the parallelization efficiency, see [49].

Figure 3 illustrates an adaptive 3D computation of a bumblebee together with the locally refined grid. The corresponding strong scaling test on different machines of GENCI (Irene SKL and Irene AMD and IDRIS Zay SKL) assesses the speed-up going up to a factor 4. The WABBIT open access software package is available on Github at the following address, https://github.com/adaptive-cfd/WABBIT.

## 2.5. Mass-spring model of flexible insect wings

Insect wings are sophisticated structures which enable insects to generate sufficient lift for flying. Yet, in the past, most research usually simplified insect wings as rigid structures. The wing's dynamic shape change is essential for developing a comprehensive understanding of insect flight because it determines the direction and magnitude of fluid dynamic forces during wing flapping [51, 52]. However, the complex interaction between the surrounding unsteady flow and the anisotropic wing structures renders the analysis of flapping flexible wings demanding. One of the challenges is the computational cost which is already high for the fluid solver. Ideally, the solid model should not increase this cost significantly. Consequently, we choose a mass-spring
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Numerical simulations of a bumblebee model with the adaptive code WABBIT. Parts A and B shows the $|ω| = 50f$ isosurface of vorticity magnitude together with the block-based adaptive grid at 0.3T. Data result from simulations allowing $J_{max} = 6$ (A) and eight (B) levels of refinement. The wavelet-based approach automatically refines the grid where necessary, i.e., at the fluid-solid interface and at locations with important flow features. (C) Strong scaling test for a bumblebee simulation with $J_{max} = 8$ on three different supercomputers (from [49]). Colored labels indicate numbers of CPUs.</figcaption></figure>

model for modelling the wing structure due to its computational efficiency and its ability to deal with large deformations.

The mass-spring model is based on the discretization of an object using mass points which are connected by massless springs. At a given time $t$, the position $\underline{x}_i$ and the velocity $\underline{v}_i$ of the mass point $i$ can be obtained by solving the dynamic equations of the system, i.e., Newton's second law, given by:

$$ \begin{aligned} \underline{F}_i &= m_i \underline{a}_i \\ \underline{F}_i &= \underline{F}_i^{\text{int}} + \underline{F}_i^{\text{ext}} \quad \text{for } i = 1 \dots n \\ \underline{v}_i(t = 0) &= \underline{v}_{0,i} \\ \underline{x}_i(t = 0) &= \underline{x}_{0,i'} \end{aligned} \qquad (5) $$
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Types of springs used in mass-spring models for flexible insect wings.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Left: *Calliphora* wing. Right: Illustration of the mass-spring model which is meshed based on measured data of the real blowfly wing shown on the left. The black and white markers represent mass centers. Color codes (red, green and blue) are used for identifying veins and the membrane is represented by the black triangular mesh. The center of wing mass is shown by the black and white marker. The length is normalized by the fly wing length. The dimensionless vein diameters are displayed by real ratios in the figure. Adapted from [48].</figcaption></figure>

where *n* is the number of mass points, *F*<sub>*i*</sub> is the total force (internal force *F*<sub>*i*</sub><sup>int</sup> and external force *F*<sub>*i*</sub><sup>ext</sup>) acting on the *i*th mass point, *m*<sub>*i*</sub>, *a*<sub>*i*</sub> are mass and acceleration of the *i*th mass point, respectively. All terms in (5) are straightforward to derive, except for the forces. The external forces come from the fluid and gravity, internal forces represent the restoring forces caused by the springs. In order to model flexible insect wings, we combine linear extension springs and bending springs (Figure 4). The former is designed to operate with axial forces resisting extension and compression, while the latter is used for torques resisting bending. The complicated properties of the restoring forces render the system of (5) nonlinear. Because the equations are numerically stiff, in particular because of the torsion springs, we resort to an implicit time marching scheme. The resulting non-linear system is solved using the Newton-Raphson method. More details about the solver can be found in [33].

The mass-spring system is then used for modeling insect wings which are complex structures composed of a network of veins, partly connected with joints, and a membrane spanned in between. Studies have shown that the vein arrangement in insect wings has strong impact on their mechanical properties [51, 53]. Thus, it will be inaccurate to consider a wing as a homogeneous structure; the vein pattern as well as the difference in terms of mechanical behaviors between vein and membrane need to be taken into account. Mechanically, veins can be considered as rods that resist torsion and bending. On the other hand, a membrane is fabric-like and resists stretching. By employing different kinds of springs, these anisotropic and inhomogeneous properties of insect wings can be included [54]. The vein structure, as well as the wing contour, are adapted from experimental measurement and encoded into the mass-spring network, as shown in Figure 5 for a blowfly (a large housefly, *Calliphora*) wing.
The wing morphology allows us to estimate the mass distribution by assuming that veins are solid rods made of cuticle with density $ρ_c = 1300 \text{ kg} \cdot \text{m}^{-3}$ as explained in [54]. However, wing stiffness remains unknown since their structural properties are still poorly understood. In [48], we propose a numerical method to evaluate wing stiffness by using a genetic algorithm with covariance matrix adaptation. The model is trained to reproduce static elasticity measurements of cantilevered *Calliphora* wings under point forces.

## 2.6. Fluid-structure coupling

To perform FSI simulation of insects with flexible wings, the solid model, in our case the developed mass-spring model, is coupled with the fluid solver, either FluSI or WABBIT in our case, see e.g., [33]. Physically the coupling conditions consist in imposing the no-slip conditions of the velocity and the continuity of the mechanical stress tensor at the fluid-solid surface. The fluid solver requires the wing position, which is given by the mask function $\chi$ and the velocity $u_s$. The solid solver requires the pressure distribution at the wing surface and interpolation becomes necessary since the wing position does not coincide with the fluid grid. In the semi-implicit staggered scheme we use, the fluid is displaced first, while keeping the wing fixed. Subsequently we obtain the new forces on the wing using the fluid field at the new time step. These are then transferred as input into the mass-spring model to move the wing to the next time level and to construct the updated mask function. This staggered approach implies that the dynamic coupling condition is fulfilled but not the kinematic one, i.e., the continuity of the mechanical stress tensor. This results in a limited stability of the scheme which depends on the mass ratio and very light wing structures cannot be considered. For such configurations a strongly coupled scheme needs to be used, see e.g., [55], or iterative coupling, the latter is implemented in FluSI using fixed point iteration with Aitken relaxation. We recall that an implicit time scheme is applied for the solid model and an explicit scheme for the fluid equations, however the difference in the time discretization schemes causes some difficulty in coupling. For details on the FSI coupling we refer to [33, 36].

# 3. Numerical results

Various high resolution numerical simulations using either the FluSI or the WABBIT code developed by the authors are presented. First results for tiny insects are presented corresponding to flows at low Reynolds numbers; thereafter results for larger insects are shown. Either rigid or flexible wings are considered and for the aerodynamics we study laminar and turbulent inlet flow. Finally an adaptive simulation of a bumblebee behind a fractal tree illustrates the capabilities of the a wavelet-based solver for computing multiscale problems and the simulation of a bumblebee in the cylinder wake yields insights into the free flight dynamics.

## 3.1. Tiny insects

We may perceive all insects as small creatures, but their variability in size is tremendous. The smallest ones are smaller than unicellular organisms like amoeba, yet they possess sophisticated locomotory organs and a neural system. In a recent study [56], a minute beetle *Paratuposa placentis* was considered as a model organism. Its body is less than half a millimeter long. Like other beetles, its forewings evolved in rigid elytra that serve to protect the more fragile hindwings, which are used for flying. The hindwings of *P. placentis* have narrow membranous central blades and thin long bristles that radiate from it. This wing condition is typical of the smallest representatives of different orders of insects, but the advantages of possessing such bristled wings remained uncertain.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Flight of a tiny beetle *Paratuposa placentis*. Left: overall view of the morphological model. Right: flow visualization using vorticity magnitude iso-surfaces $|ω| = 10f$ and $30f$, where *f* is the wing flapping frequency.</figcaption></figure>

The study [56] combined thorough 3D morphological and kinematic reconstructions with fluid and solid body dynamics calculations. Geometrical models of the insect body, elytra and wings (Figure 6, left) were constructed based on scanning electron microscope data. Special attention was paid to ensure that the positions, orientations and effective diameter of the bristles in the wing model are representative of the real wings. This was important because a prior 2D analysis [57] suggested that the ratio of the inter-bristle spacing over the bristle diameter is an important aerodynamic parameter. To respect this geometric similarity condition in the CFD simulation was challenging, because the bristle diameter was three hundred times smaller than the wing length. At the same time, the external flow domain should be much larger than the wing length at this relatively low Reynolds number. Hence, the dynamic grid adaptation implemented in WABBIT [49] became crucial for the feasibility of these numerical simulations.

The peculiar features of *P. placentis*, such as the bristled wings, figure-of-eight wing tip trajectories, dorsal and ventral clapping, body pitch oscillation and elytra opening and closing are pre-conditioned by the small size of this insect and all combine to achieve flight performance competitive with larger species. Bristled wings are lightweight, therefore, inertial torques are negligible even when the wings flap with the maximum possible amplitude and clap at reversals. The necessary vertical force is produced as the wings displace almost flat-on downwards. Edge vortices (Figure 6, right) indicate that the flow through the wing is slowed down by the viscous stresses. Wing rotation and clapping produce the kinematic asymmetry needed for propulsion at low Reynolds number. A side effect of the large-amplitude flapping is a large periodic pitching torque that is counteracted by the elytra. Numerical simulation allowed to explore the parameter space beyond the insects' capabilities. Thus, the effects of varying the Reynolds number and the number of bristles were explored [58, 59].

## 3.2. Larger insects

Kolmogorov scale in the lower atmosphere may range from 0.1 to 10 mm, therefore, large insects experience turbulence as a multiscale phenomenon. Turbulence becomes a limiting factor for the direct numerical simulation of the flight of large insects. In the following discussion we focus our attention on bumblebees that are large enough to produce a wide range of vortical motions in the wake but are still amenable to Navier-Stokes simulation without statistical modelling of turbulence.

Bumblebees are relentless foragers that can fly under strong gusty wind conditions [60]. The effect of inflow turbulence on the aerodynamic performance of the flapping bumblebee wings was studied in [46]. A bumblebee model with a fixed body and flapping wings was placed
<figure><img src="image_11.png" /><figcaption>**Figure 7.** FluSI computation of a tethered bumblebee in a virtual wind tunnel and turbulent inflow with intensity $Tu = 0.99$. Shown are vorticity magnitude $|\omega|$ and the mask function. For details see [46].</figcaption></figure>

in a virtual wing tunnel. The inflow boundary condition was generated using pre-computed homogeneous isotropic turbulence data samples (Figure 7). It was found that the time-averaged aerodynamic lift, drag, torques and power were remarkably insensitive to the inflow turbulence intensity (Figure 8). However, the variance of all of those parameters increased with increasing inflow turbulence intensity. This means that, when the body is free to rotate, flying through turbulence may require additional control effort. Uncontrolled free flight simulations [47] show that the roll angle perturbations grow in time, followed by pitch and yaw reorientation. The growth rate correlates positively with the inflow turbulence intensity and length scale.

In the field, foraging bumblebees encounter turbulence generated by the vegetation canopy. A typical situation of a bumblebee approaching a flower was modelled in [61] by introducing a vertically oriented cylindrical obstacle in an otherwise unperturbed inflow. By confronting the results of numerical simulations and animal flight experiments in a wind tunnel, it was shown that bumblebees can ride through the small-scale high-frequency turbulence without any assistance from the neural system. Active control is used to compensate for the slowly growing reorientation of the body, and it acts on the same time scale as for maneuvering.

### 3.3. *Flexible wings*

The most basic and common approach to model the flapping wing motion consists in prescribing the Euler angles of rigid wings as functions of time. To obtain these data, synchronized high-speed video recordings using multiple cameras are necessary. Wing-tip trajectories are easily amenable to automated tracking, but measuring the wing rotation about its longitudinal axis is particularly demanding because it requires 3D reconstruction of the spatial orientation of the wing plane. On the other hand, this feathering rotation occurs largely due to the mechanical
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Bumblebee simulations with different levels of inflow turbulence intensity *Tu* from 0 to 1. Statistical moments of aerodynamic measures are shown as mean values with the 95% confidence intervals ± standard deviation. Forces are normalized by the body weight *m*·*g*, torques by *m*·*g*·*R*, where *g* is gravity, *m* body mass and *R* the wing length. Aerodynamic power is in Watt per kilogram body mass.</figcaption></figure>

compliance of the wing hinge and the supporting skeletomuscular elements. This has lead to passive feathering rotation approximation using lump flexibility models previously developed for robotic flappers [62] being combined with the computational fluid dynamics solver FluSI. The results of numerical simulations of bumblebee hovering show that this reduced-order model of FSI can produce realistic feathering motion and accurate time-average estimates of the aerodynamic performance [63].

The flexible wing models provide important insights into the influence of flexibility on the aerodynamic performance of insects. The study [64] compared the aerodynamic forces and the power requirements between tethered bumblebees with highly flexible, flexible and rigid wings. The visualization in Figure 9 shows a FluSI computation of a tethered bumblebee with flapping flexible wings. The wing deformation is clearly visible for the left wing, and the complex vortical structures can be observed on the right wing and in the wake. It showed that flexibility allows reducing the energetic cost of flapping flight characterized by a higher lift-to-power ratio. However, the highly flexible wing appears to be less efficient than the flexible wing. This can be interpreted as the existence of an optimized zone of wing flexibility, which is ideal for flying. Furthermore, the wing inertia contributes to the damping of fluctuations in the aerodynamic forces and hence stabilizes the insect during flight.

Experimental data on insect wing stiffness is essential for modelling correctly wing flexibility. To this purpose, an experimental setup for measuring the wing elasticity of female (*Calliphora vomitoria*) was designed by Wehmann *et al.* [65]. The measured data were used for training a data-driven model of insect wings as proposed in [48]. The new approach ensured that the model had the same behavior as the real wing. We obtained overall nine sets of stiffness parameters corresponding to the nine measured individuals. In the same computational setup at the same
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Visualization of flow generated by a tethered flapping bumblebee with flexible wings in laminar flow at $Re = 2685$ showing normalized absolute vorticity isosurfaces $|ω| = 100$ (light blue). The flow field is plotted at time $t/T = 0.45$. The vortices are only shown in the right half of the computational domain for the purpose of visualizing the deformation of the left wing. Adapted from [64].</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Wing deformation during one flapping cycle viewed from the side in the wing system. The flexible wing (white) is plotted together with the rigid wing (blue) to indicate the deformation caused by inertial and aerodynamic forces. The arrows show the flapping direction during the downstroke (top) and the upstroke (bottom). Adapted from [48].</figcaption></figure>

Reynolds number, the influence of intra-species variability of wing stiffness on their aerodynamic performance was insignificant. Furthermore, by looking at the wing deformation (cf. Figure 10), we observed that the maximum deflection of the wing leading edge occurred at the wing tip during the mid-downstroke and corresponds to 10% of the wing length. On the other hand, during the reversals, the wing was deformed only in chord-wise direction due to strong inertial force caused by the wing rotation.

### 3.4. Fractal tree

To illustrate the flexibility of WABBIT to compute flows with multiple scales, we simulate the flow of the bumblebee model, described previously, in the wake of a model flower, corresponding to a
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Bumblebee behind a fractal tree. Top part shows the setup consisting of a bumblebee and a tree-inspired fractal turbulence generator composed of rigid cylinders. Bottom part illustrates the flow field with an isosurface of the Q-criterion. Time in figure is $t/T = 8.0$, results obtained with CDF 4/0 wavelets. Adapted from [49].</figcaption></figure>

bio-inspired turbulence generator. The latter is composed of an array of rigid cylinders arranged in a fractal manner, largely motivated by plants. For a detailed description of the geometry we refer the reader to [66]. Both bumblebee and fractal tree are simulated using the volume penalization method in a large cubic domain of size $L = 64R$, where $R$ is the wing length of the insect. Figure 11 (top, left) illustrates the setup. The thinnest cylinders in the fractal tree have a diameter $d_{min} = 0.114R$ and thus we have a Reynolds number $Re = u_{\infty}d_{min}/\nu = 250$. The other numerical parameters are chosen as in the coarse simulations of the bumblebee alone, for details we refer to [49]. The simulation is performed using CDF 4/0 wavelets where $C_{\epsilon} = 4 \times 10^{-2}$. The flow field is visualized in Figure 11 with an isosurface of the Q-criterion. It shows the turbulent wake produced by the tree model and the wake of the bumblebee after eight wingbeats. The maximum number of refinement levels is $J_{max} = 9$ and the grid is composed, on average, of $N_b = 1.8 \times 10^5$ blocks with a total number of $2.2 \times 10^9$ grid points. The dynamically adapted grid allows us to use a large domain size, which would result in 11776³ points on the corresponding, uniform grid with the same $\Delta x$. We hence use on average only 0.15% of the uniform grid.

### 3.5. *Free flight*

In [61] we studied the bumblebee's motion while flying in the von Karman vortex street generated in the wake of a circular cylinder, at $Re = 4200$ (Figure 12). The high resolution numerical simulations with more than 2.2 billion grid points run on 8192 cores allowed to study the free flight dynamics of flapping bumblebees considering two degrees of freedom, lateral displacement and roll rotation about the longitudinal axis of the body. We compared the results with experimental data and assessed the presence of active versus passive flight control.

Free flight simulations of bumblebees taking into account all six degrees of freedom in a numerical wind tunnel have been performed to study the influence of isotropic turbulent inflow with variable intensity [47], a configuration similar to that in Figure 7. Active control was excluded and only the passive response of real animals was quantified. We showed that changes in body orientation and angular velocity are highly sensitive to variations in the turbulence spectrum. Moreover, we found the translation of the insect to be small compared to its rotational motion.
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Bumblebee (colored in orange) behind a circular cylinder at $Re = 4300$ (top view). Shown are isosurfaces of vorticity magnitude, light blue $|ω| = 5f$, red $|ω| = 40f$. The velocity field in the horizontal plan of the insect is represented by arrows. From [61].</figcaption></figure>

**4. Conclusions and perspectives**

This review outlines the authors' work over the last ten years on the numerical simulation of insect flight, mainly related to the French-German AIFIT project (Aerodynamics of Insect Flight in Turbulent Flow), jointly funded by ANR and DFG. After a brief motivation for biolocomotion, in particular for flight, we gave an overview of the history of computational aerodynamics of insect flight, and then mostly focused on our own work. We reviewed the fully parallel codes we have developed for solving the 3D FSI problems related to flapping insect flight: FluSI using a Fourier pseudospectral discretization, and WABBIT using wavelet-based adaptive methods. The various high-resolution computations of tiny and larger size insects illustrated the state-of-the-art and scientific capabilities of supercomputing to study the aerodynamics of flapping flight. The first adaptive computations with WABBIT using the geometry of real insects obtained from micro computed tomography data are promising and will be published soon.

Current limitations of the computations of FSI are imposed by the computational complexity in terms of memory and CPU time, even when using adaptive locally refined grids with the WABBIT code. The available computational resources limit the Reynolds number to values below 5000. This implies that the size of the insects that can be simulated remains limited, e.g., the hawk moth (*Manduca sexta*) with a typical wingspan of 100 mm and flying at 5 m/s has a Reynolds number of about 5000 [67]. Today, it is not possible to simulate bird flight without using a turbulence model to reduce the resolution requirements, e.g., by using large eddy simulations (LES), because not all active scales can yet be solved by direct numerical simulation (DNS). An interesting perspective currently being explored is the coherent vorticity simulation (CVS) using wavelet-based adaptive grids, where the incoherent background flow produced as a thermal noise during the turbulent flow evolution is eliminated, since it corresponds to turbulent dissipation [68, 69].
The solid model needed to describe the wing deformation requires input parameters, which are based on measurements and optimization. Our current approach uses mass-spring models, which are distinguished by their simplicity, and directly takes into account the wing structure consisting of veins and membranes. Data-driven methods can then be used to identify the parameters, as done in [48]. A straight forward extension is to upgrade the mass-spring model by including torsion springs. This would allow to take into account the effect of dynamic torsion which many species exhibit, however its impact on insect flight is still unclear [70]. In the future, even more sophisticated approaches may be considered, such as finite element models, but again, close interaction with biologists is needed to determine the material properties required in the governing equations and the benefits in terms of computational cost need to be assessed.

Control issues for free flight simulations are unavoidable and certainly advanced machine learning tools, in particular deep reinforcement learning approaches, shall be exploited, including for insect swarms [71].

Possible challenging next steps for future investigations are the study of aeroacoustics of insects and also computational aerodynamics of flapping flight on Mars for designing flapping-wing nano air vehicles.

## Conflicts of interest

The authors declare no competing financial interest.

## Dedication

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

## Acknowledgments

Financial support from the Agence Nationale de la Recherche (ANR) (Grant 15-CE40-0019) and Deutsche Forschungsgemeinschaft (DFG) (Grants SE824/26-1 and LE905/17-1), project AIFIT, is gratefully acknowledged. The authors were granted access to the HPC resources of IDRIS under the allocation No. A0102A01664 attributed by Grand Équipement National de Calcul Intensif (GENCI). DK gratefully acknowledges financial support from the JSPS KAKENHI Grants Number 15F15061 and 18K13693. Centre de Calcul Intensif d'Aix-Marseille, Global Scientific Information and Computing Center at Tokyo Institute of Technology and Center of Shared Research Facilities at Lomonosov Moscow State University are acknowledged for granting access to their high performance computing resources.

## References

[1] E. J. Marey, "Determination experimentale du mouvement des ailes des insectes pendant le vol", C. R. Acad. Sci. Paris 67 (1868), p. 1341-1345.

[2] T. Weis-Fogh, "Quick estimates of flight fitness in hovering animals, including novel mechanisms for lift production", J. Exp. Biol. 59 (1973), p. 169-230.

[3] M. Lighthill, "On the Weis-Fogh mechanism of lift generation", J. Fluid Mech. 60 (1973), no. 1, p. 1-17.

[4] T. Maxworthy, "Experiments on the Weis-Fogh mechanism of lift generation by insects in hovering flight. Part 1. Dynamics of the 'fling'", J. Fluid Mech. 93 (1979), no. 1, p. 47-63.

[5] D. Kolomenskiy, H. K. Moffatt, M. Farge, K. Schneider, "The Lighthill-Weis-Fogh clap-fling-sweep mechanism revisited", J. Fluid Mech. 676 (2011), p. 572-606.
[6] S. J. Lighthill, *Mathematical Biofluiddynamics*, SIAM, Philadelphia, PA, 1975.

[7] C. S. Peskin, "Numerical analysis of blood flow in the heart", *J. Comput. Phys.* **25** (1977), p. 220-252.

[8] C. S. Peskin, "Flow patterns around heart valves: a digital computer method for solving the equations of motion", PhD Thesis, The Albert Einstein College of Medicine, New York, NY, 1972.

[9] T. Maxworthy, "The fluid dynamics of insect flight", *Annu. Rev. Fluid Mech.* **13** (1981), no. 1, p. 329-350.

[10] S. P. Sane, "The aerodynamics of insect flight", *J. Exp. Biol.* **206** (2003), p. 4191-4208.

[11] Z. J. Wang, "Dissecting insect flight", *Annu. Rev. Fluid Mech.* **37** (2005), p. 183-210.

[12] Z. J. Wang, "Insect flight: from Newton's law to neurons", *Annu. Rev. Condens. Matter Phys.* **7** (2016), p. 281-300.

[13] T. A. Ward, M. Rezadad, C. J. Fearday, R. Viyapuri, "A review of biomimetic air vehicle research: 1984-2014", *Int. J. Micro Air Veh.* **7** (2015), no. 3, p. 375-394.

[14] M. Maeda, H. Liu, "Ground effect in fruit fly hovering: A three-dimensional computational study", *J. Biomech. Sci. Eng.* **8** (2013), p. 344-355.

[15] R. Dudley, C. P. Ellington, "Mechanics of forward flight in bumblebees I. Kinematics and morphology", *J. Exp. Biol.* **148** (1990), p. 19-52.

[16] H. Liu, "Integrated modeling of insect flight: From morphology, kinematics to aerodynamics", *J. Comput. Phys.* **228** (2009), p. 439-459.

[17] H. Liu, K. Kawachi, "A numerical study of insect flight", *J. Comput. Phys.* **146** (1998), p. 124-156.

[18] R. Ramamurti, W. C. Sandberg, "A three-dimensional computational study of the aerodynamic mechanisms of insect flight", *J. Exp. Biol.* **205** (2002), p. 1507-1518.

[19] R. Ramamurti, W. Sandberg, "A computational investigation of the three-dimensional unsteady aerodynamics of Drosophila hovering and maneuvering", *J. Exp. Biol.* **210** (2009), p. 881-896.

[20] M. H. Dickinson, F.-O. Lehmann, S. P. Sane, "Wing rotation and the aerodynamic basis of insect flight", *Science* **284** (1999), p. 1954-1960.

[21] C. S. Peskin, "The immersed boundary method", *Acta Numer.* **11** (2002), p. 479-517.

[22] P. Angot, C. Bruneau, P. Fabrie, "A penalization method to take into account obstacles in incompressible viscous flows", *Numer. Math.* **81** (1999), p. 497-520.

[23] R. Mittal, G. Iaccarino, "Immersed boundary methods", *Annu. Rev. Fluid Mech.* **37** (2005), p. 239-261.

[24] K. Schneider, "Immersed boundary methods for numerical simulation of confined fluid and plasma turbulence in complex geometries: a review.", *J. Plasma Phys.* **81** (2015), article no. 435810601.

[25] K. Suzuki, K. Minami, T. Inamuro, "Lift and thrust generation by a butterfly-like flapping wing-body model: immersed boundary-lattice Boltzmann simulations", *J. Fluid Mech.* **767** (2015), p. 659-695.

[26] N. Yokoyama, K. Senda, M. Iima, N. Hirai, "Aerodynamic forces and vortical structures in flapping butterfly's forward flight", *Phys. Fluids* **25** (2013), article no. 021902.

[27] L. Zheng, T. L. Hedrick, R. Mittal, "Time-varying wing-twist improves aerodynamic efficiency of forward flight in butterflies", *PLoS ONE* **8** (2013), no. 1, article no. e53060.

[28] J. Young, S. M. Walker, R. J. Bomphrey, G. K. Taylor, A. L. R. Thomas, "Details of insect wing design and deformation enhance aerodynamic function and flight efficiency", *Science* **325** (2009), no. 5947, p. 1549-1552.

[29] D. Kolomenskiy, T. Engels, K. Schneider, "Numerical modelling of flexible heaving foils", *J. Aero Aqua Bio-Mech.* **3** (2013), p. 22-28.

[30] L. A. Miller, C. S. Peskin, "Flexible flap and fling in tiny insect flight", *J. Exp. Biol.* **212** (2009), p. 3076-3090.

[31] T. Nakata, H. Liu, "A fluid-structure interaction model of insect flight with flexible wings", *J. Comput. Phys.* **231** (2012), no. 4, p. 1822-1847.

[32] T. Nakata, H. Liu, "Aerodynamic performance of a hovering hawkmoth with flexible wings: a computational approach", *Proc. R. Soc. B* **279** (2012), no. 1729, p. 722-731.

[33] H. Truong, T. Engels, D. Kolomenskiy, K. Schneider, "A mass-spring fluid-structure interaction solver: application to flexible revolving wings", *Comput. Fluids* **200** (2020), article no. 104426.

[34] E. Arquis, J.-P. Caltagirone, "Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide milieu poreux: application à la convection naturelle", *C. R. Acad. Sci. Paris, Sér. II* **299** (1984), p. 1-4.

[35] D. Kolomenskiy, K. Schneider, "A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles", *J. Comput. Phys.* **228** (2009), p. 5687-5709.

[36] T. Engels, D. Kolomenskiy, K. Schneider, J. Sesterhenn, "Two-dimensional simulation of the fluttering instability using a pseudospectral method with volume penalization", *Comput. Struct.* **122** (2012), p. 101-112.

[37] G. Carbou, P. Fabrie, "Boundary layer for a penalization method for viscous incompressible flow", *Adv. Differ. Equ.* **8** (2003), p. 1453-2480.

[38] T. Engels, D. Kolomenskiy, K. Schneider, J. Sesterhenn, "Numerical simulation of fluid-structure interaction with the volume penalization method", *J. Comput. Phys.* **281** (2015), p. 96-115.

[39] R. Nguyen van yen, D. Kolomenskiy, K. Schneider, "Approximation of the laplace and stokes operators with Dirichlet boundary conditions through volume penalization: A spectral viewpoint", *Numer. Math.* **128** (2014), p. 301-338.
[40] E. W. Hester, G. M. Vasil, K. J. Burns, "Improving accuracy of volume penalised fluid-solid interactions", *J. Comput. Phys.* **430** (2021), article no. 110043.

[41] T. Engels, D. Kolomenskiy, K. Schneider, J. Sesterhenn, "FluSI: A novel parallel simulation tool for flapping insect flight using a Fourier method with volume penalization", *SIAM J. Sci. Comput.* **38** (2016), no. 5, p. S3-S24.

[42] B. Kadoch, D. Kolomenskiy, P. Angot, K. Schneider, "A volume penalization method for incompressible flows and scalar advection-diffusion with moving obstacles", *J. Comput. Phys.* **231** (2012), p. 4365-4383.

[43] T. Sakurai, K. Yoshimatsu, N. Okamoto, K. Schneider, "Volume penalization for inhomogeneous Neumann boundary conditions modeling scalar flux in complicated geometry", *J. Comput. Phys.* **390** (2019), p. 452-469.

[44] C. Canuto, M. Y. Hussaini, A. Quarteroni, T. Zang, *Spectral Methods in Fluid Dynamics*, Springer, Berlin, Heidelberg, 1986.

[45] K. Schneider, "Numerical simulation of the transient flow behaviour in chemical reactors using a penalisation method", *Comput. Fluids* **34** (2005), p. 1223-1238.

[46] T. Engels, D. Kolomenskiy, K. Schneider, F.-O. Lehmann, J. Sesterhenn, "Bumblebee flight in heavy turbulence", *Phys. Rev. Lett.* **116** (2016), article no. 028103.

[47] T. Engels, D. Kolomenskiy, K. Schneider, M. Farge, F.-O. Lehmann, J. Sesterhenn, "Impact of turbulence on flying insects in tethered and free flight: High-resolution numerical experiments", *Phys. Rev. Fluids* **4** (2019), article no. 013103.

[48] H. Truong, T. Engels, H. Wehmann, D. Kolomenskiy, F.-O. Lehmann, K. Schneider, "An experimental data-driven mass-spring model of flexible Calliphora wings", *Bioinspir. Biomim.* **17** (2022), no. 2, article no. 026003.

[49] T. Engels, K. Schneider, J. Reiss, M. Farge, "A wavelet-adaptive method for multiscale simulation of turbulent flows in flying insects", *Commun. Comput. Phys.* **30** (2021), p. 1118-1149.

[50] A. Cohen, I. Daubechies, J. C. Feauveau, "Biorthogonal bases of compactly supported wavelets", *Commun. Pure Appl. Math.* **45** (1992), p. 485-560.

[51] S. A. Combes, T. L. Daniel, "Flexural stiffness in insect wings I. Scaling and the influence of wing venation", *J. Exp. Biol.* **206** (2003), p. 2979-2987.

[52] L. Zhao, Q. Huang, X. Deng, S. P. Sane, "Aerodynamic effects of flexibility in flapping wings", *J. R. Soc. Interface* **7** (2009), no. 44, p. 485-497.

[53] S. A. Combes, T. L. Daniel, "Flexural stiffness in insect wings II. Spatial distribution and dynamic wing bending", *J. Exp. Biol.* **206** (2003), p. 2989-2997.

[54] H. Truong, T. Engels, D. Kolomenskiy, K. Schneider, "Fluid-structure interaction using volume penalization and mass-spring models with application to flapping bumblebee flight", in *Cartesian CFD Methods for Complex Applications* (R. Deiterding, M. O. Domingues, K. Schneider, eds.), SEMA SIMAI Springer Series, vol. 3, Springer, Cham, 2021, p. 19-35.

[55] C. Foerster, W. Wall, E. Ramm, "Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows", *Comput. Methods Appl. Mech. Eng.* **196** (2007), p. 1278-1293.

[56] S. E. Farisenkov, D. Kolomenskiy, P. N. Petrov, T. Engels, N. A. Lapina, F.-O. Lehmann, R. Onishi, H. Liu, A. A. Polilov, "Novel flight style and light wings boost flight performance of tiny beetles", *Nature* **602** (2022), p. 96-100.

[57] S. H. Lee, M. Lee, D. Kim, "Optimal configuration of a two-dimensional bristled wing", *J. Fluid Mech.* **888** (2020), article no. A23.

[58] D. Kolomenskiy, S. Farisenkov, T. Engels, N. Lapina, P. Petrov, F.-O. Lehmann, R. Onishi, H. Liu, A. Polilov, "Aerodynamic performance of a bristled wing of a very small insect", *Exp. Fluids* **61** (2020), article no. 194.

[59] T. Engels, D. Kolomenskiy, F.-O. Lehmann, "Flight efficiency is a key to diverse wing morphologies in small insects", *J. R. Soc. Interface* **18** (2021), no. 183, article no. 20210518.

[60] J. D. Crall, J. J. Chang, R. L. Oppenheimer, S. A. Combes, "Foraging in an unsteady world: bumblebee flight performance in field-realistic turbulence", *Interface Focus* **7** (2017), no. 1, article no. 20160086.

[61] S. Ravi, D. Kolomenskiy, T. Engels, K. Schneider, C. Wang, J. Sesterhenn, H. Liu, "Bumblebee minimize control challenges by combining active and passive modes in unsteady winds", *Sci. Rep.* **6** (2016), article no. 35043.

[62] J. P. Whitney, R. J. Wood, "Aeromechanics of passive rotation in flapping flight", *J. Fluid Mech.* **660** (2010), p. 197-220.

[63] D. Kolomenskiy, S. Ravi, R. Xu, K. Ueyama, T. Jakobi, T. Engels, T. Nakata, J. Sesterhenn, K. Schneider, R. Onishi, H. Liu, "The dynamics of passive feathering rotation in hovering flight of bumblebees", *J. Fluids Struct.* **91** (2019), article no. 102628.

[64] H. Truong, T. Engels, D. Kolomenskiy, K. Schneider, "Influence of wing flexibility on the aerodynamic performance of a tethered flapping bumblebee", *Theor. Appl. Mech. Lett.* **10** (2020), no. 6, p. 382-389.

[65] H.-N. Wehmann, L. Heepe, S. N. Gorb, T. Engels, F.-O. Lehmann, "Local deformation and stiffness distribution in fly wings", *Biol. Open* **8** (2019), no. 1, article no. bio038299.

[66] M. Dreissigacker, "Turbulence generated by fractal Trees - PIV measurements and comparison with numerical data", PhD Thesis, echnische Universität Berlin, 2017.

[67] A. P. Willmott, C. P. Ellington, "The mechanics of flight in the Hawkmoth Manduca sexta II. Aerodynamic consequences of kinematic and morphological variation", *J. Exp. Biol.* **200** (1997), p. 2723-2745.
[68] M. Farge, K. Schneider, N. Kevlahan, "Non-Gaussianity and coherent vortex simulation for two-dimensional turbulence using an adaptive orthonormal wavelet basis", *Phys. Fluids* **11** (1999), p. 2187-2201.

[69] M. Farge, K. Schneider, "Coherent vortex simulation (CVS), a semi-deterministic turbulence model using wavelets", *Flow Turbul. Combust.* **66** (2001), p. 393-426.

[70] S. A. Combes, "Materials, structure, and dynamics of insect wings as bioinspiration for MAVs", in *Encyclopedia of Aerospace Engineering*, vol. 7, John Wiley & Sons, Ltd, Chichester, West Sussex, UK, 2010, Part 34.

[71] S. Hong, S. Kim, D. You, "Control of a fly-mimicking flyer in complex flow using deep reinforcement learning", 2021, arXiv preprint, arXiv:211103454.