Computational evaluation of flight performance of flapping-wing nano air vehicles using hierarchical coupling of nonlinear dynamic and fluid-structure interaction analyses

This study endeavours to introduce an inventive hierarchical coupling methodology for evaluating the flight capabilities of polymer micromachined flapping-wing nano air vehicles (FWNAVs) using advanced computational tools. Furthermore, our objective is to provide a practical demonstration of FWNAVs in both tethered and airborne scenarios. These dual objectives represent the distinctive and pioneering aspects of this research. Such FWNAVs, which are insect-inspired robots, have a 2.5-dimensional structure. An FWNAV comprises a micro piezoelectric drive system and a pair of micro thin flexible wings. The drive system includes a micro transmission and a piezoelectric bimorph actuator. The flight performance of the designed FWNAV is evaluated using a hierarchical coupling approach, where the whole system is decomposed into a micro wing and a micro piezoelectric drive system. The coupling between these parts is modelled as one-way coupling and that between the micro wings and the surrounding air is modelled as strong coupling. In the one-way coupling analysis, nonlinear structural dynamic analysis is conducted for the micro piezoelectric drive system, where the dynamic response is transmitted to the micro wings via the Dirichlet boundary condition. In the strong coupling analysis, strongly coupled fluid-structure interaction analysis is conducted for the micro wings and the surrounding air to consider their strong coupling. The optimisation of flight performance is conducted using fluid-structure interaction analysis to achieve sufficient lift force to support the weight of an FWNAV. The design of a tethered and flyable FWNAV with a size of 10 mm or smaller is demonstrated. This FWNAV can be easily fabricated using polymer micromachining.


Introduction
Rotor-wing and fixed-wing micro/nano air vehicles have limited applicability for indoor inspection (Galinski, 2005) because of their relatively poor maneuverability and hovering characteristics.Other types of flight systems that have good flight characteristics are thus required and so Flapping-wing air vehicles have become increasingly popular.Depending on the characteristics of the flight system, flapping-wing air vehicles are either insect-inspired or bird-inspired.Insects can perform vertical take-off and hover more easily than can birds and hence have better maneuverability and hovering efficiency (Petricca et al., 2011).Depending on size, insectinspired flapping-wing air vehicles can be divided into flapping-wing micro air vehicles and flapping-wing nano air vehicles (FWNAVs) (Hylton et al., 2012).
There are two main challenges in developing FWNAVs at the insect scale, namely design problems due to multisystem and multiphysics coupling (e.g.electroaeromechanical design) (Ishihara et al., 2017;Ramegowda et al., 2020) and manufacturing problems due to the very small structures required to realise the complicated motion of the wings and the assembly of micro size components.Another challenge is the production of a sufficient lift force to support the weight of an FWNAV due to the complex fluid mechanics at low Reynolds numbers (Re < 1000) (Liu & Aono, 2009).The design challenges can be overcome using either computational biomimetics (Ishihara, 2022;Ishihara, Liu et al., 2022a;Ishihara, Takei et al., 2022b), which can be used to define strong multisystem and multiphysics coupling of FWNAV or the design window search methodology (Ishihara et al., 2002a).The manufacturing challenges can be overcome using microelectromechanical system technology (Ishihara et al., 2020).A sufficient lift force can be produced by optimising flight performance.
The flight performance of FWNAVs can be evaluated experimentally or computationally.Due to the coupled physics between air and the FWNAV structure, a computational evaluation is preferred.To investigate the coupled physics between flexible wing structures and an incompressible fluid and accurately predict flight performance, fluid-structure interaction (FSI) analyses, a type of numerical simulation, have been conducted (Ishihara & Yoshimura, 2005;Ishihara et al., 2002b;Liu et al., 2010;Rubio & Chakravarty, 2017).
In Figure 1, we present an overview of our FWNAV concept, which comprises several critical subsystems, including a micro piezoelectric drive system (consisting of a transmission and piezoelectric bimorph actuators), a pair of micro wings engineered to generate significant aerodynamic forces, and a robust support frame.This FWNAV inherently involves intricate multisystem interactions and multiphysics coupling, a common challenge in complex engineering problems (Felippa et al., 2001; Ishihara et al., 2015;Ramegowda et al., 2020).To address this complexity, we introduce a novel hierarchical coupling approach specifically tailored to address the coupled dynamics of FWNAVs in this study.
The primary objective of this study is to unveil a pioneering hierarchical coupling approach designed to evaluate the flight capabilities of polymer micromachined FWNAVs through the application of advanced computational tools.Additionally, our study sets out to provide tangible demonstrations of FWNAVs in both tethered and airborne conditions.It's worth noting that these dual objectives represent the unique and innovative aspects that distinguish this research.
Our approach handles the multiphysics coupling of FWNAVs through a hierarchical method.Flight performance evaluation involves partitioning the analysis into two subsystems: one-way coupling models the micro piezoelectric drive system, and strong coupling addresses interactions between the micro wing and the surrounding air.
The output from the one-way coupling analysis, which includes the response from nonlinear structural dynamics with the micro piezoelectric drive system and micro wing, serves as input for the strong coupling analysis, specifically the FSI of the micro wing and surrounding air.
Optimising flight performance via FSI analysis ensures sufficient lift force to support an FWNAV's weight.Ultimately, we demonstrate the design of a tethered and flyable FWNAV with a size of 10 mm or smaller, a highly feasible fabrication using polymer micromachining techniques.
In the pursuit of this study, we have encountered significant challenges.These encompass the intricate task of mesh generation for both fluid and structural components in the context of FSI analysis, the necessity to maintain the feathering angle within a predefined threshold to ensure the successful completion of FSI simulations, and the complex endeavour of generating sufficient lift force to counterbalance the weight of the FWNAV, all while carefully managing Reynolds numbers.

