김경종
33 KiB
symmetry exists in the r-z plane at \theta = 0 . Hence in-plane bending of the model is permitted, while deformations such as torsion about the axis of symmetry are precluded. The symmetries of the undeformed configuration and of the deformation are exploited through the assumption of particular displacement and rotation interpolations around the circumference of the shell. Specifically, Fourier series expansions are used in the µ or circumferential direction that preserve the plane of symmetry.
Geometric description
Let ( S , \theta ) be coordinate functions parametrizing the reference surface of the shell and let S _ { 3 } \in [ - h / 2 , h / 2 ] be the coordinate function in the thickness direction, where h is the shell's initial thickness. (For a detailed account of the geometric description of the finite-strain shell formulation, see ``Finite-strain shell element formulation,'' Section 3.6.5.) Then points in the reference or undeformed configuration are identified by the normal coordinates mapping
\mathbf {X} (S, \theta , S _ {3}) = \overline {{\mathbf {X}}} (S, \theta) + S _ {3} \mathbf {T} _ {3} (S, \theta),
where X is the three-dimensional position of a material point, \overline { { \mathbf { X } } } is the shell reference surface mapping, and \mathbf { T } _ { 3 } is the unit normal to the shell reference surface. The fact that \mathbf { T } _ { 3 } is a unit vector assumes that the reference configuration is (locally) of constant thickness. Owing to the axisymmetric reference configuration, \overline { { \mathbf { X } } } can be given relative to a global Cartesian coordinate system as
\overline {{{\mathbf {X}}}} (S, \theta) = \left\{ \begin{array}{l} X ^ {1} (S, \theta) \\ X ^ {2} (S, \theta) \\ X ^ {3} (S, \theta) \end{array} \right\} = \left\{ \begin{array}{c} r (S) \cos \theta \\ r (S) \sin \theta \\ z (S) \end{array} \right\},
where r ( S ) is the radius, z ( S ) is the axial position, and ( r , z , \theta ) are the cylindrical coordinates. (Note that the usual convention for cylindrical coordinates ( r , \theta , z ) has been changed, which is consistent with the axisymmetric shell elements and the axisymmetric elements allowing nonlinear bending.) By definition the normal field to the shell reference surface is \mathbf { T } _ { 3 } = \overline { { \mathbf { X } } } _ { , S } \times \overline { { \mathbf { X } } } _ { , \theta } / \| \overline { { \mathbf { X } } } _ { , S } \times \overline { { \mathbf { X } } } _ { , \theta } \| , which by direct computation yields
\mathbf {T} _ {3} = \frac {1}{[ (r ^ {\prime}) ^ {2} + (z ^ {\prime}) ^ {2} ] ^ {\frac {1}{2}}} \left\{ \begin{array}{c} - z ^ {\prime} \cos \theta \\ - z ^ {\prime} \sin \theta \\ r ^ {\prime} \end{array} \right\} \quad \mathrm{with} \quad \| \overline {{\mathbf {X}}} _ {, S} \times \overline {{\mathbf {X}}} _ {, \theta} \| = r [ (r ^ {\prime}) ^ {2} + (z ^ {\prime}) ^ {2} ] ^ {\frac {1}{2}},
where r ^ { \prime } = \partial r / \partial S and z ^ { \prime } = \partial z / \partial S . Relative to the cylindrical coordinate system,
\mathbf {T} _ {3} = [ (r ^ {\prime}) ^ {2} + (z ^ {\prime}) ^ {2} ] ^ {- \frac {1}{2}} (- z ^ {\prime} \mathbf {e} _ {r} + r ^ {\prime} \mathbf {e} _ {z}).
The basic kinematic assumption is that for any deformed configuration, the position of a point in the body can be identified by
Equation 3.6.7-1
\mathbf {x} (S, \theta , S _ {3}) = \overline {{\mathbf {x}}} (S, \theta) + S _ {3} \overline {{f}} _ {3 3} (S, \theta) \mathbf {t} _ {3} (S, \theta),
where x is the deformed position of the material point, x is the deformed shell reference surface
Elements
mapping, \mathbf { t } _ { 3 } is the deformed unit director field, and \overline { { f } } _ { 3 3 } is the thickness change parameter. \mathrm { o f } critical importance for any shell formulation is the treatment of the rotation field; that is, the treatment of the director field \mathbf { t } _ { 3 } . The geometric description and the incremental update procedure for the director field are given in detail below.
Under the kinematic assumption above, the deformed configuration of the shell is completely determined by the reference surface mapping \overline { { \mathbf { x } } } , the deformed director field \mathbf { t } _ { 3 } and the thickness parameter \overline { { f } } _ { 3 3 } .
We define the following displacement quantities. Since \overline { { \bf x } } is an element of a (linear) vector space, we can define the reference surface displacement vector u by the difference between the deformed reference surface and the undeformed reference surface; i.e.,
\mathbf {u} (S, \theta) = \overline {{\mathbf {x}}} (S, \theta) - \overline {{\mathbf {X}}} (S, \theta).
The director field, however, is a unit vector field that is not a member of a linear vector space. The orientation of the director field is defined in terms of a rotation vector \phi as
\mathbf {t} _ {3} (S, \theta) = \exp [ \hat {\pmb {\phi}} (S, \theta) ] \cdot \mathbf {T} _ {3} (S, \theta).
Here \hat { \pmb { \phi } } is the skew-symmetric matrix with axial vector \phi , defined by the properties
\hat {\boldsymbol {\phi}} \cdot \boldsymbol {\phi} = \mathbf {0} \quad \text {and} \quad \hat {\boldsymbol {\phi}} \cdot \mathbf {v} = \boldsymbol {\phi} \times \mathbf {v} \text {for all vectors} \mathbf {v}
and \exp [ \hat { \phi } ] is an orthogonal transformation given by the closed-form expression
\exp [ \hat {\pmb {\phi}} ] = \cos \| \pmb {\phi} \| \mathbf {I} + \frac {\sin \| \pmb {\phi} \|}{\| \pmb {\phi} \|} \hat {\pmb {\phi}} + \frac {(1 - \cos \| \pmb {\phi} \|)}{\| \pmb {\phi} \| ^ {2}} \pmb {\phi} \pmb {\phi}.
Alternatively, quaternion algebra can be used to specify the orientation of the deformed director field \mathbf { t } _ { 3 } . In this case the orthogonal matrix \exp [ \hat { \phi } ] is replaced by the quaternion parameter { \pmb q } = ( q _ { 0 } , { \pmb q } ) , where
q _ {0} = \cos \left\| \phi / 2 \right\| \quad \mathrm{and} \quad \mathbf {q} = \frac {\sin \left\| \phi / 2 \right\|}{\left\| \phi \right\|} \phi .
The orientation of the unit director field then follows as
\mathbf {t} _ {3} = (2 q _ {0} ^ {2} - 1) \mathbf {T} _ {3} + 2 q _ {0} \mathbf {q} \times \mathbf {T} _ {3} + 2 (\mathbf {q} \cdot \mathbf {T} _ {3}) \mathbf {q}.
Similarly, the orthogonal transformation \exp [ \hat { \phi } ] can be extracted from the quaternion parameter as
\exp [ \hat {\pmb {\phi}} ] = (2 q _ {0} ^ {2} - 1) \mathbf {I} + 2 q _ {0} \hat {\mathbf {q}} + 2 \mathbf {q} \mathbf {q}.
Interpolations
Displacement and rotation components are given relative to the cylindrical coordinate system ( r , z , \theta ) with orthonormal basis vectors \left\{ \mathbf { e } _ { r } , \mathbf { e } _ { z } , \mathbf { e } _ { \theta } \right\} that are fixed in the reference or undeformed configuration. A general interpolation scheme for \mathbf { u } = u _ { r } \mathbf { e } _ { r } + u _ { z } \mathbf { e } _ { z } + u _ { \theta } \mathbf { e } _ { \theta } and \pmb { \phi } = \phi _ { r } \mathbf { e } _ { r } + \phi _ { z } \mathbf { e } _ { z } + \phi _ { \theta } \mathbf { e } _ { \theta } using a Fourier expansion in the \theta variable is
\left. \begin{array}{l} u _ {i} (S, \theta) = \sum_ {m = 1} ^ {M} H ^ {m} (S) \left(u _ {i} ^ {m 0} + \sum_ {p = 1} ^ {P} (\cos (p \theta) u _ {i c} ^ {m p} + \sin (p \theta) u _ {i s} ^ {m p})\right) \\ \phi_ {i} (S, \theta) = \sum_ {m = 1} ^ {M} H ^ {m} (S) \left(\phi_ {i} ^ {m 0} + \sum_ {p = 1} ^ {P} (\cos (p \theta) \phi_ {i c} ^ {m p} + \sin (p \theta) \phi_ {i s} ^ {m p})\right) \end{array} \right\} \quad (i = r, z, \theta).
Here H ^ { m } ( S ) are the polynomial interpolation functions along the generator lines of the axisymmetric reference configuration; u _ { i } ^ { m 0 } , u _ { i c } ^ { m p } , u _ { i s } ^ { m p } , \phi _ { i } ^ { m 0 } , \phi _ { i c } ^ { m p } , \phi _ { i s } ^ { m p } u ic , uis , Ám0i , are the solution amplitude values (Fourier coefficients); M is the number of terms used in the interpolation along the generator lines; and P is the number of Fourier interpolation terms used around the circumference of the reference shell. Note that an axisymmetric deformation is obtained for the choice P = 0 .
The symmetry requirement in the r-z plane at \theta = 2 n \pi , n = 0 , 1 , . . . , eliminates many of the above Fourier coefficients. For the displacement vector the only admissible terms are
\left\{ \begin{array}{l} u _ {r} (S, \theta) \\ u _ {z} (S, \theta) \\ u _ {\theta} (S, \theta) \end{array} \right\} = \sum_ {m = 1} ^ {M} H ^ {m} (S) \left(\left\{ \begin{array}{c} u _ {r} ^ {m 0} \\ u _ {z} ^ {m 0} \\ 0 \end{array} \right\} + \sum_ {p = 1} ^ {P} \cos (p \theta) \left\{ \begin{array}{c} u _ {r c} ^ {m p} \\ u _ {z c} ^ {m p} \\ 0 \end{array} \right\} + \sum_ {p = 1} ^ {P} \sin (p \theta) \left\{ \begin{array}{c} 0 \\ 0 \\ u _ {\theta s} ^ {m p} \end{array} \right\}\right).
For the rotation components, symmetry requirements switch the role of the r and z components with the \theta components:
\left\{ \begin{array}{l} \phi_ {r} (S, \theta) \\ \phi_ {z} (S, \theta) \\ \phi_ {\theta} (S, \theta) \end{array} \right\} = \sum_ {m = 1} ^ {M} H ^ {m} (S) \left(\left\{ \begin{array}{l} 0 \\ 0 \\ \phi_ {\theta} ^ {m 0} \end{array} \right\} + \sum_ {p = 1} ^ {P} \cos (p \theta) \left\{ \begin{array}{l} 0 \\ 0 \\ \phi_ {\theta c} ^ {m p} \end{array} \right\} + \sum_ {p = 1} ^ {P} \sin (p \theta) \left\{ \begin{array}{l} \phi_ {r s} ^ {m p} \\ \phi_ {z s} ^ {m p} \\ 0 \end{array} \right\}\right).
For practical reasons the values of u _ { r } , u _ { z } . , and \phi _ { \theta } are often required at specific locations around the circumference of the shell. Therefore, displacement and rotation components u _ { r } ^ { m p } , u _ { z } ^ { m p } , and \phi _ { \theta } ^ { m p } are used instead of the Fourier coefficients u _ { r c } ^ { m p } , u _ { z c } ^ { m p } , and \phi _ { \theta c } ^ { m p } . Furthermore, a negative sign is introduced in the interpolation for \phi _ { \theta } for the following reason: The ABAQUS convention for axisymmetric shell elements is that the axial tangent direction is drawn between nodes in ascending node number (the shell local 1-direction). The normal to the shell is then obtained by a 9 0 ^ { \circ } counter-clockwise rotation of the tangent (the shell local 3-direction). However, a positive rotation of the normal field (about the shell local 2-direction) is counterclockwise. This convention implies a left-handed shell local coordinate system. For the axisymmetric shell bending elements, a right-handed shell local coordinate system is required at the integration points; thus, the direction of positive rotation is reversed there.
Elements
Rearranging the Fourier series expansions and making the substitution u _ { \theta s } ^ { m p } \mapsto u _ { \theta } ^ { m p } 7! um µ , the interpolations for the displacement components are
Equation 3.6.7-2
\left\{ \begin{array}{l} u _ {r} (S, \theta) \\ u _ {z} (S, \theta) \\ 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} u _ {r} ^ {m p} \\ u _ {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).
Similarly, replacing \phi _ { r s } ^ { m p } and \phi _ { z s } ^ { m p } with \phi _ { r } ^ { m p } and \phi _ { z } ^ { m p } respectively, the interpolation for the rotation components becomes
Equation 3.6.7-3
\left\{ \begin{array}{l} \phi_ {r} (S, \theta) \\ \phi_ {z} (S, \theta) \\ \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 \\ - \phi_ {\theta} ^ {m p} \end{array} \right\} + \sum_ {p = 1} ^ {P} \sin (p \theta) \left\{ \begin{array}{c} \phi_ {r} ^ {m p} \\ \phi_ {z} ^ {m p} \\ 0 \end{array} \right\}\right).
In the above interpolations, u _ { r } ^ { m p } , u _ { z } ^ { m p } , and \phi _ { \theta } ^ { m p } are physical displacement and rotation components at \theta ^ { p } = \pi ( p - 1 ) / P and R ^ { p } ( \theta ) are trigonometric interpolation functions with the property that
R ^ {p} (\theta^ {q}) = \left\{ \begin{array}{l} 1 , p = q \\ 0 , p \neq q \end{array} \right. \text { defined by: }
P = 1 \colon \theta^ {1} = 0, \theta^ {2} = \pi ,
R ^ {1} (\theta) = \frac {1}{2} (1 + \cos \theta)
R ^ {2} (\theta) = \frac {1}{2} (1 - \cos \theta)
P = 2 \colon \theta^ {1} = 0, \theta^ {2} = \pi / 2, \theta^ {3} = \pi ,
R ^ {1} (\theta) = \frac {1}{4} (1 + 2 \cos \theta + \cos 2 \theta)
R ^ {2} (\theta) = \frac {1}{2} (1 - \cos 2 \theta)
R ^ {3} (\theta) = \frac {1}{4} (1 - 2 \cos \theta + \cos 2 \theta)
P = 3 \colon \theta^ {1} = 0, \theta^ {2} = \pi / 3, \theta^ {3} = 2 \pi / 3, \theta^ {4} = \pi ,
Elements
R ^ {1} (\theta) = \frac {1}{6} (1 + 2 \cos \theta + 2 \cos 2 \theta + \cos 3 \theta)
R ^ {2} (\theta) = \frac {1}{3} (1 + \cos \theta - \cos 2 \theta - \cos 3 \theta)
R ^ {3} (\theta) = \frac {1}{3} (1 - \cos \theta - \cos 2 \theta + \cos 3 \theta)
R ^ {4} (\theta) = \frac {1}{6} (1 - 2 \cos \theta + 2 \cos 2 \theta - \cos 3 \theta)
P = 4 \text {:} \theta^ {1} = 0, \theta^ {2} = \pi / 4, \theta^ {3} = \pi / 2, \theta^ {4} = 3 \pi / 4, \theta^ {5} = \pi ,
R ^ {1} (\theta) = \frac {1}{8} (1 + 2 \cos \theta + 2 \cos 2 \theta + 2 \cos 3 \theta + \cos 4 \theta)
R ^ {2} (\theta) = \frac {1}{4} (1 + \sqrt {2} \cos \theta - \sqrt {2} \cos 3 \theta - \cos 4 \theta)
R ^ {3} (\theta) = \frac {1}{4} (1 - 2 \cos 2 \theta + \cos 4 \theta)
R ^ {4} (\theta) = \frac {1}{4} (1 - \sqrt {2} \cos \theta + \sqrt {2} \cos 3 \theta - \cos 4 \theta)
R ^ {5} (\theta) = \frac {1}{8} (1 - 2 \cos \theta + 2 \cos 2 \theta - 2 \cos 3 \theta + \cos 4 \theta)
As with the continuum axisymmetric bending element, P = 4 is the highest-order interpolation offered with respect to µ. The element becomes significantly more expensive as higher-order interpolations are used, and it is assumed that the general-purpose finite-strain shell is less expensive than using this element with P > 4 .
Virtual work
The virtual work expression from the three-dimensional theory is
\delta \Pi = \int_ {V} \sigma^ {i j} \delta E _ {i j} d V,
where V is the current volume of the deformed body, \sigma ^ { i j } are the curvilinear components of the Cauchy stress tensor, E _ { i j } are the components of the Lagrange strain tensor, and \delta E _ { i j } are the variational or linearized strain measure components. By definition the Lagrange strain tensor components are given by
E _ {i j} \stackrel {\mathrm{def}} {=} \frac {1}{2} \left(\frac {\partial \mathbf {x}}{\partial \xi^ {i}} \cdot \frac {\partial \mathbf {x}}{\partial \xi^ {j}} - \frac {\partial \mathbf {X}}{\partial \xi^ {i}} \cdot \frac {\partial \mathbf {X}}{\partial \xi^ {j}}\right).
Note that in the statement of virtual work, no choice has thus far been made regarding the curvilinear coordinate functions \xi ^ { i } ( i = 1 , 2 , 3 ) . Furthermore, the current volume measure dV is given by the parametric relationship
Elements
d V = j d \xi^ {1} d \xi^ {2} d \xi^ {3}, \quad \mathrm{where} \quad j \stackrel {\mathrm{def}} {=} \det [ \nabla \mathbf {x} ] = \frac {\partial \mathbf {x}}{\partial \xi^ {1}} \times \frac {\partial \mathbf {x}}{\partial \xi^ {2}} \cdot \frac {\partial \mathbf {x}}{\partial \xi^ {3}}.
We now introduce the kinematic assumption Equation 3.6.7-1 into the definition of E _ { i j } to find
\begin{array}{l} E _ {\alpha \beta} = \frac {1}{2} (\overline {{\mathbf {x}}}, _ {\alpha} \cdot \overline {{\mathbf {x}}}, _ {\beta} - \overline {{\mathbf {X}}}, _ {\alpha} \cdot \overline {{\mathbf {X}}}, _ {\beta}) + \frac {1}{2} S _ {3} [ \overline {{\mathbf {x}}}, _ {\alpha} \cdot (\overline {{f}} _ {3 3} \mathbf {t} _ {3}) _ {, \beta} + \overline {{\mathbf {x}}}, _ {\beta} \cdot (\overline {{f}} _ {3 3} \mathbf {t} _ {3}) _ {, \alpha} \\ - \overline {{\mathbf {X}}} _ {, \alpha} \cdot \mathbf {T} _ {3, \beta} - \overline {{\mathbf {X}}} _ {, \beta} \cdot \mathbf {T} _ {3, \alpha} ] + \frac {1}{2} (S _ {3}) ^ {2} [ (\overline {{f}} _ {3 3} \mathbf {t} _ {3}) _ {, \alpha} \cdot (\overline {{f}} _ {3 3} \mathbf {t} _ {3}) _ {, \beta} - \mathbf {T} _ {3, \alpha} \cdot \mathbf {T} _ {3, \beta} ], \\ \end{array}
E _ {3 \alpha} = \frac {1}{2} (\overline {{f}} _ {3 3} \overline {{\mathbf {x}}} _ {, \alpha} \cdot \mathbf {t} _ {3} + S _ {3} \overline {{f}} _ {3 3} \overline {{f}} _ {3 3, \alpha}),
E _ {3 3} = \frac {1}{2} [ (\overline {{f}} _ {3 3}) ^ {2} - 1 ],
where differentiation is now with respect to the parametric coordinates, so that \alpha , \beta \in \{ S , \theta \} and 3 = S _ { 3 } . Define the following shell strain measure components and kinematic relationships:
\begin{array}{l} \epsilon_ {\alpha \beta} \stackrel {{\text {def}}} {{=}} \frac {1}{2} (\overline {{\mathbf {x}}}, _ {\alpha} \cdot \overline {{\mathbf {x}}}, _ {\beta} - \overline {{\mathbf {X}}}, _ {\alpha} \cdot \overline {{\mathbf {X}}}, _ {\beta}) \qquad \qquad \text {membrane strain}, \\ \kappa_ {\alpha \beta} \stackrel {\mathrm{def}} {=} \frac {1}{2} \big (\overline {{{\mathbf {x}}}} _ {, \alpha} \cdot \mathbf {t} _ {3, \beta} + \overline {{{\mathbf {x}}}} _ {, \beta} \cdot \mathbf {t} _ {3, \alpha} \big) - b _ {\alpha \beta} ^ {0} \qquad \mathrm{bendingstrain}, \\ \gamma_ {\alpha} \stackrel {{\text {def}}} {{=}} \overline {{\mathbf {x}}} _ {, \alpha} \cdot \mathbf {t} _ {3} \quad \text {transverse shear}, \\ \mu \stackrel {{\text { def }}} {{=}} \ln (\overline {{f}} _ {3 3}) \quad \text { logarithmic stretch. } \\ \end{array}
In the above, b _ { \alpha \beta } ^ { 0 } \ { \stackrel { \mathrm { d e f } } { = } } \ { \overline { { \mathbf { X } } } } _ { , \alpha } \cdot \mathbf { T } _ { 3 , \beta } = b _ { \beta \alpha } ^ { 0 } are the components of the second fundamental form of the undeformed reference surface.
Substituting the above definitions into the virtual work expression, we find (after some manipulation) that the volume integral reduces to the following integral over the deformed reference surface
\delta \Pi = \int_ {A} \left[ N ^ {\alpha \beta} \delta \epsilon_ {\alpha \beta} + M ^ {\alpha \beta} \delta \kappa_ {\alpha \beta} + Q ^ {\alpha} \delta \gamma_ {\alpha} + L ^ {3} \delta \mu + M ^ {3 \alpha} (\delta \mu) _ {, \alpha} \right] \overline {{j}} d S d \theta ,
where \delta \mu = \delta \overline { { f } } _ { 3 3 } / \overline { { f } } _ { 3 3 } and the current reference surface Jacobian determinant is \overline { { j } } = \| \overline { { \mathbf { x } } } _ { , 1 } \times \overline { { \mathbf { x } } } _ { , 2 } \| . In the above virtual work expression, the ( S _ { 3 } ) ^ { 2 } term in E _ { \alpha \beta } has been neglected. This term is O ( h ^ { 2 } / R ^ { 2 } ) --where h is the thickness and R is some characteristic radius of curvature--and is negligible in light of the kinematic assumption Equation 3.6.7-1. The shell stress resultant components are defined by the following integrals through the thickness of the shell:
Elements
N ^ {\alpha \beta} = \frac {1}{\dot {j}} \int_ {h} \sigma^ {\alpha \beta} j d S _ {3},
M ^ {\alpha \beta} = \overline {{f}} _ {3 3} \frac {1}{\overline {{j}}} \int_ {h} S _ {3} \sigma^ {\alpha \beta} j d S _ {3},
Q ^ {\alpha} = \overline {{f}} _ {3 3} \frac {1}{\overline {{j}}} \int_ {h} \left(\sigma^ {3 \alpha} + S _ {3} \sigma^ {\alpha \beta} \mu_ {, \beta}\right) j d S _ {3},
{\cal L} ^ {3} = \overline {{f}} _ {3 3} \frac {1}{\overline {{j}}} \int_ {h} \left[ \sigma^ {3 3} \overline {{f}} _ {3 3} + \sigma^ {3 \alpha} \gamma_ {\alpha} + 2 S _ {3} \sigma^ {3 \alpha} \overline {{f}} _ {3 3, \alpha} + S _ {3} \sigma^ {\alpha \beta} (\kappa_ {\alpha \beta} + \gamma_ {\alpha} \mu_ {, \beta}) \right] j d S _ {3},
M ^ {3 \alpha} = \overline {{f}} _ {3 3} \frac {1}{\overline {{j}}} \int_ {h} S _ {3} \left(\sigma^ {3 \alpha} \overline {{f}} _ {3 3} + \sigma^ {\alpha \beta} \gamma_ {\beta}\right) j d S _ {3}.
For thin shells the Kirchhoff-Love approximation, which states that the deformed director field \mathbf { t } _ { 3 } is (approximately) the normal field to the deformed reference surface, is introduced along with the plane stress assumption \sigma ^ { 3 3 } = 0 . Consistent with these approximations, we neglect all O ( h / R ) terms and terms proportional to the gradient of the thickness parameter. Accordingly, we set
\gamma_ {\alpha} = 0, \quad \frac {j}{\overline {{j}}} = \overline {{f}} _ {3 3} \quad \mathrm{and} \quad \overline {{f}} _ {3 3, \alpha} = \mu_ {, \alpha} = 0.
We can now summarize the virtual work expression for thin (Kirchhoff-Love) shells:
Equation 3.6.7-4
\delta \Pi = \int_ {A} \left[ N ^ {\alpha \beta} \delta \epsilon_ {\alpha \beta} + M ^ {\alpha \beta} \delta \kappa_ {\alpha \beta} + Q ^ {\alpha} \delta \gamma_ {\alpha} \right] \overline {{j}} d S d \theta ,
where the shell resultant components are defined in terms of the Cauchy stress tensor components by the integrals
N ^ {\alpha \beta} = \overline {{f}} _ {3 3} \int_ {h} \sigma^ {\alpha \beta} d S _ {3} \quad \mathrm{and} \quad M ^ {\alpha \beta} = (\overline {{f}} _ {3 3}) ^ {2} \int_ {h} S _ {3} \sigma^ {\alpha \beta} d S _ {3}.
In the expression for ±¦, \begin{array} { r } { Q ^ { \alpha } = ( \overline { { f } } _ { 3 3 } ) ^ { 2 } \int _ { h } \sigma ^ { 3 \alpha } d S _ { 3 } } \end{array} is interpreted as a constraint stress that enforces that the director field remain normal to the reference surface. Two other contributions to the virtual work expression M ^ { \alpha \beta } \kappa _ { \alpha \beta } \delta \mu and M ^ { 3 \alpha } ( \delta \mu ) _ { , \alpha } , where \begin{array} { r } { M ^ { 3 \alpha } = ( \overline { { f } } _ { 3 3 } ) ^ { 2 } \int _ { h } S _ { 3 } \sigma ^ { 3 \alpha } d S _ { 3 } } \end{array} are O ( h / R ) and, thus, neglected.
Orthonormal surface coordinate system and coordinate transformation
It is desirable to define stress resultant quantities relative to an orthonormal basis in the deformed configuration. To do this, we define a normal coordinate system \{ \mathbf { n } _ { 1 } , \mathbf { n } _ { 2 } , \mathbf { n } _ { 3 } \} , where \mathbf { n } _ { 1 } and \mathbf { n } _ { 2 } are tangent to the deformed reference surface and { \bf n } _ { 3 } is the unit normal field.
Define the following notation. Let \overline { { f } } _ { \alpha } ^ { a } ( a = 1 , 2 ) be the components of \overline { { \mathbf { x } } } _ { , \alpha } relative to the basis \{ \mathbf { n } _ { 1 } , \mathbf { n } _ { 2 } \} ; that is
\frac {\partial \overline {{\mathbf {x}}}}{\partial \xi^ {\alpha}} = \overline {{f}} _ {\alpha} ^ {a} \mathbf {n} _ {a}.
Furthermore, let the inverse of the matrix of components \textstyle { \overline { { f } } } _ { \alpha } ^ { a } ] be given by [ \overline { { h } } _ { a } ^ { \alpha } ] , such that
\overline {{f}} _ {\alpha} ^ {a} \overline {{h}} _ {b} ^ {\alpha} = \delta_ {b} ^ {a} \quad \mathrm{and} \quad \overline {{h}} _ {a} ^ {\alpha} \overline {{f}} _ {\beta} ^ {a} = \delta_ {\beta} ^ {\alpha}.
Note that the basis vectors \mathbf { n } _ { 1 } and \mathbf { n } _ { 2 } induce distance measuring coordinates s ^ { 1 } and s ^ { 2 } such that
\overline {{f}} _ {\alpha} ^ {a} = \frac {\partial s ^ {a}}{\partial \xi^ {\alpha}} \quad \mathrm{and} \quad \overline {{h}} _ {a} ^ {\alpha} = \frac {\partial \xi^ {\alpha}}{\partial s ^ {a}}.
It follows from Equation 3.6.7-5 that the orthonormal tangent vectors are given by
\mathbf {n} _ {a} = \overline {{h}} _ {a} ^ {\alpha} \overline {{\mathbf {x}}} _ {, \alpha} = \frac {\partial \overline {{\mathbf {x}}}}{\partial s ^ {a}}.
For the material calculations, it is important to express both the strain and stress quantities relative to the local orthonormal frame \{ \mathbf { n } _ { 1 } , \mathbf { n } _ { 2 } \} . Accordingly, let N ^ { a b } = N ^ { b a } and M ^ { a b } = M ^ { b a } be the membrane and bending stress resultant components relative to this local orthonormal basis. Thus, we can write
N ^ {\alpha \beta} = N ^ {a b} \overline {{h}} _ {a} ^ {\alpha} \overline {{h}} _ {b} ^ {\beta} \quad \mathrm{and} \quad M ^ {\alpha \beta} = M ^ {a b} \overline {{h}} _ {a} ^ {\alpha} \overline {{h}} _ {b} ^ {\beta},
where we recall that \overline { { h } } _ { a } ^ { \alpha } \ \overset { \mathrm { d e f } } { = } \ \partial \xi ^ { \alpha } / \partial s ^ { a } . Then the stress resultant contributions to the virtual work expression can be transformed as follows. First, the membrane contribution:
\begin{array}{l} N ^ {\alpha \beta} \delta \epsilon_ {\alpha \beta} = N ^ {\alpha \beta} \overline {{\mathbf {x}}} _ {, \alpha} \cdot (\delta \overline {{\mathbf {x}}}) _ {, \beta} \\ = N ^ {a b} \overline {{h}} _ {a} ^ {\alpha} \overline {{h}} _ {b} ^ {\beta} \overline {{\mathbf {x}}} _ {, \alpha} \cdot (\delta \overline {{\mathbf {x}}}) _ {, \beta} \\ = N ^ {a b} \frac {\partial \overline {{\mathbf {x}}}}{\partial s ^ {a}} \cdot \frac {\partial (\delta \overline {{\mathbf {x}}})}{\partial s ^ {b}}. \\ \end{array}
Recall, however, that by the definition of the coordinates s ^ { a } , \partial \overline { { \mathbf { x } } } / \partial s ^ { a } = \mathbf { n } _ { a } . Thus,
N ^ {\alpha \beta} \delta \epsilon_ {\alpha \beta} = N ^ {a b} \mathrm{sym} \big [ \mathbf {n} _ {a} \cdot \frac {\partial (\delta \overline {{\mathbf {x}}})}{\partial s ^ {b}} \big ].
Similarly, the bending contribution is:
\begin{array}{l} M ^ {\alpha \beta} \delta \kappa_ {\alpha \beta} = M ^ {\alpha \beta} \left[ \overline {{\mathbf {x}}} _ {, \alpha} \cdot (\delta \mathbf {t} _ {3}) _ {, \beta} + \mathbf {t} _ {3} \cdot (\delta \overline {{\mathbf {x}}}) _ {, \alpha} \right] \\ = M ^ {a b} \big [ \mathbf {n} _ {a} \cdot \frac {\partial (\delta \pmb {\phi} \times \mathbf {t} _ {3})}{\partial s ^ {b}} + \frac {\partial \mathbf {t} _ {3}}{\partial s ^ {b}} \cdot \frac {\partial (\delta \overline {{\mathbf {x}}})}{\partial s ^ {a}} \big ]. \\ \end{array}
Let e _ { a } ^ { c } be the two-dimensional alternator, such that e _ { 1 } ^ { 1 } = e _ { 2 } ^ { 2 } = 0 and e _ { 1 } ^ { 2 } = - e _ { 2 } ^ { 1 } = 1 . Then
Elements
\mathbf {t} _ {3} \times \mathbf {n} _ {a} = e _ {a} ^ {c} \mathbf {n} _ {c} \mathrm{and}
M ^ {\alpha \beta} \delta \kappa_ {\alpha \beta} = M ^ {a b} \left[ e _ {a} ^ {c} \mathbf {n} _ {c} \cdot \frac {\partial (\delta \pmb {\phi})}{\partial s ^ {b}} + \frac {\partial \mathbf {t} _ {3}}{\partial s ^ {b}} \cdot \left(\frac {\partial (\delta \overline {{\mathbf {x}}})}{\partial s ^ {a}} - \delta \pmb {\phi} \times \mathbf {n} _ {a}\right) \right].
Since the second term in the brackets is proportional to the bending curvature, we neglect this term relative to the first, yielding
M ^ {\alpha \beta} \delta \kappa_ {\alpha \beta} = M ^ {a b} \mathrm{sym} \big [ e _ {a} ^ {c} \mathbf {n} _ {c} \cdot \frac {\partial (\delta \pmb {\phi})}{\partial s ^ {b}} \big ].
Strain displacement operators
We now write the virtual work expression Equation 3.6.7-4 in matrix operator notation. Define the following stress resultant component vectors:
\mathbf {N} \stackrel {\mathrm{def}} {=} \left\{ \begin{array}{l} N ^ {1 1} \\ N ^ {2 2} \\ N ^ {1 2} \end{array} \right\} \quad \mathrm{and} \quad \mathbf {M} \stackrel {\mathrm{def}} {=} \left\{ \begin{array}{l} M ^ {1 1} \\ M ^ {2 2} \\ M ^ {1 2} \end{array} \right\}.
Define the matrix strain displacement operators as follows:
\mathbf {B} _ {N} = \left[ \begin{array}{c c} \mathbf {n} _ {1} ^ {T} \frac {\partial}{\partial s ^ {1}} & \mathbf {0} _ {1 \times 3} \\ \mathbf {n} _ {2} ^ {T} \frac {\partial}{\partial s ^ {2}} & \mathbf {0} _ {1 \times 3} \\ \mathbf {n} _ {1} ^ {T} \frac {\partial}{\partial s ^ {2}} + \mathbf {n} _ {2} ^ {T} \frac {\partial}{\partial s ^ {1}} & \mathbf {0} _ {1 \times 3} \end{array} \right] \quad \mathrm{and} \quad \mathbf {B} _ {M} = \left[ \begin{array}{c c} \mathbf {0} _ {1 \times 3} & \mathbf {n} _ {2} ^ {T} \frac {\partial}{\partial s ^ {1}} \\ \mathbf {0} _ {1 \times 3} & - \mathbf {n} _ {1} ^ {T} \frac {\partial}{\partial s ^ {2}} \\ \mathbf {0} _ {1 \times 3} & \mathbf {n} _ {2} ^ {T} \frac {\partial}{\partial s ^ {2}} - \mathbf {n} _ {1} ^ {T} \frac {\partial}{\partial s ^ {1}} \end{array} \right],
where { \bf 0 } _ { 1 \times 3 } is the column of zeros \left\{ \begin{array} { l l l } { 0 } & { 0 } & { 0 } \end{array} \right\} ^ { T } . Then the virtual work expression Equation 3.6.7-4 is equivalently stated
\delta \Pi = \int_ {A} \left[ \mathbf {N} \cdot \mathbf {B} _ {N} \left\{ \begin{array}{l} \delta \overline {{\mathbf {x}}} \\ \delta \pmb {\phi} \end{array} \right\} + \mathbf {M} \cdot \mathbf {B} _ {M} \left\{ \begin{array}{l} \delta \overline {{\mathbf {x}}} \\ \delta \pmb {\phi} \end{array} \right\} \right] \overline {{j}} d S d \theta .
Corotational coordinate system
Thus far \mathbf { n } _ { 1 } and \mathbf { n } _ { 2 } are any two orthonormal vectors in the tangent plane to the reference surface. We can uniquely choose these two vectors by requiring the matrix components of the incremental reference surface deformation gradient [ \Delta \overline { { f _ { b } } } ] , where
\Delta \overline {{f}} _ {b} ^ {a} \stackrel {\mathrm{def}} {=} \frac {\partial s ^ {(t + \Delta t) a}}{\partial s ^ {(t) b}} = \overline {{f}} _ {\beta} ^ {(t + \Delta t) a} \overline {{h}} _ {b} ^ {(t) \beta},
to be symmetric; i.e., \Delta \overline { { f } } _ { 2 } ^ { 1 } = \Delta \overline { { f } } _ { 1 } ^ { 2 } . Note that by definition \overline { { f } } _ { \alpha } ^ { ( t + \Delta t ) a } = \Delta \overline { { f } } _ { b } ^ { a } \overline { { f } } _ { \alpha } ^ { ( t ) b } = ¢f ab f (t)b® . This symmetry condition defines a corotational orthonormal basis in the deformed configuration. This orthonormal basis is calculated as follows.
Elements
Obtain the ABAQUS convention pair of orthonormal surface vectors \hat { \mathbf { t } } _ { a } . We then rotate these vectors in the tangent plane to the reference surface about the normal vector \mathbf { n } _ { 3 } ^ { t + \Delta t } by the angle \Delta \psi , where \Delta \psi is determined by the symmetry condition \Delta \overline { { f } } _ { 2 } ^ { 1 } = \Delta \overline { { f } } _ { 1 } ^ { 2 } . Thus, we define
Equation 3.6.7-6
\begin{array}{l} \mathbf {n} _ {1} ^ {t + \Delta t} = \cos (\Delta \psi) \hat {\mathbf {t}} _ {1} + \sin (\Delta \psi) \hat {\mathbf {t}} _ {2}, \\ \mathbf {n} _ {2} ^ {t + \Delta t} = - \sin (\Delta \psi) \hat {\mathbf {t}} _ {1} + \cos (\Delta \psi) \hat {\mathbf {t}} _ {2}. \\ \end{array}
It follows by definition that
\Delta \overline {{f}} _ {b} ^ {a} = \mathbf {n} _ {a} ^ {(t + \Delta t)} \cdot \frac {\partial \overline {{\mathbf {x}}} ^ {(t + \Delta t)}}{\partial \mathcal {S} ^ {(t) b}}.
Therefore, define the quantity
\Delta \hat {f} _ {a b} \stackrel {\mathrm{def}} {=} \hat {\mathbf {t}} _ {a} \cdot \frac {\partial \overline {{\mathbf {x}}} ^ {(t + \Delta t)}}{\partial s ^ {(t) b}}.
The symmetry condition requires that
(\cos (\Delta \psi) \hat {\mathbf {t}} _ {1} + \sin (\Delta \psi) \hat {\mathbf {t}} _ {2}) \cdot \frac {\partial \overline {{{\mathbf {x}}}} ^ {t + \Delta t}}{\partial s ^ {(t) 2}} = (- \sin (\Delta \psi) \hat {\mathbf {t}} _ {1} + \cos (\Delta \psi) \hat {\mathbf {t}} _ {2}) \cdot \frac {\partial \overline {{{\mathbf {x}}}} ^ {t + \Delta t}}{\partial s ^ {(t) 1}}.
From this it can be determined that
\Delta \psi = \tan^ {- 1} \left[ \frac {\Delta \hat {f} _ {2 1} - \Delta \hat {f} _ {1 2}}{\Delta \hat {f} _ {1 1} + \Delta \hat {f} _ {2 2}} \right],
\{ \mathbf { n } _ { 1 } ^ { t + \Delta t } , \mathbf { n } _ { 2 } ^ { t + \Delta t } \} are then determined by Equation 3.6.7-6.
Having determined the updated vectors \mathbf { n } _ { a } ^ { t + \Delta t } , we can calculate the quantities required for coordinate transformation. First the incremental deformation gradient components and the Jacobian matrix components:
\Delta \overline {{f}} _ {b} ^ {a} \equiv \Delta \overline {{f}} _ {a b} = \mathbf {n} _ {a} ^ {(t + \Delta t)} \cdot \frac {\partial \overline {{\mathbf {x}}} ^ {(t + \Delta t)}}{\partial s ^ {(t) b}} \quad \mathrm{and} \quad \overline {{f}} _ {\alpha} ^ {(t + \Delta t) a} = \Delta \overline {{f}} _ {b} ^ {a} \overline {{f}} _ {\alpha} ^ {(t) b}.
Note that due to the symmetry of the incremental deformation gradient components, there is no ambiguity in writing \Delta \overline { { f } } _ { b } ^ { a } = \Delta \overline { { f } } _ { a b } . We can now calculate the inverse components
\overline {{h}} _ {a} ^ {(t + \Delta t) \alpha} = \left[ \overline {{f}} _ {\alpha} ^ {(t + \Delta t) a} \right] ^ {- 1} = \frac {\partial \xi^ {\alpha}}{\partial s ^ {(t + \Delta t) a}}.