import sympy as sm
import sympy.physics.mechanics as me
me.init_vprinting()
L = me.ReferenceFrame('L') # The laboratory frame
C = me.ReferenceFrame('C') # The disc frame
R = me.ReferenceFrame('R') # The rod frame
q1, q2 = me.dynamicsymbols('q1, q2')
r = sm.symbols('r', real=True)
t = me.dynamicsymbols._t
C.orient(L, 'Axis', (q1, L.x))
R.orient(C, 'Axis', (q2, C.y))
O, P = sm.symbols('O, P', cls=me.Point)
P.set_pos(O, r*C.y - 3*r*R.x)
P.pos_from(O).express(L)
If you express in L and differeniate each measure number you get these measure numbers of P in L:
L_v_P = (P.pos_from(O).dot(L.x).diff(t) * L.x +
P.pos_from(O).dot(L.y).diff(t) * L.y +
P.pos_from(O).dot(L.z).diff(t) * L.z)
L_v_P
L_v_P.express(C).simplify()
Same as doing:
P.pos_from(O).dt(L).express(C).simplify()
L_a_P = (L_v_P.dot(L.x).diff(t) * L.x +
L_v_P.dot(L.y).diff(t) * L.y +
L_v_P.dot(L.z).diff(t) * L.z)
L_a_P
L_a_P.express(C).simplify()
C_v_P = (P.pos_from(O).dot(C.x).diff(t) * C.x +
P.pos_from(O).dot(C.y).diff(t) * C.y +
P.pos_from(O).dot(C.z).diff(t) * C.z)
C_v_P
C_a_P = (C_v_P.dot(C.x).diff(t) * C.x +
C_v_P.dot(C.y).diff(t) * C.y +
C_v_P.dot(C.z).diff(t) * C.z)
C_a_P
Setup the problem like 2.7.
q1, q2, q3, q4, q5 = me.dynamicsymbols('q1:6')
R = sm.symbols('R')
A = me.ReferenceFrame('A')
Y = me.ReferenceFrame('Y') # intermediate frame so that simple rotations can be used
B = me.ReferenceFrame('B')
C = me.ReferenceFrame('C')
Y.orient(A, 'Axis', (q1, A.z))
B.orient(Y, 'Axis', (sm.pi/2 - q2, Y.x))
C.orient(B, 'Axis', (q3, B.z))
Define some new variables for use in defining velocities.
u1, u2, u3, u4, u5 = me.dynamicsymbols('u1:6')
Create the points of interest.
O, P, Cs = sm.symbols('O, P, C^*', cls=me.Point)
Set the relative positions of the points.
P.set_pos(O, q4 * A.x + q5 * A.y)
Cs.set_pos(P, R * B.y)
$u_1,u_2,u_4$ are equal to the measure numbers of the angular velocity of C in A as described in Problem 2.7.
C.ang_vel_in(A).express(B)
Set the velocity of P in A using the u's.
P.set_vel(A, u4*A.x + u5*A.y)
$C^*$ and $P$ are both fixed in B, so apply the two point theory:
Cs.v2pt_theory(P, A, B).express(B).simplify()
Note that this has qdots present and we only want this velocity interms of the u's. So, solve for the qdots and substitute them in.
eqs = [C.ang_vel_in(A).dot(B.x) - u1,
C.ang_vel_in(A).dot(B.y) - u2,
C.ang_vel_in(A).dot(B.z) - u3]
qdots = sm.solve(eqs, q1.diff(), q2.diff(), q3.diff())
Cs.v2pt_theory(P, A, B).express(B).subs(qdots).simplify()
A_v_Cs = Cs.v2pt_theory(P, A, B).subs(qdots).express(B)
A_v_Cs
Time differentiate the velcity vector (use the one with only u's) and then substitute in the qdots.
A_v_Cs.dt(A).subs(qdots).simplify().to_matrix(B)
Ch = me.Point('C_h')
Ch.set_pos(Cs, -R*B.y)
A_v_Ch = Cs.vel(A) + me.cross(C.ang_vel_in(A), Ch.pos_from(Cs))
A_v_Ch.express(A)
A_v_Ch.express(A).subs(qdots).simplify()
v2pt_theory()
gives the same result:
Ch.v2pt_theory(Cs, A, C).express(A).subs(qdots).simplify()
A_a_Ch = (Cs.acc(A) + me.cross(C.ang_acc_in(A), Ch.pos_from(Cs)) +
me.cross(C.ang_vel_in(A), me.cross(C.ang_vel_in(A), Ch.pos_from(Cs))))
A_a_Ch.express(B).subs(qdots).simplify()
Ch.a2pt_theory(Cs, A, C).express(B).subs(qdots).simplify()
Key things to note:
These are constant with respect to time:
theta, b, r = sm.symbols('theta, b, r')
These three reference frames are needed.
R, S, C = sm.symbols('R, S, C', cls=me.ReferenceFrame)
Create some measure numbers for the angular velocities that we will eventually solve for.
wc, ws1, ws2, ws3 = me.dynamicsymbols('omega_c, omega_{s1}, omega_{s2}, omega_{s3}')
The cone spins about the vertical axis. Note that all velocities will be define using the C unit vectors.
C.set_ang_vel(R, wc*C.y)
The angular velocity of the sphere can be anything if the sphere slips and rolls on the surfaces (the general case).
S.set_ang_vel(R, ws1*C.x + ws2*C.y + ws3*C.z)
Create some points that are needed and defined above.
Pc, Ps, Ss, A, B, O = sm.symbols('P_c, P_s, S^*, A, B, O', cls=me.Point)
Now, the positions of A and B can be set relative to $S^*$, the sphere's center. This will allow us to use the two point theorem to calculate velocities since $S^*,A,B$ all all fixed in S.
A.set_pos(Ss, r*C.x)
B.set_pos(Ss, -r*C.y)
We know that:
$$ {}^R\mathbf{v}^A = 0 = {}^R\mathbf{v}^{S^*} + {}^R\omega^S \times \mathbf{r}^{A/S^*} \\ {}^R\mathbf{v}^B = 0 = {}^R\mathbf{v}^{S^*} + {}^R\omega^S \times \mathbf{r}^{B/S^*} $$thus also:
$$ {}^R\mathbf{v}^{S^*} = {}^R\mathbf{v}^A - {}^R\omega^S \times \mathbf{r}^{A/S^*} \\ {}^R\mathbf{v}^{S^*} = {}^R\mathbf{v}^B - {}^R\omega^S \times \mathbf{r}^{B/S^*} $$R_v_Ss__A = -S.ang_vel_in(R).cross(A.pos_from(Ss))
R_v_Ss__A
R_v_Ss__B = -S.ang_vel_in(R).cross(B.pos_from(Ss))
R_v_Ss__B
These are two expressions for the velocity of $S^*$ in R and thus they should equate. We can use that fact to solve for the measure numbers.
(R_v_Ss__B - R_v_Ss__A).to_matrix(C)
sol = sm.solve((R_v_Ss__A - R_v_Ss__B).to_matrix(C), ws2, ws3)
sol
Now substitute these measure numbers into the velocities we need.
S.ang_vel_in(R).subs(sol)
S.set_ang_vel(R, S.ang_vel_in(R).subs(sol))
S.ang_vel_in(R)
R_v_Ss__A.subs(sol)
Ss.set_vel(R, R_v_Ss__A.subs(sol))
Ss.vel(R)
$\hat{\mathbf{n}}$ is the vector normal to the cone's surface pointing from $S^*$ towards the cone's spin axis.
n = -sm.cos(theta)*C.x + sm.sin(theta)*C.y
Now the velocity of the point P on S can be found.
Ps.set_pos(Ss, r*n)
Ps.v2pt_theory(Ss, R, S)
The position to $P_c$ takes a bit of geometry to locate relative to O:
Pc.set_pos(O, (r+r*sm.sin(theta))*C.y + (b-r*sm.cos(theta))*C.x)
Now, both O and $P_c$ are fixed in C, so the two point theorem can be used.
O.set_vel(R, 0)
Pc.v2pt_theory(O, R, C)
Recall ${}^R\mathbf{v}^{P_c} = {}^R\mathbf{v}^{P_s}$. We can form the expression on the left side of: $\left({}^R\mathbf{v}^{P_c} - {}^R\mathbf{v}^{P_s}\right)\cdot\hat{\mathbf{c}}_z = 0$ (there is only a $\hat{\mathbf{c}}_z$ component).
(Ps.vel(R) - Pc.vel(R)).to_matrix(C)
Use this one equation to solve for the last measure number.
sol = sm.solve((Ps.vel(R) - Pc.vel(R)).to_matrix(C), ws1)
sol
Also recall that ${}^S\omega^C \cdot \hat{n} = 0$ and we can form the necessary angular velocity using the addition theorem:
S_w_C = C.ang_vel_in(R) - S.ang_vel_in(R)
S_w_C
Now, solve this expression for $b$:
zero = S_w_C.dot(n).subs(sol)
zero
b_sol = sm.solve(zero, b)[0]
b_sol
SymPy simplies this into a seemingly less simple solution. You can check that this is the correct solution by seeing if you solution minus the expected solution result in zero when simplified.
b_expected = r*(1+sm.sin(theta))/(sm.cos(theta)-sm.sin(theta))
b_expected
sm.trigsimp((b_sol - b_expected))