Kinematic Optimization of a Flexible Wing Undergoing Flapping and Pitching

Wu, Yue; Xie, Changchuan; Meng, Yang; Yang, Chao

doi:https://doi.org/10.1155/2021/9635584

Shock and Vibration

On this page

Abstract Introduction Results Conclusions Data Availability Conflicts of Interest References Copyright Related Articles

Research Article | Open Access

Volume 2021 | Article ID 9635584 | https://doi.org/10.1155/2021/9635584

Kinematic Optimization of a Flexible Wing Undergoing Flapping and Pitching

Yue Wu,¹Changchuan Xie,¹Yang Meng,¹and Chao Yang¹

Academic Editor: Francisco Beltran-Carbajal

Received24 Sept 2021

Revised09 Nov 2021

Accepted12 Nov 2021

Published15 Dec 2021

Abstract

In recent years, there has been widespread interest in the design of microair vehicles (MAVs) for flapping flight with high-aspect ratio wings due to their high efficiency and energy savings. However, the flexibility of a flapping wing causes the aeroelastic effect, which remains a subject of investigation. Generally, existing research simulates active bending and twisting of flexible wings under the assumption of neglecting flapping inertia. In this research, the kinematic optimization of a bionic wing with passive deformation in forward flight while undergoing flapping and pitching is considered. To this end, a computational aeroelasticity framework, which includes the three-dimensional unsteady vortex lattice method (UVLM) and the Newmark-β method, is constructed for flapping flight. Under the assumption of linear elastic deformation, this tool is capable of simulating attached flows over a thin wing and capturing unsteady effects of wakes. A bionic numerical wing with an aspect ratio of 6.5, chord Reynolds number of 1.9 × 10⁵, and reduced frequency less than 0.1 is investigated in kinematic optimization. The computational aeroelasticity framework is combined with a global optimization algorithm to identify the optimal kinematics that maximize the propulsive efficiency under the minimum average lift constraint. Two types of numerical wings, rigid wing and flexible wing, are considered here to compare the influence of deformation on the aerodynamics of the flapping wing. The results show that the aeroelastic effect, which increases the flapping amplitude, yields a significant improvement in terms of propulsive efficiency. In addition, the optimization algorithm maximizes the thrust efficiency while satisfying the required lift. Moreover, the optimal kinematics of both the rigid wing and the flexible wing reach the maximum flapping angle, which indicates that a larger range of motions is needed for optimal kinetics when loosening the boundary conditions.

1. Introduction

Recently, various ornithopters have been designed and manufactured worldwide. Such aircrafts gain lift and thrust by flapping their wings. Their weights are limited to the magnitude that flapping wings can drive, which leads to the application of flexible lightweight materials in the wing components. Based on experimental [1, 2] and numerical results [3–5], flexibility plays an important role in the aerodynamic performance of flapping wings. However, the motion mechanism of the ornithopter determines that it has a fixed kinematic mode after it is manufactured. However, there is a lack of methods to optimize the kinematics of flexible flapping wings in the design stage. To determine the coupled motion of flapping and pitching of a microflapping wing machine, Matthew et al. [6] selected optimal kinematics from dozens of motion parameter combinations through experiments. This simple way to design kinematics is an alternative method in the absence of proper optimization tools. Other studies directly used common functions to represent flapping wing motion. Nick et al. [7] used symmetrical in-plane motion on a micro bat-like ornithopter without applying a kinematic optimization method to improve its aerodynamic characteristics. Razak and Dimitriadis [8] investigated the coupling effect of flapping and pitching motions on the aerodynamic characteristics of flapping wings. Their experiment explored the aerodynamic effect caused by the phase difference between flapping and pitching motions but without considering the influence of other motion parameters. Sutthiphong and Chan [9] applied a simple sawtooth function to the variation process of the flapping angle with time when exploring the flow field and noise characteristics of large bird wings due to the lack of data on the motion process. At present, the physical phenomena caused by fluid-structure interactions are not completely understood on a flexible wing. An abstract model with simplified wing structure and kinematics is an effective approach to overcome this challenge. Thus, a fully passive dynamical model without active deformation is applied in this research. The current work combines a computational aeroelasticity framework with a global optimization algorithm to construct a design method for flapping wing kinematics. The computational aeroelasticity framework includes the three-dimensional unsteady vortex lattice method UVLM and the Newmark-β method to calculate the fluid-structure coupling effect caused by flexible deformation.

The flapping wing kinematics have a complex design space with multiple local optima when taking specific aerodynamic characteristics as the optimization goal. A widely used feature to evaluate the aerodynamic characteristics in the UVLM is the propulsion efficiency. It is defined as the ratio of the propulsion power to the total aerodynamic power [10, 11], which can be easily obtained by the UVLM. When the flapping wings create a specific degree of disturbance to the airflow by motions and deformations, the maximum propulsion efficiency can be achieved, which indicates the existence of a set of optimal kinematics to maximize the propulsion efficiency [1]. This article uses the DIRECT (dividing rectangles) optimization algorithm to move towards the maximum-efficiency goal. Motivated by a modification to Lipschitzian optimization, the DIRECT optimization algorithm was first introduced by Perttunen [12]. It attempts to solve the global optimization problems of multiple parameters with bound constraints and a real-valued objective function. However, global searching leads to a high computational cost to ensure convergence to the optimal solution. As a sampling algorithm, DIRECT, which requires no knowledge of the gradient of the objective function, is suitable for searching for the optimal value of a “black box function” or simulation result.

