Rotation variables

Even though Abaqus stores and outputs the rotation vector, quaternion parameters prove to be an efficient and convenient way to treat finite rotations computationally. Let $q_0\in R$ be a scalar, and let $\mathbf{q}\in R^3$ be a vector field. The quaternion $q$ is simply the pairing

To associate $q$ with the finite rotation vector $\mathit{\varphi }$, define the following:

By trigonometric identities it follows that the orthogonal matrix $\exp \left[\widehat{\mathit{\varphi }}\right]$ in Equation 1 is given in terms of $q$ as

By the convention introduced above, $\widehat{\mathbf{q}}$ is the skew-symmetric matrix with axial vector $\mathbf{q}$.

For a more detailed discussion of quaternion algebra and its relation to other representations of finite rotations, see the discussion by Spring (1986).

A compound rotation is the successive application of two or more rotation fields. In geometrically linear problems compound rotations are obtained simply as the linear superposition of the individual (linearized) rotation vectors. This fact follows directly from the series expansion for $\exp \left[\widehat{\mathit{\varphi }}\right]$. Let ${\mathit{\varphi }}_1$ and ${\mathit{\varphi }}_2$ be infinitesimal rotations. Thus, $\exp \left[{\widehat{\mathit{\varphi }}}_1\right]\approx \mathbf{I}+{\widehat{\mathit{\varphi }}}_1$, $\exp \left[{\widehat{\mathit{\varphi }}}_2\right]\approx \mathbf{I}+{\widehat{\mathit{\varphi }}}_2$, and

In geometrically nonlinear analysis compound rotations are no longer additive. Furthermore, they are not commutative; that is, the order of application is important. A significant exception occurs when the multiple rotations share the same rotation axis. This special case is investigated further below. A detailed example of a finite compound rotation is given in Conventions.

Let ${\mathbf{C}}_i$ be the orthogonal transformation representing the compound rotation defined as the product of a set of individual or incremental rotations $\mathrm{\Delta }{\mathbf{C}}_p$, for $p=1,2,\mathrm{\dots },i$. (For the case of specified boundary conditions ${\mathbf{C}}_i$ is the final product after i steps of all the specified rotations $\mathrm{\Delta }{\mathbf{C}}_p$; for the iterative numerical solution procedure ${\mathbf{C}}_i$ is the total rotation after i increments, where $\mathrm{\Delta }{\mathbf{C}}_p$, for $p=1,\mathrm{\dots },i$, is the converged rotation field solution at each increment.) By definition, the compound rotation is the product

or equivalently by the recursion relation,

It is important to note that $\mathrm{\Delta }\mathbf{C}\cdot \mathbf{C}$, which is interpreted as the finite rotation $\mathrm{\Delta }\mathbf{C}$ superposed on the finite rotation $\mathbf{C}$, is different from $\mathbf{C}\cdot \mathrm{\Delta }\mathbf{C}$, which is interpreted as the finite rotation $\mathbf{C}$ superposed on the finite rotation $\mathrm{\Delta }\mathbf{C}$.

Although compound rotations are defined in terms of orthogonal matrices, in a numerical context the rotation vectors (or equivalently the quaternion parameters) associated with the rotation matrices are the degrees of freedom. Compound rotations are performed as follows: Given a quaternion parametrization $q=\left(q_0,\mathbf{q}\right)$ and an incremental (finite) rotation $\mathrm{\Delta }q=\left(\mathrm{\Delta }{q}_0,\mathrm{\Delta }\mathbf{q}\right)$, where $\mathrm{\Delta }q$ is defined in terms of an incremental rotation vector $\mathrm{\Delta }\mathit{\varphi }$ by Equation 2, the total or compound rotation is given by the quaternion $r=\left(r_0,\mathbf{r}\right)$, which is calculated as

Here $\circ$ denotes the quaternion product and is defined as

Equation 4 allows for the update of rotation fields without ever calculating the orthogonal matrix from the quaternion and without performing a matrix multiplication. Furthermore, all operations are singularity free regardless of the magnitude of the incremental rotation field $\mathrm{\Delta }\mathit{\varphi }$. The final (total) rotation vector can be calculated from the quaternion $\mathbf{r}$ by inverting Equation 2.

For the special case when compound rotations share the same rotation axis, the compound rotation reduces to an additive form. Let $\mathrm{\Delta }q$ and $q$ have the same rotation axis $\mathbf{p}$. Then $q=\left(\cos \parallel \mathit{\varphi }/2\parallel ,\sin \parallel \mathit{\varphi }/2\parallel \mathbf{p}\right)$, $\mathrm{\Delta }q=\left(\cos \parallel \mathrm{\Delta }\mathit{\varphi }/2\parallel ,\sin \parallel \mathrm{\Delta }\mathit{\varphi }/2\parallel \mathbf{p}\right)$, and

