channelflow.org

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Numerical simulation of a baseball in flight ======

This labs shows how to do computer simulations of some basic physics, specifically, the flight of a baseball hit from home plate. We'll model the forces on the flying ball, including gravity and air resistance to produce the Newtonian equations of motion for the baseball in flight. The solution to these equations is the path of the baseball through the air.

Unfortunately, the equations of motion are too complicated to find an explicit formula for the solution using pencil and paper mathematics! So we'll use computer simulation techniques to approximate the solution numerically. We'll graph the path of the ball to see if it clears the outfield fence. Then we'll adjust the speed and angle of the ball off the bat to try to determine the minimum speed and optimal angle to hit a home run.

After that, we'll adjust the air density parameter to explore the difference between hitting a home run at Fenway and hitting one at Mile High Stadium in Denver.

==== Equations of motion, without air resistance====

The equations of motion for a baseball in flight with no air resistance are

$ dx/dt = v_x $

$ dy/dt = v_y $ 

$ dv_x/dt = 0$ 

$ dv_y/dt = -g$

Here $x$ and $y$ are the horizontal and vertical distances from home plate, $v_x$ and $v_y$ are the components of velocity in the $x$ and $y$ directions, and $g$ is the acceleration due to gravity.
$dx/dt$ means the rate of change of $x$, and $dv_x/dt$ means the rate of change of the $x$ velocity. 

Note that these equations relate //the rates of change// of a set of variables to the //values// of those variables. This kind of system is called an Ordinary Differential Equation, or ODE.

We'll now implement some Matlab code to solve these equations numerically. We write a function ''f_nodrag'' in file ''fnodrag.m'' which defines the above equations of motion as a function mapping $x,y,v_x,v_y$ into $dx/dt, dy/dt, dv_x/dt, dv_y/dt$. 

<code matlab>
function dudt = f_nodrag(t, u) 
% equations of motion for baseball in flight without air resistance

    g = 9.8;    % acceleration of gravity m/s^2

    % extract components of vector u into variables x, y, vx, vy
    x  = u(1);
    y  = u(2);
    vx = u(3);
    vy = u(4);

    % compute rates of change of those components
    dxdt = vx;
    dydt = vy;
    dvxdt = 0.0;
    dvydt = -g;

    % pack rates of change back into vector dudt, return value of function
    dudt = [dxdt; dydt; dvxdt; dvydt];
end
</code>

Matlab expects equations like this to be written in vector form. So we let $u$ be a vector of the positions and velocities $ u = [x, y, v_x, v_y]$ and $du/dt$ be the vector of the rates of change of those quantities $ du/dt = [dx/dt, dy/dt, dv_x/dt, dv_y/dt]$.


==== Equations of motion with air resistance modeling ====


The equations of motion for a baseball in flight with turbulent air resistance are

$ dx/dt = v_x $

$ dy/dt = v_y $ 

$ dv_x/dt = -\frac{\mu}{m} v_x \sqrt{v_x^2 + v_y^2}$ 

$ dv_y/dt = -\frac{\mu}{m} v_y \sqrt{v_x^2 + v_y^2} - g$

Here $\mu$ is a coefficient determined by the air density and the size of the baseball, $m$ is the mass of the baseball, and $g$ is the acceleration due to gravity. We encode those equations into the function ''f_withdrag'' as follows

<code matlab>
function dudt = f_withdrag(t,u) 
% equations of motion for baseball in flight with turbulent air resistance

    rho_air  = 1.196;     % kg/m^3, density of dry air, 21 C, sea level   (Fenway)
    C = 0.3;              % drag coefficient for baseball (rough sphere)
    g = 9.81;             % acceleration due to gravity in m/s^2
    r = 0.0375;           % radius of baseball in m (3.75 cm)
    A = pi*r^2;           % cross-sectional area of baseball in m^2
    m = 0.145;            % mass of baseball in kg (145 gm)
    mu = rho_air*C_D*A/2; % coefficient of nonlinear |v|^2 term, in mks units

    % extract components of vector u into variables x, y, vx, vy
    x  = u(1);
    y  = u(2);
    vx = u(3);
    vy = u(4);

    % compute rates of change of those components
    dxdt = vx;
    dydt = vy;
    dvxdt = -mu/m * vx * sqrt(vx^2 + vy^2);
    dvydt = -mu/m * vy * sqrt(vx^2 + vy^2) - g;

    % pack rates of change back into vector dudt, return value of function 
    dudt = [dxdt; dydt; dvxdt; dvydt]
end
</code>


==== Solving the equations of motion numerically ====


The file ``baseballsolve.m`` solves the above equations with some numerical solution methods, given
the initial position and speed of the ball, and makes a plot that compares the flight of the ball with 
and without air resistance.

<code matlab>
% Compute the trajectory of a baseball in flight, with and without 
% air resistance and plot results

x = 0.0;      % horizontal position of home plate, meters
y = 1.0;      % height of ball over strike zone
v = 40.0;     % initial speed of ball, meters per second
theta = 76;   % angle of ball off bat, degrees

% initial position and speed of ball
u0 = [x, y, v*cosd(theta), v*sind(theta)]; 

% solve for the first ten seconds of flight
tspan = [0 10];  

% solve the baseball equations numerically. ode45 is Matlab's numerical algorithm for ODEs
[t, u_withdrag] = ode45(@f_withdrag, tspan, u0);
[t, u_nodrag]   = ode45(@f_nodrag,   tspan, u0);

% get the x,y coordinates from the solutions and plot
figure(1); clf()
x = u_nodrag(:,1);
y = u_nodrag(:,2);
plot(x, y, 'r-')
hold on

x = u_withdrag(:,1);
y = u_withdrag(:,2);
plot(x, y, 'b-')

% plot outfield fence at x = 120 m (390 ft) and y = 0 to 5.2 m (17 ft)
plot([120, 120], [0, 5.2], 'k-', 'linewidth',2) 

legend('no drag', 'with drag', 'outfield fence')

grid on
xlim([-10, 130])
ylim([0, 110])
xlabel("x, meters")
ylabel("y, meters")
title("baseball trajectory")
</code>
==== Running the Matlab code to see the results ====

You can run the Matlab code by typing ``baseballsolve`` at the Matlab prompt.

<code matlab>
>> baseballsolve
</code>

A plot like this should appear

{{ :gibson:teaching:summer-2017:techcamp:baseballsolve.png?600 |}}

The thick black line shows the outfiled fence at 390 ft or 120 m, and 17 ft or 5.2 m high. 

**Questions**

  * How important is the effect of air resistance? 

  * Air density at Fenway park (sea level) is 1.196 kg/m^3. At one mile high in Denver, it's 0.986 kg/m^3. How much does the change in density affect the path of the ball at different angles and speeds? Try to change the above codes to find out. See if you can find a speed and angle that gives a home run in Denver but not at Fenway. 

  * Play around with the parameters in the functions, like the speed and angle of the ball. Can you figure out the minimum speed and optimal angle to just clear the outfiled fence and hit a home run?
channelflow.org

User Tools

Site Tools

Page Tools
