Fluid–structure interaction simulation on flight performance of a dragonfly wing under different pterostigma weights


 The dragonfly wings provide insights for designing an efficient biomimetic micro air vehicle (BMAV). In this regard, this study focuses on investigating the effect of the pterostigma weight loading and its spatial location on the forewings of dragonfly by using the fluid–structure interaction simulation. This study also investigates the effect of change in the wing elasticity and density on the wing performance. The forewing, which mimics the real dragonfly wing, is flat with a 47.5 mm span and a 0.4 mm thickness. The wing was set to cruise at 3 m/s with a constant flapping motion at a frequency of 25 Hz. This study shows that a small increase of pterostigma loading (11% of wing weight) at the tip of the wing significantly improves the lift to drag ratio, CL/CD, which has 129.16% increment in comparison with no loading. The lift to drag ratio depends on the pterostigma location, pterostigma loading, elastic modulus and density. The results of this study can be used as a reference in future BMAV wing optimization design.


INTRODUCTION
The dragonfly is a fascinating creature that has four wings and is known to be an efficient predator due to its advanced flight abilities. The unique dragonfly wings were studied by previous researchers such as Jongerius and Lentink [1], Azuma et al. [2] and Bomphrey et al. [3] in an effort to understand the flight mechanism that contributes to its high maneuverability. This knowledge may help biomimetic engineers to develop a novel biomimetic micro air vehicle (BMAV) wing design that can be more cost-effective than the conventional MAV design. Thus, the study on the dragonfly is worthy of further investments.
The understanding of the structural mechanics of a dragonfly wing is essential to facilitate the design of biomimetic flyer wings. There were a few notable research works such as those conducted by Chen et al. [4], Sun and Bhushan [5] and Zhang et al. [6], which focused on studying the dragonfly wing structure and mechanical properties. Chen et al. [4] conducted a laser vibrometer test on the dragonfly wings' leading-edge veins. Their findings show that the dead dragonfly wings have elastic modulus identical to that of low-density polyethylene, while the live dragonfly wings have elastic modulus similar to that of elastic hydrocarbon polymer. Sun and Bhushan [5] discussed the flight ability of dragonfly by focusing on the aspects of morphology, and structural and mechanical properties of the wings. Their findings show that the tapered vein and membrane, and the corrugated shape of the dragonfly wings affected the flight performances by reducing structure weight while having high strength and stiffness. Zhang et al. [6] conducted a tensile test on the dragonfly wings and found that veins and membrane structure have a significant contribution to prevent crack during flight. Other researchers proposed alternative structural design and material to construct biomimetic flyer wings, attempting to reproduce the wings of the dragonfly for practical applications. Sivasankaran and Ward [7] used the spatial network analysis to design a biomimetic dragonfly-inspired wing that has simplified vein patterns and is fit for laser cutting machine fabrication methods. Sivasankaran et al. [8] tested various materials such as stainless steel, balsa wood, black graphite carbon fiber and red prepreg fiberglass to fabricate dragonfly-inspired BMAV wings. Their findings show that red prepreg fiberglass and black graphite carbon fiber are more appropriate materials for the fabrication of mentioned BMAV wings due to their better performance structurally. The elastic properties of dragonfly-inspired flapping wings made of acrylonitrile butadiene styrene (ABS), polylactic acid and acrylic, respectively, were also studied by Sivasankaran et al. [9] by using a high frame rate imaging system. They observed that the ABS BMAV wing possesses significant flexibility in the chordwise direction, which may improve the flight performances. The above-mentioned findings can be beneficial to optimize the structural design of biomimetic flyer wings.
The aerodynamics and kinematics of a dragonfly wing were investigated to deduce the causes of its high flight performance and maneuverability. Broering and Lian [10] investigated the aerodynamic effects of flapping wings operating in a tandem configuration by using a numerical method. They found that the lift and thrust increased due to the positive vortex interaction between the tandem wings at zero phase lag. Chen et al. [11] used a two-dimensional unsteady aerodynamic model derived from potential flow theory for predicting wing motion that mimics the dragonfly wings. They concluded that the lift coefficient could be increased by taking advantage of the aeroelastic deformation of the wings at upstroke and downstroke positions. Chen et al. [12] utilized high-speed videography on the dragonfly to study the kinematics of figure-eight flapping motion. They found that the wing rotation was improved by the existing nodus, which hinged and halved the leading-edge spar, allowing a 40°angle of movement between the two pieces of wing sections during flapping motion. Gadde and Vengadesan [13] studied the Lagrangian coherent structures of wings flapping in a tandem configuration. They discovered that the in-phase stroke of both tandem wings generates high vertical forces and thrusts, which improve take-off while counter-stroking reduces fluttering during hovering. Xie and Huang [14] investigated the vortex interactions between the forewing and hindwing of a dragonfly, and their findings show that the lift forces vary at different interaction modes. Some researchers studied the aerodynamic effect due to the surface features of dragonfly wing; for example, Kim et al. [15] showed that the dragonfly wing corrugation increases lift coefficients and has a minimal effect on drag coefficients. Chen and Skote [16] conducted numerical simulation on the gliding performance of corrugated and profiled dragonfly wings. Their simulation results show that the corrugated wing outperforms the profiled wing by 6% higher in the lift to drag ratio. These findings can be the references for designing the optimal geometry of biomimetic flyer wings. Some researchers designed and built prototypes of microflyer and tested their flight performances. Li et al. [17] designed a flapping wing rotor that combines the insect flapping with rotating wing motions. The result of this combination is a flyer that has increased lift with efficiency between insect flapping and rotating wing motions. A prototype of the mentioned flyer was tested experimentally by Guo et al. [18].
Despite that there are plenty of previous studies on the material, structural mechanics, aerodynamics and kinematics of the dragonfly wings, there were limited studies dedicated to the pterostigma. The pterostigma is a conspicuous feature of a dragonfly wing [19]. It is defined as a concentrated nontransparent spot located near the tip and leading edge of the dragonfly wings [20]. Ennos [21] observed that the wing mass of dragonflies usually tapered, where the wing base has higher mass while the wing tip has lower mass. In addition, the dragonflies possess pterostigma, a spot containing larger mass than the rest of the wing sections, which has 5-10% of total wing mass for dragonflies, situated near the leading edge. There were several suggestions proposed by researchers on the functions of pterostigma based on observations. Needham [22] proposed that the pterostigma joins the veins of the dragonfly wing at the leading edges and tips. The pterostigma was also suggested to increase the wing tip inertia, which will eventually increase the efficiency of flapping wings [22]. Tillyard [23] exemplified the different types of pterostigma in different species of dragonflies and proposed that they have a strong relationship with flight performance. There were some experimental studies involving the investigation of the effect of pterostigma load on the dragonfly flight performances. Norberg [24] conducted wing loading experiments that used small patches of adhesive tapes to increase the pterostigma load on a dragonfly wing. The experimental results showed that different pterostigma loads were able to assist in balancing the chordwise mass distribution and may prevent fluttering of wings by regulating the pitch, especially in a gliding motion. The pterostigma was also known to increase the wing-tip amplitude of the flapping motion [25]. The increase of wing-tip amplitude enables wider flapping motion that will interact with more airflow to produce higher lift. The above-mentioned works were two related studies that investigated the effect of pterostigma load on dragonfly wings by experimental methods. Despite the above-mentioned observations by previous researchers, the function of the pterostigma in the dragonfly wing flight performances was not yet fully understood, and there were limited publications that focus on this aspect.
With the advancement of computational technologies, computational fluid dynamics (CFD) have been widely used for predicting the fluid flow behavior in many research works [26][27][28]. Coupled with the computational structural dynamics (CSD), the fluid-structure interaction (FSI) simulation became popular in predicting the insect or BMAV flapping wing motions. Takizawa et al. [29] suggested a deforming spatial domain formulation technique for studying a locust wing deformation and flapping aerodynamics. The proposed computational techniques successfully modeled the flapping motion of locust in flight with deformation taken into account. Chen et al. [30] simulated a dragonfly-inspired MAV by using commercial FSI solvers to study its flapping motions. The mentioned examples proved that FSI is a reliable computational method for this study.
Due to the complex nature of flapping wings, there were limited computational studies that utilize FSI methods for investigating the wings in addition to pterostigma load. Previous researchers such as Norberg [24] performed wing loading experiments on dragonflies and found that certain pterostigma loading will produce an optimum flight performance. However, there was no detailed investigation involving the effect of different pterostigma locations, pterostigma loading, wing elastic modulus and wing density on flight performance. Thus, the objective of this study is to predict the flight performance involving the lift and drag coefficients of the dragonfly wing model by considering the manipulating variables of pterostigma location, pterostigma loading, wing elastic modulus and wing density.

