Computational Aeroelasticity of Flying Robots with Flexible Wings

Abstract

A computational co‐simulation framework for flying robots with flexible wings is presented. The authors combine a nonlinear aerodynamic model based on an extended version of the unsteady vortex‐lattice method with a nonlinear structural model based on a segregated formulation of Lagrange’s equations obtained with the Floating Frame of Reference formalism. The structural model construction allows for hybrid combinations of different models typically used with multibody systems such as models based on rigid‐body dynamics, assumed‐modes techniques, and finite‐element methods. The aerodynamic model includes a simulation of leading‐edge separation for large angles of attack. The governing differential‐algebraic equations are solved simultaneously and interactively to obtain the structural response and the flow in the time domain. The integration is based on the fourth‐order predictor‐corrector method of Hamming with a procedure to stabilize the iteration. The findings are found to capture known nonlinear behavior of flapping-wing systems. The developed framework should be relevant for conducting aeroelastic studies on a wide variety of air vehicle systems.

1. Introduction

Natural flight is characterized by highly nonlinear and unsteady flows and wing‐deformation patterns that include time‐dependent in‐and‐out‐of‐plane bending, changing camber, spanwise twisting, expanding and contracting wing surfaces, and so on. A detailed description of the deformations of wings in different flight configurations is crucial for quantifying the influence of a wing’s flexibility on the aerodynamic forces. Flight kinematics and aerodynamics have been identified and studied for a variety of insect and bird wings [1–13]. Advances in computational power have further enabled studies of highly flexible structures such as flapping and morphing wings of micro‐air vehicles (MAVs). The primary components of numerical simulations are the aerodynamic solver, the structural‐dynamic solver, and the communication between them. Early studies on fluid‐structure interactions for vehicles with flapping wings have been based on either two‐dimensional models or other simplified aerodynamic models based on thin‐airfoil theory or unsteady panel methods or Euler methods and on linear beam finite elements [14–17].

Nakata and Liu [18] used an unsteady aerodynamic model based on a Navier‐Stokes solver [19] coupled with a finite‐element formulation specifically developed for insect flight that accounts for the distribution and anisotropy of the wing’s veins and membranes. Chimakurthi et al. [20] proposed a computational framework, in which a nonlinear finite‐shell‐element solver is coupled with a Navier‐Stokes flow solver [21]. Malhan et al. [22] analyzed flexible flapping wings by using a computational fluid dynamics‐computational structural‐dynamics (CFD‐CSD) environment. They used a Reynolds‐averaged Navier‐Stokes (RANS) solver [23] along with an open‐source multibody software [24]. Unger et al. [25] investigated the flow around a flexible airfoil based on a seagull’s wing. They used a flow model based on the incompressible unsteady RANS equations (URANS) [26]. To model the structure and obtain the response, the authors generated a finite‐element model that includes shells, solids, and contact elements for modeling the various structural members of a wing, and used an ANSYS [27] solver. Recently, Bose et al. [28] solved the incompressible two‐dimensional flow field for a range of solid‐to‐fluid added‐mass ratios (characteristic of natural flappers or MAVs) by using an arbitrary Lagrangian‐Eulerian (ALE) formulation and a nonlinear discrete model of the structure. They coupled the models through a fourth‐order Runge‐Kutta scheme implemented in OpenFOAM [29], and found that the structural response experienced a supercritical Hopf bifurcation [30].

Computational difficulties and costs associated with the use of CFD models have led to alternative approaches. For example, unsteady vortex‐lattice methods (UVLMs) are an effective and a simple choice [31–33]. Related to flapping wings, Taha et al. [34] identified five main contributors to flow characteristics during hover: the wing’s translation and rotation, the leading‐edge vortices (LEVs), wake capture, viscosity, and added‐mass effects. UVLMs can be used to model all aspects except viscous and LEV effects. As shown by the experiments of Dickinson et al. [1], over the range of Reynolds numbers (75–4000) of hovering insects, viscous effects can be neglected; this makes the use of UVLMs suitable for modeling the aerodynamics of flapping wings and morphing wings.

The authors extend the approach proposed by Preidikman and Mook [32] and Preidikman [35]; a computational aeroelastic framework is presented. In these prior studies, the structures were much simpler than those considered here, and they were represented by means of traditional beam formulations based on assumed‐modes/finite‐element methods. In the authors’ novel methodology, the dynamical system is partitioned into two subsystems (the structural and aerodynamic models) that communicate with each other across the boundary of the flow field (the surface of the structure) in a strong way. The computational environment can be considered a strong co‐simulation framework [36]. Initial steps taken in this direction are presented in reference [37]. The underlying methodology is detailed in Section 2. The aerodynamic and structural models are presented in the third and fourth sections, respectively. The communication between the two models and the numerical integration scheme are presented in Section 5. In the next section, results obtained for systems with flapping and morphing wings are presented. Finally, concluding remarks are collected together and presented.

2. Methodology

The structural model (Simulator 1) and the aerodynamic model (Simulator 2) are different subsystems, which are constructed to exchange information bi‐directionally in an iterative sequence that continuously improves the estimation of the aerodynamic loads (see Figure 1) and, consequently, the structure’s response. The numerical scheme used by Simulator 2 is based on the UVLM, and the one used by Simulator 1 to solve the vehicle’s equations of motion is based on Hamming’s fourth‐order predictor‐corrector method [35]. This procedure was chosen for two reasons: (i) performance of Simulator 2 is better if the loads are evaluated at integral time steps and (ii) aerodynamic loads depend on the structure’s acceleration. To combine the simulators, a technique for transferring information (TTI) based on radial basis functions is used.

Figure 1.

Schematic diagram of fluid‐structure interaction framework [37].

The substantial deformations of the lifting surfaces induce significant changes in the aerodynamic loads, which in turn induce further changes in the deformation of the wings. This feedback between the aerodynamic loads and the movement of structural members (flowfield boundaries) generates strong coupling between the aerodynamic and structural models. In order to capture these complex interactions with a numerical model, a strong coupling method is needed; the authors have chosen a so‐called two‐way non‐monolithic strategy. Although Simulators 1 and 2 are computational implementations of independently modeled physical components, the coupling procedure is strong because information can be bi‐directionally exchanged, and the chosen time step advance is the same for both models. Traditional approaches used to study unsteady aeroelastic characteristics may not be adequate for modeling the unusual fluid‐structure interactions that characterize aerial vehicles with flexible flapping wings. Classical linear analyses can, in some instances, predict when aeroelastic systems become unstable, but they can neither detect subcritical instabilities nor describe post‐instability motions. As the solutions are obtained in the time domain here, the authors have the capability for predicting both subcritical instabilities and post‐instability motions [30].

3. Aerodynamic model