Flapping motion has a low frequency, which means that the calculation time lasts for several seconds in several flapping wing cycles to achieve a stable state. In this paper, the global optimization method needs to sample large numbers of calculation states, and the method of medium fidelity, the UVLM, is used in aerodynamic calculations. The UVLM has been widely used in flight dynamics and aeroelastic analysis for aerodynamic modelling [13, 14]. Mehdi et al. [15, 16] applied the UVLM to estimate the aerodynamic force in flapping wing research. Even in flow at low Reynolds numbers [17] and separation flow [18–20], the UVLM maintains a reasonable trend correlation with CFD (computational fluid dynamics) results.

The objective of the current paper is to investigate the optimal kinematics of a specific flexible wing undergoing flapping and pitching at the wing root. First, a flexible high-aspect ratio wing model is constructed according to the data of natural avians. Then, a computational aeroelasticity framework is formulated to predict the aerodynamic effect and propulsion efficiency of the preset flexible wing. After that, the global optimization problem is presented for the DIRECT algorithm. To prevent the effect of the flow separation caused by the flapping motions, the flapping range is selected according to the relevant experimental results. Finally, the kinematic and aerodynamic results of the flexible model are discussed and compared with those of a rigid model.

2. Wing Kinematics and Dynamics

In this section, the kinematics of a bionic slender wing flapping and pitching at the wing root are presented. The material properties and wing size, including the camber of the airfoil, are designed according to the results from research on a kind of high-aspect ratio seabird, the black-browed albatross. Although the exact kinematics of the wings of albatrosses in nature remain unknown, the characteristic angles of flapping and pitching in this article are expressed by sine functions. The Lagrange equation under rigid-elastic coupled motion is established here to solve the computational structural dynamics for the simulation of passive wing deformation.

2.1. Wing Model

Natural avian flight with a high-aspect ratio wing is commonly seen on albatrosses that migrate thousands of miles over the seas. Although albatrosses apply unpowered flight modes, such as gliding and soaring, to save energy, flapping flight by this species has been observed under calm conditions in wind fields [21]. When the speed of sea wind stays too low to supply energy for unpowered flight, flapping flight is required to maintain flight altitude. Based on the available morphological data in the literature, a kind of high-aspect ratio seabird, the black-browed albatross, is considered. A numerical wing model is designed to investigate the flexibility based on the morphological data, as shown in Table 1. The reduced frequency and Reynolds number are also calculated according to the data in Table 1. Notably, the density of the wing is approximately equal to the density of insect wings [23]. The same treatment of the estimation of the wing density was also adopted by Kodali [24].

Further consideration of the numerical wing model focuses on the wing planform and airfoil camber. The real wing planform abstracted from natural observation [21], including the shoulder spacing, is plotted in Figure 1. The fact that albatrosses commonly have a sharp shape at the wingtip is equivalent to having a high taper ratio. To simplify the wing model into linear leading and linear trailing edges, the sweep angle is set to 0^°, and a taper ratio of 7.244 is chosen for the numerical wing. Combined with the half span length and wing area data in Table 1, the exact planform of the numerical wing can be obtained, as shown in Figure 1. Because a wing with a thin airfoil needs to be replaced by a curved surface in the UVLM, the middle camber of the airfoil is chosen to model the wing. In Figure 2, schematic views of the dimensionless albatross airfoil (GOE 174) [27] and the camber are shown. Due to the lack of weight data for a single wing of the black-browed albatross, the mean thickness h_s = 8.52 mm is obtained according to a GOE 174 airfoil with the length of a mean chord shown in Table 1 to construct a wing with numerically equal thickness. Although the numerical wing is overweight compared with a real bird wing due to the sparseness of feathers, the redundancy of the material enhances the stiffness of the numerical wing, which linearizes the deformation caused by motions. This simplified wing model assumes that the wing is composed of isotropic materials, which limits the aerodynamic mechanisms of feathers, such as interlocking. Moreover, the taper ratio and camber have a positive effect on the stiffness of the model. The effects of gravity are considered in the results. In addition, the weight of the flyer in Table 1 is used only as a lift constraint in kinematic optimization.

2.2. Wing Kinematics

This article considers a three-dimensional wing in constant wind speed U_∞ that undergoes flapping and pitching at the wing root. As the shoulder joint of the bird is located at the root of the leading edge of the wing, the kinematic centre of rotation is placed at the root of the leading edge. In addition, the y-axis and x-axis of the local coordinate system are defined as the spanwise and chordwise directions of the wing, respectively. A schematic of the present case is shown in Figure 3. To describe the wing location, the flapping angle γ and the pitching angle θ are defined as functions of time t as follows:

