A method for validating rigid-body models of compliant mechanisms under dynamic loading conditions using motion tracking cameras and genetic algorithms is presented. The compliant mechanisms are modeled using rigid-body mechanics as compliant joints (CJ): spherical joints with distributed mass and three-axis torsional spring dampers. This allows compliant mechanisms to be modeled using computationally efficient rigid-body dynamics methods, thereby allowing a model to determine the desired stiffness and location characteristics of compliant mechanisms spatially distributed into a structure. An experiment was performed to validate a previously developed numerical dynamics model with the goal of tuning unknown model parameters to match the flapping kinematics of the leading edge spar of an ornithopter with contact-aided compliant mechanisms (CCMs), compliant mechanisms that feature self-contact to produce nonlinear stiffness, inserted. A system of computer motion tracking cameras was used to record the kinematics of reflective tape and markers placed along the leading edge spar with and without CCMs inserted. A genetic algorithm was used to minimize the error between the model and experimental marker kinematics. The model was able to match the kinematics of all markers along the spars with a root-mean-square error (RMSE) of less than 2% of the half wingspan over the flapping cycle. Additionally, the model was able to capture the deflection amplitude and harmonics of the CCMs with a RMSE of less than 2 deg over the flapping cycle.

## Introduction

In this paper, we present a method for validating a rigid-body mechanics model of contact-aided compliant mechanisms (CCMs) in a dynamic structure. A numerical rigid-body dynamics model of a flapping wing structure called the dynamic spar numerical model (DSNM) was previously developed as an efficient method for modeling compliant mechanism dynamics using rigid-body techniques in an ornithopter wing structure [1,2]. This model uses kinematic constraints to formulate the equations of motion for the wing structure of an avian-scale ornithopter. The focus of the current paper is a method developed to tune model parameters using a benchtop experiment and test the assumptions made during the development of the compliant mechanism component of the DSNM.

Compliant mechanisms are flexible mechanisms designed to allow localized large deformations, typically by using thin, flexible members. Unlike conventional rigid-link mechanisms, compliant mechanisms have no relative moving parts, do not require lubrication, and can be monolithic devices that require no assembly. CCMs are a class of compliant mechanisms where the compliant members are designed to come in contact with one another to perform a specific task or to improve the performance of the mechanism itself. These mechanisms can vary from a simple single-point contact to multipoint curved contact surfaces in separate parts of the mechanism [3–6]. Beyond gripping, CCMs can have nonlinear stiffness and stress relief capabilities, and be designed for path generation [3,7–15]. Several examples of ongoing efforts to model compliant mechanism deformation and kinematics using the finite element method can be found in the literature [16–18]. These models include tensegrity mechanisms to alleviate stress, distributed-compliance grippers, and translational joints. The pseudo-rigid body model was developed to model compliant mechanisms as rigid links with spherical joints and torsional springs to represent compliance [19–24]. Examples using the technique primarily focus on modeling thin beams undergoing large static deformations with bistability, composite beams, small flexible hinges, extension mechanisms, and energy storage during motion, all of which would be computationally expensive to model relative to a rigid-body mechanics model.

Flapping wing unmanned aerial vehicles, otherwise known as ornithopters, are of interest for both civilian and military applications where tight maneuverability at low speeds is necessary [25–27]. Passive wing morphing via compliant mechanisms reduces complexity, the number of moving parts, and power consumption as compared to active morphing using electronic motors or actuation. We are developing contact-aided compliant mechanisms to passively mimic biological flight kinematics for avian-scale ornithopters during steady level flight. These passive compliant mechanisms respond to the aerodynamic loads encountered naturally during flight, allowing morphing without additional sensors and actuators. Previous work has demonstrated that implementing a single degree-of-freedom (DOF) compliant mechanism into the leading edge spar was able to improve the performance of a test ornithopter by increasing thrust and decreasing power consumption, while not producing any drag penalties [3,10,28,29]. More complex wing morphing, specifically simultaneous deformation in all three spatial directions, is required to achieve the biologically inspired continuous vortex gait [10,30]. A CCM called the bend-twist-and-sweep compliant mechanism has been designed and optimized for maximized deflections during upstroke with minimized stress and mass [7]. This compliant mechanism is tailored for simultaneous bending, twist, and sweep. These three-axis deflections are strongly coupled. For example, a bending tip load passively induces a twist angle, and applying simultaneous bending and sweep loads produces a different twist angle profile as compared to independent bending and sweep loads.

The goal of the work described herein was to tune the rigid-body model using data collected in a single experimental setup to record CCM kinematics in a real ornithopter structure, then implementing a machine learning search algorithm to find model parameters that minimize the discrepancy between the model and the experiment. The tuning parameters, experimental setup, and model tuning algorithm are described in Sec. 2. A summary of the ornithopter wing structure rigid-body model is presented in Sec. 3. The results of the experiment and model comparison are presented and discussed in Sec. 4. Finally, conclusions and recommendations for future work are presented in Sec. 5.

## Model Comparison Methodology and Experimental Setup

