Matlab projectile motion with drag

Dynamics MATLAB Project

Computing Project 1 Projectile

Introduction

The computing project Projectile concerns the theory of particle dynamics. The concept is to model projectile motion including the nonlinear drag resistance caused by motion in a fluid with finite density (which causes drag that is proportional to the square of the ve- locity). The drag forces can have a significant impact on the motion of the particle and that impact can affect the outcome in practical problems (e.g., in targeting and homerun hitting).

The most widely accepted model of velocity related drag forces give a force with magni- tude equal to

𝐹𝐹𝑑𝑑 = 12 𝜌𝜌𝐶𝐶𝑑𝑑𝐴𝐴𝑝𝑝 𝑣𝑣 2 ≡ 𝑐𝑐𝑣𝑣2 (1)

where ρ is the density of the surrounding fluid, Cd is the coefficient of drag (which de- pends upon the shape of the object—something you can look up online or in a handbook), Ap is the projected area of the object perpendicular to the direction of the velocity, and 𝑣𝑣 = √𝐯𝐯 ∙ 𝐯𝐯 is the speed (the norm of the velocity v). For convenience in the following discussion we lump the coefficients into the single parameter 𝑐𝑐 = 𝜌𝜌𝐶𝐶𝑑𝑑𝐴𝐴𝑝𝑝/2. In your pro- gram you will want to input each of these as physical properties and then compute c for use throughout the computation.

The forces that act on a particle are vectors. The drag force must, therefore, be a vector (with magnitude and direction). The expression above gives the magnitude of the vector. The direction of the force is in the direction of the velocity and acts to oppose motion in the direction of the velocity. The simplest way to develop a unit vector in the direction of the velocity is to divide the velocity by its magnitude. The direction m of the drag force is then

𝐦𝐦 = −𝐯𝐯/𝑣𝑣 (2)

Putting these together the drag force can be expressed as

𝐅𝐅𝑑𝑑 = 𝑐𝑐𝑣𝑣2𝐦𝐦 = −𝑐𝑐 𝑣𝑣 𝐯𝐯 (3)

The force opposes the motion, it is proportional to the constant c, and it is proportional to the “square” of the velocity. The presence of this force makes the particle motion prob- lem nonlinear because the velocity is part of the state that we are trying to determine by integrating the equations of motion.

We can solve the problem numerically by computing the velocity and position incremen- tally based upon the state at the previous time step and the equation of motion at the cur- rent time step. Using the generalized trapezoidal rule for numerical integration we have

Arizona State University The Mechanics Project CEE 212—Dynamics CP 1—Projectile

© 2013 Keith D. Hjelmstad 2

𝐯𝐯𝑛𝑛+1 = 𝐯𝐯𝑛𝑛 + ∆𝑡𝑡[𝛽𝛽𝐚𝐚𝑛𝑛 + (1 − 𝛽𝛽)𝐚𝐚𝑛𝑛+1] (4)

𝐱𝐱𝑛𝑛+1 = 𝐱𝐱𝑛𝑛 + ∆𝑡𝑡[𝛽𝛽𝐯𝐯𝑛𝑛 + (1 − 𝛽𝛽)𝐯𝐯𝑛𝑛+1] (5)

where ∆𝑡𝑡 = 𝑡𝑡𝑛𝑛+1 − 𝑡𝑡𝑛𝑛 is the time increment and β is a numerical integration parameter, which we generally take to be β=0.5 (more on that in the Course Notes entitled Numeri- cal Methods). The subscripts on the state variables x, v, and a indicate the time—e.g., xn = x(tn) and xn+1 = x(tn+1), etc. Even though the position vector is a continuous function of time we will be finding values at certain discrete values of time that are meant to approx- imate the actual function at that time.

We get acceleration from the equations of motion F=ma. In the present case that would be

−𝑚𝑚𝑚𝑚𝐧𝐧− 𝑐𝑐𝑣𝑣𝑛𝑛+1𝐯𝐯𝑛𝑛+1 = 𝑚𝑚𝐚𝐚𝑛𝑛+1 (6)

where m is the mass of the particle, g is the acceleration of gravity (i.e., 9.81 m/s2), and n is a unit vector pointing in the direction of gravity (usually we will set our coordinate sys- tem such that n = e3).

The knowledge you need to solve this problem within the context of a MATLAB code is contained in the Course Notes entitled Numerical Methods. The state is completely de- termined at 𝑡𝑡𝑛𝑛 (we got it from doing the same calculation at the last time step). So, we can look at Eqns. (4), (5), and (6) as three (vector) equations in three (vector) unknowns.1

The task is to solve those equations at each time step and then move on to the next time step. As in all dynamics problems we must specify initial condition on position and ve- locity (i.e., xo and vo must be known as part of the problem specification). To get ready to start the time marching algorithm you must also know the acceleration ao at time zero. The acceleration cannot be specified, but rather must be computed from the equation of motion, in this case Eqn. (6) so