Hierarchical coupling approach
Our FWNAV consists of a piezoelectric drive (bimorph actuator and micro transmission) and a pair of micro wings.Similar to an insect flying in the sky, the surrounding air comes into effect during flight, making the FWNAV a complicated multiphysics system.In this study, the multiphysics coupling of the FWNAV with the surrounding air is handled using a proposed hierarchal coupling approach based on the hierarchical decomposition.
Figure 2 shows the hierarchical decomposition of this coupled problem.As shown, the hovering FWNAV is first partitioned into two subsystems, namely the piezoelectric drive and the micro wings with surrounding airflow.Each subsystem consists of each coupled problem i.e. a micro piezoelectric drive system is itself a coupled problem (electro-mechanical) and micro wing with air is also a coupled problem (aero-mechanical), and these subsystems are further partitioned, as shown in Figure 2. In  this study, the first partitioning step, highlighted by the dashed line, is used to evaluate the flight performance of the FWNAV.
In this partitioning step, the whole system is decomposed into the micro wing and the micro piezoelectric drive system.Their coupling is modelled as one-way coupling.The coupling between the micro wing and the surrounding air is modelled as strong coupling.
In the one-way coupling analysis, nonlinear structural dynamic analysis is conducted for the micro piezoelectric drive system, where the dynamic response is transmitted to the micro wing via the Dirichlet boundary condition.In the strong coupling analysis, strongly coupled FSI analysis is conducted for the micro wing and the surrounding air to consider their strong coupling.
Figure 3 shows the details of the flight performance evaluation method for the FWNAV.As shown, flight performance is evaluated using two computational methods, namely a nonlinear structural dynamic analysis of the FWNAV and an FSI analysis of the micro wing.The former provides the flapping amplitude Φ 0 as an output of the nonlinear dynamic behaviour of the FWNAV.The sinusoidal flapping angular displacement Φ = Φ 0 sin2π ft is used as an input for the FSI analysis of the flapping micro wing.

Evaluation and optimization of flight performance of FWNAV Using hierarchical coupling approach
The previous section described a hierarchal coupling approach for obtaining the flight performance of the FWNAV.This approach is now discussed in detail.
The flight performance of FWNAVs is evaluated (see Section 3) and optimised (see Section 4) using a hierarchal approach, as shown in Figure 4.As shown, the output of the nonlinear structural dynamic analysis of the FWNAV is used as the input variable for the FSI analysis of the micro wing and the output (numerical lift force) of the FSI analysis of the micro wing is used as the input variable (lift coefficient) in the calculation of the theoretical lift force for the nonlinear dynamic analysis of the FWNAV.The outputs (wing motion and lift forces) of the two analyses (nonlinear dynamic analysis of the FWNAV and FSI analysis of the micro wing) are then compared (see Section 3.3).
This comparison validates that the FSI analysis of the micro wing is sufficiently accurate for determining the flight performance of FWNAV and also diluting the computational complexity associated with the FSI analysis of FWNAV design if the latter would have been considered for determining the flight performance of FWNAV.Hence, the FSI analysis of the micro wing is used to optimise the flight performance (i.e.maximise the lift force) of a tethered and flyable FWNAV.

Nonlinear dynamic analysis
To consider geometrical nonlinearity, which is the basic characteristic of our FWNAV, a nonlinear finite element transient analysis is simulated using Marc solver, where the basic nonlinear dynamic equation, given below (Bathe, 2006), is solved.
where a, v, and u are the acceleration, velocity, and displacement vectors, respectively; M and C are the mass and diffusive matrices, respectively; and Q and F are the interior and external force vectors, respectively.The structural damping is negligible compared with the surrounding air damping.Hence, C is set to 0, and Equation ( 1) is reduced to The above equation is solved using the Newmark integration method in the Marc solver, where the integration parameters α and β are taken as 1/2 and 1/4, respectively.

Fluid-Structure interaction analysis
The multiphysics of a flexible micro wing flapping in the airflow can be accurately predicted using a threedimensional strongly coupled FSI analysis.The FSI analysis was carried out using a self-developed computer programme, which was validated in previous studies (Ishihara, 2018;Ishihara et al., 2014;Ishihara & Horie, 2014).The main equations and the finite element modelling related to FSI are described below.