The motions with constant frequency f contain the adjustable kinematic parameters of flapping amplitude γ_a, middle pitching angle θ₀, pitching amplitude θ_a, and pitch-leading phase φ_θ.

In Figure 3, the inertial coordinate system O-X-Y-Z is represented by E and the local coordinate system o-x-y-z is represented by B. As shown in Figure 3, the transformation from E to B begins with rotation γ around the X-axis first and then proceeds with rotation θ around the y-axis. Since rotation about z is relatively rare on existing ornithopters, kinematic optimization including this motion seems uncommon for aircraft design. Thus, the rotation about z is always set to 0^°. Therefore, the vector describing the relationship of the relative angles between E and B is defined as (here, the bold font indicates a vector or matrix). The vector composed of the angle velocities of the wing described in the local coordinate is defined as . According to the projecting relationships, ω_B is expressed by the derivative of ϕ with time as follows:where the mark with one dot above the symbol represents the first derivative of time and D_B is the transition matrix as follows:

Thus, equation (3) has the following expansion form:

2.3. Structural Dynamics

The passive deformation caused by inertia and aerodynamics is solved in the local coordinate system using the Lagrange equation. To express the kinetic energy T of an arbitrary microscale unit of the wing, the nondeformation position vector and the displacement vector of a microscale unit are defined as and in the local coordinate system, respectively. If the mass of the microscale unit is dm, then the kinetic energy of the microscale unit is as follows:where the marks with a curve over the symbols represent the second-order tensors that are defined as follows:

Assuming that k is the stiffness matrix and U_e is the elastic potential energy of the microunit, we have

According to the Lagrange equation, the kinetic energy T, the elastic potential energy U_e, and the force vector f on the microscale unit are related via the following equations:which is then reduced to the simple form by substituting (6) and (9) into (10).where the marks with two dots above the symbols represent the second derivative with respect to time. If the wing model is discretized into n elements in the finite element method, the total elastodynamic equation of the finite element model is given as follows:

The symbols in the total elastodynamic equation are defined as follows:

Based on the modal method, the displacement U_B is the sum of the modes Φ weighted by the generalized coordinate q:

Substituting equation (17) into dynamic equilibrium equation (12) yields the following result:

The modal mass matrix is normalized by , and the modal stiffness matrix is defined as . The replacements of the terms in equation (18) are performed as follows:

Thus, equation (18) is arranged as follows:

As an implicit method, the Newmark-β method is often used for solving structural dynamics problems discretized by time steps. In the step n_t + 1, the dynamic equation is shown as follows with the generalized coordinate and its time derivatives:

By assuming the stepping relation between time steps, the Newmark-β method transforms the dynamic eqaution (22) into a simple form as follows:where

When the length of the time step is determined, the coefficients a_i,i ₌ _1,…,6 become a set of constants defined by the Newmark-β method. If the current state and the generalized force D are known at each time step, then equation (23) can be solved to yield the solution of the next state. Further instructions of the Newmark-β method are found in reference [28]. By substituting the following equations into equations (24) and (25), the generalized coordinate q in equation (25) is solved in the time-stepping method:

To introduce the modal method into the structural dynamics, the vibration modes and natural frequencies of the numerical wing in which the grid has 7 × 27 nodes are analysed by the finite element method in the software Nastran. Although each node has six degrees of freedom, three for translation and three for rotation, the structural model extracts three translational degrees of freedom for the dynamic calculation of elastic deformation. The mesh of the wing is shown in each plot of Figure 4, and the x and y axes are along with the wing chord and span direction. On the whole wing, each region surrounded by four corner nodes is modelled with CQUAD4 element, and the solver of SOL 103 is used to generate the structural modes Φ. The mass matrix M can be calculated by Nastran when adding “PARAM, EXTOUT, DMIGPCH” to the calculation card. Under the assumption of linear elastic deformation [29], the effect of flapping and pitching motions is equivalent to the inertial force and Coriolis force on the cantilever wing in the local coordinate system, in which the y-axis and x-axis are defined as the spanwise and chordwise directions, respectively [30, 31]. To this end, the clamped boundary condition at the wing root is used to constrain the deformation of the model in the local coordinate system. The analysed results of the first four modes are shown in Figure 4. The flapping frequency in Table 1 is lower than the frequency of the first mode, suggesting that the elastic deformation of the flapping wing is dominated by the first bending mode. Although it cannot be known for certain whether the natural frequencies and vibration mode shapes are similar to those of the real wing, the experimental [1] and simulation [24] data indicate that the propulsion efficiency is improved by using the first-order bending mode.

3. Unsteady Vortex Lattice Method

