Influence of Area Distribution Along the Span Direction on Flapping Wing Aerodynamics in Hover Based on Numerical Modeling Analysis

This study performed a numerical investigation to explore the effect of span-wise area distribution on flapping wing aerodynamics in hovering configuration using an in-house developed flow solver. The span-wise area distribution was defined using a morphological beta function and the flapping motion was set to a sinusoidal movement manner. Results show that moving the area distribution toward the tip region can generate more lift and simultaneously expense more power, whereas an optimum area distribution (${\rm{\bar{r}}}_1$ = 0.45) was observed because of its aerodynamic efficiency. In general, the temporal profiles of aerodynamic forces are slightly sensitive concerning span-wise area distribution rather than the peak and mean force. Detail flow structure visualization illustrated that the flapping wing locomotion produces complex spatial and temporal vortex structures, including vortex generation, development, and shedding of leading-edge vortex, trailing edge vortex, and root vortex. For flapping wing with a larger area on the tip is in principle capable of enhancing the vortex strength, particularly for the leading-edge vortex which dominates the lift generation during flapping motion. Meanwhile, the smoother profile bounded by the tip and leading edge is beneficial for stabilizing the leading edge and tip vortex.

be able to produce enough lift to stay aloft or forward fly.The mechanism whereby they achieve flight must benefit from unsteady flows interacting with the dynamically changing wing surfaces.Till now, five unsteady high-lift mechanisms associated with flapping wing aerodynamics have been well-revealed [1], viz., Wagner effect, leading-edge vortex (LEV), delayed stall, wing rotation, wake capture, and clap-and-fling.Using the aforementioned collected knowledge related to flapping wing aerodynamics, kinematics, control strategy, etc., some flapping wing prototypes have been successfully developed with examples such as DelFly [2], Robobee [3], Microbat [4], etc.
The wing geometry, generally regarded as the morphological planform or outline, and kinematics of flapping wings are crucial parameters that dominate the aerodynamic performance, particularly in hover configuration.Over the years, extensive research has been conducted on kinematic parameters such as flapping frequency, stroke amplitude, pitching angle, tip trajectory, and the phase offset between pitch and stroke angles.Alternatively, the effect of flapping geometry, such as spanwise/chordwise area distribution, total area, location of the pivoting point, and aspect ratio (AR) was seldom reported.Looking back into nature there are millions of flying insects with diverse wing shapes, therefore, it is difficult to find normalized morphological shape factors for statistically describing the wing planform.The mechanical importance of the moments of the wing area is emphasized, which depends on the shape of the wing [5].His pioneering study approximated nearly all flying insects' wings by a semiellipse, the approximations were quite suitable for an initial investigation.Later, Ellington [6] of the University of Cambridge claimed that by comparing the measured distribution's values to those of an appropriate analytical function, one can derive the precise form of a distribution from its shape characteristics.He then defined a beta distribution (see later in Section II) compared with the wing area distribution, i.e., the planform of the wing, and a very good agreement was revealed.The rules of shape that relate the higher radii to the first radius may be used to further compress the beta distribution to a coefficient for just one parameter, giving rise to a useful tool for approximating the planform of a wing utilizing only its centroid of the surface.The air pressure distribution of the flapping wing under varied Reynolds numbers and the variation trend of pressure distribution on the wing surface with different Reynold numbers was observed, while the influence of Reynold numbers on the force production ability was also investigated [6].And they tried to optimize the morphological parameters of rectangular flapping wings depending on the simulation cases.Their results show that the adaptive amendment methods based on the air pressure distribution improve the aerodynamic performance which is consistent with the evolutionary characteristics of forewings of some large insects.Later, the LEV of hovering artificial flapping wing in different AR, Reynold numbers, and Rossby numbers by solving unsteady incompressible Navier-Stokes equations [7].In their study, four variations of AR and two orders of magnitude variations of Reynolds number were set to research.In computational fluid dynamics (CFD) solvers, the Navier-Stokes equation is pivotal.It governs fluid flow behavior, allowing CFD simulations to predict velocity, pressure, and other flow properties.Accurate solving of this equation is vital for modeling complex fluid dynamics, such as those involved in flapping wing locomotion.Results show that the best AR of the flapping wing was consistent with the results observed in natural flapping insects at different Reynolds numbers, while it was found that decreasing the Rossby numbers with a certain range would improve the spanwise flow, leading to maintaining the intensity of LEV, inducing a higher force production.Stewart et al. [8] developed a set of multiobjective optimization methods based on gradient optimization solver and error constraint.By modeling the optimization characteristics of the flapping wing.The aerodynamic performance of rigid flapping wings in forward flight and hovering state was investigated.Results show that the aerodynamic performance of the flapping wing surface was improved at the same Reynolds number through the optimization solver.Shahzad et al. [9] tried to find a collaborative solution to solve the fluid structural interaction problems of flapping wings, it contains two solvers, the immersed boundary method (IBM) is designed for solving the three-dimensional (3-D) Navier-Stokes equations, and the finite-element method solver is for structural deformation which was used to provide boundary conditions for the flow field.The results show that the absolute force of the flexible flapping wing with a larger AR is smaller than that of the rigid flapping wing, but the aerodynamic efficiency is higher.This difference is due to the deformation of the flexible wing during the pitching motion, which leads to a decrease in the equivalent pitching angle.When the AR is larger (AR = 6), the aerodynamic efficiency of the flapping wing surface area distribution near the wing tip is higher than that of the wing root, but when the AR is smaller (AR = 3), the opposite conclusion will be obtained.
It is also interesting to see from the archived document that the aerodynamic performance related to the wing planform or shape has contradictory results.In the experimental studies from Phillips et al. [10] and Wang et al. [11], similar flow topology for a variety of shapes of flapping wings have been reported, however, there is no force, particularly the temporal force variation during flapping cycles provided.Some studies [12], [13], [14], [15] reported that there is a less than 5% difference in the instantaneous lift coefficient within 10 wing planforms based on a fruit fly's wing.Luo and Sun [13] and Wilkins [14] conducted a comprehensive investigation involving both quasi-steady and computational analyses of flapping wings.Their findings revealed that wings with larger outboard areas generate increased lift but at the expense of higher power consumption.Harbig et al. [16] carried out a computational study focusing on four distinct flapping wing planforms: rectangle, ellipse, reverse ellipse, and four ellipses.Three different kinematics were applied-thrips kinematics for a Reynolds number of 12, honeybee kinematics for a Reynolds number of 1134, and 2-angle insect-inspired kinematics at Re = 13 500.The results indicated that as the Reynolds number increased from 10 to 10 000, inertial effects became the predominant factor in force generation.This phenomenon resulted from both pressure differences and the formation of vortical structures.Ke et al. [17] propose a solution for the intricate dynamics of highly coupled nonlinear hovering wingbeats with two degrees of freedom.They employ an extended quasi-steady aerodynamic and inertial forces/moments model and a numerical approach involving ordinary differential equations to achieve a solution.In addition, the study investigates the variable phase offset rule between wing pitch angle and flapping angle, enabling bioinspired flapping wing microaerial vehicles to maintain a high variable angle of attack (AoA).In a related study, Lang et al. [18] examine the influence of wing geometry and kinematic factors on the aerodynamic performance of flapping wing microair vehicles (FWMAVs).The findings emphasize that wing area has the most significant impact on lift and power loading, while a moderate sweeping amplitude with advanced rotation enhances lift.The study provides conceptual recommendations for the development of FWMAVs, highlighting the importance of these factors in optimizing aerodynamic performance.
Besides, the planform of the flapping wing shape, the effect of AR, which also can be regarded as an important morphological parameter of the flapping wing has also drawn considerable attention [19], [20], [21], [22], [23].By enabling researchers to control and manipulate these features in a controlled environment, numerical models and simulations play a crucial role in studying the relationship between morphological features and aerodynamic performance without being constrained by physical experiments.The variation in the lift coefficient for different ARs is insignificant, say less than 10% while for some cases the difference is a noticeable trend.Generally, a low AR was recommended for producing a higher lift coefficient, and the efficiency for different AR wings was barely discussed in detail.
Because morphological characteristics of flapping wings have significant influences on the aerodynamic performance and there is no systematic study to identify the aerodynamic mechanism associated with flapping wing planform.Achieving appropriate surface roughness or texture to simulate natural wing characteristics is difficult when trying to replicate specific morphological traits of engineering flapping wing designs.Other difficulties include ensuring the scaling is exact and the structural flexibility is replicated.This study is dedicated to investigating the effect of shape parameters, i.e., span-wise area distribution (SWAD) on the aerodynamic performance of flapping wing locomotion utilizing high-fidelity numerical simulation.The aerodynamic performance and flying behavior of the wing are impacted by the SWAD of the flapping wing movement, which is crucial.With variations in wing width or form along their span, insect wings frequently have intricate SWADs.Lift, drag, and stability during flight are all impacted by this distribution.To help with takeoff and maneuverability, tapered wings with broader bases produce greater lift at the wing root, while wings with narrowing tips produce less generated drag.Overall, the distribution of span-wise area is a crucial adaptation that improves the effectiveness and maneuverability of insect flying systems.Accurate wing geometry, precise fluid dynamics representation, quality meshing, suitable boundary conditions, numerical solvers, vortex modeling, time integration, force calculations, validation through comparisons, and postprocessing for analysis are essential elements of high-fidelity numerical simulations of flapping wing locomotion.Together, these elements produce a simulation that accurately captures the intricate aerodynamics during movement.
This article is organized as follows.In Section II, the flapping wing model and numerical methodology are introduced and subsequently validated in Section III.Section IV will present the results and discussion followed by a conclusion given in Section V.

