% SIMIODE-TeacherStudent-Template.tex and simiode.cls for LaTeX2e
% Version 22 November 2016.
% 
% This template is to be used with the  class file ``simiode.cls". 
%  
% The latest version of this class and template files 
% can be found at
% 
%       <http://www.SIMIODE.org> in the Author section.
%
% A ``Quick Start Guide" appears below.
%
% This file must be processed with version LaTeX2e or higher.
% Appreciation to Mark Frantz, IUPUI; John Thoo, Yuba College; and Eric Sullivan, Carroll College.
%
\documentclass{simiode}
%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario only with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resourcing, but nothing in the bibliography which would ``give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with   TEACHER VERSION, but also full author  identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation  ideas, ``solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review, but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.   
% Rename the pdf after you compile each final version.
%
\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
%\usepackage{hyperref}
\usepackage{url}
%
\renewcommand{\baselinestretch}{1.3}
%
% NOTE regarding macros %%
%
% Any macros defined for your paper should be contained in 
% the top matter. Likewise, any environment definitions such 
% as \newtheorem or \newenvironment should also go in the 
% top matter. 
% 
% Do not \input your macros as separate files within your 
% final version because more files make it harder for the 
% editors and users to keep track of and/or modify and 
% append your paper.
%
% Top matter is above and the beginning of the document below.
%
\begin{document}

\newpage
%
% You only need to change the NAME OF SCENARIO in the next line and delete the footnote. 
% See above for instructions on toggling from Teacher Version to Student Version.
%
\title{\ifnum\TeacherEdition=1 TEACHER VERSION \else STUDENT VERSION \fi\\Solid Particle Erosion}




\markboth{Solid Particle Erosion}{Solid Particle Erosion}

\author{Rich Laverty\\
Philadelphia, PA\\
rich.laverty@gmail.com
%Philadelphia Design Center\\ 
%The Boeing Company\\
%Ridley Park, PA \teacher{19078} USA\\
%\teacher{richard.r.laverty@boeing.com}
}

\makeStitlePDFLaTex


\begin{abstract}
By applying Newton's second law, and making a collection of reasonable assumptions, students will derive a system of differential equations that model the path of a rigid particle as it gouges material from a more ductile surface.  Examination of the solution will yield a formula for the volume of the gouge, that is, the amount of surface material eroded from a single impact.
\end{abstract}


\section*{SCENARIO DESCRIPTION}
The process of erosion is familiar to most of us from Earth science, where water or wind removes material from one location and transports it to another.  However, it is also an important industrial process.  We can think of erosion more generally as solid particles, entrained by a fluid, attacking and removing material from a surface.  This process can adversely affect the performance and life of such equipment as pipes, turbines, and helicopter rotorblades.  It is also the cause of the corona effect, sometimes called the Kopp-Etchells effect, experienced by helicopters where small pieces of eroded metal can oxidize and illuminate the rotorblades - a dramatic, but obviously undesirable condition for military operations.  On the other hand, erosion can have useful applications such as sand blasting, abrasive deburring, and drilling hard materials.  For a process of such significance, both good and bad, an understanding of the basic physical processes involved is required to take measures to either mitigate or take advantage of its effects.

In this project we will develop a mathematical model of the cutting process by which a relatively rigid particle removes material from a ductile surface.  The canonical situation would be a piece of sand striking a metal surface.  The mathematical model will be in the form of a system of differential equations resulting from application of Newton's Second Law of Motion.  Our model is based on the first published attempt at an analytical description of solid particle erosion \cite{Finnie1958}.  From the solution we will be able to determine the conditions for maximum material removal, and the role played by the different physical parameters of the model.  In particular, we will be interested in the angle of attack at which maximum material removal occurs.

\subsection*{PHYSICS BACKGROUND}
The physics knowledge required to formulate this problem is covered early in any calculus textbook aimed at science and engineering students.  So the concepts should be review, but taking the time to go through the details will be time well spent.

We start with kinematics, the subfield of mechanics that describes motion without concern for the causes of that motion.  For bodies that are rigid and large enough to ignore quantum effects, kinematics amounts to accounting for the position and orientation (and their rates of change) of the body in space.  Position is typically given as a vector specifying the coordinates of the Center of Mass (CM) of the body with respect to some fixed coordinate system.  Orientation, also a vector, is most often quantified by specifying the rotation of a body about coordinate axes whose origin moves with the CM.