Based on inviscid flow theory, the UVLM is widely used to predict the flow field evolution and aerodynamic force of a thin wing. Under the condition of a certain air velocity, the three-dimensional UVLM discretizes the free wake that forms at the trailing edge into interconnected vortex elements, whose instantaneous velocities are governed by the Biot–Savart law [14]. In this work, 5 vortices in chordwise and 16 vortices in the spanwise are assigned as bound vortices. The wake is held to generate velocity induction for 20 chord lengths. It is important to note that the vortex line lies on a 1/4 chord length of each panel, and the panel control point lies on the middle point of a 1/2 chord length of each vortex ring as shown in Figure 5. The green lines represent the vortex lines, and the red dots represent the control points. The zoom view clearly shows the vortices and direction of vortex rings.

In the coordinate system shown in Figure 3, the induced velocity caused by the attached and free vortices is represented by the vector . In addition, each vortex in the wake must move with the free-stream velocity U_∞. Thus, the displacement of a free vortex in each time step Δt is determined via the following scheme:

In a flapping wing problem, the simulation process is generally solved for multiple flapping periods to ensure constancy of the aerodynamic results. In this study, at least five flapping periods and twenty chord lengths of the wake are calculated for each example of a flapping problem. The data of the last complete flapping period are saved to represent the features of a specific case. Based on inviscid flow theory, the Kutta condition is applied to handle the formation of vortices at the trailing edge of the wing model. The numerical simulation of the UVLM begins with the grid generation phase, during which the wing model is divided into N subpanels with respective attached vortices. The equations of the circulation strengths Γ₁, Γ₂, …, Γ_n are solved in every time step.

The induction coefficient from the j^th vortex to the i^th vortex is expressed as a_{i, j}. The symbol RHS_i represents the component of the local stream velocity normal to the i^th subpanel, which is composed of the flow velocity and wake induced velocity. Therefore, the local stream velocity U_i of the i^th subpanel is as follows:where , , and are the velocities induced by the wake in the three axis directions. Due to the relative motion between the flapping wing and the flow, the local velocity also includes the convection velocity expressed as , , and , which makes the boundary condition different from a conventional UVLM. According to the UVLM, the pressures Δp_i on each subpanel are solved by (40) as follows:where τ_X and τ_Y are the unit vectors of the i^th subpanel along the chordwise and spanwise directions. The terms and are used to calculate the projection components of the local stream velocity. and are the circulation strengths of the vortices attached to the subpanels that are connected to the i^th subpanel in the directions of τ_X and τ_Y, respectively. In the situation where there is no connected subpanel in the directions of τ_X and τ_Y, such as solving the pressures at the leading-edge vortices, the terms of and are simply replaced by . By summing up the product of pressures Δp_i and areas ΔS_i of all subpanels, the total lift and drag of the wing are consequently obtained:where the unit vector normal to the i^th subpanel is and the unit vectors in the X-axis and Z-axis directions are and . The corresponding force coefficients are calculated by the following equations:

The average forces of lift and drag are represented by and , which are time averages of the force integrations in a total flapping period. Furthermore, the average force coefficients and are obtained by replacing the L and D in (38) and (39) with and . The average thrust coefficient is defined as .

The convergence of the vortex quantity is determined before solving the fluid-structure interaction. As shown in Table 2, an extreme kinematic parameter is chosen in this case. The flapping frequency and velocity follow the settings in Table 1. Figure 6(a) shows the lift coefficient and thrust coefficient for a complete flapping cycle, while 60, 90, and 120 vortices are assigned to the boundary. Figure 6(b) shows the effect of the length at which the wake is held. Wakes of 15, 20, and 25 chord lengths are tested in Figure 6(b). Before it reaches stabilization, 5 flapping cycles are calculated, and the data in the last cycle are shown in Figure 5. From the result, the convergence of vortex quantity is reached in different cases.

(a)

(b)

4. Fluid-Structure Interaction

The aeroelastic loosely coupled solution is based on a time-stepping solution process in which the equations modelling the dynamic behaviour of both fluid and structure are solved independently with shared boundary information. An interface module based on the surface spline interpolation method [32] is used to transmit the force and displacement information between the flow and the structure. As the behaviour of the fluid and structure is handled in their respective coordinate systems, the aerogrids and the aeroforce need to be expressed correctly before solving another dynamic equation. Figure 7 shows how the information is transmitted in fluid-structure interactions in different coordinate systems.

To solve equation (21), the preparation of modes and kinematics should be done before fluid-structure interaction as shown in Figure 7. Nastran provides the first ten modes in which translations in the x, y, and z directions are used to form Φ. In addition, Nastran provides the mass of all nodes in the x, y, and z directions. According to the result of Nastran, program is prepared to rebuild the mass matrix M of all nodes. K_p is a diagonal matrix including the first ten modal stiffnesses, and they can also be extracted from the results of Nastran. In a simulation case with a given kinematic parameter, such as Table 2, the flapping and pitching motions become explicit in the time domain. The Newmark-β method takes the time step as a known quantity. Thus, the angular velocities and acceleration are computed to construct the replacement items of L and J.

