<!-- source-page: 341 -->

# Elements

The bicurvature defines the axial strain variation in the section due to the twist of the beam. The expression for the curvature and the twist can be combined to yield

$$
\frac {d \mathbf {n} _ {\alpha}}{d S} = \epsilon_ {\alpha} ^ {\beta} (- b _ {\beta} \mathbf {t} + b \mathbf {n} _ {\beta}).
$$

Before we derive strain measures from these expressions, we will consider in detail how the above quantities and their first and second variations are obtained for a typical beam finite element.

# Beam element in undeformed configuration

In the undeformed configuration, we use capital letters for all quantities. We assume that the undeformed state has no warping, so the position of a material point is given by

$$
\hat {\mathbf {X}} (S, S ^ {\alpha}) = \mathbf {X} (S) + S ^ {\alpha} \mathbf {N} _ {\alpha} (S)
$$

and the curvatures and twist are defined by

$$
B _ {\alpha} = \epsilon_ {\alpha} ^ {\beta} \mathbf {T} \cdot \frac {d \mathbf {N} _ {\beta}}{d S},
$$

$$
B = \mathbf {N} _ {2} \cdot \frac {d \mathbf {N} _ {1}}{d S} = - \mathbf {N} _ {1} \cdot \frac {d \mathbf {N} _ {2}}{d S},
$$

where T is the unit vector orthogonal to ${ \bf N } _ { 1 }$ and $\mathbf { N } _ { 2 } ;$ ; i.e.,

$$
\mathbf {T} = \mathbf {N} _ {1} \times \mathbf {N} _ {2}.
$$

We assume that the section normal T coincides with the beam tangent S.

In the element the position of a point on the axis is interpolated from nodal positions ${ \mathbf { X } } ^ { N }$ with standard interpolation functions as

$$
\mathbf {X} = g _ {N} (\boldsymbol {\xi}) \mathbf {X} ^ {N},
$$

where $\xi$ is a parametric coordinate, typically varying between ¡1 and 1 along the element. The beam axis tangent is readily computed as

$$
\mathbf {S} = \left. \frac {d \mathbf {X}}{d \xi} \right/ \frac {d S}{d \xi} = \frac {d g _ {N}}{d \xi} \mathbf {X} ^ {N} \bigg / \left| \frac {d g _ {M}}{d \xi} \mathbf {X} ^ {M} \right|.
$$

The section normals are interpolated from user-defined nodal normals ${ { \bf N } ^ { N } }$ . However, we cannot use a simple interpolation since this would not create integration point normals orthogonal to the beam tangent S. Hence, we use a two-step approach. First, we create approximate normals by the interpolation

$$
\overline {{\mathbf {N}}} _ {\alpha} = g _ {N} (\boldsymbol {\xi}) \mathbf {N} _ {\alpha} ^ {N}.
$$

<!-- source-page: 342 -->

# Elements

Then, we orthonormalize these vectors with respect to S by

$$
\tilde {\mathbf {N}} _ {\alpha} = (\overline {{{\mathbf {N}}}} _ {\alpha} - \mathbf {S} \mathbf {S} \cdot \overline {{{\mathbf {N}}}} _ {\alpha}) / \left| (\overline {{{\mathbf {N}}}} _ {\alpha} - \mathbf {S} \mathbf {S} \cdot \overline {{{\mathbf {N}}}} _ {\alpha}) \right| \qquad \mathrm{(nosummation)}
$$

and subsequently, with respect to each other by

$$
\mathbf {N} _ {\alpha} = (\tilde {\mathbf {N}} _ {\alpha} + \epsilon_ {\alpha} ^ {\beta} \tilde {\mathbf {N}} _ {\beta} \times \mathbf {S}) / | (\tilde {\mathbf {N}} _ {\alpha} + \epsilon_ {\alpha} ^ {\beta} \tilde {\mathbf {N}} _ {\beta} \times \mathbf {S}) | \qquad \mathrm{(nosumon} \alpha),
$$

where it has been assumed that $\overline { { \mathbf { N } } } _ { \alpha }$ and S form a right-handed system. This provides $\mathbf { T } = \mathbf { S }$ . The curvature and the twist in the initial configuration are calculated directly from $\overline { { \mathbf { N } } } _ { \alpha }$ as

$$
\begin{array}{l} B _ {\alpha} = \epsilon_ {\beta} ^ {\alpha} \mathbf {T} \cdot \frac {d \overline {{\mathbf {N}}} _ {\beta}}{d S} = \epsilon_ {\alpha} ^ {\beta} \mathbf {T} \cdot \frac {d g _ {N}}{d \xi} \mathbf {N} _ {\alpha} ^ {N} \bigg / \frac {d S}{d \xi}, \\ { B } { = \frac { 1 } { 2 } \epsilon _ { \alpha } ^ { \beta } \overline { { { \mathbf { N } } } } _ { \beta } \cdot \frac { d \overline { { { \mathbf { N } } } } _ { \alpha } } { d S } = \frac { 1 } { 2 } \epsilon _ { \alpha } ^ { \beta } \overline { { { \mathbf { N } } } } _ { \beta } \cdot \frac { d g _ { N } } { d \xi } \mathbf { N } ^ { N } \Big / \frac { d S } { d \xi } . } \\ \end{array}
$$

The "average" twist is taken, since, in general,

$$
\overline {{\mathbf {N}}} _ {2} \cdot \frac {d \overline {{\mathbf {N}}} _ {1}}{d S} \neq - \overline {{\mathbf {N}}} _ {1} \cdot \frac {d \overline {{\mathbf {N}}} _ {2}}{d S}.
$$

The gradient of the normal vectors is then obtained as

$$
\frac {d \mathbf {N} _ {\alpha}}{d S} = \epsilon_ {\alpha} ^ {\beta} (- B _ {\beta} \mathbf {T} + B \mathbf {N} _ {\beta});
$$

and, therefore, the curvature and twist are also equal to