The DSNM is a numerical rigid-body dynamics model developed to simulate compliant mechanisms, specifically contact-aided compliant mechanisms, in a dynamic structure. The goal of the model is to provide a computationally efficient method of determining the optimal CCM stiffness and spatial locations such that they improve the performance of a structure by passively inducing shape change. Furthermore, modeling the effect of adding CCMs to a wing structure on the free flight agility of the ornithopter as a whole requires the entire wing structure to be included to account for structural reactions. The DSNM allows both the CCMs and wing structure to be modeled simultaneously in an efficient manner. The results of this model provide guidelines for developing a CCM in terms of a kinematic objective. The mathematical development of the DSNM is described in detail in Refs. [1] and [2]; the model is summarized briefly Sec. 3. The values of several parameters were unknown in the development of the DSNM. The nominal values chosen for these parameters needed to be tuned using experimental data. In this section, first, the unknown model parameters are described. Then, the experimental configurations and unknown CCM parameters are presented. Finally, the algorithm used to tune the model is presented.

### Model Parameters.

*β*. Thus, the two unknown model parameters for equivalent structural joints are the spar modulus of elasticity and spar damping coefficient term

_{s}*β*

_{s}The pseudo-rigid-body model is a common method used to model the flexible members of compliant mechanisms using rigid-body kinematics. However, the model requires a well-defined area moment of inertia to calculate the stiffness of the joints. Furthermore, the stiffness is assumed to be constant for all angles. Contact-aided compliant mechanisms such as the bend-twist-and-sweep compliant mechanism have complex, three-dimensional geometries and nonlinear stiffness due to self-contact, where the CCM has a different stiffness before and after self-contact [7]. The stiffness of such CCMs is more appropriately calculated using finite element modeling or experimental data. In the DSNM, CCMs are modeled as compliant joints: spherical joints with mass and three-axis nonlinear torsional spring dampers (Fig. 1). The damping of the CCMs is not known; the stiffness component of Rayleigh damping is assumed to represent the damping. Thus, there are four unknown parameters associated with compliant joints: their stiffness when not in self-contact and when in self-contact, the rotation required to cause self-contact, and their damping coefficient. The mathematical descriptions and nominal values are discussed in Sec. 2.2.

### Configurations Fabricated for Experiment and Compliant Joint Model.

A benchtop experiment was conducted to tune the unknown model parameters. There were four ornithopter configurations tested: the solid spar with no CCMs inserted (solid), a single degree-of-freedom CCM called the compliant spine at two spatial locations along the leading edge spar (compliant spine A (CSA)), and a modified version of the same compliant spine (compliant spine B (CSB)) (Fig. 2). A summary of the configurations and CCMs is shown in Table 1. The solid spar configuration served as a baseline for the ornithopter structure, as well as to determine suitable values for the spar properties such as modulus of elasticity and structural damping coefficient. The CCM designs were based on a compliant spine featuring two compliant joints, which has been tested previously (COMP04PM) [10,32,33]. CSA has the same design parameters and dimensions as the previously tested compliant spine, and CSB has the same design parameters but has been shortened to reduce the mass, and thereby reducing the rotational inertia penalty. The single degree-of-freedom compliant spine was chosen to consider only in-plane deflections, rather than simultaneously considering out of plane and coupled deflections, which simplified the experiment and parameter tuning in Sec. 4. Figure 3 shows CSA attached to spars via collars on either side, with seven reflective markers used in the experiment. A local coordinate system is shown in Fig. 3, where *x, y,* and *z* represent forward chord, outward span, and downward, respectively, with rotation angles *ϕ*, *θ*, and *ψ* representing bending, twist, and sweep, respectively. They were manufactured from a sheet of black Delrin using a waterjet cutter [34]. The minimum contact gap achievable by the waterjet is approximately 0.5–1 mm (0.02–0.04 in); this gap is made smaller by attaching shim stock to the CCM via double-sided adhesive VHB™ tape. A set of each configuration was created for each wing for symmetric flapping in an effort to reduce body motion from asymmetry. CSA was placed at two spatial locations: 40% of the distance from the wing root to wing tip, denoted CSA-40, and 35% of the distance, denoted CSA-35. CSB was placed at 40% of the distance. The 40% location is based on the insertion location during the previous benchtop and free flight experiments using COMP04PM [10,33], and was chosen previously based on the location of a typical bird elbow joint (40% of the half wingspan).

Configuration name | Solid | CSA-40 | CSA-35 | CSB-40 |
---|---|---|---|---|

CCM | CSA | CSB | ||

CCM length (mm (in)) | 114.3 (4.5) | 76.2 (3) | ||

CCM mass (g (oz)) | None | 14.2 (0.50) | 9.5 (0.34) | |

Total CCM mass/base ornithopter mass | 11.5% | 8.0% | ||

CCM location along half wingspan | 40% | 35% | 40% | |

Number of markers | 23 | 22 | 20 | 20 |

Flapping frequency | 5.3 Hz | 4.2 Hz | 4.1 Hz | 4.5 Hz |

Configuration name | Solid | CSA-40 | CSA-35 | CSB-40 |
---|---|---|---|---|

CCM | CSA | CSB | ||

CCM length (mm (in)) | 114.3 (4.5) | 76.2 (3) | ||