Governing equation for FSI
Based on the total Lagrangian method, micro-wing motion can be formulated as follows: where ρ is the mass density, d/dt is the Lagrangian time derivative operator, v i is the i-th component of the velocity vector, and σ ij is the ij-th component of the Cauchy stress tensor.Note that i corresponds to a Cartesian coordinate direction (x, y, or z), and the characteristic quantity of the micro wing is denoted by the superscript w.
In FSI, fluid motion with a moving boundary can be described using the incompressible Navier-Stokes equations and the arbitrary Lagrangian-Eulerian (ALE) method, which can be stated as follows: where ∂/∂t is the ALE time derivative operator.Note that the characteristic quantity for the ALE coordinates is denoted by the superscript m and that for air is denoted by the superscript a.
In FSI, the interaction between the micro wing and surrounding airflow at the interface can be formulated as the following coupling conditions: where n i is the i-th component of the outward unit normal vectors on the interface between the micro wings and the surrounding air.

Monolithic equation system for fluid-structure interaction
Equations ( 3)-( 7) are discretized using the finite element method and combined at the fluid-structure interface in a matrix-vector form.These combined equations are called the monolithic equation system for FSI, which can be written as ) where M, C, and G are the mass, diffusive, and divergence operator matrices, respectively, N is the convective term, a, v, and u are the acceleration, velocity, and displacement, respectively, p is the pressure vector, and q and g are the internal and external forces, respectively.Note that matrix lumping and the matrix transpose are denoted by subscripts L and τ , respectively.Equations ( 8) and ( 9) can be further simplified using the linearisation method: where M * is the generalised mass matrix, which consists of mass and tangential stiffness matrices, represents an increment, t is time, G is the pressure stabilisation term due to PSPG (pressure-stabilizing/Petrov-Galerkin) contribution (Tezduyar et al., 1992), g and h are the residual vectors of Equations ( 8) and ( 9), respectively, and a is the acceleration increment, which can be evaluated using Newmark's method, where u = β t 2 a and v = γ t a.Here, time integration is performed using a predictor-multicorrector algorithm and the pressure and elastic interior forces are evaluated using an implicit method.

Projection method using algebraic splitting
While solving the simplified monolithic system of equations given by Equations ( 10) and ( 11) using the monolithic method at the fluid-structure interface to satisfy the interface conditions ( 4) and ( 5), an ill-conditioned equation system is obtained.A projection method based on algebraic splitting (Ishihara & Horie, 2014) is derived to solve this ill-conditioned equation system.According to this projection method, Equation ( 8) can be further simplified to evaluate the state variable for the known pressure using the linearisation method, which can be expressed as follows: where â is the intermediate acceleration.
Subtracting Equation ( 12) from Equation ( 10) gives where v is the intermediate velocity.
Multiplying both sides of Equation ( 13) by τ G L M −1 yields the desired result: where M * = M * -L M. Equation ( 14) can be divided into two parts, namely the pressure Poisson equation given in Equation ( 15) and the reduced equation given in Equation ( 16).
If the pressure Poisson equation is solved using convergent coupling iteration, the incompressibility constraint given by Equation (9) (i.e.Equation ( 16)) will be satisfied.
The following steps are used to solve the monolithic equation system: (1) The state variables are predicted by solving Equation ( 12). ( 2) The pressure is obtained by solving Equation ( 15).
These steps are repeated until convergence is achieved.

Parallel computation environment
As mentioned, Equations ( 10), ( 12), and ( 15) are iteratively solved until convergence to evaluate the matrixvector products.Hence, the computational cost is very high.It can be reduced by parallel computation based on the mesh decomposition method (Ishihara & Horie, 2014).
Figure 5 shows a schematic diagram of the parallel computation.In the parallel computation, since the analysis degrees of freedom of the fluid mesh are very high compared to those of the micro wing, the fluid mesh (air mesh a ) is decomposed into N submeshes ( a1 , a2 , . . .aN ).To avoid the computational control difficulty at the interface, micro wing mesh w is kept inside air mesh a1 , as shown in Figure 5.The parallel computation is divided into two steps.
Step 1: At computational node P i (i = 1,2, 3 . . ., N d ), the matrix-vector product is evaluated using the elementby-element method: where A denotes the global matrix, p denotes the global vector, and q denotes the matrix-vector product.Note that the element number is denoted by e and the element quantity is denoted by the superscript e.
Step 2: In Equation ( 17), at the interface between ai and aj (j = i), the matrix-vector product computed at P j is transferred to P i .The products at P j and P i are added to complete the corresponding nodal data computed in step 1.
The parallel computation was executed on a computer with two 10-core Xeon 2.8-GHz CPUs and 32 GB of memory (Ishihara, 2022).Note that the Newton-Raphson method has been applied for the convergence criterion.

Finite element modelling of micro wings and surrounding fluid domain
Linear equal-order-interpolation velocity-pressure elements (Tezduyar et al., 1992) are used for the fluid (air) and fluid-structure interface (Zhang & Hisada, 2001) and mixed interpolation of tensorial components (MITC) shell elements (Bathe & Dvorkin, 1985;Noguchi & Hisada, 1993) are used for the micro wings.In addition, a pair of overlapping air nodes, per one wing node, is used at the interface to describe the pressure discontinuity between the front and back surfaces of the micro wing.
The coupling conditions given by Equations ( 6) and ( 7) are rewritten in matrix form as where T is the interpolation matrix and Q is the internal force vector that includes all effects.Note that the coupled degrees of freedom are denoted by the subscript c.The above equations are rewritten in nodal form as where v, Q, and g denote the nodal components of the velocity and the internal and external force vectors, respectively.Note, the front and back sides of the micro wing are denoted by the subscripts f and b, respectively.