A. Physical Model and Kinematics
As aforementioned, the platform, i.e., the morphology of the insects in nature varied significantly and thus required some shape parameter to quantify the morphology of flapping wings.In the work of Ellington [6], it is stated that the SWAD of the insect wings can be statistically quantified by using a polynomial distribution and a beta distribution.Shape parameters are essential for understanding and representing the diverse morphologies of insect wings.They quantify features like AR, camber, and wing shape, offering insights into flight performance and evolutionary adaptations.The formulations to generate wing planform shape from Ellington [6] are written as follows: where c and r are the chord length of the flapping wing and the distance between the wing root in a spanwise direction, which are normalized by the mean chord length c and wing half-wing span R, respectively.r2 is the nondimensional radii of the first and second moments of the flapping wing area.B (p, q) is the beta function which is defined by p and q, computed from r1 and r2 .As documented in Ellington [6],r 1 generally fall within the range [0.43, 0.563] for natural flying insects.Swing wing aerial drone (SWAD) team developed a drone with a dynamic wing design.This drone may be launched and operated by a single person.The SWAD is used to evaluate wing attributes at various sweep angles to determine if a two-position wing is appropriate for a small aircraft.To evaluate the effect morphology in terms of SWAD on the aerodynamic characteristics, five beta distributions, i.e., r1 = 0.4, 0.45, 0.5, 0.55, and 0.6 are considered as shown in Fig. 1.The AR of the tested flapping wing stayed constant at 3.0 throughout the study.The half-span and mean chord of the flapping wing is 10 cm and 3.33 cm, respectively, and thus resulted in the wing area equal to 33.3 cm 2 .The kinematics of flying insects and vertebrates in nature temporally varied in flapping in sweeping planes and pitching as shown in Fig. 2. The kinematics in terms of sweeping and pitching of the flapping wing in this study is defined as follows: where θ and α are the sweeping angle and pitching angle of the investigated flapping wing, respectively.The sweeping strokes θ 0 and pitching stroke α 0 are set to 45°(π/4) throughout the study.Keep in mind that the sweeping axis is at the wing root, whereas the pitching axis is at the leading edge.The pitching angle (θ) and sweeping angle (σ) are crucial in determining the movement of a flapping wing during a cycle of motion.The pitching angle balances lift and thrust by tilting the wing relative to airflow, while the sweeping angle represents the lateral movement of the wing as it flaps.These variations enable the wing to generate lift, thrust, and maneuver effectively during flapping-wing flight.The leading edge, the part of the wing that first contacts the air and is the foremost edge of an airfoil section, has a maximum curvature and minimum radius.Moreover, the thickness of the flapping wing is set to 0.02c.The use of a single AR (AR = 3) and fixed flapping and pitching amplitudes of 45°simplifies analysis, improves comparability between experiments, relates to biological examples, enables sensitivity analysis, and accommodates practical constraints.These parameters may be selected based on research objectives and the need to isolate specific factors to develop a better understanding of aerodynamic principles.However, they may change depending on the objectives and constraints of the study.The pitch axis in flapping wing aerodynamics, representing the rotational axis during pitching motion, is revealed to be positioned not at the leading edge but at approximately 1/4 of the wing's chord, as depicted in kinematic schematics illustrating pitch (θ) and sweep (σ) angle variations within a flapping cycle.The chord, defined as the straight-line distance between the leading and trailing edges, places the pitch axis a quarter of the way from the leading edge toward the trailing edge.This distinctive placement holds significance in the study of flapping wing aerodynamics, influencing the understanding of aerodynamic performance in both biological flyers like birds and insects and the design of biomimetic flapping-wing robots.Precise knowledge of the pitch axis location is crucial for analyzing the complex motion of flapping wings and for designing efficient and maneuverable flapping-wing systems.

