Hung Truong
email: dinh-hung.truong@univ-amu.fr

Thomas Engels
email: thomas.engels@ens.fr

Dmitry Kolomenskiy
email: dkolomenskiy@jamstec.go.jp

Kai Schneider
email: kai.schneider@univ-amu.fr

A mass-spring fluid-structure interaction solver: application to flexible revolving wings

Keywords: insect flight, wing elasticity, mass-spring model, fluid-structure interaction, spectral method, volume penalization method

The secret to the spectacular flight capabilities of flapping insects lies in their wings, which are often approximated as flat, rigid plates. Real wings are however delicate structures, composed of veins and membranes, and can undergo significant deformation. In the present work, we present detailed numerical simulations of such deformable wings. Our results are obtained with a fluidstructure interaction solver, coupling a mass-spring model for the flexible wing with a pseudo-spectral code solving the incompressible Navier-Stokes equations. We impose the no-slip boundary condition through the volume penalization method; the time-dependent complex geometry is then completely described by a mask function. This allows solving the governing equations of the fluid on a regular Cartesian grid. Our implementation for massively parallel computers allows us to perform high resolution computations with

Introduction

A fundamental characteristics of insect flight are flexible wings, which play an important role for their aerodynamics [START_REF] Young | Details of insect wing design and deformation enhance aerodynamic function and flight efficiency[END_REF]2,3], requiring lower forces than their rigid counter parts and producing reduced sound. Numerical simulation of insect flight is itself a sophisticated task, because it involves the solution of fluid-solid interaction problems. Thus, we have to model simultaneously the fluid part and the solid part by using two coupled solvers.

The fluid solver [START_REF] Engels | Flusi: A novel parallel simulation tool for flapping insect flight using a Fourier method with volume penalization[END_REF] we are using in the present work has been developed previously and is called FLUSI 1 , a fully parallel software dedicated for modeling three-dimensional flapping flight in viscous flows. The heart of this software is the Fourier pseudospectral method with adaptive time stepping used for the discretization of the 3D incompressible Navier-Stokes equations. Moreover, the volume penalization method is used to take into account the no-slip boundary conditions on the interfaces between the fluid and the solid part.

In [START_REF] Kolomenskiy | Leading-edge vortex shedding from rotating wings[END_REF], we performed numerical simulations of rotating triangular rigid wings at Reynolds number Re = 250 to investigate the Leading-Edge Vortices (LEVs) as a function of the wing aspect ratio and the angle of attack. High resolution direct numerical simulations of rotating and flapping bumblebee wings were presented in [START_REF] Engels | Helical vortices generated by flapping wings of bumblebees[END_REF] using likewise the FLUSI code with rigid wings focusing again on the role of LEVs and the associated helicity production.

The interaction of flapping bumblebees with turbulent inflow in free and tethered flight was studied in [START_REF] Engels | The impact of turbulence on flying insects in tethered and free flight: high-resolution numerical experiments[END_REF] using once again massively parallel computations with FLUSI. We found that the fluctuations of aerodynamic observables significantly grow with increasing turbulence intensity, even if the mean values are almost unaffected. Changing the length scale of the turbulent inflow, while keeping the turbulence intensity fixed, showed that the fluctuation level of forces and moments can be significantly reduced. We also found that the scale-dependent energy distribution in the surrounding turbulent flow is a relevant factor conditioning how flying insects control their body orientation.

Nevertheless, a solver for solving the deformation of the structure was not fully developed in this software FLUSI. Hence, all previous simulations of insect flight have been performed with the essential assumption that the insect is composed of linked rigid parts including the wings.

In the current work we aim at investigating the role of wing elasticity on the flight performance. Consequently, a wing model is required for simulating the deformation of the solid part under the impact from the fluid. There are many models based on continuum mechanics theory, which are used in many well-known solid solvers. Among these, Finite Element Methods (FEM) are mostly used in both research and industrial fields due to their reliability and effectiveness. Despite this dominance, the requirement for faster and more efficient methods motivates the development of alternative models. One of these is the mass-spring system, which is known for its computational efficiency and straightforward implementation [START_REF] Jarrousse | Modified mass-spring system for physically based deformation modeling[END_REF]. As part of this work, a solver using a network of masses and springs is developed to model a flexible insect wing and coupled with FLUSI.

The story of mass-spring models started back at the end of the 20 th century when people were dealing with observable deformations of flexible objects in computer graphics applications, such as soft tissue, skin, hair, ball, cloth, textiles, etc. Being considered as one of the pioneers on this field, Terzopoulos et al. [START_REF] Terzopoulos | Elastically deformable models[END_REF] proposed using elasticity theory for modeling deformable objects instead of previously-used kinematic models, where the motion and the deformation of materials were prescribed. During this period, most of simulations of flexible objects had been calculated using finite element methods until the requirement of faster models was claimed by Eischen et al. [START_REF] Eischen | Finite element modeling and control of flexible fabric parts[END_REF].

Consequently, the development of physically based deformable models started to grow, especially in the field of computer graphics [START_REF] Cai | Modeling and dynamics simulation for deformable objects of orthotropic materials[END_REF] and biomedical engineering [START_REF] Jarrousse | Modified mass-spring system for physically based deformation modeling[END_REF]. While classic solvers, based on finite elements or finite difference methods, are generally employed in the static case for computing stress and strain in a structure, new solvers for deformable objects must have the ability to deal with large deformations and large displacements, i.e. the nonlinear regime. Furthermore, these models need to be fast enough since they are usually employed for real-time applications or coupled with other models which are already time-consuming. Among all the deformable models proposed, mass-spring networks stand out as the most intuitive and simplest due to their computational efficiency [START_REF] Nealen | Physically based deformable models in computer graphics[END_REF]. Mass-spring systems have already been employed in many fields such as medical applications [START_REF] Jarrousse | Modified mass-spring system for physically based deformation modeling[END_REF] (muscle, red blood cells and virtual surgery), computer graphics, fluid-structure interaction and insect flight. Miller and Peskin [START_REF] Miller | Flexible clap and fling in tiny insect flight[END_REF] to capture all the dynamic motion at small time scale and length scale. Instead, a revolving wing model is usually employed to study the aerodynamics of flapping wings thanks to its simplicity. Accordingly, the flow fields and force generation aspects of revolving wings have been analyzed for a wide range of parameters, as reported in the literature [START_REF] Lentink | Rotational accelerations stabilize leading edge vortices on revolving fly wings[END_REF][START_REF] Jardin | Spanwise gradients in flow speed help stabilize leading-edge vortices on revolving wings[END_REF][START_REF] Jones | Biomimetic flexible plate actuators are faster and more efficient with a passive attachment[END_REF]. Di et al. [START_REF] Di | Forewings match the formation of leading-edge vortices and dominate aerodynamic force production in revolving insect wings[END_REF] studied the role of forewings in generating LEVs of three revolving insect including some perspectives for future investigations.

Flexible wing

Mass-spring model

The very basic idea of the mass-spring model is to discretize an object using mass points m i (i = 1, ..., n) connected by springs. There exist many kinds of springs for different purposes but in the limit of our work, we have used only linear extension and bending springs to model insect wings. The dynamic behavior of the mass-spring system, at a given time t, is defined by the position x i and the velocity v i of the mass point i. For this, we need to solve the dynamic equations of the system, which govern their motion in time under certain external forces. This is one of the elegant advantages of the mass-spring network where these governing equations are simply the corresponding classical well-known Newton's second law, given in eqns. [START_REF] Young | Details of insect wing design and deformation enhance aerodynamic function and flight efficiency[END_REF].