Here, UVLM is used to compute the aerodynamic loads. It can be applied to lifting surfaces of any planform, camber, and twist. The lifting surface may undergo any time‐dependent deformation and execute any maneuver in moving air. The flow surrounding the lifting surface is assumed to be inviscid, incompressible, and irrotational over the entire flowfield, except at the solid boundaries and in the wakes. Due to the relative motion between the wing and the fluid and due to the viscous effects, vorticity is generated in the fluid in a thin region next to the wing’s surface (the boundary layer). The boundary layers on the upper and lower surfaces are merged into a single vortex sheet. Vorticity in this sheet must be shed from the sharp edges (trailing edge, leading edge, and wing tips) so that the pressures are equal in the merging flows coming off the upper and lower sides of the lifting surface, the so‐called Kutta condition. This shed vorticity forms the wake as it moves with the fluid in order for the pressure to remain continuous across the wake, which is force‐free. The calculation of the velocity field associated with the vorticity is expedited by replacing the continuous vortex sheet with a lattice of discrete vortex lines. The boundary layers are represented by bound‐vortex lattices, which are attached or bound to the lifting surface and move with it. The wakes are represented by free‐vortex lattices, which move freely with the flow so that their position and vorticity distribution are determined as part of the solution [35, 37–39].

3.1. Mathematical formulation

The velocity of a fluid particle that occupies position, $$\mathit{r}$$ , at instant $$t$$ is denoted by $$\mathit{V}{\displaystyle \left(\mathit{r};t\right)}$$ . Since the flow is irrotational outside the boundary layers and wakes and considered to be incompressible, the velocity field can be expressed as the gradient of a scalar potential $$\varphi {\displaystyle \left(\mathit{r};t\right)},$$ which is governed by the continuity equation:

${\nabla }^2\varphi {\displaystyle \left(\mathit{r};t\right)}=0.$ E1

In a three‐dimensional flow, the velocity field satisfying Eq. (1) and associated with a straight, finite‐length segment of a vortex line with circulation $\Gamma {\displaystyle \left(t\right)}$ is given by the Biot‐Savart law:

$\mathit{V}{\displaystyle \left(\mathit{r},t\right)}=\frac{\Gamma {\displaystyle \left(t\right)}}{4\pi }\frac{\mathit{L}\times {\mathit{r}}_1}{{\Vert \mathit{L}\times {\mathit{r}}_1\Vert }_2^2+{{\displaystyle \left(\delta {\Vert \mathit{L}\Vert }_2\right)}}^2}{\displaystyle \left[\mathit{L}\cdot {\displaystyle \left({\widehat{\mathit{e}}}_1-{\widehat{\mathit{e}}}_2\right)}\right]}$ E2

Here, ${\mathit{r}}_1$ and ${\mathit{r}}_2$ are the position vectors of the point of interest relative to the ends of the vortex segment; ${\widehat{\mathit{e}}}_1$ and ${\widehat{\mathit{e}}}_2$ are the unit vectors parallel to ${\mathit{r}}_1$ and ${\mathit{r}}_2$ , respectively, and $\mathit{L}={\mathit{r}}_{\text{1}}-{\mathit{r}}_{\text{2}}$ is the vector along the vortex segment. The term $\delta {\Vert \mathit{L}\Vert }_2$ in Eq. (2) is added in order to avoid the singularity that appears when the point approaches the vortex line or its extension. The influence of the cutoff radius $\delta$ on the velocity is strongly felt in the immediate vicinity of the vortex line but is hardly noticeable elsewhere. A linear cutoff radius is another option [40].

3.2. Unsteady vortex‐lattice method

In UVLM, the vortex sheets are replaced by lattices of short straight vortex segments with spatially constant circulation. These segments divide the wing surface into a finite number of typically nonplanar, quadrilateral elements of area with straight edges often called panels. The model is completed by joining the free‐vortex lattices (wakes) to the bound‐vortex lattice (lifting surface) along the edges where separation occurs. The separation location is user supplied and typically based on experience. Each quadrilateral element of the bound lattice has a single unknown circulation $G{\displaystyle \left(t\right)}$ instead of the four unknown circulations around each of the short, vortex‐line segments along its edges. Consequently, the requirement of spatial conservation of circulation [41] is automatically satisfied throughout the lattices. Once all the $G{\displaystyle \left(t\right)}\text{'}\text{s}$ are known, the $\Gamma {\displaystyle \left(t\right)}\text{'}\text{s}$ around all of the straight vortex segments can be easily determined. In Figure 2, a representative example of the resulting bound and free‐vortex sheets and vortex lattices after the discretization is shown.

Figure 2.

Vortex sheets and vortex lattices.

The governing equation is complemented with the following boundary conditions: (i) the velocity field associated with the disturbance decays away from the body and its wakes. Hence,

$\underset{{\Vert r\Vert }_2\to \infty }{\lim }\Vert {\mathit{V}}_B{\displaystyle \left(\mathit{r};t\right)}+{\mathit{V}}_W{\displaystyle \left(\mathit{r};t\right)}+{\mathit{V}}_{LW}{\displaystyle \left(\mathit{r};t\right)}\Vert =\Vert {\mathit{V}}_{\infty }\Vert ,$ E3

where ${\mathit{V}}_B{\displaystyle \left(\mathit{r};t\right)}$ , ${\mathit{V}}_W{\displaystyle \left(\mathit{r};t\right)}$ , and ${\mathit{V}}_{LW}{\displaystyle \left(\mathit{r};t\right)}$ are the velocity fields associated with the bound‐vortex lattice, the free‐vortex lattice shed from the wing’s trailing edge (including the tip), and leading edge, respectively, and ${\mathit{V}}_{\infty }$ is the free‐stream velocity. The velocity field obtained from the Biot‐Savart law satisfies this condition. (ii) The normal component of the fluid velocity relative to the body’s surface must be zero:

${\displaystyle \left({\mathit{V}}_{\infty }+{\mathit{V}}_B{\displaystyle \left(\mathit{r};t\right)}+{\mathit{V}}_W{\displaystyle \left(\mathit{r};t\right)}+{\mathit{V}}_{LW}{\displaystyle \left(\mathit{r};t\right)}-{\mathit{V}}_P\right)}\cdot \widehat{\mathit{n}}=0,$ E4

where ${\mathit{V}}_P$ is the velocity of the body’s surface, and $\widehat{\mathit{n}}$ is a unit vector normal to the surface. Since there is a finite number of panels, and hence unknown circulations, Eq. (4) is only imposed at one point in each panel. Here, ${\mathit{V}}_B$ is expressed in terms of the unknown values for the loop‐vortex circulations $G_j{\displaystyle \left(t\right)}$ and the aerodynamic influence coefficients $A_{ij}{\displaystyle \left(t\right)}$ [35, 39]. To satisfy the unsteady Kutta condition at each time step, the vortex rings along the edges are shed into the flow where they have the same order as they had on the wing’s surface. The vortex rings are moved downstream with the flow by moving the end points of their vortex segments, called nodes, with the local fluid‐particle velocity $\mathit{V}$ to new positions, denoted $\mathit{r}{\displaystyle \left(t+\Delta t\right)}$ , according to the first‐order approximation:

$\mathit{r}{\displaystyle \left(t+\mathrm{\Delta }t\right)}=\mathit{r}{\displaystyle \left(t\right)}+\mathit{V}(\mathit{r},t)\mathrm{\Delta }t.$ E5

3.3. Leading‐edge separation model

