Motion of Falling Objects With Resistance: Doing Physics With Matlab Mechanics

DOING PHYSICS WITH MATLAB
MECHANICS
MOTION OF FALLING OBJECTS WITH RESISTANCE
Ian Cooper
School of Physics, University of Sydney
[email protected]
DOWNLOAD DIRECTORY FOR MATLAB SCRIPTS
mec_fr_mg_bv.m
Computation of the displacement, velocity and acceleration for the motion of an object
acted upon by a resistive force FR    v and its weight FG  m g . The equation of
motion is solved by analytical means (integration of the equation of motion) and by a
finite difference numerical method.
mec_fr_mg_bv2.m
acted upon by a resistive force FR   v 2 and its weight FG  m g . The equation of
motion is solved by analytical means (integration of the equation of motion) and by a
finite difference numerical method.
mec_stokes.m
acted upon by a resistive force FR  6 F R v (Stokes Law), the buoyance force Fbuoy
and its weight FG  m g . The equation of motion is solved by a finite difference
numerical method.
mec_drag.m
acted upon by a resistive force FR  12 CD 6 A  F v 2 and its weight FG  m g . The
equation of motion is solved by a finite difference numerical method.
mec_tt_ball.m
Experimental data for a table tennis ball falling from rest. Used to plot the actual
measurements of displacement and time from a video recording with the predicted
values of displacement using a finite difference approach.
http://www.physics.usyd.edu.au/teach_res/mp/mphome.htm 1
INTRODUCTION
We will consider the vertical motion of objects through fluids near the Earth’s
surface where the acceleration due to gravity is assumed to be constant
g = 9.80 m.s-2.
The motion of falling objects is usually described with constant acceleration. This
is only approximately true. For example, in introductory physics textbooks, two
objects of different mass when dropped simultaneously from rest will hit the
ground at the same time. This is an idealized situation and ignores the effects of
the air resisting the motion of the falling objects. Air resistance, a friction which
increases with increasing speed, acts against gravity, so the speed of falling
objects tends toward a limit called terminal velocity (terminal speed).
Motion through any real fluid (liquid or gas) gives rise to forces resisting the
motion. To a reasonable approximation, fluid resistance tends to depend on
either the first power of the speed (a linear resistance) or the second power (a
quadratic resistance). Our two models for the resistive force FR are
Model (1) FR    v low speeds
Model (2) FR   v 2 high speeds
where  and  are constants of proportionality.
Model (1) for linear resistance is often applicable when the object is moving with
low speeds. In the motion through a fluid, the resistive force FR    v is often
called the viscous drag and it arises from the cohesive forces between the layers
of the fluid. The S.I. units for the constant  are N.m-1.s-1 or kg.s-1.
Model (2) for quadratic resistance is more applicable for higher speeds. In the
motion through fluids, the resistive force FR   v 2 is usually called the drag
and is related to the momentum transfer between the moving object and the
fluid it travels through. The S.I. units for the constant are N.m-2.s-2 or kg.m-1.
Many problems in the mathematical analysis of particles moving under the

influence of resistive forces, you start with the equation of motion. To find
velocities and displacements as functions of time you must integrate the
equation of motion.
The equation of motion for an object can be derived from Newton’s Second Law
1
a
m
F
i
i Newton’s Second Law
For the vertical motion of an object through a fluid, the forces acting on the
object are the gravitational force FG (weight) and the resistive force FR. In our
frame of reference, we will take down as the positive direction.
Therefore, the vertical acceleration a of an object through a fluid is
a  g    / m v Model (1)
a  g   / m  v 2  v / v  Model (2)
In Model 2
If the object is moving down then v > 0
a  g   / m  v 2  v / v   g   / m  v 2  1  g   / m  v 2
If the object is moving up then v < 0
a  g   / m  v 2  v / v   g   / m  v 2  1  g   / m  v 2
MODEL 1 FR    v ma  m g   v
The mscript used for the plots is mec_fr_mg_bv.m
Analytical Approach FR    v ma  m g   v
frame of reference, we will take down as the positive direction.
The equation of motion of the object is determined from Newton’s Second Law.
dv 
ma  m  FR  mg   v ag v
dt m
where a is the acceleration of the object at any instance.
The initial conditions are t  0 v  v0 x  0 a     / m  v0
When a = 0, the velocity is constant v = vT where vT is the terminal velocity
0  m g   vT
mg
vT  terminal velocity

