Reanalyze the data from Example E9.11, concerning Millikan's oil-drop experiment, to use a more accurate approximation for the effective air viscosity:
$$
\eta = \frac{\eta_0}{1 + \frac{b}{ap}},
$$
where $p = 100.82\;\mathrm{kPa}$ is the air pressure, $\eta_0 = 1.869 \times 10^{-5}\;\mathrm{kg\,m^{-1}\,s^{-1}}$, $b = 7.88\times 10^{-3}\;\mathrm{Pa\, m}$, and $a$ is the droplet radius.
Solution P9.4.2
Substituting the more accurate expression for $\eta$ into the expression $-m'g - 6\pi\eta a v_g$ gives:
$$
-\frac{4}{3}\pi a^3\rho'g = \frac{6\pi\eta_0av_g}{1+\frac{b}{ap}},
$$
which rearranges to a quadratic equation in the droplet radius, $a$:
$$
a^2 + \frac{b}{p}a + \frac{9\eta_0v_g}{2\rho'g} = 0.
$$
Taking the positive root gives the following expression for $a$ (recall that $v_g < 0$):
$$
a = -\frac{b}{2p} + \frac{1}{2}\sqrt{\frac{b^2}{p^2} - \frac{18\eta_0v_g}{\rho'g}},
$$
which is used in the function get_a below.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
df = pd.read_csv("eg10-millikan-data.txt", sep=r"\s+", index_col=[0, 1])
# Air viscosity parameters, eta0 (kg.m-1.s-1) and b (Pa.m).
eta0 = 1.869e-5
b = 7.88e-3
# Air density, oil density, density difference (kg.m-3).
rho_air = 1.17
rho_oil = 917.3
rhop = rho_oil - rho_air
# Air pressure (Pa), acceleration due to gravity (m.s-2).
p = 100.82e3
g = 9.803
# Fall/rise distance (m), electric field (V.m-1).
d = 11.09e-3
E = -322.1e3
def get_eta(a):
"""Return the viscosity adjusted for droplet radius."""
return eta0 / (1 + b / a / p)
def get_a(tg):
"""Get the droplet radius from fall time."""
# Droplet fall terminal velocity.
vg = -d / tg
c = b / p
a = (-c + np.sqrt(c**2 - 18 * eta0 * vg / rhop / g)) / 2
return a
for drop in df.index.levels[0]:
drop_df = df.loc[drop].T
# Rename the index to simple 'tg' and 'te'.
drop_df.index = drop_df.index.str.slice(0, 2)
# Average over all tg values.
tg = drop_df.loc["tg"].values.mean()
# Average te for each experiment.
te = drop_df.loc["te"].mean()
a = get_a(tg)
eta = get_eta(a)
q = 6 * np.pi * eta * a * d / E * (1 / tg + 1 / te)
df.loc[drop, "q"] = q.values
sorted_q = sorted(-df.loc[:, "q"])
dq = np.diff(sorted_q)
e_estimate = dq[(dq > 1.0e-19) & (dq < 2.0e-19)].mean()
df["N"] = (df["q"] / e_estimate).astype(int)
e = (df["q"] / df["N"]).mean()
print(f"Estimated value for e = {e} C")