Combined co-rotational beam/shell elements for fluid–structure interaction analysis of insect-like flapping wing

Flapping wing micro-air vehicles are biologically inspired by nature flyers, specifically insects and birds. Specifically, insect wings generally consist of veins and membrane components. In this study, a structural analysis considering the vein/membrane components of an insect-like flapping wing is presented. Co-rotational (CR) finite elements are adopted in order to consider the complex wing configuration including both vein and membrane. The CR beam elements with warping degrees of freedom are employed for veins and CR shell elements for the wing membrane. The present structural analysis is verified against the analytical results obtained by an existing software, and it is validated by comparison to existing results from the literature. A fluid–structure interaction analysis is then performed. In the procedure, an aerodynamic analysis based on three-dimensional preconditioned Navier–Stokes equations is employed. Finally, a comparative study with respect to the structural characteristics is conducted. As a result, an efficiency of the present structural analysis is confirmed by comparing with the existing software. It is found that the present FSI results are in good agreement with the existing experimental and numerical results. Moreover, the passive wing twist may have a significant influence on the hover performance.


Introduction
Flapping wing micro-air vehicles (FWMAVs) are specified as being smaller than 15 cm with a maximum weight of 100 g and flying at a slow flight speed range of 10-15 m/s. Owing to their small sizes, flight regime in a low Reynolds number, and modes of operation (hover and maneuver), they require a significant scientific knowledge and engineering technology. Those vehicles Fig. 1 Cross-sectional configuration of vein (hind wing of Schistocerca gregaria) [9] are biologically inspired by nature flyers, specifically insects and birds [1][2][3].
Moreover, insect wings generally consist of veins and membrane components, and various relevant configurations are composed of assemblages of veins and membranes. Typical insect wings, composed of a membrane supported by a vein, are driven by the wing root. The wings exhibit a large amount of elastic deflection during the flapping motion. Therefore, an accurate structural dynamics modeling on insect wings will require a precise consideration of the wing kinematics and flexibility of the multi-components, i.e., the veins and wing membrane. The vein is a beam-like structure that supports the membrane. It typically has an open cross section or thin-walled section. Therefore, a warping effect in the vein may become significant when the wing is under a certain stroke. Figure 1 shows an example of the cross-sectional configuration of an insect. The membrane is an extremely thin and weak structure under a large rotational motion so that it may exhibit a geometrically nonlinear behavior. The shape of wings is complicated and may be different among the biological species. In addition, the vein arrangement is unique in each flyer [4][5][6]. Such variety in wing configuration will induce different aerodynamic performances. Moreover, resonant vibration between the natural modes of the wing structure and flapping motion is employed to produce insect wings kinematics and generate lift [7,8]. The relevant design of the wing is strongly related to understanding the fluid-structure interaction, and such resonant phenomena increase the flapping amplitude of the wings.
Wu [10] conducted experiments to investigate the influence of the type of wing configuration. In Wu's experiment, a Zimmerman planform, in which two elliptical planform contours intersected at the quarter chord, was employed. An elaborate investigation by varying the arrangement and flexibility of the vein was conducted. In Wu's experiment, a number of assembling methods of the vein/membrane were considered by varying the arrangement of the vein, e.g., parallel arrangement and radial diagonal arrangement. Moreover, the degrees of the flexibility of the vein were varied by using the number of layup (Carbon fiber). From Wu's experiment, a significant impact owing to the assembling methods of the vein/membrane and their flexibility was found to affect the wing performance. Heathcote et al. [11] conducted experiments on the effect of flexible plunging wings having an NACA0012 cross section. In Heathcote's experiment, three rectangular wings with different spanwise flexibility are employed to verify an effect on the aerodynamic performance. In their experiment, a moderate amount of wing flexibility was found to improve the wing performance.
In their experiments, a significant influence was discovered by fluid-structure interaction (FSI) phenomena due to both the assembling pattern of the vein/ membrane and the flexibility of the wing. Such FSI phenomena were found to affect the wing aerodynamic performance significantly.
Moreover, a number of numerical investigations have been conducted to investigate the relevant physics between the flexible wing structures and complex aerodynamic environment [12][13][14][15][16][17][18][19][20][21][22][23]. Chimakurthi et al. [15,16] built a numerical framework to facilitate the FSI analysis of flexible flapping wings at various fidelity levels. Their framework involved a finite volume Navier-Stokes fluid dynamics analysis and a finite element structural dynamic analysis based on geometrically nonlinear composite beams. The structural dynamic analysis was then extended into the co-rotational (CR) shell finite element, which is suitable for simulating planar flapping wings [16]. Gordnier et al. [17] developed the FSI analysis by combining a computational fluid dynamics (CFD) solution based on the higher-order governing equation and geometrically nonlinear composite beam analysis. On the other hand, we developed a three-dimensional FSI framework for FWMAVs by implicitly coupling a fluid dynamics solver with a nonlinear finite element analysis (planar and beam analysis) based on the CR element [18,19]. Its validation was achieved by considering a spanwise flexible wing and comparing with the experimental results reported by Heathcote et al. [11].
Those works [15][16][17][18][19] were focused on the interactions between the fluids and flapping elastic wings. Geometrically nonlinear composite beam or CR planar elements were developed and applied for a simple wing configuration, e.g., the NACA0012 rectangular wing.
Recently, Gogulapati et al. [20] attempted to consider both vein and membrane by using the anisotropic shell element provided in MSC.Marc; then, a comparison with the experimental result was conducted. However, a number of shell elements were employed to construct the vein, which was not quite efficient. Specifically, the vein configuration is slender. Thus, sole usage of the shell elements requires significant discretization. Masarati et al. [21] applied the nonlinear shell element for the flapping wing. In their approach, the multi-body concept was employed so that both beam and shell elements could be used for the vein and membrane at a time. However, no explicit validation results for a realistic flapping wing have been reported in the literature. More recently, Farhat and Lakshminaryan developed an FSI framework based on an embedded boundary method, and an application for a flexible flapping wing was conducted [22,23]. In their study, the wings used in Wu's experiment and the numerical analysis by Gogulapati et al. were analyzed, and those were discretized by using nonlinear shell elements.
Currently, it is found that the relevant disciplines were well developed and employed in the previous studies on FSI analysis of FWMAVs. Those attempts provided a good understanding of the physics involved in flapping wings. However, the structural analysis used in the previous studies has a limitation in considering the detailed wing configuration, i.e., the various vein arrangements. The structural characteristics of the vein and wing membrane are different. The vein for insect wings is slender and has an open cross section [5,6], whereas the membrane has a flexible thin configuration. Therefore, relevant structural components, specifically beam and shell elements, will be required for the vein and wing membrane, respectively. Moreover, such a strategy will be more efficient than that by a sole usage of nonlinear shell elements [20,22,23] requiring a number of elements to discrete slender vein.
In this study, CR finite elements will be adopted to consider the wing configurations in a precise manner. Hence, the CR beam elements with warping degrees of freedom will be employed for the veins, and the CR shell elements for the wing membrane. To have interconnection among the finite elements, a global-ized Lagrange multiplier [24] will be employed. The present structural analysis will then be validated using the results from the existing studies. In the validation, both frequency and time domain results will be considered. As a result, a precision structural analysis, exclusively developed for the flapping wing simulation, will be proposed. Then, an FSI analysis on an insect-like flapping wing will be performed by combining the three-dimensional preconditioned Navier-Stokes equations. Finally, a comparative study will be conducted to verify the effect of the passive wing twist behavior in hover.

