Motion of a rigid body in Abaqus/Standard

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.