Motion of a rigid body in Abaqus/Standard

Force-free motion of a rigid body

Problem description

The problem is shown in Figure 1. An arbitrary symmetric body whose rotary inertia about its axis of symmetry is different from its value along the two other principal axes spins around its axis of symmetry with an initial angular velocity $ω$ The body is modeled with a ROTARYI element whose second moments of inertia along its principal axes, $e_{α}$ ( $α =$ 1, 2, 3), have the values $I_{11} = I_{22}$ and $I_{33} = 2 I_{11}$ . The axis of symmetry is $e_{3}$ . Dummy nodes are attached rigidly to the ROTARYI element along the principal axes by using a BEAMMPC so that their displacements can be tracked. Since ROTARYI elements have only rotational degrees of freedom, a MASS element is needed on top of the ROTARYI element to activate translational degrees of freedom at these dummy nodes. Initial conditions are taken from the analytical solution presented below.

For the force-free symmetric body $N_{1} = N_{2} = N_{3} = 0$ and $I_{11} = I_{22}$ ; therefore, the Euler's equations reduce to

\begin{matrix} I_{11} {\ddot{ϕ}}^{1} + (I_{33} - I_{11}) {\dot{ϕ}}^{2} {\dot{ϕ}}^{3} = 0 \\ I_{11} {\ddot{ϕ}}^{2} + (I_{11} - I_{33}) {\dot{ϕ}}^{3} {\dot{ϕ}}^{1} = 0 \\ I_{33} {\ddot{ϕ}}^{3} = 0. \end{matrix}

The last equation can be integrated to give ${\dot{ϕ}}^{3} = ω$ , where $ω$ is a constant, defined as an initial condition of the problem.

To determine ${\dot{ϕ}}^{1}$ , we take the time derivative of the first equation:

I_{11} \frac{d {\ddot{ϕ}}^{1}}{d t} + (I_{33} - I_{11}) {\ddot{ϕ}}^{2} ω = 0.

Using the second equation to solve for ${\ddot{ϕ}}^{2}$ gives

\frac{d {\ddot{ϕ}}^{1}}{d t} + {[\frac{(I_{33} - I_{11}) ω}{I_{11}}]}^{2} {\dot{ϕ}}^{1} = 0.

Similarly,

\frac{d {\ddot{ϕ}}^{2}}{d t} + {[\frac{(I_{11} - I_{33}) ω}{I_{11}}]}^{2} {\dot{ϕ}}^{2} = 0.

These equations describe simple harmonic motion with angular frequency

Ω = \frac{(I_{33} - I_{11})}{I_{11}} ω .

With appropriate initial conditions the solution is ${\dot{ϕ}}^{1} = - A \sin Ω t$ , where A is a constant. The corresponding solution for ${\dot{ϕ}}^{2}$ can be found by substituting this solution for ${\dot{ϕ}}^{1}$ into the first of the Euler equations, giving ${\dot{ϕ}}^{2} = A \cos Ω t$ . The corresponding initial conditions are ${{\dot{ϕ}}^{1} |}_{t = 0} = 0$ , ${{\dot{ϕ}}^{2} |}_{t = 0} = A$ , ${{\dot{ϕ}}^{3} |}_{t = 0} = ω$ . We also choose ${ϕ^{1} |}_{t = 0} = - A / Ω$ , ${ϕ^{2} |}_{t = 0} = 0$ , and ${ϕ^{3} |}_{t = 0} = 0$ . These initial rotation conditions give rise to the local orientation indicated in Figure 1. The directions of the principal axes of inertia of the body are defined. We choose $A =$ 0.25, $I_{11} = I_{22} =$ 1, $I_{33} = 2 I_{11} =$ 2, so that

Ω = \frac{(I_{33} - I_{11})}{I_{11}} ω = ω = 2 π .

Initial angular velocities, ${\dot{ϕ} |}_{t = 0}$ , must be applied to node 1, and translational velocities, ${v |}_{t = 0}$ , must be applied to the dummy nodes lying along the legs of the axes of the body. The translational velocity components are obtained from

v = \dot{ϕ} \times r,

where $r$ is the vector connecting the center of the body (node 1) to one of the nodes along the principal axes (node 2, 3, or 4). This latter initial velocity calculation is performed internally by Abaqus for each dummy node as a result of applying the BEAMMPCs mentioned previously. The model is shown in rigidbodymotion_free.inp.

The dynamic response of the body subjected to the above initial conditions is tracked for two seconds. Large-rotation theory is used, so the principal axes of inertia rotate with the rotation of the ROTARYI element. Rigid body rotary inertia contributes nonsymmetric terms to the system matrix when the motion is in three dimensions. Therefore, we use the unsymmetric equation solver. Numerical damping is removed from the implicit dynamic operator.

