Flight Simulation of an Ornithopter - Semantic Scholar
Flight Simulation of an Ornithopter
Cameron Rose Ronald S. Fearing
Electrical Engineering and Computer Sciences University of California at Berkeley Technical Report No. UCB/EECS-2013-60 http://www.eecs.berkeley.edu/Pubs/TechRpts/2013/EECS-2013-60.html
May 14, 2013
Copyright © 2013, by the author(s). All rights reserved. Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior specific permission.
ii
Contents 1 Simulation Environments 1.1 Aerodynamic and Dynamic Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Yaw Dynamics Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2 3 5
2 Sagittal Plane Experimental Data 2.1 Wind Tunnel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Free Flight . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Comparison of Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7 7 8 8
3 Simulation Results 15 3.1 Free Flight Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4 Yaw Plane Simulation 19 4.1 Yaw Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.2 Yaw Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 5 Conclusions and Analysis 5.1 Sagittal Plane Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Yaw Plane Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23 23 24 24
Bibliography
25
2
Chapter 1
Simulation Environments The ornithopter simulation software package includes two modules: a sagittal plane model and a yaw plane model. Each module is implemented separately and models the H2 Bird hardware [7]. The robot and its control surfaces are shown in Figure 1. The H2 Bird is built around the Silverlit i-Bird RC flier power train and has clap-fling wings with a wingspan of 26.5 cm. For yaw and pitch control, it uses a tail propeller and a servo-controlled tail elevator. In terms of control and sensing, the H2 Bird uses an on-board ImageProc 2.4 controller that includes a 40 MIPS microprocessor, 6 DOF IMU, IEEE 802.15.4 radio, and motor drivers, powered by a 90 mAh lithium polymer battery [1].
Wings
Tail prop ImageProc Elevator servo Tail elevator
Figure 1.1: H2 Bird ornithopter MAV. The sagittal plane ornithopter simulation is modeled using rigid body dynamics with force and moment information collected from wind tunnel and free flight tests as inputs. The top-level layout of the simulator is illustrated in Figure 1.2 and is implemented in MATLAB. The duty cycle of the motor wing drive is the input, and the state of the ornithopter is updated at 1 kHz. The yaw plane simulation is also implemented in MATLAB. The yaw dynamics are determined using system identification on the step response of the ornithopter to a 90 degree right turn. The simulation computes the yaw motion of the ornithopter with a desired heading as the input.
CHAPTER 1. SIMULATION ENVIRONMENTS
3
Body Parameters
Aerodynamics Duty Cycle Controller
Body Dynamics Forces, Moments
State(t)
State(∆t)
∫ 1kHz
Figure 1.2: Simulator block diagram.
1.1
Aerodynamic and Dynamic Models
The aerodynamics of flapping winged flight are nonlinear and complex. The flapping wings disturb the airflow and create vortices that interact with the wings and tail of the ornithopter to augment the thrust forces produced by the flapping of the wings. In some cases, simply the thrust provided by the wings is not sufficient by itself to keep the robot aloft, although it still flies freely [2]. Several approaches to modeling the flight paths of ornithopters have been used, but they often involve variations of blade element theory or other methods to model the dynamics of the wings in an effort to predict the flight path [3, 5, 10]. In my model, I simplify the process by using averages of measured forces and moments to predict the motion of the ornithopter. The method of using averaged data has been previously shown to be a valid approximation [8, 11]. This strategy enables a simulation that is less computationally complex, which is helpful for simulating multiple ornithopters or working with more intricate control strategies. The previous control strategies used for earlier iterations of the ornithopter were based upon Proportional-Integral-Differential (PID) control schemes for target tracking or height regulation [1, 2]. Information about the aerodynamic interactions of the ornithopter can produce a more robust model that can be used in more sophisticated control schemes such as Linear-Quadratic Regulator (LQR) or Model Predictive Control (MPC). The aerodynamic forces and moments used in the simulation are collected using a force and torque sensor affixed to the bottom of the ornithopter in a wind tunnel. The forces and moments were measured at the center of mass of the entire robot to attempt to capture the interactions between the wings and the rest of the robot without separately modeling each interaction. These measurements were corrected by free flight data collected using a Vicon motion capture system. The modification of the fixed wind tunnel data is necessary, as the raw fixed data are not sufficient to capture the complex nonlinear dynamics of the ornithopter [4]. The pitching motion that the ornithopter exhibits in free flight is important to the flight dynamics, although this could not be captured in the wind tunnel. The horizontal and vertical forces and the pitching moment are used as inputs into the body
CHAPTER 1. SIMULATION ENVIRONMENTS
4
dynamics in the simulation. Each force and moment is stored at set data points and interpolation is used to calculate the points in between. The forces and moments are all functions of the set duty cycle, angle of attack, and velocity magnitude of the ornithopter at a given point in time. The body dynamics of the ornithopter are treated as a rigid body. The kinematic equations of motion are as follows [6]: x˙ = u cos θ + w sin θ z˙ = u sin θ − w cos θ 1 u˙ = −wQ + (T − mg sin θ) m 1 w˙ = uQ + (−L + mg cos θ) m θ˙ = Q M − damping Q˙ = Iyy
(1.1)
The state vector has six elements: x is the horizontal position in world coordinates, z is the vertical position in world coordinates, u is the forward velocity in body coordinates, w is the vertical velocity in body coordinates, θ is the pitch, and Q is the angular velocity. Additionally, T is the thrust and L is the lift in body coordinates, and M is the pitch moment. The forces are measured by the sensor in body coordinates, so the velocities are in body coordinates for ease of use. In terms of the body parameters, the pitch moment of inertia, Iyy , is 3.0×10−5 kg·m2 and was computed using SolidWorks, which computes the mass distribution of the model of the ornithopter and returns it in the form of an inertial matrix. The mass of the ornithopter, m, is 13.6 g, and was measured using a scale. The body diagram is shown in Figure 1.3.
zw
xb
L u T
𝛂
𝚹
V mg
w zb
Figure 1.3: The rigid body diagram for the model.
xw
CHAPTER 1. SIMULATION ENVIRONMENTS
1.2
5
Yaw Dynamics Model
The turning dynamic model is part of a larger model used to simulate the motion of the ornithopter in the horizontal plane in a cooperative control scenario. This model is used to predict the feasibility of achieving a goal successfully for a given initial pose. In the scenario, a ground station is used to provide guidance to the ornithopter by transmitting desired headings to a window calculated from the visual input provided by an off-the-shelf USB web camera. The layout of the experiment is depicted in Figure 1.4. The total cooperative system software layout is shown in Figure 1.5.
H2Bird
0.66 m
Window
0.94 m
Camera 3.66 m Ground Station
Figure 1.4: Conceptual sketch of experimental setup. The yaw dynamics of the ornithopter are modeled in the horizontal plane using a one dimensional model to limit the complexity. This format also matches the physical turning implementation on the actual ornithopter. In terms of control, the ornithopter is sent a desired heading and it controls to that heading using an on-board PID controller. The turning dynamics and the internal PID controller are modeled with a single one-dimensional model: θ(s) =
K u(s) 1 + Rs
(1.2)
where K and R are constants determined through system identification and the input, u(s), is the desired heading and the output, θ(s), is the physical heading of the ornithopter. A form of Dubin’s car model is used to model the translational dynamics in the horizontal
CHAPTER 1. SIMULATION ENVIRONMENTS
6
plane [9]: x˙ = us sin θ y˙ = us cos θ θ˙ = u
(1.3)
u ∈ U = [− tan φmax , tan φmax ] In the car model, φmax is the maximum steering angle and us is the forward speed of 1.2 m/s. Here, the angular position θ is an input rather than a state. The overall motion of the system can be modeled as shown in Equation 1.3, where θ is the time-domain output of Equation 1.2, θ(t).
External Control Loop – 10Hz Window position
Desired Internal Control Loop – 300Hz Heading PID P-D Ornithopter Controller Controller Dynamics
-
Translational Position Pose Estimation
Angular position
Attitude Estimation
Camera
Figure 1.5: Overall block diagram of the cooperative control system.
Pose
7
Chapter 2
Sagittal Plane Experimental Data The experimental data necessary to determine the aerodynamic properties of the ornithopter are comprised of both wind tunnel force-torque data and free flight data collected from a Vicon motion capture system.
2.1
Wind Tunnel
The wind tunnel data were collected over a series of trials using an ATI Nano17 force-torque sensor. The ornithopter was affixed to an acrylic mount attached to the sensor and facing into the wind, as shown in Figure 2.1. The ornithopter was attached at approximately its center of mass, although this point fluctuates as the wings open and close. This fluctuation is minimal, however, so it was discounted. Each trial consisted of 3 seconds of data collected for wind speeds between 0 and 3 m/s in 0.5 m/s increments and longitudinal body axis angles relative to the horizontal, or angles of attack, between -70 and +90 degrees in 10 degree increments. Each of these trials was conducted with the wings of the robot closed and stationary. Each trial was averaged over the 3 seconds of collection, and only the thrust and lift forces in body coordinates and the pitch moment were stored. The results of the experiments are in Figures 2.2 through 2.4. Each force is expressed in the world reference frame with the weight of the robot absent from the measurements. In addition to the stationary wing experiments, data were collected in still air for duty cycles between 50 percent and 100 percent in 10 percent increments using the same wind tunnel setup as the previous experiments. This range corresponds to wing flap speeds between 12 Hz and 20 Hz. These experiments required some damping around the connection of the robot to the acrylic mount with 4 mm of foam, as the oscillations from the wings moving created significant noise in the pitch moment signal. The data collected from the experiments were averaged over 5 seconds, and the thrust and lift forces in body coordinates and the pitch moment were stored. The results of the experiments are in Figures 2.5 and 2.6. Both of these results were combined to form a 4 dimensional data set dependent upon duty cycle, angle of attack, and wind speed that outputs the forward and vertical body frame forces and pitch moment. This data set is used as a look-up table for the instantaneous forces and moments for a given pose in the simulation.
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
8
L T
Figure 2.1: H2 Bird mounted to sensor
2.2
Free Flight
The free flight data was captured over variable amounts of time for circular paths in an empty room for duty cycles between 70 and 100 percent in 5 percent increments. For each trial, the initial launch transient was cut from the data, because only the behavior around equilibrium was of interest. For each trial, the position, orientation, angle of attack, and velocity magnitude and direction of the ornithopter were computed. Data from a typical trial are shown in Figures 2.7 and 2.8. The relationship between duty cycle and velocity, angle of attack, and pitch angle are depicted in Figures 2.9 through 2.11. The trends and locations of the equilibrium points in terms of the angle of attack and velocity are used in conjunction with the wind tunnel data to determine an accurate picture of the flight characteristics for simulation.
2.3
Comparison of Approaches
The two approaches to determining the flight forces of the ornithopter both have their merits and limitations. The wind tunnel data give a picture of the aerodynamics of the ornithopter over a broad range of wind conditions. Although all of the motions that contribute to the flight of the ornithopter are not present in this fixed measurement, the wind tunnel data still provides an indication of the relationship of the forces and moments to changes in wind speed and angle of attack. While these relationships can be used, the values of the individual data points may diverge, as some of the aerodynamic interactions are not present. The free flight data provide information about the exact values of the angle of attack and flight speed at equilibrium for given duty cycles. This information can be used to correct the values of the forces and moments produced by the wind tunnel experiments, but the relationships of the wind tunnel data can still be used. One limitation of the free flight data is that it only gives an indication of the flight behaviors at equilibrium, which is a narrow snapshot of the flight dynamics.
0.2
Vertical Force [N]
Horizontal Force [N]
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
0.1 0 −0.1 −0.2 0 1 2 3
0.1 0 −0.1 −0.2 0 1 2 3
AoA [Degrees]
Wind Speed [m/s]
Wind Speed [m/s]
Figure 2.2: Aerodynamic horizontal force surface in world coordinates. Pitch Moment [N*mm]
0.2
50
0
−50
9
−50
0
50
AoA [Degrees]
Figure 2.3: Aerodynamic vertical force surface in world coordinates.
5
0
−5
−10 0 1 2 3
Wind Speed [m/s]
−50
0
50
AoA [Degrees]
Figure 2.4: Aerodynamic pitch moment surface.
Additionally, the flapping wings produce noise in the Vicon tracking system, so it is difficult to capture valid data. To combine the two sets of data into one composite data set to use with the simulator, several adjustments were made to the wind tunnel data that were inspired by the free flight data. One of these adjustments was to the equilibrium points indicated by the wind tunnel data, which did not match with the values of the equilibrium points indicated by the free flight data. This issue was expected, since the ornithopter is fixed in the wind tunnel. To account for the difference, the zero force and moment point of the wind tunnel data was shifted to the point supplied by the free flight data. This change was implemented by subtracting, from the entire wind tunnel data set, the values provided by the wind tunnel at the free flight provided equilibrium point: 85 percent duty cycle, 1.3 m/s wind speed, and 53 degrees angle of attack. Therefore, the values of the data were shifted, but all of the relationships remain the same. Another discrepancy between the wind tunnel data set and the free flight data, is that in the wind tunnel data, the relationship between the pitch moment and the wind speed was incorrect. As shown in Figure 2.12, the measured pitch moment decreases as the wind speed increases. If this were the case, the robot would pitch down as the speed increases, which would cause a further speed increase. This is counterintuitive, and also contradictory of the actual flight behavior. This
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
0.1
−1
Thrust Force Lift Force
Pitch Moment [N*mm]
0.12
Forces [N]
0.08 0.06 0.04 0.02 0 −0.02 50
60
70
10
80
90
Duty Cycle [Percent]
100
Figure 2.5: Thrust (blue) and lift (green) forces as a function of the duty cycle.
−1.1 −1.2 −1.3 −1.4 −1.5 50
60
70
80
Duty Cycle [Percent]
90
100
Figure 2.6: Pitch moment as a function of the duty cycle.
deficiency could be caused by some aerodynamic interaction that is not captured in the wind tunnel that causes the robot to pitch up in free flight. To fix this problem, a further additive term was added to the pitch moment dependent on the wind speed: 7 ∗ windspeed . This term is meant to 3.0 behave like an elevator would on the end of the tail. This approximation is overly simplistic, as, in actuality, the moment provided by an elevator would also be dependent on angle of attack. However, this term fixes the relationship between wind speed and pitch moment, while preserving the pre-existing correct relationship between angle of attack and pitch moment. The progressive effects of these changes to the pitch moment surface are illustrated in Figures 2.14 through 2.16. The shift illustrated in Figure 2.15 is analogous to a similar shift in the forward and vertical forces based upon the equilibrium points determined in the free flight experiments. Figure 2.13 gives a closer look at the adjusted relationship between wind speed and the pitch moment for a particular angle of attack.
Pitch [Deg]
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
100 50 0 0
Z [m]
1
0.5
2
0 −2 1
−1 0
0 1 −1
X [m]
Y [m]
Figure 2.7: The Vicon tracked path for 80 percent duty cycle with velocity vector (blue) and angular orientation vector (green) indicated.
1
2
3
4
5
2
3
4
5
2
3
4
5
Time [s]
100 50
Velocity [m/sec]
AoA [Deg]
2
1.5
11
0 0
1
Time [s]
4 2 0 0
1
Time [s]
Figure 2.8: The Vicon tracked pitch angle, angle of attack, and velocity magnitude for 80 percent duty cycle.
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
70
Angle of Attack [Degrees]
54
Pitch [Degrees]
52 50 48 46 44 42 60
12
70
80
90
100
60 55 50 45 40 35 30 60
110
Duty Cycle [Percent]
65
Figure 2.9: Free flight pitch angle as a function of the duty cycle.
70
80
90
Velocity [m/s]
2
1.5
70
110
Figure 2.10: Free flight angle of attack as a function of the duty cycle.
2.5
1 60
100
Duty Cycle [Percent]
80
90
100
Duty Cycle [Percent]
110
Figure 2.11: Free flight velocity magnitude as a function of the duty cycle.
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
2
0.5
Pitch Moment [N*mm]
Pitch Moment [N*mm]
1
0 −0.5 −1 −1.5 −2 −2.5 0
13
0.5
1
1.5
2
Wind Speed [m/s]
2.5
3
Figure 2.12: Original pitch moment for 85 percent duty cycle and 53 degrees angle of attack as a function of wind speed.
1
0
−1
−2
−3 0
0.5
1
1.5
2
Wind Speed [m/s]
2.5
3
Figure 2.13: Adjusted pitch moment for 85 percent duty cycle and 53 degrees angle of attack as a function of wind speed.
Pitch Moment [N*mm]
Pitch Moment [N*mm]
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
5
0
−5 0 1 2 3
−5 0 2 3
Wind Speed [m/s]
Figure 2.14: Initial total pitch moment surface.
Pitch Moment [N*mm]
0
1
AoA [Degrees]
Wind Speed [m/s]
5
50
0
−50
14
0
−5
1
Wind Speed [m/s]
2
0
50
AoA [Degrees]
Figure 2.15: Equilibrium shifted total pitch moment surface.
5
0
−50
3
−50
0
50
AoA [Degrees]
Figure 2.16: Final adjusted total pitch moment surface.
15
Chapter 3
Simulation Results To determine the accuracy of the simulation, open loop flight tests of the sagittal plane model were conducted. The simulation results were compared to analogous free flight tests.
3.1
Free Flight Comparisons
For the wind tunnel inspired model, the motion of the ornithopter was simulated for duty cycles between 70 and 100 percent in 5 percent increments. The initial conditions chosen were 0.85 m/s forward velocity in the horizontal plane and a 40 degree pitch angle. For each trajectory, the path of the flight and also the angle of attack, velocity magnitude, and pitch angle were plotted, as shown in Figures 3.1 through 3.4. To determine the overall steady-state accuracy of the simulation in comparison to the real flight data, the simulation results were compared to the previously collected flight data. As mentioned in the previous chapter, this flight data was collected for duty cycles between 70 and 100 percent in 5 percent increments for circular paths. The paths are circular because the room was not big enough to collect multiple seconds of data after the initial launch transient before equilibrium flight is reached. Because of this difference, the steady state values are not exact for straight forward flight. Since the flapping wings are not modeled in the simulation, the oscillatory behavior of the true flight path is not captured. However, the steady state average angle of attack, pitch, and velocity can be compared in the free flight and simulation data. The results are shown in Figures 3.5 through 3.7. There is a marked difference in the trend for the steady state velocity between the Vicon and simulation data. Figures 3.8 through 3.10 compare the force and moment values from the look-up table for the equilibrium points reached at steady state by the ornithopter in true flight and in simulation. As shown in the plots, the pitch moment is close to zero for the equilibrium points in simulation, as expected, but nonzero for the equilibrium points above 85 percent duty cycle in true flight. A zero pitch moment is important for passive stability, both in reality and simulation, so it is clear from the plots that some interaction is missing in the wind tunnel collected data that is present in reality and contributes to the steady increase in velocity at higher duty cycles.
CHAPTER 3. SIMULATION RESULTS
70 75 80 85 90 95 100
5
z [m]
4 3 2
65
1
5
10
15
x [m]
45
5
Time [s]
10
2
1.5
1
5
Time [s]
10
15
Figure 3.3: Simulated velocity of the ornithopter.
Angle of Attack [Deg]
70 75 80 85 90 95 100
15
Figure 3.2: Simulated pitch angle of the ornithopter. 80
2.5
Velocity [m/s]
50
35 0
20
Figure 3.1: Simulated paths of the ornithopter for various duty cycles.
0.5 0
55
40
0 −1 0
70 75 80 85 90 95 100
60
Pitch [Deg]
6
16
70 75 80 85 90 95 100
70
60
50
40
30 0
5
Time [s]
10
15
Figure 3.4: Simulated angle of attack of the ornithopter.
CHAPTER 3. SIMULATION RESULTS
2.5
Vicon Simulation
70
Vicon Simulation
65
Velocity [m/s]
Angle of Attack [Degrees]
75
17
60 55 50
2
1.5
45 40 35 70
75
80
85
90
Duty Cycle [Percent]
95
1 70
100
Figure 3.5: Comparison between Vicon (red) and simulation (blue) angle of attack.
Pitch [Degrees]
55
75
80
85
90
Duty Cycle [Percent]
95
Figure 3.6: Comparison between Vicon (red) and simulation (blue) velocity magnitude.
Vicon Simulation
50
45 70
75
80
100
85
90
Duty Cycle [Percent]
95
100
Figure 3.7: Comparison between Vicon (red) and simulation (blue) pitch angle.
CHAPTER 3. SIMULATION RESULTS
0.015
Vicon Simulation
0
−0.02
0.005 0 −0.005
−0.04
−0.06 70
Vicon Simulation
0.01
0.02
Vertical Force [N]
Forward Force [N]
0.04
18
−0.01
75
80
85
90
Duty Cycle [Percent]
95
−0.015 70
100
Figure 3.8: Comparison between Vicon (red) and simulation (blue) horizontal force at equilibrium.
Pitch Moment [N*mm]
2
75
80
85
90
Duty Cycle [Percent]
95
100
Figure 3.9: Comparison between Vicon (red) and simulation (blue) vertical force at equilibrium.
Vicon Simulation
1.5
1
0.5
0
−0.5 70
75
80
85
90
Duty Cycle [Percent]
95
100
Figure 3.10: Comparison between Vicon (red) and simulation (blue) pitch moment at equilibrium.
19
Chapter 4
Yaw Plane Simulation The yaw plane dynamics of the ornithopter are modeled using a one-dimensional model to reduce computational complexity. The model is used in a simulation framework that is used to determine the feasibility of the cooperative control experiment described in Section 1.2.
4.1
Yaw Dynamics
The turning dynamics of the ornithopter in the ”Internal Control Loop” block in Figure 1.5 are modeled using a one-dimensional model (Equation 1.2). To determine the constants K and R in the turning model, the step response of the ornithopter for an input of a clockwise 90 degree turn was determined. The response of the system was recorded using the attitude estimates calculated onboard the ornithopter using its gyroscope. The results of the experiment are shown in Figure 4.1. Using MATLAB, the low order model in Equation 1.2 was fitted to the step response. The resulting model is: 0.95 θ(s) = u(s) (4.1) 1 + 0.62s Additionally, the step response rise time was 1.4 seconds and the ornithopter flies at an average forward velocity of 1.2 m/s for these experiments, thus the estimated minimum turning radius of the robot is 1.07 m. The remainder of the parameters, such as the camera and transmission latencies and the external PID constants, are determined through methods detailed in [7]. This system layout is used to predict the success of the single-view tracking and guidance of the ornithopter to the window.
4.2
Yaw Model Validation
The turning model in conjunction with the total cooperative system model, is meant to simulate the cooperative control problem of guiding the ornithopter through a window of 0.66 m by 0.94 m dimensions. To validate the turning model defined in Section 4.1, several experiments were conducted to determine the feasible and infeasible initial condition locations for successful navigation to the window. To determine the success rates in various locations of the view plane, 80 experimental trials were conducted in which the ground station tracks the location of the ornithopter with the camera
CHAPTER 4. YAW PLANE SIMULATION
1.6
Yaw Rotation [rad]
1.4
20
Reference Yaw Measured Yaw Fitted Model
1.2 1 0.8 0.6 0.4 0.2 0 −0.2
−2
0
Time [sec]
2
Figure 4.1: The step response of the ornithopter to a clockwise 90 degree turn.
and attempts to guide it to the window by transmitting a desired heading to the ornithopter. For each trial, the ornithopter was launched by hand and oriented perpendicular to the window plane, although human error and the initial launch transient added some variation to this heading. Of the 80 trials, 20 were initiated from 0.8 m directly in front of the camera, and 60 were initiated in 5 cm increments from the window plane on the edges of the view space. The results of the experiments are in Figures 4.2 and 4.3. The success rate in front of the camera was 80 percent, and the success rate around the edges of the view plane was 50 percent. To validate the model, two experiments were conducted in simulation. The first simulation experiment was to calculate the geometrically infeasible region in which it is impossible for the ornithopter to reach the window due to constraints on the minimum turning radius of 1.07 m. This region is overlayed with the infeasible region calculated by the simulation. For each scenario, the ornithopter is assumed to begin with an initial heading facing the window plane. The results are shown in Figure 4.4, with each marker representing a starting location in a 5 cm grid in the camera view space. Although this result provides some information about the regions in which successful control is possible, it does not provide any explanation about why a higher success rate in the middle of the view space opposed to the edges of the view space was observed, even though both locations are technically feasible. To explain this result, a Monte Carlo simulation of the system was performed. The trajectory of the ornithopter was simulated using initial translational positions in a 10 cm grid in the view space, and the initial heading was randomly sampled from a uniform distribution of headings between -90 degrees (left, parallel to window plane) and +90 degrees (right, parallel
CHAPTER 4. YAW PLANE SIMULATION
Success Failure
4
3
y [meters]
y [meters]
4
21
2
3 2
1
1
0
0 −2
−1
0
x [meters]
1
2
Figure 4.2: Plot of camera field of view and experimental trials within the feasible region.
Success Failure
−2
−1
0
x [meters]
1
2
Figure 4.3: Plot of camera field-of-view and perimeter experimental trials overlayed.
to window plane). Forty trials were conducted in each location within the grid. The results of the experiments are illustrated in Figure 4.5. The circled region represents the initial location of the experimental trials shown in Figure 4.2.
CHAPTER 4. YAW PLANE SIMULATION
5
y [meters]
4
22
Feasible Region Geometric and Simulation Infeasible Region Simulation Infeasible Region
3 2 1 0 −2
−1
0
x [meters]
1
2
Figure 4.4: Plot of backwards reachable set for successful window traversal.
Figure 4.5: Probability of success, resulting from a Monte Carlo simulation of the reachable set model, for all starting locations.
23
Chapter 5
Conclusions and Analysis This report has presented two software tools implemented in MATLAB to simulate a 13.6 g ornithopter in the sagittal plane and horizontal plane. Each of these simulations were compared to true flight data collected using a Vicon motion capture system to determine their accuracy.
5.1
Sagittal Plane Simulation
The sagittal plane simulation computes the steady-state motion of the ornithopter with the duty cycle of the motor wing drive as an input. The motion is computed based upon wind tunnel aerodynamics adjusted using known true flight equilibrium points. The results of the simulation are shown in Figures 3.1 through 3.4, and the steady-state comparison between the simulation and true flight is shown in Figures 3.5 through 3.7. As indicated in the angle of attack comparison between the true flight and simulation data, the simulation closely follows the trend exhibited by the true flight of the ornithopter. The same cannot be said for the relationship between the velocities. The steady-state velocity follows the same trend for both the true and simulation flight within 0.3 m/s for duty cycles up to 85 percent. For duty cycles above 85 percent, however, the steady-state velocity decreases in simulation but increases in reality for increasing duty cycle. The limitation in the accuracy of the velocity of the ornithopter in simulation could be caused by several factors. The equilibrium shift used for the wind tunnel data only exactly defines the equilibrium point at 85 percent duty cycle. This fact is illustrated in Figures 3.8 through 3.10. The two points are very close for the forward and vertical forces and pitch moment for 85 percent duty cycle. For the duty cycles above 85 percent, the pitch moment is very close to zero in steady state for the simulation, but greater than 0.8 N·mm for the equilibrium points collected from the Vicon. A zero pitch moment is important for the passive stability in the simulation and the pitch moment is dependent on the angle of attack, duty cycle, and velocity. Since the duty cycles at these equilibrium points are equal for the Vicon and simulation and the angles of attack are close, the velocity measured using the Vicon does not correspond to a stable equilibrium point in the wind tunnel surfaces. This variation is expected, since the aerodynamic interactions caused by the pitching motion of the ornithopter in true flight are absent from the data used for the simulation. Additionally, the flapping forces and moments in still air and the aerodynamic forces and moments for still wings in moving air were measured separately and added together. This method causes
CHAPTER 5. CONCLUSIONS AND ANALYSIS
24
some error in estimating the affects of the forces and moments caused by the airflow when it is disturbed by the flapping wings.
5.2
Yaw Plane Simulation
The yaw model simulation computes the motion of the ornithopter in the horizontal plane with a desired heading as an input. The model encompasses both the internal PID controller that computes the necessary yaw motor input as well as the yaw dynamics of the ornithopter. The true flight results of the ornithopter in the cooperative control system corroborates the model for the turning dynamics. As shown in Figures 4.2 and 4.3, the true flight tests yielded an 80 percent success rate in front of the camera and a 50 percent success rate along the edges of the view space. While these locations in the view space are technically feasible, as supported by the geometric and simulation backwards reachable sets in Figure 4.3, the Monte Carlo experiments provide more insight about the reduced success rate along the edges of the view space. There are several factors that contribute to the lower success rate experienced on the edges of the view space. The monocular vision control scheme has no notion of depth, so the control algorithm provides the same desired heading to the ornithopter whether it is close to the camera, where small changes in heading are adequate, or far away from the camera, where larger changes in heading are necessary to navigate to the window. Another factor that contributes to the lower success rate is noise in the launch by hand of the ornithopter at the beginning of the experiment. Small changes in initial heading, or disturbances caused by a hand can cause the ornithopter to fly out of the view plane or begin flight in an erratic manner making it impossible to reach the window successfully. Both of these causes of error are captured in the Monte Carlo experiment with the randomly selected initial headings and the lack of depth considerations in the calculation of the control inputs.
5.3
Future Work
In the future, the sagittal plane and horizontal plane simulations will be combined into a three dimensional simulation for the ornithopter. In terms of possible uses, the simulation could be used to determine what sorts of flight maneuvers are possible for the ornithopter. This information could lead to model-inspired redesigns of the platform to increase maneuverability. Additionally, the development of this model allows for more sophisticated control schemes to be implemented on the ornithopter such as LQR and MPC to generate optimal control paths to a goal. Although the simulator is an accurate predictor of the flight path of the ornithopter, there are some differences between the true flight and simulated flight, especially in terms of steady state velocity. In the future, it would be useful to examine the airflow around the ornithopter with the wings flapping and without using Particle Image Velocimetry (PIV) to determine exactly how the air moves over the control surfaces. This experiment could provide some insight as to what types of interactions are missing in the wind tunnel data collection.
25
Bibliography [1]
S.S. Baek, F.L. Garcia Bermudez, and R.S. Fearing. “Flight Control for Target Seeking by 13 gram Ornithopter.” In: IEEE Int.l Conf. on Intelligent Robots and Systems. 2011, pp. 2674– 2681.
[2]
S.S. Baek and R.S. Fearing. “Flight forces and altitude regulation of 12 gram I-Bird”. In: IEEE Int.l Conf. on Biomedical Robotics and Biomechatronics. 2010, pp. 454–460.
[3]
B. Cheng and X. Deng. “Translational and Rotational Damping of Flapping Flight and Its Dynamics and Stability at Hovering”. In: IEEE Transactions on Robotics Vol. 27.5 (2011), pp. 849–864.
[4]
G. C. H. E. de Croon, K. M. E. de Clercq, R. Ruijsink, B. Remes, and C. de Wagter. “Design, aerodynamics, and vision-based control of the DelFly”. In: International Journal of Micro Air Vehicles Vol. 1.2 (June 2009), pp. 71–97.
[5]
JH Han, JY Lee, and DK Kim. “Ornithopter Modeling for Flight Simulation”. In: Intl. Conf. on Control, Automation and Systems. 2008, pp. 1773–1777.
[6]
D.G. Hull. Fundamentals of Airplane Flight Mechanics. 1st. Springer Publishing Company, Inc., 2007.
[7]
R. Julian, C. Rose, H. Hu, and R.S. Fearing. “Cooperative Control for Window Traversal with an Ornithopter MAV and Lightweight Ground Station”. In: 12th Intl. Conf. on Autonomous Agents and Multiagent Systems, 2013. Proceedings. 2013.
[8]
Z.A. Khan and S.K. Agrawal. “Control of Longitudinal Flight Dynamics of a Flapping-Wing Micro Air Vehicle Using Time-Averaged Model and Differential Flatness Based Controller”. In: American Control Conference. 2007, pp. 5284–5289.
[9]
S. M. LaValle. Planning Algorithms. Available at http://planning.cs.uiuc.edu/. Cambridge, U.K.: Cambridge University Press, 2006.
[10]
A. T. Pfeiffer, JS Lee, JH Han, and H. Baier. “Ornithopter Flight Simulation Based on Flexible Multi-Body Dynamics”. In: Journal of Bionic Engineering Vol. 7.1 (2010), pp. 102 –111.
[11]
L. Schenato, D. Campolo, and S. Sastry. “Controllability issues in flapping flight for biomimetic micro aerial vehicles (MAVs)”. In: 42nd IEEE Conf. on Decision and Control. Vol. 6. 2003, pp. 6441–6447.
Cameron Rose Ronald S. Fearing
Electrical Engineering and Computer Sciences University of California at Berkeley Technical Report No. UCB/EECS-2013-60 http://www.eecs.berkeley.edu/Pubs/TechRpts/2013/EECS-2013-60.html
May 14, 2013
Copyright © 2013, by the author(s). All rights reserved. Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, to republish, to post on servers or to redistribute to lists, requires prior specific permission.
ii
Contents 1 Simulation Environments 1.1 Aerodynamic and Dynamic Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Yaw Dynamics Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2 3 5
2 Sagittal Plane Experimental Data 2.1 Wind Tunnel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Free Flight . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Comparison of Approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7 7 8 8
3 Simulation Results 15 3.1 Free Flight Comparisons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4 Yaw Plane Simulation 19 4.1 Yaw Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.2 Yaw Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 5 Conclusions and Analysis 5.1 Sagittal Plane Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Yaw Plane Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23 23 24 24
Bibliography
25
2
Chapter 1
Simulation Environments The ornithopter simulation software package includes two modules: a sagittal plane model and a yaw plane model. Each module is implemented separately and models the H2 Bird hardware [7]. The robot and its control surfaces are shown in Figure 1. The H2 Bird is built around the Silverlit i-Bird RC flier power train and has clap-fling wings with a wingspan of 26.5 cm. For yaw and pitch control, it uses a tail propeller and a servo-controlled tail elevator. In terms of control and sensing, the H2 Bird uses an on-board ImageProc 2.4 controller that includes a 40 MIPS microprocessor, 6 DOF IMU, IEEE 802.15.4 radio, and motor drivers, powered by a 90 mAh lithium polymer battery [1].
Wings
Tail prop ImageProc Elevator servo Tail elevator
Figure 1.1: H2 Bird ornithopter MAV. The sagittal plane ornithopter simulation is modeled using rigid body dynamics with force and moment information collected from wind tunnel and free flight tests as inputs. The top-level layout of the simulator is illustrated in Figure 1.2 and is implemented in MATLAB. The duty cycle of the motor wing drive is the input, and the state of the ornithopter is updated at 1 kHz. The yaw plane simulation is also implemented in MATLAB. The yaw dynamics are determined using system identification on the step response of the ornithopter to a 90 degree right turn. The simulation computes the yaw motion of the ornithopter with a desired heading as the input.
CHAPTER 1. SIMULATION ENVIRONMENTS
3
Body Parameters
Aerodynamics Duty Cycle Controller
Body Dynamics Forces, Moments
State(t)
State(∆t)
∫ 1kHz
Figure 1.2: Simulator block diagram.
1.1
Aerodynamic and Dynamic Models
The aerodynamics of flapping winged flight are nonlinear and complex. The flapping wings disturb the airflow and create vortices that interact with the wings and tail of the ornithopter to augment the thrust forces produced by the flapping of the wings. In some cases, simply the thrust provided by the wings is not sufficient by itself to keep the robot aloft, although it still flies freely [2]. Several approaches to modeling the flight paths of ornithopters have been used, but they often involve variations of blade element theory or other methods to model the dynamics of the wings in an effort to predict the flight path [3, 5, 10]. In my model, I simplify the process by using averages of measured forces and moments to predict the motion of the ornithopter. The method of using averaged data has been previously shown to be a valid approximation [8, 11]. This strategy enables a simulation that is less computationally complex, which is helpful for simulating multiple ornithopters or working with more intricate control strategies. The previous control strategies used for earlier iterations of the ornithopter were based upon Proportional-Integral-Differential (PID) control schemes for target tracking or height regulation [1, 2]. Information about the aerodynamic interactions of the ornithopter can produce a more robust model that can be used in more sophisticated control schemes such as Linear-Quadratic Regulator (LQR) or Model Predictive Control (MPC). The aerodynamic forces and moments used in the simulation are collected using a force and torque sensor affixed to the bottom of the ornithopter in a wind tunnel. The forces and moments were measured at the center of mass of the entire robot to attempt to capture the interactions between the wings and the rest of the robot without separately modeling each interaction. These measurements were corrected by free flight data collected using a Vicon motion capture system. The modification of the fixed wind tunnel data is necessary, as the raw fixed data are not sufficient to capture the complex nonlinear dynamics of the ornithopter [4]. The pitching motion that the ornithopter exhibits in free flight is important to the flight dynamics, although this could not be captured in the wind tunnel. The horizontal and vertical forces and the pitching moment are used as inputs into the body
CHAPTER 1. SIMULATION ENVIRONMENTS
4
dynamics in the simulation. Each force and moment is stored at set data points and interpolation is used to calculate the points in between. The forces and moments are all functions of the set duty cycle, angle of attack, and velocity magnitude of the ornithopter at a given point in time. The body dynamics of the ornithopter are treated as a rigid body. The kinematic equations of motion are as follows [6]: x˙ = u cos θ + w sin θ z˙ = u sin θ − w cos θ 1 u˙ = −wQ + (T − mg sin θ) m 1 w˙ = uQ + (−L + mg cos θ) m θ˙ = Q M − damping Q˙ = Iyy
(1.1)
The state vector has six elements: x is the horizontal position in world coordinates, z is the vertical position in world coordinates, u is the forward velocity in body coordinates, w is the vertical velocity in body coordinates, θ is the pitch, and Q is the angular velocity. Additionally, T is the thrust and L is the lift in body coordinates, and M is the pitch moment. The forces are measured by the sensor in body coordinates, so the velocities are in body coordinates for ease of use. In terms of the body parameters, the pitch moment of inertia, Iyy , is 3.0×10−5 kg·m2 and was computed using SolidWorks, which computes the mass distribution of the model of the ornithopter and returns it in the form of an inertial matrix. The mass of the ornithopter, m, is 13.6 g, and was measured using a scale. The body diagram is shown in Figure 1.3.
zw
xb
L u T
𝛂
𝚹
V mg
w zb
Figure 1.3: The rigid body diagram for the model.
xw
CHAPTER 1. SIMULATION ENVIRONMENTS
1.2
5
Yaw Dynamics Model
The turning dynamic model is part of a larger model used to simulate the motion of the ornithopter in the horizontal plane in a cooperative control scenario. This model is used to predict the feasibility of achieving a goal successfully for a given initial pose. In the scenario, a ground station is used to provide guidance to the ornithopter by transmitting desired headings to a window calculated from the visual input provided by an off-the-shelf USB web camera. The layout of the experiment is depicted in Figure 1.4. The total cooperative system software layout is shown in Figure 1.5.
H2Bird
0.66 m
Window
0.94 m
Camera 3.66 m Ground Station
Figure 1.4: Conceptual sketch of experimental setup. The yaw dynamics of the ornithopter are modeled in the horizontal plane using a one dimensional model to limit the complexity. This format also matches the physical turning implementation on the actual ornithopter. In terms of control, the ornithopter is sent a desired heading and it controls to that heading using an on-board PID controller. The turning dynamics and the internal PID controller are modeled with a single one-dimensional model: θ(s) =
K u(s) 1 + Rs
(1.2)
where K and R are constants determined through system identification and the input, u(s), is the desired heading and the output, θ(s), is the physical heading of the ornithopter. A form of Dubin’s car model is used to model the translational dynamics in the horizontal
CHAPTER 1. SIMULATION ENVIRONMENTS
6
plane [9]: x˙ = us sin θ y˙ = us cos θ θ˙ = u
(1.3)
u ∈ U = [− tan φmax , tan φmax ] In the car model, φmax is the maximum steering angle and us is the forward speed of 1.2 m/s. Here, the angular position θ is an input rather than a state. The overall motion of the system can be modeled as shown in Equation 1.3, where θ is the time-domain output of Equation 1.2, θ(t).
External Control Loop – 10Hz Window position
Desired Internal Control Loop – 300Hz Heading PID P-D Ornithopter Controller Controller Dynamics
-
Translational Position Pose Estimation
Angular position
Attitude Estimation
Camera
Figure 1.5: Overall block diagram of the cooperative control system.
Pose
7
Chapter 2
Sagittal Plane Experimental Data The experimental data necessary to determine the aerodynamic properties of the ornithopter are comprised of both wind tunnel force-torque data and free flight data collected from a Vicon motion capture system.
2.1
Wind Tunnel
The wind tunnel data were collected over a series of trials using an ATI Nano17 force-torque sensor. The ornithopter was affixed to an acrylic mount attached to the sensor and facing into the wind, as shown in Figure 2.1. The ornithopter was attached at approximately its center of mass, although this point fluctuates as the wings open and close. This fluctuation is minimal, however, so it was discounted. Each trial consisted of 3 seconds of data collected for wind speeds between 0 and 3 m/s in 0.5 m/s increments and longitudinal body axis angles relative to the horizontal, or angles of attack, between -70 and +90 degrees in 10 degree increments. Each of these trials was conducted with the wings of the robot closed and stationary. Each trial was averaged over the 3 seconds of collection, and only the thrust and lift forces in body coordinates and the pitch moment were stored. The results of the experiments are in Figures 2.2 through 2.4. Each force is expressed in the world reference frame with the weight of the robot absent from the measurements. In addition to the stationary wing experiments, data were collected in still air for duty cycles between 50 percent and 100 percent in 10 percent increments using the same wind tunnel setup as the previous experiments. This range corresponds to wing flap speeds between 12 Hz and 20 Hz. These experiments required some damping around the connection of the robot to the acrylic mount with 4 mm of foam, as the oscillations from the wings moving created significant noise in the pitch moment signal. The data collected from the experiments were averaged over 5 seconds, and the thrust and lift forces in body coordinates and the pitch moment were stored. The results of the experiments are in Figures 2.5 and 2.6. Both of these results were combined to form a 4 dimensional data set dependent upon duty cycle, angle of attack, and wind speed that outputs the forward and vertical body frame forces and pitch moment. This data set is used as a look-up table for the instantaneous forces and moments for a given pose in the simulation.
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
8
L T
Figure 2.1: H2 Bird mounted to sensor
2.2
Free Flight
The free flight data was captured over variable amounts of time for circular paths in an empty room for duty cycles between 70 and 100 percent in 5 percent increments. For each trial, the initial launch transient was cut from the data, because only the behavior around equilibrium was of interest. For each trial, the position, orientation, angle of attack, and velocity magnitude and direction of the ornithopter were computed. Data from a typical trial are shown in Figures 2.7 and 2.8. The relationship between duty cycle and velocity, angle of attack, and pitch angle are depicted in Figures 2.9 through 2.11. The trends and locations of the equilibrium points in terms of the angle of attack and velocity are used in conjunction with the wind tunnel data to determine an accurate picture of the flight characteristics for simulation.
2.3
Comparison of Approaches
The two approaches to determining the flight forces of the ornithopter both have their merits and limitations. The wind tunnel data give a picture of the aerodynamics of the ornithopter over a broad range of wind conditions. Although all of the motions that contribute to the flight of the ornithopter are not present in this fixed measurement, the wind tunnel data still provides an indication of the relationship of the forces and moments to changes in wind speed and angle of attack. While these relationships can be used, the values of the individual data points may diverge, as some of the aerodynamic interactions are not present. The free flight data provide information about the exact values of the angle of attack and flight speed at equilibrium for given duty cycles. This information can be used to correct the values of the forces and moments produced by the wind tunnel experiments, but the relationships of the wind tunnel data can still be used. One limitation of the free flight data is that it only gives an indication of the flight behaviors at equilibrium, which is a narrow snapshot of the flight dynamics.
0.2
Vertical Force [N]
Horizontal Force [N]
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
0.1 0 −0.1 −0.2 0 1 2 3
0.1 0 −0.1 −0.2 0 1 2 3
AoA [Degrees]
Wind Speed [m/s]
Wind Speed [m/s]
Figure 2.2: Aerodynamic horizontal force surface in world coordinates. Pitch Moment [N*mm]
0.2
50
0
−50
9
−50
0
50
AoA [Degrees]
Figure 2.3: Aerodynamic vertical force surface in world coordinates.
5
0
−5
−10 0 1 2 3
Wind Speed [m/s]
−50
0
50
AoA [Degrees]
Figure 2.4: Aerodynamic pitch moment surface.
Additionally, the flapping wings produce noise in the Vicon tracking system, so it is difficult to capture valid data. To combine the two sets of data into one composite data set to use with the simulator, several adjustments were made to the wind tunnel data that were inspired by the free flight data. One of these adjustments was to the equilibrium points indicated by the wind tunnel data, which did not match with the values of the equilibrium points indicated by the free flight data. This issue was expected, since the ornithopter is fixed in the wind tunnel. To account for the difference, the zero force and moment point of the wind tunnel data was shifted to the point supplied by the free flight data. This change was implemented by subtracting, from the entire wind tunnel data set, the values provided by the wind tunnel at the free flight provided equilibrium point: 85 percent duty cycle, 1.3 m/s wind speed, and 53 degrees angle of attack. Therefore, the values of the data were shifted, but all of the relationships remain the same. Another discrepancy between the wind tunnel data set and the free flight data, is that in the wind tunnel data, the relationship between the pitch moment and the wind speed was incorrect. As shown in Figure 2.12, the measured pitch moment decreases as the wind speed increases. If this were the case, the robot would pitch down as the speed increases, which would cause a further speed increase. This is counterintuitive, and also contradictory of the actual flight behavior. This
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
0.1
−1
Thrust Force Lift Force
Pitch Moment [N*mm]
0.12
Forces [N]
0.08 0.06 0.04 0.02 0 −0.02 50
60
70
10
80
90
Duty Cycle [Percent]
100
Figure 2.5: Thrust (blue) and lift (green) forces as a function of the duty cycle.
−1.1 −1.2 −1.3 −1.4 −1.5 50
60
70
80
Duty Cycle [Percent]
90
100
Figure 2.6: Pitch moment as a function of the duty cycle.
deficiency could be caused by some aerodynamic interaction that is not captured in the wind tunnel that causes the robot to pitch up in free flight. To fix this problem, a further additive term was added to the pitch moment dependent on the wind speed: 7 ∗ windspeed . This term is meant to 3.0 behave like an elevator would on the end of the tail. This approximation is overly simplistic, as, in actuality, the moment provided by an elevator would also be dependent on angle of attack. However, this term fixes the relationship between wind speed and pitch moment, while preserving the pre-existing correct relationship between angle of attack and pitch moment. The progressive effects of these changes to the pitch moment surface are illustrated in Figures 2.14 through 2.16. The shift illustrated in Figure 2.15 is analogous to a similar shift in the forward and vertical forces based upon the equilibrium points determined in the free flight experiments. Figure 2.13 gives a closer look at the adjusted relationship between wind speed and the pitch moment for a particular angle of attack.
Pitch [Deg]
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
100 50 0 0
Z [m]
1
0.5
2
0 −2 1
−1 0
0 1 −1
X [m]
Y [m]
Figure 2.7: The Vicon tracked path for 80 percent duty cycle with velocity vector (blue) and angular orientation vector (green) indicated.
1
2
3
4
5
2
3
4
5
2
3
4
5
Time [s]
100 50
Velocity [m/sec]
AoA [Deg]
2
1.5
11
0 0
1
Time [s]
4 2 0 0
1
Time [s]
Figure 2.8: The Vicon tracked pitch angle, angle of attack, and velocity magnitude for 80 percent duty cycle.
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
70
Angle of Attack [Degrees]
54
Pitch [Degrees]
52 50 48 46 44 42 60
12
70
80
90
100
60 55 50 45 40 35 30 60
110
Duty Cycle [Percent]
65
Figure 2.9: Free flight pitch angle as a function of the duty cycle.
70
80
90
Velocity [m/s]
2
1.5
70
110
Figure 2.10: Free flight angle of attack as a function of the duty cycle.
2.5
1 60
100
Duty Cycle [Percent]
80
90
100
Duty Cycle [Percent]
110
Figure 2.11: Free flight velocity magnitude as a function of the duty cycle.
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
2
0.5
Pitch Moment [N*mm]
Pitch Moment [N*mm]
1
0 −0.5 −1 −1.5 −2 −2.5 0
13
0.5
1
1.5
2
Wind Speed [m/s]
2.5
3
Figure 2.12: Original pitch moment for 85 percent duty cycle and 53 degrees angle of attack as a function of wind speed.
1
0
−1
−2
−3 0
0.5
1
1.5
2
Wind Speed [m/s]
2.5
3
Figure 2.13: Adjusted pitch moment for 85 percent duty cycle and 53 degrees angle of attack as a function of wind speed.
Pitch Moment [N*mm]
Pitch Moment [N*mm]
CHAPTER 2. SAGITTAL PLANE EXPERIMENTAL DATA
5
0
−5 0 1 2 3
−5 0 2 3
Wind Speed [m/s]
Figure 2.14: Initial total pitch moment surface.
Pitch Moment [N*mm]
0
1
AoA [Degrees]
Wind Speed [m/s]
5
50
0
−50
14
0
−5
1
Wind Speed [m/s]
2
0
50
AoA [Degrees]
Figure 2.15: Equilibrium shifted total pitch moment surface.
5
0
−50
3
−50
0
50
AoA [Degrees]
Figure 2.16: Final adjusted total pitch moment surface.
15
Chapter 3
Simulation Results To determine the accuracy of the simulation, open loop flight tests of the sagittal plane model were conducted. The simulation results were compared to analogous free flight tests.
3.1
Free Flight Comparisons
For the wind tunnel inspired model, the motion of the ornithopter was simulated for duty cycles between 70 and 100 percent in 5 percent increments. The initial conditions chosen were 0.85 m/s forward velocity in the horizontal plane and a 40 degree pitch angle. For each trajectory, the path of the flight and also the angle of attack, velocity magnitude, and pitch angle were plotted, as shown in Figures 3.1 through 3.4. To determine the overall steady-state accuracy of the simulation in comparison to the real flight data, the simulation results were compared to the previously collected flight data. As mentioned in the previous chapter, this flight data was collected for duty cycles between 70 and 100 percent in 5 percent increments for circular paths. The paths are circular because the room was not big enough to collect multiple seconds of data after the initial launch transient before equilibrium flight is reached. Because of this difference, the steady state values are not exact for straight forward flight. Since the flapping wings are not modeled in the simulation, the oscillatory behavior of the true flight path is not captured. However, the steady state average angle of attack, pitch, and velocity can be compared in the free flight and simulation data. The results are shown in Figures 3.5 through 3.7. There is a marked difference in the trend for the steady state velocity between the Vicon and simulation data. Figures 3.8 through 3.10 compare the force and moment values from the look-up table for the equilibrium points reached at steady state by the ornithopter in true flight and in simulation. As shown in the plots, the pitch moment is close to zero for the equilibrium points in simulation, as expected, but nonzero for the equilibrium points above 85 percent duty cycle in true flight. A zero pitch moment is important for passive stability, both in reality and simulation, so it is clear from the plots that some interaction is missing in the wind tunnel collected data that is present in reality and contributes to the steady increase in velocity at higher duty cycles.
CHAPTER 3. SIMULATION RESULTS
70 75 80 85 90 95 100
5
z [m]
4 3 2
65
1
5
10
15
x [m]
45
5
Time [s]
10
2
1.5
1
5
Time [s]
10
15
Figure 3.3: Simulated velocity of the ornithopter.
Angle of Attack [Deg]
70 75 80 85 90 95 100
15
Figure 3.2: Simulated pitch angle of the ornithopter. 80
2.5
Velocity [m/s]
50
35 0
20
Figure 3.1: Simulated paths of the ornithopter for various duty cycles.
0.5 0
55
40
0 −1 0
70 75 80 85 90 95 100
60
Pitch [Deg]
6
16
70 75 80 85 90 95 100
70
60
50
40
30 0
5
Time [s]
10
15
Figure 3.4: Simulated angle of attack of the ornithopter.
CHAPTER 3. SIMULATION RESULTS
2.5
Vicon Simulation
70
Vicon Simulation
65
Velocity [m/s]
Angle of Attack [Degrees]
75
17
60 55 50
2
1.5
45 40 35 70
75
80
85
90
Duty Cycle [Percent]
95
1 70
100
Figure 3.5: Comparison between Vicon (red) and simulation (blue) angle of attack.
Pitch [Degrees]
55
75
80
85
90
Duty Cycle [Percent]
95
Figure 3.6: Comparison between Vicon (red) and simulation (blue) velocity magnitude.
Vicon Simulation
50
45 70
75
80
100
85
90
Duty Cycle [Percent]
95
100
Figure 3.7: Comparison between Vicon (red) and simulation (blue) pitch angle.
CHAPTER 3. SIMULATION RESULTS
0.015
Vicon Simulation
0
−0.02
0.005 0 −0.005
−0.04
−0.06 70
Vicon Simulation
0.01
0.02
Vertical Force [N]
Forward Force [N]
0.04
18
−0.01
75
80
85
90
Duty Cycle [Percent]
95
−0.015 70
100
Figure 3.8: Comparison between Vicon (red) and simulation (blue) horizontal force at equilibrium.
Pitch Moment [N*mm]
2
75
80
85
90
Duty Cycle [Percent]
95
100
Figure 3.9: Comparison between Vicon (red) and simulation (blue) vertical force at equilibrium.
Vicon Simulation
1.5
1
0.5
0
−0.5 70
75
80
85
90
Duty Cycle [Percent]
95
100
Figure 3.10: Comparison between Vicon (red) and simulation (blue) pitch moment at equilibrium.
19
Chapter 4
Yaw Plane Simulation The yaw plane dynamics of the ornithopter are modeled using a one-dimensional model to reduce computational complexity. The model is used in a simulation framework that is used to determine the feasibility of the cooperative control experiment described in Section 1.2.
4.1
Yaw Dynamics
The turning dynamics of the ornithopter in the ”Internal Control Loop” block in Figure 1.5 are modeled using a one-dimensional model (Equation 1.2). To determine the constants K and R in the turning model, the step response of the ornithopter for an input of a clockwise 90 degree turn was determined. The response of the system was recorded using the attitude estimates calculated onboard the ornithopter using its gyroscope. The results of the experiment are shown in Figure 4.1. Using MATLAB, the low order model in Equation 1.2 was fitted to the step response. The resulting model is: 0.95 θ(s) = u(s) (4.1) 1 + 0.62s Additionally, the step response rise time was 1.4 seconds and the ornithopter flies at an average forward velocity of 1.2 m/s for these experiments, thus the estimated minimum turning radius of the robot is 1.07 m. The remainder of the parameters, such as the camera and transmission latencies and the external PID constants, are determined through methods detailed in [7]. This system layout is used to predict the success of the single-view tracking and guidance of the ornithopter to the window.
4.2
Yaw Model Validation
The turning model in conjunction with the total cooperative system model, is meant to simulate the cooperative control problem of guiding the ornithopter through a window of 0.66 m by 0.94 m dimensions. To validate the turning model defined in Section 4.1, several experiments were conducted to determine the feasible and infeasible initial condition locations for successful navigation to the window. To determine the success rates in various locations of the view plane, 80 experimental trials were conducted in which the ground station tracks the location of the ornithopter with the camera
CHAPTER 4. YAW PLANE SIMULATION
1.6
Yaw Rotation [rad]
1.4
20
Reference Yaw Measured Yaw Fitted Model
1.2 1 0.8 0.6 0.4 0.2 0 −0.2
−2
0
Time [sec]
2
Figure 4.1: The step response of the ornithopter to a clockwise 90 degree turn.
and attempts to guide it to the window by transmitting a desired heading to the ornithopter. For each trial, the ornithopter was launched by hand and oriented perpendicular to the window plane, although human error and the initial launch transient added some variation to this heading. Of the 80 trials, 20 were initiated from 0.8 m directly in front of the camera, and 60 were initiated in 5 cm increments from the window plane on the edges of the view space. The results of the experiments are in Figures 4.2 and 4.3. The success rate in front of the camera was 80 percent, and the success rate around the edges of the view plane was 50 percent. To validate the model, two experiments were conducted in simulation. The first simulation experiment was to calculate the geometrically infeasible region in which it is impossible for the ornithopter to reach the window due to constraints on the minimum turning radius of 1.07 m. This region is overlayed with the infeasible region calculated by the simulation. For each scenario, the ornithopter is assumed to begin with an initial heading facing the window plane. The results are shown in Figure 4.4, with each marker representing a starting location in a 5 cm grid in the camera view space. Although this result provides some information about the regions in which successful control is possible, it does not provide any explanation about why a higher success rate in the middle of the view space opposed to the edges of the view space was observed, even though both locations are technically feasible. To explain this result, a Monte Carlo simulation of the system was performed. The trajectory of the ornithopter was simulated using initial translational positions in a 10 cm grid in the view space, and the initial heading was randomly sampled from a uniform distribution of headings between -90 degrees (left, parallel to window plane) and +90 degrees (right, parallel
CHAPTER 4. YAW PLANE SIMULATION
Success Failure
4
3
y [meters]
y [meters]
4
21
2
3 2
1
1
0
0 −2
−1
0
x [meters]
1
2
Figure 4.2: Plot of camera field of view and experimental trials within the feasible region.
Success Failure
−2
−1
0
x [meters]
1
2
Figure 4.3: Plot of camera field-of-view and perimeter experimental trials overlayed.
to window plane). Forty trials were conducted in each location within the grid. The results of the experiments are illustrated in Figure 4.5. The circled region represents the initial location of the experimental trials shown in Figure 4.2.
CHAPTER 4. YAW PLANE SIMULATION
5
y [meters]
4
22
Feasible Region Geometric and Simulation Infeasible Region Simulation Infeasible Region
3 2 1 0 −2
−1
0
x [meters]
1
2
Figure 4.4: Plot of backwards reachable set for successful window traversal.
Figure 4.5: Probability of success, resulting from a Monte Carlo simulation of the reachable set model, for all starting locations.
23
Chapter 5
Conclusions and Analysis This report has presented two software tools implemented in MATLAB to simulate a 13.6 g ornithopter in the sagittal plane and horizontal plane. Each of these simulations were compared to true flight data collected using a Vicon motion capture system to determine their accuracy.
5.1
Sagittal Plane Simulation
The sagittal plane simulation computes the steady-state motion of the ornithopter with the duty cycle of the motor wing drive as an input. The motion is computed based upon wind tunnel aerodynamics adjusted using known true flight equilibrium points. The results of the simulation are shown in Figures 3.1 through 3.4, and the steady-state comparison between the simulation and true flight is shown in Figures 3.5 through 3.7. As indicated in the angle of attack comparison between the true flight and simulation data, the simulation closely follows the trend exhibited by the true flight of the ornithopter. The same cannot be said for the relationship between the velocities. The steady-state velocity follows the same trend for both the true and simulation flight within 0.3 m/s for duty cycles up to 85 percent. For duty cycles above 85 percent, however, the steady-state velocity decreases in simulation but increases in reality for increasing duty cycle. The limitation in the accuracy of the velocity of the ornithopter in simulation could be caused by several factors. The equilibrium shift used for the wind tunnel data only exactly defines the equilibrium point at 85 percent duty cycle. This fact is illustrated in Figures 3.8 through 3.10. The two points are very close for the forward and vertical forces and pitch moment for 85 percent duty cycle. For the duty cycles above 85 percent, the pitch moment is very close to zero in steady state for the simulation, but greater than 0.8 N·mm for the equilibrium points collected from the Vicon. A zero pitch moment is important for the passive stability in the simulation and the pitch moment is dependent on the angle of attack, duty cycle, and velocity. Since the duty cycles at these equilibrium points are equal for the Vicon and simulation and the angles of attack are close, the velocity measured using the Vicon does not correspond to a stable equilibrium point in the wind tunnel surfaces. This variation is expected, since the aerodynamic interactions caused by the pitching motion of the ornithopter in true flight are absent from the data used for the simulation. Additionally, the flapping forces and moments in still air and the aerodynamic forces and moments for still wings in moving air were measured separately and added together. This method causes
CHAPTER 5. CONCLUSIONS AND ANALYSIS
24
some error in estimating the affects of the forces and moments caused by the airflow when it is disturbed by the flapping wings.
5.2
Yaw Plane Simulation
The yaw model simulation computes the motion of the ornithopter in the horizontal plane with a desired heading as an input. The model encompasses both the internal PID controller that computes the necessary yaw motor input as well as the yaw dynamics of the ornithopter. The true flight results of the ornithopter in the cooperative control system corroborates the model for the turning dynamics. As shown in Figures 4.2 and 4.3, the true flight tests yielded an 80 percent success rate in front of the camera and a 50 percent success rate along the edges of the view space. While these locations in the view space are technically feasible, as supported by the geometric and simulation backwards reachable sets in Figure 4.3, the Monte Carlo experiments provide more insight about the reduced success rate along the edges of the view space. There are several factors that contribute to the lower success rate experienced on the edges of the view space. The monocular vision control scheme has no notion of depth, so the control algorithm provides the same desired heading to the ornithopter whether it is close to the camera, where small changes in heading are adequate, or far away from the camera, where larger changes in heading are necessary to navigate to the window. Another factor that contributes to the lower success rate is noise in the launch by hand of the ornithopter at the beginning of the experiment. Small changes in initial heading, or disturbances caused by a hand can cause the ornithopter to fly out of the view plane or begin flight in an erratic manner making it impossible to reach the window successfully. Both of these causes of error are captured in the Monte Carlo experiment with the randomly selected initial headings and the lack of depth considerations in the calculation of the control inputs.
5.3
Future Work
In the future, the sagittal plane and horizontal plane simulations will be combined into a three dimensional simulation for the ornithopter. In terms of possible uses, the simulation could be used to determine what sorts of flight maneuvers are possible for the ornithopter. This information could lead to model-inspired redesigns of the platform to increase maneuverability. Additionally, the development of this model allows for more sophisticated control schemes to be implemented on the ornithopter such as LQR and MPC to generate optimal control paths to a goal. Although the simulator is an accurate predictor of the flight path of the ornithopter, there are some differences between the true flight and simulated flight, especially in terms of steady state velocity. In the future, it would be useful to examine the airflow around the ornithopter with the wings flapping and without using Particle Image Velocimetry (PIV) to determine exactly how the air moves over the control surfaces. This experiment could provide some insight as to what types of interactions are missing in the wind tunnel data collection.
25
Bibliography [1]
S.S. Baek, F.L. Garcia Bermudez, and R.S. Fearing. “Flight Control for Target Seeking by 13 gram Ornithopter.” In: IEEE Int.l Conf. on Intelligent Robots and Systems. 2011, pp. 2674– 2681.
[2]
S.S. Baek and R.S. Fearing. “Flight forces and altitude regulation of 12 gram I-Bird”. In: IEEE Int.l Conf. on Biomedical Robotics and Biomechatronics. 2010, pp. 454–460.
[3]
B. Cheng and X. Deng. “Translational and Rotational Damping of Flapping Flight and Its Dynamics and Stability at Hovering”. In: IEEE Transactions on Robotics Vol. 27.5 (2011), pp. 849–864.
[4]
G. C. H. E. de Croon, K. M. E. de Clercq, R. Ruijsink, B. Remes, and C. de Wagter. “Design, aerodynamics, and vision-based control of the DelFly”. In: International Journal of Micro Air Vehicles Vol. 1.2 (June 2009), pp. 71–97.
[5]
JH Han, JY Lee, and DK Kim. “Ornithopter Modeling for Flight Simulation”. In: Intl. Conf. on Control, Automation and Systems. 2008, pp. 1773–1777.
[6]
D.G. Hull. Fundamentals of Airplane Flight Mechanics. 1st. Springer Publishing Company, Inc., 2007.
[7]
R. Julian, C. Rose, H. Hu, and R.S. Fearing. “Cooperative Control for Window Traversal with an Ornithopter MAV and Lightweight Ground Station”. In: 12th Intl. Conf. on Autonomous Agents and Multiagent Systems, 2013. Proceedings. 2013.
[8]
Z.A. Khan and S.K. Agrawal. “Control of Longitudinal Flight Dynamics of a Flapping-Wing Micro Air Vehicle Using Time-Averaged Model and Differential Flatness Based Controller”. In: American Control Conference. 2007, pp. 5284–5289.
[9]
S. M. LaValle. Planning Algorithms. Available at http://planning.cs.uiuc.edu/. Cambridge, U.K.: Cambridge University Press, 2006.
[10]
A. T. Pfeiffer, JS Lee, JH Han, and H. Baier. “Ornithopter Flight Simulation Based on Flexible Multi-Body Dynamics”. In: Journal of Bionic Engineering Vol. 7.1 (2010), pp. 102 –111.
[11]
L. Schenato, D. Campolo, and S. Sastry. “Controllability issues in flapping flight for biomimetic micro aerial vehicles (MAVs)”. In: 42nd IEEE Conf. on Decision and Control. Vol. 6. 2003, pp. 6441–6447.