F i = m i a i F i = F int i + F ext i for i = 1...n v i (t = 0) = v 0,i x i (t = 0) = x 0,i (1)
where n is the number of mass points, F i is the total force (internal force

F int i
and external force F ext i ) acting on the i th mass point, m i , a i are mass and acceleration of the i th mass point, respectively.

Here, all terms are quite simple to derive except for the forces. The external forces come from outside of the system and depend on the surrounding field and the problem we are dealing with. On the other hand, the internal forces represent the restoring forces caused by the springs. The complicated properties of these forces make the system (1) become nonlinear. Hence, we have a nonlinear system of 3n second order ordinary differential equations (ODEs) corresponding to three dimensions x, y, z and n mass points. In the general case, this system (1) needs to be solved numerically since its analytical solution cannot be derived. In practice, it is more convenient to convert a system of 3n second order equations into a system of 6n first order equations by using the relations a i = dv i /dt and v i = dx i /dt. Eqns. (1) can then be rewritten as below:

dx i dt = v i m i dv i dt = F int i + F ext i for i = 1...n v i (t = 0) = v 0,i x i (t = 0) = x 0,i (2)
Let us call q = x i , v i the phase vector containing positions and velocities of all mass points and f (q

) = v i , m -1 i (F int i + F ext i
) the right hand side function, eqns.( 2) can be rewritten again under the familiar form of a system of first order ODEs as follows:

dq dt = f (q, t) q(t = 0) = q 0 (3)
For time stepping, eqns. (3) need to be discretized using an appropriate numerical scheme. The choice for this scheme is not trivial since it depends on many factors. Due to the size of the wings, the mass-spring network contains a lot of very small springs, which make the system very stiff and we need an implicit scheme for time marching. For this reason, either a centered scheme or a backward scheme can be used. Although centered schemes are usually in favor for their conservation property without numerical diffusion, a centered scheme, for example the trapezoidal scheme, can lead to numerical instability at some points because the eigenvalues of the operator of the time discretization lie exactly on the imaginary axis, the boundary of the stability zone [START_REF] Engels | Numerical modeling of fluid-structure interaction in bioinspired propulsion[END_REF]. Furthermore, the coupling between the fluid and the solid solver will require an adaptive time stepping scheme. For all these reasons, a second order backward differentiation scheme with variable time steps [START_REF] Berger | A second order backward difference method with variable steps for a parabolic problem[END_REF] is used to discretize eqns. (3) in time.

q n+1 i - (1 + ξ) 2 1 + 2ξ q n i + ξ 2 1 + 2ξ q n-1 i = 1 + ξ 1 + 2ξ ∆t n f (q n+1 i ) (4)
where ξ = ∆t n /∆t n-1 is the ratio between the current ∆t n and the previous time step ∆t n-1 . Eqn. (4) is a system of nonlinear equations with the variable q n+1 , the phase vector of the system at the current time step, which needs to be solved. The Newton-Raphson method, a powerful iterative method, is employed to solve this nonlinear system of equations. With an initial guess, which is reasonably close to the true root of the equations, Newton-Raphson helps to find approximations of the root with the rate of convergence estimated to be quadratic. For our mass-spring solver, the initial guess chosen is the phase vector q n solved at the previous time step; this allows the Newton-Raphson method to converge quickly since the structure is advanced slowly and smoothly in time, which assures that q n+1 remains close to q n . In most cases, the Newton-Raphson method in the solver needs three to four iterations to converge within a relative or absolute L 2 norm error of 10 -6 as the stopping criterion.

Extension springs and bending springs

Extension springs (figure 1a) and bending springs (figure 1b) are common mechanical devices, which resist against the external forces to get back to their resting positions. The former is designed to operate with axial forces, while the latter is used for torques. The relations between the displacement and the restoring force are given by:

• Linear extension spring

F i = k e ij (||x ij || -||x 0,ij ||) e ij F j = -F i (5)
where k e ij is the extension stiffness, e ij = (x j -x i )/||x j -x i || is the unit position vector and F i and F j are the restoring forces of the extension spring acting on two points i and j, respectively;

• Linear bending spring

M ijk = -k b ijk (θ ijk -θ 0,ijk ) (6)
or in terms of forces

F i = k b ijk (θ ijk -θ 0,ijk ) (e ij × e jk ) × e ij F k = k b ijk (θ ijk -θ 0,ijk ) (e ij × e jk ) × e jk F j = -F i -F k (7)
where M ijk is the restoring moment, k b ijk is the bending stiffness, θ 0,ijk is the initial bending angle among three points i, j and k, θ ijk the current bending angle and F i , F j , F k are the restoring force vectors (as shown in figure 1) of the bending spring acting on three points i, j and k, respectively.

However, the calculation of θ ijk in eqn. [START_REF] Engels | The impact of turbulence on flying insects in tethered and free flight: high-resolution numerical experiments[END_REF] is not trivial since it involves the geometrical definition of the angle in 3D space with respect to a fixed coordinate system. Firstly, we consider a simpler case when the three points are in the same plane, a 2D coordinate system Oxy. This leads to x i = (x i , y i ) T , x j = (x j , y j ) T and x k = (x k , y k ) T . The angle is now determined by:

θ ijk = atan2 (y k -y j )(x j -x i ) -(y j -y i )(x k -x j ), (x k -x j )(x j -x i ) + (y k -y j )(y j -y i ) (8)
Here, atan2 (usually known as two-argument arctangent) is a special function first introduced in computer programming languages to give a correct and unambiguous value for the angle by taking into account the sign of both arguments. This function helps us to calculate on the whole space when the angle can vary in the range of (-π, π], instead of the range of (-π/2, π/2) when using the standard arctangent function.

For a problem in 3D space, only one bending angle will not be sufficient.

This can be easily seen by considering a simple case of one bending spring.

At an instant t, the spring is deformed and has a bending angle θ ijk . Nevertheless, corresponding to this θ ijk , there is an infinite number of solutions x i , x j and x k that can satisfy this deformation and the set of all these solutions forms a cone, like shown in figure 2. Consequently, one more angle is needed to obtain a unique solution. To define these two bending angles, the same bending spring as the 2D case above is considered but now in a 3D coordinate system Oxyz. The bending spring is projected on the Oxy and the Oxz planes which gives us two 2D bending angles θ y ijk and θ z ijk on the Oxy and the Oxz planes, respectively. Then, like in the 2D problem, these two bending angles are calculated as below:

θ y ijk = atan2 (x j -x i )(y k -y j ) -(x k -x j )(y j -y i ), (x j -x i )(x k -x j ) + (y j -y i )(y k -y j ) θ z ijk = atan2 (x j -x i )(z k -z j ) -(x k -x j )(z j -z i ), (x j -x i )(x k -x j ) + (z j -z i )(z k -z j ) (9)

Functional approach -vein and membrane models