Results and discussion

The harmonic response for the angular velocity relative to the global coordinate system is obtained in the Abaqus solution and is plotted in Figure 2, Figure 3, and Figure 4. These angular velocity values are obtained from node 1. Noting that ${{\dot{ϕ}}^{z} |}_{t = 0} = {{\dot{ϕ}}^{3} |}_{t = 0} \cos (A / Ω) - {{\dot{ϕ}}^{2} |}_{t = 0} \sin (A / Ω)$ , ${{\dot{ϕ}}^{z} |}_{t = 0}$ can be calculated as 6.268. This is shown accurately in Figure 4.

The solutions for ${\dot{ϕ}}^{1}$ and ${\dot{ϕ}}^{2}$ obtained above indicate that the vector ${\dot{ϕ}}^{1} e_{1}$ + ${\dot{ϕ}}^{2} e_{2}$ is of constant magnitude and precesses about the body 3-axis with the angular frequency $Ω = 2 π$ . The evolution of this vector with respect to the global coordinate system is plotted in Figure 5 as an X–Y plot of the history of ${\dot{ϕ}}^{x}$ versus the history of ${\dot{ϕ}}^{y}$ for node 1. As expected, the result traces a circle of diameter A. Figure 6 shows a similar plot of $u^{x}$ versus $u^{y}$ for node 4, viewed by looking down the global z-axis.

The precession described by Goldstein is relative to the body axes, which are themselves rotating in space at a frequency of $2 π$ . In large-displacement analysis in Abaqus (geometric nonlinearities considered in the step) the principal axes of inertia rotate with the rotation of the node to which the ROTARYI element is attached. This explains why the period of the motion observed in the figures is 0.5 and not 1.0.

The analysis is completed in 200 increments, with each increment requiring only 1 iteration to satisfy the moment equilibrium criterion.

Input files

rigidbodymotion_free.inp: Implicit force-free motion analysis.

Figures

Figure 1. Rigid body rotation example.

Figure 2. Angular velocity response (node 1).

Figure 3. Angular velocity response (node 1).

Figure 4. Angular velocity response (node 1).

Figure 5. Precession of angular velocity (node 1).

Figure 6. Precession of rotation (node 4).

Forced motion of a rigid body

Problem description

The top is modeled with a ROTARYI element, and a local coordinate system is used to prescribe the second moments of inertia along the principal axes $(e_{1}, e_{2}, e_{3})$ . A 2-node rigid beam element RB3D2 is used to connect the fixed point of the top, O, with its center of mass, C. The effect of the gravitational field is considered by applying a concentrated load of magnitude $M g$ in the z-direction at point C. The initial conditions for the angular velocity, $\dot{ϕ}$ , are prescribed in the global system of coordinates. For postprocessing and visualization purposes only, a second RB3D2 element is added at point C in a direction perpendicular to the axis of rotation. The Abaqus solution is compared to the analytical solution, which is outlined in the next section.

The problem is also solved using connector elements. A CONN3D2 element of type BEAM is used to model the top. A CONN3D2 element of type EULER is used to obtain the Euler angles.

Analytical solution

The solution for the motion of the symmetric top is described in Goldstein (1980), Whittaker (1988), and Macmillan (1936).

The analytical solution is described in terms of the Euler angles: $θ, γ, ψ$ , where $θ$ measures the inclination of the $e_{3}$ -axis from the vertical, $γ$ measures the azimuth of the top about the vertical, and $ψ$ is the rotation angle of the top around its own $e_{3}$ -axis. Since the system is conservative, the total energy is constant in time. By denoting $ω_{1} = {\dot{ϕ}}^{1}, ω_{2} = {\dot{ϕ}}^{2}$ , and $ω_{3} = {\dot{ϕ}}^{3}$ , the energy conservation equation gives

E = K E + V = \frac{1}{2} (ω_{1}^{2} + ω_{2}^{2}) I_{11} + \frac{1}{2} I_{33} ω_{3}^{2} + M g l \cos θ = c o n s t a n t,

where $I_{11} = I_{22}$ since the body is symmetrical.

The energy equation can be arranged in the following form:

{\dot{u}}^{2} = (1 - u^{2}) (α - β u) - {(b - a u)}^{2},

where $u = \cos θ$ and the constants $a, b, α, β$ have the following form:

a = \frac{I_{33}}{I_{11}} ω_{3}; b = \frac{K}{I_{11}}; α = \frac{2}{I_{11}} (E - \frac{I_{33}}{2} ω_{3}^{2}); β = \frac{2 M g l}{I_{11}} .

