Modeling & optimization of projectile motion and simulation using FlightGear

Contents

• Team
• Motivation
• Proposal/Requirements
• Design/Models
  • Dynamics models
  • Optimization algorithms
  • FlightGear simulations
• Implementation
• Experiments & Results
• Conclusions
• References
• Presentation

Team

• Elijah Johnson
• elijah.johnson@mail.mcgill.ca

Motivation

The goal of this project was to further the study of projectile motion. Essentially, the well-known 2 Degree-Of-Freedom (DOF) point-mass system was extended to gain more insight into the effects of wind resistance/drag forces, lift forces and the Coriolis effect on a body.

To aid in visualizing the model simulations, output was sent to FlightGear, a powerful flight simulation software package.

Finally, the models were used to solve a simple ballistics problem: determining the initial conditions necessary to hit a specified target.

Proposal/Requirements

The equations governing the motion of a projectile were to be modeled using causal block diagrams in MATLAB's Simulink software. MATLAB's scripting environment was to be used to implement the optimization algorithms.

The modeling data was then to be fed into FlightGear to visualize the data.

Design/Models

Dynamics Models

• 2 Degree Of Freedom Systems
  
  
  
  Free Body Diagram of 2DOF syst.
  
  The dynamics equation for a 2 DOF system involves a simple point mass, depicted in the figure, an initial force Fi at angle theta, a "wind" force, and the gravitational force. The effect of wind on the projectile's motion was simplified to a force.
  
  The equations governing the body's motion are as follows:
  
  
  
  To facilitate modeling using Causal Block Diagrams in MATLAB, the equations were converted to their Laplace representations:
  
  
  
  These equations represent the solution to the dynamics equations.
• 3 Degree Of Freedom Systems
  
  
  
  Figure illustrating 3DOF free body diagram
  
  
  The dynamics equations in this case are very similar to the previous 2DOF case, and are expressed here in the Laplace-domain:
  
  
  
  Induced drag, i.e. drag on the projectile due to the body's own movement, was then added to the equations. It is given by
  
  
  
  where Cd is the drag coefficient, A is the projected area of the body that induces drag, rho is the density of air and v is the magnitude of the velocity of the body. This resulted in the following equation of motion:
  
  
  
  The normalized term represents the direction of the velocity vector. To compensate for the effect of wind, the induced drag term has as one of its factors the square of the magnitude of the relative velocity of the projectile w.r.t. wind velocity.
  
  Note that the equation is left in the time-domain representation. Since the velocity term is squared, the equation would have to be linearized in order to allow classical Laplace-domain analysis. Instead, a non-linear approach was implemented using the CBD. See the implementation section for details.
  
  Subsequently, Coriolis effects were added to the model. It is given by:
  
  a_c = where
  
  where omega is the angular acceleration of the rotating body, in this case earth, v is the velocity vector of the body in its own frame, and l is the latitude of the body on earth.
  
  
  
  Figure illustrating relation between
  local coordinate system and the rotating frame
  
  The resulting equation is as follows:
  
  