We start with the equation of motion then integrate this equation where the limits
of the integration are determined by the initial conditions (t = 0 and v = v0) and
final conditions (t and v)
dv      mg
a  g   v   v 
dt m  m   
mg   du
uv du  dv    dt 
 m u
 mg 
  du   v
mg
 v  
   t  log e  u   m g
t u
    dt    log e  
m 0 u 0 u m v0 
  v0  m g 
  
 
 mg 
 v     / m t
  e
 v0  m g 
  
 
mg  m g    / m t mg
v   v0  vT 
 
e
  
Therefore, we can express the velocity v as
v  vT   v0  vT  e 
  / m t
velocity

v0  0  v  vT 1  e    / m  t 
v0  vT  v  vT
v0  vT  v increases to vT
v0  vT  v decreases to vT
In every case, the velocity v tends towards the limiting value vT.
Plots of the velocity v as functions of time t
m = 2.00 kg
 = 5.00 kg.s-1
g = 9.80 m.s-2
vT = 3.92 m.s-1
Initial values for velocity v0 [m.s-1]
blue: 10 red: vT magenta: 2 cyan: 0
The acceleration a as a function of time t is
v  vT   v0  vT  e 
  / m t
a
dv d

dt dt

vT   v0  vT  e 
  / m t

 
a   vT  v0    e    / m t
m
  vT    / m t
v0  0  a   e
 m 
 a g e
  / m t
v0  vT  a  0
v0  vT  a  0 and decreases to 0
v0  vT  a  0 and a increases to 0
t a0
We can now calculate the displacement x as a function of velocity t
v  vT   v0  vT  e    / m  t
dx
v dx  v dt
dt

dx   vT   v0  vT  e    / m  t dt 
x t
0 0
t
 m 
x   vT t     v0  vT  e    / m  t 
   0
m m
x  vT t     v0  vT  e    / m  t     v0  vT 
   
m

x  vT t     v0  vT  1  e    / m  t
 

t  x 
 m 
v0  0  x  vT  t    e    / m  t  1   
   
So far we have only considered the case where the initial velocity was either zero
or a positive quantity ( v0  0 ), i.e., the object was released from rest or
projected downward. We will now consider the case where the object was
project vertically upward (v0 < 0). Note: in our frame of reference, the origin is
taken as x = 0, the position of the object at time t = 0; down is the positive
direction and up is the negative direction.
When the object is launched upward at time t = 0, the initial velocity has a
negative value. Let u be the magnitude of the initial velocity v0
v0  0 v0  u u  0
Therefore, the equation for the velocity v as a function of time t can be

expressed as
v  vT   v0  vT  e 
  / m t
v  vT   u  vT  e 
  / m t
We can now find the time tup it takes for the object to rise to its maximum height
xup above the origin (remember: up is negative). At the highest point v = 0,
therefore,
   / m  tup
0  vT   u  vT  e
m  u
tup    log e  1  
   vT 
The maximum height xup reached by the object in time t = tup is
m

x  vT t     v0  vT  1  e    / m  t
 

m
xup  vT tup     u  vT  1  e

 

   / m  tup

For the parameters
m = 2.00 kg  = 5.00 kg.s-1 g = 9.8 m.s-2 u = 10 m.s-1 vT = 3.92 m.s-1
The time to reach maximum height is tup = 0.507 s
The max height hup reached is hup = 2.013 m xup = - 2.013 m
We can find the displacement x as a function of velocity v
dv v dv
av  g    / m v dx 
dx g    / m v
 1  v dv v dv
  / m  dx  mg /   v vT  mg /    / m  dx 
    vT  v
We can integrate this equation by a substitution method or an algebraic

manipulation method.
Substitution Method
u  vT  v du   dv v  vT  u dv   du v0  vT  u0 u0  vT  v0
  vT  u 
  / m  dx  du
u
u  v  u u  v 
0   / m  dx  u0 u du  u 0 1  uT
x
T
 du

u
  / m  x  u  vT log e  u  u   u  u0   vT log e  
u
0
 u0 
 vT  v0 
  / m  x   v0  v   vT log e  
 vT  v 
 m   v  v0  
x      v0  v   vT log e  T 
 b   vT  v  
Algebraic manipulation
v dv
  / m  dx 
vT  v
v dv  1 1 
 vT   
vT  v  vT vT  v 
  1 1  
   / m  dx  v  
x vu
   dv 
 vT vT  v  
T