METHODOLOGY 2.1 Geometry model
For this study, the primary focus is on the forewings of the dragonfly and the dynamic behavior of flapping wings under a uniform thrust flight condition. The leading and trailing edges of the wing were truthfully presented, as shown in Fig. 1, based on the dragonfly planform shape given in [24]. The cross-section A-A denotes the section mesh of Fig. 3. The wing material properties for the wing were interpreted from [1]. See Table 1 for the dimensions and material properties of the model and dragonfly wings. For this study, the hindwing of the dragonfly was not considered to avoid any downstream complex flow  affecting the forewing behavior. Tandem wing simulation is not necessary as the main objective is to test the effect of pterostigma on a wing. The forewing is chosen because the forewing has no wake interference due to its location at the upstream. Thus, the basis for this study is valid without the consideration of tandem wing simulation. The thickness of the wing is assumed to be uniform for the simplicity of FSI investigations. The basis for using a constant thickness is to ease the fabrication process of the MAV. Young's modulus is not a fixed parameter in this study as it is manipulated to observe its effect on the performances of the wings. Thus, the material properties of a real dragonfly are unnecessary to be adopted in this study.

Governing equations
The theoretical equations of the CFD and CSD software were referred to ANSYS Fluent [31] and ANSYS Mechanical APDL theory guides [32]. The mass conservation equation used in the CFD simulation is as follows: where ∇ is the divergence, ρ is the density, j is the mass flux, v is the flow velocity and t is the time.
The momentum conservation equations were adopted as well to simulate the field pattern of incompressible flow, as shown below: where v is the flow velocity, p is the pressure, τ is the stress tensor and F s represents the external forces. The transport equations of the shear stress transport (SST) kω turbulence model include the kinetic energy and the specific dissipation rate, as shown below: where E k and E ω are turbulence kinetic energy and specific rate generation by velocity gradient, I k and I ω are the k and ω diffusivity, and F k and F ω are external forces. The SST k-ω turbulence model was adopted in this study as it is known to be accurate in predicting turbulent flow, converges faster and consumes minimal computational time [33]. The SST k-ω turbulence model is also able to predict the viscous boundary layer when the meshes near the wing surfaces are having dimensionless wall distance y + < 1. The research work of Radermacher [34] is an example where the SST k-ω turbulence model was adopted for flappingwing simulation. For structural dynamics, the stress-strain relationship is defined as where σ n is the stress vector (n = x, y, z, xy, yz, xz), K is the elasticity and ε is the elastic strain vector. Equations (1)-(5) show the fluid dynamics equations, which are based on the Navier-Stokes fluid flow equations. These fluid flow equations were solved in the ANSYS Fluent, which is the main fluid dynamics code of the ANSYS Workbench. Equation (6) shows the stiffness matrix equations for the structural part of the FSI. This structural dynamics equation was solved in the AN-SYS Mechanical ADPL, which is the main structural dynamics code of the ANSYS Workbench. A system coupling code available in the ANSYS package was used to act as intermediary software that exchanges the data between the fluid and structure solvers. The exchanged data include the forces and torques from the fluid solver and the deformation of the wing from the structure solver. The system coupling code updates the mentioned exchanged data as new inputs in both the fluid and structure solvers at alternating turns for each time step, thus allowing the modeling of FSI for the dragonfly wing simulation. The mentioned system coupling method is known as the partitioned type of FSI coupling, as shown in [35].
For the airfoil analysis of the wing model, the lift and drag coefficients were defined as follows [36]: where C L is the lift coefficient, C D is the drag coefficient, F L is the lift force, F D is the drag force, ρ is the density of air, V is the velocity of airflow and A r is the reference area.
A dimensionless constant consisting of the elastic modulus and density as manipulating variables was introduced as shown in the following equation by referring to the natural frequency equation of a vibrating cantilever beam from [37]: where D c is the dimensionless constant, f is the frequency, E is the elastic modulus and A is the cross-sectional area of the wing base.