Co-rotational finite elements
The CR formulation is one of the approaches applicable to geometrically nonlinear problems. It has been established and investigated by a number of researchers in the field. Based on the assumptions of a small degree of strain and large displacement, the CR formulation allows an accurate geometrically nonlinear structural analysis. Its main advantage is that it leads to an artificial separation between the material and geometrical nonlinearity. Therefore, a local formulation will be required for the small deformational component, and this will be performed by using the existing finite element hypothesis.
Regarding the three-dimensional CR beam element, there are a number of previous studies. The existing elemental hypothesis based on a linear strain definition, i.e., the Euler-Bernoulli or Timoshenko beam analysis, was employed. However, Battini et al. [25] established a CR beam element using the local beam element considering the Green-Lagrangian strain components to improve the accuracy. Le et al. [26,27] established a consistent dynamic formulation for the threedimensional beam element by adding a seventh degree of freedom to describe the warping of the cross section. Regarding the CR shell element, Pascote [28] suggested the shell formulation following the EICR approach. Khosravi et al. [29] extended the formulation by using a facet shell element involving the optimal triangular membrane [30] and discrete Kirchhoff triangular bending plate (OPT-DKT facet shell) [31]. A relevant dynamic formulation extended from it was suggested by Chimakurthi et al. [16].
In this study, the beam and shell element based on the CR formulation will be developed and used. Each formulation is extended to consider the nonlinear inertial terms defined in accordance with the CR framework for both the two-node beam and three-node triangular shell. An additional coordinate will be considered in the fixed frame for the rigid body motion. Thus, the existing coordinates in the CR framework will be considered as the dynamic frame maintained in the structure that undergoes motion. A relevant description regarding the CR framework will be presented in the following subsection.

