Formulation

Weak form

The configuration of a beam model with shear deformation is defined by a position vector \(\boldsymbol{x}\) and an associated orientation tensor \(\boldsymbol{\alpha}\in\mathrm{SO}(3)\). Let identify the position and orientation in the deformed configuration with an appended prime, \(\boldsymbol{x}'\) and \(\boldsymbol{\alpha}'\), with \(\boldsymbol{x}\) and \(\boldsymbol{\alpha}\) the corresponding position and orientation in the reference configuration.

The generalized beam strain measures are

\[\begin{split}\boldsymbol{\varepsilon} &= \boldsymbol{\alpha}'^T\boldsymbol{x}'_{/s}- \boldsymbol{\alpha}^T\boldsymbol{x}'_{/s}\\ \boldsymbol{\beta}& = \boldsymbol{\alpha}'^T\mathrm{ax}(\boldsymbol{\alpha}'_{/s}\boldsymbol{\alpha}'^T)- \boldsymbol{\alpha}^T\mathrm{ax}(\boldsymbol{\alpha}_{/s}\boldsymbol{\alpha}^T)\\ &= \boldsymbol{\alpha}^T\boldsymbol{\varPhi}^T\mathrm{ax}(\boldsymbol{\varPhi}_{/s}\boldsymbol{\varPhi}^T)\end{split}\]

where the axial operator \(\boldsymbol{a}=\mathrm{ax}(\boldsymbol{A})\) is the operator extracting the vector \(\boldsymbol{a}\) that characterize the skew symmetric part of tensor \(\boldsymbol{A}\), i.e. \(\boldsymbol{a}\times=1/2(\boldsymbol{A}^T+\boldsymbol{A})\) and \(\boldsymbol{\varPhi}=\boldsymbol{\alpha}'\boldsymbol{\alpha}^T\) is the rotation tensor bringing \(\boldsymbol{\alpha}\) into \(\boldsymbol{\alpha}'\). Remember also that, being \(\boldsymbol{\alpha}\) orthogonal, tensor \(\boldsymbol{\alpha}_{/s}\boldsymbol{\alpha}^T\) is skew-symmetric.

The internal forces and moments, \(\boldsymbol{T}\) and \(\boldsymbol{M}\) are functions of the generalized strains \(\boldsymbol{\varepsilon}\) and \(\boldsymbol{\beta}\). A linear constitutive law can be written as

\[\begin{split}\boldsymbol{T} &=\overline{\boldsymbol{K}}_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}} \boldsymbol{\varepsilon}+ \overline{\boldsymbol{K}}_{\boldsymbol{\varepsilon}\boldsymbol{\beta}}\boldsymbol{\beta}\\ \boldsymbol{M} &=\overline{\boldsymbol{K}}_{\boldsymbol{\beta}\boldsymbol{\varepsilon}} \boldsymbol{\varepsilon}+ \overline{\boldsymbol{K}}_{\boldsymbol{\beta}\boldsymbol{\beta}}\boldsymbol{\beta}\\\end{split}\]

where \(\overline{\boldsymbol{K}}_{\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}}\), \(\overline{\boldsymbol{K}}_{\boldsymbol{\varepsilon}\boldsymbol{\beta}}\), \(\overline{\boldsymbol{K}}_{\boldsymbol{\beta}\boldsymbol{\varepsilon}}\) and \(\overline{\boldsymbol{K}}_{\boldsymbol{\beta}\boldsymbol{\beta}}\), are constant tensors defined in beam section reference system.

The Principle of Virtual Work (i.e. the linear form) reads

