Equivalent rigid body dynamic motion

ProductsAbaqus/Standard

It is often useful to obtain the equivalent rigid body motion of part of a model (or of the whole model): the position and translational velocity of the part's center of mass and its angular rotation and velocity about the same center of mass. Abaqus/Standard provides such output, based on equivalent momentum.

Let V be the volume of a part for which the equivalent rigid body motion values are requested. The density of the part in its initial configuration is $ρ_{0} (S_{i})$ , where $S_{i}$ , $i = 1, 2, 3$ , are material coordinates in the part. The spatial position of a material particle in its initial configuration is $X (S_{i})$ and in the current configuration is $x (S_{i})$ , resulting in a displacement of $u (S_{i})$ . We wish to compute the current spatial position of the center of mass of the part, $x_{\circ}$ ; the translational velocity of an equivalent rigid body, ${\dot{x}}_{\circ}$ ; the angular velocity of this equivalent rigid body, $ω_{\circ}$ ; the translational displacement of an equivalent rigid body motion, $u_{\circ}$ ; and the rotation of an equivalent rigid body motion around the center of mass, $ϕ$ . In these definitions an “equivalent rigid body” means a rigid body with the same mass distribution and the same translational and angular momentum as the actual deforming part in the current configuration.

For simplicity of notation we define some quantities. The mass of the part is

M \overset{def}{=} \int_{V^{0}} ρ_{0} d V^{0} .

Its first mass moment about the origin is

J_{\circ} \overset{def}{=} \int_{V^{0}} ρ_{0} x d V^{0} .

Its second mass moment about the origin is

I_{\circ} \overset{def}{=} \int_{V^{0}} ρ_{0} (x \cdot x I - x x) d V^{0},

where $I$ is a unit matrix.

It is convenient to invoke the relation $x = N^{N} x^{N}$ (the summation convention is assumed), where $N^{N}$ are the finite element interpolation functions associated with each degree of freedom and $x^{N}$ is the vector of current nodal positions. We can now write

I_{\circ} \overset{def}{=} \int_{V^{0}} ρ_{0} (N^{N} N^{M} x^{N} \cdot x^{M} I - N^{N} N^{M} x^{N} x^{M}) d V^{0} .

Recognizing that the primitive mass matrix is

M^{N M} \overset{def}{=} \int_{V^{0}} ρ_{0} N^{N} N^{M} d V^{0},

we have

I_{\circ} \overset{def}{=} M^{N M} (x^{N} \cdot x^{M} I - x^{N} x^{M}) .

We can immediately obtain

x_{\circ} = J_{\circ} / M

and, by equating the translational momentum of the equivalent and the actual body,

{\dot{x}}_{\circ} = \int_{V^{0}} ρ_{0} \dot{x} d V^{0} / M .

The angular velocity of the part is defined by equating the angular momentum of the part and of the equivalent rigid body about the center of mass:

\int_{V^{0}} ρ_{0} (x - x_{\circ}) \times \dot{x} d V^{0} = I \cdot ω_{\circ},

where

I = I_{\circ} + M (x_{\circ} x_{\circ} - x_{\circ} \cdot x_{\circ} I)

is the second mass moment of the part about its center of mass.

Abaqus uses the lumped mass formulation for low-order elements. As a consequence, the second mass moments of inertia can deviate from the theoretical values, especially for coarse meshes. Certain Abaqus elements produce lumped or structural contributions to this second mass moment (rotary inertias) not shown in these equations.

This provides

ω_{\circ} = I^{- 1} \cdot [H_{\circ} - J_{\circ} \times {\dot{x}}_{\circ}],

where

H_{\circ} \overset{def}{=} \int_{V^{0}} ρ_{0} x \times \dot{x} d V^{0} = M^{N M} x^{N} \times {\dot{x}}^{M}

is the angular momentum of the part about the origin.

The perceived translational motion of the center of mass in an equivalent rigid body motion is calculated as

u_{\circ} = \int_{V^{0}} ρ_{0} u d V^{0} / M .