$$
\begin{array}{l} B _ {\alpha} = \epsilon_ {\alpha} ^ {\beta} \mathbf {T} \cdot \frac {d \mathbf {N} _ {\beta}}{d S}, \\ { B } { = \mathbf { N } _ { 2 } \cdot \frac { d \mathbf { N } _ { 1 } } { d S } = - \mathbf { N } _ { 1 } \cdot \frac { d \mathbf { N } _ { 2 } } { d S } . } \\ \end{array}
$$

The procedure followed above to derive $\mathbf { N } _ { \alpha }$ and dN® $\frac { d \mathbf { N } _ { \alpha } } { d S }$ is not unique but provides values that satisfy the proper orthonormality conditions. This procedure is followed only for the undeformed configuration. For subsequent configurations $\mathbf { n } _ { \alpha }$ and $\frac { d { \bf n } _ { \alpha } } { d S }$ are obtained individually by forward integration of the kinematic equations.

# Change of position, warping, and normal direction

We assume that the position of the beam axis and the orientation of the normals can undergo (independent) changes. The change in position of the axis is described by the velocity vector $\mathbf { v } = \mathbf { v } ( { \boldsymbol { \xi } } )$ , which can be obtained from nodal velocities $\mathbf { v } ^ { N }$ with the standard interpolation functions as

<!-- source-page: 343 -->

# Elements

$$
\mathbf {v} (\xi) = g _ {N} (\xi) \mathbf {v} ^ {N}.
$$

The change in orientation of the normals is described by the spin vector $\omega = \omega ( \xi )$ , which is obtained from nodal spin vectors $\omega ^ { N }$ with the same interpolation functions $g _ { N } ( \xi )$ , as

$$
\dot {\mathbf {n}} _ {\alpha} = \pmb {\omega} (\pmb {\xi}) \times \mathbf {n} _ {\alpha}, \pmb {\omega} (\pmb {\xi}) = g _ {N} (\pmb {\xi}) \pmb {\omega} ^ {N}.
$$

Rigid body motion is included since the original position $\mathbf { X } ( \xi )$ is obtained by the same interpolation as $\mathbf { v } ( \boldsymbol { \xi } )$ . The rate of change of warping is also defined in terms of the rate of change of nodal warping $w ^ { N }$ with the standard interpolation functions as

$$
\dot {w} (\xi) = g _ {N} (\xi) \dot {w} ^ {N}.
$$

The velocity and spin describe the rate of change of the position and orientation. Finite changes in position are obtained by integration of the velocities over a finite time increment as

$$
\Delta \mathbf {x} (\xi) = \int_ {t} ^ {t + \Delta t} \mathbf {v} (\xi) d t = g _ {N} (\xi) \int_ {t} ^ {t + \Delta t} \mathbf {v} ^ {N} d t = g _ {N} (\xi) \Delta \mathbf {X} ^ {N}.
$$

Similarly, for the warping,

$$
\Delta w (\xi) = \int_ {t} ^ {t + \Delta t} \dot {w} (\xi) d t = g _ {N} (\xi) \int_ {t} ^ {t + \Delta t} \dot {w} ^ {N} d t = g _ {N} (\xi) \Delta w ^ {N}.
$$

The spin is related to the rate of change of the rotation quaternion q by

$$
2 \pmb {q} ^ {\dagger} \dot {\pmb {q}} = \pmb {\omega}, \quad \pmb {q} = \left(\cos \frac {\theta}{2}, \sin \frac {\theta}{2} \mathbf {n}\right), \quad \pmb {q} ^ {\dagger} = \left(\cos \frac {\theta}{2}, - \sin \frac {\theta}{2} \mathbf {n}\right),
$$

where ${ \pmb \theta } = \theta { \bf n }$ is the total or Euler rotation.

The relation between spin and quaternion can be integrated exactly if it is assumed that the spin is constant over the time increment $\Delta t .$ . We then define

$$
\Delta \phi = \omega \Delta t, \Delta \phi = | \Delta \phi |, \mathbf {n} = \Delta \phi / \Delta \phi ,
$$

and the incremental rotation quaternion follows from

$$
\Delta \pmb {q} = \left(\cos \frac {\Delta \phi}{2}, \sin \frac {\Delta \phi}{2} \mathbf {n}\right).
$$

This allows us to update the normal directions at the nodes and at stress points with

$$
\mathbf {n} _ {\alpha} = \Delta \pmb {q} \mathbf {N} _ {\alpha} \Delta \pmb {q} ^ {\dagger}.
$$

<!-- source-page: 344 -->

# Elements

In this equation we use the notation that $\mathbf { N } _ { \alpha }$ is $\mathbf { n } _ { \alpha }$ at the beginning of the current increment; i.e.,

$$
\mathbf {N} _ {\alpha} \stackrel {\mathrm{def}} {=} \mathbf {n} _ {\alpha} (t).
$$

For the entire motion the new normal directions can be formally expressed as

$$
\mathbf {n} _ {\alpha} = \pmb {q} \mathbf {N} _ {\alpha} \pmb {q} ^ {\dagger},
$$

where q is defined by the product rule

$$
\boldsymbol {q} = \prod_ {i = 1} ^ {n} \Delta \boldsymbol {q} ^ {(i)}.
$$

Here i is an increment and $\Delta \pmb q ^ { ( i ) }$ is the rotation quaternion for that increment. The section normal t is updated the same way.

# Change of curvature and twist

The curvature and the twist involve the derivative of the normal vector with respect to S. From the update rule for $\mathbf { n } _ { \alpha }$ ,

$$
\frac {d \mathbf {n} _ {\alpha}}{d S} = \frac {d \Delta \pmb {q}}{d S} \mathbf {N} _ {\alpha} \Delta \pmb {q} ^ {\dagger} + \Delta \pmb {q} \mathbf {N} _ {\alpha} \frac {d \Delta \pmb {q} ^ {\dagger}}{d S} + \Delta \pmb {q} \frac {d \mathbf {N} _ {\alpha}}{d S} \Delta \pmb {q} ^ {\dagger}.
$$

The second term can be written as