Vortex shedding from the leading edge (LEV) depends on the angle between the local fluid velocity and the velocity of the edge of the moving wing’s plane (the effective angle of attack). Several studies on conventional aircraft wings have reported that flow attached to the wing starts to separate when the effective angle of attack exceeds a critical value of 12–15°. Several numerical tools based on vortex‐lattice methods account for LEV on highly swept delta wings [42]. Here, the effects of leading‐edge separation are included through an on/off mechanism that consists mainly of computing the value of the effective angle of attack ${\alpha }_e$ at each time step and comparing it with a reference value ${\alpha }_c$ . Only when ${\alpha }_e\ge {\alpha }_c$ is leading‐edge separation included [38].

3.4. Aerodynamic loads

The force on each element of the bound lattice is determined based on the pressure jump across the lifting surface at the control point; this calculation is carried out by using the unsteady Bernoulli’s equation:

${\partial }_t\varphi (\mathit{r},t)+\frac{1}{2}\mathit{V}(\mathit{r};t)\cdot \mathit{V}(\mathit{r};t)+\frac{p(\mathit{r},t)}{\rho }=H{\displaystyle \left(t\right)}$ E6

Here, ${\partial }_t$ denotes the partial time derivative at a fixed location in an inertial reference frame, ρ is the fluid density, p is the pressure, and H(t) is the total energy per unit mass, which only depends on time and has the same value at every point in the whole domain of the flow [35].

4. Structural model

The equations governing the structure were developed for large rotations and displacements, called primary motions, and small/moderate rotations and displacements with respect to a moving reference frame, called secondary motions. The primary motions describe the position and orientation of each body as a whole, and the secondary motions describe deformations. The Floating Frame of Reference (FFR) formalism and Lagrange’s method for constrained systems [43] were used to derive the equations of motion. The flying robot system is modeled as a collection of ${\mathit{n}}_b$ rigid and deformable bodies subjected to ${\mathit{n}}_c$ constraints, which are used to impose the connections and predefined motions. The fuselage (called the central body) is rigid, and the wings attached to it are assumed to be flexible. Three reference‐point coordinates that typically coincide with either the body’s center of mass or its centroid represent the rigid‐body’s position. The rigid body’s orientation is described by Euler angles. In what follows, bold italic letters and bold letters are used for tensor notation and matrix notation, respectively.

4.1. Reference frames

In the FFR formalism, the rigid‐body motion of the kth body is the primary motion, and the motion of points in the body’s reference frame is the secondary motion. If the displacement field contains rigid‐body modes, a set of conditions has to be imposed to define a unique displacement field with respect to the selected body reference. Consequently, $n_b+1$ reference frames are used in order to describe the kinematics of the multibody system: (i) an inertial system $N={\displaystyle \left\{{\widehat{\mathbf{n}}}_1,{\widehat{\mathbf{n}}}_2,{\widehat{\mathbf{n}}}_3\right\}}$ and (ii) a reference system fixed to each body k, $B^k={\displaystyle \left\{{\widehat{\mathbf{b}}}_1^k,{\widehat{\mathbf{b}}}_2^k,{\widehat{\mathbf{b}}}_3^k\right\}}$ . The set of vectors ${\widehat{\mathbf{n}}}_i$ and ${\widehat{\mathbf{b}}}_i^k$ for i = 1, 2, 3 form a dextral orthonormal basis.

4.2. Velocity for an arbitrary point in a body

Let us consider the kth flexible body, which moves and deforms in space. The components of the absolute translational velocity at an arbitrary point in the body P, ${\dot{\mathit{R}}}_p^k,$ , can be expressed in the segregated form:

${\dot{\mathit{R}}}_p^k={\dot{\mathit{R}}}^k+{\mathit{Q}}_{N{B}^k}{\displaystyle \left[{\dot{\mathit{r}}}_0^k+{\dot{\mathit{u}}}^k+{\mathit{\omega }}^k\times {\displaystyle \left(r_0^k+{\mathit{u}}^k\right)}\right]},\text{ for }k=1,\mathrm{\dots },n_b$ E7

Here, ${\dot{R}}^k$ is the velocity of the origin of the body‐fixed frame $B^k$ ; $r_0^k$ is the position vector of an arbitrary point P on the kth body in the undeformed configuration, and ${\dot{\mathit{r}}}_0^k$ is identically zero because the differentiation is performed with respect to the ${\mathit{B}}^k$ frame; ${\mathit{u}}^k$ is the elastic displacement vector of the point P and ${\dot{\mathit{u}}}^k$ represents its time derivative with respect to the body‐fixed frame; ${\mathbf{\omega }}^k$ is the angular velocity vector of $B^k$ relative to the inertial frame N; and ${\mathit{Q}}_{N{B}^k}:B^k\to N$ represents a rotation tensor from the body‐fixed frame $B^k$ to the inertial frame N. After algebraic manipulations, Eq. (7) can be written as

${\dot{\mathit{R}}}_p^k={\dot{\mathit{R}}}^k+{\mathit{Q}}_{N{B}^k}{\dot{\mathit{u}}}^k-Q_{N{B}^k}{\displaystyle \left[{\displaystyle \left({\overline{\mathit{r}}}_0^k+{\overline{\mathit{u}}}^k\right)}{\mathit{\omega }}^k\right]},\text{ for }k=1,\mathrm{\dots },n_b,$ E8

where ${\overline{\mathit{r}}}_0^k$ and ${\overline{\mathit{u}}}^k$ are skew‐symmetric tensors associated with ${\mathit{r}}_0^k$ and ${\mathit{u}}^k$ , respectively. Their actions on vectors are equivalent to the cross products. The tensor ${\mathit{Q}}_{N{B}^k}$ is a function of $q^k$ , the vector of absolute coordinates that describes the motion of the Floating Frame of Reference $B^k$ attached to body k, or primary motions. The vector $q^k$ is given by

$q^k={{\displaystyle \left[{{\displaystyle \left(R^k\right)}}^T{{\displaystyle \left({\theta }^k\right)}}^T\right]}}^T,$ E9

where $R^k$ is a set of coordinates for the origin of the body reference frame $B^k$ , and ${\theta }^k$ is a set of rotational parameters that orient the body reference frame (e.g., Euler angles, Rodriguez parameters, etc.). When high‐order terms are not considered, the elastic displacement vector solved in the body‐fixed frame is ${\mathbf{u}}^k=N^k{p}^k$ , where $N^k$ is a matrix of shape functions for interpolating the elastic displacement across the body domain, and $p^k$ is a time‐dependent vector of generalized coordinates that is used to describe secondary motions.

4.3. Equations of motion

Once the absolute vector of translational velocity for any point on the body is known, the equations of motion can be derived. Additionally, it is important to consider that a flexible body can be constrained [43, 44], by either linking it with other bodies or by choosing the parameters to describe primary motions; only holonomic constraints are considered. To simplify the manipulation, Lagrange's equations are separated into two groups: one corresponding to $q^k$ and the second to $p^k$ . The constrained Lagrange’s equations are given by

