Access the full text.
Sign up today, get DeepDyve free for 14 days.
Loading next page...
/lp/arxiv-cornell-university/a-mass-spring-fluid-structure-interaction-solver-application-to-1pAcZsC9v5
- ISSN
- 0045-7930
- eISSN
- ARCH-3341
- DOI
- 10.1016/j.compfluid.2020.104426
- Publisher site
- See Article on Publisher Site
Abstract
The secret to the spectacular ight capabilities of apping insects lies in their wings, which are often approximated as at, rigid plates. Real wings are how- ever delicate structures, composed of veins and membranes, and can undergo signi cant deformation. In the present work, we present detailed numerical simulations of such deformable wings. Our results are obtained with a uid{ structure interaction solver, coupling a mass{spring model for the exible wing with a pseudo-spectral code solving the incompressible Navier{Stokes equations. We impose the no-slip boundary condition through the volume pe- nalization method; the time-dependent complex geometry is then completely described by a mask function. This allows solving the governing equations of the uid on a regular Cartesian grid. Our implementation for massively parallel computers allows us to perform high resolution computations with Email addresses: dinh-hung.truong@univ-amu.fr (Hung Truong), thomas.engels@ens.fr (Thomas Engels), dkolomenskiy@jamstec.go.jp (Dmitry Kolomenskiy), kai.schneider@univ-amu.fr (Kai Schneider) Preprint submitted to Computers & Fluids February 21, 2020 arXiv:2002.08642v1 [physics.flu-dyn] 20 Feb 2020 up to 500 million grid points. The mass{spring model uses a functional ap- proach, thus modeling the di erent mechanical behaviors of the veins and the membranes of the wing. We perform a series of numerical simulations of a exible revolving bumblebee wing at a Reynolds number Re = 1800. In order to assess the in uence of wing exibility on the aerodynamics, we vary the elasticity parameters and study rigid, exible and highly exible wing models. Code validation is carried out by computing classical benchmarks. Keywords: insect ight, wing elasticity, mass{spring model, uid{structure interaction, spectral method, volume penalization method 1. Introduction A fundamental characteristics of insect ight are exible wings, which play an important role for their aerodynamics [1, 2, 3], requiring lower forces than their rigid counter parts and producing reduced sound. Numerical sim- ulation of insect ight is itself a sophisticated task, because it involves the solution of uid-solid interaction problems. Thus, we have to model simul- taneously the uid part and the solid part by using two coupled solvers. The uid solver [4] we are using in the present work has been developed previously and is called FLUSI , a fully parallel software dedicated for mod- eling three-dimensional apping ight in viscous ows. 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. FLUSI: freely available for noncommercial use from GitHub (https://github.com/pseudospectators/FLUSI). 2 Moreover, the volume penalization method is used to take into account the no-slip boundary conditions on the interfaces between the uid and the solid part. In [5], 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 apping bumblebee wings were presented in [6] using likewise the FLUSI code with rigid wings focusing again on the role of LEVs and the associated helicity production. The interaction of apping bumblebees with turbulent in ow in free and tethered ight was studied in [7] using once again massively parallel com- putations with FLUSI. We found that the uctuations of aerodynamic ob- servables signi cantly grow with increasing turbulence intensity, even if the mean values are almost una ected. Changing the length scale of the tur- bulent in ow, while keeping the turbulence intensity xed, showed that the uctuation level of forces and moments can be signi cantly reduced. We also found that the scale-dependent energy distribution in the surrounding tur- bulent ow is a relevant factor conditioning how ying 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 ight 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 ight performance. Consequently, a wing model is required for simulating 3 the deformation of the solid part under the impact from the uid. 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 elds due to their reliability and e ectiveness. Despite this dominance, the requirement for faster and more ecient methods motivates the development of alternative models. One of these is the mass-spring system, which is known for its computational eciency and straightforward implementation [8]. As part of this work, a solver using a network of masses and springs is developed to model a exible insect wing and coupled with FLUSI. th The story of mass-spring models started back at the end of the 20 century when people were dealing with observable deformations of exible objects in computer graphics applications, such as soft tissue, skin, hair, ball, cloth, textiles, etc. Being considered as one of the pioneers on this eld, Ter- zopoulos et al. [9] 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 exible objects had been calculated using nite element meth- ods until the requirement of faster models was claimed by Eischen et al. [10]. Consequently, the development of physically based deformable models started to grow, especially in the eld of computer graphics [11] and biomedical engi- neering [8]. While classic solvers, based on nite elements or nite di erence 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 4 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 pro- posed, mass-spring networks stand out as the most intuitive and simplest due to their computational eciency [12]. Mass-spring systems have already been employed in many elds such as medical applications [8] (muscle, red blood cells and virtual surgery), computer graphics, uid-structure inter- action and insect ight. Miller and Peskin [13] used mass-spring networks to model insect wings in two-dimensional numerical simulations. A mass- spring network was used by Yeh and Alexeev [14] to model a exible plate swimmer and performed uid-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 exibility of insect wings in the three-dimensional case. The goal of this paper is to move from rigid to exible wings and to present a uid-structure interaction solver for apping ight, 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 di- cult to obtain, since the measurements usually require high tech equipment to capture all the dynamic motion at small time scale and length scale. In- stead, a revolving wing model is usually employed to study the aerodynamics of apping wings thanks to its simplicity. Accordingly, the ow elds and force generation aspects of revolving wings have been analyzed for a wide range of parameters, as reported in the literature [15, 16, 17]. Di et al. [18] studied the role of forewings in generating LEVs of three revolving insect 5 wing models: hawkmoth (Manduca sexta), bumblebee (Bombus ignitus) and fruit y (Drosophila melanogaster). Van de Meerendonk et al. [19] inves- tigated experimentally the ow eld and uid-dynamic loads of a exible revolving wing to quantify the in uence of exibility on the force generation performance of the wing. In our study, we also consider revolving exible bumblebee wings and assess the in uence of the wing deformation. The outline of this paper is the following: In sec. 2 we present a mass- spring model for describing the exible insect wings. The wing structure of the considered bumblebee wings and its mass distribution are detailed in sec. 3. The numerical artillery of uid-structure interaction is explained in sec. 4 and some validation tests are given. The numerical results as well as the discussion about the in uence of the exibility on the aerodynamic performance are presented in sec. 5. Finally, conclusions are drawn in sec. 6, including some perspectives for future investigations. 2. Flexible wing 2.1. Mass-spring model The very basic idea of the mass-spring model is to discretize an object using mass points m (i = 1; :::; n) connected by springs. There exist many kinds of springs for di erent 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 de ned by the position x and the velocity v of the mass point i. For this, we need i i 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 6 of the mass-spring network where these governing equations are simply the corresponding classical well-known Newton's second law, given in eqns. (1). F = m a i i i int ext F = F + F for i = 1:::n i i (1) v (t = 0) = v i 0;i x (t = 0) = x i 0;i where n is the number of mass points, F is the total force (internal force int ext th F and external force F ) acting on the i mass point, m , a are mass i i i i th and acceleration of the i mass point, respectively. Here, all terms are quite simple to derive except for the forces. The exter- nal forces come from outside of the system and depend on the surrounding eld 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 di erential 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 analyti- cal solution cannot be derived. In practice, it is more convenient to convert a system of 3n second order equations into a system of 6n rst order equations by using the relations a = dv =dt and v = dx =dt. Eqns. (1) can then be i i i i 7 rewritten as below: dx = v dt dv int ext m = F + F for i = 1:::n i i dt (2) v (t = 0) = v i 0;i x (t = 0) = x i 0;i Let us call q = x ; v the phase vector containing positions and velocities i i 1 int ext of all mass points and f (q) = v ; m (F + F ) the right hand side i i i function, eqns.(2) can be rewritten again under the familiar form of a system of rst order ODEs as follows: dq = f (q; t) dt (3) q(t = 0) = q 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 sti 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 di usion, 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 [20]. Furthermore, the coupling between the uid and the solid solver will require an adaptive time stepping scheme. For all these reasons, a second 8 order backward di erentiation scheme with variable time steps [21] is used to discretize eqns. (3) in time. 2 2 (1 + ) 1 + n+1 n n 1 n n+1 q q + q = t f (q ) (4) i i i i 1 + 2 1 + 2 1 + 2 n n 1 n where = t =t is the ratio between the current t and the pre- n 1 vious time step t . Eqn. (4) is a system of nonlinear equations with n+1 the variable q , the phase vector of the system at the current time step, which needs to be solved. The Newton{Raphson method, a powerful itera- tive 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 nd approximations of the root with the rate of convergence estimated to be quadratic. For our mass-spring solver, the ini- tial guess chosen is the phase vector q solved at the previous time step; this allows the Newton{Raphson method to converge quickly since the structure n+1 is advanced slowly and smoothly in time, which assures that q remains close to q . In most cases, the Newton{Raphson method in the solver needs three to four iterations to converge within a relative or absolute L norm error of 10 as the stopping criterion. 2.2. Extension springs and bending springs Extension springs ( gure 1a) and bending springs ( gure 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: 9 Linear extension spring F = k (jjx jj jjx jj) e i ij 0;ij ij ij (5) F = F j i where k is the extension sti ness, e = (x x )=jjx x jj is the unit ij j i j i ij position vector and F and F are the restoring forces of the extension i j spring acting on two points i and j, respectively; Linear bending spring M = k ( ) (6) ijk ijk 0;ijk ijk or in terms of forces F = k ( ) (e e ) e i ijk 0;ijk ij jk ij ijk (7) F = k ( ) (e e ) e k ijk 0;ijk ij jk jk ijk F = F F j i k where M is the restoring moment, k is the bending sti ness, is ijk 0;ijk ijk the initial bending angle among three points i; j and k, the current ijk bending angle and F , F , F are the restoring force vectors (as shown i j k in gure 1) of the bending spring acting on three points i; j and k, respectively. However, the calculation of in eqn.(7) is not trivial since it involves ijk the geometrical de nition of the angle in 3D space with respect to a xed 10 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 = T T T (x ; y ) , x = (x ; y ) and x = (x ; y ) . The angle is now determined by: i i j j j k k k = atan2 (y y )(x x ) (y y )(x x ); ijk k j j i j i k j (8) (x x )(x x ) + (y y )(y y ) k j j i k j j i Here, atan2 (usually known as two-argument arctangent) is a special func- tion rst 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 sucient. 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 . Nev- ijk ertheless, corresponding to this , there is an in nite number of solutions ijk x , x and x that can satisfy this deformation and the set of all these so- i j k lutions forms a cone, like shown in gure 2. Consequently, one more angle is needed to obtain a unique solution. To de ne 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 and on the ijk ijk Oxy and the Oxz planes, respectively. Then, like in the 2D problem, these two bending angles are calculated as below: 11 h = atan2 (x x )(y y ) (x x )(y y ); j i k j k j j i ijk (x x )(x x ) + (y y )(y y ) j i k j j i k j (9) = atan2 (x x )(z z ) (x x )(z z ); j i k j k j j i ijk (x x )(x x ) + (z z )(z z ) j i k j j i k j 2.3. 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 elastic- ity properties are still poorly understood. Certain studies have shown that the vein arrangements in insect wings have strong impact on their mechan- ical properties [3, 22]. Thus, it will be inaccurate to consider a wing as a homogeneous structure; the vein pattern as well as the di erence 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 speci c 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 nd a thin layer called mem- brane. 12 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 e ect of all the external loads applied on a mechanical structure results in deformation which can usually be classi ed into three main types: bending, twisting and stretching. Although the whole wing is observed to be twisted signi cantly in many studies using high speed cameras or the digital particle image velocimetry [23, 24, 25, 26, 27], 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 gure 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 ) 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 ). i 1 When apping, most of external forces will act in the direction perpen- dicular to the wing surface. As a result, the stretching deformation is negli- 13 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 sti ness of these extension springs is arti cial and they do not need to re ect the mechanical property of the vein itself. They should be chosen sti enough to make the rod unstretched but not too sti to avoid problems with numerical stability when integrating the dynamical system in time. Compared to veins, membranes are totally di erent in terms of geomet- rical 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 [28, 29]. 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 gure 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 [30] since the random orientations of the springs will average out the forces. 2.3.1. Correlation between mass-spring network models and continuum con- stitutive 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 14 object is made. There are two main parameters needed to be assigned for a mass-spring model: the masses and the spring sti ness. Although a Voronoi diagram can be used to nd the masses in a proper way [32], selecting spring sti ness is still an open question and there are two common solutions to overcome this [33]. The rst approach is based on optimization methods to minimize the di erence between the results solved by the mass-spring model and the reference data. These reference data can come from the measure- ments, the visual appearance of real objects [34] or numerical solutions using validated methods such as nite element methods [35, 36]. 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 sti ness by using optimization can give satisfying results with reasonable computational cost. The second way is about deriving a relation between spring sti ness and other continuum mechanic properties, such as Young modulus, the Poisson ratio and the exural rigidity. In contrast to the discrete models, the elas- ticity parameters have been obtained for many materials and can be used to calculate the corresponding spring sti ness. Omori et al. [31] 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 gure 5. For the vein model, a relation between the exural rigidity EI and the sti ness of the bending springs k is needed. To derive this relation, we consider a classical problem of a cantilever beam length l , under a point 15 force F at the free end ( gure 6). In the limit of small displacement, the principle of energy yields the value of the bending spring sti ness k as a function of the exural 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. (10). 2 3 F l E = (10) beam 6EI 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. Let us now study a mass-spring network discretized into n+2 mass points connected by bending and extension springs as shown in gure 6. All the bending and b e extension springs are the same with a sti ness k and k , respectively and e b k k . The rst two points are totally xed to represent the boundary condition of the xed end of the beam. Writing the equilibrium equations for the remaining n points, we have: F (n + 1 i) = k ( ) for i = 1:::n (11) i+1 i (n + 1) Considering the deformations of extension springs are very small, the total potential energy of all the bending springs of the system is b 2 E = k ( ) (12) mass spring i+1 i i=1 With eqn. (11) and eqn. (12) we obtain: 16 n 2 2 F l b 2 E = i mass spring b 2 2k (n + 1) i=1 2 2 F l n(n + 1)(2n + 1) (13) b 2 2k (n + 1) 6 2 2 F l n(2n + 1) 12k n + 1 By comparing eqn. (10) and eqn. (13), we can derive an analytical relation between k and EI as following: EI n(2n + 1) k = (14) l 2(n + 1) 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 ight, 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 [37] and then the superposition principle can be applied. Thus, it is sucient to analyze only one point force case to nd the relation between EI and k , since it does not depend on the point force F. However, as mentioned at the beginning of this section, insect wings are deformed signi cantly to create lift for ying. 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 de ection of a cantilever beam cannot be obtained analytically [38]. This problem involves calculating elliptical integrals of the second kind [39] 17 and needs to be solved numerically. Consequently, the relation between EI and k 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. 2.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, exural 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 de ection. All the parameters here are dimensionless. The vertical displacement y and the horizontal displacement x of the free end of the beam at equilibrium state can be calculated by using the fundamental Bernoulli{Euler theorem [40] 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 arti cial 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 , the relation in eqn. (14) is still valid even in very large de ection 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 can then be generalized for nonlinear, large de ection cases with good agreement between the continuum theory and the discrete mass-spring network. 18 Table 1: Comparison between the continuum theory and the discrete mass-spring network in the static large de ection case. Point force Nonlinear beam [39] Mass-spring network Relative error 2 2 2 2 F x [10 ] y [10 ] x[10 ] y[10 ] err [%] err [%] ref ref x y 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 Dynamic case. The vein model will now be compared with another solid solver developed by Engels et al. [41]. It is based on the classical nonlinear beam equation, the Euler{Bernoulli theory. All details about this solver can be found in [20]. We study the case when we have a 2D cantilever beam ( gure 7) of length l = 1, density = 0:0571, exural rigidity EI = 0:0259. The beam is in vacuum, subjected to a gravity eld g = 0:7 strong enough to cause large de ection. All the parameters here are dimensionless. Both computations are performed for the same numerical parameters with the time step dt = 10 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 ecient since it is designed to deal with 3D problems. For the computation, the mass-spring solver handles a system of 3n degrees of 19 freedom, corresponding to 3 dimensions, while the nonlinear-beam solver only solves for 2n variables. The de ection line of the two models at a given time t and the oscilla- tion of the trailing edge y (t) are shown in gure 8. The dashed blue line te 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%. 2.4.2. Membrane model - Uniaxial and isotropic deformations of a two-dimensional sheet We consider here the same test case proposed by Omori et al. [31] where a square 2D sheet with an initial side length l = 1 is extended by a uniaxial tension T = 0:005 and has a nal length l in the x-direction at the equilibrium state, as shown in gure 4. This tension must be small enough to cause small deformation on the sheet. The Young modulus E is de ned by: E = (15) where is the strain. The sheet is then discretized by using three types of meshes, illustrated in gure 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 re ne the mesh resolution. Since we are only interested in the equilibrium state of the model, the masses will not have any e ect 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 20 equations of the system with arti cial damping forces to make the system go quickly towards its balanced position. Last but not least, all the spring sti nesses 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. [31]. When the mesh is re ned, the ratio k=E 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 aws of the mesh at the border. But these can be neglected when the mesh is ne enough and we can consider it as a regular triangular mesh. Indeed, for high resolution, we nd again an excellent agreement with the analytical ratio k=E = 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 is larger than the one of Omori et al., but identical to the analytical solution for the regular triangle. This nding 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 dicult to explain the di erence. Using our mass-spring model, we are capable of reproducing the results from the references, which indicates the reliability of the solver. These re- sults 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 ne enough. However, this needs to be compromised with the computational eciency which is the main reason 21 that we choose mass-spring network in the rst place. 3. Wing structure The simulation of insects with exible wings is extremely complicated not only because it involves solving for both uid 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 ight: venation pattern, mass distribution and vein sti ness. 3.1. Venation pattern The venation architecture is claimed to a ect the anisotropy of the wing. Throughout measurements from di erent insects, Combes et al. [3, 22] sug- gest that wing exural sti ness declines exponentially towards the tip and trailing edge. This is explained by the common venation patterns of insect wings: most insect wings have thick, sti 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 [22]. Nakata and Liu [2] modeled the anisotropy caused by wing veins. To this end they took into account the variation of wing thickness and intro- duced 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 22 adapted from the data from [42] and encoded into the mass-spring network, as shown in gure 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 modi cations are made to add bending sti ness to the tip of the wing and thus to obtain a more realistic behavior during the simulation. The meshing is done by identifying rstly the contour of the wing and all the veins (green, red and blue curves in gure 11). The membrane is then discretized by a triangular mesh using SALOME , an open-source integration platform for numerical simulation and mesh generation. 3.2. 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 ight. The mass distribution is calculated based on the measured wing mass data from [42]. For our numerical simulations, the total wing mass is chosen as the same used by Kolomenskiy et al. [42], m = 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, = 1300 kg=m [42], with a circular cross section of constant diameter d c v [42] and length l , 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 L , air https://www.salome-platform.org/ 23 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 [42]. An optimization is employed to nd the parameters for the regression using the mass center from the measured data as an objective function. For a mass point m belonging to the membrane at position [x ; y ] (the z coordinate is neglected here because we assume that i i the membrane is a planar sheet), we get: 4 4 4 m = 1:75 10 2:83 10 x + 4:91 10 y (16) i i i This yields a di erence, 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 = 1. 3.3. Vein sti ness estimation To study the in uence of wing exibility on the aerodynamics perfor- mance, the exural rigidity of veins will be changed to alter the bending sti ness of the wing. Consequently, only Young modulus E will be varied. Insect cuticles are reported to have a Young modulus about 1kPa to 50MPa [43]. For our study, we are choosing two values of E (corresponding to exible and highly exible wings): 3:5 kPa and 35 kPa, which are in this range and give realistic deformations comparing to those observed in real life. Then, the exural 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. 24 Forewing Hindwing Nominal Nominal Nominal Nominal # # diameter mass diameter mass 1 0.0070 0.0209 1 0.0065 0.0180 2 0.0074 0.0237 2 0.0043 0.0071 3 0.0055 0.0076 3 0.0046 0.0024 4 0.0070 0.0063 4 0.0011 0.0001 5 0.0040 0.0031 5 0.0038 0.0043 6 0.0048 0.0094 6 0.0037 0.0005 7 0.0040 0.0019 7 0.0020 0.0012 8 0.0038 0.0009 9 0.0041 0.0023 10 0.0048 0.0064 11 0.0045 0.0017 12 0.0038 0.0018 13 0.0042 0.0010 14 0.0038 0.0020 15 0.0034 0.0008 16 0.0032 0.0005 17 0.0032 0.0004 18 0.0044 0.0009 19 0.0015 0.0001 20 0.0018 0.0001 21 0.0020 0.0009 Table 2: Dimensionless vein diameter d (adapted from [42]) and their corresponding dimensionless mass m . 25 4. Fluid-structure interaction 4.1. Numerical method The ultimate goal of this work is the uid-structure interaction simulation of insects with exible wings. To study the air ow as well as the e ect of exibility on the aerodynamic performance of the wing, the developed mass-spring model needs to be coupled with a uid solver. This is done by combining the volume penalization method [44] with a Fourier pseudospectral discretization [45], for which we developed the parallel opensource solver FLUSI, freely available on Github [4]. The code solves the incompressible, penalized Navier-Stokes equations 1 ( !) sp @ u + ! u = r + r u (u u ) r (17) t s C C r sp | {z } | {z } penalization term sponge term r u = 0 (18) u(x; t = 0) = u (x) x 2 ; t > 0 (19) where u is the uid velocity, ! is the vorticity, = p + u u is the total pressure, is the kinematic viscosity. We nd here again all the familiar terms of the classical Navier{Stokes equations except for the sponge and pe- nalization terms. The former is a vorticity damping term used to gradually damp vortices and alleviate the periodicity inherent to the Fourier discretiza- tion. The last term is used to impose the no-slip boundary conditions on the uid-solid interface as explained in [4]. All the geometry information of the https://github.com/pseudospectators/FLUSI 26 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 exible obstacle and the discontinuous mask function need to be replaced by a smooth one, eqn. (20), to avoid oscillations in the hydrodynamic forces [20]. Thus, we employ a mask function as shown below: >1 t h ( t +h) 1 w () = (20) 1 + cos t h < < t + h w w 2 2h 0 t + h where h is the semi-thickness of the smoothing layer, t is the semi-thickness of the wing and is the distance function which gives us the distance from Eulerian uid 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 ; x and x ). The distance function is computed by s;i s;j s;k cycling over all these triangles. Since we are only interested in the uid nodes near the uid-solid interface, a bounding box is used to save computing time. For each triangle, the distance from it to all the uid nodes belonging to its bounding box will be computed by using the algorithm from [46]. The distance function at one uid node is nally the minimum distance from this uid node to all the triangles nearby. (x; t) = min(jjx triangle(x ; x ; x )jj ) (21) s;i s;j s;k 2 The solid velocity eld u is calculated in the same way as the distance function . If the triangle (x ; x ; x ) is the one closest to the uid node s;i s;j s;k 27 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 exibility 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 uid node. Nevertheless, the solid velocity eld is de ned in the global reference frame for the uid 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. (22) where V (w) and are the (w) (w) translating and rotating velocity of the wing reference frame, v and x are the velocity and the position computed by the solid solver in the wing reference frame, respectively. (w) (w) u = V (w) + v + x (22) Moreover, the uid also interacts with the wing via the pressure and viscous force. However, at Re = O(10 ) 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 eld u , the pressure force is interpolated from the Eulerian uid grid onto the Lagrangian mass points. The pressure interpolation is quite straightforward because the pressure eld solved by the penalized Navier{Stokes equations is following the Darcy law and continuous even inside the wing [44]. Consequently, the pressure at any mass point can always be determined, using delta interpolation proposed by Yang et al. [47], from pressure values at neighboring Eulerian grid points [20]. Then for each 28 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 uid-solid system is advanced by employ- ing a semi-implicit staggered scheme as proposed in [20]. At time step t , n n+1 the uid velocity eld is advanced to new time level u ! u from the old n n mask function and the old solid velocity eld u by using the Adams{ n+1 Bashforth scheme. Then, the pressure eld at the new time step p is n+1 calculated from the uid velocity eld u . Finally, the solid is advanced n+1 n+1 to the new time step and u and the whole process is repeated until the nal time. 4.2. Validation For the validation of the uid-solid coupling, we consider two test cases: the Turek benchmark test case FSI3 [48, 49] and the rigid revolving bumble- bee wing test case [6]. 4.2.1. Turek benchmark FSI3 The Turek benchmark FSI3 involves a exible 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 in ow of mean ow U = 2, as shown in gure 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. 29 The uid solver, as well as the uid-solid coupling, are handled by the software FLUSI [4] and the setup remains the same as the test case FSI3 with a resolution of 5200 1152 whose details can be found in [20]. 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 mass- spring network for validation. The results of this simulation are presented in table 3 for the comparison with the reference solutions presented in the literature [20, 48, 49]. For the oscillation of the trailing edge y , the result is in excellent agree- te ment with all three references when the maximum relative error, for both maximum and minimum values of y , is only 1:76% and the relative error te 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 gure 13 to compare with the reference [49]. The two lines are almost super- posed 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 [20], but only 1% comparing with [48]. From gure 13, the curves of the two solutions appear to have the same shape but have some o set. This o set is explained in [20] 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 di erence is around 20% with respect to [48, 49] for both max and min values. Like in Engels [20], the amplitude of the lift force is over-predicted by coupling FLUSI and the mass-spring solver. In conclusion 30 3 References y [10 ] Drag Lift f te 0 max min max min max min Mass-spring network 36.22 -32.93 503.02 442.12 189.94 -186.23 5.56 (1) T. Engels [20] 35.63 -32.71 481.20 432.50 188.52 -181.30 5.44 (2) S. Turek [48] 36.37 -33.45 487.81 432.79 156.13 -151.31 5.47 (3) S. Turek [49] 36.46 -33.52 488.24 432.76 156.40 -151.40 5.47 Table 3: Results of the FSI3 benchmark. we nd 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 uid solver in 2D. 4.2.2. Rigid revolving wing Prior to studying the exibility of the wing, a common test case of a rigid revolving wing is considered to validate the coupling between the uid and the solid solver in 3D, i.e. the mask function generation and the velocity eld of the solid. The setup is taken the same as the one used by Engels et al. [6] as shown in gure 14. The angle of attack is xed at = 45 while the rotation angle (t) is given by t= (t) = e + t (23) The wing is rotated around the center of the computational domain of 31 size 4 4 2, which is discretized using a mesh of 1024 1024 512 grid points. To be consistent with the reference simulation [6], 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 gure 14b. However, the vein structure will not be taken 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 eld 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 = L = 2:817 mg and the time scale is chosen in the way that air wing tip velocity is unity, thus T = 1 s . The Reynolds number is then de ned as in [6] Re = u c = where the mean wingtip velocity u = 1 [LT ] by tip m tip 4 2 1 de nition from eqn. 23, the uid viscosity = 1:477 10 [L T ] and the mean chord, the ratio between the wing surface area A and the wing length R , c = A=R = 0:304 [L]. This yields Re = 2060. Additionally, the lift w m w and drag coecient are de ned as below F F L D C = ; C = (24) L D 2 2 MLT MLT where the lift F is the force in the vertical direction Oz and the drag F L D is the force perpendicular to the plane formed by the vertical and the wing spanwise axes, as shown in gure 14a. The computed lift and drag coecients are shown in gure 15 along with the reference values from [6]. 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 32 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 agree- ments can give us the con dence about the solid solver, based on mass-spring system, as well as the coupling with the ow solver FLUSI. Any di erence between all the numerical studies carried out can be explained by the di er- ence between the continuum model and the discrete model together with the way of generating the mask function. 5. Results and discussion In the following we present results of high resolution computations of revolving bumblebee wings which are either rigid, exible or highly exi- ble. First we perform computations for di erent resolutions to check the mesh convergence for both uid and solid solvers. Then a comparison of the exible wings with the rigid case allows to assess the in uence 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 di erence 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 = L = 4:13 mg while the time scale air remains the same T = 1 s. The corresponding Reynolds number is 1800 4 2 1 where the uid viscosity is assumed to be = 1:477 10 [L T ], the wing tip velocity u = 1 [LT ] and the mean chord calculated from the new tip wing surface area is c = A=R = 0:266 [L]. m w 33 5.1. Study of mesh convergence 5.1.1. Fluid mesh The following mesh convergence study for the uid solver is performed considering ve di erent resolutions: 128 128 64, 256 256 128, 512 512 256, 768 768 384 and 1024 1024 512. 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 nest 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 gure 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 [20]) and equal to 0:052. The drag obtained for each simulation ( gure 16) shows the convergence to the nest resolution solution when we re ne the mesh. The spatial convergence exhibits a rst to second order behavior when we plot the error versus the mesh size. 5.1.2. 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 exible wing at resolution 768 384 are run to compare between a medium-mesh and a ne-mesh wing which are discretized by 465 and 1065 mass points, respec- tively. As shown here in gure 18, although the number of mass points is more than doubled, the forces remain almost unchanged with an increase of 34 1:1% and 0:8% in average lift and drag coecients during the steady state, respectively. Since the uid solver is itself already costly in term of CPU time, the medium-mesh wing with 465 mass points is sucient and can be chosen for the following study in section 5.2. 5.2. In uence of wing exibility To examine the in uence of vein sti ness on the aerodynamic performance of the wing, the exural rigidity of veins will be varied by changing the Young modulus E. Two values of the Young modulus are used: E = 1:25 8 1 2 7 1 2 10 [ML T ] and E = 1:25 10 [ML T ], corresponding to the exible and highly exible cases, respectively. Lift and drag coecients at resolution 1024 512 for the rigid, exible and highly exible cases are presented in gure 19. During the transition phase (rotation angle 40 ), the lift generated by the rigid wing increases instantly and then decreases before going up again. The drag follows the same trend as the lift, but is larger in magnitude. When the exibility of the wing is taken into account, the rapid rise at the beginning of the forces for both exible and highly exible wings disappear. Instead, the forces increase gradually and the more exible the wing is, the lower the lift and the higher the drag are. At steady state, similar behaviors between the rigid and the exible 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 nd that the lift-to-drag ratio at the steady state of the exible 35 wing is 1:2, 14:5% higher than the one of the rigid case, which is only 1:05 ( gure 20). This nding is consistent with conclusions found in literature [50, 19]. A exible wing generates less lift and drag than a rigid one. However, due to the exibility of the wing, the bending in the chordwise direction makes the e ective geometric angle of attack decrease and alters the direction of the total resultant force upward [19]. This makes the lift-to-drag ratio raise and allows better ight performance. On the contrary, the highly exible wing acts di erently. 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 ( gure 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 gure 21 for comparison at three time instants t = 2, t = 4 and t = 6. By applying the functional approach, the di erence between the vein and the membrane is visible in the visualization. At the nest resolution of the mesh (1024 512), the ows generated by the exible wing are shown ( gure 22) by plotting their vorticity magnitude at four time instants of the simulation. The formations of the leading edge 36 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. 6. Conclusions We presented a numerical approach for uid-structure interaction in the open source framework FLUSI, which is based on a mass-spring model de- scribing 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 vol- ume 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 clas- sical benchmarks, a ow past a cylinder with a exible appendage and the ow generated by a rigid revolving wing, is likewise presented. Considering the exible wing, the exibility reduces the buildup of the aerodynamic force during the beginning of motion. Nevertheless, after the start-up phase, the wing yields a steady state con guration, and no signi cant oscillation nor unsteady deformation of the wing are observed. A better aerodynamic performance of the exible wing, characterized by the increase of the lift-to-drag ratio during the steady state, is explained by the decrease of the e ective angle of attack caused by the deformation of the exible wing. On the other hand, the highly exible wing appears to be less ecient than 37 the rigid wing. This can be interpreted that there is an optimized zone of wing exibility, which is ideal for ying. For apping wings we anticipate that exibility will become important because of the dynamic wing deformation. In the near future we are planing high resolution numerical simulations of apping insect ight with exible 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 simula- tions 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 re nement is currently being developed to be able to re- duce memory and CPU time requirements. High resolution numerical sim- ulation of apping ight 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 uid-structure interac- tion on adaptive grids at reduced computational cost. 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 38 were granted access to the HPC resources of IDRIS under the Allocation No. 2018-91664 attributed by GENCI (Grand Equipement National de Calcul Intensif ). For this work, we were also granted access to the HPC resources of Aix-Marseille Universit e nanced by the project Equip@Meso (No. ANR- 10-EQPX- 29-01). The authors thankfully acknowledge nancial support granted by the minist eres des A aires etrang eres et du d eveloppement inter- national (MAEDI) et de l'Education nationale et l'enseignement sup erieur, de la recherche et de l'innovation (MENESRI), and the Deutscher Akademis- cher Austauschdienst (DAAD) within the French-German Procope project FIFIT. D.K. gratefully acknowledges nancial 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. References [1] 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 ight eciency, Science 325 (2009) 1549{1552. [2] T. Nakata, H. Liu, A uid-structure interaction model of insect ight with exible wings, Journal of Computational Physics 231 (2012) 1822{ [3] S. A. Combes, T. L. Daniel, Flexural sti ness in insect wings I. Scaling and the in uence of wing venation, The Journal of Experimental Biology 206 (2003) 2979{2987. [4] T. Engels, D. Kolomenskiy, K. Schneider, J. Sesterhenn, Flusi: A novel 39 parallel simulation tool for apping insect ight using a Fourier method with volume penalization, SIAM J. Sci. Comp. 38 (2016) S03{S24. [5] D. Kolomenskiy, Y. Elimelech, K. Schneider, Leading-edge vortex shed- ding from rotating wings, Fluid Dyn. Res. 46 (2014) 031421. [6] T. Engels, D. Kolomenskiy, K. Schneider, M. Farge, F.-O. Lehmann, J. Sesterhenn, Helical vortices generated by apping wings of bumble- bees, Fluid Dyn. Res. 50 (2018) 011419. [7] T. Engels, D. Kolomenskiy, K. Schneider, M. Farge, F. Lehmann, J. Ses- terhenn, The impact of turbulence on ying insects in tethered and free ight: high-resolution numerical experiments, Phys. Rev. Fluids 4 (2019) 013103. [8] O. Jarrousse, Modi ed mass-spring system for physically based defor- mation modeling, KIT Scienti c Publishing, 2014. [9] D. Terzopoulos, J. Platt, A. H. Barr, K. Fleischer, Elastically deformable models, ACM Siggraph Computer Graphics 21 (1987) 205{214. [10] J. W. Eischen, S. Deng, T. G. Clapp, Finite element modeling and con- trol of exible fabric parts, IEEE Computer Graphics and Applications 16 (1996) 71{80. [11] J. Cai, F. Lin, Y. Lee, Modeling and dynamics simulation for deformable objects of orthotropic materials, Springer Berlin Heidelberg, 2016. [12] A. Nealen, M. Muller, R. Keiser, E. Boxerman, M. Carlson, Physically 40 based deformable models in computer graphics, Computer Graphics forum 25 (2006) 809{836. [13] L. A. Miller, C. S. Peskin, Flexible clap and ing in tiny insect ight, J. Exp. Biol 212 (2009) 30763090. [14] P. Yeh, A. Alexeev, Biomimetic exible plate actuators are faster and more ecient with a passive attachment, Acta Mechanica Sinica 32 (2016) 1001{1011. Doi: 10.1007/s10409-016-0592-0. [15] D. Lentink, M. H. Dickinson, Rotational accelerations stabilize leading edge vortices on revolving y wings, Journal of Experimental Biology 212 (2009) 2705{2719. Doi: 10.1242/jeb.022269. [16] T. Jardin, L. David, Spanwise gradients in ow speed help stabilize leading-edge vortices on revolving wings, Phys. Rev. E 90 (2014) 013011. Doi: 10.1103/PhysRevE.90.013011. [17] A. Jones, A. Medina, H. Spooner, et al., Biomimetic exible plate actuators are faster and more ecient with a passive attachment, Exp Fluids 57 (2016) 1{16. Doi: 10.1007/s00348-016-2143-7. [18] C. Di, D. Kolomenskiy, T. Nakata, H. Liu, Forewings match the forma- tion of leading-edge vortices and dominate aerodynamic force production in revolving insect wings, Bioinspiration & Biomimetics 13 (2017). Doi : 10.1088/1748-3190/aa94d7. [19] R. van de Meerendonk, M. Percin, B. W. van Oudheusden, Three- dimensional ow and load characteristics of exible revolving wings, Experiments in Fluids 59 (2018) 59:161. 41 [20] T. Engels, Numerical modeling of uid-structure interaction in bio- inspired propulsion, PhD Thesis (2015). [21] J. Berger, A second order backward di erence method with variable steps for a parabolic problem, BIT 38 (1998) 644{662. [22] S. A. Combes, T. L. Daniel, Flexural sti ness in insect wings II. Spatial distribution and dynamic wing bending, The Journal of Experimental Biology 206 (2003) 2989{2997. [23] S. M. Walker, A. L. R. Thomas, G. K. Taylor, Deformable wing kinemat- ics in the desert locust: how and why do camber, twist and topography vary through the stroke?, J. R. Soc. Interface 6 (2009) 735{747. [24] S. M. Walker, A. L. R. Thomas, G. K. Taylor, Photogrammetric recon- struction of high-resolution surface topographies and deformable wing kinematics of tethered locusts and free- ying hover ies, J. R. Soc. In- terface 6 (2009) 351366. [25] S. M. Walker, A. L. R. Thomas, G. K. Taylor, Deformable wing kine- matics in free- ying hover ies, J. R. Soc. Interface 7 (2010) 131142. [26] J. R. Bomphrey, N. J. Lawson, N. J. Harding, G. K. Taylor, T. A. L. R., The aerodynamics of manduca sexta: digital particle image velocimetry analysis of leading-edge vortex, J. Exp. Biol 208 (2005) 1079{1094. [27] A. M. Mountcastle, T. L. Daniel, Aerodynamic and functional conse- quences of wing compliance, Exp. Fluids 46 (2009) 873{882. 42 [28] R. T. Fenner, Engineering Elasticity: Application of Numerical and An- alytical Techniques, New York: John Wiley, 1986. [29] R. E. White, An Introduction to the Finite Element Method with Ap- plications to Nonlinear Problems, New York: John Wiley, 1985. [30] M. Chen, F. Boyle, Investigation of membrane mechanics using spring networks: application to red-blood-cell modelling, Materials Science and Engineering 43 (2014) 506{516. [31] T. Omori, T. Ishikawa, D. Barth es-Biesel, A.-V. Salsac, J. Walter, Y. Imai, T. Yamaguchi, Comparison between spring network models and continuum constitutive laws: Application to the large deformation of a capsule in shear ow, Physical Review E 83 (2011) 041918. [32] O. Deussen, L. Kobbelt, P. Tucke, Using simulated annealing to obtain good nodal approximations of deformable bodies, Springer-Verlag (1995) 30{43. [33] B. A. Lloyd, G. Sz ekely, M. Harders, Identi cation of spring parameters for deformable object simulation, IEEE Transactions on Visualization and Computer Graphics 13 (2007). [34] J. Louchet, X. Provot, D. Crochemore, Evolutionary identi cation of cloth animation models, Springer-Verlag (1995) 44{54. [35] G. Bianchi, M. Harders, G. Sz ekely, Mesh topology identi cation for mass-spring models, Proc. Medical Image Computing and Computer- Assisted Intervention (MICCAI 03) 1 (2003). 43 [36] G. Bianchi, B. Solenthaler, G. Sz ekely, M. Harders, Simultaneous topol- ogy and sti ness identi cation for mass-spring models based on fem ref- erence deformations, Proc. Medical Image Computing and Computer- Assisted Intervention (MICCAI 04) 2 (2004). [37] D. L. Logan, A rst course in the Finite Element Method, 5th Revised edition, CL Engineering, 2010. [38] H. J. Barten, On the de ection of a cantilever beam, Quarterly of Applied Mathematics 2 (1944) 168{171. [39] K. E. Bisshopp, D. C. Drucker, Large de ections of cantilever beams, Quart. Appl. Math (1945) 272{275. [40] J. Wang, J. K. Chen, S. Liao, An explicit solution of the large deforma- tion of a cantilever beam under point load at the free tip, J. Comput. Appl. Math. 212 (2006). [41] T. Engels, D. Kolomenskiy, K. Schneider, J. Sesterhenn, Numerical sim- ulation of vortex-induced drag of elastic swimmer models, Theoretical and Applied Mechanics Letters 27 (2017) 280{285. [42] D. Kolomenskiy, S. Ravi, R. Xu, K. Ueyama, T. Jakobi, T. Engels, T. Nakata, J. Sesterhenn, M. Farge, K. Schneider, R. Onishi, H. Liu, The dynamics of passive feathering rotation in hovering ight of bumble- bees., J. Fluids. Struc. (2019). Doi: 10.1016/j.jfluidstructs.2019. 03.021. [43] J. F. V. Vincent, U. G. K. Wegst, Design and mechanical properties of insect cuticle, Arthropod Structure and Development 33 (2004) 187{199. 44 [44] P. Angot, C. Bruneau, P. Fabrie, A penalization method to take into account obstacles in incompressible viscous ows, Numer. Math. (1999) 81:497{520. [45] D. Kolomenskiy, K. Schneider, A Fourier spectral method for the Navier{Stokes equations with volume penalization for moving solid ob- stacles, J. Comput. Phys. (2009) 228:5687{5709. [46] D. Eberly, Distance between point and triangle in 3d, Geometric Tools, LLC (1999). http://www.magic- software.com/Documentation/pt3tri3.pdf. [47] X. Yang, X. Zhang, Z. Li, G.-W. He, A smoothing technique for dis- crete delta functions with application to immersed boundary method in moving boundary simulations, Journal of Computational Physics 228 (2009) 7821{7836. [48] S. Turek, J. Hron, Proposal for numerical benchmarking of uid- structure interaction between an elastic object and laminar incompress- ible ow. In: H.J. Bungartz and M. Sch afer (eds) Fluid-Structure In- teraction. Lecture Notes in Computational Science and Engineering, Springer Berlin Heidelberg 53 (2006) 371{385. [49] S. Turek, J. Hron, M. Razzaq, H. Wobker, M. Sch afer, Numerical bench- marking of uid-structure interaction: A comparison of di erent dis- cretization and solution approaches. In: H.J. Bungartz, M. Mehl and M. Sch afer (eds) Fluid-Structure Interaction II. Lecture Notes in Compu- 45 tational Science and Engineering, Springer Berlin Heidelberg 73 (2010) 413{424. [50] L. Zhao, Q. Huang, X. Deng, S. P. Sane, Aerodynamic e ects of exi- bility in apping wings, J R Soc Interface 7 (2009) 485{497. 46 (a) extension spring (b) bending spring Figure 1: Illustration of the restoring forces corresponding to the deformation of extension and bending springs. Figure 2: All possible shapes of a bending spring corresponding to one bending angle . ijk 47 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: Di erent mesh structures for 2D discretization [31]. Figure 5: Relationship between the spring constant k, Young modulus E , and Poisson ratio 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]. 48 (a) Continuous beam (b) Discrete mass-spring network Figure 6: Illustration of deformation corresponding to forces applied on extension and bending springs. 49 Figure 7: Cantilever beam subjected to gravity eld modeled by continuous nonlinear beam and mass-spring network. (a) (b) Figure 8: The oscillation of the trailing edge y (t) (a) and the de ection lines at t = 2 te (b) calculated by the nonlinear beam theory [20] (dashed blue line) and the mass-spring network (solid red line). 50 Figure 9: Deformation of a 2D sheet along the x-axis under the uniaxial tension T . Figure 10: E ect of mesh topology on the relation between the spring sti ness k and the corresponding Young modulus E 51 Figure 11: Illustration of the mass-spring model which is meshed based on measured data of real bumblebee wings [42]. The black and white markers represent mass centers. Color codes (red, green and blue) are used for identifying veins and the membranes are represented by cyan circles. Figure 12: Computational domain of the FSI3 Turek benchmark and dimensions of the solid part [49]. 52 (a) (b) (c) Figure 13: The oscillation of the trailing edge y (t) (a), the drag (b) and the lift (c) from te the references [49] (blue dashed line), [20] (black dash-dotted line) and the mass-spring network (red continuous line). 53 (a) Scheme of the revolving wing, adapted from [6]. (b) Corresponding discrete mass-spring network. Figure 14: Schematic diagram of the revolving wing simulation (a) and the wing mass- spring model (b) adapted from [6]. The wing is rotated around the hinge point with an angle of attack = 45 . 54 Figure 15: Comparison of lift and drag coecients for a rigid wing, calculated using the coupling between FLUSI and the developed mass-spring solver, with the reference data from [6]. Figure 16: Drag coecient generated by a exible wing, calculated at di erent resolutions: 2 2 2 2 2 128 64, 256 128, 512 256, 768 384 and 1024 512 . 55 Figure 17: The error of the mean drag versus mesh size. The dashed lines represent rst and second order convergence. Figure 18: Lift and drag coecients generated by a revolving exible wing discretized by 465 and 1065 mass points at Re = 1800. 56 (a) Lift coecients. (b) Drag coecient. Figure 19: In uence of wing exibility on lift and drag coecients of a revolving wing at Re = 1800. 57 Figure 20: Lift-to-drag ratio for the three wings during the steady state, rotation angle 120 320 . 58 (a) Top view. (b) Side view. Figure 21: Wing deformation corresponding to rigid (dark blue), exible (blue) and highly exible (light blue) wings at three time instants, t = 2, t = 4 and t = 6. 59 Figure 22: Flows generated by a exible revolving wing, visualized by their vorticity magnitude j!j at four time instants t = 1; 2; 4 and 6. The simulation is performed with resolution 1024 512.
Journal
Physics – arXiv (Cornell University)
Published: Feb 20, 2020