which reduces to

For output purposes it is necessary to extract the rotation vector corresponding to a given quaternion. The extraction procedure is as follows: Let $r=\left(r_0,\mathbf{r}\right)$ be the quaternion, and let $\mathit{\varphi }$ be the rotation vector. Thus,

It is important to note that the extraction of the rotation vector from the quaternion is not unique. The magnitude $\parallel \mathit{\varphi }\parallel$ is determined only up to the addition of $2n\pi$, $n=\mathrm{\dots },-1,0,1,\mathrm{\dots }$Abaqus will always choose that rotation vector such that $\parallel \mathit{\varphi }\parallel <2\pi$.

As an example of the utility of the quaternion parameters, consider the incremental update of a director field for either a beam or shell analysis. At some stage of the solution the director field ${\mathbf{t}}^i$, the quaternion parametrization of the rotation field $q^i$, and the incremental rotation field $\mathrm{\Delta }{\mathit{\varphi }}^i$ are known at increment i. To update the director field by the incremental rotation to increment $i+1$, proceed as follows: First calculate the quaternion parametrization of the incremental rotation:

The director field at $i+1$ is then defined as ${\mathbf{t}}^{i+1}=\exp \left[{\widehat{\mathrm{\Delta }\mathit{\varphi }}}^i\right]\cdot {\mathbf{t}}^i$, where $\exp \left[{\widehat{\mathrm{\Delta }\mathit{\varphi }}}^i\right]$ is calculated with Equation 3. Thus, the director is calculated directly from the quaternion as

Furthermore, the update of the rotation field is obtained by quaternion multiplication $q^{i+1}=\mathrm{\Delta }{q}^i\circ q^i$ and is defined by

In the development of the balance equations, it is necessary to calculate the variation of the rotation field. Consider the vector field $\mathbf{a}$, which is obtained by rotation of the reference vector field $\mathbf{A}$:

Variations $\delta \mathbf{a}$ in this field are obtained as

where $\delta \mathbf{C}$ is the linearized rotation matrix; that is, the variation of the orthogonal tensor $\mathbf{C}$. On the other hand, the variation can be defined in terms of the linearized rotation field $\delta \mathit{\theta }$:

Consequently, it follows that

It is important to note that the linearized rotation $\delta \mathit{\theta }$, which is analogous to the angular velocity in dynamics, is not the variation of the rotation vector $\mathit{\varphi }$. By a straightforward (but involved) calculation, it can be shown that the variation of the rotation vector ($\delta \mathit{\varphi }$) is related to the linearized rotation $\delta \mathit{\theta }$ by

where

The inverse of $\mathbf{H}\left(\mathit{\varphi }\right)$ is

Let $\mathrm{d}\mathit{\theta }$ represent an infinitesimal change in the rotation field. A direct calculation of the variation of $\mathrm{d}\mathit{\theta }$, which is equivalent to calculation of the second variation of either $\mathbf{C}$ or $\mathbf{a}$, leads to an expression that is not symmetric in the variations $\delta \mathit{\theta }$ and the changes $\mathrm{d}\mathit{\theta }$. However, it is shown in Simo (1992) that the correct definition of the Hessian operator—that is, the “covariant” derivative of the weak form of the balance equations—requires only the symmetric part (with respect to the variations) of the second variation. Thus, without loss of generality, we can write

Similarly, the second variation of the vector field rotated by $\mathbf{C}$ can be written as

Taking the time derivative of the rotation matrix, we find with the same arguments as used in the calculation of the variations that

where $\widehat{\mathit{\omega }}$ is the angular velocity matrix. Equivalently, the first and second time derivative of $\mathbf{a}$ are written as

The instantaneous angular velocity vector $\mathit{\omega }$ is related to the time rate of change of the rotation vector by the relation

where $\mathbf{H}\left(\mathit{\varphi }\right)$ is given by Equation 6.

In the linearization of the dynamic balance equations, it is necessary to calculate the variation of the angular velocity, $\mathrm{d}\mathit{\omega }$. This quantity, however, can be calculated only by linearizing the specific algorithm used for the time integration of the dynamic equations.

Next, a more rigorous treatment of the two-dimensional constant velocity joint described in MPC is presented. This derivation exemplifies some of the issues associated with the treatment of finite rotations. Uniform collapse of straight and curved pipe segments deals with a different finite rotation constraint and tackles additional complications.

Let a, b, c (see Figure 1) be the nodes making up the joint, with a the dependent node.

Figure 1. Nonlinear MPC example—constant velocity joint.