B. Numerical Methodology and Setup
Employing a custom CFD solver, the unstable viscous fluid pattern of the flapping wing is quantitatively calculated.When determining a low Reynolds number for flapping wing locomotion, one compares the characteristic size (usually the length of the wing chord) and the speed of the flapping motion to the kinematic viscosity of the fluid.Small values of this ratio indicate low Reynolds number flow, which is frequently characterized by laminar, viscous flow characteristics.CFD simulation software solves complicated flow equations using numerical techniques.In CFD, there are two kinds of solvers: pressure-based and density-based.For low-speed flows, Mach numbers less than 0.3, and incompressible flows, pressure-based solutions are implemented.Density-based solutions are employed for highspeed and compressible flows.Advanced techniques such as high-order spatial and temporal discretization, adaptive grid refinement, turbulence modeling (such as RANS or LES), and aeroelastic coupling are used in the in-house numerical method for modeling the turbulent airflow around a flapping wing.These characteristics allow for a highly accurate and detailed simulation of the wing's unsteady motion and interaction with the surrounding fluid.With the help of these properties, complex flow behavior connected to flapping wing locomotion can be accurately and effectively simulated.A low-mach-number preconditioning methodology is used to widen the code's use in the low-speed realm, enabling the acquisition of relevant convergency, security, and precision findings.The code resolves the Navier-Strokes equation via a finite volume approach.In CFD solvers, the low-mach-number preconditioning technique is used to improve the stability and accuracy of simulations in situations where fluid velocities are low relative to the sound speed.It is frequently used to simulate incompressible flows in low-speed aerodynamics, interior pipe flows, and heat transfer issues.This method guarantees effective and trustworthy outcomes in a variety of engineering and scientific simulations.Furthermore, there is the dual-time-stepping method, which involves introducing a pseudo-time to move the outcome along at every physical time layer.The dual-time-stepping method divides the time integration process into two steps: pseudo-time and physical time, improving simulation stability and accuracy.Equations are iteratively solved in the pseudo-time step to arrive at a steady-state solution inside each physical time step.This decoupling enables more steady convergence, particularly in transient or complex flows, which ultimately enhances the accuracy and dependability of simulations.It is possible to write the formulas that govern with a preconditioned pseudo-timedependent component as follows: in which τ and t are the pseudo-and physical time, respectively.W and Q are the conservative variables and primitive-variable vectors, whereas F(W) and FV are the convective and viscous fluxes.ν g is the contra-variant velocity of the boundary S(t) of the control volume Ω(t) and Γ is the low Mach preconditioning matrix.Due to the relatively low Reynolds number of the flapping wing movement, the simulations successfully resolve the laminar flow.The decision to solve for totally laminar flow is influenced by the comparatively low Reynolds number in flapping wing movement by making it a meaningful factor.At lower Reynolds numbers, viscous effects predominate, and laminar flow may continue over a larger area of the wing's surface.Because laminar flow behavior greatly affects aerodynamic performance and lift generation, it is essential to precisely characterize and anticipate it when simulating flapping wing action.
It should be noted that a convergence condition is established (residual down to 10-3) to put an end to the computation.Our earlier work [24] has further information about the CFD solver.The present CFD solver is combined with a hierarchical overset grid method to instantly simulate the movement of the flapping wings.In the computational setup, overset grid schematics are essential because they allow for dynamic grid interactions.They simulate realistic motion by allowing the flapping wing grid to move with the background grid's fixed position.The simulations are more reliable because of this dynamic capability, which guarantees a precise portrayal of intricate aerodynamic interactions during the flapping process.The overset-grid-schematics used throughout the calculation are shown in Fig. 3 together with the computational grid zones.As can be seen in Fig. 3(a), there are two meshes in total, including the base mesh and a submesh for the flapping wing.It is crucial for obtaining accurate simulations to use two meshes in the computational setup-one for the flapping wing and another as a backdrop mesh.While the background mesh serves as a stable frame of reference for computations involving computational stability and boundary conditions, the wing mesh enables accurate modeling of the wing's geometry and aerodynamics.This results in a more accurate depiction of the flapping motion and its interactions with the surrounding fluid.The global size of the computational domain is 30×20 * 30c and meshed with structured hexahedron, whereas the boundary layer is refined to Y+ = 1.0 to guarantee computational accuracy.The grid cell amount of the flapping wing and background meshes for all the test cases are about 2 and 7 million, respectively.The overset mesh in CFD is a method used to simulate complex fluid flow problems like flapping wings.It involves dividing the problem domain into separate grids and determining the overlap between them.The overset approach uses two unstructured meshes that simulate the displacement of a flapping wing: the background mesh representing the wind tunnel and the overset mesh containing the wing surfaces.The computational setup consists of a primary solution domain and one or more secondary grids that combine the primary domain.At each time step, the overset mesh is determined based on the position of the flapping wings and other moving bodies.Information is then exchanged between these grids in the overlapping regions, and the CFD solver computes the fluid flow.This approach enables the simulation of dynamic and unsteady flow phenomena accurately.
In this study, the most representative parameter that can provide important information on the effect of the SWAD of flapping wing is the lift coefficient CL, power coefficients CP, a, and propulsive efficiency η (the ratio of lift to power) defined, respectively, as follows: where Ūtip is the mean tip velocity within one flapping cycle and S wing is the projected area of the flapping wing.