Computational evaluation of flight performance of FWNAVs
In this section, the flight performance of two previously designed FWNAVs is evaluated using the proposed hierarchal coupling approach.In our previous study (Rashmikant, 2022;Rashmikant & Ishihara, 2023), the iterative design window search method was used to design two FWNAVs with different sizes, namely a first prototype and its miniaturised design, as shown in Figure 6.The plan and sectional views of the two FWNAV designs are the same.The length of the support frame L s = 26.9mm and the hinge gap h = 160 μm for the first prototype; L s = 10 mm and h = 240 μm for the miniaturised design.

Nonlinear structural dynamic analysis of FWNAVs
The nonlinear structural dynamic analysis of the FWNAVs was carried out using the Marc solver.The corresponding material distribution and problem setup are shown in Figures 7 and 8, respectively.As shown in Figure 7, two polyimide materials were used for the analysis and fabrication.Their properties are shown in Table 1.As shown in Figure 8, an alternating force u x = U x × sin2π f t (U x = 81 μm is the alternating force amplitude and f is the flapping frequency), given by the force balance between the actuator and transmission (Rashmikant et al., 2021a), is applied to the upper part of the FWNAV, where the free end of the piezoelectric bimorph actuator is attached.The lower part of the FWNAV is fixed where the fixed end of the actuator is attached.In the problem's initial configuration, the flapping frequency, denoted as 'f,' is chosen for a FWNAV with meticulous consideration.The primary objective behind this selection is to ensure that the elastic hinges  involved in the FWNAV's flapping mechanism do not adhere or stick together during the flapping motion.A comprehensive explanation of this crucial design requirement has been thoroughly examined in a prior research conducted by us (Rashmikant & Ishihara, 2023).It is worth noting that the problem setup in this earlier study exhibited some variations, with the main distinction being the absence of a wing membrane in the FWNAV being investigated.Thus, the flapping frequency for the first prototype is f = 50 Hz and that for the miniaturised design is f = 100 Hz.The nonlinear structural dynamic analysis of the FWNAVs was conducted using the Marc solver with a 20-node solid hexahedral element.The detailed finite element mesh setup is given in Table 2.The analysis of the structural aspects is presented in relation to the flapping performance characteristics, which have been scrutinised for the FWNAVs.The flapping angular displacement was measured at two locations, namely the leading edge and the wing attachment.The leading edge has a large bend (see Figure 9) and thus there is a significant difference in the stroke angle between these two locations.The flapping velocity is used to calculate lift forces.The time histories of the flapping velocity for the two designs are shown in Figure 11; they were computed at the position of r 2 from the wing base (see Figure 8).As shown, for the first prototype, the mean flapping velocity P U Dynamic = 0.69 m/s and for the miniaturised design M U Dynamic = 2.13 m/s.The computed mean flapping velocity U Dynamic was compared with the quasisteady mean flapping velocity U to identify the effect of nonlinear deformation.In a previous study (Ishihara et al., 2014), U was based on a sinusoidal flapping motion and expressed as where r 2 = 0.6 is the dimensionless radius of the second moment of the wing area (Ishihara et al., 2014), L w is the spanwise length of the micro wing, and Φ and f are computed values from the nonlinear dynamic analysis (for the first prototype, P Φ L = 43.80°,P Φ w = 35.10°,and f = 50 Hz; for the miniaturised design, M Φ L = 87.02°,M Φ w = 39.25°, and f = 100 Hz).Based on the data and Equation ( 22), the calculated quasi-steady mean flapping velocity U and computed flapping velocity U Dynamic are compared in Table 3.There is a difference between U and U Dynamic because U is based on the sinusoidal flapping motion.The distinctions between U and U Dynamics give rise to differences in the theoretical lift forces derived from nonlinear dynamic analysis with the full FWNAV model and the numerical lift force obtained through FSI analysis of the flapping wing model, applied to both FWNAV models, which can be seen in Subsection 3.3.
Regarding the computational performance dependency on mesh size for nonlinear dynamic analysis, we have previously conducted an exhaustive investigation and documented our findings (Rashmikant, 2022).Our prior research firmly established that numerical results demonstrate negligible sensitivity to alterations in mesh size.It is important to note, however, that the problem and analysis setup employed in this manuscript differ from that study.