$d_t{\displaystyle \left({\partial }_{{\dot{\mathbf{q}}}^k}{T}^k\right)}-{\partial }_{q^k}{T}^k+{\mathit{B}}_{q^k}^T{\mathit{\lambda }}^k={{\displaystyle \left({\mathit{Q}}_q^k\right)}}^T,\text{ and }{d}_t{\displaystyle \left({\partial }_{{\dot{\mathbf{p}}}^k}{T}^k\right)}-{\partial }_{p^k}{T}^k+{\mathit{B}}_{p^k}^T{\mathit{\lambda }}^k+{\partial }_{p^k}{U}^k={{\displaystyle \left({\mathit{Q}}_p^k\right)}}^T,$ E10

which are complemented by a set of algebraic‐constraint equations expressed as

${\mathbf{\Phi }}^k{\displaystyle \left(q^{k-1},q^k,q^{k+1};t\right)}=0.$ E11

Here, the $q^k{\displaystyle \left(p^k\right)}$ were defined previously, and $d_t,{\partial }_{q^k}{\displaystyle \left({\partial }_{p^k}\right)}$ as well as ${\partial }_{{\dot{\mathbf{q}}}^k}{\displaystyle \left({\partial }_{{\dot{\mathbf{p}}}^k}\right)}$ denote time derivatives, the partial derivatives with respect to $q^k{\displaystyle \left(p^k\right)}$ and with respect to ${\dot{\mathbf{q}}}^k{\displaystyle \left({\dot{\mathbf{p}}}^k\right)}$ , respectively. ${\mathbf{\lambda }}^k$ is the vector of Lagrange multipliers, ${\mathit{B}}_{q^k}{\displaystyle \left({\mathit{B}}_{p^k}\right)}$ is the Jacobian tensor of constraints associated with the coordinates $q^k{\displaystyle \left(p^k\right)}$ , $T^k$ is the kinetic energy of the k‐th body, $U^k$ is the elastic potential energy, and ${\mathbf{Q}}_q^k{\displaystyle \left({\mathbf{Q}}_p^k\right)}$ is the generalized load vector corresponding to the vector $q^k{\displaystyle \left(p^k\right)}$ . The kinetic energy $T^k$ follows from ${\dot{\mathit{R}}}_p^k$ in Eq. (8), and then the equations for primary and secondary motions for the k‐th body are as follows:

$\begin{array}{c}{M}^k{\ddot{q}}^k+{\left({\mathfrak{M}}^k\right)}^T{\ddot{p}}^k+B_{q^k}^T{\mathit{\lambda }}_k={\left({\mathbf{Q}}_q^k\right)}^T,\\ m^k{\ddot{p}}^k+{\mathfrak{M}}^k{\ddot{q}}^k+K^k{p}^k+B_{p^k}^T{\mathit{\lambda }}_k={\left({\mathbf{Q}}_p^k\right)}^T,\\ {\mathbf{\Phi }}^k\left(q^{k-1},q^k,q^{k+1},p^{k-1},p^k,p^{k+1};t\right)=0,\text{ for }k=1,\mathrm{\dots },n_b.\end{array}$ E12

Here, $M^k$ is the mass matrix for primary motions, which is differentiable, symmetric, and at least positive‐semi‐definite; ${\mathfrak{M}}^k$ is the mass matrix that couples primary and secondary motions; $m^k$ is the metric tensor for secondary motions, which is constant, symmetric, positive‐definite, and its matrix representation is the “elastic mass matrix”; $K^k$ is the stiffness matrix for the secondary motions; and ${\Phi }^k$ is the set of holonomic‐rheonomic constraint equations associated with the k‐th body. It should be noted that $q^T={\displaystyle \left\{{{\displaystyle \left(q^1\right)}}^T,\mathrm{\dots },{{\displaystyle \left(q^m\right)}}^T\right\}}$ for $m\le n_b$ and $p^T={\displaystyle \left\{{{\displaystyle \left(p^1\right)}}^T,\mathrm{\dots },{{\displaystyle \left(p^n\right)}}^T\right\}}$ for $n\le n_b.$ . The total number of coordinates is calculated through $6n_b+{\displaystyle {\sum }_{k=1}^{n_b}{N}_k}=n_{\text{coord}}$ (N_k is the number of elastic generalized coordinates of the k‐th body). Thus, the number of degrees of freedom is $n_{\text{dof}}=n_{\text{coord}}-n_c$ . A beam element based on Euler‐Bernoulli/Rayleigh theory is used. Hermite polynomials are used to interpolate displacement/rotation fields in each finite element from nodal values. First‐order polynomials are used to interpolate elongation and torsion, and third‐order polynomials are used to interpolate bending [45]. Finally, the equations of motion for the complete multibody system are obtained by assembling the equations of motion for each body:

$\begin{array}{l}\mathbf{M}\ddot{x}+B_x^T\mathit{\lambda }=F,\\ \mathbf{\Phi }=0.\end{array}$ E13

Here, $M\in {\mathbb{R}}^{n_{\text{coord}}\times n_{\text{coord}}}$ is the global mass matrix, $B_x\in {\mathbb{R}}^{n_{\text{coord}}\times n_c}$ the constraint Jacobian matrix, $x\in {\mathbb{R}}^{n_{\text{coord}}\times 1}$ is the global vector of generalized coordinates, $\lambda \in {\mathbb{R}}^{n_c\times 1}$ is the global vector of Lagrange’s multipliers, $F\in {\mathbb{R}}^{n_{\text{coord}}\times 1}$ is the global vector of forces, and $\Phi \in {\mathbb{R}}^{n_c\times 1}$ is the set of all constraints for the multibody system.

5. Communication between models and numerical integration

In this section, the schemes for exchanging information (loads, displacements, and velocities) between the structural and aerodynamic models and integrating the resulting algebraic‐differential equation are described.

5.1. Information transfer between simulators for aerodynamics and structural dynamics

At each time step, the aerodynamic and structural models exchange information bi‐directionally in an iterative sequence to improve the estimates of both the structure’s response and the aerodynamic loads (see Figure 3). To this end, the authors propose the use of an interpolation procedure based on radial basis functions (RBFs) [46].

Figure 3.

Schematic diagram of transfer technique between simulators.

A linear relationship between the generalized displacements associated with both the aerodynamic mesh and the structural mesh is established through the following coupling matrix:

$w=G_{AE}^{NA}u,\text{ and}v=G_{AE}^{CP}u$ E14

Here, w and v are the displacements of the nodes and control points of the aerodynamic mesh, and u is the generalized nodal displacements of the structural mesh. $G_{AE}^{NA}$ is the coupling matrix that relates w with u, and $G_{AE}^{CP}$ relates the control‐point and the structural‐node displacements. As a first step, a continuous scalar function $s{\displaystyle \left(x\right)}$ is generated by using the radial basis functions $\phi {\displaystyle \left(\text{|}x-x_i^e\text{|}\right)}$ [47, 48], where x is the location for the evaluation of s, and $x_i^e$ is the positions of centers (structural nodes in the present case). It is required that the function s recovers the generalized structural displacements when it is evaluated at the centers. For the component z of u (denoted by $u_i^z$ ),

