Modeling rigid body dynamics: Trouble with equations, encountering nonsense answers

Question

Hi, I am working on modeling a simple rigid body dynamic system through a certain range of motion. Here is the picture of the system.

I am trying to solve the system for the variables $x_1'', y_1'', x_2'',y_2''$, the acceleration variables for points A and B respectively, the angular acceleration which we will call \ddot{\theta} corresponding to the angle seen in the diagram, and the force in the rod which is given as Fab. F is the input force that is held constant.

I could have used a typical approach of basing coordinates and equations around the center of mass but I decided I wanted to try it this way. ALSO, the rod connecting the two masses is considered massless. The masses are given value 1, so is r the rod length for mathematical simplification.

Here are the general equations of motion used. SUM OF FORCES AT $A(x_1,y_1)$

$$x_1'' = -F + F_{AB}\sin(\theta) \tag{eq.1}$$

$$y_1'' = F_{AB}\cos(\theta) - N_{A} = 0 \tag{eq.2}$$

SUM OF FORCES AT B(x_2,y_2)

$$x_2'' = Nb - F_{AB}\sin(\theta) = 0\tag{eq.3}$$

$$y_2'' = -F_{AB}cos(\theta)\tag{eq.4}$$

So far we have 5 variables and 4 equations, we will proceed to the moment equation about B

$$Mb = R_{BA} \times N_{ab} - R_{BA} \times F = I_{A}\ddot{\theta}\tag{eq.5}$$

$$Mb = rN_{A}sin(\theta) - rFcos(\theta) = I_{A}\ddot{\theta} = m_{1}r^2\ddot{\theta}\tag{eq.6}$$

As stated earlier, r and m are equal to 1 for mathematical simplification, so we get

$$N_{A}sin(\theta) - Fcos(\theta) = \ddot{\theta}\tag{eq.7}$$

The last equation is the relative acceleration for point A

$$R_{A} = R_{B} + R_{BA}\tag{eq.8}$$

$$R_{A}'' = R_{B}'' + \ddot{\theta} \times R_{BA} - w^2R_{BA}\tag{eq.9}$$

$$R_{A}'' = [0\vec i - F_{AB}cos(\theta)\vec j] + [-\ddot{\theta}\vec k \times (-rsin(\theta)\vec i - rcos(\theta)\vec j)] - \dot{\theta}^2[-rsin(\theta)\vec i - rcos(\theta)\vec j]\tag{eq.10}$$

$$R_{A}'' = [0\vec i - F_{AB}cos(\theta)\vec j] + [-r\ddot{\theta}cos(\theta)\vec i + r\ddot{\theta}sin(\theta)\vec j] + [r\dot{\theta}^2sin(\theta)\vec i + r\dot{\theta}^2cos(\theta)\vec j]\tag{eq.11}$$

Setting r = 1 as given from before and organizing terms under the correct components we get

$$R_{A}'' = [-\ddot{\theta}cos(\theta) + \dot{\theta}^2sin(\theta)]\vec i + [\dot{\theta}^2cos(\theta) + \ddot{\theta}sin(\theta) - F_{AB}cos(\theta)]\vec j\tag{eq.12}$$

Comparing these components to the components for point A

$$x_1'' = -\ddot{\theta}cos(\theta) + \dot{\theta}^2sin(\theta)\tag{eq.13}$$

$$y_1'' = 0 = \dot{\theta}^2cos(\theta) + \ddot{\theta}sin(\theta) - F_{AB}cos(\theta)\tag{eq.14}$$

Solving for $\dot{\theta}^2$ in the second equation we have

$$\dot{\theta}^2 = -\ddot{\theta}tan(\theta) + F_{AB}\tag{eq.15}$$

Plugging this into the new equation for $x_1''$ we get

$$x_1'' = -\ddot{\theta}cos(\theta) + F_{AB}sin(\theta) - \ddot{\theta}tan(\theta)sin(\theta)\tag{eq.15}$$

$$x_1'' = F_{AB}sin(\theta) + \ddot{\theta}[-tan(\theta)sin(\theta) - cos(\theta)]\tag{eq.16}$$