• 4 Degree Of Freedom System
  
  A pseudo-4DOF system was finally modeled. The time variation of theta was measured from the 3DOF model (simply the pitch given by the body's velocity vector) and used to calculate a lift force.
  
  Lift force was modeled using the equation:
  
  &
  
  where C_l is the lift coefficient, A_w is the wing area that generates the lift force, rho is the density of air and v is the magnitude of the velocity of the body. It was then directed in the vertical direction, effectively ignoring the lift-induced-drag (see lift force diagram).
  
  
  
  Lift force
  
  This resulted in the following equation of motion:
  
  
  where
  
  

Optimization algorithms

Given the models that plot the trajectory of a projectile given a set of model parameters, determining those necessary for the projectile to reach a specified target was achieved by using a root-finding method. Given the initial impulse force and direction, an script was written to measure the distance between the projectile and target after running a simulation, and this distance, or error, was used as input into the bisection root-finding method. The algorithm assumes a root, i.e. a point in the function where f(x) = 0, is between two supplied intervals, and iteratively halves the search interval until the root is found or the search space is small (see figure). In this fashion, the optimal initial direction was able to be found.

Figure illustrating the bisection method halving
the search space until the root/zero crossing is found

Equipped with an algorithm giving the angle needed to hit a target given a certain initial force, the issue of minimizing this force was investigated. If one were to plot the distribution of force vs. angle required to hit a target, an inverse-square distribution would result (see figures). A minimization algorithm was therefore used to determine the minimum force in the distribution, specifically MATLAB's built-in fminbnd function. It determines a function's local minimum given a search interval. It's algorithm is "based on golden section search and parabolic interpolation" (see http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fminbnd.html for more details).

Figure depicts distribution of force vectors
capable of hitting a target.

FlightGear simulations

In order to simulate projectile trajectories in flightGear, the typical XYZ position specified in a body frame must be translated to geodetic coordinates. This function was accomplished via the use of MATLAB's Flat Earth to LLA block. See here for more details: http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/flatearthtolla.html

Implementation

Source files

• 2DOF Projectile motion: launched_ball_simple.mdl
• 2DOF Projectile motion with logarithmic wind resistance: launched_ball.mdl
• 3DOF Projectile motion: launched_ball_3D.mdl
• 3DOF Projectile motion with drag effects: launched_ball_drag_3D.mdl
• 3DOF Projectile motion with drag, Coriolis effects: launched_ball_drag_3D_coriolis.mdl
• 4DOF Projectile motion with drag, Coriolis, lift effects: launched_ball_drag_4DOF_coriolis_simple_lift.mdl
• Optimization scripts: launched_ball_opt.m, launched_ball_opt_angle.m, sim_launched_ball.m, bisection.m
• FlightGear simulation models: launched_ball_w_flightgear_sim_2D.mdl, launched_ball_w_flightgear_sim_3D.mdl
• FlightGear run batch file: runfg.bat
• MATLAB .mat files containing default model parameters: launched_ball_def_params.mat, launched_ball_def_params3d.mat

Models

2DOF Projectile motion

MATLAB launched_ball_simple.mdl Causal Block Diagram

• F_i: magnitude of the impulse force.
• R: magnitude of the wind resistance force. Can represent a head- or tailwind.
• theta: angle at which projectile is launched
• g: gravitational force. 9.8 m/s² for all models
• m: projectile mass

Blocks:

• Impulse (builtin MATLAB block)
  • Provides a discrete impulse signal after a specified delay period. See here for more details: http://www.mathworks.com/access/helpdesk/help/toolbox/dspblks/ref/discreteimpulse.html
• sim_x & sim_y blocks: store the signal outputs and allow for further analysis after a MATLAB simulation is run.

2DOF Projectile motion with logarithmic wind resistance

MATLAB launched_ball.mdl Causal Block Diagram

• None

3DOF Projectile motion:

MATLAB launched_ball_3D.mdl Causal Block Diagram

Parameters:

• Phi: See Design/Models: 3DOF systems

3DOF Projectile motion with drag effects

MATLAB launched_ball_drag_3D.mdl Causal Block Diagram

With the addition of drag effects, a non-linear model had to be employed, as specified in the design section. Note that the velocity output signal is routed through the drag blocks and squared, then added to the output. This technique was also used in subsequent models to allow modeling of non-linear components.

Parameters:

• C_d: Drag coefficient
• A: Projected area of projectile body causing drag

Blocks:

• COESA Atmosphere Model: (built-in MATLAB block)
  • Description: Given a certain altitude, the block computes values for absolute temperature, pressure, density, and speed of sound. More details here: http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/coesaatmospheremodel.html
  • Parameters:
    • Units (metric/English)
• Wind shear model (built-in MATLAB block)
  • Description: Implements a wind shear given a certain altitude http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/windshearmodel.html
  • Parameters:
    • Wind Speed at 6m altitude
    • Wind direction at 6m altitude (degrees clockwise from north)*
      • *In case of model, north is always considered the "forward" direction of the projectile. Ex: 90 degrees would always result in a crosswind, despite the direction the projectile is launched
• Clean signal noise
  • Description: Considers the absolute value of the signal less than a certain tolerance to be zero. Prevents near-zero values from appearing in the final output
    • Parameters: Threshold
• Normalize (built-in MATLAB block)
  • Description: Determines the unit vector in the direction of the supplied vector
  • Parameters:
    • None

3DOF Projectile motion with drag, Coriolis effects:

MATLAB launched_ball_drag_3D_coriolis.mdl Causal Block Diagram


Coriolis effect subsystem

Additional Parameters:

• lat: latitude at which the projectile is located
• omega: earth's angular speed

4DOF Projectile motion with drag, Coriolis, lift effects

MATLAB launched_ball_drag_4DOF_coriolis_simple_lift.mdl Causal Block Diagram


Lift effects subsystem

Additional Parameters:

• A_w: the wing area generating lift

Additional Blocks

• Saturation (builtin MATLAB block)
  • limits the incoming signal to specific minimum and maximum values
  • parameters:
    • minimum value
    • maximum value

FlightGear Simulation:

2 Degree of Freedom simulation model

MATLAB launched_ball_w_flightgear_sim_2D.mdl Causal Block Diagram

Parameters:

• IA: The initial altitude at which the FlightGear model will be placed

Blocks:

• Ballistics simulation subsystem
  • Description: Contains the 2DOF projectile motion with logarithmic wind resistance model
• Flat Earth to LLA: (built-in MATLAB block)
  • Description: Translates the projectile trajectory from flat earth/local body coordinates to geodetic coordinates (see http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/flatearthtolla.html for more details)
  • Parameters:
    • Initial longitude & latitude array
    • Direction of flat earth x-axis (degrees clockwise from north)
• Simulation Pace: (built-in MATLAB block)
  • Description: sets the pace at which animation will be run w.r.t. the MATLAB simulation so that animations are smooth (see http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/simulationpace.html for more details)
  • Parameters:
    • Simulation pace (simulation sec per clock sec)
    • Sample time
• Generate Run Script: (built-in MATLAB block)
  • Description: Generates the batch file that will setup FlightGear for control by the MATLAB model (see http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/generaterunscript.html for more details)
  • Parameters: see http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/generaterunscript.html

4 Degree of Freedom simulation model

MATLAB launched_ball_w_flightgear_sim_3D.mdl Causal Block Diagram

Parameters as previously specified

Blocks:

• Ballistics simulation subsystem
  • Description: Contains the 4DOF projectile motion with drag, Coriolis, lift effects model

Experiments & Results

2 DOF systems

Parameters:

• Fi= 1500 N
• R= -5 N
• theta = pi/5 rad
• m = 10 kg
• g = 9.8 m/s²

3 & 4 DOF Systems:

Experiments compared to headwind-type wind shear:

Parameters:

• Fi = 1500 N
• R = 5 N
• theta = pi/5 rad
• phi = 0 rad
• m = 10 kg
• g = 9.8 m/s²
• A = 0.00785 m²
• A_w = 0.002 m²
• Cd = 0.295
• lat = 37.98°
• omega = 7.292E-05 rad/s
• wind speed = 30m/s
• wind direction = 0°

Experiments compared to crosswind-type wind shear:

Parameters:

• All parameters as specified in previous case except for crosswind trajectory (red), where wind direction is directed -90 degrees from north

Optimization experiments:

The above plot illustrates the solution to finding the optimal projectile launch direction. In black is the solution trajectory. The target location is marked by the red crosshair. The bisection method's progress is demonstrated by the blue trajectories.

Script command:
launched_ball_opt_angle(m, g, R, target_x, target_y, Fi_init, t, TOL, theta_min, theta_max, interm_plots)
For full argument descriptions, see launched_ball_opt_angle.m

Script command for above plot:
launched_ball_opt_angle(10, 9.8, 0, 1500, 100, 1500, 0:.1:30, 1E-7,0,pi/2, true)

The above plot ilustrates the solution to the minimization problem, where the black trajectories depict the intermediate solutions found by the minimization algorithm, and the blue line represents the optimal force & direction trajectory. The target is marked by the red crosshair.

Script command:
launched_ball_opt(m, g, R, target_coords, t, TOL, theta_min, theta_max, force_range, interm_plots)
For full argument descriptions, see launched_ball_opt.m

Script command for above plot:
launched_ball_opt(10, 9.8, 0, [1200 0], 0:.1:20, 1E-7, 0, pi/2, [1100, 1200], true)

FlightGear simulations:

• 2 DOF FlightGear simulation quicktime movie: flightgear_2DOF_simulation.m4v
  • The simulation depicts the path of a 2 degree of freedom system (launched_ball model) with the parameters as specified in the 2DOF experiment
• 4 DOF FlightGear simluation quicktime movie: 4DOF_flightgear_simulation.m4v
  • The simulation depicts the path of a 4 degree of freedom system (launched ball with drag, lift & coriolis effects) with all parameters specified as in 4DOF experiment, with the following exceptions:
    • wind direction: 90 degrees
    • wind shear speed: 50 m/s

Conclusions

The dynamics models employed where able to provide a qualitative description of the behaviour of projectiles in different environments. For true quantitative descriptions, more complex equations of motion would have to be employed. Improvements could be made in modelling the fluid/body interactions, lift and induced drag forces, and atmospheric conditions.

Nevertheless, the proposed models and experiment results provide motivation for further investigations.

References

Clancy, L. J. (1975). Aerodynamics. Toronto, ON: Copp Clark Publishing Company.

Coriolis Acceleration, http://scienceworld.wolfram.com/physics/CoriolisAcceleration.html

Lift Coefficient, http://scienceworld.wolfram.com/physics/LiftCoefficient.html

Laplace Transform, http://mathworld.wolfram.com/LaplaceTransform.html

fminbnd, http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fminbnd.html

Flat Earth to LLA, http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/flatearthtolla.html

COESA Atmosphere Model, http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/coesaatmospheremodel.html

Wind Shear Model, http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/windshearmodel.html

Working with the Flight Simulator Interface, http://www.mathworks.com/access/helpdesk/help/toolbox/aeroblks/f3-19546.html

Presentation

Modeling simple aerospace dynamics using CBDs.pptx