Boundary conditions and properties
The outer domain of the simulation case consists of velocity inlet, pressure outlet and symmetry surfaces at the hinge side of the wing model, with the rest considered as walls. The wing model surfaces were set to wall boundary condition. A steady wind was applied in the positive Y-axis direction, where the leading edges of the wing model were at the upstream. The relative wind was modeled as turbulent flow with a Reynolds number of ∼2000 in the vicinity of the wing model. Despite the Reynolds number of 2000, the turbulence model is still necessary for the flow computation due to the coarse mesh adopted in this study. This is so because the current computational resource is not capable of performing simulation work of the direct numerical simulation, which does not require turbulence models. The definition of Reynolds number depends on free stream velocity. The wing was subjected to a simple harmonic flapping motion with the input, as shown in the following equation: where A t is the angle of the flapping wing tip, A m is the amplitude, f is the frequency and t is the time. The maximum flapping angle is set to 15°for all simulation cases. The maximum flapping angle of 15°is adopted due to the limitation of the simulation model where the mesh deforms too much, and negative cell volume occurs as a result of that. However, the flapping angle of 15°i s adequate for a comparative study as the results show a significant effect of changing pterostigma loadings despite the applied maximum flapping angle of 15°. The location of the hinge was at   Table 2 shows the numerical simulation settings and initial conditions for the FSI simulation of the wing model. The flight speed, flapping frequency and pterostigma loading of a typical dragonfly can be referred to in [38,39]. The real dragonfly wing geometry is complex and made up of different materials. The geometry itself is too complex to be modeled, given the limitation of the available computational resources. Thus, a simplified geometry with constant thickness is modeled instead. Since the fabrication of a constant thickness wing is feasible for any MAV, this study contributes to investigating the effect of pterostigma on the mimicked wings of MAVs. Thus, real dragonfly conditions are not necessary to be modeled here in this study. The FSI cases have to be made simple due to the limitation of the available computational resources. The FSI simulation solves the fluid flow and structural dynamics equations. It is more computationally intensive to solve multiphysics problems than single-field problems. The interactions between the fluid flow and structural dynamics solvers also require additional computational power. The FSI model involves transient mesh deforming at every time step, which requires higher computational power than that of the CFD alone. As the primary focus of this study is to get a better idea of the effect of the pterostigma point load itself, we believed that the best reference condition would be at zero pitch angle. This zero pitch angle simplifies the mesh deformation pattern, which will reduce the remeshing work.
Another reason why zero pitch angle was our focus in this study can be related to our present and ongoing effort on developing an MAV based on a flapping-type drone (see Fig. 2).
The MAV has flapping wings at almost zero pitch and the model can flap different speeds. The flapping action of the wings creates wobbling behavior of the wing sheets, especially at the back. Thrust generated is inclined downward (−45°to −60°) from the horizontal axis. This is a workable model. We are still investigating the possibility of improving its flight performance considering adding a localized point mass on the wings. We hope that our results in this manuscript can be useful for the people working with a flapping-type drone, wind turbine, etc.