Elemental kinematics
The CR formulation is established by tracking an elemental motion based on the coordinates. The coordinates, attached on the element, include the elemental fixed, undeformed, the CR, and deformed frame. The CR frame is an additional intermediate configuration between the undeformed and deformed frame. In this study, such concept is extended to consider the rigid body motion for the flapping wing structure. Thus, an additional coordinate is created in the inertial frame attached to the elemental fixed frame. For such extension to apply, a simultaneously prescribed motion is combined with the governing equation. The relevant coordinates maintaining the elemental kinematics for the beam and shell element are illustrated in Fig. 2.
In Fig. 2, R o and R r are the rotational operators from the elemental fixed frame to the undeformed frame and the deformed frame, respectively. And, R G and R L denote the rotational operators defined by the rotations with respect to the undeformed frame and the CR frame, respectively. Using the CR frame, it is possible to define the local system between the CR frame and deformed frame, and the relevant local quantities are based on the small strains [32].
To attain transformation between the coordinates, a rotational operator will be required using the parameterization of the finite three-dimensional rotations [33]. In this study, the Rodrigues formula is employed for both beam and shell elements. Thus, a consistent expression of the rotations between the segregated elements is used. The resulting expression of the rotational operator is as follows: where a tilde denotes the skew-symmetric matrix, and I is a unit matrix. In Eq. (1), θ can be obtained by θ . The spatial variation can be described by using the following relation: where δψ and δφ, respectively, denote the material and spatial angular variations, i.e., the infinitesimal rota-tions imposed onto the rotational matrix. Then, the relation between two variations can be δ R = δφ R = Rδψ, and T s can be expressed as By using the coordinates shown in Fig. 2, the local displacement of the element is given by the following relation: where r i and r i e denote the position vector of Node i with respect to the undeformed frame and the vector from the origin of the deformed frame to the deformed position of Node i, respectively. And, u i G and u i L denote the displacement of Node i with respect to the deformed frame and the undeformed frame, respectively.
The rotational operator R i G with respect to θ i G is defined by using the relation in Eq. (1). θ i G is the rotations at Node i with respect to the undeformed frame. Then, the orientation of the deformed frame can be obtained by the product R r R i L . Simultaneously, the orientation can also be obtained by the product R i G R o . With the orthogonality of the operators, the local rotational operator can be obtained as The local rotation is then evaluated through the matrix logarithm of the local operator, i.e., log R i L . Then, the off-diagonal components in log R i L are expressed in terms of the local rotation. The geometric and rigid body rotations consisting of R o and R r are defined with respect to the elemental reference frame. Let e k denote the orthonormal basis vectors on the element. Then, the orthogonal matrix defining the rigid rotations is simply given by In order to impose a prescribed motion, the joint components, i.e., the prismatic and revolute joints, are introduced. Let u p and θ p be the corresponding degrees of freedom for the prismatic and revolute joint. In this paper, an inertial frame is attached to the elemental fixed frame. When a joint is attached to Node i, the relevant motion is prescribed to the joint. Then, the motion can be considered as where R p is the additional rotational operator defined by Eq. (1) with θ p .
In the following subsection, the relevant definition of the rotational operators for the beam and shell elements will be described.

Elemental coordinates in beam element
The first coordinate axis referring to the CR frame e 1 is defined by the line connecting Nodes 1 and 2 of the element.
The remaining two axes are determined with the help of a vector b: where the vector b is directed along the local e 2 direction in the initial configuration, whereas in the deformed configuration its orientation is obtained as where R av is the rotational operator defined by Eq. (1) with

Elemental coordinates in shell element
The first coordinate axis referring to the CR frame e 1 is selected as the base of the vector parallel to the side connected between Node 1 and 2.
Using an auxiliary vector, r 13 = r 3 + u 3 G −r 1 − u 1 G , the third coordinate axis e 3 is assumed to be orthogonal to the elemental plane.
The second coordinate axis e 2 is then obtained by the following expression: e 3 × e 1 . Given this, it is possible to establish the rotational operators R o and R r .