We will be idealizing our physical world with two dimensions, so the position vector will be a 2-tuple.  If the body is to remain in the plane, orientation will be reduced to a single number - the angle of rotation about an axis perpendicular to the plane of the body and passing through the CM. Therefore, the kinematics of our rigid erodent (the body causing the material removal) will be quantified with just three numbers, $\left(x, y, \phi \right)$, two representing translation in the plane, and the third describing its rotation about the CM, see Figure \ref{fig:coordinate_system}.

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{coordinate_system.jpg}
\caption{The position $\bar{r}$ of a rigid body in two dimensions and it's orientation $\phi$ with respect to an axis through the CM and perpendicular to the page}
\label{fig:coordinate_system}
\end{figure}

The second subfield of mechanics is kinetics.  Kinetics relates applied forces to the motion they cause.  Newton's famous three laws of motion are the foundation of Kinetics.  Specifically, Newton's Second Law states that the vector sum of the applied forces acting on a particle is equal to the particle's rate of change of momentum.  For many problems the mass of the particle under consideration is constant (or close enough), so that the rate of change of momentum is equal to mass times acceleration.

\begin{equation}
\sum\bar{F} = \frac{d}{dt}\left(m\bar{v}\right) = m\bar{a}
\end{equation}

The direct application of Newton's Second Law will give us a vector differential equation relating the acceleration of the CM to the forces acting on the body.  If we let $F_x$ and $F_y$ denote the sum of the force components in the $x$ and $y$ directions, respectively, then the component equations are:  

\begin{eqnarray}
F_x & = & m\ddot{x}\label{eqn:x_eqn}\\
F_y & = & m\ddot{y}\label{eqn:y_eqn}
\end{eqnarray}

\noindent where overdots are used as a short-hand for time derivatives, in this case the second derivative.

To complete the model we need to supplement these equations with a relationship that allows us to solve  for the orientation of the body.  Toward this end, consider a coordinate system anchored at the rigid body's CM, and partition the rigid body into a finite number of subsets.  Consider the external force $d\bar{F}_i$ applied to the $i^{th}$ subset of the partition.  Let the mass of of the $i^{th}$ subset be denoted $dm_i$ (we will use $dm_i$ to refer to the $i^{th}$ subset itself, and its mass).  First, consider $dm_i$ as a particle and apply Newton's Second Law, $d\bar{F}_i = dm_i\bar{a}_i$.  Now, take the cross product of each side of this equation with the position vector $\bar{r}_i$ that locates $dm_i$ with respect to the coordinate system at the CM.

\begin{equation}
\bar{r}_i\times d\bar{F}_i = \bar{r}_i\times dm_i\bar{a}_i
\end{equation}

Since the body is rigid, the subset $dm_i$ will remain at a constant distance from the CM.  As the rigid body rotates its orientation changes.  From the point of view of the CM, this will look like the subset $dm_i$ is moving along a circular path centered at the CM, see Figure \ref{fig:partition}.  This suggests a convenient decomposition for the vectors $d\bar{F}_i$ and $\bar{a}_i$ with components normal, and tangent to the trajectory of $dm_i$.  We denote the components with superscripts $n$ and $t$, see Figure \ref{fig:decomposition}.

\begin{figure}[h]
\centering
\includegraphics[width=0.5\textwidth]{rigid_body_partitioned.jpg}
\caption{The force $d\bar{F}_i$ applied to the subset $dm_i$ of the partitioned rigid body.  As the body rotates, the subset $dm_i$ will remain at a fixed distance from the CM; that is, viewed from the CM, it will travel in a circle.}
\label{fig:partition}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{decompositions.jpg}
\caption{The components of force and acceleration that are normal, and parallel to the position vector from the CM to the subset $dm_i$.}
\label{fig:decomposition}
\end{figure}

The computational convenience of this decomposition is that the cross product of parallel vectors is zero.  Since the normal components are always parallel to $\bar{r}_i$, their cross products with $\bar{r}_i$ vanish.

\begin{eqnarray}
\bar{r}_i\times \left( d\bar{F}_{i}^{n} + d\bar{F}_{i}^{t} \right) & = & \bar{r}_i\times dm_i\left( \bar{a}_{i}^{n} + \bar{a}_{i}^{t} \right)\\
\bar{r}_i\times d\bar{F}_{i}^{t} & = & \bar{r}_i\times dm_i \bar{a}_{i}^{t}
\end{eqnarray}

This is still a vector equation; however, the vectors that we are crossing lie the plane of the two dimensional rigid body.  We know that the direction of their cross product will be perpendicular to the page.  With the direction known (and unchanging throughout the motion) we need only find the magnitude of the vectors involved.  We will use the notation that a vector's magnitude is denoted with the same symbol as the vector, but without the overbar.  In general, the magnitude of a cross product is $\left|\bar{a}\times\bar{b}\right| = ab\sin\theta$, where $\theta$ denotes the angle between the vectors $\bar{a}$ and $\bar{b}$.  Due to our choice of decomposition, the cross products only involve vectors at right angles, therefore since $\sin\left(90^{\circ}\right) = 1$ their magnitudes satisfy

\begin{equation}
r_i dF_{i}^{t} = r_i dm_i a_{i}^{t}
\end{equation}

The product on the left hand side of this equation is commonly known as {\em torque}, and we will denote it by $d\tau_i$.  On the right hand side we have the tangential acceleration.  This is the rate of change of the magnitude of $dm_i$'s velocity vector (the normal component of the acceleration, which does not play a role for our analysis, is due to the change in the {\em direction} of $dm_i$'s velocity vector).  With respect to the CM-based coordinate system, $dm_i$ is traveling in a circle, denoting its angular displacement by $\phi_i$ the magnitude of its velocity is $r_i\dot{\phi}_i$.  Furthermore, since $r_i$ is constant

\begin{equation}
d\tau_i = r_i dm_i \frac{d}{dt}\left(r_i \dot{\phi}_i\right) = r_{i}^{2}dm_i\ddot{\phi}_i
\end{equation}

There is one more point to make about this equation for the subset $dm_i$.  Since the body is rigid, each point will rotate with the same angular displacement, thus $\phi_i$ is the same for all $i$.  So we can drop the subscript and refer to $\phi$ as the angular displacement, or rotation, of the entire body.

Now, sum over the partition:

\begin{equation}
\sum\limits_i d\tau_i = \ddot{\phi}\sum\limits_i r_{i}^2 dm_i
\end{equation}

\noindent As we refine the partition $i\to\infty$ and $dm_i\to 0$, and the sums become integrals over the rigid body

\begin{equation}
\int\limits_{body}d\tau = \ddot{\phi}\int\limits_{body}r^2 dm
\end{equation}

The integral on the right hand side of the equation has a special name, the {\em mass moment of inertia} of the body, typically denoted by $I$.  It encodes not just the mass of the rigid body, but how that mass is distributed with respect to the origin of the coordinate system.  We can simplify the last expression by substituting $\tau$ for the total torque applied to the body, and $I$ for the mass moment of inertia:

\begin{equation}
\tau = I\ddot{\phi}\label{eqn:phi_eqn}
\end{equation}

The three equations are now established.  The first two are the component equations of Newton's Second Law, relating the $x$ and $y$ acceleration of the CM to the forces applied to the rigid body.  The third equation relates the rotation to the applied torque.  Now we are prepared to model the trajectory of a rigid erodent attacking a ductile surface.
 