III. SOLVER VALIDATION
The proposed numerical methodology has been validated on a 3-D wing heaving in a sinusoidal manner, where the incoming freestream-based Reynolds number is 10 000.On an annular wing with a NACA0012 cross-section as well as a semi-AR of 2.0, their studies [24] were carried out.Both aerodynamic force and flow structure were measured through force sensors and 2-D particle image velocimetry (PIV) techniques.PIV, which tracks moving particles, calculates the fluid's velocity.The two primary PIV techniques are laser-based PIV and Doppler PIV.Calibration and other methods, such as phase-Doppler anemometry and laser extinction, can be used to measure the particle density in PIV.In the current investigation, the case Strouhal number (St = fc/UÝ) with a 20°AoA has been selected for validation.Validation challenges are tackled by calibrating instruments, reducing errors, conducting repeated trials, comparing with established benchmarks, and refining methods to enhance accuracy and reliability.The grid cells for the airfoil, the backdrop meshes, and the computational domain as depicted in Fig. 4(a) comprise around two million.The reliability and accuracy of the validation process are greatly improved by grid cell refining for the backdrop and airfoil meshes.In areas of interest, such as those closest to the wing surface and vortex shedding zones, it enables greater resolution, enabling more accurate modeling of flow behavior.This enhancement guarantees that the simulation accurately depicts the aerodynamics, improving the validity and realism of the process.The grid on the heaving airfoil is refined at Y+ = 1.0.Fig. 4(b) shows the overset scheme at a certain time.Note that, during the simulation, the k-ω SST turbulence model was employed to close the governing equation.The k-ω SST turbulence model was selected to close the governing equations during the simulation because it combines the benefits of the k-ε and k-ω models, providing accurate predictions for a wide range of turbulent flows.This choice improves the model's ability to capture the complex turbulence effects seen in flapping wing locomotion, enhancing the accuracy of the simulation results.
Table I compares the aerodynamic force coefficients among the computational and experimental data, showing that our computational model undervalues both lift and drag while being within a suitable range.
In investigating aerodynamic forces, the methodology typically involves setting up an experiment with a scaled model, Fig. 5. Flow structure comparison between the PIV measurement [24] and computational results.
equipping it with sensors, and collecting data during flapping motion.Analysis and visualization techniques are then used to study forces, flow patterns, and vortical structures.Numerical simulations and comparisons to theoretical models provide additional insights, allowing for a comprehensive understanding of aerodynamics.The process concludes with summarizing findings and making recommendations for further research or practical applications.Similar numerical results were also given by Tay et al. [24] using the IBM.Except for the aerodynamic forces, the flow detail in terms of vortex topology around the heaving wing is illustrated in Fig. 5.As can be seen, the numerical method in this study can accurately capture the generation, translation, and shed of both leading and trailing edge vortices, moreover, the computed vortex pattern has good agreement with the experimental result.