0 v0
v
 v 
  / m  x  vT   log e  vT  v  
 vT  v0
  vT  v0  
  / m  x    v0  v   vT log e  
  vT  v  
 m   v  v0  
x      v0  v   vT log e  T 
 b   vT  v  
QED
For the object projected up with an initial velocity v0 = - u where u > 0, the
maximum height reached xup occurs when v = 0
 m   v  v0  
x      v0  v   vT log e  T 
 b   vT  v  
 m   v  u 
xup      u   vT log e  T 
 b   vT  
 m   u 
xup     vT log e  1    u 
 b   vT  
Note: up is negative and down is positive in our frame of reference.
Numerical Approach FR    v ma  m g   v
We can also find the velocity and displacement of the object by solving Newton’s
Second Law of motion using a finite difference method.
We start with
dv mg
a  g    / m v vT 
dt 
In the finite difference method we calculate the velocity v and displacement x at

N discrete times tk at fixed time intervals t
t1, t2, … , tk, …, tN t =t2 - t1 t1 = 0 tk = (k-1) t k = 1, 2, 3, … , N
The acceleration a is approximated by the difference formula
dv  tk 1  v  tk 2   v  tk 

dt 2 t
Therefore, the velocity v  tk 2  at time tk+2 is
v  tk  2   v  t k 
 g    / m  v  tk 1 
2 t
v  tk 2   v  tk    2 t  g   2 t  / m  v  tk 1 
 v  tk 1  
v  tk 2   v  tk    2 t g   1  
 vT 
Hence, to calculate the velocity v  tk 2  we need to know the velocity at the two
previous time steps tk 1 and tk . We know t1 = 0 and v(t1) = v(0) = v0.
We estimate the velocity at the second time step t2
v(t2 )  v(t1 )  a(t1 ) t  v(t1 )    / m  v  t1  t
where we have assumed a constant acceleration in the first time step. We can
improve our estimate of v(t2 ) by using an average value of the acceleration in
the first time step
 a ( t )  a ( t2 ) 
v(t1 )   1  t  v (t2 )
 2
v(t1 )  12   / m   v  t1   v  t2   t  v (t2 )
We can now calculate the velocity v(t) at all times from t = t1 to t = tN.
The acceleration at each time step is
 
a  tk   g    v  t k 
m
The displacement at each time step is
dx x
v 
dt t
x  tk  2   x  t k 
v  tk 1  
2 t
x  tk 2   x  tk    2 t  v  tk 1 
EXAMPLE FR    v ma  m g   v
The mscript mec_fr_mg_bv.m can be used for simulations for the motion of an
object acted upon a resistive force of the form FR    v (Model 1).
Input parameters
m = 2 kg  = 10 kg.s-1 v0 = 10 m.s-1 t = 1x10-4 s tmax = 2.00 s
Outputs N numerical approach A analytical approach
terminal velocity vT = 1.96 m.s-1
acceleration t  a 0
velocity t   v  vT
displacement t  x
For the input parameters used in this simulation, there is excellent agreement
between the values calculated using the numerical and analytical approaches.
However, you always need to be careful in using numerical approaches to solve

problems. In this instance, you need to check the convergence of results by
progressively making the time step t smaller.
When the time step is

t = 5x10-2 s the numerical and
analytical results do not agree. The
time step is too large for accurate
results using the numerical
approach.
MODEL 2 FR   v 2 m a  m g   v2
The mscript used for the plots is mec_fr_mg_bv2.m
Analytical Approach FR   v 2 m a  m g   v2
frame of reference, down is the positive direction.
The equation of motion of the object is determined from Newton’s Second Law.

 FG  FR  mg   v 2  v / v  v 2 v / v 
dv
ma  m ag
dt m
where a is the acceleration of the object at any instance.
v 
The initial conditions are t  0 v  v0 x  0 a  g   / m  v02  0 
 v0 
When a = 0, the velocity is constant v = vT where vT is the terminal velocity
mg
0  m g   vT 2 vT 2 

mg
vT 

We start with the equation of motion then integrate this equation where the
limits of the integration are determined by the initial conditions (t = 0 and v = v0)
and final conditions (t and v).
Since the acceleration depends upon v2 its a more difficult problem then for the
linear resistive force example. We have to do separate analytical calculations for
the motion when the object is falling or rising.
Velocity of the object is always positive (falling object) v0  0 v  0
Equation of motion