Then, the Newmark-β method guides the process of solving equation (21). Because the initial position and velocity of all nodes are known, the UVLM and the interpolation method calculate the aerodynamic force following the step from interpolate aerogrid displacements to express the aeroforce in the local coordinate in Figure 7. According to equations (24)–(29), and are obtained when the aerodynamic force is solved. Finally, equation (23) can be solved to yield the solution of the next state, and it can update the structural deformation for the next interaction.

5. Global Optimization

Formally, a global optimization problem of the DIRECT algorithm is stated as follows with the assumption that f is the objective function and f_c is the constraint function:where x is the multiparameter vector of the sample point in state space and the bounds of the domain are represented by x_low and x_up. For an optimization problem with given constraints, the penalty-function method is used to handle the infeasible points beyond constraints. Specifically, a high value is artificially assigned to these points. The way the algorithm moves towards the optimum is described in detail in reference [33].

To investigate the kinematic optimization of flapping wing, the propulsion efficiency is chosen as the real-valued objective function corresponding to a set of specific kinematic parameters. The propulsion efficiency η of a flapping wing based on the vortex lattice method is defined as the ratio of the average propulsion power to the average total aerodynamic power :

The average powers are time averages of the power integrations in a total flapping period. The instantaneous powers based on the vortex lattice method are defined as follows:where is the instantaneous velocity of the i^th subpanel divided by the UVLM. By using light weight material, low inertia becomes one of the prominent features of MAVs, including ornithopters. To provide a kinematic strategy for MAV design, inertial power is neglected in the optimization.

The flapping mechanisms of the natural flyer generate enough lift to support the body. Thus, a desirable design specification is that the average lift of the flapping wings is not less than the weight of the flyer. In this article, the average lift of a single wing satisfies the following constraint:

Assigning to the infeasible case beyond constraints is chosen as the strategy of the penalty-function method. The multiparameter vector of the optimization problem consists of the kinematic parameters γ_a, θ₀, θ_a, and φ_θ in equations (1) and (2), in which the bounds of the domain are discussed in the next section. Therefore, the optimization problem is formulated as follows:

6. Results

The variations in the amplitudes and phase angles of the flapping and pitching are considered in the optimization process. Firstly, the upper and lower bounds of these parameters are designed to ensure that the optimal kinematics fall within the range as much as possible. In existing studies [8, 34, 35], the optimal propulsion efficiency is often associated with the Strouhal number St, which is defined as , where h₀ is the amplitude at the wingtip. Schouveiler [34] showed that the maximum propulsion efficiency is achieved at Strouhal numbers from 0.21 to 0.25, while a range from 0.25 to 0.40 was suggested by Anderson et al. [35]. In nature, flying animals [36] operate within a range of Strouhal numbers between 0.2 and 0.4. However, significant discrepancies [35] between experimental measurements and inviscid model predictions occur due to leading-edge separation when the Strouhal number exceeds 0.35. Because the UVLM is an inviscid method, the Strouhal numbers of all simulations are limited to less than 0.35 by restricting the flapping amplitude in the current research. The flapping angles of three other kinds of birds such as seagull, crane, and goose are collected for comparison. Referring to Tianshu’s [37] work, the wings of these birds are modelled as a two-jointed arm with two sections near the wing root and wing tip, in which the rotation angles change within ±40°, as shown in Figure 8.

Moreover, the conditions of the optimal kinematics are from 15° to 25° for the maximum angle of attack and 75° for the phase angle between flapping and pitching in reference [34]. Experimental results [8] show that the pitch-leading phase near 90° is conducive to obtaining the best propulsion efficiency and preventing the occurrence of dynamic stalls. Then, the bounds of and are used to ensure that the pitching angle contains the recommended range, and the constraint of the pitch-leading phase is . Finally, the upper and lower bounds of these parameters are presented in Table 3.

To compare the influence of the deformation on the aerodynamics, a rigid wing and a flexible wing are considered here. In the cases of the flexible model, two simulations are performed, with and without gravity. Twenty iterations in the DIRECT optimization are calculated for each case, and the variation of the maximum propulsion efficiencies with the number of iterations is shown in Figure 9. A total of 411 sets of different kinematic parameters are calculated for the rigid model, and the sets for the flexible model with and without gravity are 401 and 417, respectively. In Figure 9, the maximum propulsion efficiencies reach a plateau after 15 iterations of optimization, which shows that the DIRECT algorithm can quickly converge to the optimal solution in flapping design space.

