Computational Approach for the Fluid-Structure Interaction Design of Insect-Inspired Micro Flapping Wings

: A ﬂight device for insect-inspired ﬂapping wing nano air vehicles (FWNAVs), which consists of the micro wings, the actuator, and the transmission, can use the ﬂuid-structure interaction (FSI) to create the characteristic motions of the ﬂapping wings. This design will be essential for further miniaturization of FWNAVs, since it will reduce the mechanical and electrical complexities of the ﬂight device. Computational approaches will be necessary for this biomimetic concept because of the complexity of the FSI. Hence, in this study, a computational approach for the FSI design of insect-inspired micro ﬂapping wings is proposed. This approach consists of a direct numerical modeling of the strongly coupled FSI, the dynamic similarity framework, and the design window (DW) search. The present numerical examples demonstrated that the dynamic similarity framework works well to make different two FSI systems with the strong coupling dynamically similar to each other, and this framework works as the guideline for the systematic investigation of the effect of characteristic parameters on the FSI system. Finally, an insect-inspired micro ﬂapping wing with the 2.5-dimensional structure was designed using the proposed approach such that it can create the lift sufﬁcient to support the weight of small insects. The existing area of satisfactory design solutions or the DW increases the fabricability of this wing using micromachining techniques based on the photolithography in the micro-electro-mechanical systems (MEMS) technology. Hence, the proposed approach will contribute to the further miniaturization of FWNAVs.


Introduction
The advantages of flight over other forms of locomotion have resulted in the dispersal of insects on Earth [1]. They can perform extreme aerial maneuvers with robustness using flapping wings [2]. Their flight abilities have become increasingly refined through a long period of natural selection [3]. Hence, emulating nature's time-tested forms, functions, and strategies in flying insects with understanding their underlying mechanisms leads to finding sustainable solutions [4,5], which are engineering alternatives to nature's solutions. This biomimetic concept is multidisciplinary. Hence, computational approaches for this biomimetic concept will be essential because of their interdisciplinary feature, leading to the computational biomimetics.
This study focuses on insect-inspired flapping wing nano air vehicles (FWNAVs) whose sizes range from approximately 1 to 10 cm [6], on which various research groups are working as one of most popular topics in the biomimetic research area. One engineering concern in the development of FWNAVs is the reduction of their size [7,8], potentially down to the size scale of very small insects, of which wingspan length is less than 1 mm as the result of the evolution [9].
Similar to the structural design of micro-electro-mechanical systems [10], the conceptual, basic, and detailed design processes of FWNAVs will be indistinguishable from each other, because of the novelty of the biomimetic concept, and the formulation as a single design problem is difficult in each design process. Furthermore, the feasibility of the design solution is strongly constrained by the manufacturability, which is hardly formulated as the design problem.
The flight ability of the insects seems to be determined by the characteristic motions of their flapping wings such as the pitching, the cambering, and the wing's tip path [11]. In particular, the pitching motion will be the kinematical basis of the aerodynamics of the insect flapping flight [12]. The insect wings and their base are flexible [13][14][15][16][17][18][19]. Hence, their flapping motion with the large stroke angle [20] causes their deformation due to the aerodynamic force from the surrounding air, and their deformation causes the change of the surrounding air flow, which results in the change of the aerodynamic force. Recent advances on this topic are given mainly by the computational approaches using the computational fluid dynamics [21][22][23][24][25] and the numerical fluid-structure interaction (FSI) analysis [26][27][28][29][30][31][32][33]. Some studies used the experimental approaches [22,29,34]. Some of them have shown that the wing's characteristic motions can be caused by the strong coupling between the wing and the surrounding air [26,[28][29][30]32,34].
In the FSI design [35], a flight device for FWNAVs, which consists of the micro wings, the actuator, and the transmission [6], actively uses the FSI to create the characteristic wing's motions. It can reduce the mechanical and electrical complexities of this device [36]. Hence, it is essential for further miniaturization of FWNAVs. Computational approaches are necessary for this biomimetic concept because of the complexity of the FSI. However, there are few studies focusing on this topic.
Hence, in this study, a computational approach for the FSI design of insect-inspired micro flapping wings is proposed. This approach consists of a direct numerical modeling of the strongly coupled FSI, the dynamic similarity framework, and the design window (DW) search. Here, the term "direct" means that, prior to the finite element discretization, the governing equations are not transformed into a different form [37].
The key of the FSI design is the accurate prediction for three-dimensional (3-D) strong coupling of a thin flexible flapping wing and the surrounding air flow. Here, a projection method using the algebraic splitting [38,39] is used to solve this problem directly.
The dynamic similarity law for the FSI [26,29,30,32,34] is used as the framework for the systematic treatment of FSI complexities. The validity of this framework is demonstrated in the first numerical example. As shown in the second numerical example, the simulation based on the direct numerical modeling for the strong coupling of the lumped flexibility model [26,32] and the surrounding fluid is guided by this framework, and the underlying mechanism is effectively investigated, which is the basis of the FSI design.
DWs are the existing areas of satisfactory solutions in a design parameter space. They can adapt to changes in the design problem. Hence, the DW search is effective for the above-mentioned issues in the design of FWNAVs.
Finally, in the last numerical example, an insect-inspired micro flapping wing with the 2.5-dimensional (2.5-D) structure is designed using the proposed approach. The 2.5-D structure is multi-layered, and each layer is given by laminating the membrane or the thin plate on the flat layer already developed and etching it to form the in-plane shape. This structure can be fabricated using the micromachining based on the photolithography in the micro-electro-mechanical systems (MEMS) technology. However, there exists the strong restriction on the shape in the out-of-plane direction. The DW is determined such that each solution in the DW can create the lift sufficient to support the weight of small insects. This wing can be fabricated using the micromachining based on the photolithography [7,8,40,41]. Hence, the proposed approach will contribute to the further miniaturization of FWNAVs.