\subsection*{EQUATIONS OF MOTION}
To estimate the amount of material removed from the surface, we will use the trajectory of the cutting tip of the erodent.  Since the erodent is assumed rigid, the position of the cutting tip is calculable from the knowledge of the $x$ and $y$ coordinates of the CM, the tip's distance $r$ from the CM, and the body's rotation $\phi$.  Therefore, we can determine the trajectory for the tip by solving the differential equations (\ref{eqn:x_eqn}, \ref{eqn:y_eqn}, \ref{eqn:phi_eqn}) for $x$, $y$, and $\phi$.  We only need to account for all of the forces (and torques) applied to the erodent during its motion.

Figure \ref{fig:initial_conditions} shows an arbitrary erodent at the moment of contact with the ductile surface.  The ductile surface is assumed stationary - even though in practice it may be in motion, e.g. turbine blades.  By assuming the ductile surface is stationary we are observing the impact as if we were located on the ductile surface.  The angle of the impact $\alpha$ (the angle of attack) is measured from the horizontal; a grazing impact would have $\alpha = 0$, and $\alpha = \pi/2$ for a normal impact.  The speed of the impact is $v_0$.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{initial_conditions.jpg}
\caption{The initial conditions of the impact of the erodent upon the ductile surface.}
\label{fig:initial_conditions}
\end{figure}

The only force acting on the erodent is contact with the surface material, see Figure \ref{fig:force_on_erodent}.  Forces due to electromagnetic interactions, air resistance, gravity, etc. are all negligible compared to the contact force.  To simplify the analysis we assume that the surface material is perfectly plastic, so that when the surface material is disturbed it is in a state of hydrostatic stress.  That is, when the erodent is cutting through the perfectly plastic surface material, the surface material will exert a constant pressure on the erodent that will act normal to the contact area between the two materials.  This pressure is called the {\em flow stress}.  It is a material property of the surface material and we will denote it by $\sigma$.  Therefore, the force acting on the erodent has magnitude $\sigma A$, where $A$ is the area of contact between the materials; the direction of the force is normal to the contact surface and opposes the motion.  To simplify our derivation we will assume that the erodent has a uniform thickness $b$ out of the plane of motion.  Therefore, the contact area is $A = bh$, see Figure \ref{fig:contact_area}.  The magnitude of the total force on the erodent is $F = \sigma bh$.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{force_on_erodent_CM.jpg}
\caption{The contact force acting on the erodent opposes the motion and is normal to the contact surface.}
\label{fig:force_on_erodent}
\end{figure}

We first derive an expression for the $x$ component of the force.  Consulting with Figure \ref{fig:contact_area}, $F_x = -\sigma bl$; the minus sign on the right hand side is consistent with our coordinate system.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{contact_area.jpg}
\caption{The contact area between the erodent and surface is $A = bh$, therefore the magnitude of the contact force is $F = \sigma bh$. When considering the $x$-component only $F_x = -\sigma bl$ .}
\label{fig:contact_area}
\end{figure}

Ahead of the erodent the surface material builds up into a lip so that the contact area between the erodent and surface extends above the undisturbed height of the surface, Figure \ref{fig:contact_area}.  We will assume that the height $l$, which takes into account this lip, is proportional to the depth of the tip $y_{tip}$ with constant of proportionality $\psi$.  Thus, $F_x = -\sigma b\psi y_{tip}$, and so

\begin{equation}
m\ddot{x} = -\sigma b\psi y_{tip}
\end{equation}

This is unsatisfactory because it introduces a new variable, $y_{tip}$, but provides no additional equations.  Therefore, to solve our system we will need to find an additional equation, or relate $y_{tip}$ to the preferred kinematic variables $x$, $y$, and $\phi$.  The most convenient approach is to replace $y_{tip}$ with the other kinematic variables through a careful examination of their meaning.  First, we need to distinguish between {\em translation} and {\em position}.  The position of any point of the body is its location in the plane, it is the composition of a translation and rotation.  Second, in a rigid body all points experience the same translation (each point just has a unique starting point), and the common translation is $x$ and $y$ of the CM.  So to find the position $y_{tip}$ we find $y$ starting from the tip's initial position (the origin), and then apply the rotation about the CM.  Note that we use the initial condition $\phi\left( 0 \right) = \pi/2$ to evaluate the contribution from rotation; this will be explained when we consider the equation for $\phi$.  Thus, we now have, 

\begin{equation}
% This is the old equation, using the old way of measuring $\phi$
% m\ddot{x} = -\sigma b\psi \left[y +r\left(1-\cos\phi\right)\right]
% This is the new equation, with $\phi$ changing sign and zero-angle position
m\ddot{x} = -\sigma b\psi \left( y - r\sin\phi \right)
\end{equation}

For the $y$ equation we assume that the $y$ component of force is proportional to the $x$ component, $F_y = kF_x$.  This is a bold assumption that is justified by the resulting model's ability to predict observed patterns of erosion.  It is certainly something that needs to be examined further, but not here.  Thus,

\begin{equation}
% This is the old equation, using the old way of measuring $\phi$
% m\ddot{y} = -k\sigma b\psi \left[y +r\left(1-\cos\phi\right)\right]
% This is the new equation, with $\phi$ changing sign and zero-angle position
m\ddot{y} = -k\sigma b\psi \left( y - r\sin\phi \right)
\end{equation}

%in this paragraph, too, something changed when we changed the way we measure $\phi$.  Specifially, the sign of $\tau$.
Finally, we tackle the $\phi$ equation (\ref{eqn:phi_eqn}).  If $r$ is the distance from the CM to the cutting tip of the erodent, then the torque applied to the erodent by the surface material is approximately $\tau\approx -r F_x$ (negative because a positive $F_x$ would induce a negative rotation in our system).  Implicit in this expression for the torque are some assumptions and simplifications that we will explain below.