$$
\Delta \pmb {q} \mathbf {N} _ {\alpha} \frac {d \Delta \pmb {q} ^ {\dagger}}{d S} = \left[ \frac {d \Delta \pmb {q}}{d S} \mathbf {N} _ {\alpha} ^ {\dagger} \Delta \pmb {q} ^ {\dagger} \right] ^ {\dagger} = - \left[ \frac {d \Delta \pmb {q}}{d S} \mathbf {N} _ {\alpha} \Delta \pmb {q} ^ {\dagger} \right] ^ {\dagger}.
$$

Hence, the scalar parts of the first two terms cancel each other, and the vector parts are the same. Therefore,

$$
\frac {d \mathbf {n} _ {\alpha}}{d S} = 2 \mathbf {V} \left(\frac {d \Delta \pmb {q}}{d S} \mathbf {N} _ {\alpha} \Delta \pmb {q} ^ {\dagger}\right) + \Delta \pmb {q} \frac {d \mathbf {N} _ {\alpha}}{d S} \Delta \pmb {q} ^ {\dagger}.
$$

Because $\Delta \pmb q$ is a rotation quaternion, its inverse is equal to its conjugate $( \Delta \pmb q ^ { - 1 } = \Delta \pmb q ^ { \dag } )$ , and, hence, we can write

$$
\begin{array}{l} \frac {d \mathbf {n} _ {\alpha}}{d S} = 2 \mathbf {V} \left(\frac {d \Delta \mathbf {q}}{d S} \Delta \mathbf {q} ^ {\dagger} \Delta \mathbf {q} \mathbf {N} _ {\alpha} \Delta \mathbf {q} ^ {\dagger}\right) + \Delta \mathbf {q} \frac {d \mathbf {N} _ {\alpha}}{d S} \Delta \mathbf {q} ^ {\dagger} \\ = \mathbf {V} \left(2 \frac {d \Delta \pmb {q}}{d S} \Delta \pmb {q} ^ {\dagger} \mathbf {n} _ {\alpha}\right) + \Delta \pmb {q} \frac {d \mathbf {N} _ {\alpha}}{d S} \Delta \pmb {q} ^ {\dagger}, \\ \end{array}
$$

where $\mathbf { p } = \mathbf { V } ( \pmb { p } )$ denotes the vector part of a quaternion. For the first term we use the relation

<!-- source-page: 345 -->

# Elements

$$
\frac {d \Delta \pmb {q}}{d S} = \left(- \frac {1}{2} \sin \frac {\Delta \phi}{2} \frac {d \Delta \phi}{d S}, \frac {1}{2} \cos \frac {\Delta \phi}{2} \mathbf {n} \frac {d \Delta \phi}{d S} + \sin \frac {\Delta \phi}{2} \frac {d \mathbf {n}}{d S}\right).
$$

This leads to

$$
2 \frac {d \Delta \pmb {q}}{d S} \Delta \pmb {q} ^ {\dagger} = \mathbf {n} \frac {d \Delta \phi}{d S} + \sin \Delta \phi \frac {d \mathbf {n}}{d S} + (1 - \cos \Delta \phi) \mathbf {n} \times \frac {d \mathbf {n}}{d S} = \Delta \pmb {\varrho},
$$

which is a vector. From the definitions of $\Delta \phi$ and n it follows that

$$
\frac {d \Delta \phi}{d S} = \Delta \phi^ {- 1} \Delta \phi \cdot \frac {d \Delta \phi}{d S} = \mathbf {n} \cdot \frac {d \Delta \phi}{d S},
$$

$$
\frac {d \mathbf {n}}{d S} = \Delta \phi^ {- 1} \frac {d \Delta \pmb {\phi}}{d S} - \Delta \phi^ {- 2} \Delta \pmb {\phi} \mathbf {n} \cdot \frac {d \Delta \pmb {\phi}}{d S} = \Delta \phi^ {- 1} \left[ \frac {d \Delta \pmb {\phi}}{d S} - \mathbf {n} \mathbf {n} \cdot \frac {d \Delta \pmb {\phi}}{d S} \right].
$$

These results provide

$$
\begin{array}{l} \Delta \pmb {\varrho} = \textbf {n n} \cdot \frac {d \Delta \pmb {\phi}}{d S} + \frac {\sin \Delta \phi}{\Delta \phi} \frac {d \Delta \pmb {\phi}}{d S} - \frac {\sin \Delta \phi}{\Delta \phi} \textbf {n n} \cdot \frac {d \Delta \pmb {\phi}}{d S} + \frac {(1 - \cos \Delta \phi)}{\Delta \phi} \textbf {n} \times \frac {d \Delta \pmb {\phi}}{d S} \\ = \frac {\sin \Delta \phi}{\Delta \phi} \frac {d \Delta \pmb {\phi}}{d S} + \frac {(1 - \cos \Delta \phi)}{\Delta \phi} \mathbf {n} \times \frac {d \Delta \pmb {\phi}}{d S} + (1 - \frac {\sin \Delta \phi}{\Delta \phi}) \mathbf {n} \mathbf {\nabla} \mathbf {n} \cdot \frac {d \Delta \pmb {\phi}}{d S}. \\ \end{array}
$$

We see that $\Delta \pmb { \varrho } \Rightarrow \frac { d \Delta \phi } { d S } \mathrm { i f } \Delta \phi  0 .$ .

For the second term we express $\frac { d \mathbf { N } _ { \alpha } } { d S }$ in terms of the curvature and twist at the beginning of the increment, which yields

$$
\Delta \pmb {q} \frac {d \mathbf {N} _ {\alpha}}{d S} \Delta \pmb {q} ^ {\dagger} = \Delta \pmb {q} \epsilon_ {\alpha} ^ {\beta} (- B _ {\beta} \mathbf {T} + B \mathbf {N} _ {\beta}) \Delta \pmb {q} ^ {\dagger} = \epsilon_ {\alpha} ^ {\beta} (- B _ {\beta} \mathbf {t} + B \mathbf {n} _ {\beta}).
$$

Combining these terms,

$$
\frac {d \mathbf {n} _ {\alpha}}{d S} = \Delta \pmb {\varrho} \times \mathbf {n} _ {\alpha} + \epsilon_ {\alpha} ^ {\beta} (- B _ {\beta} \mathbf {t} + B \mathbf {n} _ {\beta}).
$$