$u_i^z=s{\displaystyle \left(x_i^e\right)}={\alpha }_1\phi {\displaystyle \left(\text{|}{x}_i^e-x_1^e\text{|}\right)}+\cdots +{\alpha }_i\phi {\displaystyle \left(\text{|}{x}_i^e-x_i^e\text{|}\right)}+\cdots +{\alpha }_{ne}\phi {\displaystyle \left(\text{|}{x}_i^e-x_{ne}^e\text{|}\right)}fori=1,2,\mathrm{\dots },n_e.$ E15

This constitutes a set of n_e (number of structural nodes) linear algebraic equations, which have the matrix form $A_{ee}\alpha =u^z$ , where $A_{ee}\in {\mathbb{R}}^{n_e\times n_e}$ contains the functions ϕ evaluated at the structural nodal positions, $\alpha$ is a column vector that contains the coefficients α_i, and u^z is an array which collects the z‐components associated with the nodal translations of the structural mesh.

Solving for the unknown $\alpha$ , makes it possible to approximate the aerodynamic nodal displacements by evaluating the function s(x) at the aerodynamic nodal positions $x_i^a$ . Again, this procedure leads to a linear algebraic system $A_{ae}\alpha =w^z$ , where $A_{ae}\in {\mathbb{R}}^{n_a\times n_e}$ (n_a is the number of aerodynamic nodes) contains the functions φ evaluated at the aerodynamic nodal positions. Combining the aforementioned two linear systems, a direct relation between w^z and u^z is obtained, $w^z=H{u}^z$ , where $H=A_{ae}{A}_{ee}^{-1}$ . If the same RBFs are used to interpolate the x(w^x) and y(w^y) components of w, a similar relation is found. Therefore,

$w^r=H{u}^r,\text{ for }r=x,y,z.$ E16

In fact, matrix H is a linear map that relates the aerodynamic and structural nodal displacements. According to prior studies [35, 49], the procedure to obtain each one of the matrices introduced above consists of the following: (i) computing H based on the data associated with the structural and aerodynamic meshes, (ii) obtaining $G_{AE}^{NA}$ from H, and (iii) computing $G_{AE}^{CP}$ [49]. The last matrix is used to transform the aerodynamic loads $F_A$ , acting at the control points of the aerodynamic grid, to an equivalent set of forces $F_E$ acting on the finite‐element mesh, $F_E={{\displaystyle \left(G_{AE}^{CP}\right)}}^T{F}_A$ . The structural‐aerodynamic force relation is obtained by establishing some kind of abstract equivalence; that is, both systems of forces must perform the same virtual work for any given virtual displacement [35]. In particular, this procedure has been developed for combining UVLM grids with structural meshes based on beam finite elements. As the coupling matrix only depends on the kind of RBFs used and the clouds of points to be interpolated, the same formulation is valid for other types of finite‐element meshes.

5.2. Numerical integration of the index‐1 DAEs

The set of equations of motion for the entire multibody system represents an index‐3 system of DAEs. These dynamic equations are, in general, nonlinear. To solve them by means of standard solvers for ordinary differential equations (ODEs), an index reduction is required for the set of DAEs [44]. The methodology adopted in this work includes differentiation of the constraint equations twice with respect to time. This new set of equations is often called constraint acceleration level and is given by

$\ddot{\mathbf{\Phi }}{\displaystyle \left(x,\dot{\mathbf{x}};t\right)}=B_x\ddot{\mathbf{x}}+2{\partial }_x{\displaystyle \left({\partial }_t\mathbf{\Phi }\right)}\dot{\mathbf{x}}+{\partial }_{tt}\Phi +{\partial }_x{\displaystyle \left(B_x\dot{\mathbf{x}}\right)}\dot{\mathbf{x}}=0,$ E17

where $x^T={\displaystyle \left\{q^T,p^T\right\}}$ , ${\partial }_x$ denotes the partial derivative with respect to the vector $x$ , and ${\partial }_{tt}$ denotes the second partial time derivative. The equations of motion for the multibody system (13) can be rewritten together with the acceleration level constraint equations, as an index‐1 system of DAEs as follows:

${\displaystyle \left[\begin{array}{cc}M& B_x^T\\ B_x& 0\end{array}\right]}{\displaystyle \left\{\begin{array}{c}\ddot{\mathbf{x}}\\ \nu \end{array}\right\}}={\displaystyle \left\{\begin{array}{c}Q\\ \kappa \end{array}\right\}},$ E18

where $\kappa =-2{\partial }_x{\displaystyle \left({\partial }_t\Phi \right)}\dot{\mathbf{x}}-{\partial }_{tt}\Phi -{\partial }_x{\displaystyle \left(B_x\dot{\mathbf{x}}\right)}\dot{\mathbf{x}}$ .

Numerical integration of Eq. (18) is susceptible to instabilities as a consequence of truncation procedures and round‐off errors. The evident one is that the position and velocity constraints are no longer exactly satisfied; that is, a drift‐out of the constraints does exist. In the literature, several stabilization methods to correct this numerical drift can be found, among which the most widely used due to its simplicity is Baumgarte’s technique [50]. Another technique currently used consists in discretizing numerically the ODE and projecting the approximate solution onto the selected constraints manifold (coordinate projection) [51]. In this work, the authors have adopted the coordinate projection method to control/eliminate the numerical drift. As mentioned in Section 2, the numerical scheme used by Simulator 2 is well known and can be found in prior work [35, 39]. On the other hand, the numerical procedure adopted for Simulator 1 to solve the equations of motion of the aerial robot is based on Hamming’s fourth‐order predictor‐corrector method [35, 52].

5.2.1. Co‐simulation strategy

During a time step $\Delta t$ , the wakes are convected to their new positions with the local fluid‐particle velocity. Simultaneously, the structure of the aerial robot moves to its new position as a result of the acting forces and constraint equations. This concept is implemented by performing the following sequence of steps to calculate the solution at time $t+\Delta t$ as follows:

1. Simulator 2 is used to predict the new position of the wakes. A fluid particle in the wake moves from its current position $r{\displaystyle \left(t\right)}$

2. to its new position $r{\displaystyle \left(t+\Delta t\right)}$

3. according to Eq. (5). During the rest of the procedure for this time step, the wake is frozen.

4. The current loads computed by Simulator 2 are used in Simulator 1 to predict the response of the structure.

5. The current state of the structure is used as input to Simulator 2 and the loads are recalculated. Then, these loads are used as input to Simulator 1 and the state of the structure is updated again. This step is repeated until convergence. Usually, three to seven iterations are required to reduce the error to less than 10⁻¹⁰.

6. The final position and velocity of the structure are evaluated by using Simulator 1, and used by Simulator 2 to recalculate the flow field and obtain the final estimate for the aerodynamic loads.

The procedure described above needs information from four previous time steps. At the beginning of the procedure, this information does not exist; so, the authors have used a special starting scheme: at $t=0$ , the initial conditions are used by Simulator 2 to calculate the aerodynamic loads ignoring the contribution of ${\partial }_t\varphi$ . Further studies on this issue are needed to definitively state that ${\partial }_t\varphi$ can be neglected at $t=0$ .

6. Numerical results