\begin{itemize}
\item Even though the torque is the result of a contact force distributed over the entire contact surface, we will approximate the area of contact with a single point.
\item The geometric center of the contact area would be the ideal point of application (since the contact force is uniform over the contact area), but for convenience we choose $r$ - the distant edge of the contact region (see Figure \ref{fig:initial_conditions}).  Assuming that the contact area is small compared to the size of the erodent then this should not introduce significant error into our model.
\item The $y$ component of the force does not contribute to the torque.  This is equivalent to assuming that the cutting tip starts directly below the CM, $\phi \left( 0 \right) = \pi/2$, and remains there. Thus, the $y$ component of the force passes through the CM and contributes no torque. This assumption is based on the reasonable expectation that the cutting tip will be ahead of the CM in some impacts and behind it in others, but {\em on average} the cutting tip will be directly below the CM.
\end{itemize}

We also need the mass moment of inertia $I$ of the erodent. As a representative of a typical mass moment of inertial we choose that of a flat cylinder of radius $r$ and height $b$.  Note that we are not assuming that this is a representative shape (no tip!), but that this mass moment of inertia is representative.  For such a shape the mass moment of inertia is $I = mr^2/2$.

All of these assumptions and approximations result in the $\phi$ equation:

\begin{equation}
% This is the old equation, using the old way of measuring $\phi$
% \frac{mr^2}{2}\ddot{\phi} = -r\sigma b\psi y
% This is the new equation, with $\phi$ changing sign and zero-angle position
 \frac{mr^2}{2}\ddot{\phi} = r\sigma b\psi y
\end{equation}

To complete the model we need initial conditions for each equation.  The initial conditions for the translation ($x$ and $y$) are the initial position of the cutting tip (the origin),  and the initial translational velocity components (the components of $\bar{v}_0$).  The intial rotation is the average value $\phi \left( 0 \right) = \pi/2$, and we will leave the initial rotational velocity as an undetermined constant $\dot{\phi}_0$ for now.  In summary:

\begin{eqnarray}
% These are the old equations, using the old way of measuring $\phi$
% m\ddot{x} & = & -\sigma b\psi \left[y +r\left(1-\cos\phi\right)\right]\\
% x\left( 0 \right) = 0 &,& \dot{x}\left( 0 \right) = v_0\cos\alpha\\
% m\ddot{y} & = & -k\sigma b\psi \left[y +r\left(1-\cos\phi\right)\right]\\
% y\left( 0 \right) = 0 &,& \dot{y}\left( 0 \right) = v_0\sin\alpha\\
% These are the new equations, with $\phi$ changing sign and zero-angle position
m\ddot{x} & = & -\sigma b\psi \left(y - r\sin\phi \right)\\
x\left( 0 \right) = 0 &,& \dot{x}\left( 0 \right) = v_0\cos\alpha\\
m\ddot{y} & = & -k\sigma b\psi \left(y - r\sin\phi\right)\\
y\left( 0 \right) = 0 &,& \dot{y}\left( 0 \right) = v_0\sin\alpha\\
% the $\phi$ equation did not change
\frac{mr^2}{2}\ddot{\phi} & = & -r\sigma b\psi y\\
\phi \left( 0 \right) = \pi/2 &,& \dot{\phi}\left( 0 \right) = \dot{\phi}_0
\end{eqnarray}

\subsubsection*{Linearized Equations}
It is not far-fetched to fisrt examine what happens if the erodent rotates very little during the cutting process.  Since $\phi \left( 0 \right) = \pi/2$ this would amount to considering the model when $\phi$ remains small throughout, that is $\sin\phi \approx 0$.  The relevant system of equations becomes:

\begin{eqnarray}
m\ddot{x} & = & -\sigma b\psi y\label{eqn:linear x-eqn}\\
x\left( 0 \right) = 0 &,& \dot{x}\left( 0 \right) = v_0\cos\alpha\\
m\ddot{y} & = & -k\sigma b\psi y\label{eqn:linear y-eqn}\\
y\left( 0 \right) = 0 &,& \dot{y}\left( 0 \right) = v_0\sin\alpha
\end{eqnarray}

\noindent Note that this assumption had two consequences; first, it uncoupled $\phi$ from the $x$ and $y$ equations as well as made the $\phi$ equation irrelevant (why do we need a differential equation to model what is essentially a constant?).  Second, it renders the equations for $x$ and $y$ linear, making their solution considerably simpler.  It is this linearized system that will be used to answer the questions in this project.

\subsection*{PROJECT}
Our ultimate goal is to determine the relationship between the angle of attack $\alpha$ and the volume of surface material removed and how this relationship is affected by the paramters in our model.  Therefore, we will have quantitative information that we can use to deploy mitigating material if our goal is to minimize the effects of erosion.  Or, in the case that we wish to use erosion to accomplish some desired result, such as removing rust from a ferrous surface material, then we can know how best to direct a stream of erodents to maximize material removal.

To begin, let's establish some values for the constants characterizing the erodent and surface material.  Table \ref{table:parameters} contains numbers that are typical for a quartz sand erodent and an aluminum surface material - with the exception that $\sigma$ is almost four times larger than a typical flow stress for aluminum.  Finnie \cite{Finnie1958} attempts an explanation of this, proclaiming that strain rate effects (which raise $\sigma$ for high speed events), and other factors not considered in this simplified model account for the enlarged flow stress.  Note that the units of mass are quite unusual.  {\em Dozens of Slugs} is the natural unit of mass in the English system with inches as the length unit and seconds as the time unit.

