# Kinematic Optimization of a Flexible Wing Undergoing Flapping and Pitching

Authors: Yue Wu; Changchuan Xie; Yang Meng; Chao Yang
Journal: Shock and Vibration
 (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 × 105, 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.


---

## Body

## 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].Table 1
Wing morphological parameters of the numerical albatross-wing model.
ParameterThe albatross-wing modelHalf span length [21], l (m)1.08Wing area [21], S (m2)0.178Mean chord length,c = S/l (m)0.1648Aspect ratio for one wing [21], AR6.55Weight of the flyer [22], W (N)30.184Density of the wing [23], ρs (kg/m3)1200Density of air [24], ρf (kg/m3)1.23Young’s modulus [25], E (GPa)5.2Flapping frequency [22], f (Hz)2.97Velocity [26], U∞ (m/s)16.7Reduced frequency,k = πfc/U∞0.092Reynolds number,Re1.9 × 105Further 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 hs = 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.Figure 1
Layout of the numerical and natural wings.Figure 2
Views of the dimensionless albatross airfoil.


### 2.2. Wing Kinematics

This article considers a three-dimensional wing in constant wind speedU∞ 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:(1)yt=yacos2Πft,(2)θt=θ0+θacos2Πft+θ0.Figure 3
Coordinate systems of a flapping and pitching wing.The motions with constant frequencyf contain the adjustable kinematic parameters of flapping amplitude γa, middle pitching angle θ0, pitching amplitude θa, and pitch-leading phase φθ.In Figure3, 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 ϕ=γθ0°T (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 ωB=pqrT. According to the projecting relationships, ωB is expressed by the derivative of ϕ with time as follows:(3)ωB=DBϕ˙,where the mark with one dot above the symbol represents the first derivative of time and DB is the transition matrix as follows:(4)DB=cosθ0−cosysinθ01sinysin0cosycosθ.Thus, equation (3) has the following expansion form:(5)P=y˙cosθ,q=θ˙,r=y˙sinθ.


### 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 energyT of an arbitrary microscale unit of the wing, the nondeformation position vector and the displacement vector of a microscale unit are defined as rB=r1r2r3T and uB=u1u2u3T 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:(6)T=12u˙BTu˙dm−12ωBTr˜B+u˜Br˜B+u˜BωBdm+ωBTr˜B+u˜Bu˙Bdm,where the marks with a curve over the symbols represent the second-order tensors that are defined as follows:(7)r˜B=0−r3r2r30−r1−r2r10,(8)u˜B=0−u3u2u30−u3−u2u20.Assuming thatk is the stiffness matrix and Ue is the elastic potential energy of the microunit, we have(9)∂∪e∂uB=kuB.According to the Lagrange equation, the kinetic energyT, the elastic potential energy Ue, and the force vector f on the microscale unit are related via the following equations:(10)ddt∂T∂u˙B−∂T∂uB+∂ue∂uB=f,which is then reduced to the simple Ax¨+Bx˙+Cx=D form by substituting (6) and (9) into (10).(11)u¨Bdm+2ω˜Bu˙Bdm+ω˙˜Bdm+ω˜Bω˜Bdm+kuB=f−ω˙˜+ω˜Bω˜BrBdm,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:(12)MU¨B+2MΩU˙B+MΩ˙+MΩΩ+kUB=F−MΩ˙+ΩΩRB.The symbols in the total elastodynamic equation are defined as follows:(13)UB=uB1…uBnT,(14)RB=rB1…rBnT,(15)M=m1⋯0⋮⋱⋮0…mn,mi=mi100010001,i=1,…,n,(16)K=k1…0⋮⋱⋮0⋯kn,Ω=ω˜B⋯0⋮⋱⋮0⋯ω˜B,F=f1…fnT.Based on the modal method, the displacementUB is the sum of the modes Φ weighted by the generalized coordinate q:(17)UB=Φq.Substituting equation (17) into dynamic equilibrium equation (12) yields the following result:(18)MΦq¨+2MΩΦq+MΩ˙+MΩΩ+KΦq=F−MΩ˙+ΩΩRB.The modal mass matrix is normalized byMp=ΦTMΦ=I3n×3n, and the modal stiffness matrix is defined as Kp=ΦTKΦ. The replacements of the terms in equation (18) are performed as follows:(19)MΩ=L,(20)MΩ˙+ΩΩJ.Thus, equation (18) is arranged as follows:(21)q¨+2ΦTLΦq˙+ΦTJΦ+Kpq=ΦTF−JRB.As an implicit method, the Newmark-β method is often used for solving structural dynamics problems discretized by time steps. In the step nt + 1, the dynamic equation is shown as follows with the generalized coordinate and its time derivatives:(22)Aq¨n+1+Bq˙nt+1+Cqnt+1=D.By assuming the stepping relation between time steps, the Newmark-β method transforms the dynamic eqaution (22) into a simple form as follows:(23)C^qnt+1=D^,where(24)C^=C+a0A+a1B,(25)D^=D+Aa3q¨nt+a2q˙nt+a0qnt+Ba5q¨nt+a4q˙nt+a1qnt.When the length of the time step is determined, the coefficientsai,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:(26)A=I3n×3n,(27)B=2ΦTLΦ,(28)C=ΦTJΦ+Kp,(29)D=ΦTF−JRB.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 Figure4, 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.Figure 4
Vibrations and frequencies of the first four modes.


## 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].Table 1
Wing morphological parameters of the numerical albatross-wing model.
ParameterThe albatross-wing modelHalf span length [21], l (m)1.08Wing area [21], S (m2)0.178Mean chord length,c = S/l (m)0.1648Aspect ratio for one wing [21], AR6.55Weight of the flyer [22], W (N)30.184Density of the wing [23], ρs (kg/m3)1200Density of air [24], ρf (kg/m3)1.23Young’s modulus [25], E (GPa)5.2Flapping frequency [22], f (Hz)2.97Velocity [26], U∞ (m/s)16.7Reduced frequency,k = πfc/U∞0.092Reynolds number,Re1.9 × 105Further 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 hs = 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.Figure 1
Layout of the numerical and natural wings.Figure 2
Views of the dimensionless albatross airfoil.


## 2.2. Wing Kinematics

This article considers a three-dimensional wing in constant wind speedU∞ 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:(1)yt=yacos2Πft,(2)θt=θ0+θacos2Πft+θ0.Figure 3
Coordinate systems of a flapping and pitching wing.The motions with constant frequencyf contain the adjustable kinematic parameters of flapping amplitude γa, middle pitching angle θ0, pitching amplitude θa, and pitch-leading phase φθ.In Figure3, 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 ϕ=γθ0°T (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 ωB=pqrT. According to the projecting relationships, ωB is expressed by the derivative of ϕ with time as follows:(3)ωB=DBϕ˙,where the mark with one dot above the symbol represents the first derivative of time and DB is the transition matrix as follows:(4)DB=cosθ0−cosysinθ01sinysin0cosycosθ.Thus, equation (3) has the following expansion form:(5)P=y˙cosθ,q=θ˙,r=y˙sinθ.


## 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 energyT of an arbitrary microscale unit of the wing, the nondeformation position vector and the displacement vector of a microscale unit are defined as rB=r1r2r3T and uB=u1u2u3T 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:(6)T=12u˙BTu˙dm−12ωBTr˜B+u˜Br˜B+u˜BωBdm+ωBTr˜B+u˜Bu˙Bdm,where the marks with a curve over the symbols represent the second-order tensors that are defined as follows:(7)r˜B=0−r3r2r30−r1−r2r10,(8)u˜B=0−u3u2u30−u3−u2u20.Assuming thatk is the stiffness matrix and Ue is the elastic potential energy of the microunit, we have(9)∂∪e∂uB=kuB.According to the Lagrange equation, the kinetic energyT, the elastic potential energy Ue, and the force vector f on the microscale unit are related via the following equations:(10)ddt∂T∂u˙B−∂T∂uB+∂ue∂uB=f,which is then reduced to the simple Ax¨+Bx˙+Cx=D form by substituting (6) and (9) into (10).(11)u¨Bdm+2ω˜Bu˙Bdm+ω˙˜Bdm+ω˜Bω˜Bdm+kuB=f−ω˙˜+ω˜Bω˜BrBdm,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:(12)MU¨B+2MΩU˙B+MΩ˙+MΩΩ+kUB=F−MΩ˙+ΩΩRB.The symbols in the total elastodynamic equation are defined as follows:(13)UB=uB1…uBnT,(14)RB=rB1…rBnT,(15)M=m1⋯0⋮⋱⋮0…mn,mi=mi100010001,i=1,…,n,(16)K=k1…0⋮⋱⋮0⋯kn,Ω=ω˜B⋯0⋮⋱⋮0⋯ω˜B,F=f1…fnT.Based on the modal method, the displacementUB is the sum of the modes Φ weighted by the generalized coordinate q:(17)UB=Φq.Substituting equation (17) into dynamic equilibrium equation (12) yields the following result:(18)MΦq¨+2MΩΦq+MΩ˙+MΩΩ+KΦq=F−MΩ˙+ΩΩRB.The modal mass matrix is normalized byMp=ΦTMΦ=I3n×3n, and the modal stiffness matrix is defined as Kp=ΦTKΦ. The replacements of the terms in equation (18) are performed as follows:(19)MΩ=L,(20)MΩ˙+ΩΩJ.Thus, equation (18) is arranged as follows:(21)q¨+2ΦTLΦq˙+ΦTJΦ+Kpq=ΦTF−JRB.As an implicit method, the Newmark-β method is often used for solving structural dynamics problems discretized by time steps. In the step nt + 1, the dynamic equation is shown as follows with the generalized coordinate and its time derivatives:(22)Aq¨n+1+Bq˙nt+1+Cqnt+1=D.By assuming the stepping relation between time steps, the Newmark-β method transforms the dynamic eqaution (22) into a simple form as follows:(23)C^qnt+1=D^,where(24)C^=C+a0A+a1B,(25)D^=D+Aa3q¨nt+a2q˙nt+a0qnt+Ba5q¨nt+a4q˙nt+a1qnt.When the length of the time step is determined, the coefficientsai,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:(26)A=I3n×3n,(27)B=2ΦTLΦ,(28)C=ΦTJΦ+Kp,(29)D=ΦTF−JRB.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 Figure4, 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.Figure 4
Vibrations and frequencies of the first four modes.


## 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.Figure 5
Vortex arrangement on the wing model.In the coordinate system shown in Figure3, the induced velocity caused by the attached and free vortices is represented by the vector X˙,Y˙,Z˙. 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:(30)ΔXΔYΔZT=x˙y˙z˙X˙Y˙Z˙T+U∞Δt.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 intoN subpanels with respective attached vortices. The equations of the circulation strengths Γ1, Γ2, …, Γn are solved in every time step.(31)a11a12⋯a1Na21a22⋯a2Na31a32⋯a3N⋮⋮⋮⋮aN1aN2⋯aNNΓ1Γ2Γ3⋮ΓN=RHS1RHS2RHS3⋮RHSN.The induction coefficient from thejth vortex to the ith vortex is expressed as ai, j. The symbol RHSi represents the component of the local stream velocity normal to the ith subpanel, which is composed of the flow velocity and wake induced velocity. Therefore, the local stream velocity Ui of the ith subpanel is as follows:(32)Ui=X˙+X˙w+U∞−y˙+y˙w+Z˙+Z˙wT,where X˙w, Y˙w, and Z˙w 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 −X˙, −Y˙, and −Z˙, which makes the boundary condition different from a conventional UVLM. According to the UVLM, the pressures Δpi on each subpanel are solved by (40) as follows:(33)Δpi=ρfΓi−ΓXΔτXUi.τX+Γi−ΓYΔτYUi.τY+ddtΓi,where τX and τY are the unit vectors of the ith subpanel along the chordwise and spanwise directions. The terms Ui,τX and Ui,τY are used to calculate the projection components of the local stream velocity. ΓX and ΓY are the circulation strengths of the vortices attached to the subpanels that are connected to the ith 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 Γi−ΓX and Γi−ΓY are simply replaced by Γi. By summing up the product of pressures Δpi and areas ΔSi of all subpanels, the total lift and drag of the wing are consequently obtained:(34)L=∑i=1NΔpiΔSini⋅nXOY,(35)D=∑i=1NΔpiΔSini⋅nZOY,where the unit vector normal to the ith subpanel is ni and the unit vectors in the X-axis and Z-axis directions are nZOY and nXOY. The corresponding force coefficients are calculated by the following equations:(36)CL=2LρfU∞2S,(37)CD=2DρfU∞2S.The average forces of lift and drag are represented byL¯ and D¯, which are time averages of the force integrations in a total flapping period. Furthermore, the average force coefficients C¯L and C¯D are obtained by replacing the L and D in (38) and (39) with L¯ and D¯. The average thrust coefficient is defined as C¯T=−C¯D.The convergence of the vortex quantity is determined before solving the fluid-structure interaction. As shown in Table2, 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.Table 2
Kinematic parameter in the UVLM convergence.
Wing motionParameterAngle (deg)Flapping amplitudeγa40°Middle pitching angleθ010°Pitching amplitudeθa20°Pitch-leading phaseφθ90°Figure 6
The lift coefficient and thrust coefficient in a complete cycle: (a) 60, 90, and 120 bound vortices; (b) 15, 20, and 25 chord-length wakes.
(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.Figure 7
A schematic of the aeroelastic framework for fluid-structure interaction.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. Kp 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), C^ and D^ 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 thatf is the objective function and fc is the constraint function:(38)minfx,subject toxlow≤x≤xupandfcx≤0,where x is the multiparameter vector of the sample point in state space and the bounds of the domain are represented by xlow and xup. 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 P¯out to the average total aerodynamic power P¯in:(39)η=p¯outp¯in.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:(40)pout=−DUα,(41)pin−∑i=1NΔpiΔsini⋅vmotioni,where vmotioni is the instantaneous velocity of the ith 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 liftL¯ of a single wing satisfies the following constraint:(42)L¯≥W2.Assigningη=−1 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, θ0, θ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:(43)minηx,subject toxlow≤x≤xupandL¯≥W2.


## 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 St=2h0f/U∞, where h0 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 γa≤40° 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.Figure 8
Rotations at the wing root and wing tip of seagull, crane, and goose.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 0°≤θ0≤10° and 0°≤θa≤20° are used to ensure that the pitching angle contains the recommended range, and the constraint of the pitch-leading phase is 75°≤φθ≤105°. Finally, the upper and lower bounds of these parameters are presented in Table 3.Table 3
Kinematic parameter constraints.
Wing motionParameterLower bound (deg)Upper bound (deg)Flapping amplitudeγa0°40°Middle pitching angleθ00°10°Pitching amplitudeθa0°20°Pitch-leading phaseφθ75°105°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 Figure9. 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.Figure 9
Maximum propulsion efficiency vs. the number of the optimization iterations.Table4 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 C¯L, average coefficient of thrust C¯T, 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 θ0, θ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 W/2=15.092N, 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.Table 4
Optimal kinematics results from the DIRECT algorithm.
Rigid modelFlexible model (with gravity)Flexible model (without gravity)Kinematic parameter (deg, °)γa39.97°39.48°39.75°θ03.340°3.395°3.395°θa11.07°13.29°13.32°φθ86.91°83.83°83.74°Aerodynamic characteristicsC¯L0.50930.49420.4942C¯T0.02930.05710.0589Η0.13520.17730.1792Average force (N)L¯15.555 N15.095 N15.095 ND¯－0.893 N－1.745 N－1.799 NFigure 10
Wingtip deflection vs. time in a period.Thez 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
The deflections at the leading-edge points of the 1/3 span, 2/3 span, and wingtip.Figure11 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 Figures12 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.Figure 12
Degrees of flapping angleγ and pitching angle θ varying with time in a period: (a) lift coefficient; (b) drag coefficient.Figure 13
Lift coefficient and drag coefficient vs. time in a period.
(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 Figure14. 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.Figure 14
(a) Pressure distributions (N·m2) and (b) force distributions (N) on rigid and flexible wings in a period.
(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.



---

*Source: 9635584-2021-12-15.xml*