CCM mass (g (oz)) | None | 14.2 (0.50) | 9.5 (0.34) | |

Total CCM mass/base ornithopter mass | 11.5% | 8.0% | ||

CCM location along half wingspan | 40% | 35% | 40% | |

Number of markers | 23 | 22 | 20 | 20 |

Flapping frequency | 5.3 Hz | 4.2 Hz | 4.1 Hz | 4.5 Hz |

To determine their nonlinear stiffness, the CCMs were modeled using COMSOL^{®} Multiphysics, a commercial finite element analysis (FEA) program [35]. The root of the CCM was prescribed zero displacement, and a quasi-static moment was applied to the tip surface via a rigid-connector boundary condition. Contact was modeled using COMSOL's augmented Lagrangian method. The augmented Lagrangian method incorporates contact in the model by segregating the model: first, the displacement field is solved for assuming no contact; then a contact pressure term is inserted in the virtual work equation and solved for using the displacement field. The model then iteratively solves the segregated model until it fulfills the error convergence tolerance. It was assumed that there was no friction acting at the contact surfaces; the motion of the BTSCM across the angled contact gap causes perpendicular contact and therefore little slip would occur. The simulations were completed with an increasing applied moment from zero until the yield stress of the material was reached; this was performed for positive moments and negative moments separately from zero initial conditions. The applied moment verses deflection angles for CSA and CSB as predicted by the FEA are shown in Fig. 4, where positive *ϕ* is defined in Fig. 3. For moments greater than approximately −0.5 N·m, both CCMs are relatively flexible, whereas for moments less than −0.5 N·m, the CCMs come in contact and are much stiffer. The angle that the CCM comes into self-contact is called the self-contact angle, $\varphi sc$. Additionally, while both CCMs are similarly flexible when not in contact, CSB comes into self-contact at a smaller angle than CSA; this is due to the compliant hinges being closer together. The finite element analysis results shown assume a contact gap thickness of 1 mm; modifying the contact gap shifts the self-contact angle, but not the stiffnesses before or after contact. For example, CSA with a 1 mm contact gap has a contact angle of −14 deg; in contrast, CSA with a 0.5 mm gap has a self-contact angle offset of −6 deg. The width of the gap with the shim stock added proved difficult to measure accurately, so a conservative estimate of a 1 mm gap is assumed.

The response shown in Fig. 4 is approximately bilinear, where it has one stiffness when it is not in self-contact, and one relatively larger stiffness when it is in self-contact. The nonlinear stiffness introduces asymmetry into the structure's kinematics. For this application specifically, the wing morphing during upstroke is asymmetric to the wing morphing during downstroke. Therefore, the stiffnesses are referred to as symmetric, $ks$, and asymmetric, $ka$, stiffnesses, respectively. Referring to the postcontact stiffness as the asymmetric stiffness also allows the postcontact stiffness to be stated relative to the noncontact stiffness. For example, if the mechanism has no self-contact for any angle, it has a nominal symmetric stiffness for all angles and 0% asymmetry. In contrast, if the same mechanism features self-contact, and the postcontact stiffness is twice the noncontact stiffness, the nominal symmetric stiffness remains the same, but the mechanism now exhibits 100% asymmetric stiffness.

The reaction moment of a compliant joint is modeled using the piecewise bilinear function in Eq. (2), where the function is written to describe the applied moment verses rotation angle shown in Fig. 4. Note that this function describes the FEA results in Fig. 4. Equation (2) can be modified to fit any reaction moment verses rotation angle. In Eq. (2), the spring reaction torque, $Tkj$, acting on joint *j* as a function of the rotation angle, $\varphi j$, is the sum of two Hooke's law terms. The first term is the symmetric stiffness of the joint, $ksj$, multiplied by the rotation angle, i.e., linear Hooke's law. The second term accounts for the asymmetric stiffness of the joint by adding the additional torque caused by self-contact. The first component of the asymmetric term is a Heaviside step function, $u$. The step function is zero if either the self-contact angle, $\varphi sc$, is positive and the joint angle, $\varphi j$, is less than the self-contact angle, or if the self-contact angle is negative and the joint angle is greater than the self-contact angle. The self-contact angle represents the angle that the compliant joint rotates from the reference configuration before it comes in self-contact. In the rigid-body model, this represents the angle where the value of the stiffness changes. The asymmetric stiffness, $kaj$, is expressed numerically as a relative percentage of the symmetric stiffness. For example, CSA is 388% stiffer post self-contact compared to when it is not in self-contact. The second component calculates the stiffness of the joint postcontact relative to the precontact stiffness using the asymmetric and symmetric stiffnesses. The third component calculates the change in angle postcontact. For example, if the self-contact angle is −12 deg and the joint has a rotation angle of −17 deg, then the third term equals −5 deg. The nominal values of $ksj$ and $kaj$ were calculated using a forward difference method to determine the instantaneous slopes in Fig. 4, then averaged on either side of $\varphi sc$. Using the nominal stiffnesses and self-contact angles for CSA and CSB in Eq. (2), the coefficients of determination with the FEA results are 0.995 and 0.991, respectively.

