김경종
35 KiB
Elements
\bar {B} _ {c c d} = \frac {1}{2} \left[ \begin{array}{l l l l} c _ {\xi} (\mathbf {t} ^ {A} - \mathbf {t} ^ {B}) ^ {T} & c _ {\xi} (\mathbf {t} ^ {B} - \mathbf {t} ^ {C}) ^ {T} & c _ {\xi} (\mathbf {t} ^ {C} - \mathbf {t} ^ {D}) ^ {T} & c _ {\xi} (\mathbf {t} ^ {D} - \mathbf {t} ^ {A}) ^ {T} \\ c _ {\eta} (\mathbf {t} ^ {A} - \mathbf {t} ^ {B}) ^ {T} & c _ {\eta} (\mathbf {t} ^ {B} - \mathbf {t} ^ {C}) ^ {T} & c _ {\eta} (\mathbf {t} ^ {C} - \mathbf {t} ^ {D}) ^ {T} & c _ {\eta} (\mathbf {t} ^ {D} - \mathbf {t} ^ {A}) ^ {T} \end{array} \right]
and
\bar {B} _ {c c r} = \frac {1}{2} \left[ \begin{array}{l l l l} c _ {\xi} (\pmb {\eta} _ {1} ^ {B} - \pmb {\eta} _ {2} ^ {A}) ^ {T} & c _ {\xi} (\pmb {\eta} _ {1} ^ {B} + \pmb {\eta} _ {2} ^ {C}) ^ {T} & c _ {\xi} (\pmb {\eta} _ {2} ^ {C} - \pmb {\eta} _ {1} ^ {D}) ^ {T} & c _ {\xi} (- \pmb {\eta} _ {2} ^ {A} - \pmb {\eta} _ {1} ^ {D}) ^ {T} \\ c _ {\eta} (\pmb {\eta} _ {1} ^ {B} - \pmb {\eta} _ {2} ^ {A}) ^ {T} & c _ {\eta} (\pmb {\eta} _ {1} ^ {B} + \pmb {\eta} _ {2} ^ {C}) ^ {T} & c _ {\eta} (\pmb {\eta} _ {2} ^ {C} - \pmb {\eta} _ {1} ^ {D}) ^ {T} & c _ {\eta} (- \pmb {\eta} _ {2} ^ {A} - \pmb {\eta} _ {1} ^ {D}) ^ {T} \end{array} \right].
The stabilization term has a similar formulation. The variation of the hourglass strain is
\delta \pmb {\gamma} _ {h g} = \bar {B} _ {d} \delta \overline {{\mathbf {x}}} + \bar {B} _ {r} \delta \pmb {\phi},
where
\bar {B} _ {d} = \frac {1}{2} \left[ \begin{array}{c c c c} (\Gamma_ {A} \pmb {\mathsf {t}} ^ {A} - \Gamma_ {B} \pmb {\mathsf {t}} ^ {B}) ^ {T} & (\Gamma_ {B} \pmb {\mathsf {t}} ^ {B} - \Gamma_ {C} \pmb {\mathsf {t}} ^ {C}) ^ {T} & (\Gamma_ {C} \pmb {\mathsf {t}} ^ {C} - \Gamma_ {D} \pmb {\mathsf {t}} ^ {D}) ^ {T} & (\Gamma_ {D} \pmb {\mathsf {t}} ^ {D} - \Gamma_ {A} \pmb {\mathsf {t}} ^ {A}) ^ {T} \\ (\pmb {\mathsf {t}} ^ {A} - \pmb {\mathsf {t}} ^ {B}) ^ {T} & (\pmb {\mathsf {t}} ^ {B} - \pmb {\mathsf {t}} ^ {C}) ^ {T} & (\pmb {\mathsf {t}} ^ {C} - \pmb {\mathsf {t}} ^ {D}) ^ {T} & (\pmb {\mathsf {t}} ^ {D} - \pmb {\mathsf {t}} ^ {A}) ^ {T} \end{array} \right]
and
\bar {B} _ {r} = \frac {1}{2} \left[ \begin{array}{c c c c} (\Gamma_ {B} \pmb {\eta} _ {1} ^ {B} - \Gamma_ {A} \pmb {\eta} _ {2} ^ {A}) ^ {T} & (\Gamma_ {B} \pmb {\eta} _ {1} ^ {B} + \Gamma_ {C} \pmb {\eta} _ {2} ^ {C}) ^ {T} & (\Gamma_ {C} \pmb {\eta} _ {2} ^ {C} - \Gamma_ {D} \pmb {\eta} _ {1} ^ {D}) ^ {T} & (- \Gamma_ {A} \pmb {\eta} _ {2} ^ {A} - \Gamma_ {D} \pmb {\eta} _ {1} ^ {D}) ^ {T} \\ (\pmb {\eta} _ {1} ^ {B} - \pmb {\eta} _ {2} ^ {A}) ^ {T} & (\pmb {\eta} _ {1} ^ {B} + \pmb {\eta} _ {2} ^ {C}) ^ {T} & (\pmb {\eta} _ {2} ^ {C} - \pmb {\eta} _ {1} ^ {D}) ^ {T} & (- \pmb {\eta} _ {2} ^ {A} - \pmb {\eta} _ {1} ^ {D}) ^ {T} \end{array} \right].
The hourglass force components h ^ { 1 } and h ^ { 2 } are given by the constitutive relations
\mathrm{h} ^ {1} = \frac {C _ {s 0} ^ {e f f}}{1 2} \gamma_ {b f} \quad \mathrm{and} \quad h ^ {2} = 0. 0 0 1 \frac {C _ {s 0} ^ {e f f}}{1 2} \gamma_ {c c}.
Comments on stabilization
(1) The butterfly mode \gamma _ { b f } is applied with a "large" or physical hourglass stiffness. For a reference geometry with constant Jacobian, the butterfly stabilization term can be derived from an exact integration of the assumed strain formulation of the transverse shear energy. It is important to apply this constraint with a high stiffness to prevent overly flexible response for quadrilateral elements. The crop circle mode is applied with a "small" or weak stiffness. Although this mode can propagate, it is rarely problematic and is often prevented with boundary conditions.
(2) As the quadrilateral element is degenerated to a triangle, the two hourglass constraints converge into a single constraint: the crop circle constraint. However, as is well-known, for a constant strain triangle the element will lock for certain meshes with three transverse shear constraints per element. Therefore, in the case of a triangular element, the (strong) butterfly mode stabilization is not applied. Only the (weak) crop circle mode stabilization is applied. Thus, in addition to the two homogeneous transverse shear strains, the triangle has a weak constraint to prevent spurious zero energy modes, yet avoids locking in most situations.
Elements
The initial stress contribution from the stabilization terms takes the following form:
A _ {0} (h ^ {1} \mathrm{d} \delta \gamma_ {b f} + h ^ {2} \mathrm{d} \delta \gamma_ {c c}) = A _ {0} \delta \mathbf {u} \cdot \mathbf {K} _ {h} \cdot \mathrm{d} \mathbf {u},
where { \bf K } _ { h } is the (symmetric) transverse shear stabilization contribution to the initial stress. Define the symmetric matrices
\mathbf {H} ^ {A} = \left(h ^ {1} \Gamma_ {A} + h ^ {2}\right) \left[ \mathrm{sym} \{\mathbf {t} ^ {A} (\overline {{\mathbf {x}}} _ {, 2} ^ {A}) ^ {T} \} - \gamma_ {2} ^ {A} \mathbf {I} \right], \quad \mathbf {H} ^ {B} = \left(h ^ {1} \Gamma_ {B} + h ^ {2}\right) \left[ \mathrm{sym} \{\mathbf {t} ^ {B} (\overline {{\mathbf {x}}} _ {, 1} ^ {B}) ^ {T} \} - \gamma_ {1} ^ {B} \mathbf {I} \right],
\mathbf {H} ^ {C} = \left(h ^ {1} \Gamma_ {C} + h ^ {2}\right) \left[ \mathrm{sym} \{\mathbf {t} ^ {C} (\overline {{\mathbf {x}}} _ {, 2} ^ {C}) ^ {T} \} - \gamma_ {2} ^ {C} \mathbf {I} \right], \quad \mathbf {H} ^ {D} = \left(h ^ {1} \Gamma_ {D} + h ^ {2}\right) \left[ \mathrm{sym} \{\mathbf {t} ^ {D} (\overline {{\mathbf {x}}} _ {, 1} ^ {D}) ^ {T} \} - \gamma_ {1} ^ {D} \mathbf {I} \right].
Also define the skew-symmetric matrices
\hat {\mathbf {a}} _ {h} = (h ^ {1} \Gamma_ {A} + h ^ {2}) \hat {\mathbf {t}} ^ {A}, \quad \hat {\mathbf {b}} _ {h} = (h ^ {1} \Gamma_ {B} + h ^ {2}) \hat {\mathbf {t}} ^ {B},
\hat {\mathbf {c}} _ {h} = (h ^ {1} \Gamma_ {C} + h ^ {2}) \hat {\mathbf {t}} ^ {C}, \quad \hat {\mathbf {d}} _ {h} = (h ^ {1} \Gamma_ {D} + h ^ {2}) \hat {\mathbf {t}} ^ {D}.
Then { \bf K } _ { h } is written
\mathbf {K} _ {h} =
\frac {1}{4} \left[ \begin{array}{c c c c c c c} \mathbf {0} & - \hat {\mathbf {a}} _ {h} + \hat {\mathbf {b}} _ {h} & \mathbf {0} & \hat {\mathbf {b}} _ {h} & \mathbf {0} & \mathbf {0} & \mathbf {0} & - \hat {\mathbf {a}} _ {h} \\ \hat {\mathbf {a}} _ {h} - \hat {\mathbf {b}} _ {h} & - \mathbf {H} ^ {A} + \mathbf {H} ^ {B} & \hat {\mathbf {b}} _ {h} & \mathbf {H} ^ {B} & \mathbf {0} & \mathbf {0} & - \hat {\mathbf {a}} _ {h} & - \mathbf {H} ^ {A} \\ \mathbf {0} & - \hat {\mathbf {b}} _ {h} & \mathbf {0} & - \hat {\mathbf {b}} _ {h} + \hat {\mathbf {c}} _ {h} & \mathbf {0} & \hat {\mathbf {c}} _ {h} & \mathbf {0} & \mathbf {0} \\ - \hat {\mathbf {b}} _ {h} & \mathbf {H} ^ {B} & \hat {\mathbf {b}} _ {h} - \hat {\mathbf {c}} _ {h} & \mathbf {H} ^ {B} + \mathbf {H} ^ {C} & \hat {\mathbf {c}} _ {h} & \mathbf {H} ^ {C} & \mathbf {0} & \mathbf {0} \\ \mathbf {0} & \mathbf {0} & \mathbf {0} & - \hat {\mathbf {c}} _ {h} & \mathbf {0} & - \hat {\mathbf {c}} _ {h} + \hat {\mathbf {d}} _ {h} & \mathbf {0} & \hat {\mathbf {d}} _ {h} \\ \mathbf {0} & \mathbf {0} & - \hat {\mathbf {c}} _ {h} & \mathbf {H} ^ {C} & \hat {\mathbf {c}} _ {h} - \hat {\mathbf {d}} _ {h} & \mathbf {H} ^ {C} - \mathbf {H} ^ {D} & \hat {\mathbf {d}} _ {h} & - \mathbf {H} ^ {D} \\ \mathbf {0} & \hat {\mathbf {a}} _ {h} & \mathbf {0} & \mathbf {0} & \mathbf {0} & - \hat {\mathbf {d}} _ {h} & \mathbf {0} & - \hat {\mathbf {d}} _ {h} + \hat {\mathbf {a}} _ {h} \\ \hat {\mathbf {a}} _ {h} & - \mathbf {H} ^ {A} & \mathbf {0} & \mathbf {0} & - \hat {\mathbf {d}} _ {h} & - \mathbf {H} ^ {D} & \hat {\mathbf {d}} _ {h} - \hat {\mathbf {a}} _ {h} & - \mathbf {H} ^ {D} - \mathbf {H} ^ {A} \end{array} \right].
Note that once the matrix entries in { \bf K } _ { h } are defined, { \bf K } _ { h } is filled just as { \bf K } _ { s } .
The initial stress contribution from the homogeneous part consists of two terms, one from the assumed strain formulation (evaluated at \xi = \eta = 0 ) as detailed earlier, and the other from the crop circle mode addition. These two terms can be written
A _ {0} \mathbf {q} | _ {0} \cdot \mathrm{d} \delta \bar {\pmb {\gamma}} _ {0} = A _ {0} \delta \mathbf {u} \cdot \mathbf {K} _ {s} | _ {0} \cdot \mathrm{d} \mathbf {u} + A _ {0} (q ^ {1} c _ {\xi} + q ^ {2} c _ {\eta}) \mathrm{d} \delta \gamma_ {c c},
where { \bf q } | _ { 0 } and { \bf K } _ { s } | _ { 0 } are the shear force and matrix { \bf K } _ { s } evaluated at the element center. The matrix expression for ( q ^ { 1 } c _ { \xi } + q ^ { 2 } c _ { \eta } ) \mathrm { d } \delta \gamma _ { c c } is analogous to { \bf K } _ { h } from the stabilization terms.
In-plane displacement hourglass control
The in-plane displacement hourglass control is applied in the same way as in the ABAQUS membrane elements. The hourglass strains are defined by
Elements
z _ {\alpha} = \frac {\partial \overline {{\mathbf {x}}}}{\partial S _ {\alpha}} \cdot \overline {{\mathbf {x}}} ^ {I} \gamma^ {I} - \frac {\partial \overline {{\mathbf {X}}}}{\partial S _ {\alpha}} \cdot \overline {{\mathbf {X}}} ^ {I} \gamma^ {I} = \frac {\partial \overline {{\mathbf {x}}}}{\partial S _ {\alpha}} \cdot \overline {{\mathbf {x}}} ^ {I} \gamma^ {I} - \mathbf {T} _ {\alpha} \cdot \overline {{\mathbf {X}}} ^ {I} \gamma^ {I},
where \gamma ^ { I } is the hourglass mode. This mode is obtained by making the "regular" hourglass mode { \Gamma } ^ { I } = \{ 1 , - 1 , 1 , - 1 \} orthogonal to the homogeneous deformation mode in the undeformed shape of the element. This last condition can be written as
\gamma^ {I} = \Gamma^ {I} - \Gamma^ {J} \mathbf {X} ^ {J} \cdot \mathbf {T} _ {\beta} \frac {\partial N ^ {I}}{\partial S _ {\beta}}.
Observe that
\mathbf {T} _ {\alpha} \cdot \overline {{\mathbf {X}}} ^ {I} \gamma^ {I} = \mathbf {T} _ {\alpha} \cdot \overline {{\mathbf {X}}} ^ {I} \Gamma^ {I} - \Gamma^ {J} \overline {{\mathbf {X}}} ^ {J} \cdot \mathbf {T} _ {\beta} \frac {\partial \overline {{\mathbf {X}}}}{\partial S _ {\beta}} \cdot \mathbf {T} _ {\alpha} = \mathbf {T} _ {\alpha} \cdot \overline {{\mathbf {X}}} ^ {I} \Gamma^ {I} - \Gamma^ {J} \overline {{\mathbf {X}}} ^ {J} \cdot \mathbf {T} _ {\alpha} = 0,
and consequently
z _ {\alpha} = \frac {\partial \mathbf {x}}{\partial S _ {\alpha}} \cdot \mathbf {x} ^ {I} \gamma^ {I}.
This expression can be worked out further. We define the projected nodal coordinates
S _ {\alpha} ^ {I} \stackrel {\mathrm{def}} {=} \mathbf {T} _ {\alpha} \cdot \overline {{\mathbf {X}}} ^ {I}
and the projected element area
A = \frac {1}{2} \big [ (S _ {1} ^ {3} - S _ {1} ^ {1}) (S _ {2} ^ {4} - S _ {2} ^ {2}) + (S _ {1} ^ {2} - S _ {1} ^ {4}) (S _ {2} ^ {3} - S _ {2} ^ {1}) \big ].
The hourglass mode can then be written in the form
\begin{array}{l} \gamma^ {1} = \frac {1}{A} \big [ S _ {1} ^ {2} (S _ {2} ^ {3} - S _ {2} ^ {4}) + S _ {1} ^ {3} (S _ {2} ^ {4} - S _ {2} ^ {2}) + S _ {1} ^ {4} (S _ {2} ^ {2} - S _ {2} ^ {3}) \big ], \\ \gamma^ {2} = \frac {1}{A} \big [ S _ {1} ^ {3} (S _ {2} ^ {1} - S _ {2} ^ {4}) + S _ {1} ^ {4} (S _ {2} ^ {3} - S _ {2} ^ {1}) + S _ {1} ^ {1} (S _ {2} ^ {4} - S _ {2} ^ {3}) \big ], \\ \gamma^ {3} = \frac {1}{A} \big [ S _ {1} ^ {4} (S _ {2} ^ {1} - S _ {2} ^ {2}) + S _ {1} ^ {1} (S _ {2} ^ {2} - S _ {2} ^ {4}) + S _ {1} ^ {2} (S _ {2} ^ {4} - S _ {2} ^ {1}) \big ], \\ \gamma^ {4} = \frac {1}{A} \big [ S _ {1} ^ {1} (S _ {2} ^ {3} - S _ {2} ^ {2}) + S _ {1} ^ {2} (S _ {2} ^ {1} - S _ {2} ^ {3}) + S _ {1} ^ {3} (S _ {2} ^ {2} - S _ {2} ^ {1}) \big ]. \\ \end{array}
The hourglass stiffness is chosen equal to
K _ {h} = (r _ {F} G) \frac {\partial N ^ {I}}{\partial S _ {\beta}} \frac {\partial N ^ {I}}{\partial S _ {\beta}} h A,
where G is the shear modulus and r _ { F } is a small number chosen to be 0.005 in ABAQUS/Standard and
Elements
0.05 in ABAQUS/Explicit. The hourglass force Z conjugate to z is then equal to
Z _ {\alpha} = K _ {h} z _ {\alpha}.
For virtual work we need the first variation of the hourglass strain. From the expression for the strain follows immediately
\delta z _ {\alpha} = \frac {\partial \overline {{\mathbf {x}}}}{\partial S _ {\alpha}} \cdot \delta \overline {{\mathbf {x}}} ^ {I} \gamma^ {I} + \frac {\partial \delta \overline {{\mathbf {x}}}}{\partial S _ {\alpha}} \cdot \overline {{\mathbf {x}}} ^ {I} \gamma^ {I}.
Note that the second term vanishes in the initial configuration since \overline { { \mathbf { X } } } ^ { I } \gamma ^ { I } = 0 . The second variation is needed for the Jacobian. From the first variation follows right away
\mathrm{d} \delta z _ {\alpha} = \frac {\partial \mathrm{d} \overline {{\mathbf {x}}}}{\partial S _ {\alpha}} \cdot \delta \overline {{\mathbf {x}}} ^ {I} \gamma^ {I} + \frac {\partial \delta \overline {{\mathbf {x}}}}{\partial S _ {\alpha}} \cdot \mathrm{d} \overline {{\mathbf {x}}} ^ {I} \gamma^ {I}.
The second variation does not contribute in the initial configuration since initially Z _ { \alpha } = 0
Rotational hourglass control
The expressions for the curvature change, the transverse shear constraints, and the drilling mode constraints still leave three nonhomogeneous rotational modes unconstrained. These modes correspond to zero rotation at the midedges and zero gradient at the centroid. Hence, they correspond to the familiar \Gamma ^ { I } = \{ 1 , - 1 , 1 , - 1 \} hourglass pattern. To pass curvature patch tests exactly, it is necessary to use orthogonalized hourglass patterns as derived for in-plane hourglass control.
This last aspect implies that the rotational hourglass mode corresponds to the mixed derivative of the rotation at the centroid:
\gamma^ {I} \Delta \phi^ {I} = 4 \times \frac {\partial^ {2} \Delta \phi}{\partial \xi_ {1} \partial \xi_ {2}}.
We cannot use the above formulation directly in a formulation suitable for multiple finite rotation increments. Hence, we use the same approach as for the calculation of the curvature change. For the purpose of the calculation we define the updated shell direction vectors
\tilde {\mathbf {t}} _ {i} ^ {t + \Delta t} = \Delta \boldsymbol {q} \mathbf {t} _ {i} ^ {t} \Delta \boldsymbol {q} ^ {\dagger}.
The updated shell direction vectors do not actually have to be stored: they are used only for the derivation of the expression for the hourglass strain. We now formally define the hourglass strain tensor as
\vartheta_ {i j} = \mathbf {t} _ {i} \cdot \frac {\partial^ {2} \mathbf {t} _ {j}}{\partial \xi_ {1} \partial \xi_ {2}}.
Elements
Observe that
\frac {\partial^ {2} (\mathbf {t} _ {i} \cdot \mathbf {t} _ {j})}{\partial \xi_ {1} \partial \xi_ {2}} = \mathbf {t} _ {i} \cdot \frac {\partial^ {2} \mathbf {t} _ {j}}{\partial \xi_ {1} \partial \xi_ {2}} + \frac {\partial^ {2} \mathbf {t} _ {i}}{\partial \xi_ {1} \partial \xi_ {2}} \cdot \mathbf {t} _ {j} + \frac {\partial \mathbf {t} _ {i}}{\partial \xi_ {1}} \cdot \frac {\partial \mathbf {t} _ {j}}{\partial \xi_ {2}} + \frac {\partial \mathbf {t} _ {i}}{\partial \xi_ {2}} \cdot \frac {\partial \mathbf {t} _ {j}}{\partial \xi_ {1}} = 0.
For the purpose of hourglass strain calculation we assume that all products of first-order derivatives with respect to \xi _ { 1 } and \xi _ { 2 } can be neglected. Consequently,
\frac {\partial^ {2} (\mathbf {t} _ {i} \cdot \mathbf {t} _ {j})}{\partial \xi_ {1} \partial \xi_ {2}} = \mathbf {t} _ {i} \cdot \frac {\partial^ {2} \mathbf {t} _ {j}}{\partial \xi_ {1} \partial \xi_ {2}} + \frac {\partial^ {2} \mathbf {t} _ {i}}{\partial \xi_ {1} \partial \xi_ {2}} \cdot \mathbf {t} _ {j} = 0;
and, hence,
\vartheta_ {i j} = - \vartheta_ {j i}
is skew-symmetric. Observe that the mixed derivative of \mathbf { t } _ { i } can be expressed in terms of the hourglass strain tensor with
\mathbf {t} _ {i} \vartheta_ {i j} = \mathbf {t} _ {i} \mathbf {t} _ {i} \cdot \frac {\partial^ {2} \mathbf {t} _ {j}}{\partial \xi_ {1} \partial \xi_ {2}} = \frac {\partial^ {2} \mathbf {t} _ {j}}{\partial \xi_ {1} \partial \xi_ {2}}.
In the undeformed configuration, we assume that \vartheta _ { i j } = 0 . Subsequent values of \vartheta _ { i j } are obtained incrementally. From the expression for \tilde { \mathbf { t } } _ { i } ^ { t + \Delta t } we obtain
\begin{array}{l} \frac {\partial^ {2} \tilde {\mathbf {t}} _ {i} ^ {t + \Delta t}}{\partial \xi_ {1} \partial \xi_ {2}} = \frac {\partial^ {2} \Delta \boldsymbol {q}}{\partial \xi_ {1} \partial \xi_ {2}} \mathbf {t} _ {i} ^ {t} \Delta \boldsymbol {q} ^ {\dagger} + \Delta \boldsymbol {q} \mathbf {t} _ {i} ^ {t} \frac {\partial \Delta \boldsymbol {q} ^ {\dagger}}{\partial \xi_ {1} \partial \xi_ {2}} + \Delta \boldsymbol {q} \frac {\partial \mathbf {t} _ {i} ^ {t}}{\partial \xi_ {1} \partial \xi_ {2}} \Delta \boldsymbol {q} ^ {\dagger} + \\ \frac {\partial \Delta \boldsymbol {q}}{\partial \xi_ {1}} \mathbf {t} _ {i} ^ {t} \frac {\partial \Delta \boldsymbol {q} ^ {\dagger}}{\partial \xi_ {2}} + \frac {\partial \Delta \boldsymbol {q}}{\partial \xi_ {2}} \mathbf {t} _ {i} ^ {t} \frac {\partial \Delta \boldsymbol {q} ^ {\dagger}}{\partial \xi_ {1}} + \frac {\partial \Delta \boldsymbol {q}}{\partial \xi_ {1}} \frac {\partial \mathbf {t} _ {i} ^ {t}}{\partial \xi_ {2}} \Delta \boldsymbol {q} ^ {\dagger} + \\ \frac {\partial \Delta \pmb {q}}{\partial \xi_ {2}} \frac {\partial \mathbf {t} _ {i} ^ {t}}{\partial \xi_ {1}} \Delta \pmb {q} ^ {\dagger} + \Delta \pmb {q} \frac {\partial \mathbf {t} _ {i} ^ {t}}{\partial \xi_ {1}} \frac {\partial \Delta \pmb {q} ^ {\dagger}}{\partial \xi_ {2}} + \Delta \pmb {q} \frac {\partial \mathbf {t} _ {i} ^ {t}}{\partial \xi_ {2}} \frac {\partial \Delta \pmb {q} ^ {\dagger}}{\partial \xi_ {1}}. \\ \end{array}
In this expression we also ignore all terms with products of derivatives with respect to \xi _ { 1 } and \xi _ { 2 } . Hence, the above expression simplifies to
\frac {\partial \tilde {\mathbf {t}} _ {i} ^ {t + \Delta t}}{\partial \xi_ {1} \partial \xi_ {2}} = \frac {\partial^ {2} \Delta \pmb {q}}{\partial \xi_ {1} \partial \xi_ {2}} \mathbf {t} _ {i} ^ {t} \Delta \pmb {q} ^ {\dagger} + \Delta \pmb {q} \mathbf {t} _ {i} ^ {t} \frac {\partial \Delta \pmb {q} ^ {\dagger}}{\partial \xi_ {1} \partial \xi_ {2}} + \Delta \pmb {q} \frac {\partial \mathbf {t} _ {i} ^ {t}}{\partial \xi_ {1} \partial \xi_ {2}} \Delta \pmb {q} ^ {\dagger}.
The second term on the right-hand side can be written in the form
\Delta \pmb {q} \mathbf {t} _ {i} ^ {t} \frac {\partial \Delta \pmb {q} ^ {\dagger}}{\partial \xi_ {1} \partial \xi_ {2}} = \left[ \frac {\partial \Delta \pmb {q}}{\partial \xi_ {1} \partial \xi_ {2}} \left(\mathbf {t} _ {i} ^ {t}\right) ^ {\dagger} \Delta \pmb {q} ^ {\dagger} \right] ^ {\dagger} = - \left[ \frac {\partial \Delta \pmb {q}}{\partial \xi_ {1} \partial \xi_ {2}} \mathbf {t} _ {i} ^ {t} \Delta \pmb {q} ^ {\dagger} \right] ^ {\dagger}.
Hence, the scalar parts of the first two terms cancel each other and the vector parts reinforce each other, leading to
Elements
\frac {\partial \tilde {\mathbf {t}} _ {i} ^ {t + \Delta t}}{\partial \xi_ {1} \partial \xi_ {2}} = 2 \mathbf {V} \left(\frac {\partial \Delta \pmb {q}}{\partial \xi_ {1} \partial \xi_ {2}} \mathbf {t} _ {i} ^ {t} \Delta \pmb {q} ^ {\dagger}\right) + \Delta \pmb {q} \frac {\partial \mathbf {t} _ {i} ^ {t}}{\partial \xi_ {1} \partial \xi_ {2}} \Delta \pmb {q} ^ {\dagger}.
The inverse of a rotation quaternion such as \Delta \pmb q is equal to its conjugate ( \Delta \pmb q ^ { - 1 } = \Delta \pmb q ^ { \dag } ) ; ; hence, we can write
\begin{array}{l} \frac {\partial \tilde {\mathbf {t}} _ {i} ^ {t + \Delta t}}{\partial \xi_ {1} \partial \xi_ {2}} = 2 \mathbf {V} \left(\frac {\partial \Delta \pmb {q}}{\partial \xi_ {1} \partial \xi_ {2}} \Delta \pmb {q} ^ {\dagger} \Delta \pmb {q} \mathbf {t} _ {i} ^ {t} \Delta \pmb {q} ^ {\dagger}\right) + \Delta \pmb {q} \frac {\partial \mathbf {t} _ {i} ^ {t}}{\partial \xi_ {1} \partial \xi_ {2}} \Delta \pmb {q} ^ {\dagger} \\ = \mathbf {V} \left(2 \frac {\partial \Delta \pmb {q}}{\partial \xi_ {1} \partial \xi_ {2}} \Delta \pmb {q} ^ {\dagger} \tilde {\mathbf {t}} _ {i} ^ {t + \Delta t}\right) + \Delta \pmb {q} \frac {\partial \mathbf {t} _ {i} ^ {t}}{\partial \xi_ {1} \partial \xi_ {2}} \Delta \pmb {q} ^ {\dagger} \\ = \Delta \mathbf {R} _ {\Gamma} \times \tilde {\mathbf {t}} _ {i} ^ {t + \Delta t} + \Delta \pmb {q} \frac {\partial \mathbf {t} _ {i} ^ {t}}{\partial \xi_ {1} \partial \xi_ {2}} \Delta \pmb {q} ^ {\dagger}, \\ \end{array}
where we have formally defined the incremental hourglass update vector
\Delta \mathbf {R} _ {\Gamma} \stackrel {\mathrm{def}} {=} \mathbf {V} \left(2 \frac {\partial \Delta \pmb {q}}{\partial \xi_ {1} \partial \xi_ {2}} \Delta \pmb {q} ^ {\dagger}\right),
which must be expressed in terms of the incremental rotation hourglass mode. From the definition of the incremental quaternion \Delta \pmb q follows, while neglecting the products of \xi _ { 1 } and \xi _ { 2 } derivatives:
\frac {\partial \Delta \pmb {q}}{\partial \xi_ {1} \partial \xi_ {2}} = \left(\frac {1}{2} \sin \frac {\Delta \phi}{2} \frac {\partial \Delta \phi}{\partial \xi_ {1} \partial \xi_ {2}}, \frac {1}{2} \cos \frac {\Delta \phi}{2} \mathbf {p} \frac {\partial \Delta \phi}{\partial \xi_ {1} \partial \xi_ {2}} + \sin \frac {\Delta \phi}{2} \frac {\partial \mathbf {p}}{\partial \xi_ {1} \partial \xi_ {2}}\right);
thus, for \Delta { \bf R } _ { \Gamma } again with use of the incremental quaternion definition
\Delta \mathbf {R} _ {\Gamma} = \mathbf {p} \frac {\partial \Delta \phi}{\partial \xi_ {1} \partial \xi_ {2}} + \sin \Delta \phi \frac {\partial \mathbf {p}}{\partial \xi_ {1} \partial \xi_ {2}} + (1 - \cos \Delta \phi) \mathbf {p} \times \frac {\partial \mathbf {p}}{\partial \xi_ {1} \partial \xi_ {2}}.
From the definition of \Delta \phi and p follows, again neglecting the products of first derivatives
\frac {\partial \Delta \phi}{\partial \xi_ {1} \partial \xi_ {2}} = \frac {1}{\Delta \phi} \Delta \phi \cdot \frac {\partial \Delta \phi}{\partial \xi_ {1} \partial \xi_ {2}},
\frac {\partial \mathbf {p}}{\partial \xi_ {1} \partial \xi_ {2}} = \frac {1}{\Delta \phi} \frac {\partial \Delta \pmb {\phi}}{\partial \xi_ {1} \partial \xi_ {2}} - \frac {1}{\Delta \phi^ {3}} \Delta \pmb {\phi} \Delta \pmb {\phi} \cdot \frac {\partial \Delta \pmb {\phi}}{\partial \xi_ {1} \partial \xi_ {2}} = \frac {1}{\Delta \phi} \left[ \frac {\partial \Delta \pmb {\phi}}{\partial \xi_ {1} \partial \xi_ {2}} - \mathbf {p} \mathbf {p} \cdot \frac {\partial \Delta \pmb {\phi}}{\partial \xi_ {1} \partial \xi_ {2}} \right].
After substitution in the expression for \Delta \mathbf { R } _ { \Gamma } and some algebra one obtains
\Delta \mathbf {R} _ {\Gamma} = \frac {\sin \Delta \phi}{\Delta \phi} \frac {\partial \Delta \pmb {\phi}}{\partial \xi_ {1} \partial \xi_ {2}} + \frac {1 - \cos \Delta \phi}{\Delta \phi} \mathbf {p} \times \frac {\partial \Delta \pmb {\phi}}{\partial \xi_ {1} \partial \xi_ {2}} + \left(1 - \frac {\sin \Delta \phi}{\Delta \phi}\right) \mathbf {p} \mathbf {p} \cdot \frac {\partial \Delta \pmb {\phi}}{\partial \xi_ {1} \partial \xi_ {2}}.
Note that ¢R¡ ! \Delta { \bf R } _ { \Gamma } \frac { \partial \Delta \phi } { \partial \xi _ { 1 } \partial \xi _ { 2 } } @¢Á when \Delta \phi \to 0 . For the updated hourglass tensor one readily obtains
Elements
\begin{array}{l} \vartheta_ {i j} ^ {t + \Delta t} = \tilde {\mathbf {t}} _ {i} ^ {t + \Delta t} \cdot \frac {\partial^ {2} \tilde {\mathbf {t}} _ {j} ^ {t + \Delta t}}{\partial \xi_ {1} \partial \xi_ {2}} = \tilde {\mathbf {t}} _ {i} ^ {t + \Delta t} \cdot \left[ \Delta \mathbf {R} _ {\Gamma} \times \tilde {\mathbf {t}} _ {j} ^ {t + \Delta t} + \Delta \pmb {q} \frac {\partial^ {2} \mathbf {t} _ {j} ^ {t}}{\partial \xi_ {1} \partial \xi_ {2}} \Delta \pmb {q} ^ {\dagger} \right] \\ = - \Delta \mathbf {R} _ {\Gamma} \cdot (\tilde {\mathbf {t}} _ {i} ^ {t + \Delta t} \times \tilde {\mathbf {t}} _ {j} ^ {t + \Delta t}) + \tilde {\mathbf {t}} _ {i} ^ {t + \Delta t} \cdot \left(\Delta \pmb {q} \mathbf {t} _ {k} ^ {t} \vartheta_ {k j} ^ {t} \Delta \pmb {q} ^ {\dagger}\right) \\ = - \Delta \mathbf {R} _ {\Gamma} \cdot (\tilde {\mathbf {t}} _ {i} ^ {t + \Delta t} \times \tilde {\mathbf {t}} _ {j} ^ {t + \Delta t}) + \tilde {\mathbf {t}} _ {i} ^ {t + \Delta t} \cdot \tilde {\mathbf {t}} _ {k} ^ {t + \Delta t} \vartheta_ {k j} ^ {t} = - \Delta \mathbf {R} _ {\Gamma} \cdot (\tilde {\mathbf {t}} _ {i} ^ {t + \Delta t} \times \tilde {\mathbf {t}} _ {j} ^ {t + \Delta t}) + \vartheta_ {i j} ^ {t}. \\ \end{array}
This expression simplifies further with the introduction of the hourglass vector
\vartheta_ {i} = - \frac {1}{2} e _ {i j k} \vartheta_ {j k},
which yields the update formula
\vartheta_ {i} ^ {t + \Delta t} = \Delta \mathbf {R} _ {\Gamma} \cdot \tilde {\mathbf {t}} _ {i} ^ {t + \Delta t} + \vartheta_ {i} ^ {t}.
The first and second variation are obtained in entirely the same way as the first and second variation of the curvature change. For the first variation we neglect terms of order \Delta \mathbf { R } _ { \Gamma } and obtain
\delta \vartheta_ {i} = \delta \Delta \mathbf {R} _ {\Gamma} \cdot \mathbf {t} _ {i} + \Delta \mathbf {R} _ {\Gamma} \cdot \delta \mathbf {t} _ {i} = \delta \Delta \mathbf {R} _ {\Gamma} \cdot \mathbf {t} _ {i} = \mathbf {t} _ {i} \cdot \frac {\partial^ {2} \delta \pmb {\phi}}{\partial \xi_ {1} \partial \xi_ {2}}.
For the second variation we ignore in addition the terms of order \mathbf { d R } _ { \Gamma } and \delta \mathbf { R } _ { \Gamma } with as final result
\mathrm{d} \delta \vartheta_ {i} = 0.
Degenerate elements
In general meshes it will be desirable to collapse at least some of the quadrilateral elements to triangles or to use the triangular element S3 or S3R, which is in fact an internally collapsed S4R element. For this case the calculation of the membrane strains and the curvature changes proceeds along the same lines as before. The transverse shears will now be zero at the degenerate edges. Finally, calculation of all hourglass constraints will be omitted.
Rotary inertia scaling for explicit dynamics
For numerical efficiency in explicit dynamic analysis, it is desirable to have the stable time increment determined by the membrane response of the structure. For this reason scaling of the rotary inertia based on the element's reference geometry is included in ABAQUS/Explicit.
In explicit dynamic analyses the stable time increment is proportional to the inverse of the highest frequency of the element. Therefore, we must ensure that the highest frequency associated with the transverse shear response does not exceed the highest frequency associated with the membrane response. For thick elements (that is, for elements whose thickness is order unity relative to a characteristic length in the element), the membrane frequencies are dominant. The primary consideration in choosing appropriate scalings is that in the limit as the element's thickness goes to
Elements
zero, the transverse shear frequencies remain below the membrane ones. Recall that for a one-dimensional spring-mass oscillator, the natural frequency ! can be written in terms of the stiffness K and the mass M as
\omega = \sqrt {\frac {K}{M}}.
For the transverse shear response the rotary inertia, which is proportional to the cube of the thickness, plays the role of the mass of the system. All other quantities--the membrane stiffness, the mass associated with membrane deformation, and the transverse shear stiffness--are proportional to the thickness. Hence as the thickness of the element goes to zero, the frequencies associated with transverse shear go to infinity proportional to the inverse of the thickness, while the membrane frequencies remain constant. Without scaling, the stable time increment would go to zero as the thickness becomes small.
Rotary inertia scaling
For thin elements the rotary inertia is small (negligible) relative to the rotational inertia of the mass at the nodes rotating about an axis through the element center. Therefore, we choose a scaling on the rotary inertia such that it never becomes smaller than a fixed (small) percentage of the rotational inertia of the mass at the nodes rotating about an axis through the center of the element.
Let R be the nondimensional rotary inertia scaling, where R \geq 1 . When R = 1, the true rotary inertia is used. Consider a lumped mass matrix for a 4-node element, and let the element be flat and square. For rotations about an axis in the plane of the element, parallel to an element edge, passing through the center, the contribution to the rotational inertia of the mass at the nodes is
\frac {1}{4} \rho A h L ^ {2},
where A is the area of the element, L is the edge length, and \rho is the mass density. The sum of the rotary inertia at the four nodes is
R \rho A \frac {h ^ {3}}{1 2}.
The ratio of the in-plane contribution to the rotary contribution is
\frac {3 L ^ {2}}{R h ^ {2}}.
For the rotary inertia to remain a fixed fraction--say ²--of the mass contribution as the thickness goes to zero, asymptotically R must be proportional to L ^ { 2 } / h ^ { 2 } ; that is,
R \rightarrow \epsilon \frac {3 L ^ {2}}{h ^ {2}}.
Elements
For planar geometries with element directors along the normal direction, closed-form expressions are possible for the highest membrane and transverse shear frequencies. In such cases the length parameter L is interpreted as a characteristic element length that depends on the element distortion. To handle arbitrarily shaped curved elements, exact calculation of the element frequencies becomes difficult. However, we can safely bound the frequencies by an appropriate choice of L in the following scaling:
R = \left[ 1 6 + \epsilon^ {4} \left(\frac {3 L ^ {2}}{h ^ {2}}\right) ^ {4} \right] ^ {\frac {1}{4}},
where \epsilon = 1 / 4 . For triangular elements the characteristic length, L; is the minimum element edge length; and for quadrilateral elements,
L = \frac {2 \| \mathbf {x} _ {, \xi} \times \mathbf {x} _ {, \eta} \|}{\max (\| \mathbf {x} _ {, \xi} \| , \| \mathbf {x} _ {, \eta} \|)}.
The factor 16 in the definition of R is used to protect against bending frequencies determining the stable time increment in very fine meshes subjected to loads that cause an increase in thickness of the shell.
Rotational bulk viscosity for explicit dynamics
For the displacement degrees of freedom, bulk viscosity introduces damping associated with volumetric straining. Linear bulk viscosity or truncation frequency damping is used to damp the high frequency ringing that leads to unwanted noise in the solution or spurious overshoot in the response amplitude. For the same reason, in shells we need to damp the high frequency ringing in the rotational degrees of freedom with linear bulk viscosity acting on the mean curvature strain rate. This damping generates a bulk viscosity "pressure moment," m, which is linear in the mean curvature strain rate:
m = b \frac {h _ {0} ^ {2}}{1 2} \rho c _ {d} L \frac {\Delta \kappa}{\Delta t},
where b is a damping coefficient (default = 0.06), h _ { 0 } is the original thickness, \rho is the mass density, c _ { d } is the current dilatational wave speed, L is the characteristic length used for rotary inertia scaling, and \Delta \kappa = \Delta \kappa _ { 1 1 } + \Delta \kappa _ { 2 2 } is twice the increment in mean curvature. The dilatational wave speed is given in terms of the effective Lamé's constants as
c _ {d} = \sqrt {\frac {(\lambda + 2 \mu)}{\rho}}.
The resultant pressure moment m h , , where h is the current thickness, is added to the direct components of the moment resultant.
Fully integrated finite-membrane-strain shell formulation
Elements
Element S4 is a fully integrated finite-membrane-strain shell element. Since the element's stiffness is fully integrated, no spurious membrane or bending zero energy modes exist and no membrane or bending mode hourglass stabilization is used. Drill rotation control, however, is required. Element S4 uses the same drill stiffness formulation as used for element S4R. Similarly, element S4 assumes that the transverse shear strain (and force, since the transverse shear treatment is elastic based on the initial elastic modulus of the material) is constant over the element. Therefore, all four stiffness integration locations will have the same transverse shear strain, transverse shear section force, and transverse shear stress distribution. The transverse shear treatment for S4 is identical to that for S4R.
It is well known that a standard displacement formulation will exhibit shear locking for applications dominated by in-plane bending deformation. However, a standard displacement formulation for the out-of-plane bending stiffness is not subject to similar locking response. Hence, S4 uses a standard displacement formulation for the element's bending stiffness, and the theory presented above for the rotation kinematics and bending strain measures applies to S4. The primary difference between the element formulations for S4 and S4R is the treatment of the membrane strain field. This formulation is the topic of the following discussion.
The membrane formulation used for S4 does not rely on the fact that S4 is a shell element. Hence, the discussion below details the formulation from the point of view that the membrane response is governed by the equilibrium for a three-dimensional body in a state of plane stress.
Consider an enhancement \tilde { \mathbf { D } } to the rate of deformation tensor D. We introduce the enhanced rate of deformation tensor, \bar { \mathbf { D } } , as
\bar {\mathbf {D}} \stackrel {\mathrm{def}} {=} \mathbf {D} + \tilde {\mathbf {D}},
where \tilde { \bf D } is defined subsequently.
Admissible variations in the rate of deformation are also introduced as
\bar {\delta} \mathbf {D} \stackrel {\text { def }} {=} \delta \mathbf {D} + \delta \tilde {\mathbf {D}},
where
\delta \mathbf {D} \stackrel {\mathrm{def}} {=} \mathrm{sym} \left(\frac {\partial \delta \mathbf {x}}{\partial \mathbf {x}}\right).
We now introduce constraints on the enhancements D and \tilde { \mathbf { D } } :
\int_ {V} \left(\bar {\mathbf {D}} - \mathbf {D}\right): \stackrel {*} {\pmb {\sigma}} d V = \int_ {V} \left(\delta \bar {\mathbf {D}} - \delta \mathbf {D}\right): \stackrel {*} {\pmb {\sigma}} d V = 0
so that the modified virtual work statement can be written in the form
\int_ {V} \pmb {\sigma} (\bar {\mathbf {D}}): \delta \bar {\mathbf {D}} d V + \int_ {V} \mathbf {b} \cdot \delta \mathbf {x} d V - \int_ {S _ {T}} \mathbf {t} \cdot \delta \mathbf {x} d S _ {T} = 0,