Table 4 provides a summary of the optimal results for all models as identified by the DIRECT algorithm. The results include the kinematic parameters and average coefficient of lift , average coefficient of thrust , and propulsion efficiency η. Compared with the rigid model, the flexibility of the wing yields a more than 4% increase in the propulsive efficiency and nearly doubled thrust under the same kinematic constraints. In addition, the results of the flexible models with and without gravity are close to each other in both the kinematic parameters and aerodynamic characteristics. The consistency between the two flexible cases is also reflected in the process towards the optimum in Figure 9. Clearly, gravity has hardly any notable impact on the deformation of the current model, which is mainly caused by forced motions. Considering the kinematic parameters, the optimal pitch-leading phases of both the rigid and flexible models are approximately 90°, which again proves that the pitch-leading strategy reported in reference [8] is conducive to obtaining the best propulsion efficiency. Evidence for the pitch-leading phase can also be found in owl flight [38, 39]. The phase shift between flapping and pitching is π/2, which is observed on free-flying barn owls in flapping flight. Although the optimal pitching parameters including θ₀, θ_a, and φ_θ are identified within the design domain, the optimal flapping amplitudes of all cases reach the upper bound of 40°. The Strouhal number is related to the flapping amplitude, which also has a great effect on the wing deformation. In Figure 10, the deflection at the wingtip in a flapping period T is represented by the ratio of displacement d to half span l. According to the results, the maximum deflection at the wingtip in the current cases is 0.23, which reaches the lower bound for nonlinear structural analysis [40]. Thus, the considerations of the leading-edge separation and nonlinear structural dynamics are required when a larger domain of the flapping amplitude is investigated. In addition, the results of the average lift are close to the lift constraint , half the weight of the flyer. This indicates that the optimal motions do not generate excess lift and extract as much thrust as possible from the flow field.

The z deflections at the leading-edge points of the 1/3 span, 2/3 span, and wingtip are plotted in Figure 11 to show the effect of inertia and aerodynamic forces on wing deformation. The flexible model without gravity is chosen in this case. The deflection is smaller at positions closer to the clamped boundary at the wing root.

Figure 11 shows that inertia dominates the wing deformation, while aerodynamic forces have a slight effect on the deformation and a certain advance for the phase. The deformation curves at different spans intersect twice in the downstroke and upstroke stages. The advance phases measured in Figure 11 are approximately 15° during the downstroke and 4° during the upstroke.

The angle-time and force-time plots are shown in Figures 12 and 13, respectively. Most of the positive lift and thrust are produced during the downstroke, and the force peak occurs near the 0.3 period when the pitching angle reverses in the downstroke. In the upstroke, the optimal strategy in all cases keeps the lift and thrust at lower values. Moreover, the lift remains positive during almost the entire period while the thrust returns negative in the upstroke. Clearly, the thrust peaks of the flexible models are higher than those of the rigid model, which enables a significant increase in the propulsive efficiency.

(a)

(b)

The pressure and force distributions for the optimal flight paths predicted by the DIRECT algorithm are exhibited for the rigid and flexible wings in Figure 14. The passive deformation dominated by the first bending mode is clearly observed near the wingtip in the flexible case. The pressure variation between the rigid and flexible wings is also shown near the wingtip during the downstroke. An increase in wingtip amplitude improves the speed to slice the air flow near the middle plane of the downstroke. Based on equations (38)–(40), higher circulation strengths are obtained when the local stream velocity increases, which results in higher pressure on the wing. This indicates that the flexible wing has the potential to increase aerodynamics by exploiting spanwise deformation although chordwise flexibility is not reflected in the current model. Furthermore, an interesting relation is found between the force distributions in Figure 14 and the planform of the natural bird wing in Figure 1. During the downstroke of the numerical model, the positive and negative forces centre at the leading edge of the outer wing section and wing root, respectively. The planform of the natural wing shows preservation and the absence of material in these wing sections. This indicates that the aerodynamic characteristics of the wing can be further improved by designing the planform.

(a)

(b)

7. Conclusions

In this study, kinematic global optimization of a flexible wing model is performed on a loosely coupled aeroelastic framework to investigate the aerodynamic effect of passive deformation caused by flapping and pitching motion. The propulsion efficiency is selected as the optimization goal to solve the flight strategy of flapping flyers with a high-aspect ratio configuration. According to the results of the simulation and optimization, the following conclusions are obtained:(1)The flapping wing with proper spanwise flexibility has more aerodynamic potential than the rigid model in a wide design domain of kinematics. Although the advantage of proper spanwise flexibility has been proven in existing experiments, the current work excludes the influence of the selection of kinematic strategy on this conclusion.(2)The elastic deformation of the flexible model doubles the thrust with the same flapping angle as the rigid model. The variation in the force distributions is mainly present at the leading edge of the outer wing section during the downstroke. A small area of negative pressure appears at the wing root while pitching at a negative angle. Therefore, a layout with a trim near the leading edge at the wing root is suggested.(3)Considerations of the leading-edge flow separation and nonlinear structural dynamics are required if a larger range of flapping angles and a more flexible model are used, which goes beyond the application scope of the present aeroelastic analytical framework.(4)The DIRECT algorithm exhibits fast optimizing characteristics on the multiparameter kinematic problem with constraints. Suitable for searching the optimal value of the simulation result, the DIRECT algorithm maximizes the thrust efficiently while satisfying the required lift.(5)The effect of gravity is negligible on the model constructed following the present data, and the deformation can be explained by a linear model. The large taper ratio and proper camber provide sufficient stiffness for large-scale forced movement. The excitation frequency of flapping and pitching is half of the first-order bending frequency, which is far from the resonant frequency.

Data Availability

The graphic and tabular data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

References