𝐚𝐚𝑜𝑜 = −𝑚𝑚𝐧𝐧 − 𝑐𝑐 𝑚𝑚 𝑣𝑣𝑜𝑜𝐯𝐯𝑜𝑜 (7)

Once the initial state is completely established then the time stepping algorithm can be done. Note that Eqn. (6) is nonlinear (because of the velocity term, which is sort of a vec- tor version of velocity squared). Therefore, we will need to use Newton’s method at each time step to solve our system of equations (again, see the course notes).

The time-stepping algorithm for projectile motion without drag has been implemented in the MATLAB program projectile.m which we have provided for you and is available on the Blackboard website. This program shows how to set up the time-stepping part of the MATLAB code. Note that without the drag term the equation of motion—i.e., Eqn. (6)

1 If you write out all of these three vector equations in components then it is nine scalar equations in nine scalar unknowns. It will be convenient to keep the equations in vector notation because we can implement these equations in MATLAB as vector equations.

Arizona State University The Mechanics Project CEE 212—Dynamics CP 1—Projectile

© 2013 Keith D. Hjelmstad 3

without drag—are not nonlinear so we do not have a Newton “while” loop in this code. You will need to put that in as part of this project. If you think about it, what you really need to do to get from my code to your code is to deal with what is different between my equation of motion and your equation of motion (yours has the drag term, mine does not). Note that the numerical integration equations—Eqns. (5) and (6) above—are the same in both cases because they only implement the fact that acceleration is the time rate of change of velocity which is the time rate of change of position.

The MATLAB program projectile.m is the program framework that we will use for all of the Computing Projects in this course. Therefore, understanding how this one is set up will help you to get oriented for all future CPs. The descriptions of some of the features are noted in the course notes Numerical Methods. For future computing projects you will have this basic code for a launching point (the time stepping algorithm will be basically the same in each case, the equations of motion will change). Study the code carefully.

What you need to do

For this computing project you can get going in two steps: (1) get the program to com- pute correctly using an explicit approximation (see the course notes) and (2) do the code for the real thing—i.e., the implicit implementation.

As mentioned in Numerical Methods, the explicit approximation uses the simplifying trick of using the old velocity 𝐯𝐯𝑛𝑛 in the equation of motion instead of the correct new val- ue 𝐯𝐯𝑛𝑛+1. That allows us to compute the acceleration at the new state (approximately) as

𝐚𝐚𝑛𝑛+1 = −𝑚𝑚𝐧𝐧 − 𝑐𝑐𝑣𝑣𝑛𝑛𝐯𝐯𝑛𝑛/𝑚𝑚

This equation can be implemented in projectile.m by simply adding the drag term. Be- cause it is the old velocity (which is known from the last time step) you don’t need New- ton’s method to solve the equation of motion. You should be able to get this working by changing a few lines of projectile.m!2

Once you get the explicit version of the code working move on to the implicit version. This part is the main challenge and purpose of this project. In essence, you need to im- plement Newton’s method (in a “while” loop) to satisfy the equations of motion—Eqn. (6)—at each time step. The course notes Numerical Methods shows you how to get the equations for the Newton iteration.

2 This approach—solving a simpler problem before you get into the complexities of the real prob- lem at hand—is a good programming practice. It allows you to sort things out one step at a time and at each step you can get a working code that might help to verify the next version. When you get the explicit version working you will want to save that for posterity (i.e., start the next code in a new file so you can go back to that one as needed).

Arizona State University The Mechanics Project CEE 212—Dynamics CP 1—Projectile

© 2013 Keith D. Hjelmstad 4

Once you get your code working, use it to study the phenomenon of terminal velocity. You can, in fact, use terminal velocity as one of the ways to verify the code because it is possible to derive an analytical expression from the equations of motion for terminal ve- locity (assuming that the direction does not change, which would be true for a particle traveling in exactly the same direction as gravity, velocity parallel to n). Examine what happens when the projectile velocity is not in a constant vertical direction—i.e., the gen- eral projectile-motion problem, wherein there is both a horizontal and vertical component to the motion. How does horizontal motion affect the vertical terminal speed? Will a ball that hits water in a swimming pool at an angle hit the bottom at the same time as a ball that hits going straight down if the vertical component of the velocity is the same?

You are, of course, free (and encouraged) to explore anything you want within the con- text of particle motion with drag forces. Note that your report will be evaluated on the basis of the interest and insights of your study. This part of the assignment invites you to discover and to explore. You could change your code to account for ambient wind (i.e., the fluid moving on its own independent of the particle) and then you could explore why it is harder to hit a homerun in baseball hitting into the wind compared to with the wind. You can, and should study how the numerical analysis parameters affect your solution. What happens if you change the time step? What happens if you change the parameter β?

Write a report documenting your work and the results (in accord with the specification given in the document Guidelines for Doing Computing Projects). Post it to BlackBoard prior to the deadline. Consult the document Evaluation of Computing Projects to see how your project will be evaluated to make sure that you can get full marks.