\begin{table}[ht]
\centering % used for centering table
\begin{tabular}{c c c} % centered columns (4 columns)
\hline\hline %inserts double horizontal lines
Parameter & Value & Units\\ [0.5ex] % inserts table
%heading
\hline % inserts single horizontal line
$k$ & 2 & dimensionless\\ % inserting body of the table
$\psi$ & 0.04 & dimensionless\\
$r$ & 0.04 & in\\
$b$ & 0.04 & in\\
$m$ & $2.855\times 10^{-7}$ & dozens of slugs\\
$v_0$ & 4400 & in/s\\ 
$\sigma$ & $200\times 10^3$ & lb/in$^2$\\
$\alpha$ & $20\times\frac{\pi}{180}$ & radians\\[1ex] % [1ex] adds vertical space
\hline %inserts single line
\end{tabular}
\caption{Model Parameters} % title of Table
\label{table:parameters} % is used to refer this table in the text
\end{table}

\subsubsection*{Question 1:  For the parameters provided, what are the solutions $x\left( t\right)$ and $y\left( t\right)$?}

Notice that even though the system is coupled, the $y$ equation can be solved independently by traditional methods, and its solution substituted into the $x$ equation for its subsequent solution.\newline

\noindent {\bf Hint: } Equation (\ref{eqn:linear y-eqn}) is uncoupled from equation (\ref{eqn:linear x-eqn}), so it can be solved by itself.  Once the solution $y\left( t\right)$ is known equation (\ref{eqn:linear x-eqn}) can be solved by integrating twice and using the intial conditions to determine the arbitrary constants. 

As a sanity check, Figure \ref{fig:q1a} shows the time histories of $x$, $y$, $\dot{x}$, and $\dot{y}$;  Figure \ref{fig:q1b} shows the trajectory of the particle tip parameterized by $x\left( t\right)$ and $y\left( t\right)$.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{q1a.jpg}
\caption{The $x$ and $y$ components of position and velocity for the particle tip.}
\label{fig:q1a}
\end{figure}

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{q1b.jpg}
\caption{Trajectory of the particle tip.}
\label{fig:q1b}
\end{figure}

\teacher{
\noindent {\bf Answer 1:} First, the equation for $y$ can be solved directly since it is uncoupled from the $x$ equation:
\begin{equation}
y\left( t\right)=\frac{v_0\sin\alpha}{\beta}\sin\left(\beta t\right)\textrm{, where }\beta^2 = \frac{\sigma\psi b k}{m}
\end{equation}

\noindent Once we have solved for $y$ then $x$ can be determined by direct integration:

\begin{equation}
x\left( t\right) = \frac{v_0\sin\alpha}{k\beta}\sin\left(\beta t\right)+\left(v_0\cos\alpha-\frac{v_0}{k}\sin\alpha\right) t
\end{equation}
}

\subsubsection*{Question 2:  For the parameters provided, when does the cutting stop?}

When we make our calculation of the volume of material removed we will need to calculate the area between the tip's position and the positive $x$-axis.  The solution to the differential equations (\ref{eqn:linear x-eqn}, \ref{eqn:linear y-eqn})  provide a parametric representation of the tip's trajectory.  So the missing information needed to calculate the volume is when to stop our simulation. That is, what is the time $t_c$ when the tip stops cutting?  There are two conditions that will end the cutting.

\begin{enumerate}
\item The erodent rebounds off of the surface, that is $y\left( t_c\right)=0$ for $t_c>0$.
\item The erodent loses all of its forward momentum before it can rebound, that is $y<0$, but $\dot{x}\left( t_c\right)=0$ for $t_c>0$.
\end{enumerate}

\noindent Whichever condition occurs first will determine $t_c$ and end material removal.\newline

\noindent {\bf Hint: } Figure \ref{fig:q2} shows the candidate $t_c$ values as the graphs of $0 = y\left( t_c\right)$ and $0 = \dot{x}\left( t_c\right)$.  The actual value of $t_c$ will be the lesser of the two candidates, indicated by a bold curve in the figure.  Even though these are smooth functions there is a cusp in the graph of equation $0 = \dot{x}\left( t_c\right)$.  This is because the figure plots only the real part of the complex solution to the equation - only when the solution is purely real (no imaginary part) is it of physical significance.  It transitions from a complex to a real solution precisely where the two curves intersect.\newline

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{q2.jpg}
\caption{The value of $t_c$ (bold):  Where $t_c$ is constant cutting ends because the erodent has rebounded off of the surface, the other values of $t_c$ occur when the erodent has lost its forward momentum.}
\label{fig:q2}
\end{figure}

\teacher{
\noindent {\bf Answer 2:}
The equation for $y$ was determined in answering Question 1, and in the process of solving for $x$, we would have first found $\dot{x}$.  Thus, the two possible conditions ending the cutting are:

\begin{eqnarray}
0 & = & \frac{v_0\sin\alpha}{\beta}\sin\left(\beta t_c\right)\label{eqn:tc1}\\
0 & = & \frac{v_0\sin\alpha}{k}\cos\left(\beta t_c\right) +v_0\cos\alpha-\frac{v_0}{k}\sin\alpha\label{eqn:tc2}
\end{eqnarray}

\noindent Solving these equations for $t_c$: 

\begin{equation}
t_c = \min_{t>0}\{\frac{\pi}{\beta}, \frac{1}{\beta}\arccos\left(1-\frac{k}{\tan\alpha}\right)\}
\end{equation}

\noindent A short computation will show that:

\begin{eqnarray}
k > 2\tan\alpha & \implies & t_c = \frac{\pi}{\beta}\\
k < 2\tan\alpha & \implies & t_c = \frac{1}{\beta}\left( 1-\frac{k}{\tan\alpha}\right)
\end{eqnarray}

Figure \ref{fig:q2} shows the candidate $t_c$ values as the graphs of equations (\ref{eqn:tc1}) and (\ref{eqn:tc2}); the actual value of $t_c$ that will end the cutting is bold in the figure.  Note the cusp in the graph of equation (\ref{eqn:tc2}).  This is because the figure plots only the real part of the complex solution.  At the intersection of the two curves, where the cusp occurs and $t_c$ transitions from a constant value, the imaginary part of the solution to equation (\ref{eqn:tc2}) vanishes and the solution becomes real-valued.
}