Numerical method
The FSI two-way system coupling between the ANSYS Fluent (version) and ANSYS Mechanical APDL was adopted. The performance of the wing model was predicted by applying a constant frequency of the flapping motion with the weight of the pterostigma as the manipulating variable. The PISO (pressureimplicit with the splitting of operators) algorithm was used in this study, together with the SST k-ω turbulence model. The PISO algorithm is well suited for transient analysis as it is the scheme designed to solve pressure-velocity fluid flow equations. The FSI involves a two-way system coupling of transient structural and fluid solvers. The Courant number, derived from the Courant-Friedrichs-Lewy condition, was controlled to be <1 following the recommendation of the user guide by setting the time step based on the following equation: where Co represents the Courant number, V m is the velocity magnitude, t is the time step and x is the mesh node interval. It was tested that for the current mesh density, the Courant number of <1 could be achieved by selecting the adaptive time step function. The default convergence criteria of the solution were adopted, which scale the residuals to <10 −3 for every fluid flow equation. It was suggested that the mentioned target residual value is generally sufficient for this type of simulation case. Under the spatial discretization schemes, the least-squares cell-based method was applied in gradient. The second-order discretization scheme was used for pressure. The second-order upwind scheme was used for momentum. However, to avoid complications in the complex flow, which will diverge the numerical solution of system coupling, first-order upwind schemes were used to discretize the turbulent kinetic energy and specific dissipation rate of the fluid flow equations. The transient formulation was set to first-order implicit as this formulation is sufficient for this study, as recommended by the user guide.