### Test Platform and Camera Array.

The computer motion tracking system at the Morpheus Lab at the University of Maryland was used for the experiment (Fig. 5). The system consisted of 8 Vicon^{®} Vantage V5 cameras, with a resolution of 5 megapixels and recording speed of up to 420 frames per second [36]. The cameras are connected and powered via Ethernet cables; a software called nexus is used to record the images for each camera and then calculate the 3D positions of each marker for each frame. The positions are calibrated in the control volume using a T-shaped wand with infrared light emitting diodes spaced at known distances along the frame. The inertial coordinate system is set using the wand, which can be translated to any point in the control volume. For this experiment, the marker displacements relative to one another during flapping with the ornithopter body clamped to the ground are of interest. The cameras were placed in a circular array on tripods (Fig. 5). The cameras were recorded at 400 Hz, and were calibrated to be accurate within 0.89 mm (0.04 in), or 0.2% of the half wingspan, with the coordinate system origin approximately near the wing root.

The body of the ornithopter was clamped to a wooden platform using a combination of bar clamps and straps. Figure 6 shows the setup with the body markers and solid spar configuration. A level was used on the body near the wing roots to ensure the flapping motion was perpendicular to the ground. Four markers were placed on the body, and one marker was placed on each of the wing roots for every configuration. The markers on the body and wing root were placed to determine the orientation of the body with respect to the inertial coordinate system, and to record any body motion due to the flexibility of the body. The flapping input for the ornithopter is a four-bar mechanism attached to the shoulder joint. For the solid spar configuration, reflective tape markers were distributed along the spar. For the CCM configurations, reflective tape reflective markers were similarly distributed along the spars, with spherical reflective markers placed on top of the CCMs.

### Model Tuning Algorithm.

The method used to validate the rigid-body model of CCMs in a dynamic structure is summarized in Fig. 7. In this method, first, a rigid-body model of the vehicle is developed, modeling the compliant mechanisms using traditional rigid-body spherical joints with torsional spring dampers. Note that this algorithm is not unique to this problem; this procedure can be used to tune any rigid-body dynamics model with compliant joints using experimental kinematic data. Finite element models of the compliant mechanism under simplified boundary conditions are used to find the approximate nonlinear stiffness functions of the compliant joints, and the pseudo-rigid body model is used to approximate structural flexibility. In the case of ornithopter free flight, the aerodynamics and free flight motion are very complex, and therefore may be beyond the scope of the rigid body model. Thus, an experiment was conducted, which allows the leading edge spar to flap freely without the wing membrane or diagonal spar, but constrains the body to the ground. This allows the spar to change shape dynamically in a controlled setting with minimal aerodynamic interactions. Genetic algorithms allow “black box” system identification, where the correlation between the design variables and objective function is highly nonlinear, and therefore the best way to find the correct design variables is via a heuristic search algorithm. This is particularly useful for numerical models, where the relationship between the tuning parameters input to the model and the kinematics output from the model is very complicated. The objective function of the genetic algorithm is to minimize the error between the recorded experimental marker kinematics and the rigid-body model kinematics.

*lb*and

*ub*, respectively (Eq. (5)). The RMSE function is shown in Eq. (6), where $[rm]$ is the model kinematics, $[re]$ is the experimental kinematics,

*t*is time,

*t*is the length of the cycle,

_{f}*Y*is the distance along the spar from the wing root in the reference configuration, and

*Y*

_{tip}is the distance of the spar tip from the wing root. The integrals were calculated using trapezoidal numerical integration in MATLAB

*subject to*

## Summary of the Dynamic Spar Numerical Model

The dynamic spar numerical model was developed as a computationally efficient rigid-body model of the ornithopter wing structure with CCMs inserted. CCM information such as stiffness and mass are input in the model, geometric and inertial properties are calculated, and then a forward time integrator is used to solve a constraint-driven system of first-order differential-algebraic equations with user-input boundary conditions to solve for the state of the system and constraint reaction forces through time. The error is calculated;, then a plotting function is used to plot the system response dynamically through the flapping cycle.

Several assumptions were made in the development of the rigid-body model. The assumptions used to model CCMs were described in Sec. 2.2. All components and links—mainly the body, spars, and CJs, in the model—are assumed to be rigid, and the only relative motion is at joints and CJs. It was assumed that all connections between links are revolute and spherical joints with no friction and no joint clearance. Assumptions for the wing flapping angle as a function of time, denoted as the driving link angle, are described in Sec. 4.1. All external loading is applied at the centers of mass of the links. This focus of the current paper is on experimentally validating the structural components of the model, i.e., structural stiffness and structural damping. Previous benchtop experiments performed with the leading edge spar of the test ornithopter flapping in vacuum and ambient conditions showed very little aerodynamic effects on the spar kinematics without the wing membrane attached [32,37]. Thus, it is assumed that in the experiment, the aerodynamic forces acting on the spar and CCM without the wing membrane attached are small due to low Reynold's number, and therefore the tuned model parameters using ambient data would be very similar to using vacuum data.