The current curvature and twist are, therefore,

$$
b _ {\alpha} = \epsilon_ {\alpha} ^ {\beta} \mathbf {t} \cdot \frac {d \mathbf {n} _ {\beta}}{d S} = \epsilon_ {\alpha} ^ {\beta} \Delta \pmb {\varrho} \cdot (\mathbf {n} _ {\beta} \times \mathbf {t}) - \epsilon_ {\alpha} ^ {\beta} \epsilon_ {\beta} ^ {\gamma} B _ {\gamma} = \Delta \pmb {\varrho} \cdot \mathbf {n} _ {\alpha} + B _ {\alpha}
$$

and

$$
{ b } { = \mathbf { n } _ { 2 } \cdot \frac { d \mathbf { n } _ { 1 } } { d S } = \Delta \pmb { \varrho } \cdot ( \mathbf { n } _ { 1 } \times \mathbf { n } _ { 2 } ) + \epsilon _ { 1 } ^ { 2 } B = \Delta \pmb { \varrho } \cdot \mathbf { t } + B . }
$$

<!-- source-page: 346 -->

Hence, the current curvature and twist are updated by a summation over all increments as given by

$$
b _ {\alpha} ^ {(n)} = \sum_ {i = 1} ^ {n} \Delta \pmb {\varrho} ^ {(i)} \cdot \mathbf {n} _ {\alpha} ^ {(i)} + B _ {\alpha},
$$

$$
{b ^ {(n)} =} {\sum_ {i = 1} ^ {n} \Delta \pmb {\varrho} ^ {(i)} \cdot \mathbf {t} _ {\alpha} ^ {(i)} + B.}
$$

# First variations

The first variations of the geometric quantities are readily obtained. Recall that

$$
\lambda = \left| \frac {d \mathbf {x}}{d S} \right|,
$$

$$
\gamma_ {\alpha} = \mathbf {s} \cdot \mathbf {n} _ {\alpha} \qquad (\mathrm{with} \mathbf {s} = \lambda^ {- 1} \frac {d \mathbf {x}}{d S}),
$$

$$
b _ {\alpha} = \epsilon_ {\alpha} ^ {\beta} \mathbf {t} \cdot \frac {d \mathbf {n} _ {\beta}}{d S},
$$

$$
{ b } { = \mathbf { n } _ { 2 } \cdot \frac { d \mathbf { n } _ { 1 } } { d S } , }
$$

$$
\chi = \frac {d w}{d S}.
$$

It follows that

$$
\delta \lambda = \mathbf {s} \cdot \frac {d \delta \mathbf {x}}{d S},
$$

$$
\delta \gamma_ {\alpha} = \delta \mathbf {s} \cdot \mathbf {n} _ {\alpha} + \mathbf {s} \cdot \delta \mathbf {n} _ {\alpha} = - \lambda^ {- 1} \gamma_ {\alpha} \delta \lambda + \lambda^ {- 1} \frac {d \delta \mathbf {x}}{d S} \cdot \mathbf {n} _ {\alpha} + \mathbf {s} \cdot (\delta \pmb {\phi} \times \mathbf {n} _ {\alpha})
$$

$$
\approx \lambda^ {- 1} \mathbf {n} _ {\alpha} \cdot \frac {d \delta \mathbf {x}}{d S} - \epsilon_ {\alpha} ^ {\beta} \mathbf {n} _ {\beta} \cdot \delta \pmb {\phi},
$$

$$
\delta \chi = \frac {d \delta w}{d S},
$$

where in the expression for $\delta \gamma _ { \alpha }$ we have assumed that $\gamma _ { \alpha } \simeq 0$ and $\mathbf { s } \simeq \mathbf { t }$ . For the variations in the curvature and twist, we note that $\delta \mathbf { t } = \delta \phi \times \mathbf { t }$ . Hence, it follows that

$$
\begin{array}{l} \delta b _ {\alpha} = \epsilon_ {\alpha} ^ {\beta} \left(\delta \mathbf {t} \cdot \frac {d \mathbf {n} _ {\beta}}{d S} + \mathbf {t} \cdot \frac {d \delta \mathbf {n} _ {\beta}}{d S}\right) = \epsilon_ {\alpha} ^ {\beta} \left[ (\delta \pmb {\phi} \times \mathbf {t}) \cdot \frac {d \mathbf {n} _ {\beta}}{d S} + \mathbf {t} \cdot (\delta \pmb {\phi} \times \frac {d \mathbf {n} _ {\beta}}{d S} + \frac {d \delta \pmb {\phi}}{d S} \times \mathbf {n} _ {\beta}) \right] \\ = \epsilon_ {\alpha} ^ {\beta} \frac {d \delta \pmb {\phi}}{d S} \cdot (\mathbf {n} _ {\beta} \times \mathbf {t}) = \mathbf {n} _ {\alpha} \cdot \frac {d \delta \pmb {\phi}}{d S}, \\ \end{array}
$$

$$
\delta b = \delta \mathbf {n} _ {2} \cdot \frac {d \mathbf {n} _ {1}}{d S} + \mathbf {n} _ {2} \cdot \frac {d \delta \mathbf {n} _ {1}}{d S} = (\delta \pmb {\phi} \times \mathbf {n} _ {2}) \cdot \frac {d \mathbf {n} _ {1}}{d S} + \mathbf {n} _ {2} \cdot (\delta \pmb {\phi} \times \frac {d \mathbf {n} _ {1}}{d S} + \frac {d \delta \pmb {\phi}}{d S} \times \mathbf {n} _ {1})
$$

$$
= \frac {d \delta \pmb {\phi}}{d S} \cdot (\mathbf {n} _ {1} \times \mathbf {n} _ {2}) = \mathbf {t} \cdot \frac {d \delta \pmb {\phi}}{d S}.
$$

<!-- source-page: 347 -->

# Elements

These terms will be used in the virtual work equation that will be discussed later. In using the above expressions the rotational quantities are obtained by interpolation from nodal variational quantities that are assumed to be valid for the velocity fields as