IV. RESULTS AND DISCUSSION
Except for the effect SWAD on the aerodynamic performance of the flapping wing, three flapping frequencies, i.e., f = 7, 14, and 28 Hz are also considered, corresponding to the tip-velocity-based Reynolds number around 5000, 10 000, and 20 000, respectively.In investigations of flapping motion, the Reynolds numbers (5000, 10 000, and 20 000) are essential for identifying the various flow modes.These numbers are chosen by researchers to accurately match particular flapping frequencies, ensuring a realistic depiction of aerodynamic or hydrodynamic forces.This part addresses and discusses the calculated aerodynamic forces and the flow framework, such as the surface-pressure gradient and vortex morphology.

A. Aerodynamic Forces in Hover
In the majority of the archived articles investigating the flapping wing shapes, the aerodynamic forces have been examined by plotting the averaged values of aerodynamic forces.To better understand the mechanics of flapping flight, researchers are examining the aerodynamic forces generated by flapping wings.These forces, including as lift and drag, are important because they influence the wing's capacity to provide lift for flight and control its trajectory.Designing effective biomimetic flying devices and gaining insights into the intricate aerodynamics of flapping flight depends on an understanding of these forces.However, the temporal variation in force generation during flapping cycles also requires a deep insight in terms of its performance.The impact of SWAD and flapping frequency on the force generation of a flapping wing, with a focus on the computation of lift forces for various test cases.The temporal variation of lift generation throughout a flapping cycle is illustrated in Fig. 6.In each cycle, two distinct lift force peaks are observed: the first peak occurs around 3/4 of the forward stroke, while the second, with a lower magnitude, occurs around 3/4 of the backward stroke.The findings reveal a nuanced pattern in lift generation, highlighting specific timing and magnitudes during the wing's flapping motion under the tested conditions.The flapping cycle having two lift maxima improves the wing's aerodynamic efficiency.These two peaks work together to provide lift more steadily and smoothly, lowering the possibility of stalling and increasing overall efficiency during flapping action.It is also logical to see that the flapping locomotion in this study can generate positive lift throughout the flapping cycle.Similar force variation for different SWADs indicates that moving the area along a spanwise direction is not sensitive to the force variation topology but only influences the magnitude.The research is centered on assessing the combined impact of SWAD and the Utip velocity-based Reynolds number on lift generation throughout a single flapping cycle.SWA and SWD are key parameters determining wing motion and, subsequently, impacting aerodynamics and airflow patterns.Simultaneously, the Utip velocity-based Reynolds number, denoted by the hues red, green, and blue, acts as a measure of variations in flapping frequency.To understand the complex equilibrium between enhanced lift and heightened drag with more SWA, research evaluates lift generation across various SWAD situations.The temporal history of lift across SWAD clusters is graphically represented with unique hues, allowing for the study of lift patterns across time.The varied hues correlate to different Reynolds numbers, providing a visual representation of how flapping frequency influences the dynamics of lift production.
Because higher Reynolds numbers signify a change from laminar to turbulent flow, increasing the Reynolds number typically results in enhanced lift generation.Reduced drag and increased lift are produced by turbulent flow, which produces more vortices and less airflow separation.However, there can be outliers at particular Strouhal numbers (SWAD), such as 0.4 and 0.6, where the interaction between vortex shedding and wing flapping frequency may disturb lift generation, leading to fluctuations or reductions in lift.Fig. 7 shows the variation of the period-averaged lift generation (herein denoted by CL), power consumption (CP, a), and propulsive efficiency (CL/CP, a) of the flapping wings with different SWAD flapped at different Reynolds numbers.As seen, the lift generation increases with increasing ther 1 , i.e., moving the wing area toward the wing tip, which is as expected.Another point should be noted that a higher Reynolds number, defined by the mean tip velocity, can increase the lift generation for some cases, however, it is not the case for r1 = 0.4 and 0.6.It is also surprising to see that the lift generation is not sensitive concerning the Reynolds number, or flapping frequency.The presence of flow separation, boundary layer effects, and complex interactions between vortices that can lessen the expected impact of these parameters on lift production are factors that contribute to the surprising observation of insensitivity in lift generation to Reynolds number or flapping frequency.A similar tendency is also found in the power consumption coefficients that different SWADs significantly influence the power consumption during flight, whereas moving the area outside requires more power input.The main determinants of power consumption in flapping-wing flight are wing kinematics (flapping frequency, amplitude, and stroke pattern), wing shape, airfoil properties, and the aerodynamic forces produced, which together determine the energy needed to maintain and control flight.Moving to the propulsive efficiency as seen in Fig. 7(b), the propulsive efficiency peaks at r1 = 0.45, Re = 20 000.Besides, increasing the flapping frequency can augment the propulsive efficiency, i.e., generating more lift with identical power input.To conclude on the effect of SWAD on flapping wings, moving the wing area toward the tip can significantly generate more lift and require more power input, while there exists an optimum propulsive efficiency.