The equivalent rigid body rotation of the part with respect to its center of mass requires some conceptual approximations as follows. Denote the relative positions of a material particle with respect to the center of mass in the undeformed configuration and in the deformed configuration $\bar{X}$ and $\bar{x}$ , respectively. Consider that the configurations are known and that the axis of rotation of the body is denoted by the unit vector $p$ . A material particle sees such rotation relative to the center of mass as

ϕ = ϕ p = \frac{{\bar{X}}_{p} \times {\bar{x}}_{p}}{| {\bar{X}}_{p} \times {\bar{x}}_{p} |} \tan^{- 1} (\frac{| {\bar{X}}_{p} \times {\bar{x}}_{p} |}{{\bar{X}}_{p} \cdot {\bar{x}}_{p}}),

where the subscript p denotes the projection of a vector into a plane normal to $p$ . We now generalize this concept by integration of the constituent parts. Define

a = \int_{V} ρ {\bar{X}}_{p} \times {\bar{x}}_{p} d V, b = \int_{V} ρ {\bar{X}}_{p} \cdot {\bar{x}}_{p} d V .

The average Euler rotation then follows with the equation

ϕ_{a v} = \frac{a}{| a |} \tan^{- 1} (\frac{| a |}{b}) .

These integrals are not easily calculated, but with the modifications below they can be expanded in such a way that the (initially unknown) current center of mass, $x_{\circ}$ , only appears in products with the position of particles in the known current configuration.

A necessary condition for the validity of the intuitive generalization above is that if the part undergoes an arbitrary rigid body rotation, the formula returns the rotation. That can easily be proved as follows. Every material particle rotates exactly an angle $ϕ = ϕ_{a v}$ in such a way that

a = p \sin ϕ \int_{V} ρ | {\bar{X}}_{p} | | {\bar{x}}_{p} | d V, b = \cos ϕ \int_{V} ρ | {\bar{X}}_{p} | | {\bar{x}}_{p} | d V,

and, therefore,

\frac{a}{| a |} = p, \frac{| a |}{b} = \tan ϕ .

In all of these equations the direction vector $p$ is unknown. To determine $p$ we consider the characteristics of the displacement field of a rigid body rotation. For such a field, 1) the displacement of a particle is orthogonal to the rotation vector, and 2) the displacement is orthogonal to the position vector at half the motion. In a deformable body context we try to determine $p$ by forcing these two statements to be true in an average sense. Considering that we are looking strictly at a rotation with respect to the center of mass, its definition automatically ensures that the first statement is satisfied. The second condition can be written in the form

\int_{V} ρ [p \times (\bar{x} + \bar{X})] \times [\bar{x} - \bar{X}] d V = 0,

which is a homogeneous set of equations in the components of $p$ with coefficients made out of integrals of known quantities, from which $p$ can be solved. We can then calculate the projection of the old and new position onto the plane normal to $p$ with

{\bar{X}}_{p} = \bar{X} - p p \cdot \bar{X}, {\bar{x}}_{p} = \bar{x} - p p \cdot \bar{x},

and by simple substitution one then obtains

a = \int_{V} ρ [\bar{X} - (\bar{X} \cdot p) p] \times [\bar{x} - (\bar{x} \cdot p) p] d V = p p \cdot \int_{V} ρ \bar{X} \times \bar{x} d V,

where the vector $\int_{V} ρ \bar{X} \times \bar{x} d V$ is easily calculated from available quantities. With $p$ known, $a$ becomes determined. The quantity b can be calculated using with the same expressions, which yields

b = \int_{V} ρ [\bar{X} - (\bar{X} \cdot p) p] \cdot [\bar{x} - (\bar{x} \cdot p) p] d V = (I - p p) : \int_{V} ρ \bar{X} \bar{x} d V .

Once $a$ and b are known, $ϕ$ is readily determined.

The determination of the equivalent rigid body rotation is based on average particle translations. Rotational degrees of freedom are ignored in the calculation of this variable; it is assumed that such rotations will produce motions of points that will measurably contribute to the calculation. However, it is possible to find pathological cases in which that would not be the case; for instance, if the part of the model considered consists of rotary inertia elements only, the calculated average rigid body rotation will be calculated as zero, even if the elements have indeed rotated.