The joint is operated by prescribing an axial rotation ${\mathit{\varphi }}^c={\varphi }^c{\mathbf{e}}_x$ at c and an out-of-plane rotation ${\mathit{\varphi }}^b={\varphi }^b{\mathbf{e}}_z$ at b. The compounding of these two prescribed rotation fields will determine the total rotation at a. We can formally write this constraint as follows:

The constraint can be written in terms of the rotation matrices as

With the previously defined variational expressions, the constraint can be linearized as

This expression can be simplified by right-multiplying the expression by ${\mathbf{C}}^T\left({\mathit{\varphi }}^a\right)$ and by making use of the constraint Equation 7, which yields

which can be written in vector form as

Since

the linearized constraint is indeed identical to the one derived based on simple linear considerations in Linear constraint equations.

The linearized constraint is used for the calculation of equilibrium. It can also be used for the recovery of the dependent rotation, ${\mathit{\varphi }}^a$, as is done in Linear constraint equations. The resulting rotation will satisfy the constraint approximately (unless one of the angles ${\varphi }^b$ or ${\varphi }^c$ is constant, in which case the constraint is linear and the recovery is exact).

For an exact enforcement of the constraint, user subroutine MPC must define the components of the total rotation vector ${\mathit{\varphi }}^a$ exactly. To do so, ${\mathit{\varphi }}^a$ must be updated based on the current values of ${\mathit{\varphi }}^b$ and ${\mathit{\varphi }}^c$. This is most easily accomplished with the aid of the quaternion parameters. Let $q^b=\left(\cos \left({\varphi }^b/2\right),\sin \left({\varphi }^b/2\right){\mathbf{e}}_z\right)$ and $q^c=\left(\cos \left({\varphi }^c/2\right),\sin \left({\varphi }^c/2\right){\mathbf{e}}_x\right)$ be the quaternion parameterizations associated with the finite rotation vectors ${\mathit{\varphi }}^b$ and ${\mathit{\varphi }}^c$, respectively. The total compound rotation ${\mathit{\varphi }}^a$ is given by the quaternion $q^a=\left(q_0^a,{\mathbf{q}}^a\right)$, where

according to the quaternion compound formula Equation 4. The rotation vector ${\mathit{\varphi }}^a$ is extracted from the quaternion $q^a$ as follows:

where $\parallel {\mathbf{q}}^a\parallel$ is the norm of the vector ${\mathbf{q}}^a$.

MPC shows the implementation of the linearized form of the constraint in user subroutine MPC. The implementation of the exact nonlinear constraint is shown below:

SUBROUTINE MPC(UE,A,JDOF,MDOF,N,JTYPE,X,U,UINIT,MAXDOF,LMPC,
      * KSTEP,KINC,TIME,NT,NF,TEMP,FIELD)
 C
       INCLUDE 'ABA_PARAM.INC'
 C
       DIMENSION UE(MDOF), A(MDOF,MDOF,N), JDOF(MDOF,N), X(6,N),
      * U(MAXDOF,N), UINIT(MAXDOF,N), TIME(2), TEMP(NT,N),
      * FIELD(NF,NT,N)
       PARAMETER( SMALL = 1.E-14 )
 C
       IF ( JTYPE .EQ. 1 ) THEN
         A(1,1,1) =  1.
         A(2,2,1) =  1.
         A(3,3,1) =  1.
         A(3,1,2) = -1.
         A(1,1,3) = -COS(U(6,2))
         A(2,1,3) = -SIN(U(6,2))
 C
         JDOF(1,1) = 4
         JDOF(2,1) = 5
         JDOF(3,1) = 6
         JDOF(1,2) = 6
         JDOF(1,3) = 4
 C
         CPHIB = COS(0.5*U(6,2))
         SPHIB = SIN(0.5*U(6,2))
         CPHIC = COS(0.5*U(4,3))
         SPHIC = SIN(0.5*U(4,3))
 C
         QA0 = CPHIB*CPHIC
         QAX = CPHIB*SPHIC
         QAY = SPHIB*SPHIC
         QAZ = CPHIB*SPHIC
 C
         QAMAG = SQRT( QAX*QAX + QAY*QAY + QAZ*QAZ )
         IF ( QAMAG .GT. SMALL ) THEN
           PHIA  = 2.*ATAN2( QAMAG , QA0 )
           UE(1) = PHIA*QAX/QAMAG
           UE(2) = PHIA*QAY/QAMAG
           UE(3) = PHIA*QAZ/QAMAG
         ELSE
           UE(1) = 0.
           UE(2) = 0.
           UE(3) = 0.
         END IF
       END IF
 C
       RETURN
       END