. Three-dimensional numerical simulations of flexible flapping wings based on the fluid-structure interaction in biological and bioinspired flow have become a vibrant and challenging research topic. The present paper focuses on a parametric study of the aerodynamic performance of a bionic flexible clapping wing. The proposed model deforms the wing in spanwise and chordwise directions based on the six rigid bodies connected along the wing veins using ball links and springs. Unsteady effects of flapping wing micro air vehicles with a double clap-fling configuration are investigated using an air-solid interaction model based on smoothed particle hydrodynamics and rigid multi-body dynamics. A validation experiment determined the convergence conditions and computational model accuracy. The proposed numerical model is evaluated in terms of flexible variation law and aerodynamic performance. The results indicate that the flapping frequency, angle of attack, and wind velocity significantly influence the lift. Furthermore, increasing the frequency will monotonically expand the maximum and time-averaged lift curve values. When the angle of attack is less than 30◦, the influence on the time-averaged and maximum lift is proportional to the angle of attack. When the angle of attack is larger than 45◦, a stall-like condition is detected. To broaden the applicability of the present findings, a dimensionless parameter, reduced frequency, is defined, and its influence on the maximum and time- averaged lift is investigated. This parametric study shows that as the reduced frequency increases, the maximum and time-averaged lift increases and then decreases. The present study could reach a modeling framework that better explains the clapping wing aerodynamics. Keywords: flapping-wing micro air vehicle; double clap-fling configuration; smoothed particle hydrodynamics; fluid-structure interaction; experimental measurement; three- dimensional flexible clapping wing. Submitted to: Simulation Modelling Practice and Theory Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 2 1. Introduction In nature, several flying and swimming species, including birds, insects, fish, and marine mammals, have evolved over several hundred million years to hunt, avoid predators, reproduce, and execute cluster motions. These creatures have undergone a lengthy process of environmental adaptation, and natural selection [1]. The research on flying biomechanics has a multifaceted significance. On the one hand, scientists must comprehend how flight mechanics interact with biological physiology, ecology, animal behavior, and evolution. On the other hand, engineers must master biomechanical technology to develop bionic machines. Moreover, the fast-growing field of biomimetics requires specialized support from its disciplines. Bionic machine development has advanced at a breakneck pace in the last few decades, with numerous significant advancements in designing bionic birds, insects, and fish. Numerous prototypes have been built over time, with several developed into products for national defense and commercial use [2, 3]. Amongst these bionic machines, the flapping-wing micro air vehicle (FWMAV) is the most representative. The flapping motion of the FWMAV wings comprises three parts: first, the sweeping motion driven by the gear and hinge mechanism, which can vary in frequency but is constrained in amplitude; second, the small amplitude heaving motion occurring due to the bending of the leading edge; third, the passive pitching motion of the wings during a stroke. FWMAV with a clapping- wing, called a clapping-wing micro air vehicle (CWMAV), is considered to have better aerodynamic efficiency. The fluid-structure interaction generated by flapping motion is an important aspect of flapping flight mechanics. Changing the wing geometry and kinematics parameter can positively influence the aerodynamic, inertial, and elastic forces [4, 5]. Owing to the complexity of flapping-wing propulsion and its association with unsteady aerodynamics, FWMAV and CWMAV were designed and studied experimentally and empirically. Particle image velocimetry (PIV) has been the primary experimental approach for quantitative visualization of unsteady effects and is an effective tool for characterizing flapping-wing propulsion systems. Ristroph et al. [6] observed insects such as Drosophila using high-speed cameras. Bross et al. [7] examined the flow field of a vortex system created by the leading edge vortex (LEV) and tip vortex (TV) of a wing with rotation and revolution motion. Percin et al. [8] employed tomographic PIV technology (suited for rotating flat plates with a Reynolds number of less than 10,000) to investigate the 3D flow field and unsteady force experimental investigation using flapping and rotating flat plates under impact conditions. In addition to the aforementioned experimental approaches, developing compu- tational models that can simulate the unsteady characteristics of the fluid-structure interaction involved in flexible flapping-wing propulsion is important. Such models would contribute towards identifying the key high-lift mechanisms behind the flexible flapping-wing propulsion, designing and optimizing the flapping-wing, and providing aerodynamic control strategies. Furthermore, well-developed computational models can Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 3 evaluate the aerodynamic forces at significantly reduced costs compared to experimental designing, manufacturing, and testing [9, 10, 11, 12, 13, 14]. To accomplish a 3D simulation of the flapping wing, Liu et al. [15] solved the whole three-dimensional N-S equation on the structural grid for the first time, confirming the critical function of span-wise flow over the stabilizing LEV. Later, other methods for solving the flow field of 3D unsteady incompressible N-S equations emerged, including the virtual compression method developed by Sun and Tang [16], the finite element method = by Ramamurti and Sandberg [17], and the immersed boundary method by Mittal and Iaccarino [18]. Kruppa’s unsuccessful validation test of the Kasper wing [19] had encouraged scientists to recognize the importance of performing a more detailed analysis of flexible wings, which presents the importance of flapping wings. Ifju et al. [20] compared the relationship between the flexible and rigid lift curves under different angles of attack (AoAs) and successfully adapted the flexible wing to FWMAV. Combes and Daniel [21] examined the wing’s bending condition and degree of bending and discovered that the wing’s deformation is determined by the inertial force. Tang et al. [22] determined the effect of the Reynolds number on the aerodynamic force of a flexible wing. Shyy et al. [23] investigated the TV effect in flexible wings. Tay et al. [14] studied CWMAV and discovered that higher lift occurs when the aircraft’s pitch angle is increased. They also discovered that thrust increases with increasing span- wise flexibility but increases and then decreases with increasing chord-wise flexibility. Therefore, a combination of span-wise and chord-wise flexibility is considered the best layout for aerodynamics. Fluid-structure interaction (FSI) is a core problem in flapping wing dynamics. The meshless Lagrangian smoothed particle hydrodynamics (SPH) method, first proposed by [24], is well-suited for solving fluid impact problems involving free-surface fluid- fluid or fluid-structure interaction, even with large deformations and nonlinear material behavior [25]. The SPH approach partitions the continuous medium into particles and solves the fluid’s governing equations on each particle using a Lagrangian methodology. A particle’s physical properties are estimated by interpolating the values of its nearest neighbors using continuous and differentiable interpolation functions. These functions interpolate fluid velocity, pressure, density, temperature, or any other field parameters requiring interpolation. This type of interpolation eliminates the need for a fixed mesh but requires nearest-neighbor searches to discover interpolating particles at each material point [26]. A single-point Gauss quadrature approach integrates the fluid dynamics equations of motion. The particles are transported using simple implicit or explicit time-marching algorithms after computing the accelerations. The pressure field for fluid flow is derived by solving a weakly compressible state equation [27], or a pressure Poisson equation [28]. Because the SPH approach is meshless, it is easier to describe the violent fluid deformations exhibiting free surface behavior [29]. Furthermore, the Lagrangian frame eliminates the spatial discretization of the convection factor in the governing equations, thereby eliminating diffusion. Several open-source SPH software packages have been developed, including the SPHysics [30], which is currently being developed and Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 4 has been effectively used for environmental, oceanic, and coastal engineering challenges [31, 32]. However, most SPH models developed in the literature fail to consider the air- solid problem, especially the unsteady air-solid interaction of flapping wing propulsion [33]. The SPH method has been developed as a possible tool in aeronautical applications. Additional publications [34, 35] that connect SPH with the structural model opened up new avenues for investigating structure–air interactions considering structure deformation during an aerodynamic problem. This feature can solve the aerodynamic characteristics of flapping wings. In this study, we develop and validate an efficient computational model with a coupled air-solid interaction model for a 3D flexible clapping-wing. The structural domain representing the clapping-wing is modeled using rigid-body dynamics and flexible links, while the fluid domain is modeled using SPH. The flapping frequencies (F) adopted for the present study varied from 4 Hz to 20 Hz, based on the mechanical capacity of the drive mechanism and the theoretical frequency [8]. The incoming flow velocities (U) are based on three wind intensities - light or quasi- hovering (0.5 m /s - 1.5 m /s), moderate (2.5 m /s - 3.5 m /s), and heavy (4.5 m /s) winds. Furthermore, to evaluate unsteady mechanisms under different flight attitudes, -15 to 45◦angles of attack (AoAs) were adopted. The novelty of the present paper is threefold: (1) a validated and efficient coupled air-solid interaction computational model is proposed for incoming airflow impact on 3D flexible clapping-wing; (2) a parametric analysis is presented, where input variables represent different flight conditions that a clapping-wing is expected to face during its flight task; and (3) an understanding of the high lift mechanisms behind flapping-wing propulsion to derive new design requirements for CWMAV. The remaining paper is organized as follows. Section 2 describes the background of the proposed FSI model, experimental setup, and modeling scheme. Section 3 presents the results and discussions. Lastly, Section 4 summarizes the conclusions, limitations of the current work, and recommendations for future studies. 2. Method 2.1. Experiment 2.1.1. Test model Given the difficulties in accurately controlling CWMAV, the current investigation employs a tethered arrangement for our experiment. Qualitative smoke visualization is effective in exposing specific flow patterns [36]. A unique PIV flow field measurement system for CWMAV is proposed, consisting of a test model, clamping regulator, wind tunnel, mechanical data gathering system, and PIV system (figure 1). Synchronizers ensure that the signals from all measuring equipment have consistent time sampling. The test model comprises a biplane wing configuration with wings made of PE film reinforced with carbon rods as veins. According to the flight time sequence, the entire Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 5 Figure 1: PIV system of the CWMAV. flapping duration is divided into six phases filmed by a high-speed camera in Figure 2. The wingspan is 280 mm, and the weight is 17 g. The kinematic profile of the wing is shown in Figure 3. The instroke and outstroke phases indicate the folding and moving apart of the clapping wings, respectively. A clamping regulator placed at the bottom of the CWMAV secures the MAV, regulates the AoA and location, and provides input for control signals. 2.1.2. Experimental setup The experiments are conducted in a low-speed wind tunnel comprising a closed-circuit facility. The rectangular cross-section is 1.2 m wide, 1.2 m high, and 3 m long, as depicted in Figure 4. The entrance features a regular hexagonal honeycomb structure and a five-layer 24 mesh /inch damping net. Eight baffles are positioned at each of the wind tunnel’s four corners. The eight-blade fan in the power section has a diameter of 0.8 m and is operated by a variable frequency speed motor rated at 986 RPM. The tunnel contains an entering section with screens, and a 5:1 contraction ratio to ensure that the wind deflection angles are less than 0.5◦and turbulence levels are less than 0.5 %. The tunnel’s flow rate is monitored using a Pitot tube, and the free stream speed is regulated using a manual speed controller that drives the tunnel’s fan, generating an experimental flow rate ranging from 0.1 to 5 m /s, also recorded using a Pitot tube. The mechanical data acquisition system comprises a tiny six-dimensional mechanical sensor (FC-Nano17, ATI Industrial Automation, Fx, Fy= 8 N, Fz= 14 N, Mx, My, Mz= 0.05 Nm, 10.1 g, a sampling rate of 7200 Hz, a resolution of 1/160 N, and a measurement uncertainty of 1 % of full-scale load at 95 % confidence level), a data acquisition card (National Instruments, NI) and a terminal mechanical data acquisition software. This acquisition system offers high precision in a small range, suitable for the aerodynamic force acquisition of CWMAV. Additionally, a PIV image capture system consisting of a complementary Metal- Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 6 I IV II V III VI Figure 2: Visualization of the flapping cycle represented by high-speed camera recordings. oxide-semiconductor (CMOS) high-speed camera (FASTCAM Mini AX, Photron, 6400 Hz, 1024 ×1024 pixels, ISO 40000, Figure 5(a)), full E-mount lens (Grand Master, Sony, 85 mm, F1.4), middle energy level Q-switched nanosecond Nd:YAG laser (Nimma series, eamtech, 532 nm, 250 mJ, 1-10 Hz, 2560 ×2160 pixels, Figure 5(b)) is used. PIV’s synchronization controller coordinates and controls the CCD high-speed camera, high- frequency laser, and image acquisition software. The PIV and CWMAV synchronization controllers are operated independently via a personal computer. The computer stores and analyzes all images collected by a CCD high-speed camera. As shown in Figure 5(c), our test model is placed in the wind tunnel. The velocity and vorticity fields are recorded at five distinct spanwise locations of the imaging plane for multiple flapping cycles at different flapping frequencies. Figure 5(d) illustrates the operation of the PIV image acquisition module. These experimental setups can measure the flow field surrounding the CWMAV and the aerodynamic force/torque generated by the CWMAV during the flapping period. 2.2. Proposed modeling technique A brief background of the smooth particle hydrodynamics (SPH) method and details of the contact algorithm used in this study are presented below. Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 7 Figure 3: Flapping kinematics. 2.2.1. SPH implementation All SPH simulations are performed using the open- source code DualSPHysics [37, 38], with the air and solid domains defined as discrete particles. The SPH equations are solved in their weakly compressible form. The discrete approximation for a physical quantity, β, for particle iin DualSPHysics, is provided by βi=X j∈ΩβjWi,jVj, (1) where j∈Ω denotes the set of neighboring particles. The kernel function is represented Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 8 Figure 4: Low-speed wind tunnel. byWi,j=W(|xi,j|, h) and is calculated as a function of the particle distance ( |xi,j|= |xi−xj|) and smoothing length, h. The volume of a neighboring particle jis given by Vj=mj/ρj, where mjandρjsignify the mass and density of particle j, respectively. For all simulations, a quintic Wendland kernel [39] is used, defined as: Wi,j=αD(1−2q)4(2q+ 1)for0≤q≤2, (2) where q=|xi,j|/h, and αDis a normalisation term. For the 3D simulations presented in this paper, αD= 21/(16πh3), and his set to 1 .2√ 2dp, where dpis the particle spacing. 2.2.2. Governing equations Fluid quantities are calculated using the principles of mass (continuity) and momentum conservation, as defined by Dρ/Dt +ρ∇ ·u= 0, (3) Du/Dt=−1 ρ∇p+g+ Γ, (4) where ρrepresents the fluid density, u= (u, v, w ) is the velocity vector with components in the ( x, y, z )-directions, pis the fluid pressure, and gis the gravitational acceleration. Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 9 Figure 5: PIV instruments. (a) High-speed camera. (b) High-frequency laser generator. (c) Positioned test model. (d) Diagram of laser irradiation. The material derivative is denoted by D/Dt , while the dissipative terms are denoted by Γ. The continuity equation’s weakly compressible SPH form, incorporating the δ-SPH density diffusion term of Fourtakas et al. [40], is given by dρi/dt=X j∈Ωmjvi,j· ∇Wi,j+δhc0X j∈ΩVjΨi,j· ∇Wi,j, (5) where vi,j=vi−vjand∇Wi,jis the kernel gradient. The speed of sound c0is set to 20√ghd, where√ghdis the phase speed in fluid at a depth of ghd. Gravitational acceleration is indicated by g. The δ-SPH coefficient, δ, is assumed to be 0.1 (e.g. Kanehira et al. [41]). The diffusion term Ψi,jis given as in Molteni and Colagrossi [42]: Ψi,j= 2(ρD j−ρD i)xi,j |xi,j|= 2(ρT i,j−ρH i,j)xi,j |xi,j|, (6) which is the formulation first described in Fourtakas et al. [40]. The superscripts D, T, and Hsignify the dynamic, total, and hydrostatic densities, respectively. The Tait equation of state is used to relate the pressure and density, and thus the conservation Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 10 of mass and momentum equations, for weakly compressible SPH: p=c2 0ρ0 γ[(ρ ρ0)γ−1], (7) where ρ0= 1.29kg/m3is the reference density of air, and γ= 7 is the polytropic index. Eq. 7 is an extremely stiff equation, which when combined with particle disorder, produces large non-physical density fluctuations. In the present simulations, the aforementioned δ-SPH density diffusion factor is incorporated in Eq. 5 to dampen these density fluctuations. The momentum equation, laminar viscosity, and sub-particle scale (SPS) turbulence are as follows: (µ∇2v)i=X j∈Ωmj(4µ|xi,j| · ∇Wi,j (ρi+ρj)(x2 i,j+ 0.001h2)), (8) where µis the kinematic viscosity, which is set to 0.0000148 based on the validation studies discussed next. This treatment for the turbulence phenomenon can avoid numerical instabilities without leading to an over-diffusive behavior [43]. The shifting algorithm used in SPH is highly dependent on the presence of a complete kernel. Additionally, time-stepping is accomplished with a symplectic explicit second-order time-integration method based on prediction and correction. The dimensionless Courant number that represents the time that an SPH particle stays in the kernel domain is set to 0.2. 2.2.3. Contact algorithm for fluid structure interaction (Coupling DualSPHysics with Project Chrono) The multi-body solver used in this study is from ”Project Chrono” [44, 45], an open-source high-performance library and collection of programs that supports various applications. DualSPHysics 5.0 uses the Project Chrono library’s rigid body, constraint, and collision detection components. The current integration of DualSPHysics and Project Chrono enables using multiple rigid bodies with the possible constraints applied to each rigid body specification of restrictions on each body (from an implemented list of joints, hinges, and springs). The computation of interactions between bodies is more stable than the discrete element method (DEM) implementation. Additionally, the Project Chrono framework enables the modeling of structure–structure interaction with different contact tracing capabilities. The library is designed to handle massive systems of 3D rigid bodies. The coupling enables considering bodies of any form, and the solver incorporates externally imposed forces and torques, and the impacts of kinematic and dynamic-type constraints, and interior collisions. The implementation for the CWMAV simulation is summarized in Figure 6: DualSPHysics calculates the linear and angular accelerations of floating objects based on their interaction with fluid particles, and passes the informationdV dt,dΘ dt; ∆tSPHto the Project Chrono rigid- body dynamics module, which integrates the equations of motion in time using smaller time steps ∆ tC. The active and reactive forces exerted on bodies are calculated using the implemented list DualSPHysics using the position R, angular velocity dΘ, and linear velocity dVdetermined by Chrono to update the data for the floating particles Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 11 associated with each floating object. Therefore, the implementation technique adopted in our study is to connect the partitioning via a message-passing interface. Figure 6: Contact algorithm for fluid-structure interaction (FSI). 2.3. Parametric study A parametric study is considered to analyze the effect of (1) flapping frequency, (2) angle of attack, and (3) wind velocity on the responses of the aerodynamic forces. The following subsections describe the modeling parameters considered for the parametric analysis. 2.3.1. Computational flapping wing modeling A summary of the important parameters used in the simulations is presented in Table 1. The numerical flapping wing is a 3D representation of the actual CWMAV, whereas the flow field is a cuboid with dimensions of 0.8 m ×1 m×0.2 m. Flexibility defining has been challenging in SPH Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 12 frameworks [46], and the following strategy is adopted to overcome it. The flapping wing under examination is constructed from carbon rods and PE film, as represented in Figure 7(b), and is based on the experimental wings (figure 7(a)). The framework that connects DualSPHysics and Project Chrono can only enable fluid-structure interaction between a fluid and a rigid-body solid portion, indicating that no flexible solids are supported. Thus, a quasi-flexible wing is proposed by substituting spring connections for the obliquely placed carbon rods. To achieve a material property comparable to that of a PE film, the Young module of the wings is increased to 0.2 GPa, with a Poisson ratio of 0.3, a restitution coefficient of 0.001, and a coefficient of friction of 0.45. The link configuration is shown in Figure 7(b). We segmented the wings along the dotted line and divided them into six partitionings. The wing partitioning is fixed at the intersection between the wing root and leading edge (shown in the red circle in Figure 7(b)) using ball links and at the other end (shown in the blue circle) using linear springs. Flexibility can be controlled by changing these link parameters. The simulated flexibility of the flapping wing is shown in Figures 7(c) and (d). The developed model is available in a public data set [47]. The FSI framework in this paper is as follows: As the clapping wing moves, the solid domain receives the physical quantity by computational fluid dynamics (CFD) in the fluid domain and transmits it to the links in the six partitionings. The solid domain then completes the deformation by stretching the linear springs. Finally, the fluid domain is recalculated when the solid domain moves and deforms during the flapping flight. The stiffness and damping of the linear springs are set to 300 and 0.1, respectively, so that the numerical flexible deformation is close to the observed flexible deformation in the experimental images. The ball links are used to ensure zero displacement of the contact points. The motivation for this design is that the veins play an important role in flexible deformation, and these links are set to play a similar role. The flexibility of the wing is innovatively simulated through links. It can also be applied to other wing layouts by different segment methods and stiffness or damping values of the links. Hence, the proposed design solves our demand for building flexible wings and expands the application of the SPH method in handling flexible bodies. 2.3.2. Input parameters The flapping frequency used for the analysis is an important parameter as different flapping frequency imparts different flight patterns, wingtip velocity, and required time steps for our simulation. The flapping frequency is chosen based on different flight conditions to support CWMAV in forward and hovering flight. The flapping frequency distribution used in this study is based on the model developed by Per¸ cin et al. [8]. The adopted flapping frequencies used in the simulations and validation are 4, 8, 12, 16, and 20 Hz. The angles of attack considered in the present study are defined as the relative angle between the horizontal plane and CWMAV in Figure 8. We selected -15, 0, 7.5, 15, 30, and 45◦as the angles of attack (AoAs). These differences represent different flapping motion attitudes for our test model, i.e., AoA <0, AoA = 0, and AoA >0 Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 13 Table 1: Key parameters and formulations used in the simulations. Parameter Value SPH kernel Quintic Wendland Fluid and solid interaction Chorno Particle spacing 0.004 m Smoothing length 1 .2√ 2dp Density diffusion 0.1 Diffusion term Fourtakas et al. [40] Particle shifting Off Speed of sound 20√ghd Reference density 1.34 kg/m-3 Polytropic index 7 viscosity 0.0000148 Time integration Predictor–corrector Courant number 0.2 Simulation time 1 s Simulation output frequency 1000 Hz Dynamic boundary condition mDBC [48] represent descendent, forwarding, and ascendant flight tendencies, respectively. Considering the flight velocity and wind conditions, three classes of wind velocity are defined: light wind or quasi-hovering (0.5 m /s - 1.5 m /s), moderate wind (2.5 m /s - 3.5 m /s), and heavy wind (4.5 m /s). 2.3.3. Validation study Validation research is conducted to ascertain the computa- tional models’ convergence conditions and accuracy. For the validation cases, the start- ing particle spacing dpis adjusted between 0.008 m, 0.004 m, and 0.002 m. Figures 9 illustrates the disparities between experimental and SPH measurements by lift evolution of a periodic stroke. The raw force data is filtered with a low-pass filter having a cut-off frequency twice the flapping frequency to consider aerodynamic force fluctuations alone. The measurements are divided into two groups to reduce the influence of inertial forces: wings without PE film (only with flapping driving mechanism and carbon rod) and the complete wings. The experimental results are based on the difference between these two sets of values. The numerically predicted trends are similar to the experimentally observed trends: At the outstroke phase, there is a sharp rise at first, reaching the first lift peak. After the first peak, there is a high fluctuation, and the lift remains at a high value. At the end of the outstroke phase, the lift rapidly declines, and the value approaches 0; A similar trend is observed at the instroke phase. The total number of particles is 108987, 1576003, and 12184977 for dp= 0.008 Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 14 Figure 7: Flapping-wing structure. (a) Experimental wings, (b) Computational wings, (c) Rear view of the flexible wing, d) Front view of the flexible wing. m, 0.004 m, and 0.002 m, respectively. The comparison of lift history and flapping- cycle-averaged lift demonstrates that situations with dp= 0.004 m and dp= 0.002 m have similar trends and values as the experimental data. The main reason for the differences in the results is that the movement and deformation range of wings are extensive, and a smaller dpvalue can capture more unsteady characteristics to make lift evolution more similar to experimental measurement. Comparative analysis shows that when the total number of particles is included, setting dp= 0.002 m results in an excessively time-consuming model that takes 20 h every 1 s of simulated flapping time at a computational capability of 8.6 (with NVIDIA GeForce RTX 3060Ti). Thus, to strike a balance between sufficient precision and computational expense, dp= 0.004 m is used in all subsequent simulations. The resulting average lift discrepancies between experimental and SPH measurements are shown in Table 2. The average lift is derived during the second period of flapping flight (from instroke to downstroke). There is a deviation from the experimental result, but the difference in most cases is less than 10 %, which is within our acceptable range. Despite the considerable gap between our numerical and experimental peak values in the instantaneous force curve, the high similarity of our instantaneous forces and the low difference in the average forces prove that setting dp= 0.004 m is sufficient to meet the accuracy requirements. Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 15 Figure 8: Description of the angle of attack. 9 Figure 9: Difference between experimental and SPH FSI measurements. (F = 16 Hz, U = 3.5 m /s, AoA = 15 deg .) Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 16 Table 2: Average lift differences between experimental and SPH FSI measurements. Case Experimental measurements (N) SPH FSI measurements (N) Differences (%) F = 8 Hz, U= 3.5 m/s AoA = 15 deg. 0.1084 0.1352 24.7 F = 12 Hz, U= 3.5 m/s AoA = 15 deg. 0.1681 0.1603 4.6 F = 16 Hz, U= 3.5 m/s AoA = 15 deg. 0.1562 0.1661 6.3 F = 16 Hz, U= 0.5 m/s AoA = 15 deg. 0.0414 0.0381 7.9 F = 16 Hz, U= 1.5 m/s AoA = 15 deg. 0.1183 0.1048 11.4 F = 16 Hz, U= 2.5 m/s AoA = 15 deg. 0.1089 0.0949 12.8 F = 16 Hz, U= 4.5 m/s AoA = 15 deg. 0.2342 0.2308 1.4 F = 16 Hz, U= 3.5 m/s AoA = -10 deg. 0.0202 0.0221 9.4 F = 16 Hz, U= 3.5 m/s AoA = 0 deg. 0.0401 0.0412 2.7 F = 16 Hz, U= 3.5 m/s AoA = 7.5 deg. 0.0714 0.0842 17.9 F = 16 Hz, U= 3.5 m/s AoA = 30 deg. 0.3039 0.3040 0.0 F = 16 Hz, U= 3.5 m/s AoA = 45 deg. 0.3195 0.3326 4.1 3. Results and discussions 3.1. Parametric study: Effects of flapping frequency, angle of attack, and wind velocity To assess the relationship between the lift and flapping motion , lift evolution calculated by the proposed SPH FSI model in Figure 9, experimental images in Figure 2, and experimental vorticity contour in Figure 10 are combined. As the peeling occurs at the start of the outstroke, a substantial increase in lift is established, with a peak occurring at approximately t∗= 0.3 (phase II). The high-speed image at this instant shows that although the wing leading edges are apart, driven by the flapping mechanism at a stroke angle of approximately 30 deg ., the trailing edges get ready to be separated. At this moment, the LEV grows larger due to the low- pressure region formed by leading-edge separation. A gradual decrease in a lift during the second half of the outstroke ends in a negative lift peak when two wings are at around t∗= 0.4. Correspondingly, LEV moves toward the trailing edge along chord-wise wings and disappears at the end of the instroke phase. Thus, at the start of the instroke, particularly at t∗= 0.6, the lift is at the minimum of the complete flapping period. In phase IV, which corresponds to a stroke reversal, the upper and lower wing leading edges are almost stationary, and the wing surfaces have a small amount of deformation. As the upper and lower wings start moving in the instroke direction, the force gradually increases with increasing wing speed and deformation, and the new LEV and TEV appear. The maximum force during the instroke, a second lift peak, is attained at about t∗= 0.8, while the second lift peak value is smaller than the first peak. The size and strength of LEV and TEV in this phase are enormous as observed in the vorticity contour. For instance, the low-pressure region formed by LEV requires a vertical wing surface area to produce a suction force. In the ending period of the instroke, the stroke velocity decreases, reducing the magnitude of the resultant forces. Subsequently, the third minimum of the lift is reached (at t∗= 1, phase VI) at the end of the instroke. In this timing from the flapping images, the wings clap, and presumably Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 17 Figure 10: Vorticity contour of six phases. (AoA = 15 deg .U = 3.5 m /s, and F = 16 Hz) the momentum jet generated due to this ipsilateral wing-wing interaction contributes to the lift (vertical force) production because the FWAMV is under a certain value of angle of attack. The vorticity contour reveals that the previous vortices disappear because the leading edges touch each other, and a new LEV with weak strength on the outside of the wings occurs. The relationships between flapping motion and lift history illustrate the importance of wing flexibility and promote our further parametric study of unsteady mechanisms that do not comply with the laws of conventional aerodynamics. The influence of flapping frequency ( F = 4, 8, 12, 16, and 20 Hz) on force generation is discussed in Figure 11(a) via the lift history calculation where the wings flap with AoA = 15 deg .and U = 3.5 m /s. An increased amplification in lift generation during instroke and outstroke is observed as a direct consequence of the higher wing velocities. In contrast, the larger flapping frequencies could change the relative size relation between the first and second lift peaks. Therefore, increasing the flapping frequency could increase the average lift, but not every instant. Moreover, increasing the wings’ flapping velocity increases the dynamic pressure, proportionally increasing the lift while changing the lift history. Figure 11(b) compares the impact force history for cases with varying angles of Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 18 attack for wings flap (AoA = -15, 0, 7.5, 15, 30, and 45 deg .) at F = 16 Hz and U = 3.5 m/s. In addition to the flapping velocity, AoA plays an important role in influencing lift evolution. Among the tested configurations, AoA ≤30 deg .shows similar trends in lift; the lift increases rapidly at all time points as AoA = increases. When AoA ≤ 7.5 deg ., the lift provided is insufficient to support the flight. Moreover, there is an obviously negative lift when AoA = -15 deg .. Furthermore, a higher AoA leads to a larger amplitude, which inevitably destroys CWMAV stability. To investigate the effects of wind velocity (U = 0.5, 1.5, 2.5, 3.5, and 4.5 m /s) on the aerodynamic behavior, a comparative analysis of different wind velocities (U = 0.5, 1.5, 2.5, 3.5, and 4.5 m /s) is performed fixing AoA = 15 deg .and F = 16 Hz, with the results in Figure 11(c). The relative difference in the performance measures of these different cases is in the predicted amplitudes, while the force shape and trend are well captured. Moreover, the force coefficients follow a slightly increasing trend concerning wind velocity. In addition, a comparison of the time-average and maximum lifts for varying flapping frequencies and angles of attack are presented in Figures 12(a) and (b), respectively. In Figure 12(a), F = 20 Hz has a relatively high lift for the complete range of test flapping frequencies. The maximum force is generated by F = 20 Hz at the expense of the maximum flapping frequency. Combined with Figure 12(b), AoA is identified as another key parameter influencing average and maximum lifts. The wings with low flapping frequency have poor characteristics under our flapping conditions. The average and maximum lifts are low when AoA = 15 deg .; hence, the cases of AoA = 15 deg.has no value for further research and can only be used for parameter comparison. In addition, for cases of AoA ≥15 deg ., F≥12 Hz outperforms the lift performance, which is sufficient to support the weight of our test model with a typical mass of 17 g; however, (0.17 N) in hovering configuration produces increased negative lift. Moreover, the average lift of AoA = 30 deg .is two times larger than that of AoA = 15 deg ., and the maximum lift of AoA = 15 deg .is two thirds more than that of AoA = 30 deg .. 3.2. Parametric study: Reduced frequency and Reynolds number To extend our parametric study to other CWMAVs of the same type but different sizes or environments, the effects of varying reduced frequency ( λ) are briefly discussed. ( λ) is an important similarity parameter defined as the measure of unsteadiness in flapping- wing flight. In forward flight, the free-stream velocity (U) and angular frequency of the flapping motion (2 πF) based on the mean chord length (c = 0.082 m) is defined as follows: λ=πFc U, (9) One period of the flapping motion is given in Figure 13. We investigate the effect of changing λbetween 0.229 (F = 4 Hz and U = 4.5 m /s) and 1.648 (F = 20 Hz and Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 19 Figure 11: Comparison of lift history for flapping wings: (a) AoA = 15 deg ., U = 3.5 m/s, F = 4, 6, 8, 12, 16, 20 Hz (b) AoA = -15, 0, 7.5, 15, 30, 45 deg ., U = 3.5 m /s and F = 16 Hz (c) AoA = 15 deg ., U = 0.5, 1.5, 2.5, 3.5, 4.5 m /s, and F = 16 Hz. Figure 12: Comparing varying flapping frequencies and angles of attack of (a) time- average and (b) maximum lifts Bionic simulation of double clap-and-fling wing mechanism with SPH FSI method 20 U = 2.5 m /s) where AoA is locked at 15 deg .. The largest value of lift is λ= 0.687 ((corresponding to F = 12 Hz and U = 4.5 m /s). At λ≤0.687, the average lift is slightly smaller than 0 .687≤λ≤0.178 (corresponding to F = 12 Hz, U = 4.5 m /s and F = 16 Hz, U = 3.5 m /s, respectively). However, as λcontinues to increase, the average lift reaches a lower trough after λ≥1.178. The maximum lift also appears at λ= 0.687, which is a 60 % increase over the minimum value of the maximum lift. Moreover, there is a strong correlation between the maximum and average lift, presumably due to the mechanical characteristics of CWMAVs. Figure 13: Lift performance at different reduced frequencies: (a) lift history (b) average and maximum lifts 4. Conclusion and next steps The present study focuses on developing an air-solid interaction model to reduce the costs and time involved in designing CWMAV. The numerical model merges SPH to represent the air and rigid-body dynamics to describe a bionic wing. The proposed model is compared and validated using the experimental motion and vorticity contours. A numerical and experimental model investigates the aerodynamic characteristics of the clapping wings to understand the unsteady mechanism of CWMAV better. This FSI model is the theoretical basis of the control model, which is convenient for determining the mapping relationship between the wind speed, flapping wing frequency, angle of attack, and aerodynamic force for controlling. The vortices, flapping motions, and instantaneous lifts of each phase are obtained and compared to determine the high lift mechanism. The effects of wind speed, flapping frequency, and angle of attack on aerodynamic performance are explored by a parametric study, focusing on lift evolution, time-average lift, and maximal lift in a cycle. The results show the significant impact of frequency on the lift. The lift curve changes dramatically when the frequency increases, with the maximum and average values growing monotonically. When the angle of attack is less than 30◦., the lift evolution essentially is constant with further increase in the angle of attack, while the influence on the average and maximum values are proportional to the angle of attack. However, a