In these relations K is the moment of momentum with respect to the z-axis. Its value is constant in time and is given by

K = I_{11} ω_{1} η_{1} + I_{22} ω_{2} η_{2} + I_{33} ω_{3} η_{3},

where in terms of Eulerian angles the direction cosines are

η_{1} = \sin θ \sin ψ; η_{2} = \sin θ \cos ψ; η_{3} = \cos θ .

The equation of motion for $u = \cos θ$ is an elliptic function of time, and the integration is not straightforward since the function presents singularities.

We can arrange this equation in the following form ${\dot{u}}^{2} = β f (u)$ . The function $f (u)$ has two real roots $u_{3} = \cos θ_{3}$ and $u_{2} = \cos θ_{2}$ situated between $u =$ −1 and $u =$ +1. The third root $u_{1} > u_{2} > u_{3}$ is greater than +1. The top will move such that $θ$ always remains between the roots $θ_{2}$ and $θ_{3}$ , which are called turning angles.

The equation of motion for $θ$ can be expressed in terms of these three roots as follows:

{\dot{u}}^{2} = β (u - u_{3}) (u - u_{2}) (u - u_{1}) .

By expressing the constants of integration in terms of the three roots, one can obtain the analytical solution of this equation by reducing the elliptic integral to a normal form. This solution is given in Macmillan (1936):

u = u_{3} + 4 h^{1 / 2} \sqrt{(u_{1} - u_{3}) (u_{2} - u_{3})} {[\frac{\sin ξ - h^{2} \sin 3 ξ + h^{6} \sin 5 ξ + \dots}{1 - 2 h \cos 2 ξ + 2 h^{4} \cos 4 ξ - 2 h^{9} \cos 6 ξ + \dots}]}^{2},

where

h = e^{(- \frac{Ω_{3}}{Ω_{1}} π)}; ξ = \frac{π}{2 ω_{1}} (\frac{\sqrt{β}}{2} t + Ω_{3}) .

In the above equation $Ω_{1}$ and $Ω_{3}$ are elliptic integrals of the first kind and have the following expressions:

Ω_{1} = \frac{1}{\sqrt{(u_{1} - u_{3})}} \int_{0}^{π / 2} \frac{d ϕ}{\sqrt{1 - k^{2} \sin^{2} ϕ}},

Ω_{3} = \frac{1}{\sqrt{(u_{1} - u_{3})}} \int_{0}^{π / 2} \frac{d ϕ}{\sqrt{1 - k_{1}^{2} \sin^{2} ϕ}},

where

k^{2} = \frac{u_{2} - u_{3}}{u_{1} - u_{3}}, k_{1}^{2} = \frac{u_{1} - u_{2}}{u_{1} - u_{3}} .

The values of the elliptic integrals are usually tabulated in calculus books or in mathematical tables. As soon as the roots of the polynomial $f (u)$ are found, we know the solution for the equation of motion.

After determining $u = \cos θ$ from the above equation, the remaining Euler's angles, $γ$ and $ψ$ , can be found from

\dot{γ} = \frac{b - a \cos θ}{\sin^{2} θ}; \dot{ψ} = \frac{I_{11} a}{I_{33}} - \cos θ \frac{b - a \cos θ}{\sin^{2} θ} .

The coordinates of the center of mass of the top in the x–y plane can be obtained if the first two Euler's angles $θ$ and $γ$ are known: $x = l \sin θ \cos γ$ and $y = l \sin θ \sin γ$ .

Results and discussion

In this section we will present the comparative results between Abaqus and the analytical solution for two situations often discussed in the literature. Many different response characteristics are possible depending on the initial conditions and inertia properties.

Case 1