\subsubsection*{Question 3:  For the parameters provided, what is the volume removed?}

If we know the trajectory of the tip $\left( x_{tip}\left( t\right), y_{tip}\left( t\right)\right)$ then we could integrate to find the area between the trajectory of the tip and the $x$-axis.  This area multiplied by the width of the erodent, $b$, is the volume of removed material.  Since the linearized model ignores rotation of the erodent, we can use $x\left( t\right)$ and $y\left( t\right)$, starting at the origin, as surrogates for $x_{tip}\left( t\right)$ and  $y_{tip}\left( t\right)$

\begin{eqnarray}
vol & = & b\int\limits_{t=0}^{t_c} ydx\\
& = & b\int\limits_{0}^{t_c} y\dot{x}dt
\end{eqnarray}\newline

\noindent {\bf Hint: } I recommend using the expressions for $y$ and $\dot{x}$ found in Question 1, and the value of $t_c$ found in Question 2 to perform a numerical integration using something like midpoint rule. \newline

\teacher{
\noindent {\bf Answer 3:}
\noindent This integral is a little tedious, but can be done by hand.  The solution is

\begin{eqnarray}
vol & = & b\int\limits_{0}^{t_c} y\dot{x}dt\\
& = & b\int\limits_{0}^{t_c}\frac{v_0\sin\alpha}{\beta}\sin\left(\beta t\right)\left[ \frac{v_0\sin\alpha}{k}\cos\left(\beta t\right) +v_0\cos\alpha-\frac{v_0}{k}\sin\alpha\right] dx\\
& = & \frac{b}{k}\left( \frac{v_0}{2\beta}\right)^2\left[ 4\cos\left( \beta t_c\right)-\cos\left( 2\beta t_c\right)-3 \right]\sin^2\alpha + \frac{b}{2}\left( \frac{v_0}{\beta}\right)^2\left[ 1-\cos\left( \beta t_c\right)\right]\sin\left( 2\alpha\right)
\label{eqn:vol}
\end{eqnarray}

\noindent However, I recommend performing a numerical approximation of the volume integral.  The complicated expression doesn't yield any physical insight (that I can see), and its length makes the probability of incorrectly entering the formula high relative to the ease of programming a numerical integration that will ultimately provide the same answer.

The value for the volume removed corresponding to the parameters in Table \ref{table:parameters} is

\begin{equation}
vol = 1.8164\times 10^{-4}\textrm{ in}^3
\end{equation}
}
\subsubsection*{Question 4:  For the parameters provided (except $\alpha$), at what angle is the volume removed maximized?}

\noindent {\bf Hint: } The angle at which maximum material removal occurs, $\alpha_{max}$, will happen when

\begin{equation}
0 = \frac{d}{d\alpha}\left( vol\right)\bigg\rvert_{\alpha = \alpha_{max}}
\label{eqn:alpha_max}
\end{equation}

\noindent I highly recommend finding $\alpha_{max}$ numerically.  That is, repeat the calculation for the volume removed at {\it many} values of $\alpha$ between $0$ and $\pi/2$ (note that $\alpha>0$ is required).\newline