ag v2
m
dv  
a  g    v2
dt m
dv dv mg
dt   vT 2 
    mg 2 
g    v 2     v 
m  m    
  dv  1  1 1 
   dt  2     dv
m v  vT 2
 2vT   v  vT v  vT 
 mg     1 1 
2    dt     dv
   m   v  vT v  vT 
4 g  1 1 
 dt     dv
m  v  vT v  vT 
4 g t v 1 1 

m 0 dt   
v0 v  v
 T

v  vT
 dv

4 g
t    log e  v  vT   log e  v  vT   v
v

m 0
4 g
t  log e  v  vT   log e  v  vT   v
v
m 0
4 g   v  vT   v  vT 
t  log e    log e  
m   v0  vT   v0  vT 
4 g  v  vT   v0  vT  
 t  log e   
m  v0  vT   v  vT  
4 g
 v  vT   v0  vT   t 4 g 4 g 2 4 g2 2 g
  e   
m
 v0  vT   v  vT  m mg vT 2 vT
v v   vT t v v   vT t
2g 2g
v  vT   v  vT   0 T e   v  vT  K K  0 T e
 v0  vT   v0  vT 
1 K 
v 1  K   vT 1  K  v  vT  
1 K 
  v0  vT   vT t 
2g
1  e 
v v
v  vT   0 T  
 v v   vT t 
2g
 1   0 T e 
  v0  vT  
  t 
2g
  v  v    v0  vT  e T 
v
v  vT  0 T
 t 
2g
  v  v    v  v  e vT 
 0 T 0 T 
valid only if v0  0 v  0
We can now calculate the displacement x as a function of velocity v
dv v dv
a   g   / m  v 2
dt dx
  / m   m g /   v 2 
v dv
vT 2  m g / 
dx
 m  v dv  VT 2   2v  dv
dx    2    2
    vT  v   2 g   vT  v 
2 2
x   v 2  v  2 v 
0
dx   T  
 2 g  v0  vT  v 
2 2
dv
 v 2 
x   T  log e  vT 2  v 2 
v
 2g  v0
v 2  v 2 v 2 
x   T  log e  T 2 02 
2g   vT  v 
valid only if v0  0 v  0
We can now investigate what happens as time t  
 v  v   0 
2g
 t
v  t     vT  0 T  e vT
0
  v0  vT   0 
v  t     vT
In falling, the object will finally reach a constant velocity vT (a = 0 ) which is

known as the terminal velocity.
1
t   v  vT vT  v  0 
vT  v
V 2  v 2 v 2 
x   T  log e  T 2 02   
 2g   vT  v 
In falling, as time t increases the objects displacement x just gets larger and
larger.
Example Small rock dropped from rest:
m = 0.010 kg  = 1.0010-4 kg.m-1 v0 = 0 m.s-1  vT = 31.3 m.s-1
Example Small rock thrown vertically downward (v < vT)
m = 0.010 kg  = 1.0010-4 kg.m-1 v0 = +5.00 m.s-1  vT = 31.3 m.s-1
Example Small rock thrown vertically downward (v > vT)
m = 0.010 kg  = 1.0010-4 kg.m-1 v0 = +40.0 m.s-1  vT = 31.3 m.s-1
For problems in which the object is projected vertically upward, you have to
divide the problem into two parts. (1) Calculate the time to reach its maximum
height and calculate the maximum height reached for the upward motion. (2)
Reset the initial conditions to the position at maximum height where the initial
velocity becomes v0 = 0 and do the calculations for the downward movement of
the object.
Velocity of the object negative and moving up v0  0 and v0
Equation of motion Note: up is the positive direction

ag v2 valid only if v0  0 and v0
m
dv  
a  g    v2
dt m
dv dv mg
dt   vT 2 
    mg 2 
g    v 2     v 
m  m    
  dv
  dt  2
m v  vT 2
  t v dv
  0 dt  v0 2
m v  vT 2
dx 1 x
Standard Integral a 2
x 2
   atan    C
a a
v
   1   v   m   m g   vT 
 t    atan      
 m   vT   vT   v0   vT    vT g   g 
v
v   v   vT    v   v0 
t T  atan  v      atan    atan  
 g   T   v0  g    vT   vT 
The time tup to reach maximum height occurs when v = 0
 v   0   v0    vT   v0 
tup   T  atan    atan        atan   v0  0
 g   vT   vT   g  vT 
The velocity v as a function of time t is
 v   v0   g 