Modeling insect wings is challenging due to the fact that these wings have complex structures composed of a network of veins, partly connected through hinges, with thin membranes spanned in between and their elasticity properties are still poorly understood. Certain studies have shown that the vein arrangements in insect wings have strong impact on their mechanical properties [3,[START_REF] Combes | Flexural stiffness in insect wings II. Spatial distribution and dynamic wing bending[END_REF]. 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. However this is not an easy task, since insects are the most diverse group of animals living on Earth with more than one million known species with varying wing sizes and wing shapes. As a result, in this study, we want to limit ourself to a specific case when we examine only bumblebee (Bombus ignitus). Bumblebee wings are mainly composed of veins and membranes.

The longitudinal veins extending along the wing in the spanwise direction are big, hollow and providing conduits for nerves while the cross veins are smaller, solid and connecting the longitudinal veins to form a kind of truss structure. In the space between these veins, we find a thin layer called membrane.

From a mechanical point of view, veins can be considered as rods which resist mainly the torsion and bending deformation. On the other hand, a membrane is fabric-like, it behaves like a piece of cloth which resists again the extension deformation. Consequently, instead of considering the wing as a homogeneous structure, a functional approach is used to distinguish veins and membranes. We then propose two models using mass-spring networks to imitate the mechanical behavior of the vein and the membrane.

A vein is considered as a rod whose length is much greater than its height and width. The total effect of all the external loads applied on a mechanical structure results in deformation which can usually be classified into three main types: bending, twisting and stretching. Although the whole wing is observed to be twisted significantly in many studies using high speed cameras or the digital particle image velocimetry [START_REF] Walker | Deformable wing kinematics in the desert locust: how and why do camber, twist and topography vary through the stroke?[END_REF][START_REF] Walker | Photogrammetric reconstruction of high-resolution surface topographies and deformable wing kinematics of tethered locusts and free-flying hoverflies[END_REF][START_REF] Walker | Deformable wing kinematics in free-flying hoverflies[END_REF][START_REF] Bomphrey | The aerodynamics of manduca sexta: digital particle image velocimetry analysis of leading-edge vortex[END_REF][START_REF] Mountcastle | Aerodynamic and functional consequences of wing compliance[END_REF], it is not entirely clear that torsion happens locally at veins or the unsynchronized bending deformations between veins cause the whole wing to twist. To simplify the model, we study only the latter in which we ignore the local torsion of veins and model solely the bending deformation of veins by using extension and bending springs. Thus, we model a vein by a curve line which is discretized by n mass points as shown in figure 3. Two neighboring points are connected with each other by an extension spring (e.g. the mass points i and i + 1 are connected by the extension spring k e i ) and three neighboring points are connected with each other by a bending spring (e.g. the mass points i -1, i and i + 1 are connected by the bending spring k b i-1 ). When flapping, most of external forces will act in the direction perpendicular to the wing surface. As a result, the stretching deformation is negli-gible comparing to the bending deformation. Thus, the role of the extension springs in the model is solely to conserve the length of the vein. The stiffness of these extension springs is artificial and they do not need to reflect the mechanical property of the vein itself. They should be chosen stiff enough to make the rod unstretched but not too stiff to avoid problems with numerical stability when integrating the dynamical system in time.

Compared to veins, membranes are totally different in terms of geometrical and mechanical properties. Geometrically, a membrane is an object whose thickness is much smaller than its extent. Consequently, a membrane is usually considered as a planar two-dimensional sheet or a set of planar sheets in the case of non-planar three-dimensional membranes [START_REF] Fenner | Engineering Elasticity: Application of Numerical and Analytical Techniques[END_REF][START_REF] White | An Introduction to the Finite Element Method with Applications to Nonlinear Problems[END_REF]. On the mechanical side, a membrane is a special kind of structure compared to other structural elements, i.e. a rod, a bar, a plate or a beam. It behaves like a piece of cloth which is much easier to be bent than to be stretched or compressed. Keeping these in mind, the membrane part of the wing is modeled by a 2D sheet which is discretized by a system of mass points and extension springs. There are several ways to discretize a 2D sheet (as shown in figure 4) but an unstructured triangular mesh needs to be employed for our problem due to the complicated geometry of insect wing. Moreover, an unstructured mesh is preferred for modeling isotropic membranes [START_REF] Chen | Investigation of membrane mechanics using spring networks: application to red-blood-cell modelling[END_REF] since the random orientations of the springs will average out the forces.

Correlation between mass-spring network models and continuum constitutive laws

Besides the mesh topology, the parameter setting is another challenge that one has to solve in order to correctly model the material of which the object is made. There are two main parameters needed to be assigned for a mass-spring model: the masses and the spring stiffness. Although a Voronoi diagram can be used to find the masses in a proper way [START_REF] Deussen | Using simulated annealing to obtain good nodal approximations of deformable bodies[END_REF], selecting spring stiffness is still an open question and there are two common solutions to overcome this [START_REF] Lloyd | Identification of spring parameters for deformable object simulation[END_REF]. The first approach is based on optimization methods to minimize the difference between the results solved by the mass-spring model and the reference data. These reference data can come from the measurements, the visual appearance of real objects [START_REF] Louchet | Evolutionary identification of cloth animation models[END_REF] or numerical solutions using validated methods such as finite element methods [START_REF] Bianchi | Mesh topology identification for mass-spring models[END_REF][START_REF] Bianchi | Simultaneous topology and stiffness identification for mass-spring models based on fem reference deformations[END_REF]. In general, this approach cannot be applied if the system has too many degrees of freedom with many unknown spring constants or the mesh topology varies in time since one set of parameters works for solely one mesh structure. Otherwise, tuning the spring stiffness by using optimization can give satisfying results with reasonable computational cost.

The second way is about deriving a relation between spring stiffness and other continuum mechanic properties, such as Young modulus, the Poisson ratio and the flexural rigidity. In contrast to the discrete models, the elasticity parameters have been obtained for many materials and can be used to calculate the corresponding spring stiffness. Omori et al. [START_REF] Omori | Comparison between spring network models and continuum constitutive laws: Application to the large deformation of a capsule in shear flow[END_REF] succeeded in doing this for a planar membrane by considering a 2D sheet under small uniaxial deformation. The relation between spring networks and continuum models for three types of meshes is shown in figure 5.

For the vein model, a relation between the flexural rigidity EI and the stiffness of the bending springs k b i is needed. To derive this relation, we consider a classical problem of a cantilever beam length l b , under a point force F at the free end (figure 6). In the limit of small displacement, the principle of energy yields the value of the bending spring stiffness k b as a function of the flexural rigidity EI. The energy stored in this beam at the static state can be calculated easily using the Euler-Bernoulli beam theory as shown in eqn. [START_REF] Eischen | Finite element modeling and control of flexible fabric parts[END_REF].

E beam = F 2 l 3 b 6EI (10)
The mass-spring network is called an equivalent model of the continuous beam if under the same external loads, its mechanical behavior (in this case, it is the energy stored in the system) is the same as the one of the beam. for the remaining n points, we have:

F l b (n + 1) (n + 1 -i) = k b (θ i+1 -θ i ) for i = 1...n (11)
Considering the deformations of extension springs are very small, the total potential energy of all the bending springs of the system is

E mass-spring = 1 2 k b n i=1 (θ i+1 -θ i ) 2 (12)
With eqn. [START_REF] Cai | Modeling and dynamics simulation for deformable objects of orthotropic materials[END_REF] and eqn. [START_REF] Nealen | Physically based deformable models in computer graphics[END_REF] we obtain:

E mass-spring = F 2 l 2 b 2k b (n + 1) 2 n i=1 i 2 = F 2 l 2 b 2k b (n + 1) 2 n(n + 1)(2n + 1) 6 = F 2 l 2 b 12k b n(2n + 1) n + 1 (13)
By comparing eqn. [START_REF] Eischen | Finite element modeling and control of flexible fabric parts[END_REF] and eqn. ( 13), we can derive an analytical relation between k b and EI as following:

k b = EI l b n(2n + 1) 2(n + 1) (14)
Since eqn. ( 14) is derived based on the assumption of small displacement, we still have here a linear problem thus the principle of superposition can be applied. During the flight, the aerodynamic loads acting on insect wings can be considered to be equivalent to distributed loads on the surface of the wings. These distributed loads can be discretized into many point forces using a work-equivalent method [START_REF] Logan | A first course in the Finite Element Method[END_REF] and then the superposition principle can be applied. Thus, it is sufficient to analyze only one point force case to find the relation between EI and k b , since it does not depend on the point force F.

However, as mentioned at the beginning of this section, insect wings are deformed significantly to create lift for flying. Here, we are dealing with a large displacement problem and the question is if eqn. ( 14) still remains valid.

The technique used to derive ( 14) is no longer applicable since the solution for large deflection of a cantilever beam cannot be obtained analytically [START_REF] Barten | On the deflection of a cantilever beam[END_REF].

This problem involves calculating elliptical integrals of the second kind [START_REF] Bisshopp | Large deflections of cantilever beams[END_REF] and needs to be solved numerically. Consequently, the relation between EI and k b is put into a large displacement, nonlinear test case to check if we still get the same mechanical behaviors between the continuous beam and the mass-spring model. The results are presented in the next section. of the beam at equilibrium state can be calculated by using the fundamental Bernoulli-Euler theorem [START_REF] Wang | An explicit solution of the large deformation of a cantilever beam under point load at the free tip[END_REF] and the mass-spring network as given in table

1.The static state of the vein model (discretized by n = 64 mass points) is obtained by solving the dynamic equations of the system with artificial damping forces to make the system go quickly towards its balanced position.

Despite the small displacement assumption for deriving the relation between EI and k b , the relation in eqn. ( 14) is still valid even in very large deflection problem. For the case F = 3.92, the vertical displacement of the free end is already more than 30% of the total length of the beam and we still get very good agreements between both models with the relative error being smaller than 1%. The mapping from EI to k b can then be generalized for nonlinear, large deflection cases with good agreement between the continuum theory and the discrete mass-spring network. Dynamic case. The vein model will now be compared with another solid solver developed by Engels et al. [START_REF] Engels | Numerical simulation of vortex-induced drag of elastic swimmer models[END_REF]. It is based on the classical nonlinear beam equation, the Euler-Bernoulli theory. All details about this solver can be found in [START_REF] Engels | Numerical modeling of fluid-structure interaction in bioinspired propulsion[END_REF]. We study the case when we have a 2D cantilever beam (figure 7) of length l b = 1, density ρ = 0.0571, flexural rigidity EI = 0.0259.

The beam is in vacuum, subjected to a gravity field g = 0.7 strong enough to cause large deflection. All the parameters here are dimensionless. Both computations are performed for the same numerical parameters with the time step dt = 10 -2 and n = 64 grid points. Although both solvers require the same amount of CPU time for the same resolution, the mass-spring network is still more efficient since it is designed to deal with 3D problems. For the computation, the mass-spring solver handles a system of 3n degrees of freedom, corresponding to 3 dimensions, while the nonlinear-beam solver only solves for 2n variables.

The deflection line of the two models at a given time t and the oscillation of the trailing edge y te (t) are shown in figure 8. The dashed blue line calculated by the nonlinear beam theory and the solid red line calculated by the mass-spring network are almost coincident with each other. We have an excellent agreement between these two models with a relative error smaller than 1%.

Membrane model -Uniaxial and isotropic deformations of a two-dimensional sheet

We consider here the same test case proposed by Omori et al. [START_REF] Omori | Comparison between spring network models and continuum constitutive laws: Application to the large deformation of a capsule in shear flow[END_REF] where a square 2D sheet with an initial side length l 0 = 1 is extended by a uniaxial tension T = 0.005 and has a final length l in the x-direction at the equilibrium state, as shown in figure 4. This tension must be small enough to cause small deformation on the sheet. The Young modulus E s is defined by:

E s = T (15)
where is the strain.

The sheet is then discretized by using three types of meshes, illustrated in figure 4. The grid size ∆l is the side length of one triangle element of the mesh and inversely proportional to the number of grid points n. The grid size is varied to refine the mesh resolution. Since we are only interested in the equilibrium state of the model, the masses will not have any effect on the result and they are chosen properly for the numerical convergence. Instead of solving the static equation of the system, we still solve here the dynamic equations of the system with artificial damping forces to make the system go quickly towards its balanced position. Last but not least, all the spring stiffnesses are set to the same value k = 1. Figure 10 shows the results we get for all three mesh structures. First, for the cross-center structure, we are able to reproduce the result of Omori et al. [START_REF] Omori | Comparison between spring network models and continuum constitutive laws: Application to the large deformation of a capsule in shear flow[END_REF]. When the mesh is refined, the ratio k/E s converges to the analytical value 3/4 with the relative error being smaller than 1%. For the regular triangle, due to the shape of the square, we have some minor flaws of the mesh at the border. But these can be neglected when the mesh is fine enough and we can consider it as a regular triangular mesh. Indeed, for high resolution, we find again an excellent agreement with the analytical ratio k/E s = √ 3/2 derived by Omori et al. The relative error is also smaller than 1%. However, for the unstructured mesh, the convergent value of k/E s is larger than the one of Omori et al., but identical to the analytical solution for the regular triangle. This finding is in fact expected by Omori et al. since these two meshes are both constructed with triangles, each node being connected to six springs. Yet, the random structure of the unstructured mesh makes it difficult to explain the difference.

Using our mass-spring model, we are capable of reproducing the results from the references, which indicates the reliability of the solver. These results for both small and large deformation cases allow us to have the same conclusions as in the literature about the mass-spring system. Even though the mechanical properties of spring networks are strongly dependent on the mesh topology, a correlation between the discrete model and the continuum model can still be obtained if the mesh is fine enough. However, this needs to be compromised with the computational efficiency which is the main reason that we choose mass-spring network in the first place.

Wing structure

The simulation of insects with flexible wings is extremely complicated not only because it involves solving for both fluid and solid dynamics, but also due to the fact that insect wings are sophisticated structures. In our work, we want to take into account as much as possible all the mechanical properties of the bumblebee wing, in order to correctly model its dynamic behaviors. In our wing model there are three main factors introduced, which are considered to have the most impact on wing deformation during flight: venation pattern, mass distribution and vein stiffness.

Venation pattern

The venation architecture is claimed to affect the anisotropy of the wing.

Throughout measurements from different insects, Combes et al. [3,[START_REF] Combes | Flexural stiffness in insect wings II. Spatial distribution and dynamic wing bending[END_REF] suggest that wing flexural stiffness declines exponentially towards the tip and trailing edge. This is explained by the common venation patterns of insect wings: most insect wings have thick, stiff veins at the leading edge and cross veins are thinner as they expand toward the wing tip. This structure allows insect wings to resist against strong bending deformation in the spanwise direction, while creating camber for lift generation in the chordwise direction [START_REF] Combes | Flexural stiffness in insect wings II. Spatial distribution and dynamic wing bending[END_REF]. Nakata and Liu [2] modeled the anisotropy caused by wing veins. To this end they took into account the variation of wing thickness and introduced a "rule of mixture" of composite materials.

In our model, the functional approach is used to take into consideration the venation pattern. The vein structure, as well as the wing contour, are adapted from the data from [START_REF] Kolomenskiy | The dynamics of passive feathering rotation in hovering flight of bumblebees[END_REF] and encoded into the mass-spring network, as shown in figure 11. Comparing to the reference data, two more veins are added (vein 21 in the forewing and vein 7 in the hindwing) and two forewing veins 19 and 20 are extended toward the tip of the wing. These modifications are made to add bending stiffness to the tip of the wing and thus to obtain a more realistic behavior during the simulation. The meshing is done by identifying firstly the contour of the wing and all the veins (green, red and blue curves in figure 11). The membrane is then discretized by a triangular mesh using SALOME 2 , an open-source integration platform for numerical simulation and mesh generation.

Mass distribution

Another property which plays an essential role on the aerodynamics of the wing, is the mass distribution. It represents the inertia of the system and the position of the mass center has a strong connection with the stability of the wing during flight. The mass distribution is calculated based on the measured wing mass data from [START_REF] Kolomenskiy | The dynamics of passive feathering rotation in hovering flight of bumblebees[END_REF]. For our numerical simulations, the total wing mass is chosen as the same used by Kolomenskiy et al. [START_REF] Kolomenskiy | The dynamics of passive feathering rotation in hovering flight of bumblebees[END_REF], m w = 0.791 mg.

The mass is then distributed into vein and membrane parts based on their geometry and material. Each vein is considered as a rod composed of cuticle,

ρ c = 1300 kg/m 3 [42]
, with a circular cross section of constant diameter d v [START_REF] Kolomenskiy | The dynamics of passive feathering rotation in hovering flight of bumblebees[END_REF] and length l v , calculated directly from the model. The mass of each vein is then calculated and shown in table 2. Both diameter and mass are dimensionless quantities, normalized by wing length L and air density ρ air L 3 , respectively.

The mass distribution for the membrane is more tricky since we do not have the material properties of bumblebee membranes. A bi-linear regression is employed due to the fact that the membrane is heavier near the wing root and the leading edge [START_REF] Kolomenskiy | The dynamics of passive feathering rotation in hovering flight of bumblebees[END_REF]. An optimization is employed to find the parameters for the regression using the mass center from the measured data as an objective function. For a mass point m i belonging to the membrane at position [x i , y i ] (the z coordinate is neglected here because we assume that the membrane is a planar sheet), we get:

m i = 1.75 • 10 -4 -2.83 • 10 -4 x i + 4.91 • 10 -4 y i (16)
This yields a difference, between two mass centers, of 0.0013 in the x-direction and 0.0085 in the y-direction which are really small compared to the reference wing length R w = 1.

Vein stiffness estimation

To study the influence of wing flexibility on the aerodynamics performance, the flexural rigidity of veins will be changed to alter the bending stiffness of the wing. Consequently, only Young modulus E will be varied.

Insect cuticles are reported to have a Young modulus about 1kP a to 50M P a [START_REF] Vincent | Design and mechanical properties of insect cuticle[END_REF]. For our study, we are choosing two values of E (corresponding to flexible and highly flexible wings): 3.5 kP a and 35 kP a, which are in this range and give realistic deformations comparing to those observed in real life. Then, the flexural rigidity EI of each vein will be calculated using the second moment of inertia I of circular-section veins with diameters given in table 2. [START_REF] Kolomenskiy | The dynamics of passive feathering rotation in hovering flight of bumblebees[END_REF]) and their corresponding dimensionless mass m v .

Forewing

4. Fluid-structure interaction

Numerical method

The ultimate goal of this work is the fluid-structure interaction simulation of insects with flexible wings. To study the airflow as well as the effect of flexibility on the aerodynamic performance of the wing, the developed mass-spring model needs to be coupled with a fluid solver. This is done by combining the volume penalization method [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] with a Fourier pseudospectral discretization [START_REF] Kolomenskiy | A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles[END_REF], for which we developed the parallel opensource solver FLUSI, freely available on Github3 [START_REF] Engels | Flusi: A novel parallel simulation tool for flapping insect flight using a Fourier method with volume penalization[END_REF]. The code solves the incompressible, penalized Navier-Stokes equations

∂ t u + ω × u = -∇Π + ν∇ 2 u - χ C η (u -u s ) penalization term - 1 C sp ∇ × (χ sp ω) ∇ 2 sponge term (17)
∇ • u = 0 (18) u(x, t = 0) = u 0 (x) x ∈ Ω, t > 0 ( 19
)
where u is the fluid velocity, ω is the vorticity, Π = p + 1 2 u • u is the total pressure, ν is the kinematic viscosity. We find here again all the familiar terms of the classical Navier-Stokes equations except for the sponge and penalization terms. The former is a vorticity damping term used to gradually damp vortices and alleviate the periodicity inherent to the Fourier discretization. The last term is used to impose the no-slip boundary conditions on the fluid-solid interface as explained in [START_REF] Engels | Flusi: A novel parallel simulation tool for flapping insect flight using a Fourier method with volume penalization[END_REF]. All the geometry information of the solid is encoded in the mask function χ, which is usually taken as one inside the solid and zero otherwise. However, we are dealing with a moving flexible obstacle and the discontinuous mask function χ need to be replaced by a smooth one, eqn. [START_REF] Engels | Numerical modeling of fluid-structure interaction in bioinspired propulsion[END_REF], to avoid oscillations in the hydrodynamic forces [START_REF] Engels | Numerical modeling of fluid-structure interaction in bioinspired propulsion[END_REF].

Thus, we employ a mask function χ as shown below:

χ(δ) =              1 δ ≤ t w -h 1 2 1 + cos π (δ-tw+h) 2h t w -h < δ < t w + h 0 δ ≥ t w + h ( 20
)
where h is the semi-thickness of the smoothing layer, t w is the semi-thickness of the wing and δ is the distance function which gives us the distance from Eulerian fluid nodes to the center surface of the wing. As presented in section 3, an unstructured triangular mesh is employed for our wing model. Thus, the discretized wing surface is composed of triangles constructed by three vertices (e.g. x s,i , x s,j and x s,k ). The distance function δ is computed by cycling over all these triangles. Since we are only interested in the fluid nodes near the fluid-solid interface, a bounding box is used to save computing time. For each triangle, the distance from it to all the fluid nodes belonging to its bounding box will be computed by using the algorithm from [START_REF] Eberly | Distance between point and triangle in 3d, Geometric Tools[END_REF]. The distance function at one fluid node is finally the minimum distance from this fluid node to all the triangles nearby.

δ(x, t) = min(||x -triangle(x s,i , x s,j , x s,k )|| 2 ) (21)
The solid velocity field u s is calculated in the same way as the distance function δ. If the triangle (x s,i , x s,j , x s,k ) is the one closest to the fluid node

x, x will be projected onto the plane of the triangle and the solid velocity of the projected point is interpolated from the velocities of the three vertices by using barycentric interpolation. Because we do not consider the flexibility of the wing in the direction perpendicular to the wing surface, the velocity of the projected point should be the same as the one of the fluid node.

Nevertheless, the solid velocity field is defined in the global reference frame for the fluid solver while the velocity solved by the solid solver is in the local wing reference frame. These velocities are needed to be transformed back to the global reference frame using eqn. [START_REF] Combes | Flexural stiffness in insect wings II. Spatial distribution and dynamic wing bending[END_REF] where V O (w) and Ω are the translating and rotating velocity of the wing reference frame, v (w) and x (w) are the velocity and the position computed by the solid solver in the wing reference frame, respectively.

u s = V O (w) + v (w) + Ω × x (w) (22)
Moreover, the fluid also interacts with the wing via the pressure and viscous force. However, at Re = O(10 3 ) in our study, the viscous force is considered to be very small and only the pressure force is transferred into the mass-spring system as external force. Contrary to the previous calculation of the solid velocity field u s , the pressure force is interpolated from the Eulerian fluid grid onto the Lagrangian mass points. The pressure interpolation is quite straightforward because the pressure field solved by the penalized Navier-Stokes equations is following the Darcy law and continuous even inside the wing [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF]. Consequently, the pressure at any mass point can always be determined, using delta interpolation proposed by Yang et al. [START_REF] Yang | A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations[END_REF],

from pressure values at neighboring Eulerian grid points [START_REF] Engels | Numerical modeling of fluid-structure interaction in bioinspired propulsion[END_REF]. Then for each triangle element of the wing, the pressure forces at the three vertices are perpendicular to the triangle and have magnitudes equaling to the pressure multiplied by one third of the triangle area. This is done for all the triangles and then accumulated to obtain the overall external pressure forces acting on the mass-spring system.

For time-stepping, the coupled fluid-solid system is advanced by employing a semi-implicit staggered scheme as proposed in [START_REF] Engels | Numerical modeling of fluid-structure interaction in bioinspired propulsion[END_REF]. At time step t n , the fluid velocity field is advanced to new time level u n → u n+1 from the old mask function χ n and the old solid velocity field u n s by using the Adams-Bashforth scheme. Then, the pressure field at the new time step p n+1 is calculated from the fluid velocity field u n+1 . Finally, the solid is advanced to the new time step χ n+1 and u n+1 s and the whole process is repeated until the final time.

Validation

For the validation of the fluid-solid coupling, we consider two test cases: the Turek benchmark test case FSI3 [START_REF] Turek | Proposal for numerical benchmarking of fluidstructure interaction between an elastic object and laminar incompressible flow[END_REF][START_REF] Turek | Numerical benchmarking of fluid-structure interaction: A comparison of different discretization and solution approaches[END_REF] and the rigid revolving bumblebee wing test case [START_REF] Engels | Helical vortices generated by flapping wings of bumblebees[END_REF].

Turek benchmark FSI3

The Turek benchmark FSI3 involves a flexible appendage of length l = 0.35 and thickness h = 0.02 placed right behind a circle cylinder of radius R = 0.05; the whole obstacle is immersed in a channel of size H × L = 0.41 × 2.5 with a Poiseuille inflow of meanflow Ū = 2, as shown in figure 12.

The center of the cylinder is placed a bit lower to the centerline at (0.2, 0.2) to trigger the instability and to make the appendage oscillate.

The fluid solver, as well as the fluid-solid coupling, are handled by the software FLUSI [START_REF] Engels | Flusi: A novel parallel simulation tool for flapping insect flight using a Fourier method with volume penalization[END_REF] and the setup remains the same as the test case FSI3 with a resolution of 5200 × 1152 whose details can be found in [START_REF] Engels | Numerical modeling of fluid-structure interaction in bioinspired propulsion[END_REF]. Only the solid solver based on the nonlinear beam equation, which is used for validation in section 2, is now replaced by the new solver using the massspring network for validation. The results of this simulation are presented in table 3 for the comparison with the reference solutions presented in the literature [START_REF] Engels | Numerical modeling of fluid-structure interaction in bioinspired propulsion[END_REF][START_REF] Turek | Proposal for numerical benchmarking of fluidstructure interaction between an elastic object and laminar incompressible flow[END_REF][START_REF] Turek | Numerical benchmarking of fluid-structure interaction: A comparison of different discretization and solution approaches[END_REF].

For the oscillation of the trailing edge y te , the result is in excellent agreement with all three references when the maximum relative error, for both maximum and minimum values of y te , is only 1.76% and the relative error for the frequency of oscillation is 1.65%. The vertical displacement of the trailing edge with respect to time in the periodic state is also plotted in figure 13 to compare with the reference [START_REF] Turek | Numerical benchmarking of fluid-structure interaction: A comparison of different discretization and solution approaches[END_REF]. The two lines are almost superposed on each other. Nevertheless, the computed drag is less accurate with a relative error which can go up to 4.57% comparing the maximum value with the reference [START_REF] Engels | Numerical modeling of fluid-structure interaction in bioinspired propulsion[END_REF], but only 1% comparing with [START_REF] Turek | Proposal for numerical benchmarking of fluidstructure interaction between an elastic object and laminar incompressible flow[END_REF]. From figure 13, the curves of the two solutions appear to have the same shape but have some offset. This offset is explained in [START_REF] Engels | Numerical modeling of fluid-structure interaction in bioinspired propulsion[END_REF] to be due to the smoothing layer in the mask function, which plays a role as surface roughness. This leads to the over-prediction of the drag force. Concerning the lift force, the mass-spring model yields results very close to the one calculated by Engels with the error of 2.76%, and the difference is around 20% with respect to [START_REF] Turek | Proposal for numerical benchmarking of fluidstructure interaction between an elastic object and laminar incompressible flow[END_REF][START_REF] Turek | Numerical benchmarking of fluid-structure interaction: A comparison of different discretization and solution approaches[END_REF] for both max and min values. Like in Engels [START_REF] Engels | Numerical modeling of fluid-structure interaction in bioinspired propulsion[END_REF], the amplitude of the lift force is over-predicted by coupling FLUSI and the mass-spring solver. In conclusion we find satisfying agreement with the results from the literature, for the solid solver alone, as well as for the FSI algorithm coupling the solid solver with the fluid solver in 2D.

Rigid revolving wing

Prior to studying the flexibility of the wing, a common test case of a rigid revolving wing is considered to validate the coupling between the fluid and the solid solver in 3D, i.e. the mask function generation and the velocity field of the solid. The setup is taken the same as the one used by Engels et al. [START_REF] Engels | Helical vortices generated by flapping wings of bumblebees[END_REF] as shown in figure 14. The angle of attack is fixed at α = 45 • while the rotation angle φ(t) is given by

φ(t) = τ e -t/τ + t -τ (23)
The wing is rotated around the center of the computational domain of size 4 × 4 × 2, which is discretized using a mesh of 1024 × 1024 × 512 grid points. To be consistent with the reference simulation [START_REF] Engels | Helical vortices generated by flapping wings of bumblebees[END_REF], the wing shape is not the one presented in section 3 but adapted from the wing planform taken from the reference. The wing shape is then discretized by a triangular mesh as shown in figure 14b. However, the vein structure will not be taken

C L = F L M LT -2 ; C D = F D M LT -2 (24)
where the lift F L is the force in the vertical direction Oz and the drag F D is the force perpendicular to the plane formed by the vertical and the wing spanwise axes, as shown in figure 14a.

The computed lift and drag coefficients are shown in figure 15 along with the reference values from [START_REF] Engels | Helical vortices generated by flapping wings of bumblebees[END_REF]. To evaluate quantitatively the error, the average lift and drag during the steady state (for rotation angles φ varying from 160 • to 320 • ) are computed and compared with the reference. A very good agreement is obtained with the relative error of 1.3% for the drag and 1.6%

for the lift.

From the results obtained from these two test cases, the satisfactory agreements can give us the confidence about the solid solver, based on mass-spring system, as well as the coupling with the flow solver FLUSI. Any difference between all the numerical studies carried out can be explained by the difference between the continuum model and the discrete model together with the way of generating the mask function.

Results and discussion

In the following we present results of high resolution computations of revolving bumblebee wings which are either rigid, flexible or highly flexi- The mean drag generated during the second half cycle of the rotation is chosen for the evaluation of the mesh convergence (160

• ≤ φ ≤ 320 • ).
Because it is impossible to obtain the exact values for the mean drag in this case, we use here the result obtained with the finest mesh as a reference value. The relative error of the mean drag with respect to the reference drag for all the mesh size is shown in figure 17. In all the simulations, the penalization parameter C η is chosen to satisfy that the number of points per thickness of the penalization boundary layer K η = νC η /∆x is always constant (as recommended in [START_REF] Engels | Numerical modeling of fluid-structure interaction in bioinspired propulsion[END_REF]) and equal to 0.052. The drag obtained for each simulation (figure 16) shows the convergence to the finest resolution solution when we refine the mesh. The spatial convergence exhibits a first to second order behavior when we plot the error versus the mesh size.

Solid mesh

As mentioned above, the dynamics of the mass-spring system depends strongly on the mesh size. Thus another convergence test on the number of mass points is performed. Two simulations of a revolving flexible wing at resolution 768 2 × 384 are run to compare between a medium-mesh and a fine-mesh wing which are discretized by 465 and 1065 mass points, respectively. As shown here in figure 18, although the number of mass points is more than doubled, the forces remain almost unchanged with an increase of 1.1% and 0.8% in average lift and drag coefficients during the steady state, respectively. Since the fluid solver is itself already costly in term of CPU time, the medium-mesh wing with 465 mass points is sufficient and can be chosen for the following study in section 5.2.

Influence of wing flexibility

To examine the influence of vein stiffness on the aerodynamic performance of the wing, the flexural rigidity of veins will be varied by changing the The drag follows the same trend as the lift, but is larger in magnitude. When the flexibility of the wing is taken into account, the rapid rise at the beginning of the forces for both flexible and highly flexible wings disappear. Instead, the forces increase gradually and the more flexible the wing is, the lower the lift and the higher the drag are.

At steady state, similar behaviors between the rigid and the flexible wings can be observed. When the rotation angle reaches 160 • , the forces generated by these two wings are stabilized. This can be explained by the fact that no dynamic deformation of the wings takes place and just the shape plays a role.

We also find that the lift-to-drag ratio at the steady state of the flexible wing is 1.2, 14.5% higher than the one of the rigid case, which is only 1.05 (figure 20). This finding is consistent with conclusions found in literature [START_REF] Zhao | Aerodynamic effects of flexibility in flapping wings[END_REF][START_REF] Van De Meerendonk | Threedimensional flow and load characteristics of flexible revolving wings[END_REF]. A flexible wing generates less lift and drag than a rigid one. However, due to the flexibility of the wing, the bending in the chordwise direction makes the effective geometric angle of attack decrease and alters the direction of the total resultant force upward [START_REF] Van De Meerendonk | Threedimensional flow and load characteristics of flexible revolving wings[END_REF]. This makes the lift-to-drag ratio raise and allows better flight performance.

On the contrary, the highly flexible wing acts differently. Both the lift and the drag increase gradually to attain their maximum values at the rotation angle φ = 120 • and then decline instead of being stabilized as in the other simulations. The lift-to-drag ratio is surprisingly much less than the one of the rigid case at the beginning of the steady state but then increases and keeps up with the rigid wing. This can be explained by the fact that the bending of the wing in the spanwise direction (figure 21) prevents the development of the LEV growing further toward the wing tip and makes the LEV burst sooner at mid-span of the wing.

The change of aerodynamic forces compared to the rigid case is linked to the deformation of the wing, which is modeled by the mass-spring solver.

The wing deformation for all three cases is shown in the same figure 21 for comparison at three time instants t = 2, t = 4 and t = 6. By applying the functional approach, the difference between the vein and the membrane is visible in the visualization.

At the finest resolution of the mesh (1024 2 × 512), the flows generated by the flexible wing are shown (figure 22) by plotting their vorticity magnitude at four time instants of the simulation. The formations of the leading edge the rigid wing. This can be interpreted that there is an optimized zone of wing flexibility, which is ideal for flying.

For flapping wings we anticipate that flexibility will become important because of the dynamic wing deformation. In the near future we are planing high resolution numerical simulations of flapping insect flight with flexible wings where the dynamical deformation plays an important role.

The limitations of the current approach are the resolution and CPU time requirements imposed by the use of an uniform grid. Hence large scale simulations become prohibitively expensive. Moreover, the thickness of bumblebee wing is much smaller than the spatial mesh size in our present simulations.

Consequently, the virtual thickness of wings studied here is set to 4 times the mesh size, necessary for the usage of the volume penalization method.

An adaptive version of the FLUSI code, likewise fully parallel, using wavelet-based grid refinement is currently being developed to be able to reduce memory and CPU time requirements. High resolution numerical simulation of flapping flight for larger species and large Reynolds numbers will thus become possible.

Implementing the solid solver presented in the current paper into the adaptive Navier-Stokes solver will allow to perform fluid-structure interaction on adaptive grids at reduced computational cost.

were granted access to the HPC resources of IDRIS under the Allocation No.

  https://github.com/pseudospectators/FLUSI).

  used mass-spring networks to model insect wings in two-dimensional numerical simulations. A massspring network was used by Yeh and Alexeev [14] to model a flexible plate swimmer and performed fluid-structure interaction simulation with Lattice-Boltzmann methods in 3D. The development of our solid solver is motivated by their mass-spring network approach, aiming to model the flexibility of insect wings in the three-dimensional case. The goal of this paper is to move from rigid to flexible wings and to present a fluid-structure interaction solver for flapping flight, based on the open access software FLUSI, where we integrated a solid model based on mass-spring systems. However, the wing kinematics of insects is very difficult to obtain, since the measurements usually require high tech equipment

  wing models: hawkmoth (Manduca sexta), bumblebee (Bombus ignitus) and fruitfly (Drosophila melanogaster). Van de Meerendonk et al.[START_REF] Van De Meerendonk | Threedimensional flow and load characteristics of flexible revolving wings[END_REF] investigated experimentally the flow field and fluid-dynamic loads of a flexible revolving wing to quantify the influence of flexibility on the force generation performance of the wing. In our study, we also consider revolving flexible bumblebee wings and assess the influence of the wing deformation.The outline of this paper is the following: In sec. 2 we present a massspring model for describing the flexible insect wings. The wing structure of the considered bumblebee wings and its mass distribution are detailed in sec. 3. The numerical artillery of fluid-structure interaction is explained in sec. 4 and some validation tests are given. The numerical results as well as the discussion about the influence of the flexibility on the aerodynamic performance are presented in sec. 5. Finally, conclusions are drawn in sec. 6,

  Let us now study a mass-spring network discretized into n+2 mass points connected by bending and extension springs as shown in figure 6. All the bending and extension springs are the same with a stiffness k b and k e , respectively and k e k b . The first two points are totally fixed to represent the boundary condition of the fixed end of the beam. Writing the equilibrium equations

2. 4 .

 4 Validations of the mass-spring model 2.4.1. Vein model -Cantilever beam under gravitational force Static case. Firstly, we consider a static case of a cantilever beam with length L = 0.3, flexural rigidity EI = 0.24 and loaded by a point force F at the free end. The force F varies from 0.39 to 11.76 and it must be strong enough to cause a large deflection. All the parameters here are dimensionless. The vertical displacement δy and the horizontal displacement δx of the free end

  into account in this model because we are considering a rigid wing. The triangular mesh is solely exploited for the creation of the mask function and the solid velocity field by using the algorithm presented at the beginning of this section. Here, all the quantities are normalized. The wing length is chosen as the length scale L = 13.2 mm; the mass scale is based on the air density M = ρ air ×L 3 = 2.817 mg and the time scale is chosen in the way that wing tip velocity is unity, thus T = 1 s . The Reynolds number is then defined as in[START_REF] Engels | Helical vortices generated by flapping wings of bumblebees[END_REF] Re = ūtip c m /ν where the mean wingtip velocity ūtip = 1 [LT -1 ] by definition from eqn. 23, the fluid viscosity ν = 1.477 • 10 -4 [L 2 T -1 ] and the mean chord, the ratio between the wing surface area A and the wing length R w , c m = A/R w = 0.304 [L]. This yields Re = 2060. Additionally, the lift and drag coefficient are defined as below

5 . 1 .

 51 ble. First we perform computations for different resolutions to check the mesh convergence for both fluid and solid solvers. Then a comparison of the flexible wings with the rigid case allows to assess the influence of the wing deformation on the aerodynamic forces. The setup is exactly the same as described in the revolving wing test case of the previous section. The only difference is the wing shape which is now changed back to the one presented in section 3. As a result, the length scale and the mass scale are changed as follows: L = 15 mm and M = ρ air × L 3 = 4.13 mg while the time scale remains the same T = 1 s. The corresponding Reynolds number is 1800 where the fluid viscosity is assumed to be ν = 1.477 • 10 -4 [L 2 T -1 ], the wing tip velocity u tip = 1 [LT -1 ] and the mean chord calculated from the new wing surface area is c m = A/R w = 0.266 [L]. Study of mesh convergence 5.1.1. Fluid mesh The following mesh convergence study for the fluid solver is performed considering five different resolutions: 128 × 128 × 64, 256 × 256 × 128, 512 × 512 × 256, 768 × 768 × 384 and 1024 × 1024 × 512.

  Young modulus E. Two values of the Young modulus are used:E = 1.25 • 10 8 [M L -1 T -2 ] and E = 1.25 • 10 7 [M L -1 T -2 ],corresponding to the flexible and highly flexible cases, respectively. Lift and drag coefficients at resolution 1024 2 × 512 for the rigid, flexible and highly flexible cases are presented in figure 19. During the transition phase (rotation angle φ ≤ 40 • ), the lift generated by the rigid wing increases instantly and then decreases before going up again.

  2018-91664 attributed by GENCI (Grand Équipement National de Calcul Intensif). For this work, we were also granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (No. ANR-10-EQPX-29-01). The authors thankfully acknowledge financial support granted by the ministères des Affaires étrangères et du développement international (MAEDI) et de l'Education nationale et l'enseignement supérieur, de la recherche et de l'innovation (MENESRI), and the Deutscher Akademischer Austauschdienst (DAAD) within the French-German Procope project FIFIT. D.K. gratefully acknowledges financial support from the JSPS KAK-ENHI Grant No. JP18K13693. K.S. thanks the organizers of ICCFD10, in particular C.H. Bruneau, for inviting him for a plenary talk.

Figure 1 :

 1 Figure 1: Illustration of the restoring forces corresponding to the deformation of extension and bending springs.

Figure 2 :

 2 Figure 2: All possible shapes of a bending spring corresponding to one bending angle θ ijk .

Figure 3 :

 3 Figure 3: Illustration of a vein modeled by mass points, extension springs and bending springs. White circles represent mass points, solid lines represent extension springs, black circles represent both mass points and bending springs.

Figure 4 :

 4 Figure 4: Different mesh structures for 2D discretization [31].

Figure 5 :

 5 Figure 5: Relationship between the spring constant k, Young modulus E s , and Poisson ratio ν s in the small deformation limit for a 2D membrane under uniaxial deformation. The membrane is discretized by three types of meshes and subjected to homogeneous deformation. Figure adapted from [31].

Figure 9 :

 9 Figure 9: Deformation of a 2D sheet along the x-axis under the uniaxial tension T .

Figure 10 :Figure 11 :

 1011 Figure 10: Effect of mesh topology on the relation between the spring stiffness k and the corresponding Young modulus E s

Figure 12 :Figure 13 :

 1213 Figure 12: Computational domain of the FSI3 Turek benchmark and dimensions of the solid part [49].

Figure 22 :

 22 Figure 22: Flows generated by a flexible revolving wing, visualized by their vorticity magnitude |ω| at four time instants t = 1, 2, 4 and 6. The simulation is performed with resolution 1024 2 × 512.

Table 1 :

 1 Comparison between the continuum theory and the discrete mass-spring network in the static large deflection case. ] δy ref[10 -2 ] δx[10 -2 ] δy[10 -2 ] err x [%] err y[%]

	Point force	Nonlinear beam [39]	Mass-spring network	Relative error
	F δx ref [10 -2 0.39 29.96	-1.46	29.96	-1.46	0	0
	1.96	29.02	-6.93	29.01	-6.92	0.03	0.14
	3.92	26.87	-12.14	26.85	-12.1	0.07	0.33
	7.84	22.53	-17.93	22.53	-17.85	0	0.45
	11.76	19.37	-20.69	19.40	-20.57	0.15	0.58

Table 2 :

 2 Dimensionless vein diameter d v (adapted from

	Hindwing

Table 3 :

 3 Results of the FSI3 benchmark.

https://www.salome-platform.org/

https://github.com/pseudospectators/FLUSI

Acknowledgments

Financial support from the Agence Nationale de la Recher-che (ANR Grant No. 15-CE40-0019) and Deutsche Forschungsgemeinschaft (DFG Grant No. SE 824/26-1), project AIFIT, is gratefully acknowledged. The authors

vortex as well as the tip vortex can be observed clearly at the beginning of the rotation (t = 1.0 and t = 2.0). Then, the vortex burst happens and a region of inhomogeneous vorticity forms at the wing tip. However, the LEV remains attached to the wing surface and this results in constant lift and drag.

Conclusions

We presented a numerical approach for fluid-structure interaction in the open source framework FLUSI, which is based on a mass-spring model describing the structure of the insect wings and a pseudospectral method for solving the incompressible Navier-Stokes equations. For imposing no-slip boundary conditions in the complex time-changing geometry we used the volume penalization technique. The solver has been implemented on massively parallel supercomputers using MPI and allows high resolution computations, here with more than half a billion grid points. Code validation for two classical benchmarks, a flow past a cylinder with a flexible appendage and the flow generated by a rigid revolving wing, is likewise presented.

Considering the flexible wing, the flexibility reduces the buildup of the aerodynamic force during the beginning of motion. Nevertheless, after the start-up phase, the wing yields a steady state configuration, and no significant oscillation nor unsteady deformation of the wing are observed. A better aerodynamic performance of the flexible wing, characterized by the increase of the lift-to-drag ratio during the steady state, is explained by the decrease of the effective angle of attack caused by the deformation of the flexible wing.

On the other hand, the highly flexible wing appears to be less efficient than