$$
\delta \mathbf {x} (\xi) = g _ {N} (\xi) \delta \mathbf {x} ^ {N},
$$

$$
\delta w (\xi) = g _ {N} (\xi) \delta w ^ {N},
$$

$$
\delta \pmb {\phi} (\xi) = g _ {N} (\xi) \delta \pmb {\phi} ^ {N}.
$$

# Solution with Newton's algorithm

Newton's algorithm involves linearization of the incremental equations. The equations must be linearized around the current (latest) state. At the integration points these equations simply take the form

$$
\mathrm{d} \lambda = \mathbf {s} \cdot \frac {d \mathrm{d} \mathbf {x}}{d S},
$$

$$
\mathrm{d} \gamma_ {\alpha} \approx \lambda^ {- 1} \mathbf {n} _ {\alpha} \cdot \frac {d \mathrm{d} \mathbf {x}}{d S} - \epsilon_ {\alpha} ^ {\beta} \mathbf {n} _ {\beta} \cdot \mathrm{d} \pmb {\phi},
$$

$$
\mathrm{d} b _ {\alpha} = \frac {d \mathrm{d} \phi}{d S} \cdot \mathbf {n} _ {\alpha},
$$

$$
{\mathrm{d} b} {= \frac {d \mathrm{d} \phi}{d S} \cdot \mathbf {t},}
$$

$$
\mathrm{d} \chi = \frac {d \mathrm{d} w}{d S}.
$$

The corrections dx and dw are readily obtained from the nodal corrections with the interpolation functions $g ^ { N } ( \xi )$ as

$$
\mathrm{d} \mathbf {x} = g ^ {N} (\xi) \mathrm{d} \mathbf {x} ^ {N}, \mathrm{d} w = g ^ {N} (\xi) \mathrm{d} w ^ {N}.
$$

The corrections $\mathrm { d } \phi$ are defined with respect to the end of the increment: we call such corrections "compound" corrections. However, it has been assumed that the total incremental rotation $\Delta \phi$ would be interpolated with $g ^ { N } ( \xi )$ as

$$
\Delta \phi (\xi) = g ^ {N} (\xi) \Delta \phi^ {N}.
$$

Linearization of this equation yields

$$
\mathrm{d} \Delta \boldsymbol {\phi} (\xi) = g ^ {N} (\xi) \mathrm{d} \Delta \boldsymbol {\phi} ^ {N},
$$

where $\boldsymbol { \mathrm { l } } \Delta \phi$ are corrections in $\Delta \phi$ in an additive sense such that

$$
\Delta \phi^ {n e w} = \Delta \phi^ {o l d} + \mathrm{d} \Delta \phi .
$$

<!-- source-page: 348 -->

# Elements

To relate the additive correction d $1 \Delta \phi$ to the compound correction $\mathrm { d } \phi ,$ we use the formula obtained for $\Delta \phi$ to find

$$
\mathrm{d} \pmb {\phi} = \Delta \pmb {q} ^ {\dagger} \mathrm{d} \pmb {q} = \frac {\sin \Delta \phi}{\Delta \phi} \mathrm{d} \Delta \pmb {\phi} + \frac {(1 - \cos \Delta \phi)}{\Delta \phi} \mathbf {n} \times \mathrm{d} \Delta \pmb {\phi} + (1 - \frac {\sin \Delta \phi}{\Delta \phi}) \mathbf {n} \mathbf {n} \cdot \mathrm{d} \Delta \pmb {\phi}.
$$

Similarly, at the nodes we find

$$
\begin{array}{l} \mathrm{d} \pmb {\phi} ^ {N} = \frac {\sin \Delta \phi^ {N}}{\Delta \phi^ {N}} \mathrm{d} \Delta \pmb {\phi} ^ {N} + \frac {(1 - \cos \Delta \phi^ {N})}{\Delta \phi^ {N}} \mathbf {n} ^ {N} \times \mathrm{d} \Delta \pmb {\phi} ^ {N} \\ + (1 - \frac {\sin \Delta \phi^ {N}}{\Delta \phi^ {N}}) \mathbf {n} ^ {N} \mathbf {n} ^ {N} \cdot \mathrm{d} \Delta \pmb {\phi} ^ {N} (\mathrm{nosummation}). \\ \end{array}
$$

Now we assume that the difference in incremental rotation along a beam element is small; i.e.,

$$
| \Delta \phi (\xi) - \Delta \phi^ {N} | \ll 1,
$$

which implies that either

$$
| \mathbf {n} (\xi) - \mathbf {n} ^ {N} | \ll 1 \quad \text {and} \quad | \Delta \phi (\xi) - \Delta \phi^ {N} | \ll 1
$$

or

$$
| \Delta \phi (\xi) | \ll 1 \quad \text { and } \quad | \Delta \phi^ {N} | \ll 1.
$$

In the first case we can use the approximations ${ \bf n } ( \xi ) \approx { \bf n } ^ { N } \approx \overline { { \bf n } }$ and $\Delta \phi ( \boldsymbol { \xi } ) \approx \Delta \phi ^ { N } \approx \Delta \overline { { \phi } }$ which, with the use of the interpolation functions for $\mathrm { d } \Delta \phi _ { : }$ , gives

$$
\begin{array}{l} \mathrm{d} \pmb {\phi} \approx g _ {N} (\xi) \left[ \frac {\sin \Delta \overline {{{\phi}}}}{\Delta \overline {{{\phi}}}} \mathrm{d} \Delta \pmb {\phi} ^ {N} + \frac {(1 - \cos \Delta \overline {{{\phi}}})}{\Delta \overline {{{\phi}}}} \overline {{{\mathbf {n}}}} \times \mathrm{d} \Delta \pmb {\phi} ^ {N} + (1 - \frac {\sin \Delta \overline {{{\phi}}}}{\Delta \overline {{{\phi}}}}) \overline {{{\mathbf {n}}}} \overline {{{\mathbf {n}}}} \cdot \mathrm{d} \Delta \pmb {\phi} ^ {N} \right], \\ \mathrm{d} \pmb {\phi} ^ {N} \approx \frac {\sin \Delta \overline {{\phi}}}{\Delta \overline {{\phi}}} \mathrm{d} \Delta \pmb {\phi} ^ {N} + \frac {(1 - \cos \Delta \overline {{\phi}})}{\Delta \overline {{\phi}}} \overline {{\mathbf {n}}} \times \mathrm{d} \Delta \pmb {\phi} ^ {N} + (1 - \frac {\sin \Delta \overline {{\phi}}}{\Delta \overline {{\phi}}}) \overline {{\mathbf {n}}} \overline {{\mathbf {n}}} \cdot \mathrm{d} \Delta \pmb {\phi} ^ {N}. \\ \end{array}
$$