In this section, results obtained with the implementation of the proposed methodology written in Fortran 90 are presented. The implementation is structured in a modular organization; that is, so that each part can be individually removed and replaced, as well as the capability of adding new models without modifying the general structure of the program. Although the parallelization of the code is not implemented, acceptable computing times, for full nonlinear simulations, were achieved. Simulations involving only structural models may take only a few minutes, but aeroelastic simulations including unsteady, nonlinear aerodynamics may take above 1 day on average in a desktop computer with an i7 processor, and RAM DDR3 of 4 GB.

6.1. Aerial vehicles with flapping wings

In this subsection, a series of numerical results obtained with the developed computational tool for a mechanical system based on flapping wings is presented. The first set of results consists of the aerodynamic analysis of MAVs with flapping wings, including some validations against experimental data. The second set of results is related to the study of the unsteady aerodynamics and nonlinear dynamics of MAVs, which explains some basic features of the involved flight mechanisms.

6.1.1. Aerodynamics of flapping wings

In this first study, the lift forces obtained from numerical simulations and comparisons of them with the experimental data reported by Dickinson et al. [1] are presented. The dynamically scaled experimental model imitated a Drosophila melanogaster, dubbed Robofly, with a wing span of 25 cm (from the force sensor to the wingtip). The wing executed an insect‐like flapping motion at a frequency of 0.145 Hz with the wing tip tracing out a figure of eight. The viscosity of the oil, the length of the wing, and the flapping frequency were chosen in order to match the Reynolds number (Re) typical of the flight of a fruit fly (Re = 136). The kinematic pattern employed by Dickinson’s team consisted of a stroke amplitude of 160°, and an angle of attack at midstroke of 40° for both upstroke and downstroke. In Figure 4, the authors show the three kinematic patterns used: (i) wing rotation precedes the reversal stroke by 8% of the wing‐beat cycle (advanced pattern), (ii) wing rotation occurs symmetrically with respect to the reversal stroke (symmetrical pattern), and (iii) wing rotation is delayed with respect to the stroke reversal by 8% of the stroke cycle (delayed pattern). Details of the kinematic model used in this work to prescribe the motion of the Robofly’s wing, as well as an extensive study of kinematical parameters can be found in reference [53].

Figure 4.

(a) Stroke position angle and angle of rotation of the wing for three different rotational timings. (b) Time derivatives of the stroke position angle and the rotation angle for three rotational timings. Here, T_f is the flapping period.

In Figure 5, numerical results for the lift force (advanced pattern) with and without leading‐edge separation and the comparison with those obtained by Dickinson et al. are shown. The results obtained are very encouraging because they show better agreement than those reported in previously published comparisons such as the CFD study by Sun and Tang [54], and the two‐dimensional aerodynamic model developed by Ansari et al. [55, 56]. These results are significant because they justify the use of the nonlinear unsteady vortex‐lattice method to study the three‐dimensional aerodynamics of insects executing different maneuvers. Similar results were found for the symmetrical and delayed patterns [38]. These results suggest that the present model will also serve well for the study of aerial robots.

The Fast Fourier Transforms (FFTs) of the force data from numerical simulations and experiments were also computed and the results are shown in Figure 5b. The frequency content for the experimental data clearly shows the flapping frequency of the motion (point A in Figure 5b, $n_f=0.1446\text{ Hz}$ ) and twice this frequency (point B in Figure 5b, $2n_f=0.2893\text{ Hz}$ ) along with a number of harmonics. These harmonics appear because there are two half‐strokes per wing‐beat cycle (the wing‐passing frequency). The frequency content of the predicted lift force closely matches the flapping frequency of the motion ( $n_f=0.145\text{ Hz}$ and $2n_f=0.29\text{ Hz}$ ). Another indication of the quality of the numerical prediction is given by the mean lift, $\overline{L}$ . The square and circular symbols in Figure 5b represent the experimental and predicted mean lift force, respectively. The difference between the predicted mean lift force and the experimental mean lift force for the advanced pattern is about 4.5%. For a detailed study of the aerodynamic of flapping wings using a modified version of the unsteady vortex‐lattice method, the reader is referred to the work of Roccia et al. [38].

Figure 5.

Comparison of numerical results with experimental data for the advanced pattern (first‐stroke cycle). Dotted line for Dickinson’s data, solid line for numerical results without LEV, and broken line for numerical results with LEV.

6.1.2. Dynamics of flapping wings

In this section, the numerical results obtained for the dynamic behavior of a fruit fly (D. melanogaster) in hovering flight are presented. The general equations of motion were adapted to consider only imposed deformations. These include torsion and bending in the directions normal and parallel to the wing’s chord. Details of the kinematics and the dynamic models that allow one to prescribe such patterns of deformation are given in reference [57, 58]. Data reported by Fry et al. [58] on the actual kinematics of a fruit fly in hover were used to describe the wing motion over a flapping cycle (see Figure 6a where the solid blue line represents the stroke deviation angle, the dotted red line represents the stroke position angle, and the broken black line represents the rotation angle). Due to the complex motion experienced by the wings during a stroke cycle, the wake shed from the leading edge during the downstroke is likely cut by the wing during the upstroke. A wake rupture model is not yet available, and therefore no leading‐edge separation was allowed in this case. For this numerical experiment, the following parameters were chosen: flapping frequency $n_f=210\text{ Hz}$ ; wing span $R=2.5\text{ mm}$ and wing area $S=2.21{\text{ mm}}^2$ ; air density ${\rho }_{air}=1.2{\text{ kg/m}}^3$ ; insect’s mass of 0.84 mg. A spatial discretization of the MAV of 1405 aerodynamic panels was used: 200 panels over each of the wings and 1005 panels over the central body (head, thorax, and abdomen of the insect). Only the first stroke cycle was simulated, which was discretized into 100 time steps. In Figure 6b, the authors show the computational model of the insect where the relative size and topology of the aerodynamic grid can be appreciated.

Figure 6.

(a) Actual kinematics of a fruit fly in hovering, circular markers indicate experimental data. (b) Computational model of the insect.

Two cases are presented next. In the first, the authors consider a rigid wing, and in the second, the authors consider a flexible wing with prescribed deformation. The prescribed deformation consists of wing torsion, which varies linearly along the wing span, combined with bending in the normal and tangential direction of the wing chord. The values used for the three different deformations are based on a detailed study carried out by the authors on the influence of spanwise twisting and bending on the lift‐force generation in flapping‐wing MAVs [58]. In order to show the role of the wing’s flexibility in the production of lift, the parameters that regulate the deformation were tuned in order to increase the lift force throughout the stroke cycle with respect to the rigid‐wing model. The set of initial conditions for the two cases studied consists of a body angle of 75° and a stroke plane angle of 15°. These values produce a horizontal stroke plane. Both the linear and angular velocities of the central body are zero at $t=t_0$ .