Unified CR formulation for stiffness matrix and internal force vector
In this section, a unified CR formulation for the beam and shell elements is briefly introduced. Regardless of the type of element, the present formulation can be extended by considering the number of nodes and the corresponding nodal degrees of freedom. The details of the beam and shell formulations can be found in previous studies [25,28], and a compact form of the formulations is as follows. The virtual energy can be expressed by equating the local/global internal load vectors and the variation of the corresponding nodal displacement vectors.
It is now possible to express the local displacement vector as a function of the global displacement vector. The relationship between the variation of the local and global displacements is By equating the relationship between the local and global internal load vectors, the internal local load vector can be expressed as follows: The global tangent stiffness matrix is then defined by taking the variation of Eq. (4). The resulting form of the global tangent stiffness matrix is as follows: The transformation matrix B will be obtained by introducing the rigid body and deformational components.
where F 1 and F 2 are obtained from P T f L , E is determined using the rotational operator composed of the rigid body rotations R r and P is the projector matrix considering the deformational components. The projector matrix in nodal form can be expressed by the following equation: Details of the derivation regarding A and G can be found in Refs. [25,28]. Using the present unified formulation, it will be possible to obtain the internal quantities of the beam and shell elements. In this study, the local beam element suggested by Battini et al. [25] is employed, and OPT-DKT facet shell element [29] is adopted for local shell formulation.

Consistent nonlinear inertial formulation for beam and shell elements
The material velocity, u Tθ T T , and acceleration, ü Tθ T T , with respect to the local coordinate can be obtained by using the relation in Eq. (15). By simplification, it is possible to neglect the projector matrix [34,35]: where H denotes the interpolation function of an element. The kinetic energy of the element is expressed by considering the velocity vector.
where A ρ and I ρ are the mass per unit length (or an area for a shell element) and the spatial inertia dyadic tensor, respectively. The variation of the kinetic energy can be expressed in terms of the material velocity and acceleration: By introducing the elemental kinematics into the resulting variation form of Eq. (22), the inertial load vector is obtained as follows: where f K L is the inertial load vector related to the local coordinate and it can be expressed as Using the linearization of the inertial load vector, the mass and gyroscopic matrices will be obtained. The compact form of the matrices is as follows: For the beam element, interdependent interpolation element formulations are employed for the beam shape functions, and the detailed expressions are described in [26,27].

Interconnecting strategy and solution methodology
In this section, an interconnecting approach between the present beam and shell elements will be presented. Then, a solution methodology for the resulting governing equation will be described. Figure 3 shows a schematic of the present interconnecting approach for the wing composed of the multicomponents. To enforce the kinematic constraints among the components, a globalized Lagrange multiplier method is employed [24,36]. This method will yield a straightforward formulation, in which the addition of the virtual work of individual components is performed to obtain the total virtual work of the entire body. In this procedure, each of the kinematic constraints between any two bodies or two adjacent nodes can be considered using the Lagrange multiplier (Fig. 3b). The N constraint equations between the beam and shell elements can be expressed using the translational and rotational degrees of freedom for adjacent nodes.

Global Lagrange multiplier technique
To solve the additional equation, a Lagrange multiplier is introduced. The variational formulation of the total energy for the multi-components can be expressed as where λ is the Lagrange multiplier for the constraint between the components. The superscripts B and S denote the physical quantities regarding the beam and shell components, respectively. By introducing K mb , the resulting stiffness matrix for the multi-components can be written in the following compact form: To solve the above equation, the displacements can be obtained: For the Lagrange multiplier vector, the following equation can be derived from the constraint equations.
When the interconnected degrees of freedom are increased, the system matrix in Eq. (29) becomes sparse. In order to solve the system of equations, a sparse linear solver (PARDISO) [37] is employed. By taking PARDISO into account, it is possible to handle the sparsity of the system matrix.