atan    atan   t
 vT   vT   vT 
v0  0 and v0
 v   g  
v  vT tan atan  0   t atan  tan  -1
  vT   vT  
The displacement x as a function of velocity v is
dv v dv
a   g   / m  v 2
dt dx
  / m   m g /   v 2 
v dv
vT 2  m g / 
dx
 m  v dv  vT 2   2v  dv
dx    2    2
    vT  v   2 g   vT  v 
2 2
x  v 2  v  2v 
0
dx   T  
 2 g  v0  vT  v 
2 2
dv v0  0 and v0
v 2 
x   T  log e  vT 2  v 2 
v
2g  v0
v 2   v 2  v2 
x   T  log e  T2 2 
2g   vT  v0 
The maximum height xup reached by the object occurs when v = 0
v 2   v 2 
xup   T  log e  2 T 2 
2g   vT  v0 
Example Small rock thrown vertically upward (v0 < 0 v0 = - u u > 0)
m = 0.010 kg  = 1.0010-4 kg.m-1 v0 = - 12.0 m.s-1  vT = 31.3 m.s-1
The terminal velocity vT is
vT 2  m g / 
vT  m g /   10   9.8 / 10  m.s
2 4 -1
vT  31.31 m.s-1
When v = 0 the object reaches its maximum height xup (up is negative)
v 2   v 2 
xup   T  log e  2 T 2 
2g   vT  v0 
xup  6.855 m
The time tup to reach the maximum height
v  v 
tup    T  atan  0 
 g  vT 
tup  1.169 s
The calculations agree with the values for tup and xup determined from the graphs.
From the graphs:
x=0 time t = 2.275 s velocity v = 11.21 m.s-1
Time to fall from max height to origin tdown = (2.375 – 1.169) s = 1.206 s
takes slight longer to fall then rise to and from origin to max height
Launch speed = 12.00 m.s-1 slightly greater than return speed = 11.21 m.s-1
In the absence of any resistive forces a  g v  v0  a t v 2  v02  2 a s

At maximum height v  0 tup  1.2245 s xup  7.3469 m
Numerical Approach FR   v 2 m a  m g   v2
We can also find the velocity and displacement of the object by solving Newton’s
Second Law of motion using a finite difference method.
We start with
dv v
(4) a  g   / m  v 2  
dt v
In the finite difference method we calculate the velocity v and displacement x at

N discrete times tk at fixed time intervals t
t1, t2, … , tk, …, tN t =t2 - t1 t1 = 0 tk = (k-1) t k = 1, 2, 3, … , N
The acceleration a is approximated by the difference formula
dv  tk 1  v  tk 2   v  tk 

dt 2 t
Therefore, the velocity v  tk 2  at time tk+2 is
v  tk  2   v  t k 
 g   / m  v  tk 1  / v
3
2 t
v  tk 2   v  tk    2 t  g   / m  v  tk 1  / v
3
dv
a  tk    g   / m  v  tk  / v  tk 
3
dt
Hence, to calculate the velocity v  tk 2  we need to know the velocity at the two
previous time steps tk 1 and tk . We know t1 = 0 and v(t1) = v(0) = v0.
We estimate the velocity at the second time step t2
(13) v(t2 )  v(t1 )  a(t1 ) t
where we have assumed a constant acceleration in the first time step. We can
improve our estimate of v(t2 ) by using an average value of the acceleration in
the first time step
 a ( t )  a ( t2 ) 
v(t1 )   1  t  v(t2 )
 2
We can now calculate the velocity v(t) at all times from t = t1 to t = tN.
The displacement at each time step is
dx x
v 
dt t
x  tk  2   x  t k 
v  tk 1  
2 t
(15) x  tk 2   x  tk    2 t  v  tk 1 
Provided the time step is small enough, there is excellent agreement between
the numerical values and analytical values for the acceleration, velocity and
displacement.
REYNOLDS NUMBER and MOTION THROUGH A FLUID
To characterize the motion of an object moving through a fluid, it is useful to

define a dimensionless quantity called the Reynolds number NR
 2 R 
NR   F  v
 F 
where the fluid is characterized by its viscosity  F and density  F . The effective
radius of the object is R and its velocity is v. Different values of the Reynolds
number NR determine the different regimes of flow in which different laws of the
resistive force are valid.
Reynolds number NR Resistive force

0 to 10 Stokes’ law Fdrag   6  F R  v
10 to 300 Transition region
300 to 300 000 constant drag coefficient CD
 