FSI analysis of micro wings
Figure 12(A) shows the design of the micro wing, which are the flexible parts of the FWNAV that produce aerodynamic forces.The problem setup for the FSI analysis of the micro wing is shown in Figure 12(B).The material properties of the polyimides for the FSI problem setup are given in Table 1.For the finite element modelling of the micro wing, 120 MITC shell elements with 143 nodes are used.
In the FSI problem setup (see Figure 12(B)), the stroke axis coincides with the wing attachment tip of the transmission (see Figure 8).At the stroke axis, a sinusoidal angular displacement Φ = Φ 0 sin2π f t is applied to flap the wing, where flapping frequency and angular amplitude Φ 0 ( = Φ L or Φ w ) are taken from the nonlinear dynamic analysis of the FWNAVs (see Section 2.2).The time increment t for the FSI analysis is T/5000, where T is the flapping period (inverse of f ).
For the FSI analysis, the structure, material properties, and finite element modelling of the fluid were as follows.The geometry of the rectangular fluid tank is shown in Figure 13(A), where the no-slip condition is imposed on all wall boundaries.The mass density and viscosity of the fluid are ρ a = 1.205 kg/m 3 and μ a = 1.837 × 10 −5 kg/m•s, respectively.Figure 13(B) shows the fluid mesh with stabilised linear equal-orderinterpolation velocity-pressure elements.The number of nodes is 50,460 and the number of elements is 271,922.Figure 13(C) shows the fluid and micro wing meshes in the xy plane.
The results of the FSI analysis are divided into two characteristic values, namely the numerical lift force and the feathering angle (wing motion).This section discusses the numerical lift force and the next section discusses the feathering motion.The time histories of the numerical lift forces for the first prototype and the miniaturised design are shown in Figure 14(A and B), respectively.The numerical lift force is used to evaluate the mean lift coefficient as shown in the flow diagram of flight performance (see Figure 4), discussed in Section 2.2.The mean lift coefficient can be evaluated using the following quasi-steady formula (Liu & Sun, 2008): where F is the numerical mean lift force, S is the area of the micro wing (S = L w c w ), ρ f is the fluid mass density, C is the mean lift coefficient, and U is the quasi-steady mean flapping velocity calculated using Equation ( 22) (see previous section).
The numerical mean lift force F for the first prototype and the miniaturised design from Figure 14 (A  and B) are respectively P F L = 0.0041 mN, P F w = 0.0025 mN, and M F w = 0.04 mN.From Equation ( 23) and the quasi-steady mean flapping velocity in Table 3, the calculated mean lift coefficient for the first prototype and the miniaturised design are P C w = 0.72, P C L = 0.79, and M C w = 2.40.
For the miniaturised design, only M C w is discussed because it was observed that the stiffness of the micro wing membrane (t w = 6.1 μm) (Rashmikant & Ishihara, 2023) is low for a high stroke angle; for example, M Φ L produces a mean lift coefficient M C L < 0.50 and a large feathering angle of more than 65°.According to a previous study (Dickinson et al., 1999), the lift force is maximum when the feathering angle is about 45°.Hence, in Section 4, the wing membrane thickness is designed to increase the stiffness of the wing membrane.
In the present quasi-steady model, the trapezoidal wing flapping velocity profile, as illustrated in Figure 15, leads to a uniform and steady fluid flow.This characteristic results in the insignificance of rotational lift and added mass forces in our quasi-steady aerodynamic force calculations, justifying their omission.In our ongoing efforts to refine the quasi-steady model and further enhance  aerodynamic lift, we plan to incorporate these factors in alignment with the work of other researchers (Chen & Zhou, 2021; Y. J. Lee et al., 2016;Sane & Dickinson, 2002).
In this context, a re-evaluation of computational performance concerning the FSI analysis mesh size has not been replicated, given its thorough investigation in prior research (Ishihara et al., 2014).This earlier work consistently demonstrates that numerical results remain unaffected by variations in mesh size.It is essential to note, though, that the problem and analysis setup in  that study differed from those employed in the present manuscript.
The mean lift coefficients are used below to evaluate the theoretical lift force from the nonlinear dynamic analysis of FWNAVs, as shown in the flow chart (see Figure 4).

Flight performance comparison between structural dynamic and fluid-structure interaction analyses
The feathering angle (wing motion) and the lift force obtained from the nonlinear structural dynamic analysis (see Section 3.1) and the FSI analysis (see Section 3.2) are compared in this section.
To validate the FSI analysis of a micro wing for the accurate prediction of flight performance, the feathering angle is first compared.The time histories of the feathering angle obtained from the nonlinear structural dynamic analysis of the two FWNAVs (first prototype and miniaturised design) are shown in Figure 16 and those of the feathering angle obtained from the FSI analysis of the micro wing are shown in Figure 17.The mean feathering angles obtained from these time histories are shown in Table 4 for comparison.As shown, θ L or θ w is similar to θ Dynamic .It follows from this result that the wing motion or flight simulation in the FSI analysis of the micro wing is similar to that in the structural dynamic analysis of the FWNAVs.This means that the aerodynamic effect on the micro piezoelectric drive (micro transmission and piezoelectric actuator) is negligible compared to the micro wing.Hence, it can be concluded that the FSI analysis of the micro wing is sufficiently accurate compared with the FSI analysis of the FWNAV if the latter would have been considered for wing motion and prediction of flight performance of FWNAV; the former reduces the computational cost and complexity associated with the latter.
The flight performance (i.e.lift force) determined from the two analyses is compared.The numerical mean lift force obtained from the FSI analysis of the micro wing for two FWNAVs was discussed in the previous section.The theoretical lift force from the nonlinear dynamic analysis of the two FWNAVs was calculated using Equation ( 23), values in Table 3 (only U Dynamic ), and the calculated mean lift coefficient ( P C w , P C L , and M C w ) in Section 3.2.The time histories of the theoretical lift force from the nonlinear dynamic analysis of FWNAVs are shown in Figure 18.The mean lift forces from the two analyses are compared in Table 5.As shown, the theoretical mean lift force from the nonlinear dynamic analysis of the FWNAV is larger than the numerical mean lift force from the FSI analysis of the micro wing.Based on this observation and flight mechanics (the fluid effect should not be ignored), the lift force is sensitive to the FSI, which should thus not be ignored in the estimation of flight performance.
In summary, the comparison of wing motion indicated that the aerodynamic effect on the micro piezoelectric drive is negligible compared to the micro wing and the comparison of the lift force indicated that the FSI should not be ignored in the flight performance evaluation.The FSI analysis of the micro wing was found to be sufficient for the accurate estimation of the wing motion and flight performance of an FWNAV.The mean lift force from the FSI analysis is less than the weight of the miniaturised design (10 mm) of the FWNAV, as shown in Table 6 (Rashmikant, 2022).To design and develop a flyable FWNAV, the flight performance obtained from the FSI analysis of the micro wing was optimised, as described in the next section.