Governing equation
In this study, the generalized-α method is employed to solve the resulting nonlinear equation. The resulting algorithm to solve the nonlinear equation is a prediction/corrector type, similar to those in the previous studies [38,39]. For nonlinear system, the resulting equation can be expressed by: where F mb G denotes the unbalanced load vector, where f mb e is the external load vector of the multicomponents.
In the generalized-α method, displacements and velocities are updated by introducing an algorithmic accelerations a, and they are obtained in Newmark integration formulae. q mb G,n+1 = q mb G,n + hq mb G,n+1 + h 2 1 2 − β a n + βa n+1 (34) q mb G,n+1 =q mb G,n + h (1 − γ ) a n + γ a n+1 (35) where the subscript n indicates an index of the time step, and h denotes the size of the time step. γ and β can be selected in order to have suitable accuracy and stability properties. a is defined by the following recurrence relation: Here, α m and α f are computed as: where the spectral radius ρ ∞ ∈ [0, 1]. And, β and γ are defined by: For the multi-components systems, the linearized equations of motion can be expressed as: where coefficients c 1 , and c 2 are defined by the following relation: Moreover, s is a scaling factor for the Lagrange multipliers, and it is chosen to be (h 2 k r +hc r +m r )/ h 2 . k r , c r and m r are characteristic stiffness, damping and stiffness coefficients of the system. In Eq. (39), the system matrix for the differential algebraic equations (DAEs) is ill-conditioned, and it becomes almost singular when employing small time steps. Hence, an optimal conditioning of the system matrix proposed by Bottasso et al. [40] can be applied in order to avoid the ill-conditioned situation. For correction step, the displacements, velocities and accelerations are updated by the following relations 4 Fluid-structure interaction framework

Fluid analysis
Three-dimensional preconditioned Navier-Stokes equations are chosen as governing equations for the simulation of the flapping wing at a low Mach number flow region. The integral form of the nondimensionalized governing equations with free stream conditions is written as where W is the conservative solution vector, Q p is the primitive solution vector, F i is the inviscid flux vector and F v is the viscous flux vector. Note that a fictitious time derivative term is added for accurate unsteady computations. For accurate and efficient computations on low Mach number flows, the preconditioning matrix of Weiss and Smith [41] was chosen. The dual-time stepping method in conjunction with the approximate factorization-alternate direction implicit method is used for the time derivative term of the governing equations. t and τ denote physical time and pseudo-time used in the time-marching procedure. Roe's approximate Riemann solver [42] and the central difference method are used to discretize the inviscid flux vector and the viscous flux vector, respectively. Details of the flow analysis are included in [43,44]. In the present aerodynamic analysis, a radial basis function (RBF) is employed to consider the structural deformations. The geometric conservation law is employed as well to alleviate the problem in the volume computation of the deforming grid.

Coupling methodology
To couple the structural with the aerodynamic analysis, an implicit coupling approach is applied. Owing to the thin configuration and large amount of displacement in the biomimetic wing, an explicit coupling between the fluid at low Reynolds number and thin structure can induce numerical instability, i.e., an added-mass effect. The added-mass effect is the instability depending upon the densities of fluid and structure and also on the geometry of the domain [45]. In the implicit coupling approach, both aerodynamic and structural solutions are determined iteratively by exchanging the results every coupled time step. A diagram of the present implicit coupling approach is presented in Fig. 4. The sub-iteration in the present FSI analysis is performed in accordance with the increments of pseudotime τ [Eq. (42)] in the fluid analysis. In this study, the parameter n f si, which controls the number of result exchanges, is introduced to improve the computational efficiency. The value of this parameter is the period of the result exchange. For example, if n f si is 3, exchange of the result will be performed every 3 sub-iterations. If n f si is 1, exchange of the result will be performed at every sub-iteration, which is the same as that of the general implicit coupling algorithm. The present coupling approach is further improved from the previous approach presented in [18,19]. In the CFD analysis,  In the case of an irregular configuration, i.e., the Zimmerman wing or realistic insect-like wing, the linear interpolation is limited owing to the nonuniform grid topology. Hence, the RBF interpolation is employed for the analysis involving a complex configuration. Figure 5 shows the present approach for interface between fluid and structure. When considering the multi-components structural analysis (Ω S = Ω B S ∪ Ω S S ), the pressures p f obtained by CFD are interpolated to the pressures p s acting on the shells (Ω S S ) that share the same boundary with CFD surface grids (Ω F ), S S = F , while enforcing the force conservation between the fluid and structure. Then, p s is transformed to the nodal force vector f S e for the shell components with respect to the elemental reference frame. The pressure is a follower load so that the relevant rela-tion to f S e can be expressed by using the relation in Eq. (15).
where H is the interpolation matrix regarding the shell element, and e 3 = 0, 0, 1 T . Thus, the external force vector f mb e in Eq. (33) is expressed as: Similarly, the displacement results of the shell components are transferred to the CFD surface grids.

