Phase 1
This phase corresponds to the rotational motion of the toast about \(O\). We take advantage of conservation of energy to find that angle \(\theta\) is given by \[ \dot{\theta}^2 = \omega^2 \sin\theta , \qquad \omega^2 = \frac{6gd}{b^2 + 3d^2} \] The friction and normal force are given by \[ F = mg B \sin\theta, \quad N= mg D \cos\theta \] with \[ B= \frac{b^2 +9d^2}{b^2 +3d^2},\qquad D= \frac{b^2}{b^2 +3d^2} \] This phase lasts as long as \(F\leq \mu N\), that is, \(\theta \leq \theta_1 = \tan^{-1}(\mu D/B)\). We can then easily find \(\theta\) vs \(t\) or rather \(t\) vs \(\theta\) by using the following quadrature: \[ t = \frac{1}{\omega} \int_0^\theta \frac{du}{\sqrt{\sin u}}, \qquad 0< \theta \leq\theta_1 \qquad (1) \] This gives \(t=t_1\) by setting \(\theta=\theta_1\).
Python Coding: This is straightforward to code by using the function quad of the module
scipy.integrate which needs to be loaded. We need to define a function to define \(f(x)=1/\sqrt{\sin{x}}\). Prior to finding time \(t_1\), we define the constants \(B\), \(D\), \(\omega\), \(\theta_1\), \(\omega_1 =\dot{\theta}_1\). To plot \(\theta\) vs \(t\), we define 2 arrays according to
equation (1).
## Import the necessary libraries/modules
import numpy as np
import matplotlib.pyplot as pypl
from scipy.integrate import quad
# We define parameters: b, d, h, g and mu
b=0.04 ; # toast half length [m]
d= 0.004 ; # offset length [m]
h=0.75 ; # table height [m]
g= 9.81 ; # gravitational acceleration [m/s^2]
mu= 1; # static/kinetic friction coefficient
###################### PHASE 1 ##############################
# Rotational phase: we take advantage of the conservation of
# energy to reduce the problem to a simple quadrature.
# First we solve for angle theta1, the incipient slip angle, the
# corresponding value of omega=omega1 and the end time t=t1
##############################################################
B= 1+ (6*d**2)/(b**2 + 3*d**2)
D= 1+ (-3*d**2)/(b**2 + 3*d**2)
theta1= np.arctan((mu*D)/B)
print ("theta1 = ",180*theta1/np.pi," degrees")
# conservation of energy gives theta_dot^2 = omega^2 * sin(theta)
# with omega:
omega = np.sqrt(6*g*d/(b**2 +3*d**2))
# find omega1 = theta_dot at t=t1 [rad/s]
omega1= omega*np.sqrt(np.sin(theta1))
print ("omega1 = ",omega1," rad/s")
# find time t=t1 such that theta(t1)=theta1
def NoSlip(theta):
return 1/(omega*np.sqrt(np.sin(theta)) )
time1, err1 = quad(NoSlip, 0, theta1);
print ("time1 = ",time1," sec"," error= ",err1)
# solve for time t, 0 < t <t1 such that t= (1/omega)* int_0^theta d phi/sqrt(sin(phi))
Theta1 = np.linspace(0, theta1, 100)
Time1= np.linspace(0, time1, 100)
for i in range(1, 100):
Time1[i], err1 = quad(NoSlip, 0, Theta1[i])
if( err1 > 0.0001):
print('warning! error = ',err1)
Theta1[i]= 180*Theta1[i]/np.pi
pypl.plot(Time1,Theta1) # plot theta vs time (PHASE 1)Phase 2: Slip Phase
To compute the motion during phase 2, we integrate the ODEs governing \((r,\theta)\):
\[
\begin{cases}
m (\ddot{r}-r\dot{\theta}^2) =-\mu N + mg \sin\theta\\
m (2\dot{r}\dot{\theta}+ r \ddot{\theta})= -N + mg \cos\theta\\
\frac{1}{3} mb^2 \ddot{\theta}= r N
\end{cases}
\qquad (2-4)
\]
with initial conditions \(r=d\), \(\dot{r}= 0\), \(\theta= \theta_1\), \(\dot{\theta}= \omega_1\). To do this numerically we use the function odeint of the module scipy.integrate. We must the system of ODE as a 4-dimensional system:
\[
\begin{cases}
\dot{r}= u \\
\displaystyle\dot{u} = r\omega^2 +g \sin\theta-{\mu b^2} \frac{g \cos\theta-2 u \omega}{3r^2+ b^2}\\
\dot{\theta}= \omega\\
\displaystyle \dot{\omega}= 3r \frac{g\cos\theta-2 u \omega}{3r^2+ b^2}
\end{cases}
\qquad (5)
\]
To find the time duration \(t=t_2\) of phase 2, we need to find the value
of \(t\) for which \(N=0\). This is equivalent to finding the
value of \(t\) for which \(\ddot{\theta}=0\) according to equation (4). We could implement a test to monitor the sign of \(\ddot{\theta}\) in the system of equations (5). This is clumsy. Instead, we define a separate function to compute \(\ddot{\theta}\) and build an interpolation of \(t\) vs \(\ddot{\theta}\) computed from arrays of \(t\) and \(\ddot{\theta}\) values.
We use function interp1d from the scipy.interpolate module: given two arrays \(x\) and \(y\), it returns a function \(y= f(x)\) by using linear interpolation. The restriction is that \(x\) must be an array of monotonically increasing or decreasing values.
This easily yields the value of \(t_2\) by seeking the value of the interpolating function at zero. We proceed the same way to find the value \((r_2, \dot{r}_2, \theta_2,\dot{\theta}_2)\) of
\((r,\dot{r},\theta, \dot{\theta})\) at time \(t_2\). We then conclude by plotting \(\theta\) vs \(t\) for phase 2.
import numpy as np
import matplotlib.pyplot as pypl
from scipy.integrate import odeint
from scipy.interpolate import interp1d
# define the system of 1st-order ODEs in the variables (r,rdot,theta,thetadot)
def SlipPhase(y, t):
[r, rdot, theta, thetadot] = y
Ftheta= (-2*rdot*thetadot +g*np.cos(theta))/(3.0*r**2+b**2)
G= r*thetadot**2 +g*np.sin(theta)
return [rdot, G-mu*b**2*Ftheta, thetadot, 3*r*Ftheta]
tmax=0.15 # [s]: max time of integration.
# This may need to be adjusted
# choose tmax sufficiently large so that theta_ddot changes sign
n=2000 # number of integration points
t = np.linspace(0, tmax, n)
values = odeint(SlipPhase, [d, 0, theta1, omega1], t)
r=values[:,0]
rdot=values[:,1]
theta=values[:,2]
thetadot= values[:,3]
# to find the duration t= t2 of phase2 we need to find the value
# of t for which N=0. This is equivalent to finding the
# value of t for which theta_ddot=0
# first we create a function which return theta_ddot
def ThetaDDot(R,Rdot,Theta,ThetaDot):
return (-2*Rdot*ThetaDot +g*np.cos(Theta))/(3.0*R**2+ b**2)
# we create array thetaddot from the arrays r,rdot,theta,thetadot
thetaddot= ThetaDDot(r,rdot,theta,thetadot)
# we interpolate t vs thetaddot
TvsThetaddot = interp1d(thetaddot,t) # interpolation of Y=time vs X=theta_ddot
# this create a function t= t(theta_ddot)
# we find t2= time2 and theta2=theta(t2), omega2= theta_dot(t2)
# r2= r(t2) and rdot2= r_dot(t2)
time2=TvsThetaddot(0) # time t2 is found by setting theta_ddot=0
TvsTheta = interp1d(t,theta) # interpolation of X=time vs Y=theta
TvsThetaDot = interp1d(t,thetadot) # interpolation of X=time vs Y=theta_dot
TvsR = interp1d(t,r) # interpolation of X=time vs Y=r
TvsRdot = interp1d(t,rdot) # interpolation of X=time vs Y=rdot
theta2= TvsTheta(time2)
omega2= TvsThetaDot(time2)
r2=TvsR(time2)
rdot2=TvsRdot(time2)
print ("time2 = ",time2," s ")
print ("theta2 = ",180*theta2/np.pi," degrees")
print ("omega2 = ",omega2," rad/s")
Theta2 = np.linspace(theta1, theta2, 100)
Time2= np.linspace(0, time2, 100)
for i in range(0, 100):
Theta2[i] = 180*TvsTheta(Time2[i])/np.pi
Time2[i]=time1+Time2[i]
pypl.plot(Time2,Theta2) # plot theta vs timePhase 3: Free Fall
To compute the motion during phase 3, we use the standard parabolic motion of mass center \(G\) and the fact that \(\ddot{\theta}=0\): this yields the following system of equations
\[
\begin{cases}
x_G = (\dot{r}_2 \cos\theta_2 -r_2 \omega_2 \sin\theta_2) t + r_2 \cos\theta_2
\\
\displaystyle y_G = g \frac{t^2}{2} + (\dot{r}_2 \sin\theta_2 +r_2 \omega_2 \cos\theta_2) t + r_2 \sin\theta_2
\\
\theta= \theta_2 +\omega_2 t
\end{cases}
\qquad(6-8)
\]
The toast hits the floor when one of its two edges reaches the level \(y=h\) of the floor: the
corresponding time
\(t=t_3\) is found by solving the equation
\[
\max\Big(g \frac{t^2}{2}+ (\dot{r}_2 \sin\theta_2 +r_2 \omega_2 \cos\theta_2) t + r_2 \sin\theta_2 \pm b \sin(\omega_2 t+\theta_2) \Big) = h\qquad (9)
\]
Since this is a nonlinear equation, we use the nonlinear solver brentq of the module
scipy.optimize. This requires to define a function according to equation (9). Once \(t_3\)
is found, we can find angle \(\theta_3 =\theta(t_3)=\theta_2 +\omega_2 t_3\). If angle \(\theta_3\) satisfies the condition
\({ \pi/2} < \theta_3 < {3\pi}/{2}\), the toast will hit the floor with the butter side down. The final step is to plot \(\theta\) vs \(t\) according to equation (8).
############### PHASE 3 ##################################
### to study phase 3, we use the parabolic motion of mass
### center G and the linear evolution of angle theta= omega2*t+theta2.
#############################################################
from scipy import optimize
# define position of G at t=t2
xGinit= r2*np.cos(theta2)
yGinit= r2*np.sin(theta2)
vGx= rdot2*np.cos(theta2) - r2*omega2*np.sin(theta2)
vGy= rdot2*np.sin(theta2) + r2*omega2*np.cos(theta2)
# define the path of G
def xG(T):
return vGx*T+xGinit
def yG(T):
return 0.5*g*T**2 + vGy*T + yGinit
# define the orientation of toast
def theta_free(T):
return theta2 + omega2*T
# find the time of impact: first define vertical coordinate
# of edges
def Target(T):
yTip1= yG(T)+ b*np.sin(theta_free(T))
yTip2= yG(T)- b*np.sin(theta_free(T))
return max(yTip1,yTip2) -h
# solve for t=time3 by solving for the root of function Target
time3 = optimize.brentq(Target,0.001,1.5)
theta3= theta_free(time3)
print ("time3 = ",time3," sec")
print ("theta3 = ",180*theta3/np.pi,"degrees")
# decide whether landing is on butter side up or down
if(theta3 < 0.5*np.pi) or (theta3 > 1.5*np.pi):
print ("Landing: butter side up")
if(theta3 > 0.5*np.pi) and (theta3 < 1.5*np.pi):
print ("Landing: butter side down")
# define time and theta arrays for phase 3
Theta3 = np.linspace(theta2, theta3, 100)
Time3 = np.linspace(0, time3, 100)
for i in range(0, 100):
Theta3[i] = 180*theta_free(Time3[i])/np.pi
Time3[i]= time1+ time2 + Time3[i]
pypl.plot(Time3,Theta3) # plot theta vs time
# plot boundary points between 3 phases
t_list = [0,time1,time2+time1,time1+time2+time3]
theta_list = [0,180*theta1/np.pi,180*theta2/np.pi,180*theta3/np.pi]
pypl.scatter(t_list, theta_list, s=10)
pypl.xlabel(r'$t$',fontsize=18)
pypl.ylabel(r'$theta$',fontsize=18)
pypl.grid()
pypl.title('mu = 1.0')
pypl.show()
pypl.savefig("toastplot-mu1.svg")
pypl.close()Results
Study 1: We use the numerical values \(b =4\) cm, \(d=4\) mm, $h =0.75 $ m and \(\mu =1\). Here is the numerical output of the code:
theta1 = 42.5 degrees
omega1 = 9.83 rad/s
time1 = 0.145 sec error= 3.623451538814493e-11
time2 = 0.041 s
theta2 = 67.5 degrees
omega2 = 11.23 rad/s
time3 = 0.362 sec
theta3 = 300.8 degrees
Landing: butter side upWe conclude that the toast lands butter-side up. The evolution of angle \(\theta\) is shown in Figure 1.
Study 2: We fix the parameters to the values \(b =4\) cm, \(d=4\) mm and \(h =0.75\) m, and we vary parameter \(\mu\) to find the range for which the toast lands butter-side down. We find the threshold to be reached at about \(\mu \approx 0.44\), that is, if \[ \mu_{min}\approx 0 \leq \mu \leq \mu_{max}\approx 0.44, \qquad \text{the toast lands butter-side down} \] \[ \mu > \mu_{max}\approx 0.44, \qquad \text{the toast lands butter-side up} \] Here are a few notable features:
The transition from “butter-side-down” to “butter-side-down” takes place when \(\theta_3 = 270^o\) is reached.
As \(\mu\) is increased, the duration of the rotational phase increases while the slip phase vanishes.
Conversely, as \(\mu \to 0\), the duration of the rotational phase vanishes while the duration of the slip phase increases.
The free fall phase remains sensibly constant in duration.
The dynamics of the toast is shown in Figure 2 for a few values of \(\mu\).
The full script is available here.