S. Heathcote, Z. Wang, and I. Gursul, “Effect of spanwise flexibility on flapping wing propulsion,” Journal of Fluids and Structures, vol. 24, no. 2, pp. 183–199, 2008.
View at: Publisher Site | Google Scholar
S. Heathcote and I. Gursul, “Flexible flapping airfoil propulsion at low Reynolds numbers,” in Proceedings of the43rd AIAA Aerospace Sciences Meeting and Exhibit, AIAA, Reno, NV, USA, January 2005.
View at: Publisher Site | Google Scholar
K. Deepa, K. K. Chang, and A. Hikaru, “An analytical aeroelastic model of spanwise flexible flapping wings in forward flight,” in Proceedings of the 55th AIAA Aerospace Sciences Meeting, Grapevine, TX, USA, January 2017.
View at: Google Scholar
Y. Wenqing, S. Bifeng, S. Wenping, and W. Liguang, “The effects of span-wise and chord-wise flexibility on the aerodynamic performance of micro flapping-wing,” Chinese Science Bulletin, vol. 57, no. 22, pp. 2887–2897, 2012.
View at: Publisher Site | Google Scholar
Z. Qiang, “Numerical simulation of a flapping foil with chordwise or spanwise flexibility,” AIAA Journal, vol. 45, no. 10, pp. 2448–2457, 2007.
View at: Publisher Site | Google Scholar
K. Matthew, K. Karl, and W. Henry, “Development of the nano hummingbird: a tailless flapping wing micro air vehicle,” in Proceedings of the 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Nashville, TN, USA, January 2012.
View at: Google Scholar
T. N. S. Pornsin, Y. C. Tai, C. M. Ho, and M. Keennon, “Microbat: a palm-sized electricallypowered ornithopter,” in Proceedings of the Nasa/jpl Workshop on Biomorphic Robotics, Pasadena, CA, USA, August 2001.
View at: Google Scholar
N. A. Razak and G. Dimitriadis, “Experimental study of wings undergoing active root flapping and pitching,” Journal of Fluids and Structures, vol. 49, no. 1, pp. 687–704, 2014.
View at: Publisher Site | Google Scholar
S. Sutthiphong and W. L. Chan, “Flow field of flapping albatross-like wing and sound at low Reynolds number,” Journal of Unmanned System Technology, vol. 1, no. 2, 2013.
View at: Google Scholar
K. S. Bret and S. B. Philip, “Analytical sensitivity analysis of an unsteady vortex-lattice method for flapping-wing optimization,” Journal of Aircraft, vol. 47, no. 2, pp. 647–662, 2010.
View at: Google Scholar
W. Shyy, H. Aono, and S. K. Chimakurthi, “Recent progress in flapping wing aerodynamics and aeroelasticity,” Progress in Aerospace Sciences, vol. 46, no. 7, pp. 284–327, 2010.
View at: Publisher Site | Google Scholar
C. D. Perttunen, D. R. Jones, and B. E. Stuckman, “Lipschitzian optimization without the lipschitz constant,” Journal of Optimization Theory and Applications, vol. 79, no. 1, pp. 157–181, 1993.
View at: Google Scholar
S. Deng, C. Jiang, Y. Wang, and H. Wang, “Acceleration of unsteady vortex lattice method via dipole panel fast multipole method,” Chinese Journal of Aeronautics, vol. 34, no. 2, pp. 265–278, 2021.
View at: Publisher Site | Google Scholar
J. Katz and A. Plotkin, Low Speed Aerodynamics: From wing Theory to Panel Methods, McGraw-Hill, New York, NY, USA, 1991.
G. Mehdi, R. H. Muhammad, and T. M. Dean, “Global optimization of actively morphing flapping wings,” Journal of Fluids and Structures, vol. 33, no. 1, pp. 210–228, 2012.
View at: Google Scholar
M. Ghommem, N. Collier, A. H. Niemi, and V. M. Calo, “Shape optimisation and performance analysis of flapping wings,” in Proceedings of the 8th International Conference on Engineering Computational Technology, Civil-Comp, Stirlingshire, Scotland, September 2012.
View at: Google Scholar
J. D. Eldredge and A. R. Jones, “Leading-edge vortices: mechanics and modeling,” Annual Review of Fluid Mechanics, vol. 51, pp. 75–104, 2019.
View at: Publisher Site | Google Scholar
Y. Wang, X. Zhao, R. Palacios, and K. Otsuka, “Unsteady aeroelasticity of slender wings with leading-edge separation,” in Proceedings of the AIAA Scitech Forum, Nashville, TN, USA, January 2021.
View at: Google Scholar
X. Xia and K. Mohseni, “A flat plate with unsteady motion effect of angle of attack on vortex shedding,” in Proceedings of the AIAA SciTech Forum: 52nd Aerospace Sciences Meeting, National Harbor, Maryland, January 2014.
View at: Publisher Site | Google Scholar
X. Xia and K. Mohseni, “Lift evaluation of a two-dimensional pitching flat plate,” Physics of Fluids, vol. 25, Article ID 091901, 2013.
View at: Publisher Site | Google Scholar
C. J. Pennycuick, “The flight of petrels and albatrosses (procellariiformes), observed in south Georgia and its vicinity,” Philosophical Transactions of the Royal Society of London - Series B: Biological Sciences, vol. 300, no. 1908, pp. 75–106, 1982.
View at: Google Scholar
C. J. Pennycuick, “Wingbeat frequency of birds in steady cruising flight: new data and improved predictions,” Journal of Experimental Biology, vol. 199, no. 7, pp. 1613–1618, 1996.
View at: Publisher Site | Google Scholar
S. A. Combes and T. L. Daniel, “Flexural stiffness in insect wings I. Scaling and the influence of wing venation,” Journal of Experimental Biology, vol. 206, no. 17, pp. 2979–2987, 2003.
View at: Publisher Site | Google Scholar
K. Deepa, M. Cory, K. K. Chang, and H. Aono, “Effects of spanwise flexibility on the performance of flapping flyers in forward flight,” Journal of The Royal Society Interface, vol. 14, no. 136, pp. 1–14, 2017.
View at: Publisher Site | Google Scholar
R. D. B. Fraser and T. P. Macrae, “Molecular structure and mechanical properties of keratins,” Symposia of the Society for Experimental Biology, vol. 34, no. 34, pp. 211–246, 1980.
View at: Google Scholar
H. Tennekes, The Simple Science of Flight - from Insects to Jumbo Jets, MIT, Cambridge, MA, USA, 1996.
M. Hassanalian, G. Throneberry, and M. Ali, “Role of wing color and seasonal changes in ambient temperature and solar irradiation on predicted flight efficiency of the Albatross,” Journal of Thermal Biology, vol. 71, pp. 112–122, 2018.
View at: Publisher Site | Google Scholar
R. Todd, Kinematic Optimization in Birds, Bats and Ornithopters, University of Toronto, Toronto, Canada, 2011, Institute for aerospace studies [dissertation].
H. J. Hung, L. J. Seong, and K. D. Kwan, “Bio-inspired flapping UAV design: a university perspective,” in Proceedings of the of SPIE Smart Structures and Materials + Nondestructive Evaluation and Health Monitoring, vol. 7295, no. 1, pp. 1–12, San Diego, CA, USA, April 2009.
View at: Publisher Site | Google Scholar
M. J. C. Smith, “Simulating moth wing aerodynamics: towards the development of flapping-wing technology,” AIAA Journal, vol. 34, no. 7, pp. 1348–1355, 1996.
View at: Publisher Site | Google Scholar
K. Isogai and Y. Harino, “Optimum aeorelastic design of a flapping wing,” Journal of Aircraft, vol. 44, no. 6, pp. 2040–2048, 2007.
View at: Publisher Site | Google Scholar
X. Changchuan and Y. Chao, “Surface splines generalization and large deflection interpolation,” Journal of Aircraft, vol. 44, no. 3, pp. 1024-1025, 2007.
View at: Publisher Site | Google Scholar
J. M. Gablonsky, Direct Version 2.0 Userguide, North Carolina State University: Center for Research in Scientific Computation, Raleigh, North Carolina, 2001, Report No.: CRSC-TR01-08.
L. Schouveilera, F. S. Hover, and M. S. Triantafyllou, “Performance of flapping foil propulsion,” Journal of Fluids and Structures, vol. 20, no. 7, pp. 949–959, 2005.
View at: Publisher Site | Google Scholar
J. M. Anderson, K. Streitlien, and D. S. Barrett, “Oscillating foils for high propulsive efficiency,” Journal of Fluid Mechanics, vol. 360, pp. 41–72, 1998.
View at: Publisher Site | Google Scholar
M. F. Platzer, K. D. Jones, Y. John, and C. S. L. Joseph, “Flapping-wing aerodynamics: progress and challenges,” AIAA Journal, vol. 46, no. 9, pp. 2136–2149, 2008.
View at: Publisher Site | Google Scholar
L. Tianshu, K. Kuykendoll, R. Rhew, and S. Jones, “Avian wing geometry and kinematics,” AIAA Journal, vol. 44, no. 5, pp. 954–963, 2006.
View at: Publisher Site | Google Scholar
T. Wolf and R. Konrath, “Avian wing geometry and kinematics of a free-flying barn owl in flapping flight,” Experiments in Fluids, vol. 56, no. 28, 2015.
View at: Publisher Site | Google Scholar
L. Xinyu, S. Bifeng, Y. Wenqing, and W. Song, “Aerodynamic performance of owl-like airfoil undergoing bio-inspired flapping kinematics,” Chinese Journal of Aeronautics, vol. 34, no. 5, pp. 239–252, 2021.
View at: Publisher Site | Google Scholar
A. Frederico, V. José, O. Éder, F. Lau, and A. Suleman, “A review on non-linear aeroelasticity of high aspect-ratio wings,” Progress in Aerospace Sciences, vol. 89, pp. 40–57, 2017.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2021 Yue Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

510

Downloads

547

Citations