Time-transient analysis of a beam/shell combined structure
In this section, the present structural analysis is verified in terms of the accuracy and efficiency. Time-transient analysis using the multi-components analysis is conducted. The results are compared to that by the existing software, and the efficiency is assessed by taking the computational time into account. In order to verify an accuracy of the present multi-component analysis, one of the existing multi-body dynamic analysis software, DYMORE [46], is used and compared. The example used in this section is depicted in Fig. 6. The distinguished flexible components consisting of beam and plate are considered. Such structure is under the prescribed flapping motion. In order to consider the prescribed flapping motion, a revolute joint component is connected to the beam root (Station R) and the motion is prescribed as − 20 circ × (1 − cos(10π t)) at Station R. In DYMORE, a quadrilateral shell element, different with the present shell analysis, is only available. Thus, additional reference solution is obtained by using ANSYS. In ANSYS prediction, the structure is discretized by using a number of three-dimensional solid elements in order to obtain the reference solution. And the same number of nodes is utilized in the present and DYMORE predictions to compare the efficiency in terms of the computational time. Table 1 shows the summary of the relevant information of discretization used in each analysis. Comparison of the displacements at Stations A is illustrated in Fig. 7. Figure 7b shows close view of   (39). Adequate handling sparsity of the resulting differential algebraic equations exhibits decrease in the computational time in the present analysis.

Validation with experimental results
In this section, the present structural and FSI analyses will be assessed by comparing the results with those of the experiment by Wu [10]. In Wu's experiment, the Zimmerman planform was employed, and a single The relevant configuration of the wings is illustrated in Fig. 8. Figure 9 shows the sketches presenting the geometrical information of the veins/wing membrane and the scheme to vary the vein components. In Fig. 9b, the lines represent the number of layers for layup-leading edge varies from 3 to 1 layers of unidirectional carbon fiber and battens vary from 1 to 2 layers. And each carbon layer is 0.8 mm wide and 0.1 mm thick. The resulting wings are named Li B j, where L stands for leading edge, i for the number of carbon fiber layups in the leading edge, B for batten along the chordwise direction and j for the number of carbon fiber layups in the batten.
To analyze the wings in the aforementioned experiment, seven flexible components are employed, that is, six-beam and one-shell components. The relevant configuration used in the present analysis is illustrated in Fig. 10. The present analysis employs a total of 6882 degrees of freedom including the flexible beam/shell and the Lagrange multipliers. The mechanical properties of these wings are described in [20]. Based on those in [20], the properties used in the present analysis are assumed to be an isotropic materials. The relevant information is summarized in Table 3.

Structural analysis
The first validation of the present analysis will be accomplished by comparing the natural frequencies and mode shapes measured in vacuum with those in the experiment by Wu [10]. In the modal analysis, all the nodes in the region of the rigid triangular base (Station R in Fig. 10) are set to be the fixed boundary condition.
Comparison of the natural frequencies shows that the present results are in good correlation, and the relevant 214 H. Cho et al.   Fig. 11. Moreover, a qualitative comparison on the mode shapes is shown in Fig. 12. As a result, validity of the wing modeling by the present structural analysis is ascertained. Moreover, it is found that the bending modes show a significant coupling with the torsion. An increase in wing rigidity will be more influential when the stiffness of the leadingedge batten is increased further than that of the chordwise battens will do.