In the second case it follows directly that

$$
\mathrm{d} \pmb {\phi} \approx g _ {N} (\xi) \mathrm{d} \Delta \pmb {\phi} ^ {N}, \mathrm{d} \pmb {\phi} ^ {N} \approx \mathrm{d} \Delta \pmb {\phi} ^ {N}.
$$

In either case

$$
\mathrm{d} \boldsymbol {\phi} (\boldsymbol {\xi}) \approx g _ {N} (\boldsymbol {\xi}) \mathrm{d} \boldsymbol {\phi} ^ {N}.
$$

This approximate relationship is used only in the creation of the Jacobian, so the approximation will at worst result in reduced speed of convergence.

<!-- source-page: 349 -->

# Elements

Once a nodal update vector $\mathrm { d } \phi ^ { N }$ is obtained, an exact update procedure is followed. This is achieved by a transformation into quaternions, use of exact quaternion update formula, and transformation of the results back into an incremental Euler rotation vector:

$$
\mathrm{d} \pmb {q} ^ {N} = \left(\mathrm{d} q ^ {N}, \mathrm{d} \pmb {q} ^ {N}\right) = \left(\cos \frac {\mathrm{d} \phi^ {N}}{2}, \frac {\mathrm{d} \pmb {\phi} ^ {N}}{\mathrm{d} \phi^ {N}} \sin \frac {\mathrm{d} \phi^ {N}}{2}\right), \qquad \mathrm{(nosummation)}
$$

$$
\Delta \pmb {q} ^ {N} = \left(\Delta q ^ {N}, \Delta \pmb {q} ^ {N}\right) = \left(\cos \frac {\Delta \phi^ {N}}{2}, \frac {\Delta \pmb {\phi} ^ {N}}{\Delta \phi^ {N}} \sin \frac {\Delta \phi^ {N}}{2}\right), \qquad \mathrm{(nosummation)}
$$

$$
\begin{array}{l} \Delta \pmb {q} _ {n e w} ^ {N} = \mathrm{d} \pmb {q} ^ {N} \Delta \pmb {q} ^ {N} \\ = \left(\mathrm{d} q ^ {N} \Delta q ^ {N} - \mathrm{d} \mathbf {q} ^ {N} \cdot \Delta \mathbf {q} ^ {N}, \mathrm{d} q ^ {N} \Delta \mathbf {q} ^ {N} + \mathrm{d} \mathbf {q} ^ {N} \Delta q ^ {N} + \mathrm{d} \mathbf {q} ^ {N} \times \Delta \mathbf {q} ^ {N}\right), \mathrm{(nosummation)} \\ \end{array}
$$

$$
\Delta \pmb {\phi} _ {n e w} ^ {N} = 2 \frac {\mathbf {q} _ {n e w} ^ {N}}{| \mathbf {q} _ {n e w} ^ {N} |} \tan^ {- 1} (| \mathbf {q} _ {n e w} ^ {N} |, q _ {n e w} ^ {N}) \qquad \mathrm{(nosummation).}
$$

The incremental rotation vector at the integration point is obtained by interpolation. Subsequently, we can calculate the updated integration point normals and the incremental curvature and twist.

# Second variations

For the calculation of the Jacobian, we also need the second variation in the generalized quantities. These follow from the first variations:

$$
\mathrm{d} \delta \lambda = \mathrm{d} \mathbf {s} \cdot \frac {d \delta \mathbf {x}}{d S} = \lambda^ {- 1} \frac {d \mathrm{d} \mathbf {x}}{d S} \cdot (\mathbf {I} - \mathbf {s s}) \cdot \frac {d \delta \mathbf {x}}{d S},
$$

$$
\begin{array}{l} \mathrm{d} \delta \gamma_ {\alpha} = \lambda^ {- 2} \gamma_ {\alpha} \mathrm{d} \lambda \delta \lambda - \lambda^ {- 1} \gamma_ {\alpha} \mathrm{d} \delta \lambda - \lambda^ {- 1} \delta \lambda \mathrm{d} \gamma_ {\alpha} - \lambda^ {- 2} \mathrm{d} \lambda \frac {d \delta \mathbf {x}}{d S} \cdot \mathbf {n} _ {\alpha} \\ + \lambda^ {- 1} \frac {d \delta \mathbf {x}}{d S} \cdot (\mathrm{d} \boldsymbol {\phi} \times \mathbf {n} _ {\alpha}) - \lambda^ {- 1} \mathrm{d} \lambda \mathbf {s} \cdot (\delta \boldsymbol {\phi} \times \mathbf {n} _ {\alpha}) + \lambda^ {- 1} \frac {d \mathrm{d} \mathbf {x}}{d S} \cdot (\delta \boldsymbol {\phi} \times \mathbf {n} _ {\alpha}) \\ + \frac {1}{2} \mathbf {s} \cdot \left[ \delta \boldsymbol {\phi} \times \left(\mathrm{d} \boldsymbol {\phi} \times \mathbf {n} _ {\alpha}\right) + \mathrm{d} \boldsymbol {\phi} \times \left(\delta \boldsymbol {\phi} \times \mathbf {n} _ {\alpha}\right) \right] \\ \approx - \lambda^ {- 1} (\delta \lambda \mathrm{d} \gamma_ {\alpha} + \mathrm{d} \lambda \delta \gamma_ {\alpha}) + \frac {1}{2} \delta \pmb {\phi} \cdot (\mathbf {n} _ {\alpha} \mathbf {t} + \mathbf {t n} _ {\alpha}) \cdot \mathrm{d} \pmb {\phi} \\ - \lambda^ {- 1} \epsilon_ {\alpha} ^ {\beta} (\delta \lambda \mathrm{d} \pmb {\phi} \cdot \mathbf {n} _ {\beta} + \mathrm{d} \lambda \delta \pmb {\phi} \cdot \mathbf {n} _ {\beta}) + \lambda^ {- 1} \epsilon_ {\alpha} ^ {\beta} (\frac {d \delta \mathbf {x}}{d S} \cdot \mathbf {n} _ {\beta} \mathbf {t} \cdot \mathrm{d} \pmb {\phi} + \frac {d \mathrm{d} \mathbf {x}}{d S} \cdot \mathbf {n} _ {\beta} \mathbf {t} \cdot \delta \pmb {\phi}) \\ \approx - \lambda^ {- 2} \frac {d \delta \mathbf {x}}{d S} \cdot (\mathbf {n} _ {\alpha} \mathbf {t} + \mathbf {t} \mathbf {n} _ {\alpha}) \cdot \frac {d \mathrm{d} \mathbf {x}}{d S} + \lambda^ {- 1} \epsilon_ {\alpha} ^ {\beta} \left(\frac {d \delta \mathbf {x}}{d S} \cdot \mathbf {n} _ {\beta} \mathbf {t} \cdot \mathrm{d} \pmb {\phi} + \delta \pmb {\phi} \cdot \mathbf {t} \mathbf {n} _ {\beta} \cdot \frac {d \mathrm{d} \mathbf {x}}{d S}\right) \\ + \frac {1}{2} \delta \pmb {\phi} \cdot (\mathbf {n} _ {\alpha} \mathbf {t} + \mathbf {t n} _ {\alpha}) \cdot \dot {\pmb {\phi}}, \\ \end{array}
$$