Optimisation of flight performances using FSI analysis of micro wing
As discussed in Section 3.2, the FSI analysis of the micro wing for stroke angle Φ = 87°and flapping frequency f = 100 Hz produces a feathering angle larger than 65°a nd low lift (value normalised using the dynamic pressure is smaller than 0.5) for a wing membrane thickness of t w = 6.1 μm, which is a design solution reported in our previous study (Rashmikant & Ishihara, 2023).This means that the stiffness of the presently designed wing membrane is insufficient.According to our previous studies (Ishihara, 2018(Ishihara, , 2022;;Ishihara et al., 2014), the passive feathering angle depends on the membrane stiffness; this dependency leads to a lift force change.Furthermore, from Equations ( 22) and ( 23), the mean lift force is a function of stroke angle and flapping frequency, which are important flight maneuvre parameters (Onishi & Ishihara, 2019).
According to a previous study (Dickinson et al., 1999), the lift force is maximum when the feathering angle reaches its optimal value, which is about 45°.In this section, we design the wing membrane thickness such that it produces a feathering angle of about 45°for an individual set of stroke angles and flapping frequency and then maximise the lift force to design a flyable FWNAV (see Table 6) (Rashmikant, 2022).

Basic problem setup for optimization of flight performance
The flight performance is maximised by performing a parametric study using the FSI analysis of the micro wing.The problem setup for the parametric study is similar to that for the FSI analysis of the micro wing, as shown in Figure 12(B).The difference is that in this parametric study, there are three variable parameters, namely the design parameter (i.e.wing membrane thickness t w ), sinusoidal input parameter (i.e.stroke angle Φ), and flapping frequency f.
The wing membrane thickness t w varies from 10 to 30 μm.The lower bound of t w was selected such that the feathering angle is less than 60°in the flapping frequency range of 100-125 Hz, as designed in our previous study (Rashmikant & Ishihara, 2023).The lift coefficient is thus in the range of 1-2 (Dickinson et al., 1999).The upper bound of t w was determined based on the actual measurement for a polyimide sheet after a curing process during polymer micromachining (Rashmikant et al., 2021a).A stroke angle of Φ 0 = 80°was selected since our previous studies showed that our FWNAV flapped in the frequency range of 100-125 Hz and produced a stroke angle in the range of 80-125° (Rashmikant, 2022;Rashmikant & Ishihara, 2023).A flapping frequency of f = 100-250 Hz and Φ 0 = 80°were chosen such that the Reynolds number Re for these sets of flapping frequency and stroke angle is comparable with that for the actual flapping frequency and produced stroke angle in the FWNAV, which is validated later.Here, Re is given as (Ishihara et al., 2014) where ρ a = 1.205 kg/m 3 is the mass density of air, μ a = 1.837 × 10 −5 kg/m•s is the dynamic viscosity of air, c m = 3.2 mm is the chordwise length of the micro wing, and U is the mean flapping velocity, which can be evaluated using Equation ( 22).From Equation ( 22

Design of wing membrane thickness for maximizing lift force
The wing motion (i.e. the amplitude of the feathering angle at the middle of each half stroke) and the mean lift force were taken from the FSI parametric study.The wing motion is first discussed.The amplitude of the feathering angle versus wing membrane thickness t w is shown in Figure 19(A) for the sets of frequency f = 100-250 Hz and stroke angle Φ 0 = 80°.As shown, satisfactory design solutions or the design window (Ishihara et al., 2002;Rashmikant & Ishihara, 2023) can be obtained for t w such that the feathering angle is about 45°( optimal value) and the lift force is maximum.Thus, the design window for t w is 12-22 μm in the frequency range of f = 100-250 Hz.This implies that within the range of membrane thickness considered, and under the conditions of stroke angle (Φ 0 ) at 80°and flapping frequency (f ) ranging from 100 to 250 Hz, the lift force exhibits an increasing trend.This augmentation in lift force is pivotal, as it contributes to the maximisation of the total lift capable of supporting the weight of the FWNAV.
The mean lift force was maximised to obtain a flyable FWNAV.The mean lift force versus t w is shown in Figure 19(B).The mean lift force is more than the weight of the FWNAV (see Table 6) for f = 200 Hz, Φ = 80°, and t w = 12-22 μm.Although f = 200 Hz and Φ = 80°don't fall in the range of designed frequency of the FWNAVs but when we consider the Re similarity as Re(80, 200) ( = 795) ≈ Re(125, 125) ( = 780), we can say that our FWNAV can fly at f = 125 Hz and Φ = 125°, which is in the range of the designed values in previous studies (Rashmikant, 2022;Rashmikant & Ishihara, 2023).Thus, it can be concluded that for tethered flight, t w should be 12-22 μm for an FWNAV with a size of 10 mm or smaller.Furthermore, the Reynolds number in our study closely aligns with that of our model insect, Eristalis Tenax, providing additional substantiation of the flight capabilities inherent in our designed FWNAV.
Upon a careful examination of Figure 19(A and B) in conjunction, it becomes evident that when the wing stiffness is low, and there is significant angular motion in  flapping (comprising both stroke angle and flapping frequency), wing deformation becomes pronounced.Consequently, this deformation significantly influences the feathering angle and lift force.This observation effectively validates the fundamental principles of aerodynamics underlying our study.
To facilitate a comprehensive discussion of FSI within the domain of fluid dynamics, let's delineate the specific parameters governing our FSI analysis.This configuration involves a 15 μm-thick membrane flapping consistently at a frequency of 200 Hz, with a fixed stroke angle of 80°.In Figure 20, we present a comprehensive time cycle depiction that includes lift, feather angle, and stroke angle within a single frame.This analysis yields several noteworthy conclusions: (1) Reduction in the feathering angle correlates with a decrease in the lift component of the aerodynamic force perpendicular to the wing surface.
(2) The feathering angle achieves its maximum immediately after the stroke motion changes direction, followed by a rapid decrease.(3) Just before the midpoint of the half-stroke, the feathering angle reaches its maximum.(4) Feathering primarily results from the inertial force generated during the stroke transition.
Figure 21 provides the wing's motion corresponding to Figure 20, offering a visual representation of the wing's deformation throughout one complete cycle, spanning from the 5th to the 6th cycle.Reiterating the FSI setup: a 15 μm-thick membrane, a flapping frequency of 200 Hz, and a fixed stroke angle of 80°. Figure 21 visually represents both the wing's motion and the accompanying flow field.In the left column, during the upstroke, the wing undergoes a right-to-left flapping motion, while in the right column, during the downstroke, it flaps from left to right.These depictions emphasise the passive pitching motion, mirroring the behaviour observed in natural insects, effectively captured by our FSI model.Specifically, the wing maintains a high angle of attack at the midpoint of each half-stroke and experiences rapid rotations during stroke reversals.The significance of this portrayal lies in the discernible formation of the leading-edge vortex, a pivotal flow structure in the context of insect flapping wings.Notably, the leading-edge vortex separates during the mid-stroke phase of the stroke motion.This phenomenon has been extensively documented in prior research, as cited in our discussion (Bhat et al., 2019;Chen et al., 2020;Eldredge & Jones, 2019;Ellington et al., 1996).Importantly, our figures illustrate the stable attachment of these vortices to the wing, a critical factor contributing to the substantial translational lift generated at the midpoint of each half-stroke.

Conclusion
The core aim of this study is to unveil an inventive hierarchical coupling methodology, meticulously crafted for conducting a comprehensive evaluation of the flight capabilities inherent in polymer micromachined FWNAVs, all through the utilisation of cuttingedge computational tools.Simultaneously, we endeavour to tangibly exhibit the practical feasibility of both tethered and autonomously flying FWNAVs.Notably, these twin objectives encapsulate the distinctive novelty of this present study.
To effectively evaluate the flight performance of two distinct FWNAV prototypes while managing the intricate multiphysics interactions, we proposed a hierarchical coupling approach.In this novel method, we divide the hovering FWNAVs into two distinct subsystems: oneway coupling, responsible for modelling the micro piezoelectric drive system, and strong coupling, addressing the interactions between the micro wing and the surrounding air.
In the one-way coupling analysis, we conduct nonlinear structural dynamic assessments of the micro piezoelectric drive system.The resulting responses are then transmitted to the micro wing through Dirichlet boundary conditions.In contrast, the strong coupling analysis involves solving for the micro wing and the surrounding air in conjunction to fully account for their intricate interactions.
The connection between one-way and two-way coupling is established by using the output (sinusoidal angular displacement) from the former as the input for the latter.Specifically, the feathering angle (wing motion) and lift force obtained from both one-way and twoway coupling analyses were compared.This was done to validate the accuracy and efficiency of employing FSI analysis specifically for micro-wing design, as opposed to a more computationally intensive FSI analysis for overall FWNAV design, which would entail higher computational costs and complexity.
Employing the proposed methodology, we conducted an extensive evaluation of the flight performance for two distinct FWNAV designs.Our findings indicate that the first prototype exhibited lower lift force and feathering angles compared to the miniaturised design.This difference was attributed to the smaller sinusoidal angular displacement observed in the former FWNAV.
Furthermore, utilising the same methodology, we optimised the wing membrane's stiffness by varying the membrane thickness within the range of 12-22 μm.This optimisation aimed to achieve the optimal feathering angle or maximum lift force.The results of this optimisation effort demonstrated a mean lift force exceeding the weight of the designed FWNAV, which is a significant milestone in enhancing flight performance.
Consequently, this study culminated in the design of a functional FWNAV with a compact size of 10 mm, which can be feasibly fabricated using polymer micromachining techniques.
In our current investigation, we have achieved lift forces slightly exceeding the weight of the FWNAV.However, it is important to note that for stable hovering, the generated lift force ideally should be in the range of 2-2.5 times the weight of the FWNAV.Addressing this limitation requires further refinements to our quasi-steady model and FSI programme.These enhancements should encompass the incorporation of rotational lift and added mass forces into our analysis, paving the way for more precise and effective lift generation.

Figure 2 .
Figure 2. Coupled problem for FWNAV.Region enclosed by dashed lines is the partition method used for the flight performance evaluation of FWNAVs.

Figure 4 .
Figure 4. Flow diagram of (A) flight performance evaluation and (B) flight performance optimisation of FWNAVs.

Figure 6 .
Figure 6.Detailed design of polymer micromachined FWNAV.(A) Plan view and (B) sectional view.

Figure 7 .
Figure 7. Plan view of material distribution in FWNAV.(A) Front and (B) rear.

Figure 8 .
Figure 8. Problem setup for nonlinear structural dynamic analysis of FWNAV.
Figure 10 shows the time history of the flapping angular displacement for the first prototype and the miniaturised design.As shown, the evaluated stroke angles for the first prototype are P Φ L = 43.80°andP Φ w = 35.10°andthose for the miniaturised design are M Φ L = 87.02°andM Φ w = 39.25°.Note that the superscripts P and M denote the characteristic values for the first prototype and the miniaturised design and the subscripts L and W denote the characteristic values due to the leading edge and wing attachment deformation, respectively.

Figure 9 .
Figure 9. Illustration of Φ L and Φ w measurements at the moment of maximum Φ L , clarifying the rationale behind the greater value of Φ L compared to Φ w .Here, Φ L represents the leading-edge stroke angle, and Φ w signifies the wing attachment stroke angle.

Figure 11 .
Figure 11.Flapping velocity time history for (A) first prototype and (B) miniaturised design.

Figure 12 .
Figure 12. (A) Plan view of micro wing design, (B) Problem setup and mesh setup for FSI analysis of micro wing.

Figure 14 .
Figure 14.Time history of the numerical lift force obtained from the FSI analysis of the micro wing, illustrating results for (A) the initial prototype and (B) the miniaturised design.It is important to highlight that the mean lift force was determined by considering the data encompassed within the dashed-line rectangle over a single cycle.

Figure 15 .
Figure 15.Time history of flapping velocity from FSI analysis of micro wing for (A) first prototype and (B) miniaturised design.Flapping velocity shows the trapezoidal behaviour.

Figure 16 .
Figure 16.Time history of feathering angle from nonlinear structural dynamic analysis of (A) first prototype and (B) miniaturised design.

Figure 17 .
Figure 17.Time history of feathering angle from FSI analysis of micro wing for (A) first prototype and (B) miniaturised design.Note that the data enclosed by the dashed-line rectangle for an individual cycle give the mean feathering angle.

Figure 18 .
Figure 18.Time history of theoretical lift force from nonlinear structural dynamic analysis of (A) first prototype and (B) miniaturised design.
), U is a function of stroke angle Φ and flapping frequency f.Hence, Re is a function of Φ and f and can be written as Re(Φ, f ).The specific values of Re for certain sets of (Φ, f ) are as follows: Re(80, 250) = 995, Re(80, 225) = 895, Re(80, 200) = 795, and Re(125, 125) = 780.These values of Re validate that Re for sets of f = 100-250 Hz and Φ 0 = 80°chosen for the present FSI parametric analysis is comparable with that for the actual flapping frequency range of f = 100-125 Hz and the produced stroke angle Φ 0 = 80°-125°in the FWNAV.

Figure 19 .
Figure 19.Optimisation of flight performance; (A) Feathering angle and (B) mean lift force.

Figure 20 .
Figure 20.A comprehensive representation of the time cycle, incorporating lift, feather angle, and stroke angle within a single frame, to facilitate an insightful discussion of FSI.Detailed information regarding the problem and analysis setup for this result can be found in the main body corresponding to this figure.

Figure 21 .
Figure 21.Illustration of the wing's flapping motion and the surrounding fluid flow during one stroke.In the left column, the wing flaps from right to left, while in the right column, the wing flaps from left to right.(A) Time instant of 5.25 cycle; (B) The time instant of 5.32 cycle; (C) The time instant of 5.38 cycle; (D) The time instant of 5.44 cycle; (E) The time instant of 5.50 cycle; (F) The time instant of 4.75 cycle; (G) The time instant of 4.82 cycle; (H) The time instant of 4.88 cycle; (I) The time instant of 4.94 cycle; (J) The time instant of 5.00 cycle.

Table 2 .
Analysis setup for nonlinear structural dynamic analysis of FWNAVs.

Table 3 .
Comparison of mean flapping velocity obtained with nonlinear structural dynamic analysis of FWNAVs and quasisteady formula.

Table 4 .
Comparison of wing motion (i.e.mean feathering angle) obtained with nonlinear dynamic analysis of FWNAVs and FSI analysis of micro wing.

Table 5 .
Comparison of flight performance (i.e.mean lift force) obtained with nonlinear dynamic analysis of FWNAVs and FSI analysis of micro wing.