The rigid-body model's geometry is based on the commercially available Park Hawk ornithopter, and is capable of handling any number of CJs in both the leading edge and diagonal spars. As an example, Fig. 8 shows a schematic for one CJ in the leading edge spar, inboard from the leading edge spar to diagonal spar interface. The inertial frame origin is located at the center of mass of the body link at the beginning of the simulation, with inertial axes *x _{I}, y_{I},* and

*z*representing nose, right wing, and underside, respectively. Angles about the inertial axes are

_{I}*ϕ*(roll),

*θ*(pitch), and

*ψ*(yaw). The body is clamped at the center of mass, and the wing root is modeled as a revolute joint about

*x*

_{1}. Black lines along the leading edge and diagonal links represent the stiff spars. The local angles between adjoining links are $\varphi j,\theta j,$ and $\psi j$. In a flexible wing model, these angles along the leading edge spar correspond to bending in the direction of flapping, twist about the span direction, and sweep in the chord direction, respectively. The driving link represents the ornithopter flapping mechanism input at the shoulder joint.

*{r}*, Euler parameters,

*{p}*, inertial linear velocity, $r\u02d9$, and local angular velocity,

*{ω}′*. The DSNM is based on kinematic constraints, {

*Φ*}. The constraints are split into five sections: driving link angle, symmetric flight, clamped body, spherical joints, and interface between the leading edge and diagonal spars. A complete description of the constraints used and model development can be found in Refs. [1] and [2]. All of the terms are then compiled into a mixed system of differential-algebraic equations of motion

The first row represents the translational component of the constrained Newton–Euler equations of motion, where the first term is the system mass matrix, $[M]$, multiplied by the translational acceleration, ${r\xa8}$, the second term is the position Jacobian, $[\Phi r]$, multiplied by Lagrange multipliers, $\lambda $, and the third term is the inertial applied forces, ${F}$. The second row represents the rotational component of the constrained Newton–Euler equations of motion. The first term is the system mass moment of inertia matrix, $[I]\u2032$, multiplied by the angular accelerations, ${\omega \u02d9}\u2032$. The second term is the orientation Jacobian, $[\Phi \pi \u2032]$, multiplied by the Lagrange multipliers. The third term is the applied torques on each link in local coordinates, ${T}\u2032$, minus the body-fixed component of the definition of the moment of a body with respect to a body-fixed coordinate system. The parameters used to tune the model appeared in Eqs. (1)–(3); these equations are used to calculate the applied torques due to the joint stiffness and damping. The final row represents the acceleration equations associated with the kinematic constraints.

## Experimental Results and Comparison to Model

In this section, the driving link angle, pseudo-rigid body joint stiffness and damping, and CCM stiffness and damping are assessed. An interpretation of the model tuning is also provided. Finally, sources of experimental uncertainty are discussed. All experimental data presented in this section are relative to the wing root marker. For visualization purposes, three markers are shown for each configuration: the inboard and outboard markers of the CCM, *M _{i}* and

*M*, and the tip marker,

_{o}*M*(Fig. 9). For the solid configuration, the CCM markers are the ones at similar spatial locations to those of CSA at 40%.

_{t}### Extraction of the Experimental Data and Driving Link Angle.

In the formulation of the DSNM, the driving link angle constraint prescribes the rotation angle as a function of time and angular acceleration of the driving (first) link along the leading edge spar. Therefore, a twice-differentiable driving link angle function is necessary. A sinusoidal input to the model was unable to replicate the rich response measured in the experiment. Thus, the experimental driving link angle was calculated and used for the model. This angle was calculated using the position of the first marker on the spar with respect to the wing root marker, effectively measuring the angle of the wing root shoulder joint with respect to the ground. Mid upstroke is denoted as the point in the flapping cycle where the shoulder joint is parallel to the ground with upward velocity (i.e., velocity in −*z*). This instant in time is not likely to occur at a precise time-step recorded by the camera array. A simple option would be to prescribe the wing root angle as interpolated values from the experimental data; however, this would not account for the camera uncertainty. The kinematics of all the cycles recorded were overlaid on top of one another; then a moving average function was used to determine an average flapping cycle. First, the midupstroke time was found by interpolating the time when the driving link angle crossed the horizontal plane in the upward direction for each cycle. Subtracting each cycle's beginning time from the following cycle's beginning time yielded the cycle's period; this value was used to determine the spar's average flapping frequency. The time for each data point in the cycle was normalized by the cycle's period; this eliminated artificial harmonics in the moving average due to small changes in the flapping frequency. A moving average was performed for each marker's *XYZ* positions and the driving link angle, where the sliding window was equal to the number of cycles for the solid configuration (107), and 40 for the CCM configurations because only 20, 10, and 25 cycles were recorded for CSA-40, CSA-35, and CSB-40, respectively, and a larger window would reduce artificial harmonics. The first and second time derivatives of the driving link angle were calculated using Klein et Morelli's global Fourier smoothing technique [39]. The averaged position and driving link data were then resampled to 800 Hz, twice the camera sampling frequency, to improve the comparison between the model and experiment. As an example, Fig. 10 shows the comparison between the resampled average driving link angle and the overlaid cycles. Figure 10 also shows a 95% confidence interval, as well as a 2.3 deg camera uncertainty. This uncertainty is discussed in Sec. 4.4.

