Question
Consider an axle which connects two wheels of radius $R$ each at distance $L$ from its center; each wheel may spin frictionlessly about the axle but they are constrained by strong static friction such that they roll with respect to the ground (See the figure below). The following equations are the Newton-Euler equations for each part of the system shown in the figure
\begin{align*} {\bf{X}}_{P_1}+{\bf{X}}_{C_1}+\,m{\bf{g}} &= m \ddot{{\bf{r}}}_{C_1} \\ {\bf{X}}_{P_2}+{\bf{X}}_{C_2}+\,m{\bf{g}} &= m \ddot{{\bf{r}}}_{C_2} \\ -{\bf{X}}_{C_1}-{\bf{X}}_{C_2}+M{\bf{g}} &= M \ddot{{\bf{r}}}_{C_3} \\ {\bf{T}}_{1}+{\bf{Y}}_{C_1}+{}_{C_1}{\bf{r}}_{P_1}\times{\bf{X}}_{P_1} &= {}_{C_1}\dot{\bf{H}}_{1} \\ {\bf{T}}_{2}+{\bf{Y}}_{C_2}+{}_{C_2}{\bf{r}}_{P_2}\times{\bf{X}}_{P_2} &= {}_{C_2}\dot{\bf{H}}_{2} \\ -{\bf{Y}}_{C_1}-{\bf{Y}}_{C_2}-{}_{C_3}{\bf{r}}_{C_1}\times{\bf{X}}_{C_1}-{}_{C_3}{\bf{r}}_{C_2}\times{\bf{X}}_{C_2} &= {}_{C_3}\dot{\bf{H}}_{3} \tag{1} \end{align*}
where we also know that
\begin{align*} {\bf{Y}}_{C_1} \boldsymbol{\cdot} {\bf{b}}_1 &= 0 \\ {\bf{Y}}_{C_2} \boldsymbol{\cdot} {\bf{b}}_1 &= 0 \tag{2} \end{align*}
which means that the connection of the wheels and the shaft is such that no constraint torque is applied in ${\bf{b}}_1$ direction.
Also, by assuming that both wheels roll with respect to the ground one can obtain the following constraint equations
\begin{align*} \dot{x} &= \frac{R}{2}(\dot{\varphi_1}+\dot{\varphi_2})\sin\theta \\ -\dot{y} &= \frac{R}{2}(\dot{\varphi_1}+\dot{\varphi_2})\cos\theta \\ \theta-\theta_0 &= \frac{R}{2L}(\varphi_1-\varphi_2) \tag{3} \end{align*}
My final aim is to find the generalized coordinates $\varphi_1,\varphi_2,\varphi_3,\theta,x,y$ as a function of time.
So I need another three scalar equations in terms of generalized coordinates and their derivatives other than the three constraint equations in $(3)$. I think that we should obtain them by manipulating equations in $(1)$; however, I am unable to obtain those. Any hint or help for finding them is appreciated. :)
Notice
- All of the quantities $\ddot{{\bf{r}}}_{C_1},\ddot{{\bf{r}}}_{C_2},\ddot{{\bf{r}}}_{C_3},\dot{\bf{H}}_{1},\dot{\bf{H}}_{2},\dot{\bf{H}}_{3},{}_{C_3}{\bf{r}}_{C_1},{}_{C_1}{\bf{r}}_{P_1}$ depend on the generalized coordinates $\varphi_1,\varphi_2,\varphi_3,\theta,x,y$ and their first and second time derivatives.
- The system is linear with respect to the constraint forces ${\bf{X}}_{P_1},{\bf{X}}_{P_2},{\bf{X}}_{C_1},{\bf{X}}_{C_2},{\bf{Y}}_{C_1},{\bf{Y}}_{C_2}$.
- ${}_{P_1}{\bf{r}}_{C_1}={}_{P_2}{\bf{r}}_{C_2}=R{\bf{a}}_3$ and $-{}_{C_3}{\bf{r}}_{C_1}={}_{C_3}{\bf{r}}_{C_2}=L{\bf{b}}_1$.
- That would be nice if we could obtain a vectorial equation without choosing a basis. So we may think of elimiating constraint forces from equations in $(1)$. As an attemp I summed the first three and second three equations in $(1)$ to obtain
\begin{align*} {\bf{X}}_{P_1}+{\bf{X}}_{P_2}+(M+2m){\bf{g}} &= (M+2m)\ddot{{\bf{r}}}_{C_3} \\ {\bf{T}}_{1}+{\bf{T}}_{2}+{}_{C_1}{\bf{r}}_{P_1}\times({\bf{X}}_{P_1}+{\bf{X}}_{P_2}) &= {}_{C_1}\dot{\bf{H}}_{1}+{}_{C_2}\dot{\bf{H}}_{2}+{}_{C_3}\dot{\bf{H}}_{3} \\ -{}_{C_3}{\bf{r}}_{C_1}\times({\bf{X}}_{C_1}-{\bf{X}}_{C_2}) & \tag{4} \end{align*} $\quad\,\,$ however, I don't know how to proceed further and make ${\bf{X}}_{C_1}-{\bf{X}}_{C_2}$ disappear.
Notation
- $m$ is the mass of each wheel and $M$ is the mass of the shaft. $R$ is the the radius of wheels and $2L$ is the length of the shaft.
- ${\bf{g}}$: accelaration due to gravity.
- ${\bf{X}}_{P_i}$: The constraint force applied at point $P_i$ which is exerted by surface to the wheel $i$.
- ${\bf{X}}_{C_i}$: The constraint force applied at point $C_i$ which is exerted by shaft to the wheel $i$.
- ${\bf{Y}}_{C_i}$: The constraint torque which is exerted by shaft to the wheel $i$.
- ${\bf{T}}_{i}$: The external torque applied by the actuator to the wheel $i$.
- ${\bf{r}}_{P}$: The position vector of point $P$ with respect to origin of frame $a$.
- ${}_{Q}{\bf{r}}_{P}$: The position vector of point $P$ with respect to piont $Q$.
- ${}_{P}{\bf{H}}_i$: The angular momentum of part $i$ with respect to point $P$.
- The overdot notation represents the time derivative with respect to the inertial frame $a$.

Hints
$1$. Combine the equations in $(4)$ to eliminate ${\bf{X}}_{P_1}+{\bf{X}}_{P_2}$.
$2$. Consider the result of step $1$ along with the third equation in $(1)$.
$3$. Take a look at this post and then solve for ${\bf{X}}_{C_1}$ and ${\bf{X}}_{C_2}$.
$4$. Think about the rest of the way. It is rather straight forward.