Simulation cases
There are, in total, four FSI simulation cases conducted in this study, which include the cases for optimizing the location of pterostigma, the amount of pterostigma load, wing elastic modulus and wing density. The optimum pterostigma location is determined first, followed by the pterostigma load. For the initial step, an 8% pterostigma load was set to different locations according to the planform view, as shown in Fig. 1. Then, after the optimized location is deduced, the pterostigma load was varied in between the range of 8-14% of wing weight. A simulation case without the pterostigma load was performed as well to infer the flight performance improvement. A typical dragonfly pterostigma load was 9% of its total wing weight [9]. In order to understand the effect of using different materials on the proposed dragonfly-inspired wing model performance, variables of elastic modulus and density were applied to the simulation. The wing elastic modulus and density were then varied with the optimized pterostigma loading and location to see their effect on flight performance.

Meshing and grid independence test
The meshing was done by using the ANSYS Meshing module for both the ANSYS Mechanical, transient structural analysis, and ANSYS Fluent solvers. Tetrahedral mesh cells were adopted for the fluid mesh domain, while the prism cells were adopted for the solid mesh domain. Five inflation layers with a growth rate of 1.2 and 0.272 transition ratios were generated from the wing surfaces to make the dimensionless wall distance y + < 1, until it reaches the viscous sublayer of the boundary flow. In order to conserve limited computational power and optimize the mesh domains, a grid independence test was conducted. A single test case of a wing model without pterostigma load, with a flapping frequency of 25 Hz and a wind speed of 3 m/s, was simulated with coarse, medium and fine meshes. The lift and drag coefficients at 0.1 s were converged at the fine mesh, as shown in Table 3. Thus, the mesh properties of the fine mesh domains were adopted for this study. t = 0.1 s was chosen for the grid independence study because the simulation is converged at t = 0.1 s, and further time period will incur an unnecessary loss in computational resources. t = 0.1 s is the most economical time period for the grid independence study tested by trial and error while giving adequate time for simulation to converge. Thus, the grid independence test is sound in this study. Figure 3a and b shows the wing chord mesh section of the fluid mesh domain, while Fig. 3c shows the mesh of the structural domain.

Solver validation
The validation of the current simulation model was done by referring to an existing dragonfly wing section study conducted by Vargas et al. [40]. The validation was done by comparing the results of lift and drag coefficients of NACA 0008 with the results from [41] , similar to the validation work of Vargas et al. [40]. The wing has a constant chord length of 0.01 m with a span of 0.0475 m. Both the chord and span of the airfoil wing were set to have a similar size to the biomimetic wing model in order to validate the model and setup. It was presumed that the validated model could be used to predict the performances of the biomimetic wing model by taking into consideration the discrepancies of the results of the airfoil wing simulation cases with the reference. The airfoil wing was in prism shape, where there was no tapering at the base and tip due to the available reference being modeled as a 2D infinite airfoil. The airfoil wing was subject to the same boundary conditions and simulation settings. The validation shows that the current simulation model has comparable accuracy, where the highest percentage difference with the reference was only 17.5%. The percentage difference may be caused by the tip loss effect as well as the spanwise interaction of the 3D model, which are not taken into account in the 2D model results from [41]. The flow separation along  the wing is not intended to be modeled in this study in the first place. Thus, the coarse mesh is adopted with the implementation of wall functions. According to Fluent, wall functions for omega turbulence models will be applied if the mesh falls outside the viscous layer of the flow boundary layer, but since in this study the meshing indeed reached the viscous boundary layer, the flow separation is considered in the simulation. The 17.5% difference between drag coefficients is quite accurate as it is impossible to match the experiment with simulation results exactly given there are inherent errors in the solver and turbulence models and also uncertainties related to the coding of software, which is beyond the responsibility of authors. Despite the 17.5% error, the simulation results are sound because they are based on comparison with the presumption that the magnitude of error will not vary too much for each simulation case. A comparative study is conducted by changing the pterostigma location, load value and material properties. The discussions are based on the compared results, which will make the errors insignificant as all results will have almost the same magnitude of the error. Thus, the results are adequate for comparative study. Figure 4 shows the closeup mesh section, while Table 4 shows the CFD validation results of the mentioned NACA 0008 airfoil wing.