Let us consider first that the symmetric top is spinning about its own axis $e_{3}$ , which is fixed in some direction $θ_{0} =$ 20°. At time $t = 0$ the symmetry or figure axis is released, and the top rotates around the $e_{3}$ -axis with angular velocity ${\dot{ϕ}}^{3} =$ 50. In addition to the angular velocity around the symmetry axis, we prescribe an angular velocity ${\dot{ϕ}}^{1} =$ 0.5 around the $e_{1}$ - or $e_{x}$ -axis. Usually the motion of the top is depicted by tracing the curve of the intersection of the $e_{3}$ -axis on a sphere of unit radius. This curve is called the “locus” of the figure axis. In our representation we will trace the projection of the locus in the x–y plane. According to the analytical solution, the ratio $u^{'} = \frac{b}{a}$ lies between the roots $u_{2}$ and $u_{3}$ , and the locus of the top axis exhibits loops (Goldstein, 1980).

We have chosen the length of the top axis $l =$ 1 and $M g =$ 20. The initial velocities in Abaqus are prescribed in the global coordinate system; therefore, the two components of the angular velocities $ω_{y} =$ −17.101 and $ω_{z} =$ 46.9846 in the global system will create a resultant angular velocity $\dot{ϕ^{3}} =$ 50 in the local system (Figure 7). The initial velocity in the global x-direction is the same as the initial velocity in the local $e_{1}$ -direction. The turning angles, obtained by solving the equation $f (u) = 0$ , are $u_{3} =$ 0.9517 and $u_{2} =$ 0.9112 or $θ_{3} =$ 24.32° and $θ_{2} =$ 17.88°, respectively. Based on the fact that the first Euler angle, $θ$ , is equal to the spherical angle used in the polar representation, the variation of this angle in time is obtained in Abaqus from the displacements. The turning angles are reproduced accurately in Abaqus, and the analytical solution is in good agreement with the Abaqus solution. This comparison is shown in Figure 8.

The numerical damping coefficient ALPHA was taken equal to zero in the direct integration scheme used in Abaqus. It is worth mentioning that the analytical solution is an approximate solution since the accuracy of this solution will depend on the number of terms taken in the expansion series and on the accuracy with which the elliptic integrals are evaluated. The projection on the x–y plane of the top's locus is depicted in Figure 9, where the analytical solution and the Abaqus solution are shown. The locus exhibits loops along with precession in the counterclockwise direction. The Abaqus solution agrees with the analytical solution; however, the analytical solution is extremely sensitive to the values of the elliptic integrals taken from the tables.

The averaged precession frequency prediction can be found from the analytical solution for a fast top; that is, a top that has a large initial kinetic energy compared to the maximum change in the potential energy. The theoretical averaged precession frequency is

f r e q = \frac{M g l}{I_{33}} ω_{3} = 0.4 r a d / s .

The total time for the complete precession in the x–y plane is 15 s, and the precession frequency given by Abaqus is, therefore, $f r e q = 2 π / 15 = 0.418 r a d / s$ .

The change in the potential energy is reflected in the external work; due to the small applied force and small displacements, the exernal work has small values. Therefore, the total energy is approximately equal to the kinetic energy of the system. The total energy and the external work obtained in Abaqus are presented in Figure 10 and Figure 11, respectively. For better visualization, the time variation of the external work is shown in Figure 11 only for the first 3 s of the spinning process.

Case 2

A second case assumes that the top is spinning only about its own axis. For this case the ratio $u^{'} = \frac{b}{a}$ coincides with one of the roots of the polynomial $f (u)$ , and the locus of the top axis exhibits cusps touching circles (Goldstein, 1980). In this case we prescribe only the angular velocity around the $e_{3}$ -axis, $\dot{ϕ^{1}} =$ 50. All of the other parameters are kept the same as before. The “turning” angles are obtained by solving again the equation $f (u) = 0$ with the new coefficients and are found to be $θ_{3} =$ 21.76° and $θ_{2} =$ 20°, respectively. The variation in time of the first Euler angle, $θ$ , is presented in Figure 12 the first 3 s of the process. The projection of the top's locus on the x–y plane, obtained in Abaqus, is presented in the Figure 13 where the analytical solution is also shown. The total energy and external work done for this case are presented in Figure 14 and Figure 15.

Abaqus/Explicit is also used to study the forced motion of the rigid top presented in this section. Due to the explicit time integration, the running time is less in Abaqus/Explicit. The top is modeled using a rigid R3D4 element and a ROTARYI element. The rigid body reference node is identical to the node of the ROTARYI element.

The problem is also solved in Abaqus/Explicit and Abaqus/Standard using connector elements. The Euler angles are obtained directly (in radians) as output variable CPR. The solution obtained using connector elements agrees well with the analytical solution.

Input files

rigidbodymotion_forced_std.inp: Implicit forced motion analysis.
rigidbodymotion_verify.f: Code used to generate the analytical solution.
rigidbodymotion_forced_xpl.inp: Forced motion analysis with Abaqus/Explicit.
rigidbodymotion_conn_f_std.inp: Forced motion analysis in Abaqus/Standard, using connector elements.
rigidbodymotion_conn_f_xpl.inp: Forced motion analysis in Abaqus/Explicit, using connector elements.

References

Goldstein, H., Classical Mechanics, Second Edition, Addison-Wesley, 1980.
Fowles, G. R., Analytical Mechanics, Third Edition, Holt, Rinehart and Winston, 1977.
MacMillan, W. D., Dynamics of Rigid Bodies, First Edition, McGraw-Hill Book, 1936.