At this point, I have enough equations to solve for my 6 variables($x_1'',y_1'',x_2'',y_2'',F_{ab},\ddot{\theta}$). I put these into a matrix form and got a nonsense answer, I was hoping to get some help. I understand I can formulate the problem in a different way around the center of mass but tvecs not what I was aiming to do. I want to formulate it with these coordinates. I feel like I got a sign wrong somewhere but I'm stuck trying to find it. Edit: format edited with matrices added

The matrix of coefficients is A $$ \begin{pmatrix} 1 & 0 & 0 & 0 & -sin(\theta) & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & cos(\theta) & 0 \\ 0 & 0 & 0 & 0 & 0.5sin(2\theta) & -1 \\ 1 & 0 & 0 & 0 & -sin(\theta) & tan(\theta)sin(\theta) + cos(\theta) \\ \end{pmatrix} $$

Where X is

$$ \begin{pmatrix} x_1'' \\ y_1'' \\ x_2'' \\ y_2'' \\ F_{AB} \\ \ddot{\theta} \\ \end{pmatrix} $$

And the output matrix B is given with the input force F = 1

$$ \begin{pmatrix} -1 \\ 0 \\ 0 \\ 0 \\ cos(\theta) \\ 0 \\ \end{pmatrix} $$

Answer 1

I usually prefer to derive the equations of motion using Lagrangian approach, or in this case its generalization D'Alembert's principle, because of the external force $F$, which may or may not be conservative. But let's go with the force description approach to both endpoints of the segment.

I think the equation of angular velocity / torque (8) is not correct. I strongly advise you to calculate all torques with respect to the center of mass of the two mass-points. Moreover, that equation might be redundant in general as it is derived from the force equations. So you can get the equation of motion and the various force components without it. Observe that this is a two mass-point system connected by a massless rod, so it is interpreted as two mass-pints bound by holonomic constraints and interpretation as a solid body dynamics is not needed.

Furthermore, equations (8) and (9) seem to be redundant and incorrectly expanded. what I mean is that by definition $\vec{R}_{BA} = \vec{R}_B - \vec{R}_A$, so I am not sure one gets anything new by reparametrizing this in terms of $\theta$. In addition to that force enters these equations where there should be acceleration. Equation (14) is also incorrect, leading to wrong consequent equations. You can see that if you carefully consider the full set of holonomic constraints of this system. Notice that $y_a = y_1 = 0$ which is not parametrized in terms of $\theta$ and therefore, you cannot write equation (14). Moreover, it can be a bit tricky (if possible at all) to eliminate the first derivative $\dot{\theta}$ in a system without some kind of conservation law, like for example conservation of energy, which you do not have if $F$ is a general non-conservative force.

Here is my derivation. I have included the torque equation, just to keep track of it.

\begin{align} &\text{Equations of motion:}\\ &m_a \, \frac{d^2}{dt^2} \big(x_a \,\vec{i} + y_a \,\vec{j}\big) \, = \, \vec{F}_{ab} \, + \, F \,\vec{i} \, - \, N_a \, \vec{j}\\ &\\ &m_b \, \frac{d^2}{dt^2} \big(x_b \,\vec{i} + y_b \,\vec{j}\big) \, = \, - \, \vec{F}_{ab} \, - \, N_b \, \vec{i}\\ &\\ &I \, \frac{d^2\theta}{dt^2} \,\vec{k} \, = \, - \, \frac{m_b}{m_a+m_b} \vec{AB} \times ({F}\, \vec{i} - N_a \, \vec{j}) \, +\, \frac{m_a}{m_a+m_b} \vec{AB} \times ( - N_b \, \vec{i})\\ &\\ &\text{Holonomic constraints: }\\ & y_a = 0\\ & x_b = 0\\ & x_a^2 + y_b^2 = r^2 \end{align}

Observe that since $\vec{F}_{ab}$ is aligned with $\vec{AB}$, one can infer that $\vec{AB} \times \vec{F}_{ab} = \vec{0}$, so that is why this force is not present in the third equations, the angular velocity / torque equation.