RESULTS AND DISCUSSION
It was presumed that the addition of pterostigma weight would increase the performance of the wing model by increasing the wing-tip amplitude. Thus, the lift and drag coefficients were the main concern in this study, as these variables will determine the flight performance of the proposed dragonfly-inspired wing model. The lift to drag ratio was also calculated to decide the optimized pterostigma location, loading, elastic modulus and density of the wing model. The lift to drag ratio is calculated by dividing the lift forces generated by the wings by the drag incurred due to the friction forces and also the pressure forces against the oncoming relative wind. The wings produce lift, as shown in the results. Thrust is not considered in this study, and the main concern is the ability for the dragonfly-inspired MAV to be able to Coefficients stay in position in forward flight. Thus, only the lift to drag ratio is taken into account as it is the most important parameter to determine whether an airfoil is efficient or not. Without the lift force, there will be no flight. The forces, stresses, lengths, times, etc. can be deduced from the graph as the normalized terms are obtained by dividing the maximum values, which can be found in the text.

General fluid flow pressure and velocity contour plots of flapping wing
The flow pressure contour was plotted at a 75% span of 0% and 11% pterostigma loads at upstroke and downstroke positions, as shown in Fig. 5. A comparison was made to understand the effect of the presence of pterostigma on the wing. The 0% and 11% pterostigma loads are the two values with distinctive features and roles in this study. The 0% pterostigma load was set as the control reference to reflect the effect in the presence of pterostigma, while the 11% pterostigma load is the load that produces the maximum lift on the wing according to the simulation results. The rest of the pterostigma loads have effects in between those for the 0% and 11% cases. Hence, they were not chosen for detailed discussion. The addition of the pterostigma weights changes the surface pressure distribution of the wing model, where the magnitude of negative and positive pressures at the leading edges of upstroke and downstroke flapping motion was reduced in the presence of pterostigma weights. This phenomenon is confirmed by the airflow pressure distribution contour at the vicinity of the chordwise section of the wing model, where the leading-edge pressure magnitudes are reduced by the presence of pterostigma weights. The pterostigma weights shifted the maximum positive pressure toward the leading edges, as shown in the chordwise flow pressure distribution contours. This pressure zone shift may be due to the slight twist of the  wing section, which aligned the leading edges to point toward the wind direction. The wake of the wing model at upstroke and downstroke is shown in Fig. 6, where the flow velocity distribution contour is plotted at the chordwise section plane. The high-velocity magnitude region at the bottom of the upstroke and at the top of the downstroke flapping motion indicates that the flapping motion altered the free stream velocity to generate lift force. The velocity contour plots confirm that the presence of the pterostigma load improves the performance of the wing. It was noted that the current contour plots are less refined due to the coarse mesh used in this study. Meshes as fine as those used in other studies were not adopted due to the limitation of the current workstation, where finer mesh density will give too much burden to the PC and consequently make the simulation unable to progress. Despite this limitation, the contour plots were still able to show the correlations between the airflow pressure and velocities.

Blade section pressure distribution of the wing model
Pressure distribution plots on the blade sections at the same location as the pressure and velocity contour plots were extracted to further study the effect of the pterostigma on the wing performance. All the pressure distribution has a distinctive shape where the peak and dip are near the leading edge at y/c = 0, while they curve downward at the trailing edges at y/c = 1.0. The upstroke and downstroke pressure distributions are different where the upstroke pressure distribution curves have lower magnitude peaks and dips, while the downstroke pressure distribution curves have higher magnitude peaks and dips at the leading edges where y/c = 0. It was observed that the pterostigma weight changes the mid-and trailing edge pressure distribution profile to become positive pressure at the lower surface and neg-ative pressure at the upper surface of the wing sections compared with the 0% pterostigma case. This may be caused by the inertial forces induced from the additional pterostigma weight, which deforms the pitch of the wing sections, especially for those near the tip to deflect more incoming wind, which results in the highly positive pressure zone. The presence of pterostigma load also increases the magnitude of pressure at the leading edge, which explains the increase of performance on the wing. Figure 7 shows the pressure distribution plot of the blade sections, where y is the chordwise distance in meters and c is the chord length in meters.