FSI analysis
To validate the present FSI analysis, Wu's experiment is considered. The wings with the semi span length (R S )  Table 4. Flapping frequency f o is varied from 5 to 40 Hz, and the flapping amplitude And for the present fluid analysis, C-H-type multi-block grid topology is employed, and a total of 0.7 million grid points (3200 cells on the wing surface) are used from grid convergence test (Fig. 13). For structural analysis, the same geometry and discretization of the wing presented in Fig. 10 are used and a fixed boundary condition is assigned to the rigid triangular base (Station R). Along X -axis with respect to Point O, the flapping motion is prescribed at the revolute joint, connected to the rigid base. The relevant prescribed motion is set to be β o × sin(2π f o × t). In the present FSI analysis, a time step is chosen to be 1/(2000 × f o ).  In the present FSI analysis, the periodic results are obtained after five periods are elapsed. In Wu's experiment, the increase in the aerodynamic coefficient, i.e., either C T or C L , was not acquired. However, the wing tip deflection and average thrust values with respect to the flapping frequency were included. In this study, Fig. 14 Measurement of the wing twist in Wu's experiment [10] L2B1 and L3B1 wings under the flapping motion for the frequencies from 15 Hz to 30 Hz and a flapping amplitude of 35 • are considered. Additionally, the wing twist increase is computed by following the procedure in Wu's experiment (Fig. 14). The first validation of the present analysis is conducted by using the wing tip deflection increase under a flapping frequency of 25 Hz (U ref = 4.58 m/s, Re = 7304). In this comparison, the previous numerical prediction obtained by Lakshminaryan and Farhat [23] is used. The relevant comparison of the tip position is illustrated in Fig. 15. The wing twist increase at 83% position of the wing span for L2B1 and L3B1 wings is compared, as shown in Fig. 16. Then, the average value of the present results for L2B1 and L3B1 wings is compared with the experimental observation. The relevant comparison from 15 Hz to 30 Hz is shown in Fig. 17.
The present results show a satisfactory correlation with those observed in the experiment. Moreover, the present analysis shows a good agreement with the result predicted by Lakshminaryan and Farhat [23]. In their computation, nonlinear shell elements with 13,000 degrees of freedom were employed. By contrast, beam and shell elements with 6882 degrees of freedom are used in the present analysis. Thus, it may be concluded that the present structural analysis is more efficient for a similar accuracy. Specifically, the leading edge and chordwise batten were discretized by shell elements in the previous numerical study. However, the discretization using beam elements results in the present improvement in terms of computational efficiency.
From the above examination, it is found that the present analysis shows a good agreement with the previous numerical results and experimental observations. Furthermore, L3B1 wing exhibits an elastic twist owing to the geometrically anisotropic feature. Such elastic twist is more significant in L2B1 wing. Regarding the aerodynamic phenomena, the flapping wing is operated in the fluid at the low Reynolds number and it is known that the viscous effect will change a formation of vortex structure [47]. When the wing is in the hover, the Reynolds number is varied by the flapping frequency, and it will result in different aerodynamic environment due to the wing elastic/inertial forces and varied viscous effects. Such complex phenomena are not fully investigated in the present analysis. Moreover, in L2B1 wing, the increasing rate of the averaged thrust is degraded when the flapping frequency is greater than 25Hz. Such situation may be related to the resonant effects. The thrust is 2 × f o excitation and the first bending mode of L2B1 wing is 43.1Hz, quite close to 50Hz. However, extensive investigation regarding the resonant effect is still required. In the next section, a comparative study among six different wings will be conducted to examine the influence of the structural behavior on the wing performance at the same flapping frequency.

