Working in 3D is not too difficult once we get familiar with proper tools. Here we shall give a thorough walk-through of the process on a non-trivial example. The hardest part is the kinematics, i.e., describing the positions of the points relative to each other. Once this is done, formulating forces and moments to do equilibrium is a trivial matter with computer support.
We shall here learn to work with a rotated coordinate system.
We start out by defining the coordinate system, the red axis is the x-axis, the green axis is the y-axis and the blue axis is the z-axis. See Figure 9.1.
\(A=(30, 60, 75)\) mm, \(r_{AB}=65\) mm, \(r_{BC}=50\) mm, \(\theta_1=35^\circ\) around the local \(x-\) axis at \(A\) and \(E=(145,40,0)\).
Let \(\theta_2\) denote the rotation of the shaft \(AB\). In the picture above, the \(\mathbb{r}_{BC}\) vector is parallel to the \(x-\) axis, and this is when \(\theta_2=0^\circ\).
The center of gravity can be defined using local coordinates as: \(\mathbb r_{AG} = 34\tilde x + 51 \tilde z\), where \(\tilde x\) and \(\tilde z\) are the unit vectors in the local coordinate system according to Figure 9.2.:
Similarly, we can define \(\mathbb r_{CD}=70 \tilde x + 40 \tilde z\) using local coordinates.
Given the spring constant \(k=0.1\) N/mm and the unloaded spring length \(L_0=120\) mm:
We are going to work with vectors to define the path from \(O\) to \(D\) and \(G\) as a function of the rotation angle \(\theta\). To do this, we need to introduce local coordinates and rotate these using rotation matrices:
\[ \left. \mathbb R_x\left(\theta\right)=\left[\begin{array}{ccc}1&0&0\\0&\cos\theta&-\sin\theta\\0&\sin\theta&\cos\theta\end{array}\right.\right] \]
\[ \left. \mathbb R_y\left(\theta\right)=\left[\begin{array}{ccc}\cos\theta&0&\sin\theta\\0&1&0\\-\sin\theta&0&\cos\theta\end{array}\right.\right] \]
\[ \mathbb R_{z}\left(\theta\right)=\begin{bmatrix}\cos\theta&-\sin\theta&0\\\sin\theta&\cos\theta&0\\0&0&1\end{bmatrix} \]
We start by getting to point \(A\) using the vector \(\mathbb r_{OA}\). From here we need to establish the direction of the shaft axis, \(AB\), which we can do by rotating the \(\mathbb e_z\) vector around the \(x-\) axis by \(-\theta_1\) degrees (using the right hand rule).
This gives us our \(\mathbb e_{AB}\) direction vector:
\[ \mathbb e_{AB} = \tilde z = \mathbb R_x(-\theta_1) \mathbb e_z \]
import sympy as sp
theta_1 = sp.symbols('theta_1', real=True)
Rx = lambda v: sp.rot_ccw_axis1(v) # x-axis rotation matrix
Ry = lambda v: sp.rot_ccw_axis2(v) # y-axis rotation matrix
Rz = lambda v: sp.rot_ccw_axis3(v) # z-axis rotation matrix
e_z = sp.Matrix([0, 0, 1]) # unit vector along z-axis
e_AB = (Rx(-theta_1) @ e_z)Sympy has built in functions for rotation matrices. Note, however, that we need to use rot_ccw_axis1, rot_ccw_axis2 and rot_ccw_axis3 which adhere to the right-hand-rule convention! See sympy documentation for more information.
Note that we rotate \(-\theta_1\) due to the right hand rule, i.e., put the thumb in the direction of \(x\) and the fingers will point in the positive direction of rotation.
The shaft \(AB\) needs to rotate around its axis, or the \(\tilde{z}-\) axis. The rotation can be achieved by rotating around the global axis, \(x,y\) and \(z\). From the initial confguration, see below, where the shaft axis is aligned with the \(z\)-axis and \(\mathbb r_{BC}\) oriented towards the \(x-\) axis. To get the correct rotation of the shaft, we need to first rotate around the \(z-\) axis and then around the \(x-\) axis.
Mathematically, this is the same as rotating our coordinate system by
\[ \tilde{\mathbb{X}}=[\tilde{x},\tilde{y},\tilde{z}]=\mathbb{R}_x\left(-\theta_1\right)\mathbb{R}_z\left(\theta_2\right)\mathbb{I} \]
where \(\mathbb I\) is the identity matrix containing the Cartesian bases, i.e.,
\[ [\mathbb e_x,\mathbb e_y,\mathbb e_z]=\begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix} \]
Let us verify this by setting \(\theta_1=0\) and \(\theta_2=0\), we expect \(\tilde{\mathbb{X}} = \mathbb I\)
theta_1 = 0;
theta_2 = 0;
Rx(theta_1) @ Rz(theta_2) @ sp.eye(3)\(\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]\)
Good! Now, we can set \(\theta_2 = 90^\circ\) and we should expect to get \(\tilde x = [0, 1, 0]^\mathsf T\)
theta_1 = 0; theta_2 = 90*sp.pi/180;
Rx(theta_1) @ Rz(theta_2) @ sp.eye(3)\(\displaystyle \left[\begin{matrix}0 & -1 & 0\\1 & 0 & 0\\0 & 0 & 1\end{matrix}\right]\)
Here, we get that the \(\tilde x\) -axis points in the positive \(y\) direction, the \(\tilde y\) -axis points in the negative \(x\) direction and the \(\tilde z\) -axis points in the positive \(z\) direction. This is correct, since we have rotated the coordinate system 90 degrees around the \(z\)-axis according to the right-hand-rule.
Now, we can also rotate around the \(x-\) axis by \(\theta_1 = -90^\circ\), we should get \(\tilde z = [0, 1, 0]^\mathsf T\)
theta_1 = -90*sp.pi/180; theta_2 = 90*sp.pi/180;
Rx(theta_1) @ Rz(theta_2) @ sp.eye(3)\(\displaystyle \left[\begin{matrix}0 & -1 & 0\\0 & 0 & 1\\-1 & 0 & 0\end{matrix}\right]\)
We can confirm that the new coordinates point in the directions we expect. The \(\tilde x\) -axis points in the negative \(z\) direction, the \(\tilde y\) -axis points in the negative \(x\) direction and the \(\tilde z\) -axis points in the positive \(y\) direction. This is correct.
Thus have confirmed that the rotation makes sense, we have built confidence in our model and can set the rotation to solve our problem:
theta_1 = -35*sp.pi/180; theta_2 = 0*sp.pi/180;
X = Rx(theta_1) @ Rz(theta_2) @ sp.eye(3)
X = X.evalf(4)
eex = X[:, 0]
eey = X[:, 1]
eez = X[:, 2]
X\(\displaystyle \left[\begin{matrix}1.0 & 0 & 0\\0 & 0.8192 & 0.5736\\0 & -0.5736 & 0.8192\end{matrix}\right]\)
Now we can easily model the rest of the points.
\[ \begin{align*} \mathbb{r}_{AB} & = r_{AB}\tilde{z}\\ \mathbb{r}_{BC} & = r_{BC}\tilde{x}\\ \mathbb{r}_{CD} & = 70\tilde{x}+40\tilde{z}\\ \mathbb{r}_{AG} & = 34\tilde{x}+51\tilde{z}\\ \mathbb{r}_{DE} & = \mathbb{r}_{OE}-\mathbb{r}_{OD} \end{align*} \]
We do a final verification of the kinematics by visualizing the geometry in Python:
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from sympy import pi
Rx = lambda v: sp.rot_ccw_axis1(v) # x-axis rotation matrix
Ry = lambda v: sp.rot_ccw_axis2(v) # y-axis rotation matrix
Rz = lambda v: sp.rot_ccw_axis3(v) # z-axis rotation matrix
theta_2 = 26;
r_AB = 65; r_BC = 50;
X = Rx(-35*pi/180) @ Rz(theta_2*pi/180) @ sp.eye(3)
eex = X[:, 0].evalf(4)
eey = X[:, 1].evalf(4)
eez = X[:, 2].evalf(4)
rr_O = sp.Matrix([0, 0, 0])
rr_OA = sp.Matrix([30, 60, 75])
rr_OE = sp.Matrix([175, 105, 15])
rr_AB = r_AB * eez
rr_BC = r_BC * eex
rr_CD = eex * 70 + eez * 40
rr_OD = rr_OA +rr_AB + rr_BC + rr_CD
rr_OB = rr_OA + rr_AB
rr_OC = rr_OB + rr_BC
rr_D = rr_OD - rr_OA
rr_DE = rr_OE - rr_OD
rr_AG = 34 * eex + 51 * eez
rr_OG = rr_OA + rr_AG
# Visualization
offset = 0.05
plt.style.use('seaborn-v0_8-whitegrid')
fig = plt.figure(figsize=(7.25, 7.25))
ax = fig.add_subplot(111,projection='3d')
ax.plot(*zip(rr_OA,rr_OB), '-k', linewidth=2)
ax.text(*rr_OA, 'A')
ax.text(*rr_OB, 'B')
ax.plot(*zip(rr_OB,rr_OC), '-k', linewidth=2)
ax.text(*rr_OC, 'C')
ax.plot(*zip(rr_OC,rr_OD), '-k', linewidth=2)
ax.text(*rr_OD, 'D')
ax.plot(*zip(rr_OD,rr_OE), '-r', linewidth=2)
ax.text(*rr_OE, 'E')
ax.plot(*zip(rr_OG), 'ok', markerfacecolor='k')
ax.text(*rr_OG, 'G')
ax.view_init(elev=3) # Default 3D view
ax.grid(True)
ax.set_xlim([0, 200])
ax.set_ylim([-100, 250])
ax.set_zlim([0, 250])
ax.view_init(elev=20, azim=-45)
ax.set_title(rf'$\theta_2={theta_2}^\circ$', fontsize=14) # Using LaTeX for title
plt.show()
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, PillowWriter
from sympy import pi
Rx = lambda v: sp.rot_ccw_axis1(v)
Rz = lambda v: sp.rot_ccw_axis3(v)
r_AB = 65
r_BC = 50
# Full sweep from -90 to 90 and back
theta_range = list(range(-90, 91)) + list(range(89, -91, -1))
#fig = plt.figure(figsize=(7.25, 7.25))
#ax = fig.add_subplot(111, projection='3d')
def update(frame_idx):
theta_2 = theta_range[frame_idx]
ax.cla()
X = Rx(-35*pi/180) @ Rz(theta_2*pi/180) @ sp.eye(3)
eex = X[:, 0].evalf()
eey = X[:, 1].evalf()
eez = X[:, 2].evalf()
rr_O = sp.Matrix([0, 0, 0])
rr_OA = sp.Matrix([30, 60, 75])
rr_OE = sp.Matrix([175, 105, 15])
rr_AB = r_AB * eez
rr_BC = r_BC * eex
rr_CD = eex * 70 + eez * 40
rr_OD = rr_OA + rr_AB + rr_BC + rr_CD
rr_OB = rr_OA + rr_AB
rr_OC = rr_OB + rr_BC
rr_AG = 34 * eex + 51 * eez
rr_OG = rr_OA + rr_AG
def plot_line(p1, p2, style='-k'):
ax.plot(
[p1[0], p2[0]],
[p1[1], p2[1]],
[p1[2], p2[2]],
style,
linewidth=2
)
plot_line(rr_OA, rr_OB)
plot_line(rr_OB, rr_OC)
plot_line(rr_OC, rr_OD)
plot_line(rr_OD, rr_OE, '-r')
ax.plot([rr_OG[0]], [rr_OG[1]], [rr_OG[2]], 'ok', markerfacecolor='k')
ax.text(*rr_OA, 'A')
ax.text(*rr_OB, 'B')
ax.text(*rr_OC, 'C')
ax.text(*rr_OD, 'D')
ax.text(*rr_OE, 'E')
ax.text(*rr_OG, 'G')
azim = -45 + frame_idx * 0.2 # Slowly rotate the camera clockwise
ax.view_init(elev=20, azim=azim)
ax.set_xlim([0, 200])
ax.set_ylim([-100, 250])
ax.set_zlim([0, 250])
ax.set_title(rf'$\theta_2={theta_2}^\circ$', fontsize=14)
ax.grid(True)
#ani = FuncAnimation(fig, update, frames=len(theta_range), interval=50)
# Save to GIF
#ani.save('graphics/ForceVectorsIn3D_anim1.gif', writer=PillowWriter(fps=20))
The spring force is modeled using
\[ L = r_{DE}, \quad \mathbb e_{DE} = \dfrac{\mathbb r_{DE}}{L}, \quad L_0 = 120 \text{ mm} \]
\[ F = k(L-L_0), \text{ and } \mathbb F = F \mathbb e_{DE} \]
The moment at \(A\) is given by
\[ \mathbb M_{A} = \mathbb r_{AD} \times \mathbb F + \mathbb r_{AG} \times \mathbb W \]
where \(\mathbb W = [0,0,-mg]^\mathsf T\), \(m=150\) g, and
\[ \mathbb r_{AD} = \mathbb r_{OD} - \mathbb r_{OA} \]
Finally the moment around the axcis of the motor is given by
\[ M_{AB} = \mathbb M_{A} \cdot \tilde z \]
Now we can plot the graphs using a symbolic \(\theta_2\)
import sympy as sp
theta_1 = -35*sp.pi/180;
theta_2 = sp.symbols('theta_2', real=True)
Rx = lambda v: sp.rot_ccw_axis1(v) # x-axis rotation matrix
Ry = lambda v: sp.rot_ccw_axis2(v) # y-axis rotation matrix
Rz = lambda v: sp.rot_ccw_axis3(v) # z-axis rotation matrix
e_z = sp.Matrix([0, 0, 1]) # unit vector along z-axis
X = Rx(theta_1) @ Rz(theta_2) @ sp.eye(3)
X = X.evalf(4)
eex = X[:, 0]
eey = X[:, 1]
eez = X[:, 2]
r_AB = 65
r_BC = 50
k = 0.1;
rr_AB = r_AB * eez
rr_BC = r_BC * eex
rr_CD = eex * 70 + eez * 40
rr_OA = sp.Matrix([30, 60, 75])
rr_OE = sp.Matrix([175, 105, 15])
rr_OD = rr_OA +rr_AB + rr_BC + rr_CD
rr_OB = rr_OA + rr_AB
rr_OC = rr_OB + rr_BC
rr_D = rr_OD - rr_OA
rr_DE = rr_OE - rr_OD
rr_AG = 34 * eex + 51 * eez
rr_OG = rr_OA + rr_AG
rr_O = sp.Matrix([0, 0, 0])
rr_AD = rr_OD - rr_OA
L = rr_DE.norm()
e_DE = rr_DE / L
L0 = 120;
F = k*(L-L0);
FF = F*e_DE;
WW = sp.Matrix([0, 0, -0.15*9.82])
import numpy as np
import matplotlib.pyplot as plt
from numpy import pi
plt.style.use('seaborn-v0_8-whitegrid')
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot()
theta_2_np_rad = np.linspace(-90*pi/180,90*pi/180, 180)
theta_2_np_deg = np.linspace(-90,90, 180)
L_f = sp.lambdify(theta_2, L)
ax.plot(theta_2_np_deg, L_f(theta_2_np_rad), 'b', linewidth=2)
ax.set_xlim([-90, 90])
ax.set_xlabel(r'$\theta_2$ [deg]', fontsize=14)
ax.set_ylabel(r'$L(\theta)$ [mm]', fontsize=14)
plt.grid(True)
plt.title(r'$L$ vs. $\theta_2$', fontsize=16)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
from numpy import pi
plt.style.use('seaborn-v0_8-whitegrid')
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot()
theta_2_np_rad = np.linspace(-90*pi/180,90*pi/180, 180)
theta_2_np_deg = np.linspace(-90,90, 180)
F_f = sp.lambdify(theta_2, F)
ax.plot(theta_2_np_deg, F_f(theta_2_np_rad), 'b', linewidth=2)
ax.set_xlim([-90, 90])
ax.set_xlabel(r'$\theta_2$ [deg]', fontsize=14)
ax.set_ylabel(r'$F(\theta)$ [N]', fontsize=14)
plt.grid(True)
plt.title(r'$F$ vs. $\theta_2$', fontsize=16)
plt.show()
MM_A = rr_AD.cross(FF) + rr_AG.cross(WW)
e_AB = eez
M_AB = MM_A.dot(e_AB)import numpy as np
import matplotlib.pyplot as plt
from numpy import pi
plt.style.use('seaborn-v0_8-whitegrid')
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot()
theta_2_np_rad = np.linspace(-90*pi/180,90*pi/180, 180)
theta_2_np_deg = np.linspace(-90,90, 180)
M_AB_f = sp.lambdify(theta_2, M_AB)
ax.plot(theta_2_np_deg, M_AB_f(theta_2_np_rad), 'b', linewidth=2)
ax.set_xlim([-90, 90])
ax.set_xlabel(r'$\theta_2$ [deg]', fontsize=14)
ax.set_ylabel(r'$M_{{AB}}(\theta)$ [Nmm]', fontsize=14)
plt.grid(True)
plt.title(r'$Nmm$ vs. $\theta_2$', fontsize=16)
plt.show()
Zooming in at the force graph, we learn that the minimum force and length of the spring occurs at \(26^\circ\).
The maximum moment occurs at around \(-78^\circ\). The moment is zero around \(29.5^\circ\)
Another approach to rotating the coordinate system is to rotate around a local coordinate system. This makes the modeling and thinking easier. One does not need to be confined to thinking about rotations around global coordinates.
This can be done by rotating the \(z\) rotation matrix, \(\mathbb R_z\) around the \(x-\) axis by
\[ \mathbb R_{\tilde{z}}(\theta_2) = \mathbb R_x(\theta_1) \mathbb R_z(\theta_2) \]
This rotation matrix rotates vectors around the \(\tilde{z}\) axis. This approach still relies on defning the rotation matrix by rotating it around global axis. The order of rotations is important, thus chaining rotations can become tedious.
It would be nicer to have a function that lets us rotate around a arbitrary vector \(\boldsymbol u\). To achieve this we visit our favorite source of goodyness, Wikipedia, and find that the very elegant way of rotating a vector \(\boldsymbol v\) around an axis \(\boldsymbol u\) by an angle \(\theta\) according to the right-hand rule is stated as
\[ \boldsymbol{R}_u\boldsymbol{v}=\boldsymbol{u}\left(\boldsymbol{u}\cdot\boldsymbol{v}\right)+\cos\theta\left(\boldsymbol{u}\times\boldsymbol{v}\right)\times u+\sin\theta\left(\boldsymbol{u}\times\boldsymbol{v}\right) \]
import sympy as sp
def cosd(a):
return sp.cos(a*sp.pi/180)
def sind(a):
return sp.sin(a*sp.pi/180)
Ruv = lambda theta, u, v: u*u.dot(v) + cosd(theta)*u.cross(v).cross(u) + sind(theta)*u.cross(v)We can test this by rotating the \(x-\) axis 90 degrees around the \(z-\) axis. We should get the \(y-\) axis, i.e., [0, 1, 0]\(^\mathsf T\).
Ruv(90, sp.Matrix([0,0,1]), sp.Matrix([1,0,0]))\(\displaystyle \left[\begin{matrix}0\\1\\0\end{matrix}\right]\)
Now we can repeat the definitions of the points from the example above using the new rotation tool.
rr_OA = sp.Matrix([30, 60, 75])
rr_OE = sp.Matrix([175, 105, 15])
r_AB = 65
r_BC = 50
k = 0.1
L0 = 120
theta_1 = 35
theta_2 = 26We formulate the direction \(\mathbb e_{AB}\) by rotating the \(z-\) axis around the \(x-\) axis.
e_AB = Ruv(-theta_1, sp.Matrix([1,0,0]), sp.Matrix([0,0,1])).evalf(4)
e_AB\(\displaystyle \left[\begin{matrix}0\\0.5736\\0.8192\end{matrix}\right]\)
Ok, nothing new, hardly easier than using standard rotation matrices. But now we can rotate around \(\mathbb e_{AB}\)
ee_x = Ruv(theta_2, e_AB, sp.Matrix([1,0,0])).evalf(4)
ee_x\(\displaystyle \left[\begin{matrix}0.8988\\0.3591\\-0.2514\end{matrix}\right]\)