Fdrag ~ 21 CD A  F v 2
> 300 000 Some variation in drag coefficient CD with v
Fdrag only roughly proportional to v2
Low Reynolds Number (0 < NR < 10)
For a small object moving through a fluid with a low velocity, the
flow around the object is essentially laminar where the fluid flows
in layers and there is no turbulence. In this situation, it is found
experimentally that the viscous drag force Fdrag acting on the
object by the fluid is directly proportional to the velocity v of the
object. For a small sphere of radius R, the drag force Fdrag is also
directly proportional to the radius R.
For a small sphere of radius R, the viscous drag force Fdrag is given by Stokes’ Law
viscous drag force Fdrag   6  F R  v STOKES LAW
where the drag force Fdrag is always in the opposite direction to the velocity v.
The quantity  F is called the viscosity of the fluid [Pa.s kg.m-1.s-1]. Viscosity is a
measure of the frictional force acting on a fluid flowing over a solid surface.
For an object falling through a fluid, we make the assumption that there are
three distinct forces acting on the object
weight FG  m g
viscous drag force Fdrag    6  F R  v STOKES’ LAW
buoyancy force Fbuoy    43  R3   F g
where in our frame of reference, down is taken as the positive direction.
Newton’s Second Law applied to the falling object gives
m a  m g   6  F R  v   43  R3   F g
The acceleration a is
a  g   6  F R / m  v   43  R3   F g / m
The terminal velocity vT (a = 0) of the object is
m g   43  R 3   F g
vT 
6  F R
The motion of tiny droplets of water falling through the air, tiny particles of rock
settling under the sea and the sedimentation of red blood cells in blood plasma,
the rise of oil drops in water are examples of the motion of objects moving
through a fluid with low Reynolds numbers.
We can use the finite difference method to numerically find the velocity v and
displacement x of the object as functions of time.
dv  tk 1  v  tk 2   v  tk 
a  tk 1   
dt 2 t
v  tk  2   v  t k 
 g   6   R / m  v  tk 1    43  R 3   g / m
2 t

v  tk 2   v  tk    2 t  g   6   R / m  v  tk 1    43  R 3   g / m 
Example: A small oil drop released in water will rise to the surface
mscript mec_stokes.m
Input parameters
radius of oil drop R = 5x10-3 m

density of oil drop  = 900 kg.m-3
density of water F = 1000 kg.m-3
viscosity of water F = 0.1 Pa.s
max time tMax = 0.4 s
No. time steps N = 10000
Output parameters
mass of oil drop m = 4.71x10-4 kg
terminal velocity |vT| = 0.054 m.s-1
time to reach vT t(vT) ~ 0.3 s
Reynolds number is in the range from 0 to 10, hence Stokes’ Law  FR  v  is

valid. The main reason the oil drop rises is because of the buoyancy force, the
viscous drag force plays only a minor role.
Example: How do small water droplets fall?
mscript mec_stokes.m
Input parameters
radius of water drop R = 4x10-5 m

density of water drop  = 1000 kg.m-3
density of air F = 1.20 kg.m-3
viscosity of air F = 1.80x10-5 Pa.s
Output parameters
mass of water drop m = 2.7x10-10 kg
Reynolds number is in the range from 0 to 1, hence Stokes’ Law  FR  v  is valid.
In the upper atmosphere, water vapour condenses on hygroscopic nuclei

(~ 1010 / m3 radii ~ 10-7 m) forming water droplets whose radii range from about
4x10-6 m to 4x10-5 m. For R < 4x10-5 m, the water droplets very quickly each
their very low values of terminal speed vT < 0.2 m.s-1. Because of the very small
values of the terminal speeds of the water droplets, they can be kept afloat by
vertical air currents. The aggregation of suspended water droplets forms clouds
or near ground level water droplets forms fogs.
In this example, the buoyance force acting on water droplets is insignificant.
Stokes’ law is not valid for raindrops. A raindrop of only 1 mm in radius has a
terminal velocity that is very large vT ~ 120 m.s-1 and NR ~ 10 000.
Intermediate Reynolds Number (300 < NR < 300 000)
For large fast moving objects such as cars, aeroplanes, cricket balls, baseballs,
hailstones, and skydivers, the drag force to a good approximation varies with the
square of its speed. The drag force Fdrag can be written as
Fdrag   1
2 
CD A  F v 2
A effective cross-sectional area of object [m2]