### Solid Spar Configuration.

The solid spar configuration test was performed primarily to assess the assumption that the spar allows only allows small deformations during flapping. A schematic of the configuration is shown in Fig. 9 (left). MATLAB's genetic algorithm from the global optimization toolbox [40] was used to solve the optimization problem, where the design variables were the elastic modulus, $Es$, and damping coefficient, $\beta s$, of the spar. Three configurations are compared: assuming the spar is rigid, using the nominal modulus of elasticity for the spar and an estimated damping coefficient, and the tuned modulus and damping coefficient. The configurations, driving link RMSE, and marker RMSE are shown in Table 2. The modulus of elasticity was reduced by 20% compared to the spar specification. All configurations had a driving link RMSE of 0.67 deg, or 1.2% of the wing stroke angle. A rigid spar assumption produced a marker RMSE of 2.5% of the half wingspan. Including equivalent structural joints with nominal values decreased the marker RMSE to 1.2%. The parameters found using the GA decreased the RMSE to 0.9% of the half wingspan. Three markers were chosen to illustrate the comparison of the experiment to the model: the markers at 25% and 58% of the half wingspan, and the spar tip. Figure 11 shows the rigid spar results. Note that only 25% of the experimental data points are shown for illustration purposes. The largest discrepancy between the model and experiment is at the stroke transitions; this is because the spar flexes due to its inertia. Figure 12 shows the three markers for the optimized configuration. Good agreement is seen at all times during the flapping cycle, and the flexibility introduced by the equivalent structural joints allows the model spar tip to deflect similar to the experiment at the transitions.

Configuration | Rigid | Normal | GA |
---|---|---|---|

Population-generations | N/A | N/A | 30–26 |

β_{s} (N·s/rad) | (Rigid) | 1.0 × 10^{−3} | 6.8 × 10^{−4} |

E_{s} (%) | (Rigid) | 100 | 80 |

Driver RMSE (deg) | 0.67 | 0.67 | 0.67 |

Marker RMSE (%) | 2.5 | 1.2 | 0.9 |

Configuration | Rigid | Normal | GA |
---|---|---|---|

Population-generations | N/A | N/A | 30–26 |

β_{s} (N·s/rad) | (Rigid) | 1.0 × 10^{−3} | 6.8 × 10^{−4} |

E_{s} (%) | (Rigid) | 100 | 80 |

Driver RMSE (deg) | 0.67 | 0.67 | 0.67 |

Marker RMSE (%) | 2.5 | 1.2 | 0.9 |

### CCM Configurations.

The goal of the CCM configurations was to assess the CJs stiffness found via finite element simulations, and to empirically develop an appropriate damping coefficient for the CJs. A schematic of the CCM configuration is shown in Fig. 9. For each configuration, a simulation using the nominal finite element derived values was performed; then the optimization problem described in Sec. 2.4 was solved using a genetic algorithm. The tuned spar modulus-damping parameter from the solid spar configuration was used for equivalent structural spar joints inboard and outboard from the CCM for both the nominal and genetic algorithm simulations. The optimization design parameters were the symmetric and asymmetric stiffnesses, $ks$ and $ka$, respectively, and the self-contact angle, $\varphi sc$. The results of the genetic algorithms are shown in Table 3. The optimized configurations reduced the CCM RMSE by between 61% and 75%, and marker RMSE by between 1.3% and 2.8% of the half wingspan. All three configurations had CCM RMSEs less than 2 deg and marker RMSEs less than 2% of the half wingspan. The finite element model underpredicted the symmetric stiffnesses, *k _{s}*, and overpredicted the asymmetric stiffnesses,

*k*. The damping coefficient of CSA,

_{a}*β*, was found to be double the damping coefficient of the spar equivalent structural joints,

_{c}*β*. The damping coefficient of CSB was nearly four times that of CSA. In the formulation of CJ damping, only the damping term from Rayleigh damping proportional to stiffness was considered, not the mass term. The difference in coefficients can be attributed to the difference in mass between CSA and CSB. The finite element model overpredicted the self-contact angle,

_{s}*ϕ*

_{sc}, by 3–5 deg. A comparison of the CSA-40 model and experiment is shown in Figs. 13 and 14. Very close agreement is seen for the driving link angle with an RMSE of 0.15 deg (Fig. 13 top), and a CCM angle RMSE of 1.4 deg (Fig. 13 bottom). The magnitude and phase of the CCM oscillations are very close (Fig. 13 bottom), where each high frequency oscillation from the experiment is present in the model, however slightly more damped. Figure 14 shows the CCM root and tip and spar tip markers from the model compared to the experiment. A small lag occurs between the inboard markers and tip markers, specifically at the transitions where the tip continues to flex after the shoulder joint has reached its maximum angle. This is due to the inertia of the CCM resisting changing flapping direction; this lag was not present during the solid spar configuration. The model was able to match the marker kinematics very well throughout the flapping cycle with the exception of the upstroke to downstroke transition. This is likely due to the complicated dynamic self-contact phenomena. Similar results were seen for CSA-35 and CSB-40. All three CCM configurations achieved similar angle magnitudes throughout the flapping cycle, with CSB having the smallest peak angle of the three configurations. Interpretation of these results is discussed in Sec. 4.4.