(1)\[\int_L \left(\delta\boldsymbol{\varepsilon}\boldsymbol{T}+\delta\boldsymbol{\beta}\boldsymbol{M}\right) \mathrm{d}s= \int_L \left(\delta\boldsymbol{x}'\boldsymbol{t} + \boldsymbol{\varphi}_\delta\boldsymbol{m}\right) \mathrm{d}s+ \sum_i\delta\boldsymbol{x}'_i\boldsymbol{F}_i + \boldsymbol{\varphi}_{\delta i}\boldsymbol{C}_i\]

where the virtual rotation vector characterizes the virtual rotation of the orientation tensor in the deformed configuration, \(\boldsymbol{\varphi}_\delta = \mathrm{ax}(\delta\boldsymbol{\alpha}'\boldsymbol{\alpha}'^T) =\mathrm{ax}(\delta\boldsymbol{\varPhi}\boldsymbol{\varPhi}^T)\) and the terms \(\delta\boldsymbol{x}'_i\boldsymbol{F}_i + \boldsymbol{\varphi}_{\delta i}\boldsymbol{C}_i\) account for concentrated forces and moments.

Rotation parametrization

The weak form sketched in the previous Section is independent from any parametrization of the orientation and rotation tensors.

Dolfin, just like most of the finite element codes, cannot approximate a field of tensors \(\in\mathrm{SO}(3)\) with an interpolation scheme that respects the peculiar properties of the special orthogonal group they belong to (see [1, 3, 4] for details).

One of the more popular solution to this problem is to resort to a parametrization of the rotation field and interpolate the corresponding parameters.

The rotation vector \(\boldsymbol{\varphi}\) is one of the possible parametrizations of tensors \(\in\mathrm{SO}(3)\):

\[\begin{split}\boldsymbol{\varPhi} &= \exp(\boldsymbol{\varphi}\times)\\ &= \boldsymbol{I}+ a\boldsymbol{\varphi}\times+ b\boldsymbol{\varphi}\times\boldsymbol{\varphi}\times\end{split}\]

where

\[\begin{split}a &= \frac{\sin(\varphi)}{\varphi}\\ b &= \frac{1-\cos(\varphi)}{\varphi^2}\end{split}\]

and \(\varphi = \sqrt{\boldsymbol{\varphi}\boldsymbol{\varphi}}\).

The virtual rotation vector \(\boldsymbol{\varphi}_\delta\) is related to the virtual change of the rotation vector \(\delta\boldsymbol{\varphi}\) by

\[\boldsymbol{\varphi}_\delta = \boldsymbol{\varGamma}\delta\boldsymbol{\varphi}\]

where

\[\begin{split}\boldsymbol{\varGamma} &= \mathrm{dexp}(\boldsymbol{\varphi}\times)\\ &= \boldsymbol{I}+ b\boldsymbol{\varphi}\times+ c\boldsymbol{\varphi}\times\boldsymbol{\varphi}\times\end{split}\]

and

\[c = \frac{1-a}{\varphi^2}\]

Note that tensor \(\boldsymbol{\varGamma}\) also allows to compute

\[\boldsymbol{\varGamma}\boldsymbol{\varphi}_{/s} = \mathrm{ax}(\boldsymbol{\varPhi}_{/s}\boldsymbol{\varPhi}^T)\]

Coefficients \(a\), \(b\) and \(c\) are numerically ill-conditioned around 0, and should be computed using a Maclaurin expansion whenever \(\varphi\) goes below a given threshold. Their value in 0 is given by their limit for \(\varphi\rightarrow 0\).

By resorting to the chosen parametrization the Principle of Virtual Work Eq. (1) becomes

\[\int_L \left(\delta\boldsymbol{\varepsilon}\boldsymbol{T}+\delta\boldsymbol{\beta}\boldsymbol{M}\right) \mathrm{d}s= \int_L \left[\delta\boldsymbol{x}'\boldsymbol{t} + \left(\boldsymbol{\varGamma}\delta\boldsymbol{\varphi}\right)\boldsymbol{m}\right] \mathrm{d}s+ \sum_i\delta\boldsymbol{x}'_i\boldsymbol{F}_i + \left(\boldsymbol{\varGamma}\delta\boldsymbol{\varphi}\right)\boldsymbol{C}_i\]

where

\[\begin{split}\boldsymbol{\varphi}_\delta &= \boldsymbol{\varGamma}\delta\boldsymbol{\varphi}\\ \boldsymbol{\beta}&= \boldsymbol{\varPhi}^T\boldsymbol{\varGamma}\boldsymbol{\varphi}_{/s}\\ &=\boldsymbol{\varGamma}^T\boldsymbol{\varphi}_{/s}\\ \boldsymbol{\alpha}'&=\exp(\boldsymbol{\varphi}\times)\boldsymbol{\alpha}\end{split}\]

are now understood.

Tensor \(\boldsymbol{\varGamma}\) plays a crucial role here, as it appears into the explicit expressions of \(\delta\boldsymbol{\varepsilon}\) and \(\delta\boldsymbol{\beta}\) as a function of \(\delta\boldsymbol{\varphi}\). Tensor \(\delta\boldsymbol{\Gamma}\) is equal to \(\boldsymbol{I}\) for \(\varphi=0\). However, it has null determinant for rotations multiple of \(2\pi\), and this reflects the fact that the chosen parametrization is not bijective.

Complementary rotation vector

One of the solutions proposed since [5] to overcome the singularity of the rotation vector parametrization is to resort to the complementary rotation vector

(2)\[\boldsymbol{\varphi}_C = \boldsymbol{\varphi}\left(1-\frac{2\pi}{\varphi}\right)\]

whenever any rotation vector has magnitude \(\varphi > \pi\). The same solution was reproposed and discussed by many authors. As shown in [5] this “trick” makes sense since

(3)\[\boldsymbol{\varGamma(\boldsymbol{\varphi}_{C})}\boldsymbol{\varphi}_{C/s} = \boldsymbol{\varGamma(\boldsymbol{\varphi})}\boldsymbol{\varphi}_{/s}\]

thus it makes no difference, at the continuum level, whether the strains are computed with \(\boldsymbol{\varphi}\) or \(\boldsymbol{\varphi}_C\). The same of course holds for \(\boldsymbol{\varphi}_\delta\).

Things change when one move to an approximated rotation field, where the rotation vector is approximated using Continuous Galerkin CG elements, and one has to compute the complemnt vector of the discrete unknowns before approximating the rotation field. Let

\[\boldsymbol{\varphi({x})}=\sum_i w_i(\boldsymbol{x}) \boldsymbol{\varphi}_i\]

be the rotation vector approximated at point \(\boldsymbol{x}\). Then, it is easy to verify that

\[\boldsymbol{\varphi}_C = \boldsymbol{\varphi}\left(1-\frac{2\pi}{\varphi}\right) \neq \sum_i w_i(\boldsymbol{x}) \boldsymbol{\varphi}_{iC}\]

whenever vectors \(\boldsymbol{\varphi}_i\) are not coaxial. Thus, after reaching an equilibrium configuration, if one rigidly translates the structure, the new configuration will not be in an equilibrium state. This is something no author emphasize, but that should be kept in mind when not resorting to approximation schems specifically designed to deal with SO(3) fields as that explained in [1, 3, 4].

That said, this error tends to zero together with the dimension of the element, so one can choose to pragmatically live with it when computing an approximate solution.

The second difficulty arises from that fact that one cannot switch to the complement rotation inside the UFL form, and leave the discrete unknwon vectors \(\boldsymbol{\varphi}_i\) unchanged. Doing so would force to Dolfin to compute \(\boldsymbol{\varphi}\) as a function of the discrete unknowns, \(\boldsymbol{\varphi}_i\), then switch it to \(\boldsymbol{\varphi}_C\), and eventually use the derivative of the interpolated complement vector \(\boldsymbol{\varphi}_C\) with respect to the intepolated vector \(\boldsymbol{\varphi}\)

(4)\[\frac{\partial\boldsymbol{\varphi}_C}{\partial\boldsymbol{\varphi}} = \left(1 + \frac{2\pi}{\varphi}\right) \boldsymbol{I} - \frac{2\pi}{\varphi^3} \boldsymbol{\varphi}\otimes\boldsymbol{\varphi}.\]

The derivative is defined in such a way that the Gâteaux differential of \(\boldsymbol{\varphi}_C\) in the direction \(\boldsymbol{\varphi}\) can be computed as

(5)\[\delta\boldsymbol{\varphi}_C = \frac{\partial\boldsymbol{\varphi}_C}{\partial\boldsymbol{\varphi}} \cdot \delta\boldsymbol{\varphi}.\]

This would actually change nothing and leave everything with the same singularity of the parametrization, just as without using the complement vector.

Thus, one should change the discrete unknowns \(\boldsymbol{\varphi}_i\) into their complement rotation and let Dolfin do its magic without even konwing that something was changed under his foots.

This way the singularity of the parametrization is really worked around. However, the elements for which ony some of the unknown vectors \(\boldsymbol{\varphi}_i\) has been switched deserve a special tratment: one cannot interpolate the rotation vector field and get meaningful result for them, so for those element (and only for those) one should modify on the fly, and transparently from Dolfin, the coefficients \(\boldsymbol{\varphi}_i\) so that they are all consistent one with each other. After that, one should assemble the resulting contribution to the residual vector accounting by hand for the fact that some of the unknows were subject to the transformation of Eq. (2) bringing back their magnitute to a value \(\varphi>\pi\). The corresponding test functions are then transformed according to (5). Thus, the contribution to the residual \(\boldsymbol{f}_i\) computed by Dolfin shuld be tranformed into

(6)\[\left(\frac{\partial\boldsymbol{\varphi}_C}{\partial\boldsymbol{\varphi}}\right)^T \boldsymbol{f}_i = \left(\left(1 + \frac{2\pi}{\varphi}\right) \boldsymbol{I} - \frac{2\pi}{\varphi^3} \boldsymbol{\varphi}\otimes\boldsymbol{\varphi}\right) \boldsymbol{f}_i\]

before assembling it.

Similarly, the Jacobian matrix should be assembled by accounting for the fact that

1. both the test and trial function are transformed according to (5). the corresponding Dolfin jacobian \(\boldsymbol{J}_i\) rows must be multiplied by \(\left(\partial\boldsymbol{\varphi}_C/\partial\boldsymbol{\varphi}\right)^T\) on the left; the columns by \(\left(\partial\boldsymbol{\varphi}_C/\partial\boldsymbol{\varphi}\right)\) on the right;

2. one should actually linearize, (6), and thus add

   \[\frac{2\pi}{\varphi^3}\left[ \left( \boldsymbol{f}_i\otimes\boldsymbol{\varphi} + \boldsymbol{\varphi}\otimes\boldsymbol{f}_i \right) + (\boldsymbol{\varphi}\boldsymbol{f}_i) \left( \boldsymbol{I} - 3\ \frac{\boldsymbol{\varphi}}{\varphi}\otimes\frac{\boldsymbol{\varphi}}{\varphi} \right) \right]\]

   to account for the linearization of tensor \(\frac{\partial\boldsymbol{\varphi}_C}{\partial\boldsymbol{\varphi}}\) in Eq. (6).