von Mises stress on wing model
The results show that the von Mises stresses were all below the yield strength of the material of the wing, which is 25 MPa. From all the von Mises stress contour plots, the maximum stresses are concentrated in the wing root region near the leading edges. This is expected as the joint was located nearby the maximum stress region. It was observed that the addition of pterostigma increases the maximum stress in the wing model root region. The maximum von Mises stress during upstroke and downstroke is 4.395 and 4.539 MPa, respectively, for 11% pterostigma weight. These von Mises stresses were slightly higher than the 0% pterostigma weight cases, which can be mitigated by simply increasing the thickness of the wing root. The shape of the wing can be further optimized to accommodate more stress capacity by tapering the tip or corrugating the wingspan. Figure 8 shows the von Mises stress distribution contours of wing models with and without pterostigma loading.

Effect of different locations of pterostigma load
The effect of pterostigma load location on the lift and drag coefficients of the flapping wing model was investigated by trial and error. Simulation cases involving various pterostigma locations according to the planform map, as shown in Fig. 1, were executed. For the initial setting, the 8% pterostigma weight was chosen for the location test. The simulation results show that the maximum lift coefficient, as well as the lift to drag ratio, was the highest when the pterostigma load was at x = 47 mm and y = −2.5 mm. A single case is adequate to represent the pattern of lift and drag coefficients. The rest of the figures show only the peak values for comparison purposes. Figure 9 shows the lift and drag coefficients for 8% pterostigma weight at the optimized location. The maximum lift to drag ratio near 0.5 s was 17.93. It was observed that the amplitude of the lift coefficient increases gradually to a constant value for the case of the optimized pterostigma location, while the rest of the locations have a decrease of the lift after a peak is reached when approaching 0.5 s. For the case of the drag coefficients, the peaks of the graph were fluctuating, and the degree of fluctuation is the highest for the optimized case. Figure 10 shows the normalized peaks of lift and drag coefficients against time at different locations. Despite having the high drag coefficient, the lift to drag ratio of the optimized case is the highest; thus, it is concluded that the location farthest from the wing model joints, which is also the wing tip, was the best place to locate the pterostigma loadings. It was observed that the changes in the lift to drag ratio by the pterostigma location along the  chordwise direction have less effect than those along the spanwise direction. The percentage difference of the highest to lowest lift to drag ratio for the chordwise direction is 4.88%, while it is 35.91% for the case of spanwise direction. Thus, it is concluded that the spanwise direction determines the efficiency of the pterostigma load on the wing model performance. The sudden drop in the lift to drag ratio along the chordwise direction may be caused by the flow separation, which is predicted by the simulation models. The flow separation may occur in the mentioned region, which reduces the lift while increasing the drag. Thus, there is a sudden drop in the lift to drag ratio along the chordwise direction. Figure 11 shows the lift to drag ratios for various pterostigma locations along spanwise, x/L, and chordwise, y/c, directions.

Effect of different amounts of pterostigma load
After the determination of the optimized location of pterostigma, different weights of the pterostigma loading were applied at x = 47 mm and y = −2.5 mm. The simulation results show that the change of pterostigma loading affects the lift and drag coefficients of the wing model. It was predicted that at pterostigma loading of 11% of wing weight, the maximum lift coefficient is highest. The trend of the lift and drag coefficients is similar to the previous section, where the lift coefficient has a smooth increment of flapping amplitude while the drag coefficient has fluctuating peaks. Figure 12 shows the lift and drag coefficients against time for 11% pterostigma weight at an optimized location. Figure 13 shows the normalized peaks of lift and drag coefficients against time for pterostigma loads of 8%, 9%, 10%, 11%, 12% and 14% of wing weight. It was observed that the highest lift to drag ratio of 18.65 occurs when the pterostigma weight is 11% of the total wing weight. This optimized pterostigma load produces a lift to drag ratio, which is 129.16% higher than the wing model without pterostigma load. Thus, this pterostigma loading of 11% of wing weight is adopted for the next section. Figure 14 shows the lift to drag ratios for various pterostigma loads.