where we have again used $\gamma _ { \alpha } \approx 0$

<!-- source-page: 350 -->

For the second variation of the curvature we find

$$
\begin{array}{l} \mathrm{d} \delta b _ {\alpha} = \epsilon_ {\alpha} ^ {\beta} \left[ \mathrm{d} \delta \mathbf {t} \cdot \frac {d \mathbf {n} _ {\beta}}{d S} + \delta \mathbf {t} \cdot \frac {d \mathrm{d} \mathbf {n} _ {\beta}}{d S} + \mathrm{d} \mathbf {t} \cdot \frac {d \delta \mathbf {n} _ {\beta}}{d S} + \mathbf {t} \cdot \frac {d \mathrm{d} \delta \mathbf {n} _ {\beta}}{d S} \right] \\ = \epsilon_ {\alpha} ^ {\beta} \left[ \frac {1}{2} \delta \boldsymbol {\phi} \cdot \left(\mathbf {t} \frac {d \mathbf {n} _ {\beta}}{d S} + \frac {d \mathbf {n} _ {\beta}}{d S} \mathbf {t}\right) \cdot d \boldsymbol {\phi} - \left(\mathbf {t} \cdot \frac {d \mathbf {n} _ {\beta}}{d S}\right) (\delta \boldsymbol {\phi} \cdot d \boldsymbol {\phi}) \right. \\ + (\delta \boldsymbol {\phi} \times \mathbf {t}) \cdot \left(\frac {d \mathrm{d} \boldsymbol {\phi}}{d S} \times \mathbf {n} _ {\beta} + \mathrm{d} \boldsymbol {\phi} \times \frac {d \mathbf {n} _ {\beta}}{d S}\right) + (\mathrm{d} \boldsymbol {\phi} \times \mathbf {t}) \cdot \left(d \delta \boldsymbol {\phi} \times \mathbf {n} _ {\beta} + \delta \boldsymbol {\phi} \times \frac {d \mathbf {n} _ {\beta}}{d S}\right) \\ + \mathbf {t} \cdot \left(\frac {1}{2} \frac {d \mathrm{d} \boldsymbol {\phi}}{d S} \delta \boldsymbol {\phi} \cdot \mathbf {n} _ {\beta} + \frac {1}{2} \mathrm{d} \boldsymbol {\phi} \frac {d \delta \boldsymbol {\phi}}{d S} \cdot \mathbf {n} _ {\beta} + \frac {1}{2} \mathrm{d} \boldsymbol {\phi} \delta \boldsymbol {\phi} \cdot \frac {d \mathbf {n} _ {\beta}}{d S}\right) \\ + \mathbf {t} \cdot \left(\frac {1}{2} \frac {d \delta \boldsymbol {\phi}}{d S} \mathrm{d} \boldsymbol {\phi} \cdot \mathbf {n} _ {\beta} + \frac {1}{2} \delta \boldsymbol {\phi} \frac {d \mathrm{d} \boldsymbol {\phi}}{d S} \cdot \mathbf {n} _ {\beta} + \frac {1}{2} \delta \boldsymbol {\phi} \mathrm{d} \boldsymbol {\phi} \cdot \frac {d \mathbf {n} _ {\beta}}{d S}\right) \\ \left. - \mathbf {t} \cdot \left(\frac {d \delta \boldsymbol {\phi}}{d S} \cdot \mathrm{d} \boldsymbol {\phi} \mathbf {n} _ {\beta} + \delta \boldsymbol {\phi} \cdot \frac {d \mathrm{d} \boldsymbol {\phi}}{d S} \mathbf {n} _ {\beta} + \delta \boldsymbol {\phi} \cdot \mathrm{d} \boldsymbol {\phi} \frac {d \mathbf {n} _ {\beta}}{d S}\right) \right]. \\ \end{array}
$$

Rewriting the second and third terms and combining the others yields