In Figure 7, the authors show the lift and horizontal forces obtained numerically for both the rigid‐wing model and the flexible‐wing model along with experimental measurements obtained by Fry et al. [58] for a Drosophila in hover. The most important requirement of hovering flight is undoubtedly that the value of the vertical force must compensate for the insect’s weight. In Figure 7a, a peak force can be observed at the center of each half‐stroke, mainly during the upstroke. The peak force produced by the flexible wing is approximately 30% higher than what the rigid‐wing model predicts. This fact is also reflected in Figure 7d, in which the degree of freedom associated with the vertical displacement of the center of mass of the central body has been plotted. When the authors use a rigid‐wing model, the lift force produced can hardly support the insect’s weight, showing a slight downward movement at the end of the stroke cycle. On the other hand, when the authors impose a suitable deformation pattern, not only the lift generated is sufficient to balance the insect’s weight, but it also rises almost continuously.

Figure 7.

Aerodynamic forces. (a) Lift force. (b) Horizontal force. (c) Longitudinal displacement. (d) Vertical displacement. (e) Pitch angle. (f) Roll, yaw, and lateral displacements.

In order to assess the quality of numerical simulations, the aerodynamic forces obtained numerically are compared with experimental measurements. As can be observed in Figure 7a, the time evolution of the lift force for the rigid and flexible wing shows a shape and trend similar to the curve reported by Fry et al. [58]. Nevertheless, lift levels predicted by the numerical model are lower than those measured experimentally (approximately 14% lower for flexible wing case). Possibly, this fact is because the LEV phenomenon was not taken into account. On the other hand, the predicted horizontal force (Figure 7b) departs quite noticeably from the results of Fry et al. However, the shape and trend is well captured. It should be noted that the aerodynamic model used is inviscid. Furthermore, UVLM does not account for leading‐edge suction. In addition, the horizontal force is always opposite to the wing’s motion, leading to a longitudinal oscillation of the insect’s body. This dynamic behavior can be observed in Figure 7c.

Unlike forward thrust and lift, sideways thrust cancels instantly due to the bilateral symmetry of the wing motion (see Figure 7f). Similarly, because each wing contributes to yaw and roll torque in opposite directions, these moments cancel out. On the other hand, the sign and magnitude of pitch torque is the same for both wings. This fact directly affects the pitch angle (see Figure 7e). Fry et al. [58] observed that pitch torque averages zero over a complete cycle, and therefore the insect’s body experiences an oscillation in the pitch angle. In real insects, this phenomenon is controlled by the wing kinematics, the flapping frequency, and the precise motion of the head and abdomen relative to the insect’s thorax. Such a control mechanism is essential to prevent the pitch torque from growing indefinitely. In this work, no control strategy was applied, and therefore the pitch angle does not behave as expected. It decreases continuously for most of the flapping cycle, only showing a slight increase toward the end of the upstroke.

6.2. Aerial vehicles based on the morphing‐wing concept

In this section, a series of numerical results for an aerial vehicle with morphing wings is presented. The first set describes the aerodynamic behavior of a bio‐inspired, morphing‐wing concept. The second set of results is related to the aeroelastic study of an aerial‐robot wing model whose planform geometry is based on the seagull wing.

6.2.1. Aerodynamics of morphing wings

The model’s geometry preserves certain morphological parameters of a gull wing [59], and a folding wing mechanism is used to change the wing’s shape. The model consists of the right and the left wings, which are joined at the wing root. Because the two parts are mirror images, the model is represented by only the right wing (see Figure 8a).

Figure 8.

(a) Wing’s geometry, (b) definition of the inner and outer wing, and (c) definition of the dihedral angles.

Each wing is represented by two parts: the inner wing and the outer wing. In Figure 8a, the boundary between these parts is pointed out. A seagull‐wing profile obtained by Liu et al. [60] is used. The center line, contained in the plane of symmetry xz, does not move during the morphing process. The inner and outer parts of the wing are connected and can move with respect to each other. The folding motion of the inner and outer wings is described by the dihedral angles ${\theta }_B$ and ${\theta }_C$ , respectively (see Figure 8b and c). During the simulation, the wings remain fully extended until a time T_I, and then the morphing process begins by varying the dihedral angles according to expression (19) until a time T_S, where the wings take a configuration that is maintained in a steady state and is determined by the angles ${\theta }_B^S$ and ${\theta }_C^S$ . The simulation ends at a time T_F.

${\theta }_B{\displaystyle \left(t\right)}={\theta }_B^S{\displaystyle \left(\frac{t-T_I}{\Delta T}-\frac{1}{2\pi }\sin {\displaystyle \left(\frac{2\pi {\displaystyle \left(t-T_I\right)}}{\Delta T}\right)}\right)}\text{and}{\theta }_C{\displaystyle \left(t\right)}={\theta }_C^S{\displaystyle \left(\frac{t-T_I}{\Delta T}-\frac{1}{2\pi }\sin {\displaystyle \left(\frac{2\pi {\displaystyle \left(t-T_I\right)}}{\Delta T}\right)}\right)}.$ E19

Here $\Delta T=T_S-T_I$ . Three passages through the morphing processes were simulated with a mesh of 8 × 32 panels in the chordwise and spanwise directions, respectively. The morphing wing was reconfigured in time by means of Eq. (19) with the values given in Table 1. In each case, two conditions were considered: (i) with shedding vorticity from the wing tips and (ii) without shedding vorticity from the wing tips.

Case	${\theta }_B^S$	${\theta }_C^S$	TI	TS	TF	V8	$\alpha$
I	−30°	30°	35 c0/V8	70 c0/V8	90 c0/V8	12.4 [m/s]	5°
II	10°	30°
III	−10°	−30°

Table 1.

Simulation parameters.

In Figure 9, the lift coefficients for the three cases are plotted as functions of time. There are only slight differences between the solutions obtained with conditions (i) and (ii); thus, the vorticity being shed from the wingtips may be ignored. For all three cases, the number of elements in the wake associated with the wingtips is 20% of the total in the wake; when the wing‐tip elements are eliminated the computational time is reduced by 35%.

Figure 9.

Lift coefficient as a function of time.

During the first part of the flight simulation, when the wing is fully extended, the lift rapidly increases until the steady state is reached at T_I. During the second part of the flight T_I < t < T_S, which is the morphing process, the lift depends on the prescribed kinematics of the dihedral angles. In the third phase of the process after the morphing is completed, t > T_S, the wings approach new steady states. In Case III, the lift smoothly decreases until it reaches a local minimum, t ≈ T_M = (T_I + T_s)/2, and then it increases to a value that corresponds to the wing’s configuration in the new steady state. In this case, it appears that during the morphing process there could be a loss in altitude or perhaps a decrease in the rate of climb or an increase in the rate of descent. In Case II, the lift increases during the morphing process to reach a local maximum, t ≈ T_M, and then it begins to decrease until it reaches a value that corresponds to the wing’s new configuration in the steady state. In this case, it appears that the morphing process would produce a gain in altitude. In Case I, the lift behavior is different; it decreases smoothly until the steady‐state solution is reached. For this last case, the lift characteristic does not have either a local maximum or a local minimum during the morphing. With some tweaking in the three cases, the wing could arrive at any (within reason) prescribed steady‐state lift, but most likely with the aircraft at different altitudes and different airspeeds. Here, although the airspeed was held constant in these three examples, it could be changed step by step. The present model provides an accurate tool to account for the following: (i) the strong dependence between the unsteady aerodynamic loads and the morphing motion, and (ii) the power needed to drive the morphing [61].