Effect of changing wing model elastic modulus and density
The percentage change of the original elastic modulus, 1.1 × 10 9 Pa, was used to test the effect of different elastic modulus values on the lift and drag coefficients. The proposed percentage change of the original elastic modulus ranged from −30% to +30%. Figure 15 shows the normalized peaks of lift and drag coefficients against time for elastic modulus E. It was observed that only the curve of 0% change of the original elastic modulus has the highest lift coefficient and also exhibits a smooth increase starting from 0 to 0.5 s. The rest of the curves were having peaks initially, followed by a gradually decreasing slope toward 0.5 s. This shows that the pterostigma weight has to be adjusted to suit the different elastic modulus values in order to obtain an optimized lift. The drag coefficient of the original elastic modulus,  however, remains the highest. Despite the high drag coefficient, the lift to drag ratio for 0% change of the original elastic modulus remains the highest with the value of 18.65. For the case of densities, similar to the case of the elastic modulus, the original density of 950 kg/m 3 was adjusted by percentage change ranging from −30% to +30%. Figure 16 shows the normalized peaks of lift and drag coefficients against the time for density ρ. The trend lines of the case of densities were similar to those for the elastic modulus where the original density values produce the highest lift, highest drag and highest lift to drag ratio. This also proved that the density needs to be adjusted along with the pterostigma weight in order to achieve optimized flight performance. The cross-sectional area A is 0.5 mm 2 , while the frequency f is 25 Hz. From Eq. (9) , the dimensionless constant D c graph was plotted. It was found that when D c is 0.0001032, the lift to drag ratio is the highest with the value of 18.65. The dimensionless constant D c graph can serve as a future guide for the users to determine the optimized elastic modulus and density for the proposed wing model design. Figure 17 shows the lift to drag ratio for various values of elastic modulus E (%), density ρ (%) and dimensionless constant D c .

CONCLUSIONS
The results show that the effect of pterostigma on the dragonflyinspired wing model has the potential to be used in future MAV wing design. The findings are summarized as follows: (1)The optimized pterostigma location was predicted to be x = 47 mm and y = −2.5 mm. This shows that the pterostigma affects the flight performance of the wing model. The best location is at the tip of the wing. (2)The optimized pterostigma load is 11% of the wing model weight. This optimized pterostigma load produces a lift to drag ratio 129.16% higher than the zero pterostigma load wing model. This implies that certain pterostigma load produces the highest lift of the wing model.

Figure 14
Lift to drag ratios for various pterostigma loads (% of wing weight).
(3)The peak lift to drag ratio is 18.65 at D c of 0.0001032. The D c graph in this study can be adopted by future users in the dragonfly wing model design.   (6)The von Mises stress contour plots on the wing show that the critical stress is located near the hinge regardless of the pterostigma loads, and upstroke and downstroke motion. This implies that the pterostigma load increases the maximum stress. This problem can be mitigated by increasing the thickness of the blade section near the wing root.

NOMENCL ATURE
A = cross-sectional area of the wing base A m = amplitude A r = reference area c = chord length C D = drag coefficient C L = lift coefficient Co = Courant number D c = dimensionless constant E = elastic modulus E k = turbulence kinetic energy generation by velocity gradient E ω = specific rate generation by velocity gradient f = frequency F D = drag force F k = external force source term for kinetic energy F L = lift force F s = the external and gravitational forces F ω = external force source term for specific rate I k = kinetic energy diffusivity I ω = specific rate diffusivity j = mass flux K = elasticity p = pressure t = time v = flow velocity V = velocity of airflow V m = velocity magnitude Greek symbols t = time step x = mesh node interval θ t = angle of the flapping wing tip ε = elastic strain vector ρ = density σ n = stress vector (n = x, y, z, xy, yz, xz) τ = stress tensor ∇ = divergence

ACKNOWLEDGE MENTS
This research was financially supported by the MOSTI-Science