B. Effect of SWAD on Flow Topology
The vortical structures at various time instants are visualized using the Q-criterion, which is displayed with pressure coefficient for the identical purpose of describing the pressure gradient as seen in Fig. 8, to obtain a comprehensive understanding of the 3-D flow patterns of the flapping wing in circling arrangement.Using the Q-criterion visualization technique, vorticity, and strain rate are used to determine the presence of vortical structures in the flow field of a flapping wing.As a way to see and understand these patterns, vortical zones are found when the Qcriterion is greater than a certain threshold.As can be observed, the vorticity map shows a variety of intricate vortex formations, with the LEV, root vortex (RV), trailing vortex (TV), and shaded vortex (SV) being the four most noticeable ones.Studying Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.complex vortex structures in vorticity plots reveals fine-scale flow details and vortex interactions critical for understanding and optimizing aerodynamics in flapping wing systems, insights often missed by other analyses.Location on an aerodynamic surface is the primary difference between the LEV, trailing edge vortex (TEV), and RV.To increase lift, LEV develops close to the leading edge.Drag is increased by TEV, which happens at the trailing edge.Spanwise flow is regulated and induced drag is decreased by RV, which arises at the base of a lifting surface.At t/T = 0.0, when it is at the onset of the new flapping stoke, the LEV, TEV, and RV that are generated during the previous outstroke are experiential to shed from the leading edge, trailing edge, and root, respectively.From t/T = 0.0 to 0.1, the flapping wings are starting to move, at this moment the LEV, TEV, and RV shedding from the previous flapping cycle were merged with the newly generated vortices on the flapping wing leading edge, trailing edge, and root.Because it may affect aerodynamic effectiveness and lift production, the merger of previously shed vortices with the fresh vortices produced during the beginning of flapping wing motion is significant.The overall flight performance and stability of flapping-wing systems can be affected by the way this process affects flow stability and vortex interactions.Moving to t/T = 0.2, the leading and trailing edges are getting stronger but still attached to the wing which results in a force lift peak as evidenced in the force variation, see Fig. 1.When the flapping wing moves to t/T = 0.3, the LEV starts to shed from the leading edge, and the trailing edge is already shed in space forming a revolving SV as indicated by SV.At this moment, the aerodynamic lift starts to decline along with the LEV shedding.At t/T = 0.4, the leading edge transverse velocity is much slower than the trailing edge, and thus the shed LEV from t/T = 0.3 and TEV at this moment are interacted and merged on the low-pressure side of the flapping wing.When t/T = 0.4, the vortices on the low-pressure side of the flapping wing interact and merge, which improves lift production.A stronger vortex is produced as a result of the merging process, which promotes lift and aids in flight or propulsion during flapping motion by maintaining low pressure on the wing's upper surface.At t/T = 0.5, when it is at the end of the half-stoke, the leading edge stays at a fixed location with only trailing edge locomotion, and the LEV and TV are shed from the leading edge and trailing edge, indicating there is no lift generation at this moment, seen in Fig. 1.From the t/T = 0.6 to 0.9, another half-stroke occurs, and the flow topology is identical to what happened during t/T = 0.0 to 0.5, this can be also revealed from the symmetric temporal force variation in Fig. 1.At the beginning of a new flapping stroke (t/T = 0.0) and the end of a half-stroke (t/T = 0.5), vortices are typically smaller and closer to the flapping object, whereas at the end (t/T = 0.5), larger vortices have formed and are further downstream due to the accumulation of fluid momentum during the stroke.
To reveal the flow structure associated with different SWADs and flapping frequencies, the flow topology in terms of Qcriterion and surface pressure coefficients distribution at a certain time instant (t/T = 0.8) are plotted in Figs. 4 and 5, respectively.In flapping wing flight, analyzing the flow topology and surface pressure coefficient distribution provides important information about vortex formation, pressure distribution affecting lift and drag, flow separation sites, control systems, and overall flight efficiency and stability.Collectively, these discoveries advance our knowledge of the aerodynamics underlying flying with flapping wings.As seen in Fig. 4, for a given flapping frequency, i.e., the same Reynolds number, the LEV becomes stronger when moving the wing area toward the wing tip.The local chord length of an airfoil or wing directly affects the LEV's strength.Typically, a stronger and more prominent LEV is produced by a greater chord length.The reason for this connection is that a longer chord offers more surface area for the airflow to interact with, increasing the vorticity and circulation close to the leading edge.As a result, a stronger LEV aids in increasing lift and enhancing aerodynamic performance in wing systems that flap.Such phenomenon can be explained by two reasons: first, increasing r1 end with a relatively smaller local chord length which allows a shorter traveling distance for an easier LEV formation, particularly for a slower root tangential velocity; second, moving the area toward the wing tip allows a comprehensive interaction between the flapping wing and the flow field in the near field thus resulting in a higher pressure difference between lower and upper surfaces.Looking into the effect of different Reynolds numbers, e.g., Fig. 4(a)-(c), the flow topologies are almost identical concerning the flapping frequency, while the only dominant difference is the vortex strength, this again can be reflected from the aerodynamic force variation in Fig. 1.Except for the leading edge, the tip vortex also resulted from the pressure gradient between the suction and pressure surfaces.The tip vortex contributes to overall lift generation during flapping motion by reducing the wing's pressure on the upper surface, creating downwash, and enhancing lift through the circulation of air around the wingtip.Since the tip velocity reaches its maximum at t/T = 0.8, it is therefore the vortex strength is correspondingly larger.Whereas for a smaller r1 , the tangential line of the root profile has a smaller angle with the leading edge and hence the interaction area between the tip vortex and LEV becomes even more serious, which accelerates the LEV shedding and energy dissipation (see the vortex sequences in Fig. 9 and Fig. 10).While the situation is altered for a larger r1 , the angle between the leading edge and tangential line w.r.t the wing trailing edge becomes bigger and declines the interaction area of LEV and TV so that the LEV can be attached to the surface more easily.When taking into account a bigger 1/r value for the wing shape (showing a more rounded wing), particular changes in flow behavior often include reduced vortex strength, lowered drag, and possibly lower lift because of smoother flow separation and reduced vorticity compared with sharper-edged wings lead to these changes in flow behavior.

