# Elements The iterative solution procedure requires the calculation of the consistent tangent stiffness. This stiffness has two parts, one resulting from the material model and the other resulting from the changingresulting from geometry. We denote the second variational quantities with d. Returning to Equation 3.6.7-4, the virtual work expression can be written as $$ \delta \Pi = \int_ {V} J \sigma^ {a b} \delta E _ {a b} d V ^ {0} = \int_ {A} \int_ {h} J \sigma^ {a b} \delta E _ {a b} j ^ {0} d S _ {3} d S d \theta , $$ where $J = j / j ^ { 0 } , j ^ { 0 }$ is the Jacobian determinant of the reference configuration parametrization given by $\boldsymbol { j } ^ { 0 } = \operatorname* { d e t } [ \boldsymbol { \nabla } \mathbf { X } ]$ , and components are relative to the corotated frame with $a , b \in \{ 1 , 2 \}$ . We assume constitutive behavior such that $$ \mathrm{d} (J \sigma^ {a b}) = J C ^ {a b c d} \mathrm{d} E _ {c d}. $$ Thus, the variation of the virtual work expression yields $$ \begin{array}{l} \mathrm{d} (\delta \Pi) = \int_ {A} \Big [ \int_ {h} \frac {j}{\overline {{j}}} \big (\delta \epsilon_ {a b} + \overline {{f}} _ {3 3} S _ {3} \delta \kappa_ {a b} \big) C ^ {a b c d} \big (\mathrm{d} \epsilon_ {c d} + \overline {{f}} _ {3 3} S _ {3} \mathrm{d} \kappa_ {c d} \big) d S _ {3} \\ \left. + N ^ {a b} \mathrm{d} (\delta \epsilon_ {a b}) + M ^ {a b} \mathrm{d} (\delta \kappa_ {a b}) \right] \overline {{j}} d S d \theta . \\ \end{array} $$ The first term in the integrand forms the material stiffness, while the second two terms form the geometric stiffness. The second variation of the strain measure components are calculated in \`\`Finite-strain shell element formulation,'' Section 3.6.5. The second variation of the membrane strain is $$ \mathrm{d} (\delta \epsilon_ {a b}) = \mathrm{sym} \Big [ \frac {\partial (\delta \overline {{{\bf x}}})}{\partial s ^ {a}} \cdot \frac {\partial (\mathrm{d} \overline {{{\bf x}}})}{\partial s ^ {b}} - 2 \delta \epsilon_ {a c} \mathrm{d} \epsilon_ {c b} \Big ]. $$ The second variation of the bending strain is zero; i.e., $$ \mathrm{d} (\delta \kappa_ {a b}) = 0. $$ # Incremental degrees of freedom: interpolations and configuration updates At the beginning of an increment we have the configuration at iteration $k ,$ denoted $\left( \overline { { \mathbf { x } } } ^ { k } , \mathbf { t } _ { 3 } ^ { k } \right)$ . In the incremental solution procedure we solve for the incremental displacement field and the incremental rotation field in components relative to the reference cylindrical coordinate system: $$ (\Delta u _ {r} ^ {m p}, \Delta u _ {z} ^ {m p}, \Delta u _ {\theta} ^ {m p}, \Delta \phi_ {r} ^ {m p}, \Delta \phi_ {z} ^ {m p}, \Delta \phi_ {\theta} ^ {m p}). $$ We use these incremental fields to update the configuration to iteration $k + 1$ . 1. Reference surface update: The displacement increments are interpolated using the same # Elements interpolation scheme as the total displacement vector in Equation 3.6.7-2: $$ \Delta \mathbf {u} (S, \theta) = \left\{ \begin{array}{l} \Delta u _ {r} (S, \theta) \\ \Delta u _ {z} (S, \theta) \\ \Delta u _ {\theta} (S, \theta) \end{array} \right\} = \sum_ {m = 1} ^ {M} H ^ {m} (S) \left(\sum_ {p = 1} ^ {P + 1} R ^ {p} (\theta) \left\{ \begin{array}{c} \Delta u _ {r} ^ {m p} \\ \Delta u _ {z} ^ {m p} \\ 0 \end{array} \right\} + \sum_ {p = 1} ^ {P} \sin (p \theta) \left\{ \begin{array}{c} 0 \\ 0 \\ \Delta u _ {\theta} ^ {m p} \end{array} \right\}\right). $$ The reference surface position map is updated from the displacement increment by $$ \overline {{{\mathbf {x}}}} ^ {k + 1} (S, \theta) = \overline {{{\mathbf {x}}}} ^ {k} (S, \theta) + \Delta \mathbf {u} (S, \theta). $$ 2. Rotation field update: The incremental rotation field is updated with the same interpolation scheme as the total rotation field in Equation 3.6.7-3: $$ \Delta \pmb {\phi} (S, \theta) = \left\{ \begin{array}{l} \Delta \phi_ {r} (S, \theta) \\ \Delta \phi_ {z} (S, \theta) \\ \Delta \phi_ {\theta} (S, \theta) \end{array} \right\} = \sum_ {m = 1} ^ {M} H ^ {m} (S) \left(\sum_ {p = 1} ^ {P + 1} R ^ {p} (\theta) \left\{ \begin{array}{c} 0 \\ 0 \\ - \Delta \phi_ {\theta} ^ {m p} \end{array} \right\} + \sum_ {p = 1} ^ {P} \sin (p \theta) \left\{ \begin{array}{c} \Delta \phi_ {r} ^ {m p} \\ \Delta \phi_ {z} ^ {m p} \\ 0 \end{array} \right\}\right). $$ This incremental rotation vector corresponds to a finite rotation, characterized by quaternion parameters as $$ \Delta q _ {0} = \cos \left\| \Delta \phi / 2 \right\| \quad \mathrm{and} \quad \Delta \mathbf {q} = \frac {\sin \left\| \Delta \phi / 2 \right\|}{\left\| \Delta \phi \right\|} \Delta \phi . $$ The total rotation quaternion parameter can be updated by the update formula $$ \boldsymbol {q} ^ {k + 1} = \left(q _ {0} ^ {k + 1}, \mathbf {q} ^ {k + 1}\right) = \left(\Delta q _ {0} q _ {0} ^ {k} - \Delta \mathbf {q} \cdot \mathbf {q} ^ {k}, \Delta q _ {0} \mathbf {q} ^ {k} + q _ {0} ^ {k} \Delta \mathbf {q} + \Delta \mathbf {q} \times \mathbf {q} ^ {k}\right). $$ Similarly, the deformed unit director field can be updated from the incremental rotation field as $$ \mathbf {t} _ {3} ^ {k + 1} = (2 (\Delta q _ {0}) ^ {2} - 1) \mathbf {t} _ {3} ^ {k} + 2 \Delta q _ {0} \Delta \mathbf {q} \times \mathbf {t} _ {3} ^ {k} + 2 (\Delta \mathbf {q} \cdot \mathbf {t} _ {3} ^ {k}) \Delta \mathbf {q} \stackrel {\mathrm{def}} {=} (\Delta q _ {0}, \Delta \mathbf {q}) \circ \mathbf {t} _ {3} ^ {k}. $$ Here we have used the notation $( q _ { 0 } , \mathbf { q } ) \circ \mathbf { v }$ to denote the rotation of the vector v by the quaternion $( q _ { 0 } , \mathbf { q } )$ . The curvatures are calculated from the gradient of the director field, which is updated by $$ \mathbf {t} _ {3, \alpha} ^ {k + 1} = (\Delta q _ {0}, \Delta \mathbf {q}) \circ \mathbf {t} _ {3, \alpha} ^ {k} + \Delta \mathbf {R} _ {\alpha} \times \mathbf {t} _ {3} ^ {k + 1}, $$ where $$ \Delta \mathbf {R} _ {\alpha} = \mathbf {H} (\Delta \boldsymbol {\phi}) \cdot \Delta \boldsymbol {\phi} _ {, \alpha}, $$ $$ \mathbf {H} (\Delta \phi) = \frac {\sin \left\| \Delta \phi \right\|}{\left\| \Delta \phi \right\|} \mathbf {I} + \frac {1}{\left\| \Delta \phi \right\| ^ {2}} \left(1 - \frac {\sin \left\| \Delta \phi \right\|}{\left\| \Delta \phi \right\|}\right) \Delta \phi \Delta \phi + \frac {\left(1 - \cos \left\| \Delta \phi \right\|\right)}{\left\| \Delta \phi \right\| ^ {2}} \widehat {\Delta \phi}. $$ For completeness, we record the values of $\Delta \phi _ { , \alpha }$ . First, along the generator lines we have # Elements $$ \Delta \phi_ {, S} = \sum_ {m = 1} ^ {M} H _ {, S} ^ {m} (S) \left(\sum_ {p = 1} ^ {P + 1} R ^ {p} (\theta) \left\{ \begin{array}{c} 0 \\ 0 \\ \Delta \phi_ {\theta} ^ {m p} \end{array} \right\} + \sum_ {p = 1} ^ {P} \sin (p \theta) \left\{ \begin{array}{c} - \Delta \phi_ {r} ^ {m p} \\ \Delta \phi_ {z} ^ {m p} \\ 0 \end{array} \right\}\right) . $$ For the circumferential derivative we must account for the derivative of the basis vectors in the µ-direction: $\mathbf { e } _ { r , \theta } = \mathbf { e } _ { \theta }$ and $\mathbf { e } _ { \theta , \theta } = - \mathbf { e } _ { r }$ . Hence, $\Delta \phi _ { , \theta } = ( \Delta \phi _ { r , \theta } - \Delta \phi _ { \theta } ) \mathbf { e } _ { r } + \Delta \phi _ { z , \theta } \mathbf { e } _ { z } + ( \Delta \phi _ { \theta , \theta } + \Delta \phi _ { r } ) \mathbf { e } _ { \theta }$ . Introducing the interpolation function, we have $$ \Delta \pmb {\phi} _ {, \theta} = \left\{ \begin{array}{c} - \Delta \phi_ {\theta} \\ 0 \\ \Delta \phi_ {r} \end{array} \right\} + \sum_ {m = 1} ^ {M} H ^ {m} (S) \left(\sum_ {p = 1} ^ {P + 1} R _ {, \theta} ^ {p} (\theta) \left\{ \begin{array}{c} 0 \\ 0 \\ - \Delta \phi_ {\theta} ^ {m p} \end{array} \right\} + \sum_ {p = 1} ^ {P} p \cos (p \theta) \left\{ \begin{array}{c} \Delta \phi_ {r} ^ {m p} \\ \Delta \phi_ {z} ^ {m p} \\ 0 \end{array} \right\}\right), $$ where $R _ { , \theta } ^ { p } ( \theta )$ are computed from the definitions of the interpolation functions $R ^ { p } ( \theta )$ as given above and $\Delta \phi _ { \theta }$ and $\Delta \phi _ { \mathrm { i } }$ r are given by the interpolation for the incremental rotation field as given in the Rotation Field Update. # Strain increment and stress resultant update Following the formulation of the finite-strain shell element formulation in \`\`Finite-strain shell element formulation,'' Section 3.6.5, the three-dimensional (finite) strain increment is calculated as $$ \Delta E _ {a b} = \Delta \epsilon_ {a b} + S _ {3} \overline {{f}} _ {3 3} ^ {t + \Delta t} \Delta \kappa_ {a b}. $$ Here the increment in the membrane strain components $\Delta \epsilon _ { a b }$ is given by the Hughes-Winget second-order approximation to the logarithmic membrane strain increment $$ \Delta \epsilon_ {a b} = 2 (\Delta \overline {{f}} _ {a c} - \delta_ {a c}) (\Delta \overline {{f}} _ {c b} + \delta_ {c b}) ^ {- 1}, $$ and the increment in the bending strain components $\Delta \kappa _ { a b }$ is given by the expression $$ \Delta \kappa_ {a b} = \mathrm{sym} \left[ e _ {a} ^ {c} \mathbf {n} _ {c} ^ {t + \Delta t} \cdot \Delta \mathbf {r} _ {b} \right], $$ where $$ \Delta \mathbf {r} _ {b} = \overline {{h}} _ {b} ^ {\alpha} \Delta \mathbf {R} _ {\alpha}. $$ As an example of the stress update procedure, consider the simple case of a Saint Venant-Kirchhoff material model. In this case, $$ \sigma^ {(t + \Delta t) a b} = \Delta \overline {{f}} _ {c} ^ {a} \sigma^ {(t) c d} \Delta \overline {{f}} _ {d} ^ {b} + \Delta \overline {{f}} _ {e} ^ {a} \Delta \overline {{f}} _ {f} ^ {b} C ^ {(t) e f g h} \Delta \overline {{f}} _ {g} ^ {c} \Delta \overline {{f}} _ {h} ^ {d} \Delta E _ {c d}, $$ where $C ^ { ( 0 ) e f g h }$ are the plane stress elastic coefficients given by # Elements $$ C ^ {(0) e f g h} = \frac {E}{(1 - \nu^ {2})} \left\{\nu \delta^ {e f} \delta^ {g h} + \frac {(1 - \nu)}{2} \big (\delta^ {e g} \delta^ {f h} + \delta^ {e h} \delta^ {f g} \big) \right\}. $$ Note that $a , b , c , d , e , f , h \in \{ 1 , 2 \}$ , and the components $\sigma ^ { ( t + \Delta t ) a b }$ are relative to the current orthonormal basis $\left\{ _ 1 ^ { t + \Delta t } , \mathbf n _ { 2 } ^ { t + \Delta t } \right\}$ f1 . # Pressure loads and load stiffness For geometrically linear problems, equivalent nodal loads due to applied surface pressure are readily calculated since the geometry is axisymmetric. For geometrically nonlinear problems, however, asymmetric deformations must be taken into consideration. The equivalent nodal loads associated with surface pressure p can be obtained by considering the virtual work contribution to the external loading $$ - \int_ {A} p \mathbf {n} \cdot \delta \overline {{\mathbf {x}}} d A = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p \left(\frac {\partial \overline {{\mathbf {x}}}}{\partial \theta} \times \frac {\partial \overline {{\mathbf {x}}}}{\partial S}\right) \cdot \delta \overline {{\mathbf {x}}} d S d \theta , $$ where S is the parametric surface coordinate in the R-Z plane and the reference surface position is $$ \overline {{\mathbf {x}}} = r \mathbf {e} _ {r} + z \mathbf {e} _ {z} + u _ {\theta} \mathbf {e} _ {\theta}, $$ with $r = R + u _ { r }$ and $z = Z + u _ { z }$ . Recall that the current position of a point can be expressed in terms of the axial interpolator $H ^ { m } ( S )$ and the circumferential interpolators $R ^ { p } ( \theta )$ by Equation 3.6.7-7 $$ \left\{ \begin{array}{l} r \\ z \\ u _ {\theta} \end{array} \right\} = \sum_ {m = 1} ^ {M} H ^ {m} (S) \left(\sum_ {p = 1} ^ {P + 1} R ^ {p} (\theta) \left\{ \begin{array}{c} r ^ {m p} \\ z ^ {m p} \\ 0 \end{array} \right\} + \sum_ {p = 1} ^ {P} \sin p \theta \left\{ \begin{array}{c} 0 \\ 0 \\ u _ {\theta} ^ {m p} \end{array} \right\}\right). $$ The terms in Equation 3.6.7-7 can be worked out as follows: $$ \frac {\partial \overline {{\mathbf {x}}}}{\partial \theta} = \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) \mathbf {e _ {r}} + \frac {\partial z}{\partial \theta} \mathbf {e} _ {z} + \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) \mathbf {e} _ {\theta}, $$ $$ \frac {\partial \overline {{\mathbf {x}}}}{\partial S} = \frac {\partial r}{\partial S} \mathbf {e} _ {r} + \frac {\partial z}{\partial S} \mathbf {e} _ {z} + \frac {\partial u _ {\theta}}{\partial S} \mathbf {e} _ {\theta}, $$ and, hence, $$ \begin{array}{l} \frac {\partial \overline {{\mathbf {x}}}}{\partial \theta} \times \frac {\partial \overline {{\mathbf {x}}}}{\partial S} = \left[ \frac {\partial z}{\partial S} \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) - \frac {\partial u _ {\theta}}{\partial S} \frac {\partial z}{\partial \theta} \right] \mathbf {e} _ {r} \\ + \left[ \frac {\partial u _ {\theta}}{\partial S} \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) - \frac {\partial r}{\partial S} \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) \right] \mathbf {e} _ {z} \\ + \left[ \frac {\partial r}{\partial S} \frac {\partial z}{\partial \theta} - \frac {\partial z}{\partial S} \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) \right] \mathbf {e} _ {\theta}. \\ \end{array} $$ # Elements The variations $\delta \mathbf { \overline { { x } } }$ are written $\delta \mathbf { \overline { { x } } } = \delta \overline { { x } } _ { r } \mathbf { e } _ { r } + \delta \overline { { x } } _ { z } \mathbf { e } _ { z } + \delta \overline { { u } } _ { \theta } \mathbf { e } _ { \theta }$ , where the components have interpolations similar to those in Equation 3.6.7-7. Therefore, the virtual work contribution becomes $$ \begin{array}{l} - \int_ {A} p \mathbf {n} \cdot \delta \overline {{{\mathbf {x}}}} d A = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p \left\{\left[ \frac {\partial z}{\partial S} \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) - \frac {\partial u _ {\theta}}{\partial S} \frac {\partial z}{\partial \theta} \right] \delta \overline {{{x}}} _ {r} \right. \\ + \left[ \frac {\partial u _ {\theta}}{\partial S} \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) - \frac {\partial r}{\partial S} \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) \right] \delta \overline {{x}} _ {z} \\ \left. + \left[ \frac {\partial r}{\partial S} \frac {\partial z}{\partial \theta} - \frac {\partial z}{\partial S} \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) \right] \delta u _ {\theta} \right\} d S d \theta . \\ \end{array} $$ With the introduction of the interpolation functions, we obtain the equivalent nodal forces: $$ \begin{array}{l} F _ {r} ^ {m p} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p H ^ {m} (S) R ^ {p} (\theta) \left[ \frac {\partial z}{\partial S} \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) - \frac {\partial u _ {\theta}}{\partial S} \frac {\partial z}{\partial \theta} \right] d S d \theta \\ F _ {z} ^ {m p} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p H ^ {m} (S) R ^ {p} (\theta) \left[ \frac {\partial u _ {\theta}}{\partial S} x \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) - \frac {\partial r}{\partial S} \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) \right] d S d \theta \\ F _ {\theta} ^ {m p} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p H ^ {m} (S) \sin p \theta \left[ \frac {\partial r}{\partial S} \frac {\partial z}{\partial \theta} - \frac {\partial z}{\partial S} \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) \right] d S d \theta . \\ \end{array} $$ For geometrically linear analysis, the equivalent nodal forces reduce to the standard axisymmetric expressions $$ \begin{array}{l} F _ {r} ^ {m p} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p H ^ {m} (S) R ^ {p} (\theta) \frac {\partial Z}{\partial S} R d S d \theta \\ F _ {z} ^ {m p} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} - p H ^ {m} (S) R ^ {p} (\theta) \frac {\partial R}{\partial S} R d S d \theta \\ F _ {\theta} ^ {m p} = 0. \\ \end{array} $$ The linearization of the pressure loading term leads to the following pressure load stiffness matrix: $$ \left\{ \begin{array}{l} \mathrm{d} F _ {r} ^ {m p} \\ \mathrm{d} F _ {z} ^ {m p} \\ \mathrm{d} F _ {\theta} ^ {m p} \end{array} \right\} = \left[ \begin{array}{l l l} K _ {r r} ^ {m p n q} & K _ {r z} ^ {m p n q} & K _ {r \theta} ^ {m p n q} \\ K _ {z r} ^ {m p n q} & K _ {z z} ^ {m p n q} & K _ {z \theta} ^ {m p n q} \\ K _ {\theta r} ^ {m p n q} & K _ {\theta z} ^ {m p n q} & K _ {\theta \theta} ^ {m p n q} \end{array} \right] \left\{ \begin{array}{l} \mathrm{d} \overline {{x}} _ {r} ^ {n q} \\ \mathrm{d} \overline {{x}} _ {z} ^ {n q} \\ \mathrm{d} u _ {\theta} ^ {n q} \end{array} \right\}, $$ with $$ \begin{array}{l} K _ {r r} ^ {m p n q} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p H ^ {m} (S) R ^ {p} (\theta) \frac {\partial z}{\partial S} H ^ {n} (S) R ^ {q} (\theta) d S d \theta \\ K _ {z r} ^ {m p n q} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p H ^ {m} (S) R ^ {p} (\theta) \left[ - \frac {\partial r}{\partial S} H ^ {n} (S) R ^ {q} (\theta) - \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) \frac {d H ^ {n}}{d S} R ^ {q} (\theta) \right. \\ \left. + \frac {\partial u _ {\theta}}{\partial S} H ^ {n} (S) \frac {d R ^ {q}}{d \theta} \right] d S d \theta \\ \end{array} $$ # Elements $$ K _ {\theta r} ^ {m p n q} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p H ^ {m} (S) \sin p \theta \left[ \frac {\partial z}{\partial \theta} \frac {d H ^ {n}}{d S} R ^ {q} (\theta) - \frac {\partial z}{\partial S} H ^ {n} (S) \frac {d R ^ {q}}{d \theta} \right] d S d \theta $$ $$ K _ {r z} ^ {m p n q} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p H ^ {m} (S) R ^ {p} (\theta) \left[ \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) \frac {d H ^ {n}}{d S} R ^ {q} (\theta) - \frac {\partial u _ {\theta}}{\partial S} H ^ {n} (S) \frac {d R ^ {q}}{d \theta} \right] d S d \theta $$ $$ K _ {z z} ^ {m p n q} = 0 $$ $$ K _ {\theta z} ^ {m p n q} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p H ^ {m} (S) \sin p \theta \left[ \frac {\partial r}{\partial S} H ^ {n} (S) \frac {d R ^ {q}}{d \theta} - \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) \frac {d H ^ {n}}{d S} R ^ {q} (\theta) \right] d S d \theta $$ $$ K _ {r \theta} ^ {m p n q} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p H ^ {m} (S) R ^ {p} (\theta) \left[ q \frac {\partial z}{\partial S} H ^ {n} (S) \cos q \theta - \frac {\partial z}{\partial \theta} \frac {d H ^ {n}}{d S} \sin q \theta \right] d S d \theta $$ $$ \begin{array}{l} K _ {z \theta} ^ {m p n q} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p H ^ {m} (S) R ^ {p} (\theta) \left[ - \frac {\partial u _ {\theta}}{\partial S} H ^ {n} (S) \sin q \theta + \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) \frac {d H ^ {n}}{d S} \sin q \theta \right. \\ \left. - q \frac {\partial r}{\partial S} H ^ {n} (S) \cos q \theta \right] d S d \theta \\ \end{array} $$ $$ K _ {\theta \theta} ^ {m p n q} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} p H ^ {m} (S) \sin p \theta \frac {\partial z}{\partial S} H ^ {n} (S) \sin q \theta d S d \theta . $$ In the case of hydrostatic pressure (p dependent on z), additional terms must be included in the pressure load stiffness. These terms appear due to the variation of the pressure magnitude and are readily obtained from the expression $$ - \int_ {A} \mathrm{d} p \mathbf {n} \cdot \delta \mathbf {u} d A = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} \mathrm{d} u _ {z} \frac {d p}{d z} \Bigg \{\left[ \frac {\partial z}{\partial S} \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) - \frac {\partial u _ {\theta}}{\partial S} \frac {\partial z}{\partial \theta} \right] \delta \overline {{x}} _ {r} + $$ $$ \left[ \frac {\partial u _ {\theta}}{\partial S} \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) - \frac {\partial r}{\partial S} \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) \right] \delta \overline {{x}} _ {z} + $$ $$ \left. \left[ \frac {\partial r}{\partial S} \frac {\partial z}{\partial \theta} - \frac {\partial z}{\partial S} \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) \right] \delta u _ {\theta} \right\} d S d \theta . $$ With use of the interpolation functions and denoting the additional contributions with an over-bar, we obtain the additional load stiffness contributions: $$ \overline {{K}} _ {r z} ^ {m p n q} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} H ^ {m} (S) R ^ {p} (\theta) \frac {d p}{d z} \left[ \frac {\partial z}{\partial S} \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) - \frac {\partial u _ {\theta}}{\partial S} \frac {\partial z}{\partial \theta} \right] H ^ {n} (S) R ^ {q} (\theta) d S d \theta $$ $$ \overline {{K}} _ {z z} ^ {m p n q} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} H ^ {m} (S) R ^ {p} (\theta) \frac {d p}{d z} \left[ \frac {\partial u _ {\theta}}{\partial S} \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) - \frac {\partial r}{\partial S} \left(r + \frac {\partial u _ {\theta}}{\partial \theta}\right) \right] H ^ {n} (S) R ^ {q} (\theta) d S d \theta $$ $$ \overline {{K}} _ {\theta z} ^ {m p n q} = \int_ {0} ^ {2 \pi} \int_ {- 1} ^ {+ 1} H ^ {m} (S) \sin p \theta \frac {d p}{d z} \left[ \frac {\partial r}{\partial S} \frac {\partial z}{\partial \theta} - \frac {\partial z}{\partial S} \left(\frac {\partial r}{\partial \theta} - u _ {\theta}\right) \right] H ^ {n} (S) R ^ {q} (\theta) d S d \theta . $$ # Penalty constraints: transverse shear and drill rotation It is necessary to enforce rotation constraints at selected points on the element surface to prevent singular modes of deformation. One axial transverse shear constraint is enforced on the $\phi _ { \theta }$ rotation field between each pair of nodes within each nodal plane. One circumferential transverse shear constraint and one drill rotation constraint on the rotation fields $\phi _ { r }$ and $\phi _ { z }$ are enforced between each pair of nodes on a circumferential line of nodes. In each instance the rotation field is constrained to follow the nodal displacements. To summarize, for element SAXAMN: - Axial transverse shear: $M ( N + 1 )$ （号 - Circumferential transverse shear: $( M + 1 ) N$ - Drill rotation: $( M + 1 ) N$ constraints are enforced, where $M \in \{ 1 , 2 \}$ is the order of integration in the axial direction and $N \in \{ 1 , 2 , 3 , 4 \}$ is the number of Fourier modes. # Transverse shear The transverse shear strain is the measure of the amount the director field $\mathbf { t } _ { 3 }$ has rotated relative to the normal to the shell surface. We define the transverse shear strain as $$ \gamma_ {3 c} \stackrel {\mathrm{def}} {=} \mathbf {t} _ {3} \cdot \overline {{\mathbf {t}}} _ {c}, \quad \text {where} \quad \overline {{\mathbf {t}}} _ {c} \stackrel {\mathrm{def}} {=} \overline {{\mathbf {x}}} _ {, \xi^ {c}} / \| \overline {{\mathbf {x}}} _ {, \xi^ {c}} \|, $$ with no sum on c and force this quantity to be zero with a penalty constraint. Note that $\overline { { \mathbf { t } } } _ { c }$ is a unit vector tangent to the shell surface defined by the displacement field along a parametric coordinate line. Recall that $\xi ^ { 1 } = S$ and $\xi ^ { 2 } = \theta$ . Hence, °31 is the axial transverse shear strain and $\gamma _ { 3 2 }$ is the circumferential transverse shear strain. For convenience, record the variation of the unit vector $\overline { { \mathbf { t } } } _ { c }$ : $$ \delta \overline {{\mathbf {t}}} _ {c} = \frac {1}{\| \overline {{\mathbf {x}}} , \xi^ {c} \|} \delta \overline {{\mathbf {x}}} _ {, \xi^ {c}} - \frac {\overline {{\mathbf {x}}} _ {, \xi^ {c}} \cdot \delta \overline {{\mathbf {x}}} _ {, \xi^ {c}}}{\| \overline {{\mathbf {x}}} , \xi^ {c} \| ^ {3}} \overline {{\mathbf {x}}} _ {, \xi^ {c}}, $$ with no sum on c. Equivalently, the definition of $\overline { { \mathbf { t } } } _ { c }$ can be used to write $$ \delta \overline {{\mathbf {t}}} _ {c} = \frac {1}{\| \overline {{\mathbf {x}}} , \xi^ {c} \|} \left[ \delta \overline {{\mathbf {x}}} _ {, \xi^ {c}} - (\overline {{\mathbf {t}}} _ {c} \cdot \delta \overline {{\mathbf {x}}} _ {, \xi^ {c}}) \mathbf {t} _ {c} \right]. $$ The linearized transverse shear strain is calculated as $$ \delta \gamma_ {3 c} = \delta \mathbf {t} _ {3} \cdot \overline {{\mathbf {t}}} _ {c} + \mathbf {t} _ {3} \cdot \delta \overline {{\mathbf {t}}} _ {c}. $$ Since by definition $\delta \mathbf { t } _ { 3 } = \delta { \pmb { \phi } } \times \mathbf { t } _ { 3 }$ , it follows that # Elements $$ \delta \gamma_ {3 c} = \delta \pmb {\phi} \cdot (\mathbf {t} _ {3} \times \overline {{\mathbf {t}}} _ {c}) + \frac {1}{\| \overline {{\mathbf {x}}} , \xi^ {c} \|} (\mathbf {t} _ {3} - \gamma_ {3 c} \overline {{\mathbf {t}}} _ {c}) \cdot \delta \overline {{\mathbf {x}}} _ {, \xi^ {c}}, $$ with no sum on c. For completeness, the second variation of the transverse shear strain is $$ \begin{array}{l} \mathrm{d} (\delta \gamma_ {3 c}) = \delta \pmb {\phi} \cdot \frac {1}{2} [ \mathbf {t} _ {3} \overline {{\mathbf {t}}} _ {c} + \overline {{\mathbf {t}}} _ {c} \mathbf {t} _ {3} ] \cdot \mathrm{d} \pmb {\phi} - \delta \overline {{\mathbf {x}}} _ {, \xi^ {c}} \cdot \frac {1}{\| \overline {{\mathbf {x}}} _ {, \xi^ {c}} \| ^ {2}} [ \mathbf {t} _ {3} \overline {{\mathbf {t}}} _ {c} + \overline {{\mathbf {t}}} _ {c} \mathbf {t} _ {3} ] \cdot \mathrm{d} \overline {{\mathbf {x}}} _ {, \xi^ {c}} \\ + \delta \pmb {\phi} \cdot \frac {1}{\| \overline {{{\mathbf {x}}}} _ {, \xi^ {c}} \|} [ \hat {\mathbf {t}} _ {3} - (\mathbf {t} _ {3} \times \overline {{{\mathbf {t}}}} _ {c}) \overline {{{\mathbf {t}}}} _ {c} ] \cdot \mathrm{d} \overline {{{\mathbf {x}}}} _ {, \xi^ {c}} + \delta \overline {{{\mathbf {x}}}} _ {, \xi^ {c}} \cdot \frac {1}{\| \overline {{{\mathbf {x}}}} _ {, \xi^ {c}} \|} [ - \hat {\mathbf {t}} _ {3} - \overline {{{\mathbf {t}}}} _ {c} (\mathbf {t} _ {3} \times \overline {{{\mathbf {t}}}} _ {c}) ] \cdot \mathrm{d} \pmb {\phi}, \\ \end{array} $$ where terms proportional to $\gamma _ { 3 c }$ have been neglected and the coupling between $\delta \phi$ and $\Delta \phi$ has been symmetrized. # Drill rotation Mathematically, the equations governing the deformation of the shell are invariant with respect to drill rotations; that is, a rotation of the director field with axis of rotation parallel to the director. It is necessary then to assign a kinematic definition to this rotation. We define the drill strain as the difference between the rotation of the circumferential tangent vector as measured by the displacement field and that measured by the rotation field. Accordingly, we define the drill strain by $$ \gamma_ {d} = \mathbf {a} \cdot \mathbf {b}. $$ Here, a is the rotated reference axial tangent vector and b is the deformed (unit) circumferential tangent vector, defined by $$ \mathbf {a} = \exp [ \hat {\boldsymbol {\phi}} ] \mathbf {A} \quad \text {and} \quad \mathbf {b} = \frac {\overline {{\mathbf {x}}} _ {, \theta}}{\| \overline {{\mathbf {x}}} _ {, \theta} \|}. $$ In the above A is a linear approximation to the axial tangent vector in the reference configuration given by $\mathbf { A } _ { 1 } = ( \mathbf { X } ^ { K } - \mathbf { X } ^ { N } ) / \| \mathbf { X } ^ { K } - \mathbf { X } ^ { N } \|$ , where $\mathbf { X } ^ { K }$ and ${ \mathbf { X } } ^ { N }$ are the position vectors to two adjacent nodes in an axial plane. The drill rotation constraint then requires that the component of rotation about the normal to the surface match the in-plane rotation of the surface as measured by the displacement field. The linearized drill strain calculation is similar to the transverse shear linearized strain calculation. Without repeating the calculation, $$ \delta \gamma_ {d} = \delta \pmb {\phi} \cdot (\mathbf {a} \times \mathbf {b}) + \frac {1}{\| \overline {{\mathbf {x}}} _ {, \theta} \|} (\mathbf {a} - \gamma_ {d} \mathbf {b}) \cdot \delta \overline {{\mathbf {x}}} _ {, \theta}. $$ Similarly, the second variation of the drill strain # Elements $$ \begin{array}{l} \mathrm{d} (\delta \gamma_ {d}) = \delta \pmb {\phi} \cdot \frac {1}{2} [ \mathbf {a b} + \mathbf {b a} ] \cdot \mathrm{d} \pmb {\phi} - \delta \overline {{\mathbf {x}}} _ {, \theta} \cdot \frac {1}{\| \overline {{\mathbf {x}}} _ {, \theta} \| ^ {2}} [ \mathbf {a b} + \mathbf {b a} ] \cdot \mathrm{d} \overline {{\mathbf {x}}} _ {, \theta} \\ + \delta \pmb {\phi} \cdot \frac {1}{\| \overline {{\mathbf {x}}} , _ {\theta} \|} [ \hat {\mathbf {a}} - (\mathbf {a} \times \mathbf {b}) \mathbf {b} ] \cdot \mathrm{d} \overline {{\mathbf {x}}} _ {, \theta} + \delta \overline {{\mathbf {x}}} _ {, \theta} \cdot \frac {1}{\| \overline {{\mathbf {x}}} , _ {\theta} \|} [ - \hat {\mathbf {a}} - \mathbf {b} (\mathbf {a} \times \mathbf {b}) ] \cdot \mathrm{d} \pmb {\phi}, \\ \end{array} $$ where terms proportional to $\gamma _ { d }$ have been neglected and the coupling between $\delta \phi$ and $\Delta \phi$ has been symmetrized. # Zero radius: collapsed edge For the case when the reference radius of any point on the shell surface goes to zero, all of the offset nodes collapse to the same point and the edge constraints along that circumferential edge become redundant. It is, therefore, necessary to treat the zero radius case separately. For the zero radius case all redundant degrees of freedom are constrained to follow the average motion of the nodes at $\theta = 0 ^ { \circ }$ and $\theta = 1 8 0 ^ { \circ }$ . The edge constraints are broken into two parts: First, a circumferential transverse shear strain is defined that requires that the rotation in the radial direction follow the circumferential rotations at the first and last nodal plane: $$ \gamma_ {3 2} ^ {0} = \phi_ {r} + \frac {1}{2} (\phi_ {\theta} ^ {1} + \phi_ {\theta} ^ {P + 1}) \sin \theta . $$ Introducing the interpolations, the linearized strain is $$ \delta \gamma_ {3 2} ^ {0} = \sum_ {p = 1} ^ {P} \sin p \theta \delta \phi_ {r} ^ {p} + \frac {1}{2} \sin \theta (\delta \phi_ {\theta} + \delta \phi_ {\theta} ^ {P + 1}). $$ Note that $\mathrm { d } ( \delta \gamma _ { 3 2 } ^ { 0 } ) = 0$ . Second, a drill rotation strain is defined that requires that the rotation about the Z-axis be zero: $$ \gamma_ {d} ^ {0} = \phi_ {z}. $$ This leads to a linearized strain given by $$ \delta \gamma_ {d} ^ {0} = \sum_ {p = 1} ^ {P} \sin p \theta \delta \phi_ {z} ^ {p}. $$ Note that $\mathrm { d } ( \delta \gamma _ { d } ^ { 0 } ) = 0$ # Mass matrix At each material point the displacement components in the three directions (radial, axial, circumferential) are dependent only on the corresponding nodal displacement components. Hence, the mass matrix does not involve any coupling between the radial, axial, and circumferential degrees of freedom, and we can write the mass matrix in the form of three separate expressions: # Elements $$ M _ {u _ {r}} ^ {m n p q} = \int_ {z} \rho \mathrm{d} z \int_ {A} N _ {r} ^ {m p} N _ {r} ^ {n q} \mathrm{d} A $$ $$ M _ {u _ {z}} ^ {m n p q} = \int_ {z} \rho \mathrm{d} z \int_ {A} N _ {z} ^ {m p} N _ {z} ^ {n q} \mathrm{d} A $$ $$ M _ {u _ {\theta}} ^ {m n p q} = \int_ {z} \rho \mathrm{d} z \int_ {A} N _ {\theta} ^ {m p} N _ {\theta} ^ {n q} \mathrm{d} A. $$ Similarly, we can write the terms associated with the rotational degrees of freedom as: $$ M _ {\phi_ {r}} ^ {m n p q} = \int_ {z} \rho z ^ {2} \mathrm{d} z \int_ {A} N _ {\theta} ^ {m p} N _ {\theta} ^ {n q} \mathrm{d} A $$ $$ M _ {\phi_ {z}} ^ {m n p q} = \int_ {z} \rho z ^ {2} \mathrm{d} z \int_ {A} N _ {\theta} ^ {m p} N _ {\theta} ^ {n q} \mathrm{d} A $$ $$ M _ {\phi_ {\theta}} ^ {m n p q} = \int_ {z} \rho z ^ {2} \mathrm{d} z \int_ {A} N _ {r} ^ {m p} N _ {r} ^ {n q} \mathrm{d} A. $$ Here the superscripts m and n refer to a particular node in the $_ { r - z }$ plane, and the superscripts $p$ and $q$ refer to a particular position along the circumference. The interpolation functions $N _ { r } ^ { m p } , N _ { z } ^ { m p }$ , and $N _ { \theta } ^ { m p }$ are the product of interpolation functions $H ^ { m } ( S )$ in the $\displaystyle r - z$ plane and interpolation functions in the µ-direction: $$ N _ {r} ^ {m p} = H ^ {m} (S) R ^ {p} (\theta) = N _ {z} ^ {m p} $$ $$ N _ {\theta} ^ {m p} = H ^ {m} (S) \sin p \theta . $$ The area integral used to form the mass matrix can be split then into an integral along the length of the element in the $_ { r - z }$ plane and an integral around the circumference. For the $r { - } r$ component of the mass matrix this yields $$ M _ {u _ {r}} ^ {m n p q} = \int_ {z} \rho \mathrm{d} z \int_ {S} H ^ {m} (S) H ^ {n} (S) 2 \pi R \mathrm{d} S \frac {1}{2 \pi} \int_ {0} ^ {2 \pi} R ^ {p} (\theta) R ^ {q} (\theta) \mathrm{d} \theta , $$ and for the $\phi _ { r }$ component of the mass matrix this yields: $$ M _ {\phi_ {r}} ^ {m n p q} = \int_ {z} \rho z ^ {2} \mathrm{d} z \int_ {S} H ^ {m} (S) H ^ {n} (S) 2 \pi R \mathrm{d} S \frac {1}{2 \pi} \int_ {0} ^ {2 \pi} R ^ {p} (\theta) R ^ {q} (\theta) \mathrm{d} \theta . $$ The matrix can be written in a convenient form by defining the primitive mass matrix as $$ M _ {p r i m} ^ {m n} = \int_ {z} \rho \mathrm{d} z \int_ {S} H ^ {m} (S) H ^ {n} (S) 2 \pi R \mathrm{d} S, $$ or for the rotational components as