Skip to main content
Library homepage
 

Text Color

Text Size

 

Margin Size

 

Font Type

Enable Dyslexic Font
Mathematics LibreTexts

3.5.4: Falling Raindrops

( \newcommand{\kernel}{\mathrm{null}\,}\)

A SIMPLE PROBLEM THAT APPEARS IN MECHANICS is that of a falling raindrop through a mist. The raindrop not only undergoes free fall, but the mass of the drop grows as it interacts with the mist. There have been several papers written on this problem and it is a nice example to explore using numerical methods. In this section we look at models of a falling raindrop with and without air drag.

First we consider the case in which there is no air drag. A simple model of free fall from Newton’s Second Law of Motion is

d(mv)dt=mg

In this discussion we will take downward as positive. Since the mass is not constant. we have

mdvdt=mgvdmdt

In order to proceed, we need to specify the rate at which the mass is changing. There are several models one can adapt.We will borrow some of the ideas and in some cases the numerical values from Sokal(2010)1 and Edwards, Wilder, and Scime (2001).2 These papers also quote other interesting work on the topic.

1

A. D. Sokal, The falling raindrop, revisited, Am. J. Phys. 78, 643-645, (2010).

2

B. F. Edwards, J. W. Wilder, and E. E. Scime, Dynamics of Falling Raindrops, Eur. J. Phys. 22, 113-118, (2001).

While v and m are functions of time, one can look for a way to eliminate time by assuming the rate of change of mass is an explicit function of m and v alone. For example, Sokal (2010) assumes the form

dmdt=λmσvβ,λ>0

This contains two commonly assumed models of accretion:

  1. σ=2/3,β=0. This corresponds to growth of the raindrop proportional to the surface area. Since mr3 and Ar2, then mA implies that ˙mm2/3.
  2. σ=2/3,β=1. In this case the growth of the raindrop is proportional to the volume swept out along the path. Thus, ΔmA(vΔt), where A is the cross sectional area and vΔt is the distance traveled in time Δt

In both cases, the limiting value of the acceleration is a constant. It is g/4 in the first case and g/7 in the second case.

Another approach might be to use the effective radius of the drop, assuming that the raindrop remains close to spherical during the fall. According to Edwards, Wilder, and Scime (2001), raindrops with Reynolds number greater than 1000 and with radii larger than i mm will flatten. Even larger raindrops will break up when the drag force exceeds the surface tension. Therefore, they take 0.1 mm<r<1 mm and 10<Re<1000. We will return to a discussion of the drag later.

It might seem more natural to make the radius the dynamic variable, than the mass. In this case, we can assume the accretion rate takes the form

drdt=γrαvβ,γ>0

Since, m=43πρdr3

dmdtr2drdtm2/3drdt

Therefore, the two special cases become

  1. α=0,β=0. This corresponds to a growth of the raindrop proportional to the surface area.
  2. α=0,β=1. In this case the growth of the raindrop is proportional to the volume swept out along the path.

Here ρd is the density of the raindrop. We will also need

vmdmdt=4πρdr243πρdr3vdrdt=3vrdrdt=3γrα1vβ+1

Putting this all together, we have a systems of two equations for v(t) and r(t) :

dvdt=g3γrα1vβ+1drdt=γrαvβ

Example 3.5.4.1

Determine v=v(r) for the case α=0,β=0 and the initial conditions r(0)=0.1 mm and v(0)=0 m/s.

In this case Equations 3.5.4.2 become

\boldsymbol{\begin{aligned} &\dfrac{d v}{d t}=g-3 \gamma r^{-1} v \\[4pt] &\dfrac{d r}{d t}=\gamma \end{aligned} \end{equation}\label{3.41}}

Noting that

dvdt=dvdrdrdt=γdvdr

we can convert the problem to one of finding the solution v(r) subject to the equation

dvdr=gγ3vr

with the initial condition v(r0)=0 m/s for r0=0.0001 m.

Rearranging the differential equation, we find that it is a linear first order differential equation,

dvdr+3rv=gγ

This equation can be solved using an integrating factor, μ=r3, obtaining

ddr(r3v)=gγr3

Integrating, we obtain the solution

v(r)=g4γr(1(r0r)4)

Note that for large r,vg4γr. Therefore, dvdtg4.

While this case was easily solved in terms of elementary operations, it is not always easy to generate solutions to Equations 3.5.4.2 analytically. Sokal (2010) derived a general solution in terms of incomplete Beta functions, though this does not help visualize the solution. Also, as we will see, adding air drag will lead to a nonintegrable system. So, we turn to numerical solutions.