Effect of flexibility on aerodynamic performance
In this subsection, four types of wing structures, i.e., a rigid wing and three different flexible wings, are used to compare the effect of the structural behavior on the wing performance. Regarding the flexible wings, a single-component (shell) wing and four multicomponent (beam and shell) wings are considered. To interpret the influence of the wing structures, the coordinate of the aerodynamic loads is established with respect to the operating situation of the flapping wing MAV (Fig. 18), and a pair of wing is considered. Thus, the thrust and lift are the direction of the vertical and sideward force, respectively. As shown in Fig. 18, the force acting along the spanwise axis will be canceled. Hence, only the lift is evaluated as the sideward force.
For structural analysis, same geometry and discretization of the wing presented in Fig. 10 are used and a fixed boundary condition is assigned to the rigid triangular base (Station R). Along the X -axis, the flapping motion is prescribed at the revolute joint. And, the joint is connected to the rigid base at Point O. The relevant prescribed motion is set to be β o × cos(2π f o × t). And, the flapping frequency f o is chosen to be 25 Hz. The single-component wing is a thin aluminum plate, and the multi-component wings are L2B1, L2C B1, L3B1  and L3C B1 Zimmerman wings, respectively. L2C B1 and L3C B1 wings are designed in the present analysis, and the curved batten is incorporated inside the wing membrane. The relevant configuration is shown in Fig. 19. For multi-component wings, the structural properties are the same as those in Table 3. Moreover, thin aluminum plate (E = 70 GPa, ν = 0.33, and ρ = 2700 kg/m 3 ) with thickness, 0.27 mm, is considered for isotropic wing. Thus, the first bending fre-quency is 67.85 Hz, similar to that of L3B1 wing. The operating condition is the same as that in Table 4.
First, the variation of aerodynamic forces is evaluated. Figure 20a shows a comparison of the thrust history, while the sideward force history is illustrated in Fig. 20b. Moreover, the average values of the thrust and sideward force are summarized in Table 5. As shown in the comparison of the sideward force history, the peak-to-peak value is significant for the L2B1 wing.   Figure 21 shows a comparison of the tip displacement and wing elastic twist. Furthermore, the peakto-peak value of the tip displacement and wing elastic twist are summarized in Table 6. It is found that a similar amount of tip bending deflection is predicted for both L2B1, L2C B1 and thin aluminum plate wings. As shown in Fig. 21b, all the flexible wings exhibit an elastic twist behavior. The elastic twist in the thin aluminum plate wing is induced by its unsymmetrical configuration and boundary condition of the wing. For the multi-component wings, the geometrically anisotropic feature induces a significant amount of elastic twist. However, L2C B1 and L3C B1 wings show decreased elastic twist when compared to those by the L2B1 and L3B1 wings. Such situation indicates that the batten arrangement upon the wing planform influence on the wing elastic twist, significantly. On the other hand, the rigidity of leading edge vein is related to the wing bending deflection.
Moreover, it is found that the wing elastic twist may be a crucial factor in terms of the thrust performance in hover. To interpret the details, pressure and vorticity  contours at 0%, 25%, 50% and 75% spanwise position in 0.25 t/T are analyzed. Figure 22 shows comparison of pressure coefficient contours among six different wings. Regarding the flexible wings, a lower pressure distribution is dominant near the wing surface, and a larger pressure difference is induced between the upper and lower region. This indicates an increase in the net aerodynamic force, which may be induced by the bending deflection. In the multi-component wings, the large low-pressure region is not sustained near the wing but it moves downward (from the leading edge to trailing edge). Such a different pressure distribution near the wing surface may be induced by the wing elastic twist, and those may result in an increase in the averaged thrust. Figure 23 shows comparison of the vorticity contours among six different wings. For the multi-component wings, a reattachment of the leading-edge vortex (LEV) along the chordwise location is clearly   shown when compared with those of the thin aluminum plate and rigid wings. Specifically, the strong LEV and its reattachment are predicted in the L2B1 wing. And such LEV reattachment follows an amount of the wing elastic twist. As a result, the elastic twist behavior reduces the effective angle of attack, and induces the LEV reattachment. Also, the attached vorticity moves from the leading edge to trailing edge. Thus, the significant difference in the pressure between the leading and trailing edge induces an increase in thrust.

Conclusion
In this paper, an improved computational approach for an insect-like flapping wing is presented. To consider the insect-like wing configuration more comprehensively, CR finite elements are adopted. Then, CR beam elements are employed for veins, and CR shell elements for the wing membrane. The constraints between the beam and shell element are established by using the global Lagrange multiplier technique. An implicit coupling scheme is employed to combine the structural and fluid analyses. Three-dimensional preconditioned Navier-Stokes equations are employed for the fluid analysis. First, the present structural analysis is verified by comparing with that obtained by DYMORE and ANSYS. Then, the present FSI analysis is validated by comparing with the results obtained from the experiment. In the structural analysis, the accurate time-transient analysis is verified when compared with ANSYS results, and 92% reduced computational time is required when compared to that by DYMORE. The present FSI results show a good correlation with the experimental observation. Moreover, an explicit investigation on the influence of the wing structural feature is conducted by using the presently improved FSI analysis. As a result, the following accomplishments are made in this study: -A nonlinear dynamic analysis for two-node beam and three-node triangular shell elements is developed, and the analysis is extended to the multi-body dynamic approach. -Verification of the present structural analysis is conducted, and its efficiency is estimated. -The FSI framework is validated by considering the Zimmerman wing experiment, and the present results show slightly better agreement than that by the previous study. -The present FSI analysis examines the effect of different wing structures and investigates the relevant influence of the flapping wing performance in hover.
Finally, the influence of the elastic twist behavior is investigated by employing a coordinate based on the FWMAV. It is found that the geometrically anisotropic feature of the multi-component wing structure induces a large amount of elastic twist, which is beneficial to generate a vertical force. Also, the leading edge rigidity mainly influences the bending deflection, and the arrangement of vein upon the wing planform may affect the wing twist behavior. Thus, the relevant arrangement or rigidity of the vein on the flapping wing should be considered to enhance the wing performance in hover.
In the future, more extensive comparative study will be conducted and optimization framework for the flapping wing will be developed.
the source, provide a link to the Creative Commons license, and indicate if changes were made.