V. CONCLUSION
Given the massive variation in flapping wing morphological characteristics in nature flyers.The effect of morphological features on the aerodynamics of flapping locomotion is still unclear.This study is dedicated to investigating the effect of SWAD on flapping wing aerodynamics in hovering configuration through an in-house developed flow solver.Five SWAD wing profiles were defined using a morphological beta function proposed by Ellington at Cambridge University.Conclusions can be drawn as follows.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 3 .
Fig. 3. Computational setup of the flapping wings.(a) CFD mesh zones of the flapping wing and (b) overset mesh schematic shown at a certain time instant.

Fig. 6 .
Fig. 6.Effect of SWAD and Utip velocity-based Reynolds number on the lift generation over one flapping cycle.

Fig. 6
shows the effect of SWAD and Reynolds number on the lift generation over one flapping cycle: (a) lift generation for different SWAD; (b) temporal lift generation for all the test 15 cases, note that five clusters (say different SWAD) are plotted, and the red, green, and blue color represents different Reynolds number, i.e., flapping frequency.

Fig. 7 .
Fig. 7. Effect of SWAD and Reynolds number on the averaged aerodynamic properties.(a) Lift generation and power consumption.(b) Propulsive efficiency.

1 )
Moving the area distribution toward the tip region can generate more lift and simultaneously expense more power, whereas an optimum area distribution (0.45) was observed because of its aerodynamic efficiency.2) The temporal variation of aerodynamic lift is slightly sensitive for SWAD rather than the peak and mean force.3) Detail flow structure visualization illustrated that the flapping wing locomotion produces complex spatial and temporal vortex structures, including vortex generation, development, and shedding of LEV, TEV, and RV.For flapping wing with a larger area on the tip is in principle capable of enhancing the vortex strength, particularly for the LEV which dominates the lift generation during flapping motion.Meanwhile, the smoother profile bounded by the tip and leading edge is beneficial for stabilizing the leading edge and tip vortex.DATA AVAILABILITYThe data underlying the results presented in the study are available within the manuscript.CODE AVAILABILITY Not applicable.

TABLE I AERODYNAMIC
FORCE COMPARISON BETWEEN THE EXPERIMENTAL AND NUMERICAL RESULTS