In MATLAB, we can use the function in raindropf.m to capture the system 3.5.4.2).Hereweputthevelocityin\(y(1) and the radius in y(2).

function dy=raindropf(t,y);

global alpha beta gamma g

dy=[g-3*gamma*y(2)^(alpha-1)*y(1)^(beta+1); ...

gamma*y(2)^alpha*y(1)^beta];

We then use the Runge-Kutta solver, ode45, to solve the system. An implementation is shown below which calls the function containing the system. The value γ=2.5×107 is based on empirical results quoted by Edwards, Wilder, and Scime (2001).

clipboard_ec7521a84e9b11d4d209f16b678f60f7d.png
Figure 3.5.4.1: The plots of position and velocity as a function of time for α=β=0.

clear

global alpha beta gamma g

alpha=0;

beta=0;

gamma=2.5e-07;

g=9.81;

r0=0.0001;

v0=0;

y0=[v0;r0];

tspan=[0 1000];

[t,y]=ode45(@raindropf,tspan,y0);

plot(1000*y(:,2),y(:,1),’k’)

clipboard_e0773d48afbb5e96a3daa2e11f6eea3f5.png
Figure 3.5.4.2: The plot the velocity as a function of position for α=β=0.
clipboard_ef53f6aa9f8cdf91af413605cd87674e2.png
Figure 3.5.4.3: The plot the velocity as a function of position for α=0,β=1.

The resulting plots are shown in Figures 3.5.4.1-3.5.4.2. The plot of velocity as a function of position agrees with the exact solution, which we derived in the last example. We note that these drops do not grow much, but they seem to obtain large speeds.

For the second case, α=0,β=1, one can also obtain an exact solution. The result is

v(r)=[2g7γr(1(r0r)7)]12

For large r one can show that dvdtg7. In Figures 333332 we see again large velocities, though about a third as fast over the same time interval. However, we also see that the raindrop has significantly grown well past the point it would break up.

In this simple model of a falling raindrop we have not considered air drag. Earlier in the chapter we discussed the free fall of a body with air resistance and this lead to a terminal velocity. Recall that the drag force given by

fD(v)=12CDAρav2

where CD is the drag coefficient, A is the cross sectional area and ρa is the air density. Also, we assume that the body is falling downward and downward is positive, so that fD(v)<0 so as to oppose the motion.

We would like to incorporate this force into our model 3.5.4.2. The first equation came from the force law, which now becomes

mdvdt=mgvdmdt12CDAρav2

Or

dvdt=gvmdmdt12mCDAρav2

The next step is to eliminate the dependence on the mass, m, in favor of the radius, r. The drag force term can be written as

fDm=12mCDAρav2=12CDπr243πρdr3ρav2=38ρaρdCDv2r

We had already done this for the second term; however, Edwards, Wilder, and Scime (2001) point to experimental data and propose that

dmdt=πρmr2v

where ρm is the mist density. So, the second terms leads to

vmdmdt=34ρmρdv2r

clipboard_e4e4c82730916c154b3a87f8f3288491b.png
Figure 3.5.4.4: The plots of position and velocity as a function of time for α=0,β=

But, since m=43πρdr3

dmdt=4πρdr2drdt

SO,

drdt=ρm4ρdv

This suggests that their model corresponds to α=0,β=1, and γ=ρm4ρd.

Now we can write down the modified system

dvdt=g3γrα1vβ+138ρaρdCDv2rdrdt=γrαvβ

Edwards, Wilder, and Scime (2001) assume that the densities are constant with values ρa=.856 kg/m3,ρd=1.000 kg/m3, and ρm=1.00×103 kg/m3. However, the drag coefficient is not constant. As described later in Section 3.5.7, there are various models indicating the dependence of CD on the Reynolds number,

Re=2rvv

where v is the kinematic viscosity, which Edwards, Wilder, and Scime (2001) set to v=2.06×105 m2/s. For raindrops of the range r=0.1 mm to 1 mm, the Reynolds number is below rooo. Edwards, Wilder, and Scime (2001) modeled CD=12Re1/2. In the plots in Section 3.5.7 we include this model and see that this is a good approximation for these raindrops. In Chapter 10 we discuss least squares curve fitting and using these methods, one can use the models of Putnam (1961) and Schiller-Naumann (1933) to obtain a power law fit similar to that used here. So, introducing

CD=12Re1/2=12(2rvv)1/2

and defining

δ=923/2ρaρdv1/2

we can write the system of equations 3.5.4.6 as

dvdt=g3γv2rδ(vr)32drdt=γv.

Now, we can modify the MATLAB code for the raindrop by adding the extra term to the first equation, setting α=0,β=1, and using δ=0.0124 and γ=2.5×107 from Edwards, Wilder, and Scime (2001).

clipboard_e8729b10254bde778386416d5e21782cb.png
Figure 3.5.4.5: The plots of position and velocity as a function of time with air drag included.
clipboard_e438058eae7b1269cefdc62aa05634249.png
Figure 3.5.4.6: The plot the velocity as a function of position with air drag included.

In Figures 3.5.4.5- 3.5.4.6 we see different behaviors as compared to the previous models. It appears that the velocity quickly reaches a terminal velocity and the radius continues to grow linearly in time, though at a slow rate.

We might be able to understand this behavior. Terminal, or constant v,

g3γv2rδ(vr)32=0

Looking at these terms, one finds that the second term is significantly smaller

δ(vr)32g

vr(gδ)2/385.54 s1

This agrees with the numerical data which gives the slope of the v vs r plot would occur when than the other terms and thus

Or as 85.5236 s1.


This page titled 3.5.4: Falling Raindrops is shared under a CC BY-NC-SA 3.0 license and was authored, remixed, and/or curated by Russell Herman via source content that was edited to the style and standards of the LibreTexts platform.

  • Was this article helpful?

Support Center

How can we help?