If a body has an insignificant rotation compared to the path its center of gravity travels, it is said to be a particle. In later chapters we will make an proper introduction to generalized Newtonian-Euler motion equations, but for now, suffice to say if either the rotational acceleration or the size of the body is small enough, then we can neglect the change of angular momentum
The equation of motion for a particle is given by Newton’s second law:
\[ \sum F = m a \]
and the rotational equation of motion is given by Euler’s second law:
\[ \sum M_G = \bar{I}\alpha \]
If either the angular acceleration, \(\alpha\) or the size of the body (and thus \(\bar I\)) is small enough, then we can neglect the change of angular momentum, and the rotational equation of motion reduces to:
\[ \alpha = 0 \text{ or } \bar{I} = 0 \Rightarrow \sum M_G = \bar{I}\alpha = 0 \tag{6.2.1}\]
This is what defines a particle.
We start with kinematics and the mathematical descriptions of a bodies position, velocity and acceleration. This is the foundation of dynamics, including kinetics.
Assume an universal coordinate system and stick to it! We can choose between cartesian, normal and polar coordinate systems.
The position of a bodies center of gravity as a function of time \(t\) can be expressed using the position vector \(\mathbb{r}(t) = [x(t), y(t), z(t)]^T\) where each component is a function of the same independent variable \(t\). The components are called coordinate functions and \(\mathbb{r}(t)\) is known as the parameter curve with the parameter \(t\). The parameter curve is also known as the curve or path.
The mean velocity between two positions \(\mathbb{r}(t)\) and \(\mathbb{r}(t+\Delta t)\) can be described as \(\bar{\mathbb{v}} = \frac{\Delta \mathbb{r}}{\Delta t}\). In the limit where \(\Delta t \to 0\) we get the time derivate of the position vector, known as the velocity vector
\[ \lim_{\Delta t \to 0} \frac{\Delta \mathbb{r}}{\Delta t} = \frac{d\mathbb{r}}{dt} =: \dot{\mathbb{r}} = \mathbb{v} \tag{6.2.2}\]
As can be seen in the animation, the velocity vector \(\dot{\mathbb{r}}\) at \(t\) is always tangential to the path.
The dot notation in \(\dot{r}\) or in \(\dot{\ddot{\ddot{\ddot{\ddot{r}}}}} = \frac{d^9}{dt^9}r\) (a.k.a. Newton’s abomination) comes from Newton who studied dynamic problems and needed to develop the mathematics to deal with these problems.
The mathematics came to be known as calculus which was expanded and is really made up by a great many contributors (Leibnitz, Lagrange). see e.g., this video for an historical overview.
The magnitude of the velocity vector, \(|\dot{\mathbb{r}}|\), is known as the speed. The difference between velocity and speed being that velocity also contains information about the direction of the object.
Similarly, the acceleration is given by the change of velocity over time, \(\frac{d\dot{\mathbb{r}}}{dt}\). The average acceleration at time \(t\) can be formulated as
\[ \lim_{\Delta t \to 0} \frac{\Delta \dot{\mathbb{r}}}{\Delta t} = \frac{d\dot{\mathbb{r}}}{dt} = \frac{d^2\mathbb{r}}{dt^2} =: \ddot{\mathbb{r}} = \mathbb{a} \tag{6.2.3}\]
From the animation we can see that the acceleration is pointing inwards, towards the center of curvature. Also, we see that the acceleration in general differs from the direction of the position vector. The acceleration vector can be decomposed into a tangential and a normal component, \(\mathbb{a} = a_t \mathbb{e}_t + a_n \mathbb{e}_n\) where \(\mathbb{e}_t\) is the unit tangent vector and \(\mathbb{e}_n\) is the unit normal vector. The tangential component \(a_t\) is responsible for changes in speed, while the normal component \(a_n\) is responsible for changes in direction. If the speed is constant, then \(a_t = 0\) and the acceleration is purely normal.
In the special case of rectilinear motion (motion without curvature), such as described in Figure 6.2.3
the position vector reduces to a single component, for example \(x(t)\), and the velocity and acceleration reduces to the time derivatives of this single component. We sometimes use the (simple) scalar notation
\[ s = |\mathbb{r}|=x(t), v = |\dot{\mathbb{r}}|=\dot x(t) \text{ and } a = |\ddot{\mathbb{r}}|=\ddot x(t) \tag{6.2.4}\]
We can express these scalar values using the Leibniz notation for time derivative
\[ v = \frac{ds}{dt} \text{ and } a = \frac{dv}{dt} \tag{6.2.5}\]
solving for the time increment \(dt\), we get
\[ dt = \frac{ds}{v} = \frac{dv}{a} \tag{6.2.6}\]
from which we get
\[ \boxed{a \, ds = v \, dv} \tag{6.2.7}\]
Which is convenient to use in special cases where we are not interested in time explicitly, but should be used with care and really be avoided all together, instead the proposed working method is described in what follows.
Working purely with the scalar notation from the previous section is troublesome. A strong recommendation is to avoid the use of derived formulas often found in classical books in mechanics and formularies, since it is often more tedious to figure out the assumptions under which these formulas are valid than to create a model using the differential equations directly. Many ready-to-go formulas assume the acceleration to be constant. We shall also refrain from using the notations \(s\), \(v\) and \(a\) for varying position, velocity and acceleration since it is harder to let go of the “thinking by formularies” and move towards “Concept Based Modeling” with “Computational Thinking”. Instead use \(\mathbb{r}, \dot{\mathbb{r}}\) and \(\ddot{\mathbb{r}}\) along with knowledge of differential equations!
Remember that the derivative and integral operator works naturally on vectors i.e.,
\[ \dot{\mathbb{r}}(t) = [\dot x(t), \dot y(t), \dot z(t)]^\mathsf T \quad \text{and} \quad \int \mathbb{r}(t) dt = \left[\int x(t) dt, \int y(t) dt, \int z(t) dt\right]^\mathsf T \]
The velocity (vector) is the time derivative of the position (vector), \(\dot{\mathbb{r}} = \frac{d\mathbb{r}}{dt}\) with respect to time, \(t\). The magnitude \(|\dot{\mathbb{r}}|\) is known as speed (scalar value). The direction is tangent to the trajectory path. The acceleration (vector) is the time derivative of the velocity (vector) \(\ddot{\mathbb{r}} = \frac{d\dot{\mathbb{r}}}{dt} = \frac{d}{dt}(\frac{d\mathbb{r}}{dt}) = \frac{d^2\mathbb{r}}{dt^2}\), or the second time derivative of the position (vector).
The acceleration (vector) is the time derivative of the velocity (vector) \(\ddot{\mathbb{r}} = \frac{d\dot{\mathbb{r}}}{dt} = \frac{d}{dt}(\frac{d\mathbb{r}}{dt}) = \frac{d^2\mathbb{r}}{dt^2}\), or the second time derivative of the position (vector).
In what follows we shall rely more on working with the vector form of \(\mathbb{r}\) and formulate the differentials as ordinary differential equations which we can easy solve using either a symbolic differential equation solve, e.g., sympy.dsolve() or numerically by utilizing e.g., the Euler method.
The modern workflow in computational dynamics can be summarized in the figure below
To get between the different quantities position, \(\mathbb{r}(t)\), velocity \(\dot{\mathbb{r}}(t)\) and acceleration \(\ddot{\mathbb{r}}(t)\) we take the derivative in one direction and integrate in the other. In practice one does not typically manually integrate expressions, one typically solves an ordinary differential equation (ODE) instead, which takes us all the way down to \(\mathbb{r}(t)\) which is the solution to an ODE. Then we can take the derivative to determine the needed quantities. This is typical for a symbolic solution. Non-linear ODEs, on the other hand, are typically solved using various numerical algorithms which will first generate velocity, \(\dot{\mathbb{r}}\) on the way down to the final position \(\mathbb{r}\).
Newtonian physics is formulated as the kinetic equation
\[ \sum \mathbb{F} = m\ddot{\mathbb{r}} \]
which is a second order ODE, but in kinematics we se equally often differential equations of first order. Regardless, it is very convenient to solve any equation using a computational approach.
The computational approach to dynamics is to formulate an initial value problem (IVP) and solve it using a differential equation solver.
The inquiries to the model can be essentially divided into two categories:
There are several types of bases (or coordinate systems) on which we can formulate kinematic expressions and work with derivatives, traditionally this is done to simplify calculations, done by hand, and the resulting symbolic expressions.
We can avoid this whole exercise of practicing to compute velocities and acceleration in various bases by just formulating relations in Cartesian coordinates and directly work with the position vector \(\mathbb{r}\) which can be easily projected onto any other bases if needed. Let the computer do the tedious mathematics!
We shall here derive the most common coordinate systems and their base vectors such that this projection can be done.
Sometimes the bases can be denoted simply \(\mathbb i, \mathbb j\) and \(\mathbb k\), but throughout this treatment we will try to be explicit and denote basis vectors and unit vectors with a specific sub-index for clarity.
The tangent vector \(\mathbb{e}_t\) is always pointing in the parameter direction, tangential to the path and is defined as the direction of the velocity
\[ \boxed{\mathbb{e}_t = \frac{\dot{\mathbb{r}}}{|\dot{\mathbb{r}}|} } \]
The normal vector \(\mathbb{e}_n\) is perpendicular to \(\mathbb{e}_t\) and always points inwards to the center of curvature
\[ \boxed{\mathbb{e}_n = \frac{\dot{\mathbb{e}}_t}{|\dot{\mathbb{e}}_t|}} \]
here, suddenly a time derivate of a unit vector appears, see this section for an explanation.
Finally the lesser known bi-normal vector is defined as the cross product between the tangent and normal vectors,
\[ \boxed{ \mathbb{e}_b := \mathbb{e}_t \times \mathbb{e}_n } \]
which defines a right-hand system. Note that for planar problems (2D) \(\mathbb{e}_b = \mathbb{e}_z\).
\(\mathbb{e}_r\) has the same direction as \(\mathbb{r}(t)\) and is given by
\[ \boxed{\mathbb{e}_r = \frac{\mathbb{r}(t)}{|\mathbb{r}(t)|}} \]
\(\mathbb{e}_\theta\) is given by
\[ \boxed{\mathbb{e}_\theta = \frac{\dot{\mathbb{e}}_r}{|\dot{\mathbb{e}}_r|}} \]
In order to ensure orthogonality of two vectors, we can take their dot product, if the resulting scalar is zero, they must be orthogonal. If two vectors are parallel then their dot product must be one, i.e., \(\mathbb{e} \cdot \mathbb{e} = 1\). Now we can take the time derivative of this expression \(\frac{d}{dt}(\mathbb{e} \cdot \mathbb{e} = 1) \stackrel{\text{product rule}}{=} \frac{d}{dt}\mathbb{e} \cdot \mathbb{e} + \mathbb{e} \cdot \frac{d}{dt}\mathbb{e} = 0 = \dot{\mathbb{e}} \cdot \mathbb{e} + \mathbb{e} \cdot \dot{\mathbb{e}} = 0 \Leftrightarrow \dot{\mathbb{e}} \cdot \mathbb{e} = 0 \Rightarrow \mathbb{e} \perp \dot{\mathbb{e}}\)
Also note that typically \(|\dot{\mathbb{e}}| \neq 1\).
So at the limit \(\Delta t \to 0\) we clearly see that \(\mathbb{e} \perp \dot{\mathbb{e}}\).
Let \(\mathbb{r}(t) = [\cos{\theta(t)}, \sin{\theta(t)}]\)
Be careful when defining angles as degrees! Angles will be functions of time in what follow, meaning implicit derivation will yield expressions where angles are directly multiplied with distances or velocities, this only works if the angles are defined as radians!
We start by defining the angle as a function of time, \(\theta(t)\) and the position vector as a function of this angle.
import sympy as sp
t = sp.symbols('t', real=True, positive=True)
theta = sp.Function('theta', real=True)(t)
rr = sp.Matrix([sp.cos(theta), sp.sin(theta)])
rr\(\displaystyle \left[\begin{matrix}\cos{\left(\theta{\left(t \right)} \right)}\\\sin{\left(\theta{\left(t \right)} \right)}\end{matrix}\right]\)
We compute the normalized position vector \(\mathbb{e}_r\)
rr.normalized()\(\displaystyle \left[\begin{matrix}\frac{\cos{\left(\theta{\left(t \right)} \right)}}{\sqrt{\sin^{2}{\left(\theta{\left(t \right)} \right)} + \cos^{2}{\left(\theta{\left(t \right)} \right)}}}\\\frac{\sin{\left(\theta{\left(t \right)} \right)}}{\sqrt{\sin^{2}{\left(\theta{\left(t \right)} \right)} + \cos^{2}{\left(\theta{\left(t \right)} \right)}}}\end{matrix}\right]\)
ee_r = sp.trigsimp(rr.normalized())
ee_r\(\displaystyle \left[\begin{matrix}\cos{\left(\theta{\left(t \right)} \right)}\\\sin{\left(\theta{\left(t \right)} \right)}\end{matrix}\right]\)
We cannot simplify the trigonometric expressions if the angle \(\theta\) is not assumed to be a real valued function.
The time derivative of the position vector, \(\dot{\mathbb{e}}_r\), is given by
ee_r_dot = ee_r.diff(t)
ee_r_dot\(\displaystyle \left[\begin{matrix}- \sin{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \theta{\left(t \right)}\\\cos{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \theta{\left(t \right)}\end{matrix}\right]\)
What we computed above is the a scaled direction vector which is orthogonal to the direction vector \(\mathbb{e}_r\). Note that the time derivate of \(\theta\) appears from the implicit differentiation, which is why \(|\dot{\mathbb{e}}|\neq 1\) necessarily. In fact this is the angular velocity, \(\omega = \dot{\theta}\). The direction of rotation is the sign of \(\dot{\theta}\) i.e., {-1,1}.
The above we re-write in Newton notation and define the basis vector as \(\mathbb e_\theta\)
\[ \boxed{\dot{\mathbb{e}}_r = \dot{\theta}\mathbb{e}_\theta} \]
Thus we can get the perpendicular vector, \(\mathbb{e}_\theta\), by normalizing \(\dot{\mathbb{e}}_r\), i.e.,
\[ \boxed{\mathbb{e}_\theta = \frac{\dot{\mathbb{e}}_r}{|\dot{\mathbb{e}}_r|}} \]
ee_theta = sp.trigsimp(ee_r_dot.normalized().subs(sp.diff(theta, t), sp.symbols(r'\dot{\theta}', real=True)))
ee_theta\(\displaystyle \left[\begin{matrix}- \frac{\dot{\theta} \sin{\left(\theta{\left(t \right)} \right)}}{\left|{\dot{\theta}}\right|}\\\frac{\dot{\theta} \cos{\left(\theta{\left(t \right)} \right)}}{\left|{\dot{\theta}}\right|}\end{matrix}\right]\)
Here we can see that \(\frac{\dot{\theta}}{|\dot{\theta}|}\) only returns the sign of \(\dot{\theta}\). Thus in this example \(\mathbb{e}_\theta = [-\sin\theta, \cos\theta]^\mathsf T\) and we get
\[ \mathbb{e}_r \cdot \mathbb{e}_\theta = [\cos\theta, \sin\theta] \cdot [-\sin\theta, \cos\theta] = -\cos\theta\sin\theta + \sin\theta\cos\theta = 0 \]
Thus we have shown that
\[ \mathbb{e}_r \perp \dot{\mathbb{e}}_r, \text{ and } \mathbb{e}_r \perp \mathbb{e}_\theta \]
□
We shall here derive the kinematic equations explicitly in polar coordinates. Note again that in practice, working with a symbolic math manipulator like sympy, we can just formulate $ $ and \(\ddot{\mathbb r}\) in Cartesian coordinates and project onto the base. There is no need to learn the derived expressions below by heart in a modern workflow, instead we explore the connections between these entities and their physical meaning.
We begin by stating the earlier established bases.
\[ \mathbb e_r = \frac{\mathbb{r}}{|\mathbb{r}|} \]
\[ \mathbb e_\theta = \frac{\dot{\mathbb{e}}_r}{|\dot{\mathbb{e}}_r|} \]
We can express the position vector using the distance and base
\[ \mathbb{r} = r\mathbb{e}_r, \]
where \(r\) is the scalar length of the position vector.
Taking the derivative of the position vector we get, using the product rule
\[ \dot{\mathbb{r}} = \frac{d}{dt}(r(t)\mathbb{e}_r(t)) = \dot{r}\mathbb{e}_r + r\dot{\mathbb{e}}_r \tag{6.2.8}\]
and with \(\dot{\mathbb{e}}_r = \dot{\theta}\mathbb{e}_\theta\) from the previous section we have
\[ \boxed{\dot{\mathbb{r}} = \dot{r}\mathbb{e}_r + r\dot{\theta}\mathbb{e}_\theta} \]
where
\[ \boxed{\omega = r\dot{\theta}} \]
is recognized as the angular velocity.
Taking the time derivative of the expression for velocity Equation 6.2.8 above we get
\[ \ddot{\mathbb{r}} = \frac{d}{dt}(\dot{\mathbb{r}}) = (\ddot{r}\mathbb{e}_r + \dot{r}\dot{\mathbb{e}}_r) + (\dot{r}\dot{\theta}\mathbb{e}_\theta + r\ddot{\theta}\mathbb{e}_\theta + r\dot{\theta}\dot{\mathbb{e}}_\theta) \]
and with \(\dot{\mathbb{e}}_r = \dot{\theta}\mathbb{e}_\theta\) and \(\dot{\mathbb{e}}_\theta = -\dot{\theta}\mathbb{e}_r\) we have the expression of the acceleration vector expressed in polar coordinates
\[ \ddot{\mathbb{r}} = (\ddot{r} - r\dot{\theta}^2)\mathbb{e}_r + (r\ddot{\theta} + 2\dot{r}\dot{\theta})\mathbb{e}_\theta \]
Arguably more important are the kinematic relation expressed in Natural coordinates and even more useful are the acceleration terms which are split into normal and tangential directions. Here we shall derive these expressions explicitly, but note that we can always get these by just projecting the total velocity or acceleration vector onto any bases. Thus, there is no need to know these formulas or use them only for the sake of computing the quantities in one of these special bases.
We begin by stating the earlier established bases:
\[ \boxed{\mathbb{e_t} = \frac{\dot{\mathbb r}}{|\dot{\mathbb r}|}} \]
\[ \boxed{\mathbb{e_n} = \frac{\dot{\mathbb e_t}}{|\dot{\mathbb e_t}|}} \]
Next, we introduce the center of curvature, \(C\), at some point \(\mathbb{r}(t)\) and the corresponding curvature radius \(R\). We then can fit a circle arc length in the neighborhood of \(\mathbb{r}(t)\) such that we get a constant radius \(R\) for an infinitesimal circle arc angle \(d\beta\). This way we can get a relation between the section distance \(ds\) and the angle
This circle arc length relation, along with the Pythagorean theorem and trigonometry is a corner stone in dynamics and the main building blocks in our modelling toolbox.
We remember the definition of velocity as \(v = \frac{ds}{dt}\) and using the arc length relation we get
\[ \dot{\mathbb{r}} = \dot{r}\mathbb{e}_t = \frac{ds}{dt} \mathbb{e}_t = R\dot{\beta}\mathbb{e}_t \]
\[ \boxed{\dot{r} = R\dot{\beta}} \]
This relation is useful for when the angular rate \(\dot{\beta}\) is sought.
As seen in the figure above, a small change in the angle \(d\beta\) causes a corresponding change \(d\mathbb{e}_t\) in the unit tangent vector. The direction of \(d\mathbb{e}_t\) is perpendicular to \(\mathbb{e}_t\) and is defined as the normal direction \(\mathbb{e}_n\). We can thus formulate the relation:
\[ d\mathbb{e}_t = \mathbb{e}_n d\beta \Rightarrow \frac{d\mathbb{e}_t}{dt} = \frac{d\beta}{dt}\mathbb{e}_n \Rightarrow \dot{\mathbb{e}}_t = \dot{\beta}\mathbb{e}_n = \frac{\dot{r}}{R}\mathbb{e}_n \]
Thus we arrive at the relation
\[ \dot{\mathbb{e}}_t = \frac{\dot{r}}{R}\mathbb{e}_n \tag{6.2.9}\]
\[ R = \frac{\dot{r}}{|\dot{\mathbb{e}}_t|} \]
Since \(|\dot{\mathbb{r}}| = R\dot{\beta}\) we get
\[ \dot{\beta} = \frac{\dot{r}}{R} = \frac{\dot{r}^2}{|\dot{\mathbb{e}}_t|} \]
We see from the figure above that the acceleration vector can be split into a tangential and a normal component, \(\ddot{\mathbb{r}} = a_t\mathbb{e}_t + a_n\mathbb{e}_n\). We can find these components by taking the time derivative of the velocity vector and using the product rule
\[ \mathbb{a} = \ddot{\mathbb{r}} = \frac{d}{dt}(\dot{r}\mathbb{e}_t) = \frac{d\dot{r}}{dt}\mathbb{e}_t + \dot{r}\frac{d\mathbb{e}_t}{dt} \] \[ \ddot{\mathbb{r}} = \ddot{r}\mathbb{e}_t + \dot{r}\dot{\mathbb{e}}_t \]
and with Equation 6.2.9 we get
\[ \ddot{\mathbb{r}} = \underbrace{\frac{\dot{r}^2}{R}\mathbb{e}_n}_{a_n} + \underbrace{\ddot{r}\mathbb{e}_t}_{a_t} \]
From this expression we note that we always have an acceleration as long as there is non-zero velocity and a curved path.
This split into the normal and tangential directions will be revisited and used once we get to the kinetics chapter.
The theory of curves is an extensive topic within the field of differential geometry. Here, we will only introduce some basic concepts, such as the arc length. For a general curve, the arc length may not have a closed-form solution and is found by integrating the speed \(|\dot{\mathbb{r}}(t)|\):
\[ s(t) = \int_0^t |\dot{\mathbb{r}}(t)| dt \]
Using the arc length \(s\) as the parameter instead of time \(t\) can sometimes lead to simpler expressions for derivatives. However, since \(s\) is itself a function of \(t\), this change of variable does not necessarily simplify computer-based calculations.
Another important concept is curvature, \(\kappa(t) = \frac{1}{R(t)}\). As expected, the curvature approaches zero as the radius of curvature \(R(t)\) tends to infinity, which corresponds to a straight path. The center of curvature, \(\mathbb{r}_C(t)\), is given by:
\[ \mathbb{r}_C(t) = \mathbb{r}(t) + R(t)\mathbb{e}_n(t) \]
This point is the center of the osculating circle (a term coined by Leibniz as Circulus Osculans), which is the circle that best approximates the curve at a given point. Natural coordinates are particularly useful for this type of analysis.
In the special case when \(R(t)\) is constant the bases of polar coordinates coincide with the bases of the natural coordinates.
The torsion of a curve, \(\tau(t)\), measures how sharply it twists out of its plane of curvature. A simple analogy is the pitch of a screw. The derivations for both curvature and torsion can be complex, see e.g., the Frenet-Serret formulas, so the formulas are stated here without proof:
\[ \kappa(t) = \frac{|\dot{\mathbb{r}}(t) \times \ddot{\mathbb{r}}(t)|}{|\dot{\mathbb{r}}(t)|^3} \]
\[ \tau(t) = - \frac{(\dot{\mathbb{r}}(t) \times \ddot{\mathbb{r}}(t)) \cdot \dot{\ddot{\mathbb{r}}}(t)}{|\dot{\mathbb{r}}(t) \times \ddot{\mathbb{r}}(t)|^2} \]
where \(\dot{\ddot{\mathbb{r}}}\) denotes the third-order time derivative of the position vector, also known as jerk or jolt. For the names of higher-order derivatives, see this list.
Let us work out the curvature and torsion on a simple example, let
\[ \mathbb{r}(t) = R [\cos(t), \sin(t), 0]^\mathsf{T} \]
import sympy as sp
t = sp.symbols('t', real=True, positive=True)
R = sp.symbols('R', real=True, positive=True)
r = R * sp.Matrix([sp.cos(t), sp.sin(t), 0])
r\(\displaystyle \left[\begin{matrix}R \cos{\left(t \right)}\\R \sin{\left(t \right)}\\0\end{matrix}\right]\)
r_dot = r.diff(t)
r_ddot = r.diff(t, 2)
r_dddot = r.diff(t, 3)
# kappa = |r_dot x r_ddot| / |r_dot|^3
kappa = sp.simplify((r_dot.cross(r_ddot)).norm() / r_dot.norm()**3)
kappa\(\displaystyle \frac{1}{R}\)
# tau = - (r_dot x r_ddot) . r_dddot / |r_dot x r_ddot|^2
# The cross product r_dot.cross(r_ddot) is not zero, so we can proceed.
tau = sp.simplify(- (r_dot.cross(r_ddot)).dot(r_dddot) / (r_dot.cross(r_ddot)).norm()**2)
tau\(\displaystyle 0\)
The results seem reasonable, torsion is zero since the path is zero in the \(z\)-direction.
□
Let’s define a path and find the necessary components for the osculating circle.
The path is given by: \[ \mathbb{r}(t) = [t, \sqrt{t}\sin(2t)]^\mathsf{T} \]
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from IPython.display import Image, display
# Define the symbolic variables
t = sp.Symbol('t')
# Define the path
r = sp.Matrix([t, sp.sqrt(t)*sp.sin(2*t)])
# Calculate derivatives
r_dot = r.diff(t)
r_ddot = r.diff(t, 2)
# Calculate curvature
kappa = (r_dot[0]*r_ddot[1] - r_dot[1]*r_ddot[0]) / (r_dot[0]**2 + r_dot[1]**2)**(3/2)
R = 1/kappa
# Calculate the normal vector
e_t = r_dot / r_dot.norm()
e_n = sp.Matrix([-e_t[1], e_t[0]])
# Calculate the center of curvature
r_C = r + R*e_n
# Lambdify the symbolic expressions for numerical evaluation
r_func = sp.lambdify(t, r, 'numpy')
r_C_func = sp.lambdify(t, r_C, 'numpy')
R_func = sp.lambdify(t, R, 'numpy')
e_t_func = sp.lambdify(t, e_t, 'numpy')
# Set up the figure and axes for the animation
fig, ax = plt.subplots(figsize=(8, 6))
fig.patch.set_facecolor('#212121')
ax.set_facecolor('#212121')
ax.set_aspect('equal')
ax.set_xlabel('x', color='white')
ax.set_ylabel('y', color='white')
ax.set_title('Osculating Circle Animation', color='white')
ax.tick_params(axis='x', colors='white')
ax.tick_params(axis='y', colors='white')
for spine in ax.spines.values():
spine.set_edgecolor('white')
# Define the time range for the animation
t_vals = np.linspace(0.01, 10, 400)
path = r_func(t_vals)
# Plot the path
ax.plot(path[0].flatten(), path[1].flatten(), 'c-', label='Path')
# Initialize the plot elements to be animated
particle, = ax.plot([], [], 'ro', label='Particle')
circle_line, = ax.plot([], [], 'g-', label='Osculating Circle')
tangent_arrow = ax.quiver(0, 0, 0, 0, color='b', scale=1, angles='xy', scale_units='xy', label='Tangent')
radius_arrow = ax.quiver(0, 0, 0, 0, color='r', scale=1, angles='xy', scale_units='xy', label='Radius of Curvature')
# Set the axis limits
ax.set_xlim(0,10)
ax.set_ylim(-3,3)
#legend = ax.legend(loc='upper right')
#legend.get_frame().set_facecolor('#424242')
#for text in legend.get_texts():
# text.set_color('white')
ax.grid(False)
# Animation update function
def update(frame):
# Get the current position, center of curvature, and radius
pos = r_func(frame)
center = r_C_func(frame)
radius = R_func(frame)
et = e_t_func(frame)
# Update the particle's position
particle.set_data(pos[0], pos[1])
# Update the tangent arrow
tangent_arrow.set_offsets(pos.flatten())
tangent_arrow.set_UVC(et[0], et[1])
# Update the radius arrow
radius_vector = center - pos
norm_radius_vector = radius_vector / np.linalg.norm(radius_vector)
radius_arrow.set_offsets(pos.flatten())
radius_arrow.set_UVC(norm_radius_vector[0], norm_radius_vector[1])
# Generate the points for the osculating circle
theta = np.linspace(0, 2 * np.pi, 100)
circle_x = center[0] + radius * np.cos(theta)
circle_y = center[1] + radius * np.sin(theta)
circle_line.set_data(circle_x, circle_y)
return particle, circle_line, tangent_arrow, radius_arrow
# Create the animation
ani = FuncAnimation(fig, update, frames=t_vals, blit=False, interval=50)
# Save and display the animation as a GIF
writer = PillowWriter(fps=20)
gif_path = "graphics/kinematics_osculating_circle.gif"
ani.save(gif_path, writer=writer)
plt.close(fig)A special case of curvilinear motion is planar circular motion.
Let us study the case of a particle in polar coordinates using the arm below
We begin by expressing the particles position
\[ \mathbb{r} = r(t)[\cos\theta(t), \sin\theta(t)]^\mathsf{T} \]
import sympy as sp
t = sp.symbols('t', real=True, positive=True)
theta = sp.Function('theta', real=True)(t)
r = sp.Function('r', real=True, positive=True)(t)
rr = r * sp.Matrix([sp.cos(theta), sp.sin(theta)])
rr\(\displaystyle \left[\begin{matrix}r{\left(t \right)} \cos{\left(\theta{\left(t \right)} \right)}\\r{\left(t \right)} \sin{\left(\theta{\left(t \right)} \right)}\end{matrix}\right]\)
We then compute the derivatives
rr_dot = rr.diff(t)
rr_dot\(\displaystyle \left[\begin{matrix}- r{\left(t \right)} \sin{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \theta{\left(t \right)} + \cos{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} r{\left(t \right)}\\r{\left(t \right)} \cos{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} \theta{\left(t \right)} + \sin{\left(\theta{\left(t \right)} \right)} \frac{d}{d t} r{\left(t \right)}\end{matrix}\right]\)
We can clean up the expression by introducing the angular velocity in Newton notation, \(\dot{\theta} = \frac{d}{dt}\theta\) and the angular acceleration \(\ddot{\theta} = \frac{d^2}{dt^2}\theta\).
rr_dot = rr_dot.subs(sp.diff(theta, t), sp.symbols(r'\dot{\theta}'))
rr_dot = rr_dot.subs(sp.diff(r, t), sp.symbols(r'\dot{r}'))
rr_dot\(\displaystyle \left[\begin{matrix}- \dot{\theta} r{\left(t \right)} \sin{\left(\theta{\left(t \right)} \right)} + \dot{r} \cos{\left(\theta{\left(t \right)} \right)}\\\dot{\theta} r{\left(t \right)} \cos{\left(\theta{\left(t \right)} \right)} + \dot{r} \sin{\left(\theta{\left(t \right)} \right)}\end{matrix}\right]\)
Now the acceleration can be computed as
rr_ddot = rr.diff(t, 2).subs(sp.diff(theta, t, 2), sp.symbols(r'\ddot{\theta}'))
rr_ddot = rr_ddot.subs(sp.diff(theta, t), sp.symbols(r'\dot{\theta}'))
rr_ddot = rr_ddot.subs(sp.diff(r, t,2), sp.symbols(r'\ddot{r}'))
rr_ddot = rr_ddot.subs(sp.diff(r, t,1), sp.symbols(r'\dot{r}'))
rr_ddot\(\displaystyle \left[\begin{matrix}- \ddot{\theta} r{\left(t \right)} \sin{\left(\theta{\left(t \right)} \right)} + \ddot{r} \cos{\left(\theta{\left(t \right)} \right)} - \dot{\theta}^{2} r{\left(t \right)} \cos{\left(\theta{\left(t \right)} \right)} - 2 \dot{\theta} \dot{r} \sin{\left(\theta{\left(t \right)} \right)}\\\ddot{\theta} r{\left(t \right)} \cos{\left(\theta{\left(t \right)} \right)} + \ddot{r} \sin{\left(\theta{\left(t \right)} \right)} - \dot{\theta}^{2} r{\left(t \right)} \sin{\left(\theta{\left(t \right)} \right)} + 2 \dot{\theta} \dot{r} \cos{\left(\theta{\left(t \right)} \right)}\end{matrix}\right]\)
In the expression above, we made substitutions to make the expression more compact. Clearly we can handle tedious expressions like this using a computer algebra system.
Now we define the polar bases
ee_r = sp.Matrix([sp.cos(theta), sp.sin(theta)])
ee_theta = sp.Matrix([-sp.sin(theta), sp.cos(theta)])where \(\mathbb{e}_\theta\) is perpendicular to \(\mathbb{e}_r\) and rotated 90 degrees counter-clockwise.
Now just project the Cartesian vectors onto the polar bases. We can just use the dot product to project onto the bases since the bases are orthonormal, i.e., \(|\mathbb e_r| = |\mathbb e_\theta| = 1\).
We get \(\mathbb{r}_{r\theta} = [\mathbb{e}_r \cdot \mathbb{r} , \mathbb{e}_\theta \cdot \mathbb{r} ]\).
rr_dot_polar = sp.Matrix([ee_r.dot(rr_dot), ee_theta.dot(rr_dot)])
rr_dot_polar = sp.simplify(rr_dot_polar)
rr_dot_polar\(\displaystyle \left[\begin{matrix}\dot{r}\\\dot{\theta} r{\left(t \right)}\end{matrix}\right]\)
rr_ddot_polar = sp.Matrix([ee_r.dot(rr_ddot), ee_theta.dot(rr_ddot)])
rr_ddot_polar = sp.simplify(rr_ddot_polar)
rr_ddot_polar\(\displaystyle \left[\begin{matrix}\ddot{r} - \dot{\theta}^{2} r{\left(t \right)}\\\ddot{\theta} r{\left(t \right)} + 2 \dot{\theta} \dot{r}\end{matrix}\right]\)
We can see that we end up with the same expressions as we have tediously derived using the classical approach in 6.2.10.
For the special case of constant radius, \(r(t) = R\), we get the well-known expressions for circular motion
rr_dot_polar.subs(sp.symbols(r'\dot{r}'), 0).subs(sp.symbols(r'\ddot{r}'), 0)\(\displaystyle \left[\begin{matrix}0\\\dot{\theta} r{\left(t \right)}\end{matrix}\right]\)
The first element is the radial velocity and the second element is the tangential velocity.
rr_ddot_polar.subs(sp.symbols(r'\dot{r}'), 0).subs(sp.symbols(r'\ddot{r}'), 0)\(\displaystyle \left[\begin{matrix}- \dot{\theta}^{2} r{\left(t \right)}\\\ddot{\theta} r{\left(t \right)}\end{matrix}\right]\)
The first term is the centripetal acceleration, always pointing towards the center of rotation, and the second term is the tangential acceleration, which is zero if the angular velocity is constant.
To summarize planar circular motion: with a constant \(r\), we have \(\dot{r} = \ddot{r} = 0\) and \(\mathbb{r} \perp \dot{\mathbb{r}}\) since \(\frac{d}{dt}(r^2) = 0 = \frac{d}{dt}(\mathbb{r} \cdot \mathbb{r}) = \dot{\mathbb{r}} \cdot \mathbb{r} + \mathbb{r} \cdot \dot{\mathbb{r}} \Leftrightarrow \mathbb{r} \cdot \dot{\mathbb{r}} = 0\). Note especially the centripetal acceleration \(-r\dot{\theta}^2\mathbb{e}_r\), it is directed in the negative direction of \(\mathbb{e}_r\), towards the center of the circle.
If the rotational velocity is constant, i.e., if \(\omega = \dot{\theta} = \text{constant}\) or \(\alpha = \ddot{\theta} = 0\) then we only have centripetal acceleration, caused by a curvilinear motion and perpendicular to the tangent direction, directed inwards, towards the center.
\[ \begin{cases} \mathbb{r}(t) = r\,\mathbb{e}_r \\ \mathbb{v}(t) = r\dot{\theta}\mathbb{e}_\theta \\ \mathbb{a}(t) = -r\dot{\theta}^2\mathbb{e}_r + r\ddot{\theta}\mathbb{e}_\theta \end{cases} \]
especially when \(\omega = \dot{\theta} = \text{constant}\)
\[ \begin{cases} \mathbb{r}(t) = r\,\mathbb{e}_r \\ \mathbb{v}(t) = r\omega\mathbb{e}_\theta \\ \mathbb{a}(t) = -r\omega^2\mathbb{e}_r = \{v = |\mathbb{v}| = r\omega\} = -\frac{v^2}{r}\mathbb{e}_r \end{cases} \]
Note that we can express bases also using the cross product
\[ \mathbb{e}_\theta = \mathbb{e}_z \times \mathbb{e}_r \]
which is very convenient since circular motion is common and we often have angular velocity as input data, thus we can formulate
\[ \mathbb{v}(t) = \dot{\theta}\mathbb{e}_z \times \mathbb{r} = r\dot{\theta}\mathbb{e}_\theta = r\omega\mathbb{e}_\theta \]
Another easy way is to explicitly rotate the base \(\mathbb{e}_r\) by 90 degrees into the direction of curvature, the drawback is we need to know which direction to rotate into, which is why the previous method is preferred, since the direction of rotation is determined by the sign of the angular velocity \(\omega\).
Many times we can express position and motion (velocities and acceleration) relative to some other point rather than the origin. This is done such that the modeling becomes significantly easier and less error prone. Besides, in real world applications, engineers need to make measurements and many times it is more accurate to make relative measurements rather than absolute. This approach is called relative motion analysis and will be used extensively in rigid body dynamics.
If we know the position of particle \(B\) as well as the relative position of particle \(A\) with respect to \(B\), that is \(\mathbb{r}_{A/B}\) (Can be read: A seen from B) we can express the absolute position of \(A\) as
\[ \mathbb{r}_A = \mathbb{r}_B + \mathbb{r}_{A/B} \]
Again it might be more convenient to express \(\mathbb{r}_{A/B}\) and \(\mathbb{r}_B\) than \(\mathbb{r}_A\) directly.
If the position vectors are expressed as variables of time, we can differentiate them with respect to time and get
\[ \frac{d}{dt}(\mathbb{r}_A = \mathbb{r}_B + \mathbb{r}_{A/B}) = \dot{\mathbb{r}}_A = \dot{\mathbb{r}}_B + \dot{\mathbb{r}}_{A/B} \quad \text{or} \quad \mathbb{v}_A = \mathbb{v}_B + \mathbb{v}_{A/B} \]
and
\[ \ddot{\mathbb{r}}_A = \ddot{\mathbb{r}}_B + \ddot{\mathbb{r}}_{A/B} \quad \text{or} \quad \mathbb{a}_A = \mathbb{a}_B + \mathbb{a}_{A/B} \]
Passengers in the A380 flying east at a speed of 800 km/h observe a stunning JAS 39 Gripen passing underneath with a heading of 45\(^\circ\) in the north-east direction, although to the passengers it appears that the JAS is moving away at an angle of 60\(^\circ\) as shown in the figure. Determine the true velocity of the JAS 39 Gripen.
Solution
We have two unknows, the absolute speed of the JAS, \(v_B\) and the relative speed observed from the A380, \(v_{B/A}\).
We formulate the relative motion equation to suit our problem \[ \mathbb{v}_B = \mathbb{v}_A + \mathbb{v}_{B/A} \]
Using the known directions we can express the matrix equation
\[ v_{B}\begin{bmatrix}\cos45^{\circ}\\ \sin45^{\circ} \end{bmatrix}=\begin{bmatrix}800\\ 0 \end{bmatrix}+v_{B/A}\begin{bmatrix}-\cos60^{\circ}\\ \sin60^{\circ} \end{bmatrix} \]
import sympy as sp
v_B, v_B_A = sp.symbols('v_B v_B_A', real=True, positive=True)The velocity of the A380 is given as \(v_A = 800\) km/h in the east direction.
vv_B = sp.Matrix([v_B * sp.cos(45*sp.pi/180), v_B * sp.sin(45*sp.pi/180)])
vv_B \(\displaystyle \left[\begin{matrix}\frac{\sqrt{2} v_{B}}{2}\\\frac{\sqrt{2} v_{B}}{2}\end{matrix}\right]\)
vv_A = sp.Matrix([800, 0])
vv_A\(\displaystyle \left[\begin{matrix}800\\0\end{matrix}\right]\)
vv_B_A = v_B_A * sp.Matrix([-sp.cos(60*sp.pi/180), sp.sin(60*sp.pi/180)])
vv_B_A\(\displaystyle \left[\begin{matrix}- \frac{v_{B A}}{2}\\\frac{\sqrt{3} v_{B A}}{2}\end{matrix}\right]\)
EQ = sp.Eq(vv_B, vv_A + vv_B_A)
EQ\(\displaystyle \left[\begin{matrix}\frac{\sqrt{2} v_{B}}{2}\\\frac{\sqrt{2} v_{B}}{2}\end{matrix}\right] = \left[\begin{matrix}800 - \frac{v_{B A}}{2}\\\frac{\sqrt{3} v_{B A}}{2}\end{matrix}\right]\)
sol = sp.solve(EQ, (v_B, v_B_A))
sol[v_B].evalf(3)\(\displaystyle 717.0\)
sol[v_B_A].evalf(3)\(\displaystyle 586.0\)
To conclude, the JAS is moving away from the passengers at 586 km/h and its absolute speed is 717 km/h.
The acceleration performance of a new solar car is being tested at the Jönköping airfield. The JTH students are measuring time against velocity and the data is presented as two vectors. Create a suitable continuous function of the speed and determine how far the car has traveled at the last measurement point. Plot the distance - time curve as well as the acceleration curve.
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt
# Given data:
T = np.array([1.1486, 2.0569, 3.1485, 4.0443, 5.0165, 6.2552, 7.1682, 8.2789, 9.209, 10.175])
V = np.array([9.1023, 16.245, 22.732, 26.446, 33.052, 36.889, 41.304, 43.496, 45.861, 48.319])
# Fit a quadratic function to the data
c = np.polyfit(T, V, 2)
t = sp.Symbol('t')
v = c[0]*t**2 + c[1]*t + c[2]
print(f"v(t) = {sp.pretty(v)}")
# Plot the data and the quadratic function
fig, ax = plt.subplots()
ax.plot(T, V, 'ro', label='Data')
tn = np.linspace(0, max(T), 100)
vn = np.polyval(c, tn)
ax.plot(tn, vn, 'b-', label='Fit')
ax.set(xlabel='Time, t [s]', ylabel='Velocity, v [m/s]', title='Velocity vs. Time')
ax.grid(True)
ax.legend()
ax.fill_between(tn, vn, color='b', alpha=0.2)
plt.show()
# Acceleration is the derivative of velocity
a = sp.diff(v, t)
print(f"a(t) = {sp.pretty(a)}")
# Plot acceleration
an = sp.lambdify(t, a, 'numpy')(tn)
fig, ax = plt.subplots()
ax.plot(tn, an, 'm-')
ax.set(xlabel='Time, t [s]', ylabel='Acceleration, a [m/s^2]', title='Acceleration vs. Time')
ax.grid(True)
plt.show()
# Total distance is the integral of velocity
s = sp.integrate(v, (t, 0, max(T)))
print(f"Total distance traveled: {s:.2f} m")v(t) = 2
- 0.309736757594532⋅t + 7.75454224391216⋅t + 1.05549793621097
a(t) = 7.75454224391216 - 0.619473515189064⋅t
Total distance traveled: 303.39 m
Study the motion of the slider-crank mechanism. Plot the motion of A, G and B as well as their velocities and accelerations as a function of \(\theta\). Let the angular velocity \(\omega=25\) rad/s.
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
t = sp.Symbol('t')
theta = sp.Function('theta')(t)
L, r = sp.symbols('L r', positive=True)
h = r*sp.sin(theta)
a = r*sp.cos(theta)
b = sp.sqrt(L**2 - h**2)
r_A = sp.Matrix([-(a+b), 0])
r_B = r*sp.Matrix([-sp.cos(theta), sp.sin(theta)])
r_BA = r_A - r_B
r_G = r_B + 100/(100+250)*r_BA
theta_dot, theta_ddot = sp.symbols('theta_dot theta_ddot', positive=True)
theta_sym = sp.Symbol('theta', positive=True)
r_A_dot = sp.diff(r_A, t).subs([(sp.diff(theta,t), theta_dot), (theta, theta_sym)])
r_B_dot = sp.diff(r_B, t).subs([(sp.diff(theta,t), theta_dot), (theta, theta_sym)])
r_G_dot = sp.diff(r_G, t).subs([(sp.diff(theta,t), theta_dot), (theta, theta_sym)])
r_A_ddot = sp.diff(r_A, t, 2).subs([(sp.diff(theta,t,2), theta_ddot), (sp.diff(theta,t), theta_dot), (theta, theta_sym)])
r_B_ddot = sp.diff(r_B, t, 2).subs([(sp.diff(theta,t,2), theta_ddot), (sp.diff(theta,t), theta_dot), (theta, theta_sym)])
r_G_ddot = sp.diff(r_G, t, 2).subs([(sp.diff(theta,t,2), theta_ddot), (sp.diff(theta,t), theta_dot), (theta, theta_sym)])
# Substitute numerical data
num_data = {L: (250+100)/1000, r: 125/1000, theta_dot: 25, theta_ddot: 0}
r_A_dot_m = r_A_dot.norm().subs(num_data)
r_B_dot_m = r_B_dot.norm().subs(num_data)
r_G_dot_m = r_G_dot.norm().subs(num_data)
r_A_ddot_m = r_A_ddot[0].subs(num_data) # Acceleration is only in x
r_B_ddot_m = r_B_ddot.norm().subs(num_data)
r_G_ddot_m = r_G_ddot.norm().subs(num_data)
# Lambdify for plotting
vA_func = sp.lambdify(theta_sym, r_A_dot_m, 'numpy')
vB_func = sp.lambdify(theta_sym, r_B_dot_m, 'numpy')
vG_func = sp.lambdify(theta_sym, r_G_dot_m, 'numpy')
aA_func = sp.lambdify(theta_sym, r_A_ddot_m, 'numpy')
aB_func = sp.lambdify(theta_sym, r_B_ddot_m, 'numpy')
aG_func = sp.lambdify(theta_sym, r_G_ddot_m, 'numpy')
theta_vals = np.linspace(0, 2*np.pi, 360)
# Plot Speed
plt.figure(figsize=(10, 6))
plt.plot(np.rad2deg(theta_vals), vA_func(theta_vals), label='$|\\dot{\\mathbb{r}}_A|$')
plt.plot(np.rad2deg(theta_vals), vB_func(theta_vals), label='$|\\dot{\\mathbb{r}}_B|$')
plt.plot(np.rad2deg(theta_vals), vG_func(theta_vals), label='$|\\dot{\\mathbb{r}}_G|$')
plt.xlabel('$\\theta$ [degrees]')
plt.ylabel('Speed [m/s]')
plt.title('Crank-slider speed')
plt.xticks(np.arange(0, 361, 45))
plt.grid(True)
plt.legend()
plt.show()
# Plot Acceleration
plt.figure(figsize=(10, 6))
plt.plot(np.rad2deg(theta_vals), aA_func(theta_vals), label='$\\ddot{r}_{Ax}$')
plt.plot(np.rad2deg(theta_vals), aB_func(theta_vals), label='$|\\ddot{\\mathbb{r}}_B|$')
plt.plot(np.rad2deg(theta_vals), aG_func(theta_vals), label='$|\\ddot{\\mathbb{r}}_G|$')
plt.xlabel('$\\theta$ [degrees]')
plt.ylabel('Acceleration [m/s$^2$]')
plt.title('Crank-slider acceleration')
plt.xticks(np.arange(0, 361, 45))
plt.grid(True)
plt.legend()
plt.show()