Configuration | CSA-40 | CSA-35 | CSB-40 | |
---|---|---|---|---|

Population-generations | 100-9 | 100-11 | 100-9 | |

K_{s} (%) | Nominal | 100 | ||

GA | 135 | 116 | 153 | |

K_{a} (%) | Nominal | 100 | ||

GA | 82 | 61 | 70 | |

β_{c} (N·s/rad) | Nominal | 1.0 × 10^{−3} | ||

GA | 1.4 × 10^{−3} | 1.5 × 10^{−3} | 5.9 × 10^{−3} | |

ϕ_{sc} (deg) | Nominal | −14 | −12 | |

GA | −10.7 | −8.1 | −9.4 | |

Driver RMSE (deg) | Nominal | 0.15 | 0.21 | 0.22 |

GA | ||||

CCM RMSE (deg) | Nominal | 5.7 | 3.5 | 5.1 |

GA | 1.4 | 1.3 | 2.0 | |

Marker RMSE (%) | Nominal | 4.3 | 2.9 | 3.8 |

GA | 1.5 | 1.6 | 1.7 |

Configuration | CSA-40 | CSA-35 | CSB-40 | |
---|---|---|---|---|

Population-generations | 100-9 | 100-11 | 100-9 | |

K_{s} (%) | Nominal | 100 | ||

GA | 135 | 116 | 153 | |

K_{a} (%) | Nominal | 100 | ||

GA | 82 | 61 | 70 | |

β_{c} (N·s/rad) | Nominal | 1.0 × 10^{−3} | ||

GA | 1.4 × 10^{−3} | 1.5 × 10^{−3} | 5.9 × 10^{−3} | |

ϕ_{sc} (deg) | Nominal | −14 | −12 | |

GA | −10.7 | −8.1 | −9.4 | |

Driver RMSE (deg) | Nominal | 0.15 | 0.21 | 0.22 |

GA | ||||

CCM RMSE (deg) | Nominal | 5.7 | 3.5 | 5.1 |

GA | 1.4 | 1.3 | 2.0 | |

Marker RMSE (%) | Nominal | 4.3 | 2.9 | 3.8 |

GA | 1.5 | 1.6 | 1.7 |

### Interpretation of Results and Experimental Uncertainty.

The VICON cameras were calibrated to record marker positions within 0.89 mm, or 0.2% of the half wingspan. The outboard marker on the CCM oscillates with a peak-to-peak amplitude of 322 mm, or 64% of the half wingspan. The spar tip oscillates with a peak-to-peak amplitude of 611 mm, or 120% of the half wingspan. Thus, the marker uncertainty is very small relative to the motion of the spar. The experimental kinematics used to tune the model are averaged over 107 cycles for the solid configuration, and at least seven cycles for the three CCM configurations. Additionally, the cameras recorded the kinematics 2 orders of magnitude faster than the flapping frequency. When considering the uncertainty propagation associated with the marker data used to calculate the wing root angle and CCM angles, the uncertainty of the two markers was calculated to be approximately 2.3 deg. This number was found using the inverse tangent of double the marker position accuracy divided by the measured distance between the markers. There is also uncertainty associated with the marker size; the cameras denote the marker position as being the center of the reflective surface; however, due to the size and shape of the markers and reflective tape used, the relative locations of the markers may not be exactly as measured. Marker dropout occurred due to the small markers moving at high speeds; however, this was eliminated by resampling the data for multiple flapping cycles.

The optimized parameters relative to the nominal values, as well as a comparison of the RMSEs, were shown in Tables 2 and 3. The driving link angle RMSE was less than 1 deg for all configurations, within the experimental uncertainty. The solid spar configuration isolated the reaction of the spar, allowing us to tune the spar parameters independently from the CCM. Equivalent structural joints reduced the solid configuration marker RMSE from 2.5% of the half wingspan using a rigid spar assumption to 1.2%, and furthermore to 0.9% using the parameters found using the GA. The rigid spar assumption case's error was primarily due to the spar flexing at the stroke transitions. Adding equivalent structural joints allowed the model to capture this flexibility, and tuning the model reduced the error further. The tuned parameters were then used for the CCM configuration simulations. The nominal CJ stiffness parameters were found using a finite element model, and the damping parameter was estimated. Tuning the stiffness and damping parameters reduced the CCM angle RMSE by 60% to 75% to less than 2 deg, and the marker RMSE from a range of 3–4% down to around 1.6%. For all CCM configurations, the nominal symmetric stiffnesses were higher than the nominal values, and asymmetric stiffnesses were lower than the nominal values. The self-contact angle was found to be between 8 deg and 10 deg, which is 3–5 deg less than predicted. This was due to the shim stock in the contact gap which was not accounted for in the model. All three configurations had similar bending angle magnitudes. This means the change in CCM insertion location and reducing the mass while maintaining the same design parameters had a relatively small effect on the CCM response. While CSB-40 would be expected to have a smaller response due to the decreased inertia, the driving link was flapping 0.5 Hz faster than CSA-40. This increased the kinetic energy from flapping, which was then converted into elastic potential energy at the transitions. The bending angle for CSA-35 was similar to CSA-40, indicating that a 5% difference in insertion location has a small effect on the CCM response under benchtop conditions. Overall, the equivalent structural joints with optimized parameters were able to account for the flexibility of the spar, and the optimized CJ models were able to match the experimental deformations, producing an error of less than 2%.