Write $$\vec{F}_{ab} \, = \, - \, F_{ab} \sin(\theta)\, \vec{i} \, + \, F_{ab} \cos(\theta)\, \vec{j}$$ $$\vec{AB} \, = \, - \, r \sin(\theta)\, \vec{i} \, + \, r \cos(\theta)\, \vec{j}$$

For simplicity, set $r_a = \frac{m_b \, r}{m_a + m_b}$ and $r_b = \frac{m_a \, r}{m_a + m_b}$

\begin{align} &\text{Equations of motion:}\\ &m_a \, \frac{d^2}{dt^2} \big(x_a \,\vec{i} + y_a \,\vec{j}\big) \, = \, \big(\, F - F_{ab}\sin(\theta) \,\big)\,\vec{i} \, + \, \big(\, F_{ab}\cos(\theta) - \, N_a\,\big) \, \vec{j}\\ &\\ &m_b \, \frac{d^2}{dt^2} \big(x_b \,\vec{i} + y_b \,\vec{j}\big) \, = \, \big(\, F_{ab}\sin(\theta) - \, N_b\,\big) \, \vec{i} \, - \, F_{ab}\cos(\theta) \, \vec{j}\\ &\\ &I \, \frac{d^2\theta}{dt^2} \,\vec{k} \, = \, r_a \big(\, \sin(\theta)\, \vec{i} - \cos(\theta)\, \vec{j} \,\big) \times ({F}\, \vec{i} - N_a \, \vec{j}) \, +\, r_b \big(\, \sin(\theta)\, \vec{i} - \cos(\theta)\, \vec{j} \,\big) \times ( N_b \, \vec{i})\\ &\\ &\text{Holonomic constraints: }\\ & y_a = 0\\ & x_b = 0\\ & x_a^2 + y_b^2 = r^2 \end{align}

Simplify the cross-products and parametrize the holonomic constraints \begin{align} &\text{Equations of motion:}\\ &m_a \, \frac{d^2}{dt^2} \big(x_a \,\vec{i} + y_a \,\vec{j}\big) \, = \, \big(\, F - F_{ab}\sin(\theta) \,\big)\,\vec{i} \, + \, \big(\, F_{ab}\cos(\theta) - \, N_a\,\big) \, \vec{j}\\ &\\ &m_b \, \frac{d^2}{dt^2} \big(x_b \,\vec{i} + y_b \,\vec{j}\big) \, = \, \big(\, F_{ab}\sin(\theta) - \, N_b\,\big) \, \vec{i} \, - \, F_{ab}\cos(\theta) \, \vec{j}\\ &\\ &I \, \frac{d^2\theta}{dt^2} \,\vec{k} \, = \, \Big( \, r_a \, F\cos(\theta) - r_a\,N_a \sin(\theta) + r_b \, N_b \cos(\theta)\,\Big)\, \vec{k}\\ &\\ &\text{Holonomic constraints: }\\ & x_a = r\sin(\theta)\\ & y_a = 0\\ & x_b = 0\\ & y_b = r\cos(\theta) \end{align}

Lets perform some of the differentiations \begin{align} &\text{Holonomic constraints: }\\ & x_a = r\sin(\theta)\\ & y_a = 0\\ & x_b = 0\\ & y_b = r\cos(\theta) \\ & \frac{d^2x_a}{dt^2} = r \cos(\theta) \frac{d^2\theta}{dt^2} - r \sin(\theta)\left(\frac{d\theta}{dt}\right)^2\\ & \frac{d^2 y_a}{dt^2} = 0\\ & \frac{d^2 x_b}{dt^2} = 0\\ & \frac{d^2y_b}{dt^2} = -\, r \sin(\theta) \frac{d^2\theta}{dt^2} - r \cos(\theta)\left(\frac{d\theta}{dt}\right)^2 &\\ &\text{Equations of motion:}\\ &m_a \, \frac{d^2x_a}{dt^2} \,\vec{i} \, = \, \big(\, F - F_{ab}\sin(\theta) \,\big)\,\vec{i} \, + \, \big(\, F_{ab}\cos(\theta) - \, N_a\,\big) \, \vec{j}\\ &\\ &m_b \, \frac{d^2y_b}{dt^2} \, \vec{j} \, = \, \big(\, F_{ab}\sin(\theta) - \, N_b\,\big) \, \vec{i} \, - \, F_{ab}\cos(\theta) \, \vec{j}\\ &\\ &I \, \frac{d^2\theta}{dt^2} \,\vec{k} \, = \, \Big( \, r_a \, F\cos(\theta) - r_a\,N_a \sin(\theta) + r_b \, N_b \cos(\theta)\,\Big)\, \vec{k} \end{align}