$$
\begin{array}{l} \mathrm{d} \delta b _ {\alpha} = \epsilon_ {\alpha} ^ {\beta} \left[ \delta \pmb {\phi} \cdot \left(\mathbf {t} \frac {d \mathbf {n} _ {\beta}}{d S} + \frac {d \mathbf {n} _ {\beta}}{d S} \mathbf {t}\right) \cdot \mathrm{d} \pmb {\phi} - 2 \left(\mathbf {t} \cdot \frac {d \mathbf {n} _ {\beta}}{d S}\right) (\delta \pmb {\phi} \cdot \mathrm{d} \pmb {\phi}) \right. \\ + \frac {1}{2} \frac {d \delta \boldsymbol {\phi}}{d S} \cdot (\mathbf {n} _ {\beta} \mathbf {t} + \mathbf {t n} _ {\beta}) \cdot \mathrm{d} \boldsymbol {\phi} + \frac {1}{2} \delta \boldsymbol {\phi} \cdot (\mathbf {n} _ {\beta} \mathbf {t} + \mathbf {t n} _ {\beta}) \cdot \frac {d \mathrm{d} \boldsymbol {\phi}}{d S} \\ - \delta \boldsymbol {\phi} \cdot \mathbf {n} _ {\beta} \mathbf {t} \cdot \frac {d \mathrm{d} \boldsymbol {\phi}}{d S} + \left(\mathbf {t} \cdot \frac {d \mathbf {n} _ {\beta}}{d S}\right) (\delta \boldsymbol {\phi} \cdot \mathrm{d} \boldsymbol {\phi}) - \delta \boldsymbol {\phi} \cdot \frac {d \mathbf {n} _ {\beta}}{d S} \mathbf {t} \cdot \mathrm{d} \boldsymbol {\phi} \\ \left. - \mathrm{d} \boldsymbol {\phi} \cdot \mathbf {n} _ {\beta} \mathbf {t} \cdot \frac {d \delta \boldsymbol {\phi}}{d S} + \left(\mathbf {t} \cdot \frac {d \mathbf {n} _ {\beta}}{d S}\right) (\mathrm{d} \boldsymbol {\phi} \cdot \delta \boldsymbol {\phi}) - \mathrm{d} \boldsymbol {\phi} \cdot \frac {d \mathbf {n} _ {\beta}}{d S} \mathbf {t} \cdot \delta \boldsymbol {\phi} \right] \\ = \frac {1}{2} \epsilon_ {\alpha} ^ {\beta} \left[ \frac {d \delta \boldsymbol {\phi}}{d S} \cdot (\mathbf {n} _ {\beta} \mathbf {t} - \mathbf {t n} _ {\beta}) \cdot \mathrm{d} \boldsymbol {\phi} - \delta \boldsymbol {\phi} \cdot (\mathbf {n} _ {\beta} \mathbf {t} - \mathbf {t n} _ {\beta}) \cdot \frac {d \mathrm{d} \boldsymbol {\phi}}{d S} \right]. \\ \end{array}
$$

The second variation of twist is

$$
\begin{array}{l} \mathrm{d} \delta b = \mathrm{d} \delta \mathbf {n} _ {2} \cdot \frac {d \mathbf {n} _ {1}}{d S} + \delta \mathbf {n} _ {2} \cdot \frac {d \mathrm{d} \mathbf {n} _ {1}}{d S} + \mathrm{d} \mathbf {n} _ {2} \cdot \frac {d \delta \mathbf {n} _ {1}}{d S} + \mathbf {n} _ {2} \cdot \frac {d \mathrm{d} \delta \mathbf {n} _ {1}}{d S} \\ = \frac {1}{2} \delta \boldsymbol {\phi} \cdot \left(\mathbf {n} _ {2} \frac {d \mathbf {n} _ {1}}{d S} + \frac {d \mathbf {n} _ {1}}{d S} \mathbf {n} _ {2}\right) \cdot d \boldsymbol {\phi} - \left(\mathbf {n} _ {2} \cdot \frac {d \mathbf {n} _ {1}}{d S}\right) (\delta \boldsymbol {\phi} \cdot d \boldsymbol {\phi}) \\ + \left(\delta \boldsymbol {\phi} \times \mathbf {n} _ {2}\right) \cdot \left(\frac {d \mathrm{d} \boldsymbol {\phi}}{d S} \times \mathbf {n} _ {1} + \mathrm{d} \boldsymbol {\phi} \times \frac {d \mathbf {n} _ {1}}{d S}\right) + \left(\mathrm{d} \boldsymbol {\phi} \times \mathbf {n} _ {2}\right) \cdot \left(\frac {d \delta \boldsymbol {\phi}}{d S} \times \mathbf {n} _ {1} + \delta \boldsymbol {\phi} \times \frac {d \mathbf {n} _ {1}}{d S}\right) \\ + \mathbf {n} _ {2} \cdot \left(\frac {1}{2} \frac {d \mathrm{d} \boldsymbol {\phi}}{d S} \delta \boldsymbol {\phi} \cdot \mathbf {n} _ {1} + \frac {1}{2} \mathrm{d} \boldsymbol {\phi} \frac {d \delta \boldsymbol {\phi}}{d S} \cdot \mathbf {n} _ {1} + \frac {1}{2} \mathrm{d} \boldsymbol {\phi} \delta \boldsymbol {\phi} \cdot \frac {d \mathbf {n} _ {1}}{d S}\right) \\ + \mathbf {n} _ {2} \cdot \left(\frac {1}{2} \frac {d \text {delta} \boldsymbol {\phi}}{d S} \mathrm{d} \boldsymbol {\phi} \cdot \mathbf {n} _ {1} + \frac {1}{2} \delta \boldsymbol {\phi} \frac {d \mathrm{d} \boldsymbol {\phi}}{d S} \cdot \mathbf {n} _ {1} + \frac {1}{2} \delta \boldsymbol {\phi} \mathrm{d} \boldsymbol {\phi} \cdot \frac {d \mathbf {n} _ {1}}{d S}\right) \\ - \mathbf {n} _ {2} \cdot \left(\frac {d \delta \boldsymbol {\phi}}{d S} \cdot \mathrm{d} \boldsymbol {\phi} \mathbf {n} _ {1} + \delta \boldsymbol {\phi} \cdot \frac {d \mathrm{d} \boldsymbol {\phi}}{d S} \mathbf {n} _ {1} + \delta \boldsymbol {\phi} \cdot \mathrm{d} \boldsymbol {\phi} \frac {d \mathbf {n} _ {1}}{d S}\right). \\ \end{array}
$$