6.2.2. Aeroelasticity of morphing wings

The wing model adopted for the numerical aeroelastic simulations presented in this section consists of an inner and outer wing as it was described in Section 6.2.1. The wing’s structural skeleton consists of a set of structural elements modeled as beams with rectangular cross sections (width = 2 cm and height = 0.6 cm) and material properties: density ρ_b = 1187 kg/m³, Young’s modulus E = 3.18 GPa, and shear modulus G = 1.35 GPa (see Figure 10).

Figure 10.

Computational wing model.

The inner and outer wing sections are oriented in space through the dihedral angles ${\theta }_B$ and ${\theta }_C$ . However, the wings do not have a prescribed motion (pitching, flapping, and lagging) and the wing’s roots are considered to be clamped to the fuselage. This dynamic configuration deactivates the large rotations and displacements associated with primary motions; therefore, the set of governing equations (12) is reduced to the well‐known set of classical ODEs for structural dynamics,

$m\ddot{\mathbf{p}}+Kp={{\displaystyle \left({\mathit{Q}}_p\right)}}^T.$ E20

After algebraic manipulations and after expressing secondary motions as a linear combination of the free-vibration modes of the structure, the set of Eq. (20) is rewritten in a dimensionless form

$d_{\widehat{t}}{}_{\widehat{t}}\widehat{\mathbf{p}}{\displaystyle \left(\widehat{t}\right)}+\widehat{\Lambda }\widehat{\mathbf{p}}{\displaystyle \left(\widehat{t}\right)}={\displaystyle \left(\frac{1}{2}{\rho }_C{L}_C^4\right)}{\tilde{\mathbf{m}}}^{-1}{\mathbf{\Psi }}^T{{\displaystyle \left(G_{AS}^{CP}\right)}}^T{\widehat{\mathbf{F}}}_A,$ E21

where $\widehat{t}$ is the Dimensionless-time (t/T_C); $d_{\widehat{t}\widehat{t}}$ denotes the second dimensionless time derivative; $\widehat{\mathbf{p}}$ is a column vector containing the modal coordinates associated with the free vibration modes; $\tilde{\mathbf{m}}$ is the diagonal modal mass matrix; $\widehat{\mathbf{\Lambda }}$ is a diagonal matrix that contains the natural frequencies in a dimensionless form (reduced frequency); $\Psi$ is the modal matrix; $L_C$ , $T_C$ and ${\rho }_C$ are characteristic properties of the problem [35]; and ${\widehat{\mathbf{F}}}_A$ is the dimensionless form of the aerodynamic loads [49]. The rest of the variables are as defined before.

The main goal of the numerical experiment is to compute the flutter velocity $V_F$ for different dihedral angles, ranging from ${\theta }_B=0°$ to ${\theta }_B=45°$ , with steps of $\Delta {\theta }_B=5°$ . ${\theta }_C=0°$ was used for all values of ${\theta }_B$ . Due to the lack of primary motions, the dihedral angles remain constant during the whole simulation. Additionally, the angle of attack is $\alpha =0°$ for all simulations. The computational model consists of the following: (i) a finite‐element mesh for the structural discretization consisting of 21 beam elements and (ii) a vortex‐lattice mesh for the aerodynamic surface consisting of 256 panels (8 in the chordwise direction and 32 in the spanwise directions). In Figure 11, the authors show the flutter speed versus the dihedral angle of the inner wing ${\theta }_B$ . It can be observed that V_F decreases linearly from the maximum value of 29.66 m/s (for ${\theta }_B=0°$ ), between 10° and 30°, to the minimum of 18 m/s for ${\theta }_B=45°$ . The curve represents a boundary for the stable and unstable regions. All flight conditions located below the curve lead to a response that tends to a static equilibrium position for any perturbation. On the contrary, the response for flight conditions located above the curve achieves a periodic motion with an amplitude that depends on the flight speed, that could lead to a collapse of the structure.

Figure 11.

Flutter speed versus dihedral angle of the inner wing.

In Figure 12a, the amplitudes of three vibration modes are shown as functions of dimensionless time for a flight condition (pointed out in Figure 11 by a) located in the stable region (free stream of 21 m/s and a configuration defined by ${\theta }_B=20°$ ). An initial perturbation on the second mode ${\widehat{p}}_2{\displaystyle \left(0\right)}=0.08$ was imposed in the simulation. The amplitudes of all vibration modes clearly decay, not to zero, but to their equilibrium positions. This is exclusively due to the aerodynamic damping, because the structural model does not include any damping. In Figure 12b, the authors show the phase portraits for two modes. Both modes have a behavior similar to a stable focus.

Figure 12.

(a) Aeroelastic response for three modes as function of the dimensionless time. (b) Phase portrait for two modes (V_∞ = 21.5 m/s and θ_B = 20°).

In Figure 13a, the amplitudes of the same three modes are also shown as functions of time but for a flight condition located in the unstable region (pointed out in Figure 11 by b, with a free stream of 24.5 m/s and a configuration defined by ${\theta }_B=20°$ ). Two initial conditions were considered in the simulation. The corresponding responses (blue-dashed line and red-solid line) are characterized by a transient behavior that evolves into periodic orbits whose amplitudes are the same. In Figure 13b, the phase portraits for two of the modes clearly reveal the presence of a limit-cycle oscillation (LCO). The trajectory associated with the first initial condition approaches the LCO from outside, while the trajectory associated with the other initial condition approaches the LCO from inside.

Figure 13.

(a) Aeroelastic response for three modes as function of dimensionless time. (b) Phase portrait for two modes (V_∞ = 24.5 m/s and θ_B = 20° for both responses, blue-dashed line and red-solid line).

7. Concluding remarks

A combined structural and aerodynamic model for studying the aeroelastic behavior of flying robots with flexible wings has been presented, and aeroelasticity studies have been carried out for systems with flexible wings. The modular approach is based on the following: (i) a general, nonlinear‐unsteady vortex‐lattice method, (ii) a segregated version of Lagrange’s equations and the Floating Frame of Reference formalism for constrained systems, (iii) tight coupling of the aerodynamic model (UVLM) with the structural model, (iv) radial basis functions to transfer information between the aerodynamic and structural grids, and (v) an efficient and novel algorithm, based on a computational platform with modular structure, to integrate the equations interactively and simultaneously in the time domain.

As illustrated by the examples, the authors’ methodology will enable fully coupled aeroelastic simulations for MAVs with flapping wings and/or morphing wings and further investigations in this direction. Their development paves the way to advance the current methodology by implementing a vortex‐particle method to improve the description of the wakes and a fast multi‐pole method for computing the velocity field, and by including various other beam and shell elements.

Acknowledgments

The authors gratefully acknowledge the partial support received from the Consejo Nacional de Investigaciones Científicas y Técnicas, Argentina, the U.S. National Science Foundation through Grant No. CMMI‐1250187, the U.S. Air Force Office of Scientific Research through Grant No. FA95501510134, and the Minta Martin Foundation. In addition, the authors would like to thank the Grupo de Electrónica Aplicada (GEA) and Grupo de Matemática Aplicada (GMA), Engineering School, Universidad Nacional de Río Cuarto, Argentina.