F density of fluid [kg.m-3]
CD drag coefficient [dimensionless]
For objects of a define shape such as spheres, the drag coefficient is nearly
constant over a wide range of speeds and sizes.
For larger Reynold numbers, the turbulence

associated with the motion of the object
through the fluid is mainly responsible for the
frictional retarding force on the object as it
moves through the fluid. This is the reason why
the resistive force depends upon v2 and it is
much greater than the resistive force for
laminar flow where the resistive force depends
upon v.
For an object falling through a fluid, we make the assumption that there are two
distinct forces acting on the object. The buoyant force can be assumed to be
negligible.
weight FG  m g
drag force Fdrag   1
2 
CD A  F v 2
where in our frame of reference, down is taken as the positive direction.
Newton’s Second Law applied to the falling object gives
v
m a  m g   21 CD A  F  v 2  
v
The acceleration a is
v
a  g   21 CD A  F / m  v 2  
v
The terminal velocity vT (a = 0) of the object is
2m g
vT 
CD A  F
We can use the finite difference method to numerically find the velocity v and
displacement x of the object as functions of time.
dv  tk 1  v  tk 2   v  tk 
a  tk 1   
dt 2 t
v  tk  2   v  tk   v  tk 1  
 g   21 CD A  F / m  v 2  
2 t  v  tk 1  
  v  tk 1  
v  tk 2   v  tk    2 t   g   21 CD A  F / m  v 2   
  v  tk 1   
Example Falling raindrops

mscript mec_drag.m
Input parameters
radius of raindrop R = 5.00x10-3 m (big raindrop)

density of water drop  = 1000 kg.m-3
drag coefficient CD = 0.500
Output parameters
mass of raindrop m = 5.24x10-4 kg
cross-sectional area A = 7.85x10-5 m2
time to reach vT t(vT) ~ 4 s
Our model in assuming Fdrag  v 2 should be OK since the Reynolds number
varies from zero to about 10 000.
A raindrop could fall a distance ~ 5000 m. With no frictional forces acting the
raindrops would reach the ground with speeds ~ 300 m.s-1 (> 1000 km.s-1) !!! This
does not happen. Raindrops quickly each their terminal speed within seconds
and reach the ground at speeds < 20 m.s-1.
Example A falling table tennis ball
mscript mec_drag.m mec_tt_ball.m
We can model a table tennis ball falling from rest using the numerical approach
where the resistive force or drag force is
Fdrag   1
2 
CD A  F v 2
Input parameters
radius of table tennis ball R = 20.0x10-3 m

mass of table tennis ball m = 12.70x10-3 kg
mass entered and not calculated from density
Output parameters
cross-sectional area A = 1.26x10-3 m2

distance to reach vT x(vT) ~ 17 m
Our model in assuming Fdrag  v 2 should be OK since the Reynolds number is
mostly in the range from 300 to 300 000.
But, how good is our model?
We can test our model against experimental data. The data for the fall of a table
tennis ball was taken from a paper by French. The data for time and displacement
of the falling table tennis ball was stored in the mscript mec_tt_ball.m and this
mscript was used to compare the video measurements with our displacements /
time predictions using the finite difference approximation to solve the equation of
motion. The following plot shows the French data and the theoretical values for
the displacement as functions of time.
The agreement between our model’s predictions and the actual measured values
is very good.
Example A falling steel ball
mscript mec_drag.m
We can model a steel ball with a radius identical to that of a table tennis ball
falling from rest using the numerical approach where the resistive force or drag
force is
Fdrag   1
2 
CD A  F v 2
We can then compare the fall of the steel ball with that of the table tennis ball.
Input parameters
radius of steel ball R = 20.0x10-3 m

density of steel ball  = 7800 kg.m-3
mass of table tennis ball m = 12.70x10-3 kg
Output parameters table tennis ball steel ball free fall

Fdrag = 0
terminal velocity vT [m.s-1] 8.79 86 
t = 3.00 s a [m.s-2] 0.049 8.75 9.80
t = 3.00 s x [m] 20.9 43.3 44.1
t = 3.00 s v [m.s-1] 8.77 28.3 29.4
t = 3.00 s Reynolds number NR 2.4x104 7.8x10 4
It is “true” that heavier objects do fall faster than lighter ones.
table tennis ball steel ball
Example A falling skydiver
mscript mec_drag.m
We can model a falling skydiver falling from rest