Additional uncertainty can be attributed to the motor being modeled as a prescribed angle, the motor heating over time, and the motion of the body potentially inducing inertia into the system. The clamps constraining the ornithopter were only able to constrain the bottom of the body due to the onboard electronics, allowing the top of the body to undergo small lateral oscillations of less than 2 mm during flapping. The body and wing root were leveled using a liquid air bubble; therefore, a small uncertainty exists in determining the inertial coordinate axis directions, causing small roll and yaw in the kinematic data. The black Delrin used to fabricate the CCMs may have different mechanical properties than the published data for white Delrin [34]. The collar to CCM connection had a small clearance, and the CCM and collars were tapered due to the waterjet manufacturing. Furthermore, all material properties used in the FEA assumed quasi-static deformation, whereas the experiment is flapping rapidly with little aerodynamic damping, only structural damping, to damp large deformations from inertia. The experimental tip marker vertical position for the CCM cases was larger in magnitude at the transitions compared to the model; this likely indicates that the CCM's stiffness decreased at larger angles, whereas in the model, the stiffness is bilinear and does not change with larger angles. With the CCM installed, the outboard spar moves much quicker than the solid spar configuration due to the CCM flexibility, making it much more difficult for a camera to track closer markers near the tip.

## Conclusions and Future Work

A numerical dynamic model of contact-aided compliant mechanisms in a flapping spar was tuned and validated using a benchtop test. The goal of this test was to tune five model parameters that were unknown or based on finite element models. A flapping angle was calculated using marker kinematics averaged over many flapping cycles. A genetic algorithm was used to minimize the root-mean-square error between the experimental marker kinematics and their corresponding model marker kinematics. Equivalent structural joints successfully accounted for the flexibility of the spar, and a damping coefficient and relative modulus of elasticity of the spar were found. The tuned stiffness of the CCMs was found relative to the FEA predicted stiffness for all three configurations. An approximate Rayleigh damping coefficient and the CCM self-contact angle offset were also calculated. The genetic algorithm successfully reduced the model marker RMSE to less than 2% of the half wingspan for all four configurations.

This experiment considered a single degree-of-freedom CCM inserted in the leading edge spar of a benchtop clamped ornithopter. The model is capable of modeling both the leading edge and diagonal spars. Ongoing work includes using the tuned rigid-body model to optimize the spatial location and stiffnesses of multiple compliant mechanisms inserted into the wing structure. The results of the optimization will give insight into the performance requirements of a new CCM, which will improve the free flight agility of the ornithopter.

## Acknowledgment

The authors gratefully acknowledge the resources of the NASA Langley Research Center, The Pennsylvania State University, University of Maryland, Star Lab, the Engineering Design and Optimization Group, and Morpheus Lab.

## Funding Data

This work was funded by the Air Force Office of Scientific Research (AFOSR) (Grant No. FA9550-13-0126 under the direction of Dr. David Stargel and Mr. James Fillerup).

## Nomenclature

- CCM =
contact-aided compliant mechanism

- CJ =
compliant joint

*F*=force

- GA =
genetic algorithm

*I*=inertial frame

- [
*I*]′ = local inertia matrix

*j, p, q*=joint number

*J′*=area moment of inertia

*k*=joint stiffness

*M*=applied Moment

- [
*M*] = mass matrix

- {
*p*} = Euler parameter vector

- {
*r*} = position vector

- RMSE =
root-mean-square error

- {
*s*} = position vector of points on a body

*t*=time

- {
*T*}′ = local torque vector

*x, y, z*=Cartesian coordinate system

- {
*y*} = state variable

*β*=Rayleigh damping coefficient

- {
*γ*} = acceleration equation

- {
*λ*} = Lagrange multiplier

*ϕ,θ,ψ*=roll, pitch, and yaw angles

*ϕ*=^{j},θ^{j},ψ^{j}local rotation angles of joint j

*ϕ*_{Dr}=driving link angle function

*ϕ*_{sc}=contact angle offset

- {
*Φ*} = constraint vector

- [
*Φ*] =_{r} constraint position Jacobian

- [
*Φ*] =_{π′} constraint orientation Jacobian

- {
*Φ*} =_{t} constraint time derivative

- {
*ω*}′ = local angular velocity vector

Braces indicate a column vector and square brackets indicate matrices. Internal numerical subscripts represent link numbers. Internal numerical superscripts represent joint numbers, and alphabetic subscripts represent partial derivatives. External subscripts represent matrix row and column references, and external superscripts represent reference frames. Dots represent time derivatives, hats represent unit vectors, tildes represent skew symmetric matrices, bars represent constant values, and primes indicate local reference frames.