\teacher{
\noindent {\bf Answer 4:}
Given that a formula for the volume removed (equation \ref{eqn:vol}) exists, it might be tempting to directly differentiate this equation with respect to $\alpha$ and proceed to solve equation \ref{eqn:alpha_max}.  This will probably fail because $t_c$ depends on $\alpha$ in a non-trivial way that is not explicit in the formula.  Furthermore, that dependence is not differentiable at the transition between the different formulas for $t_c$.

If you neglect my recommendation and attack the problem analytically with the formula above, then you might arrive at the point where you need to solve for $\alpha$ in the equation.

\begin{equation}
\tan\left( 2\alpha\right) = \frac{4k\left( \cos\left( \beta t_c\right) -1\right)}{4\cos\left( \beta t_c\right)-\cos\left( 2\beta t_c\right)-3}
\end{equation}

\noindent This is only correct where $t_c$ is independent of $\alpha$.  That is why I recommend plotting the volume removed as a function of $\alpha$ and pick out the point on the graph where the maximum occurs.  However, if the relationship between $\alpha$ and $t_c$ is ignored (in effect assuming that $t_c = constant$) then given the correct value for $t_c$ the right answer for $\alpha_{max}$ can be found - but for the wrong reason!  In the next paragraph I explain why, but the justification for using $t_c = constant$ is far from rigorous.

From physical intuition, it could be argued that maximum material removal could never occur if the cutting stops due to $\dot{x}=0$.  On mathematical grounds this is a more difficult statement to defend. I have not seen a proof of this, but would be interested in seeing one. Despite the lack of a rigorous demonstration, it seems to be true that the cutting will stop at the angle of maximum volume removal when $y=0$. If we proceed assuming this to be true then $\beta t_c = \pi$, that is, $t_c$ is indeed indpendent of $\alpha$ where equation (\ref{eqn:alpha_max}) is satisfied.  Therefore, our expressions for the derivative of volume removed is correct in a neighborhood of $\alpha_{max}$, and is greatly simplified.  So much so that the angle at which maximum material removal occurs is

\begin{equation}
\alpha_{max} = \frac{1}{2}\arctan\left( k\right)
\end{equation}

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{q4.jpg}
\caption{The volume of material removed as a function of $\alpha$.}
\label{fig:q4}
\end{figure}

However you choose to evaluate $\alpha_{max}$, for the parameters in Table \ref{table:parameters} (except $\alpha$):

\begin{equation}
\alpha_{max} = 31.7\textrm{ degrees}
\end{equation}

}
\subsubsection*{Question 5:  Does the erodent's speed $v_0$ affect the angle of maximum material removal?}

We implicitly assume that $v_0$ is large enough to remove material - that is, the stress produced from the impact is large enough that the surface material will {\it flow}.  If we maintain that this assumption is true for any $v_0$, will different values of $v_0$ change the angle of maximum material removal?  Put differently, does the initial velocity $v_0$ have an influence on the volume of material removed during an impact?\newline

\noindent {\bf Hint: } Select (at least) 3 different initial velocities that cover a relatively large, but realistic range, say 2000 to 8000 in/s.  On the same graph plot the volume removed vs. angle of attack for each initial velocity.  From these graphs is should be clear if the speed is effecting angle of maximum material removal.\newline

\teacher{
\noindent {\bf Answer 5:}
No.  It is clear from the equation for the volume of material removed that the volume of material will change, but the angle at which this maximum occurs will not.  Figure \ref{fig:q5} shows the volume removed for three values of $v_0$.  The peak in each occurs at the same angle.


\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{q5.jpg}
\caption{The volume removed at three different impact velocities - the angle of maximum material removal does not change with the impact velocity.}
\label{fig:q5}
\end{figure}
}
\subsubsection*{Question 6:  Does the choice of erodent affect the angle of maximum material removal?}

The erodent is characterized by two parameters, its size $r$ and its mass $m$.  Within any real stream of particles we can expect a considerable amount of variability in these two numbers, they could even change mid-stream if particles collide and break-up into smaller particles.  So measuring the angle of maximum removal from a stream of particles instead of a single particle impact could potentially mislead our conclusions if $r$ and $m$ have an effect of the angle of maximum material removal.  Will either of them have a major impact on our results?\newline

\noindent {\bf Hint: } First, consider the impact of $r$.  The equations that govern the material removal are equations (\ref{eqn:linear x-eqn}, \ref{eqn:linear y-eqn}).  Is $r$ even present in these equations?  If not, then it cannot affect the results.  Second, consider the mass $m$.  As in Question 5, pick (at least) 3 masses that would be typical for small particles such as sand (use the mass value in Table \ref{table:parameters} as a guide for a typical mass).  On a single graph plot the volume of material removed vs. the angle of attack and compare the angle at which maximum material removal occurs for each mass.\newline

\teacher{
\noindent {\bf Answer 6:}
Looking at the linearized equations, $r$ is not present.  Therefore, $r$ has no effect.  This point is reinforced by evidence documented in Hutchings \cite{Hutchings1976}; experiments were carried out using larger erodents (easier to conduct single impact experiments), but the results on angle of maximum volume removed were independent of the particle size.  An argument could be made that $r$ and $m$ are not independent, thus $\alpha_{max}$ independent of $r$ implies it is also independent of $m$, but this misses the point of the question since $r$ is not present in our model.  However, our model does support the conclusion that changing the mass has no effect on the angle of maximum material removal, Figure \ref{fig:q6}.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{q6.jpg}
\caption{The volume removed using three different erodent masses - the angle of maximum material removal is independent of erodent mass.}
\label{fig:q6}
\end{figure}
}

\subsubsection*{Question 7:  Does the surface material affect the angle of maximum material removal?}

The surface material is characterized only by its flow stress $\sigma$.  The flow stress for other ductile materials that fit our canonical situation are shown in Table \ref{table:flow_stress}. 

\begin{table}[ht]
\centering % used for centering table
\begin{tabular}{c c c} % centered columns (4 columns)
\hline\hline %inserts double horizontal lines
Metal & Flow Stress & Units\\ [0.5ex] % inserts table
%heading
\hline % inserts single horizontal line
Mild Steel & $37\times 10^{3}$ & lb/in$^2$\\ % inserting body of the table
Brass & $27\times 10^{3}$ & lb/in$^2$\\
Copper & $10\times 10^{3}$ & lb/in$^2$\\
Stainless Steel & $75\times 10^{3}$ & lb/in$^2$\\ 
Titanium & $136\times 10^3$ & lb/in$^2$\\[1ex] % [1ex] adds vertical space
\hline %inserts single line
\end{tabular}
\caption{Flow Stress} % title of Table
\label{table:flow_stress} % is used to refer this table in the text
\end{table}

\noindent {\bf Hint: } On a single graph plot the volume of material removed vs. the angle of attack for 3 (or more) of the materials in Table \ref{table:flow_stress} to determine if the surface material has an effect. \newline

\teacher{
\noindent {\bf Answer 7:}
The volume removed by our projectile on all of these metals is shown in Figure \ref{fig:q7}; they all share the same angle of maximum material removal.  So it seems, that our model predicts that regardless of the velocity of the impacts, or the materials involved, we expect that the maximum material will be eroded when the erodents strike the surface at an angle of 31.7 degrees.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{q7.jpg}
\caption{The volume of material removed for the surface materials of Table \ref{table:flow_stress}}
\label{fig:q7}
\end{figure}
}

\subsubsection*{Question 8:  How do the modeling parameters $\psi$ and $k$ affect the angle of maximum material removal?}
Recall that $\psi$ is the constant of proportinality between the depth of the eroding tip, and the height of the material building up in front of the erodent.  It's value is required to calculate the area of contact between the erodent and the surface material, it's value was {\em assigned} to be 0.04.  There was no physical justification given for this value - it was more wild guess than educated guess.  

Similar to $\psi$, $k$ is a constant of proportionality.  This time, between the $x$ and $y$ components of force acting on the erodent.  That this ratio between the force compenents remains constant is a bold assumption; that the force acting on the erodent does not change direction during the cutting process.  And like $\psi$, it's value was guessed ($k = 2$).

So, were the guesses made for the values of $\psi$ and $k$ consequential in determining $\alpha_{max}$?\newline

\noindent {\bf Hint: }As in Questions 5, 6, and 7, we need to select some representative values for $\psi$ and $k$ then plot the volume of material removed vs. angle of attack.  Selecting typical values of these parameters, however, is not as straightforward as with the physical quantities of velocity, mass, radius, and flow stress.  Therefore, I recommend the following:  For $\psi$ use positive values and stay less than 5.0 (this is more than 100 times our initial choice for $\psi$).  For $k$, remember it is the ratio of the $y$-component of force to the $x$-component.  So $k$ should also be positive, and somewhere ``close" to unity, $0.1 < k < 10$. \newline

\teacher{
\noindent {\bf Answer 8 for $\psi$:}
As can be seen in Figure \ref{fig:q8a} the parameter $\psi$ has no effect on the agle of material removal.  However, it has a great effect on the amount of material removed.  Nonintuitively, increasing $\psi$ increases the material removed.

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{q8a.jpg}
\caption{The volume of material removed vs. angle of attach for 3 values of $\psi$: $\psi = 0.04$, $0.4$, and $4.0$. The angle of maximum material removal is the same for each.}
\label{fig:q8a}
\end{figure}
}

\teacher{
\noindent {\bf Answer 8 for $k$:}
Figure \ref{fig:q8b} shows the effect of altering $k$.  Finally, the angle of maximum material removal has changed.  As $k$ increases, so does the angle at which maximum material removal occurs.  This makes intuitive sense because smaller $k$ values result in smaller $y$-components of the force acting on the erodent.  This will allow the erodent to cut further before rebounding off of the surface, ending the cutting process.  

\begin{figure}[h]
\centering
\includegraphics[width=0.9\textwidth]{q8b.jpg}
\caption{The volume of material removed vs. angle of attach for 3 values of $k$: $k = 1.0$, $2.0$, and $4.0$. The angle of maximum material removal is increases with increasing $k$.}
\label{fig:q8b}
\end{figure}
}

\subsubsection*{Additional questions That are considerably more difficult:}

There are three ways that this investigation could go further, each addresses our biggest assumptions in deriving the linearized model.

\begin{enumerate}
\item What if we do not make the assumption that $\phi$ remains small throughout the cutting process?  Will this change our answer for the angle of maximum material removal?  Note that this requires making an additional assumption regarding the initial angular velocity of the erodent $\dot{\phi}\left( 0\right)$.
\item We assumed that that $F_y = kF_x$ with $k$ constant, independent of $\alpha$.  A first step towards a more more realistic model would be that the force acts in the direction opposite of the erodent's intial velocity, $k = \tan\alpha$.  This will allow us to keep the linearized system of equations while reducing the severity of one of our assumptions.  Will this change our answer of the angle of maximum material removal?
\item Addressing the $k$ assumption again, what if $k$ changes with the trajectory of the erodent?  A better expression for $k$ would produce a contact force opposing the velocity of the projectile throughout the motion, not just its initial direction.  So, what if
\begin{equation}
\frac{F_y}{F_x} = \frac{\dot{y}}{\dot{x}}\textrm{ ?}
\end{equation}
\noindent Under this new expression for $k$ even the two equation system neglecting $\phi$ is nonlinear.  Will this change our answer for the angle of maximum material removal?
\end{enumerate}

\begin{thebibliography}{98}
\bibitem{Finnie1958} Finnie, I. 1958. {\it The Mechanism of Erosion of Ductile Metals}. Proc. 3rd U.S. Natl. Congr. Appl. Mech. pp.~ 527-532.
\bibitem{Hutchings1976} Hutchings, I.M., R.E.~Winter., and J.E.~Field. 1976. {\it Solid particle erosion of metals: the removal of surface material by spherical projectiles}. Proc. R. Soc. Lond. A. 348. pp. 379-392.
\end{thebibliography}

\teacher{
\section*{NOTES FOR TEACHERS}
We found that the angle at which maximum material is removed does not depend upon the erodent, or the surface being eroded.  This couldn't be literally true, but it's not difficult to imagine that this critical angle is not so different for a great deal of materials that we would be interested in.  We also found that the angle is independent of $\psi$, but not $k$.  In practice, these two parameters are used to {\em tune} the model to data.  In a typical experiment, some known mass of erodent is applied to a plate, also of known mass.  After the experiment the mass of the plate is measured again (after cleaning to remove any embedded erodents) to determine how much material was actually removed.  The erodents are also measured to determine a typical size $r$ and mass $m$.  Simple division of the total erodent mass by the typical mass determines how many particles are in the mass of erodent.  If the experiment is repeated at many different angles then the angle of maximum material removal can be determined.  The value of $k$ is then adjusted in the model so that the model predicts the correct $\alpha_{max}$.  Then $\psi$ is determined so that the correct amount of material is removed in an ideal single impact.  In fact, one alternate interpretation of the parameter $\psi$ is that it quantifies the fraction of particles that strike the surface in the ideal manner that we have modeled.  In reality, many erodents will not cut as we have modeled because many impacts will occur with a flat side of the erodent striking the surface (a belly flop), and some particles in the stream never strike the surface, having been redirected by impacts with other erodents.  In fact, it is expected that very few particles will actually erode the surface material in the manner described by our model.

However the modeling parameters are interpreted, there does seem to be a ``correct" $\psi$ and $k$ for each pair of erodent and surface material.  It is analogous to the coefficients of friction between different materials - every different pair of materials possess coefficients of static and kinetic friction that can only be determined by experiment.

One further comment: this is just one method by which material is removed from the surface.  Experiments consistently show that while this model does well at predicting material removal trends at low angles, from grazing to past $\alpha_{max}$, at angles closer to normal impingement the model under predicts material removed, \cite{Finnie1958}.  This is because as the erodents strike the surface at close to normal impact they don't ``cut" material.  The mechanism is more of a fracture where the surface material breaks brittlely, then is liberated by subsequent impacts.  Finnie used a {\em correction factor} at angles greater than $\alpha_{max}$ to account for this additional material removed at large angles so that his model could characterize erosion for the full range of impact angles.
}

\end{document}