using the numerical approach where the resistive
force or drag force is
Fdrag   1
2 
CD A  F v 2
Input parameters
radius (estimate) R = 0.500 m
mass m = 70 kg
mass entered and not calculated from density
Output parameters
cross-sectional area A = 0.7854 m2
terminal velocity |vT| = 54 m.s-1 = 192 km.h-1
time to reach vT t(vT) ~ 18 s
distance to reach vT x(vT) ~ 760 m
You can investigate how the skydiver can control their speed of fall by varying
their effective cross-sectional area.

Motion of Falling Objects With Resistance: Doing Physics With Matlab Mechanics

Uploaded by

Copyright:

Available Formats

Motion of Falling Objects With Resistance: Doing Physics With Matlab Mechanics

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Motion of Falling Objects With Resistance: Doing Physics With Matlab Mechanics

Uploaded by

Copyright:

Available Formats

DOING PHYSICS WITH MATLAB

MOTION OF FALLING OBJECTS WITH RESISTANCE

DOWNLOAD DIRECTORY FOR MATLAB SCRIPTS

Model (1) FR    v low speeds

Model (2) FR   v 2 high speeds

where  and  are constants of proportionality.

Many problems in the mathematical analysis of particles moving under the

Therefore, the vertical acceleration a of an object through a fluid is

If the object is moving down then v > 0

If the object is moving up then v < 0

The mscript used for the plots is mec_fr_mg_bv.m

where a is the acceleration of the object at any instance.

The initial conditions are t  0 v  v0 x  0 a     / m  v0

When a = 0, the velocity is constant v = vT where vT is the terminal velocity

Plots of the velocity v as functions of time t

Initial values for velocity v0 [m.s-1]

blue: 10 red: vT magenta: 2 cyan: 0

Initial values for velocity v0 [m.s-1]

blue: 10 red: vT magenta: 2 cyan: 0

Initial values for velocity v0 [m.s-1]

blue: 10 red: vT magenta: 2 cyan: 0

Therefore, the equation for the velocity v as a function of time t can be

The maximum height xup reached by the object in time t = tup is

m = 2.00 kg  = 5.00 kg.s-1 g = 9.8 m.s-2 u = 10 m.s-1 vT = 3.92 m.s-1

The time to reach maximum height is tup = 0.507 s

The max height hup reached is hup = 2.013 m xup = - 2.013 m

We can integrate this equation by a substitution method or an algebraic

Note: up is negative and down is positive in our frame of reference.

In the finite difference method we calculate the velocity v and displacement x at

t1, t2, … , tk, …, tN t =t2 - t1 t1 = 0 tk = (k-1) t k = 1, 2, 3, … , N

The acceleration a is approximated by the difference formula

Therefore, the velocity v  tk 2  at time tk+2 is

We estimate the velocity at the second time step t2

v(t2 )  v(t1 )  a(t1 ) t  v(t1 )    / m  v  t1  t

The displacement at each time step is

m = 2 kg  = 10 kg.s-1 v0 = 10 m.s-1 t = 1x10-4 s tmax = 2.00 s

Outputs N numerical approach A analytical approach

terminal velocity vT = 1.96 m.s-1

However, you always need to be careful in using numerical approaches to solve

When the time step is

The mscript used for the plots is mec_fr_mg_bv2.m

where a is the acceleration of the object at any instance.

When a = 0, the velocity is constant v = vT where vT is the terminal velocity

Velocity of the object is always positive (falling object) v0  0 v  0

We can now investigate what happens as time t  

In falling, the object will finally reach a constant velocity vT (a = 0 ) which is

m = 0.010 kg  = 1.0010-4 kg.m-1 v0 = 0 m.s-1  vT = 31.3 m.s-1

m = 0.010 kg  = 1.0010-4 kg.m-1 v0 = +5.00 m.s-1  vT = 31.3 m.s-1

m = 0.010 kg  = 1.0010-4 kg.m-1 v0 = +40.0 m.s-1  vT = 31.3 m.s-1

Equation of motion Note: up is the positive direction

The time tup to reach maximum height occurs when v = 0

The velocity v as a function of time t is

The maximum height xup reached by the object occurs when v = 0

m = 0.010 kg  = 1.0010-4 kg.m-1 v0 = - 12.0 m.s-1  vT = 31.3 m.s-1

The terminal velocity vT is

The time tup to reach the maximum height

From the graphs:

x=0 time t = 2.275 s velocity v = 11.21 m.s-1