You can express $N_a = F_{ab}\cos(\theta) \,$ and $\, N_b = F_{ab} \sin(\theta)$ and write

\begin{align} &m_a r \cos(\theta) \frac{d^2\theta}{dt^2} - m_ar \sin(\theta)\left(\frac{d\theta}{dt}\right)^2 \, = \, F - F_{ab}\sin(\theta)\\ &m_b r \sin(\theta) \frac{d^2\theta}{dt^2} + m_br \cos(\theta)\left(\frac{d\theta}{dt}\right)^2 \, = \, F_{ab}\,\cos(\theta) \\ &I \, \frac{d^2\theta}{dt^2} \, = \, r_a \, F\cos(\theta) + (r_b-r_a)\,F_{ab}\, \cos(\theta)\sin(\theta) \end{align} Multiply both sides of the first equation by $-\sin(\theta)$, both sides of the first equation by $\cos(\theta)$ and add them. Consequently, one can express $$F_{ab} \, = \, (m_br - m_ar)\cos(\theta)\sin(\theta) \frac{d^2\theta}{dt^2} \, + \, (m_a+m_b)r\left(\frac{d\theta}{dt}\right)^2 \, - \, F\sin(\theta)$$ Multiply both sides of the first equation by $\cos(\theta)$, both sides of the first equation by $\sin(\theta)$ and add them. In this case, you obtain the equation of motion, expressed in $\theta$ $$\Big(\,m_a r \cos^2(\theta) + m_b r \sin^2(\theta)\,\Big) \frac{d^2\theta}{dt^2} \, + \, r(m_b - m_a) \cos(\theta) \sin(\theta) \left(\frac{d\theta}{dt}\right)^2 \, =\, F \cos(\theta)$$

Answer 2

This post may not answer your question explicitly, but perhaps you (or others) can use the results to validate your own methods.

Assuming I don't misinterpret the drawing, you can depict the mechanics as shown

Mass 1 is constrained to move along the $x$-axis and mass 2 is constrained to move along the $y$-axis.
Also, both masses are constrained from rotation which leaves only two degrees of freedom: $x_1$ and $y_2$. However, the rod imposes a constraint between those coordinates, as $x_1=-r\sin(\theta)$ and $y_2=r\cos(\theta)$.
You can therefore describe the system using only one coordinate $\theta$.

The Lagrangian $\mathscr{L}=K-P$ is the difference between kinetic energy $$K=\frac{m_1\dot{x}_1^2}{2} + \frac{m_2\dot{y}_2^2}{2}=\frac{r^2\dot{\theta}^2}{2}\left(m_1\cos^2(\theta) + m_2\sin^2(\theta)\right)$$ and potential energy $$P = m_2g \ y_2=m_2g \ r\cos(\theta)$$ and you derive the equations of motion using the Lagrangian equation: $$\frac{d}{dt}\left(\frac{\partial\mathscr{L}}{\partial\dot{\theta}}\right)-\frac{\partial\mathscr{L}}{\partial\theta}=Q=-Fr\cos(\theta)$$ Working out the details will give you the equation $$\ddot{\theta}=\frac{r\dot{\theta}^2\sin(\theta)\cos(\theta)(m_1 -m_2 ) + m_2g \sin(\theta)-F\cos(\theta)}{r(m_1\cos^2(\theta)+m_2\sin^2(\theta))}$$

Typically when dealing with nonlinear systems, you obtain a numerical solution using an ODE solver to get a sequence of state vectors $\mathbf{x}(t)=\begin{bmatrix}\theta(t)\\\dot{\theta}(t)\end{bmatrix}$corresponding to a time interval. Then you can find the motion of the masses from $\mathbf{x}$.