Model Wing for the Fluid-Structure Interaction Design
In the FSI design [35], a flight device for FWNAVs consisting of the micro flapping wings, the actuator, and the transmission [6] actively uses the FSI to create the characteristic motions of the wings. This concept will reduce the mechanical and electrical complexities of this device [36], which leads to further miniaturization of FWNAVs.
The wings of actual insects consist of multiple and multilayered structures spanning a wide range of length scale [13,[42][43][44][45]. The wing's veins and membranes keep enough strength against loads during the flight [13], while these structures allow their characteristic deformations such as the twisting and the cambering [11,46], which improve the aerodynamic performance [22][23][24][25]27,31]. Recent advances on this topic such as References [18,19] will improve the understanding for the material mechanics of the insects' wings.
The wing base consists of steering muscles and tendons with both flexible and rigid parts that work together as a transmission to redirect power from the main flight muscles to the wing's flapping while to allow the fine-tuned wing's control [47][48][49][50].
It seems that there exist two approaches for modeling these structures [51]. One is the reduced-order modeling [26,27,29,32,[52][53][54], and the other is the realistic modeling [55][56][57]. In References [26,29,32], the lumped flexibility model is used, where the spring encapsulates the detailed structures as the macroscopic constitutive relationship that describes their elastic mechanism. In Reference [29], the continuous constitutive relationship like anisotropic materials is used. In References [52][53][54], the structural elements such as beam, membrane, and shell are used. The reduced-order modeling can be incorporated in the computational FSI analysis using the current computer resources [26,27,29,32]. Hence, in this study, the lumped flexibility model is used.
The basic model wing for the FSI design is shown in Figure 1, where the spring for the torsion corresponds to the torsional flexibility that generates the wing's pitching motion [26,29], the spring for the elevation corresponds to the elevation flexibility that generates the wing's tip path [32], and the spring for the flapping corresponds to the transmission for the wing's flapping. Furthermore, the veins are represented by three beams according to their functions to form the camber [52][53][54]. This model wing emulates the characteristic wing's motions with their elastic mechanisms and leads to the engineering alternatives to the nature's solution.

Model Wing for the Fluid-Structure Interaction Design
In the FSI design [35], a flight device for FWNAVs consisting of the micro flapping wings, the actuator, and the transmission [6] actively uses the FSI to create the characteristic motions of the wings. This concept will reduce the mechanical and electrical complexities of this device [36], which leads to further miniaturization of FWNAVs.
The wings of actual insects consist of multiple and multilayered structures spanning a wide range of length scale [13,[42][43][44][45]. The wing's veins and membranes keep enough strength against loads during the flight [13], while these structures allow their characteristic deformations such as the twisting and the cambering [11,46], which improve the aerodynamic performance [22][23][24][25]27,31]. Recent advances on this topic such as References [18,19] will improve the understanding for the material mechanics of the insects' wings.
The wing base consists of steering muscles and tendons with both flexible and rigid parts that work together as a transmission to redirect power from the main flight muscles to the wing's flapping while to allow the fine-tuned wing's control [47][48][49][50].
It seems that there exist two approaches for modeling these structures [51]. One is the reduced-order modeling [26,27,29,32,[52][53][54], and the other is the realistic modeling [55][56][57]. In References [26,29,32], the lumped flexibility model is used, where the spring encapsulates the detailed structures as the macroscopic constitutive relationship that describes their elastic mechanism. In Reference [29], the continuous constitutive relationship like anisotropic materials is used. In References [52][53][54], the structural elements such as beam, membrane, and shell are used. The reduced-order modeling can be incorporated in the computational FSI analysis using the current computer resources [26,27,29,32]. Hence, in this study, the lumped flexibility model is used.
The basic model wing for the FSI design is shown in Figure 1, where the spring for the torsion corresponds to the torsional flexibility that generates the wing's pitching motion [26,29], the spring for the elevation corresponds to the elevation flexibility that generates the wing's tip path [32], and the spring for the flapping corresponds to the transmission for the wing's flapping. Furthermore, the veins are represented by three beams according to their functions to form the camber [52][53][54]. This model wing emulates the characteristic wing's motions with their elastic mechanisms and leads to the engineering alternatives to the nature's solution. The basic model wing for the FSI design with the springs of flapping, torsion, and elevation, which correspond to the axes of flapping, torsion, and elevation, respectively. Φ is the stroke angle, φ is the stroke angular displacement, θ is the pitching angular displacement, and χ is the elevation angular displacement.  ϕ is the stroke angular displacement, θ is the pitching angular displacement, and χ is the elevation angular displacement.

Fluid-Structure Interaction Analysis
The key of the FSI design is the accurate prediction for the 3-D strong coupling of the thin flexible flapping wing and the surrounding air flow. Here, a projection method using the algebraic splitting [38,39] solves this problem directly. This method was validated using the benchmarks in References [38,39]. The computer program used in this study was also  [29] using the actual experiment of the dynamically scaled model, which is similar to the FSI problem in this study.

Governing Equations for the FSI
The equilibrium equation for the elastic wing can be expressed as where d/dt in the left-hand side is the Lagrangian time derivative, the superscript s indicates a quantity corresponding to the structure, ρ s is the mass density of the structure, v s i is the i-th component of the structural velocity vector, and σ s ij is the ij-th component of the Cauchy stress tensor of the structure. Note that i equals 1, 2, or 3, and corresponds to the x, y, or z direction, respectively, in the Cartesian coordinate system.
The incompressible Navier-Stokes (NS) equations for the fluid can be expressed as where the arbitrary Lagrangian-Eulerian (ALE) method is used, ∂/∂t in the left-hand side is the ALE time derivative, the superscript f indicates a quantity corresponding to the fluid, The interface conditions describing the interaction between the wing and the surrounding fluid can be expressed as where n f i and n s i are the i-th components of the outward unit normal vectors on the fluid-structure interface corresponding to the fluid and the structure, respectively.

Monolithic Equation System for the FSI
Applying finite element discretization to the total Lagrangian formulation of Equation (1), the nonlinear equilibrium equation system is obtained in matrix form. Similarly, applying finite element discretization to Equations (2) and (3), the nonlinear equation system is obtained in matrix form. Applying the fluid-structure interface conditions (4) and (5) to the spatially discretized governing equations, the monolithic equilibrium equations and the incompressibility constraint for the fluid can be obtained as follows: where M is the mass matrix, C is the diffusive matrix, G is the divergence operator matrix, N is the convective term vector, q is the elastic internal force vector, g is the external force vector, a is the acceleration vector, v is the velocity vector, u is the displacement vector, p is the pressure vector, the subscript L indicates the lumping of the matrix, and the subscript τ indicates the transpose of the matrix. The monolithic Equations (6) and (7) are linearized using the increments of the state variables as follows: where the pressure term and the elastic interior force term are evaluated implicitly, M* is the generalized mass matrix, ∆t denotes the time increment, ∆g and ∆h are the residual vectors of the equilibrium Equation (6) and the incompressibility constraint (7), respectively, G ε comes from the pressure stabilization term of the PSPG method. The relations between the increments of the state variables are given as ∆u = β∆t 2 ∆a and ∆v = γ∆t∆a using the Newmark's β method. The fluid convection and diffusion terms are evaluated explicitly to reduce the computational cost. Therefore, M* is given as L M + β∆t 2 K, while the Courant's and diffusion number conditions are imposed on the time increment ∆t for the stability of the time integration. The predictor-multicorrector algorithm (PMA) is used for the time integration.

Projection Method Using the Algebraic Splitting
The monolithic method solves the linearized monolithic Equations (8) and (9) simultaneously. Hence, this method can satisfy the interface conditions (4) and (5) to avoid spurious numerical power on the fluid-structure interface, which yields numerical instability. However, the formulation might lead to ill-conditioned coefficient matrices. In this study, the projection method using the algebraic splitting [38,39] is used to avoid this difficulty. This method is briefly described as follows: The state variables are predicted as the intermediate state variables from the equilibrium Equation (6) for the known pressure, which can be linearized as where a* is the intermediate acceleration or the acceleration predictor, and ∆ is the increment from the known variable. Subtracting both sides of Equation (10) from those of Equation (8) gives, after suitable rearrangement, where v* is the intermediate velocity or the velocity predictor. The relations between the increments of the intermediate state variables are given as ∆u* = β∆t 2 ∆a* and ∆v* = γ∆t∆a* using the Newmark's β method. Left multiplying both sides of Equation (11) by τ G L M −1 , the following equation is obtained: where M * is defined as M * − L M. If the following pressure Poisson equation (PPE) is solved, then Equation (12) is reduced as If the linear convergence of the state variables is assumed, v* agrees with v asymptotically in the nonlinear iterations. Therefore, the second term of Equation (14) will vanish asymptotically, and the incompressibility constraint for the unknown fluid velocity will be satisfied.
It follows from the above formulation that the monolithic equation system is split into the two equilibrium Equations (8) and (10), and the PPE (13). In the nonlinear iterations, the linearized equilibrium equation for the previous pressure (10) is solved to derive the intermediate velocity, the PPE (13) is solved to determine the current pressure such that the current velocity field satisfies the incompressibility constraint (7), and the linearized equilibrium equation for the current pressure (8) is solved to derive the current velocity.

Parallel Computation
The following matrix-vector product is considered, since it is most expensive computation in the iterative solvers for Equations (8), (10) and (13): Assuming that the number of degrees of freedoms (DOFs) of the shell structure is far smaller than that of the fluid, the following parallel solution procedure is used based on the mesh decomposition. The mesh is decomposed as shown in Figure 2. The symbols Ω S and Ω F denote the structural and fluid meshes, respectively, and Ω Fi (i = 1, 2, . . . , N d ) denotes i-th submesh of Ω F . Note that Ω S is set inside Ω F1 as shown in Figure 2 under the above-mentioned assumption. In this case, the matrix-vector product concerning the fluid-structure interface can be computed at one computational node, and the complexity of the computer implementation is reduced. The computation concerning Ω Fi is executed at the computational node P i (i = 1, 2, . . . , N d ), while the computation concerning Ω S is executed at P 1 . Using the above setup, the matrix-vector product (15) is computed as follows:

•
Step 1. The matrix-vector product (15) is independently computed at P i (i = 1, 2, . . . , N d ) using the element-by-element method as where A denotes the global matrix, p denotes the global vector, q denotes their matrixvector product, A (e) and p (e) are their elemental counterparts, respectively, and e denotes the element number.

•
Step 2. The nodal data of the matrix-vector product (16) on the interface between Ω Fi and Ω Fj (j = i) computed at P j is transferred to P i and added to complete the corresponding nodal data computed in Step 1.
asymptotically, and the incompressibility constraint for the unknown fluid velocity will be satisfied. It follows from the above formulation that the monolithic equation system is split into the two equilibrium Equations (8) and (10), and the PPE (13). In the nonlinear iterations, the linearized equilibrium equation for the previous pressure (10) is solved to derive the intermediate velocity, the PPE (13) is solved to determine the current pressure such that the current velocity field satisfies the incompressibility constraint (7), and the linearized equilibrium equation for the current pressure (8) is solved to derive the current velocity.

Parallel Computation
The following matrix-vector product is considered, since it is most expensive computation in the iterative solvers for Equation (8), (10) and (13): Assuming that the number of degrees of freedoms (DOFs) of the shell structure is far smaller than that of the fluid, the following parallel solution procedure is used based on the mesh decomposition. The mesh is decomposed as shown in Figure 2. The symbols ΩS and ΩF denote the structural and fluid meshes, respectively, and ΩFi (i = 1, 2, …, Nd) denotes i-th submesh of ΩF. Note that ΩS is set inside ΩF1 as shown in Figure 2 under the above-mentioned assumption. In this case, the matrix-vector product concerning the fluid-structure interface can be computed at one computational node, and the complexity of the computer implementation is reduced. The computation concerning ΩFi is executed at the computational node Pi (i = 1, 2, …, Nd), while the computation concerning ΩS is executed at P1. Using the above setup, the matrix-vector product (15) is computed as follows: • Step 1. The matrix-vector product (15) is independently computed at Pi (i = 1, 2, …, Nd) using the element-by-element method as where A denotes the global matrix, p denotes the global vector, q denotes their matrixvector product, A (e) and p (e) are their elemental counterparts, respectively, and e denotes the element number.

•
Step 2. The nodal data of the matrix-vector product (16) on the interface between ΩFi and ΩFj (j ≠ i) computed at Pj is transferred to Pi and added to complete the corresponding nodal data computed in Step 1.

Shell Modeling of the Thin Elastic Wing
In this study, mixed interpolation of tensorial components (MITC) shell elements [58,59] are used for a thin elastic wing, while linear equal-order-interpolation velocity-pressure elements [60] are used for the fluid. The same arrangements of the fluid and shell nodes on the fluid-structure interface [61] are used. Hence, a pair of two overlapping fluid nodes per one shell node on the fluid-structure interface is necessary to describe the pressure discontinuity between the obverse and reverse of the shell structure. The interface conditions (4) and (5)  where Q is the internal force vector including all effects of the fluid or the structure, T is the interpolation matrix, and the subscript c indicates the coupled DOFs. The nodal expressions of these equations are where v, Q, and g are the nodal components of the velocity, and the internal and external force vectors, respectively, and the subscripts o and r indicate the obverse and reverse sides of the shell structure, respectively.

Dynamic Similarity Framework for the FSI
The dynamic similarity law can be used to correctly incorporate morphological and kinematical data from actual insects into model wings. Initially, the dynamic similarity law for the fluid dynamics was used for the experimental model [12,62]. The dynamic similarity law for the FSI was first introduced in the numerical model [26], and then it was used for the experimental model [34]. In this study, it is used as the framework for the systematic treatment for the complexities of the FSI.
The nondimensional parameter α, the Reynolds number Re, the Cauchy number Ca, which represents the ratio of the fluid dynamic pressure to the structural elastic force, and the mass ratio M R are obtained from the dimensional analysis for the governing Equations (1)-(4) with respect to the characteristic length L, the characteristic velocity V, and the characteristic time T as follows: where µ f is the fluid viscosity, C is the wing's compliance, and m w is the wing's mass. The parameter (21) is the ratio of the inertial force due to the local acceleration to that due to the convective acceleration. As shown in Figure 3, these four nondimensional parameters determine all ratios among all force terms appear in each FSI system. Hence, each FSI system can be specified using these four nondimensional parameters, and the dynamic similarity between two different FSI systems can be measured as their distance in the nondimensional parameter space.

Design Window (DW) Search
Similar to the structural design of micro-electro-mechanical systems [10], the conceptual, basic, and detailed design processes of FWNAVs will be indistinguishable from each other. Hence, the formulation as a single design problem is difficult in each design process. Furthermore, the feasibility of the design solution is strongly constrained by the manufacturability, which is hardly formulated as the design problem.

Design Window (DW) Search
Similar to the structural design of micro-electro-mechanical systems [10], the conceptual, basic, and detailed design processes of FWNAVs will be indistinguishable from each other. Hence, the formulation as a single design problem is difficult in each design process. Furthermore, the feasibility of the design solution is strongly constrained by the manufacturability, which is hardly formulated as the design problem.
In this study, DWs are searched in the design parameter space using the simulation based on the direct numerical modeling of the FSI. The DW is an existing area of satisfactory solutions in a design parameter space. The final solution can be adaptively chosen from the DWs taking into account other conditions such as the change given to the design problem and the manufacturability.

Problem Setup
The lumped torsional flexibility model is used in this study. The concept of this model is illustrated in Figure 4. In this model, the characteristic quantities L, T, V, and C in Section 2.3 are given as the mean chord length c w , the mean flapping velocity V m , the inverse of the flapping frequency 1/f ϕ , and the wing's torsional compliance C θ , respectively. V m can be expressed as 2r 2 ΦL w f ϕ , where r 2 is the nondimensional radius of the second moment of the wing's area, Φ is the stroke angle, and L w is the wing's longitudinal length. Then, Equations (21)-(24) are reduced to

Design Window (DW) Search
Similar to the structural design of micro-electro-mechanical systems [10], the conceptual, basic, and detailed design processes of FWNAVs will be indistinguishable from each other. Hence, the formulation as a single design problem is difficult in each design process. Furthermore, the feasibility of the design solution is strongly constrained by the manufacturability, which is hardly formulated as the design problem.
In this study, DWs are searched in the design parameter space using the simulation based on the direct numerical modeling of the FSI. The DW is an existing area of satisfactory solutions in a design parameter space. The final solution can be adaptively chosen from the DWs taking into account other conditions such as the change given to the design problem and the manufacturability.

Problem Setup
The lumped torsional flexibility model is used in this study. The concept of this model is illustrated in Figure 4. In this model, the characteristic quantities L, T, V, and C in Section 2.3 are given as the mean chord length cw, the mean flapping velocity Vm, the inverse of the flapping frequency 1/fφ, and the wing's torsional compliance Cθ, respectively. Vm can be expressed as 2r2ΦLwfφ, where r2 is the nondimensional radius of the second moment of the wing's area, Φ is the stroke angle, and Lw is the wing's longitudinal length. Then, Equations (21)-(24) are reduced to where r A is the wing's aspect ratio defined as 2L w /c w . As an implementation of the lumped torsional flexibility model, the rectangular wing and the plate spring are used as shown in Figure 5a. r 2 of the rectangular wing is 1/3 0.5 , which is close to that of many insect wings [63]. Hence, this shape is used for the purpose of simplicity. The finite element model is shown in Figure 5b, where MITC shell elements [58,59] are used. The fluid domain is rectangular as shown in Figure 6a. The no-slip condition is imposed to all boundaries on the wall. The finite element mesh for the fluid domain is shown in Figure 6b,d, where linear equal-order-interpolation velocity-pressure elements [60] are used. Based on the observation of actual insects [11], the time history of the flapping angular velocity is given using the trapezoidal function. Following this time history, the time history of the flapping angular displacement is given as shown in Figure 7, which is close to the sinusoidal motion after the initial transient phase.
tion of actual insects [11], the time history of the flapping angular velocity is given using the trapezoidal function. Following this time history, the time history of the flapping angular displacement is given as shown in Figure 7, which is close to the sinusoidal motion after the initial transient phase.
Two different FSI systems are considered here. One is the insect wing, and the other is the dynamically scaled wing. The characteristic length of the dynamically scaled wing is 20 times larger than that of the insect wing. Following this, the values of the characteristic parameters are converted from the insect wing to the dynamically scaled wing using Equations (25)(26)(27)(28). The values of the characteristic parameters for these wings are summarized in Table 1.       Two different FSI systems are considered here. One is the insect wing, and the other is the dynamically scaled wing. The characteristic length of the dynamically scaled wing is 20 times larger than that of the insect wing. Following this, the values of the characteristic parameters are converted from the insect wing to the dynamically scaled wing using Equations (25)- (28). The values of the characteristic parameters for these wings are summarized in Table 1.

Results and Discussion
The purpose of this test is to demonstrate that the present dynamic similarity framework can make two different FSI systems dynamically similar to each other. For this purpose, the pitching angular displacement and the nondimensional lift are used, since the former is the representative characteristic quantity of the structure in the passive pitching motion, and the latter is the representative characteristic quantity of the surrounding fluid. Here, the lift force is defined as the total fluid force acting on the wing in the direction perpendicular to the stroke plane. The present FSI system is geometrically similar to that in Reference [29], and their nondimensional parameters (α, Re, Ca, and M R ) are close to those in Reference [29]. The mesh used in Reference [29] had the sufficient resolution for the lift force. The mesh used in this study has the finer resolution. Figure 8a,b show the time histories of the pitching angle and the nondimensional lift, respectively. As shown in Figure 8a, the transient large deformations in the insect wing and the dynamically scaled wing agree well with each other. Similarly, as shown in Figure 8b, the surrounding flow dynamics of these two wings agree well with each other. That is, these wings strongly coupled with the surrounding fluid flows are dynamically similar to each other. Ca represents the ratio of the fluid dynamic pressure to the structural elastic force. In Figure 8a, Ca contributes the equivalence of the pitching angles, since they are determined by the equilibrium between the wing and the surrounding fluid. It follows from these results that the present dynamic similarity framework works well to make two different FSI systems dynamically similar to each other. Figure 8a,b show the time histories of the pitching angle and the nondimensional lift, respectively. As shown in Figure 8a, the transient large deformations in the insect wing and the dynamically scaled wing agree well with each other. Similarly, as shown in Figure  8b, the surrounding flow dynamics of these two wings agree well with each other. That is, these wings strongly coupled with the surrounding fluid flows are dynamically similar to each other. Ca represents the ratio of the fluid dynamic pressure to the structural elastic force. In Figure 8a, Ca contributes the equivalence of the pitching angles, since they are determined by the equilibrium between the wing and the surrounding fluid. It follows from these results that the present dynamic similarity framework works well to make two different FSI systems dynamically similar to each other.

Effect of Flapping Frequency on Passive Wing Motion
One of key parameters in insect flapping wings will be the flapping frequency. However, multiple material and kinematical properties, and complicated fluid and structural forces might prevent us from systematic understanding of its effect on the wing's passive pitching motion. The dynamic similarity law can reduce these properties and forces to only four nondimensional parameters. Hence, the framework using this law will work as the guideline for the systematic investigation of the effect of the flapping frequency on the passive pitching motion.

Problem Setup
The basic setup is the same as that in Section 3.1. The mesh used here is the same one as Reference [29]. The nondimensional parameters for the cranefly are used as the reference values, which are α = 0.073, Re = 254, Ca = 0.19, MR = 16, Φ = 123°, and rA = 10.7 [29]. Here, the following dynamic similarity framework is used for the systematic investigation of the effect of the flapping frequency fφ on the passive pitching motion:

Effect of Flapping Frequency on Passive Wing Motion
One of key parameters in insect flapping wings will be the flapping frequency. However, multiple material and kinematical properties, and complicated fluid and structural forces might prevent us from systematic understanding of its effect on the wing's passive pitching motion. The dynamic similarity law can reduce these properties and forces to only four nondimensional parameters. Hence, the framework using this law will work as the guideline for the systematic investigation of the effect of the flapping frequency on the passive pitching motion.

Problem Setup
The basic setup is the same as that in Section 3.1. The mesh used here is the same one as Reference [29]. The nondimensional parameters for the cranefly are used as the reference values, which are α = 0.073, Re = 254, Ca = 0.19, M R = 16, Φ = 123 • , and r A = 10.7 [29]. Here, the following dynamic similarity framework is used for the systematic investigation of the effect of the flapping frequency f ϕ on the passive pitching motion: The frequency ratio F r is used as the control parameter, which is defined as f n /f ϕ , where f n is the natural frequency of the wing's torsion. F r is changed such that the product αF r is kept constant (Φ and f ϕ are changed to change α and F r , respectively), while Re, Ca, and M R are kept constant. The four cases shown in Table 2 are investigated.  Figure 9a,b show the time histories of the pitching angular displacement and the nondimensional lift, respectively. In the middle of each half stroke in insect flapping wings, the dynamic pressure and the elastic force are dominant or the fluid and structural inertial forces are significantly damped [30]. Here, Ca values are same among different F r . Hence, in Figure 9a, the pitch angles in the middle of each half stroke are almost equivalent among different F r . On the contrary, the first peak of the pitch angle in each half stroke increases following the decrease of F r . This can be explained as follows:

Results and Discussion
at the time instant of 6.65 cycle when the nondimensional lift takes the first peak in the second half stroke. Note that the flow fields are considered on the cylindrical planes, of which normal direction coincides with the wing's longitudinal direction. The objectives of these figures are to show the present analysis captured the most important flow structure and to explain why the first lift peak during each half stroke in the case of Fr = 3.3 is larger than that in the case of Fr = 5.1. The characteristic vortex on the leading edge or the leading-edge vortex, which is one of most important flow structures in insect flapping wings [62], was formed clearly in each figure. The large translational lift in the middle of each half stroke is given by the stable attachment of this vortex on the wing. As shown in Table 2, α increases as F r decreases. The increase of α corresponds to the increase of the ratio of the inertial force to the convective and elastic forces, since M R , Re, and Ca are kept constant. The inertial effect for the pitching motion is most significant just after the stroke reversal [30]. Hence, the first peak of the pitch angle in each half stroke increases following the decrease of F r .
As shown in Figure 9b, the first peak of the nondimensional lift in each half stroke increases as F r decreases. Figure 10a,b show the flow fields for F r = 3.3 and 5.1, respectively, at the time instant of 6.65 cycle when the nondimensional lift takes the first peak in the second half stroke. Note that the flow fields are considered on the cylindrical planes, of which normal direction coincides with the wing's longitudinal direction. The objectives of these figures are to show the present analysis captured the most important flow structure and to explain why the first lift peak during each half stroke in the case of F r = 3.3 is larger than that in the case of F r = 5.1. The characteristic vortex on the leading edge or the leading-edge vortex, which is one of most important flow structures in insect flapping wings [62], was formed clearly in each figure. The large translational lift in the middle of each half stroke is given by the stable attachment of this vortex on the wing. As shown in these figures, the magnitude of the downflow in front of the wing in the case of Fr = 3.3 is larger than that in the case of Fr = 5.1. That is, the downward momentum given to the fluid from the wing's flapping motion in the case of Fr = 3.3 is larger than that for Fr = 5.1 at this time instant. As the reaction from the fluid, the first lift peak during each half stroke in the case of Fr = 3.3 is larger than that in the case of Fr = 5.1. As shown in these figures, the magnitude of the downflow in front of the wing in the case of F r = 3.3 is larger than that in the case of F r = 5.1. That is, the downward momentum given to the fluid from the wing's flapping motion in the case of F r = 3.3 is larger than that for F r = 5.1 at this time instant. As the reaction from the fluid, the first lift peak during each half stroke in the case of F r = 3.3 is larger than that in the case of F r = 5.1.

Micro Flapping Wing with 2.5-D Structure
In this section, the possibility of the micro flapping wing for FWNAVs is demonstrated using the proposed approach. Figure 11 shows the concept of the proposed micro flapping wing using the FSI design. As shown in this figure, the pitching motion is caused by the FSI. Furthermore, the base excitation given by a micro actuator is amplified to obtain the sufficient amplitude of the flapping motion using the resonance. Here, the flapping frequency is chosen as the design parameter since it is a key parameter for the passive pitching motion as shown in the previous section. As illustrated in Figure 12a as well as Figure 11, the most significant feature of the proposed wing is the 2.5-D structure. This structure is the engineering alternative of the nature's solution such that it can be fabricated using the micromachining technique based on the photolithography in the MEMS technology. Figure 12b shows the prototyping of this wing using a polymer micromachining [7,8,40].
(a) (b) Figure 10. Flow fields on the cylindrical plane crossing the flapping wing at the time instant of 6.65 cycle when the nondimensional lift takes the first peak in the second half stroke: (a) Φ = 123°, Fr = 5.1, α = 0.073; (b) Φ = 80°, Fr = 3.3, α = 0.11. In each cylindrical plane, the color contour shows the magnitude of the vorticity in the wing's longitudinal direction, the arrows show the fluid velocity field, and the color of each arrow shows the magnitude of the velocity.
As shown in these figures, the magnitude of the downflow in front of the wing in the case of Fr = 3.3 is larger than that in the case of Fr = 5.1. That is, the downward momentum given to the fluid from the wing's flapping motion in the case of Fr = 3.3 is larger than that for Fr = 5.1 at this time instant. As the reaction from the fluid, the first lift peak during each half stroke in the case of Fr = 3.3 is larger than that in the case of Fr = 5.1.

Micro Flapping Wing with 2.5-D Structure
In this section, the possibility of the micro flapping wing for FWNAVs is demonstrated using the proposed approach. Figure 11 shows the concept of the proposed micro flapping wing using the FSI design. As shown in this figure, the pitching motion is caused by the FSI. Furthermore, the base excitation given by a micro actuator is amplified to obtain the sufficient amplitude of the flapping motion using the resonance. Here, the flapping frequency is chosen as the design parameter since it is a key parameter for the passive pitching motion as shown in the previous section. As illustrated in Figure 12a as well as Figure 11, the most significant feature of the proposed wing is the 2.5-D structure. This structure is the engineering alternative of the nature's solution such that it can be fabricated using the micromachining technique based on the photolithography in the MEMS technology. Figure 12b shows the prototyping of this wing using a polymer micromachining [7,8,40].

Problem Setup
The wing membrane is made of the polyimide (the mass density ρ s = 1.43 g/cm 3 , the Young's modulus E = 3GPa, and the Poisson's ratio ν = 0.4), the span length Lw = 2.5 mm, the chord length cw = 0.8 mm, and the thickness tw = 1.6μm, which are comparable to those of small flies. The leading edge is made of the single crystal silicon (the mass density ρ s = 2.383 g/cm 3 , the Young's modulus E = 180GPa, and the Poisson's ratio ν = 0.3), and the cross section is 100 μm width and 50 μm thickness. The plate spring is made of the same material of the wing membrane, and the spring length ls is set as 50μm. The surrounding fluid is considered as air (the mass density ρ f = 1.18 × 10 −3 g/cm 3 , the viscosity μ = 1.82 × 10 −4 g/(cm × s)).
The amplitude and flapping frequency of the micro actuator are assumed as follows: The amplitude of the base excitation displacement in the x-direction u0 (see Figure 11) is set as 80 μm, which should be smaller than 100 μm taking into account the available micro actuator [7,8], and the flapping frequency fφ is chosen from 100 Hz to 1000 Hz based on

Problem Setup
The wing membrane is made of the polyimide (the mass density ρ s = 1.43 g/cm 3 , the Young's modulus E = 3 GPa, and the Poisson's ratio ν = 0.4), the span length L w = 2.5 mm, the chord length c w = 0.8 mm, and the thickness t w = 1.6 µm, which are comparable to those of small flies. The leading edge is made of the single crystal silicon (the mass density ρ s = 2.383 g/cm 3 , the Young's modulus E = 180 GPa, and the Poisson's ratio ν = 0.3), and the cross section is 100 µm width and 50 µm thickness. The plate spring is made of the same material of the wing membrane, and the spring length l s is set as 50 µm. The surrounding fluid is considered as air (the mass density ρ f = 1.18 × 10 −3 g/cm 3 , the viscosity µ = 1.82 × 10 −4 g/(cm × s)). The amplitude and flapping frequency of the micro actuator are assumed as follows: The amplitude of the base excitation displacement in the x-direction u 0 (see Figure 11) is set as 80 µm, which should be smaller than 100 µm taking into account the available micro actuator [7,8], and the flapping frequency f ϕ is chosen from 100 Hz to 1000 Hz based on actual insects. The objective of the design is to find the solution that can generate the mean lift F L larger than 7 µN, which is equivalent to the weight of small flies.
MITC shell elements [58,59] are used for a thin elastic structure, which consists of the beam of the leading edge, the plate spring, and the wing membrane (number of nodes: 225, number of elements: 196). As shown in Figure 13b, the plate spring is modeled as a part of the beam of the leading edge at its base. In this modeling, the Young's modulus is used to keep the spring constant.  The fluid domain is rectangular as shown in Figure 13a. The no-slip condition is imposed to all boundaries on the wall. Stabilized linear equal-order-interpolation velocitypressure elements [60] are used for the fluid. The fluid mesh is shown in Figure 13c,d. In this mesh, the number of nodes is 46,911 and the number of elements is 254,352. Note that the resolution of the fluid mesh is almost equivalent to that in Reference [29]. The time increment ∆t is set as T ϕ /5000, where T ϕ is the flapping period. For the parametric study, the parallel computation in Section 2.2.4 is executed in a multiple core processor (10-core Xeon 2.8 GHz × 2 CPUs, 32 GB memory).

Results and Discussion
In the time histories of the x-displacement of the wing's tip and base in the case of f ϕ = 428 Hz, the wing tip displacement was 19 times larger than the wing base displacement because of the resonance as shown in Figure 14a. Figure 14b shows the time history of the flapping angular displacement. The stroke angle reaches about 73 • at the quasi-steady state, which is comparable to the stroke angle in actual insects.   Figure 15a shows the time history of the pitching angular displacement. As shown in this figure, the passive pitching motion is caused by the FSI periodically following the flapping motion, and the pitch angle reaches 62 • , which is comparable to the pitch angle in actual insects. Figure 15b shows the time history of the lift from a pair of the wings. As shown in this figure, the average lift is larger than the weight of the small fly (7 µN).   Figure 16 shows the wing's motion and the surrounding flow. In the left column of this figure, during the upstroke, the wing flaps from right to left. In the right column of this figure, during the downstroke, the wing flaps from left to right. As shown in these figures, the passive pitching motion similar to that observed in actual insects is given by the FSI adequately. That is, the wing keeps a high angle of attack in the middle of each half stroke, and it rotates quickly during the stroke reversals. The leading-edge vortex, which is one of most important flow structures in insect flapping wings [62], was formed clearly in each figure. These figures also show the stable attachment of this vortex on the wing. This stable attachment contributes the large translational lift in the middle of each half stroke.
Finally, the DW search was done based on the frequency response functions (FRFs). Using the FRF for the mean lift, which is shown in Figure 17a, the DW is given as the flapping frequency area ranging from 410Hz to 460Hz approximately. That is, the proposed wing using any frequency in this DW generates the sufficient lift to support the weight of the small fly (7 µN). Figure 17b shows the FRF for the stroke angle. As shown in this figure, the stroke angle is larger than 50 • approximately in the DW, which is comparable to the stroke angle in actual insects. It follows from these results that the possibility of the micro flapping wing for FWNAVs is demonstrated. the FSI adequately. That is, the wing keeps a high angle of attack in the middle of each half stroke, and it rotates quickly during the stroke reversals. The leading-edge vortex, which is one of most important flow structures in insect flapping wings [62], was formed clearly in each figure. These figures also show the stable attachment of this vortex on the wing. This stable attachment contributes the large translational lift in the middle of each half stroke. Finally, the DW search was done based on the frequency response functions (FRFs). Using the FRF for the mean lift, which is shown in Figure 17a, the DW is given as the flapping frequency area ranging from 410Hz to 460Hz approximately. That is, the proposed wing using any frequency in this DW generates the sufficient lift to support the weight of the small fly (7 μN). Figure 17b shows the FRF for the stroke angle. As shown in this figure, the stroke angle is larger than 50° approximately in the DW, which is comparable to the stroke angle in actual insects. It follows from these results that the possibility of the micro flapping wing for FWNAVs is demonstrated.

Conclusions
In this study, the computational approach for the FSI design of insect-inspired micro flapping wings was developed. This approach consists of the simulation based on the direct numerical modeling of the strongly coupled FSI, the dynamic similarity framework, and the design window search. The key of the FSI design is the accurate prediction for the 3-D strong coupling of the thin flexible flapping wing and the surrounding air flow. Here, a projection method using the algebraic splitting was used to solve this problem directly. It seems that there are few studies about computational approaches for the FSI design of FWNAVs, where a flight device actively uses the FSI to create the characteristic wing's motions to reduce the mechanical and electrical complexities of this device. Finally, the DW search was done based on the frequency response functions (FRFs). Using the FRF for the mean lift, which is shown in Figure 17a, the DW is given as the flapping frequency area ranging from 410Hz to 460Hz approximately. That is, the proposed wing using any frequency in this DW generates the sufficient lift to support the weight of the small fly (7 μN). Figure 17b shows the FRF for the stroke angle. As shown in this figure, the stroke angle is larger than 50° approximately in the DW, which is comparable to the stroke angle in actual insects. It follows from these results that the possibility of the micro flapping wing for FWNAVs is demonstrated.

Conclusions
In this study, the computational approach for the FSI design of insect-inspired micro flapping wings was developed. This approach consists of the simulation based on the direct numerical modeling of the strongly coupled FSI, the dynamic similarity framework, and the design window search. The key of the FSI design is the accurate prediction for the 3-D strong coupling of the thin flexible flapping wing and the surrounding air flow. Here, a projection method using the algebraic splitting was used to solve this problem directly. It seems that there are few studies about computational approaches for the FSI design of FWNAVs, where a flight device actively uses the FSI to create the characteristic wing's motions to reduce the mechanical and electrical complexities of this device.

Conclusions
In this study, the computational approach for the FSI design of insect-inspired micro flapping wings was developed. This approach consists of the simulation based on the direct numerical modeling of the strongly coupled FSI, the dynamic similarity framework, and the design window search. The key of the FSI design is the accurate prediction for the 3-D strong coupling of the thin flexible flapping wing and the surrounding air flow. Here, a projection method using the algebraic splitting was used to solve this problem directly. It seems that there are few studies about computational approaches for the FSI design of FWNAVs, where a flight device actively uses the FSI to create the characteristic wing's motions to reduce the mechanical and electrical complexities of this device.
In the first numerical example, it was demonstrated that the dynamic similarity law for the FSI makes two different lumped torsional flexibility models strongly coupled with the surrounding fluids dynamically similar with each other.
In the second numerical example, the effect of the flapping frequency on the passive pitching motion was systematically investigated using the dynamic similarity framework, which worked as the guide for the complexities of this FSI.
Finally, in the last numerical example, an insect-inspired micro flapping wing with the 2.5-D structure was designed using the proposed approach. The proposed wing was the engineering alterative for the nature's solution in actual insects such that it can be fabricated using the microfabrication based on the photolithography in the MEMS technology. Here, the flapping frequency was taken as the design parameter. The DW was found such that each solution in the DW can generate the lift sufficient to support the weight of small insects.
In future work, this wing will be fabricated using the micromachining. Hence, the proposed approach will contribute to the further miniaturization of FWNAVs. Data Availability Statement: The datasets analyzed during the current study are available from the corresponding author on reasonable request.

Conflicts of Interest:
The author declares no conflict of interest.