Elastoplastic law of Cosserat type in shell theory with drilling rotation

Abstract

Within the framework of six-parameter non-linear shell theory, with strain measures of the Cosserat type, we develop small-strain J₂-type elastoplastic constitutive relations. The relations are obtained from the Cosserat plane stress relations assumed in each shell layer, by through-the-thickness integration employing the first-order shear theory. The formulation allows for unlimited translations and rotations. The constitutive relations are expressed in terms of the following material data: Young’s modulus, Poisson’s ratio, Cosserat shear modulus, characteristic length and material parameters pertaining to plasticity. Numerical examples are presented to show the correctness of the proposed model and capacity to analyze regular and irregular shells.

Keywords

Cosserat elastoplasticity Cosserat surface drilling rotation shell finite elements

1. Introduction

Materials with internal structure constitute a vast class of bodies [1, 2] which possess, in comparison with the Cauchy continuum, a more advanced mathematical constitution that describes the three main ingredients of theory: equilibrium equations, kinematical relations and material law. It is generally acknowledged that the foundations of this class of bodies were laid down by the Cosserat brothers [3]. A comprehensive historical outline can be found in the literature [1, 2, 4 –6]. The structure of the equilibrium equations for three-dimensional (3D) Cosserat solids, 3D structural shells and 3D structural rods seems to be well known. The same is true for the kinematic relations between strains and generalized displacement fields, both in linear and non-linear formulations [1, 2, 4, 5, 7, 8 –13].

However, the material relations within the realm of the 3D Cosserat theory are still an open question. This is attributed mainly to the fact that the formulation requires additional (as compared with the Cauchy continuum) material parameters. As shown, for instance, in Nowacki [14], Cowin [15] and Eremeyev and Pietraszkiewicz [16] for linear elastic Cosserat materials, the number of elastic constants is equal to six which causes serious problems where identification of their values is concerned. Some recent developments and earlier efforts are discussed for example in Altenbach and Eremeyev [8] and Chróścielewski et al. [17 –19] and references therein.

This paper is concerned with the non-linear six-parameter shell theory, i.e. a Cosserat-type surface naturally endowed with a drilling degree of freedom. The main purpose is the formulation of the elastoplastic constitutive relation consistent with a Cosserat-type surface. The necessary theoretical background pertaining to Cosserat-type surfaces can be found in Libai and Simmonds [20], Chróścielewski et al. [21, 22] and Bîrsan and Neff [23]. Extending our previous work [24] we show details of the elastoplastic constitutive relations that take into account the Cosserat-type structure of the plane stress relation. The salient ingredients of the formulation are:

strains are small everywhere;

at each shell layer the Cosserat plane stress is assumed;

through-the-thickness integration is based on the first-order shear deformation theory (FOSD);

additive decomposition of small elastoplastic strain rate;

associative flow rule;

isotropic hardening.

2. Strain and stress measures

This section presents the framework of the shell theory used, necessary to formulate the constitutive relations. Let $x$ denote a typical point on the undeformed shell base surface, and $M$ and $y (x)$ denote the position of that point in the current configuration. The microstructure of the shell embedded in a fixed Cartesian frame $e_{i}, i = 1, 2, 3$ ( $x = x_{i} e_{i}, i = 1, 2, 3$ ) is defined by the so-called structure tensor $T_{0} \in SO (3)$ satisfying

t_{i}^{0} (x) = T_{0} (x) e_{i} .

(1)

For the purposes of the present study we limit our attention to the theory in which the initial directors $t_{β}^{0} = x,_{β}$ , $β = 1, 2$ , $(.),_{β} = \frac{d (.)}{d s^{β}}$ $(s^{β} are curvilinear arc length coordinates)$ are orthonormal vectors locally tangent to $M$ at $x$ and the normal vector is defined by

t^{0} \equiv t_{3}^{0} = \frac{t_{1}^{0} \times t_{2}^{0}}{| | t_{1}^{0} \times t_{2}^{0} | |},

(2)

where $\times$ stands for the standard vector product. Vector $t_{3}^{0}$ is crucial in the employed shell theory due to this fact the index 3 in equation (2) is sometimes omitted. The microstructure triad $t_{i}^{0} (x)$ satisfies the orthogonality conditions $〈 t_{i}^{0}, t_{j}^{0} 〉 = δ_{ij}$ , where $〈 \cdot, \cdot 〉$ is the standard inner product on $R^{3}$ and $| | \cdot | |$ is the corresponding norm $| | \cdot | | = {〈 \cdot, \cdot 〉}^{\frac{1}{2}}$ . It is assumed that $t_{i}^{0} (x)$ remains rigid in motion $u (x) = (u (x), Q (x))$ , where the field $u (x) = y (x) - x$ denotes the translation vector and the independent field $Q (x) \in SO (3)$ , $t_{i} (y (x)) = Q (x) t_{i}^{0} (x)$ describes the finite rotation of the directors $t_{i}^{0} (x)$ . We wish to emphasize that the rotations are limited only by the local impenetrability condition, i.e. $t_{3}$ is never tangent to the deformed shell base.

In the present approach $Q (x)$ is expressed in terms of the finite rotation vector $z$ through the canonical representation

Q = 1 + a Z + b Z^{2}, Z = {axl}^{- 1} (z) \in so (3) .

(3)

The coefficients in equation (3) are defined as $a = \frac{\sin z}{z}$ , $b = \frac{1 - \cos z}{z^{2}}$ , $z = | | z | |$ and the operator $axl$ defines the axial vector $z$ of the skew tensor $Z$ (i.e. $b = Za = z \times a$ ).

The vectors of Cosserat-type strain measures in the defined basis are given in, for example, Chróścielewski et al. [21, 22] as:

ε_{β} = y,_{β} - t_{β} = u,_{β} + (1 - Q) t_{β}^{0}, δ ε_{β} = v,_{β} + y,_{β} \times w,

(4)

κ_{β} = axl (Q,_{β} Q^{T}), δ κ_{β} = w,_{β} .

(5)

Here $v (x)$ and $w (x)$ are vector fields of virtual translations and virtual rotations, respectively. In components, equations (4)₁ and (5)₁ read

ε_{β} = ε_{β 1} t_{1} + ε_{β 2} t_{2} + ε_{β} t_{3}, κ_{β} = t_{3} \times (κ_{β 1} t_{1} + κ_{β 2} t_{2}) + κ_{β} t_{3} .

(6)

In the Finite Element Method (FEM) approximation they are assembled in a one-column array (vector form)

ε^{M} = {ε_{11}^{M} ε_{22}^{M} ε_{12}^{M} ε_{21}^{M} | ε_{1}^{M} ε_{2}^{M} | | κ_{11}^{M} κ_{22}^{M} κ_{12}^{M} κ_{21}^{M} | κ_{1}^{M} κ_{2}^{M}}^{T} = {ε_{m}^{M} | ε_{s}^{M} | | ε_{b}^{M} | ε_{d}^{M}}^{T},

(7)

where the labels m, s, b, d denote respectively: the membrane, transverse shear, bending/twisting and drilling parts respectively. The superscript M is introduced to draw a distinction between these strains and in-layer strains, which will be shown later.

To equations (4)₁ and (5)₁ we adjoin the energy conjugated stress and couple stress resultants

n^{β} = N^{β 1} t_{1} + N^{β 2} t_{2} + Q^{β} t_{3}, m^{β} = t_{3} \times (M^{β 1} t_{1} + M^{β 2} t_{2}) + M^{β} t_{3},

(8)

which are assembled in the one-column array as well

s = {N^{11} N^{22} N^{12} N^{21} | Q^{1} Q^{2} | | M^{11} M^{22} M^{12} M^{21} | M^{1} M^{2}}^{T} = {s_{m} | s_{s} | | s_{b} | s_{d}}^{T} .

(9)

The column arrays (vectors) in equations (7) and (9) are connected through the matrix constitutive relation of the type $s = C ε^{M}$ $(C_{(12 \times 12)})$ which is the topic of the following sections. The remaining details of the formulation can be found in the literature [19, 21, 22, 25 –27].

3. The elastoplastic constitutive relation

The constitutive relation presented in this section generalizes the idea presented in, for example, De Borst [28]. Assume that in each layer of the shell there exists the Cosserat plane stress for which we define the one-column array (vector) of generalized stresses and strains, respectively, given by

σ = {\begin{matrix} σ_{m} \\ σ_{d} \end{matrix}} = {\begin{matrix} σ_{11} \\ σ_{22} \\ σ_{12} \\ σ_{21} \\ \dots \\ m_{1} \\ m_{2} \end{matrix}}, ε = {\begin{matrix} ε_{m} \\ ε_{d} \end{matrix}} = {\begin{matrix} ε_{11} \\ ε_{22} \\ ε_{12} \\ ε_{21} \\ \dots \\ κ_{1} \\ κ_{2} \end{matrix}},

(10)

where the Cosserat media components are the couple stresses $m_{1}$ , $m_{2}$ and the microcurvatures $κ_{1}$ , $κ_{2}$ .

The material law for the elastic part is written in the following form

σ = C^{e} ε^{e},

(11)

where

C^{e} = [\begin{matrix} {Gc}_{1} & {Gc}_{2} & 0 & 0 & 0 & 0 \\ {Gc}_{2} & {Gc}_{1} & 0 & 0 & 0 & 0 \\ 0 & 0 & G + G_{c} & G - G_{c} & 0 & 0 \\ 0 & 0 & G - G_{c} & G + G_{c} & 0 & 0 \\ 0 & 0 & 0 & 0 & 2 {Gl}^{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & 2 {Gl}^{2} \end{matrix}]

(12)

with $c_{1} = (1 - ν) / (1 - 2 ν)$ and $c_{2} = ν c_{1}$ . Material parameters introduced here are $G$ , $ν$ , $G_{c}$ and $l$ which denote, respectively, the Kirchhoff moduli, Poisson’s ratio, additional shear moduli and the micropolar characteristic length. The structure of the elastic or elastoplastic constitutive matrix can be arranged in a block manner

C = [\begin{matrix} C_{mm} & C_{md} \\ C_{dm} & C_{dd} \end{matrix}] .

(13)

Equation (12) shows the structure of the constitutive matrix valid for purely elastic analysis. Bearing in mind the formulation of the elastoplastic material law, we introduce the term $C_{md}$ in equation (13), which becomes non-zero during the plastic flow.

We further assume that at each point $x$ of the shell reference surface $M$ there exists a non-zero shell thickness $h = h^{+} + h^{-}$ parameterized by the physical coordinate $ζ \in [- h^{-}, + h^{+}]$ . Here $h^{-}$ denotes the position of the bottom face of the 3D shell-like body while $h^{+}$ stands for the position of the upper face.

While the developments presented so far are free of kinematical assumptions for cross-section behavior, here (and only here) we employ the kinematical assumption of Mindlin–Reissner, known also under the name of the first-order shear deformation (FOSD) theory [29]. Under the FOSD hypothesis, the generalized strains in the 3D shell-like body can be defined as

ε_{m} = ε_{m}^{M} + ζ ε_{b}^{M} .

(14)

For the drilling strains we postulate the following relation

ε_{d} = ε_{d}^{M},

(15)

which means that the microcurvatures $κ_{1}$ , $κ_{2}$ are constant through the thickness of the shell and equal to the corresponding values from base surface.

Through-the-thickness integration yields simplified expressions for the stress resultants and stress couples defined in equation (9) as functions of stresses defined in equation (11)

\begin{matrix} s_{m} = \int_{- h^{-}}^{+ h^{+}} σ_{m} d ζ = \int_{- h^{-}}^{+ h^{+}} [C_{mm} \underset{ε_{m}}{\underset{︸}{(ε_{m}^{M} + ζ ε_{b}^{M})}} + C_{md} \underset{ε_{d}}{\underset{︸}{ε_{d}^{M}}}] d ζ \\ = A_{mm} ε_{m}^{M} + B_{mb} ε_{b}^{M} + C_{md} ε_{d}^{M} \end{matrix},

(16)

A_{mm} = \int_{- h^{-}}^{+ h^{+}} C_{mm} d ζ, B_{mb} = B_{bm}^{T} = \int_{- h^{-}}^{+ h^{+}} C_{mm} ζ d ζ, C_{md} = C_{dm}^{T} = \int_{- h^{-}}^{+ h^{+}} C_{md} d ζ,

(17)

\begin{matrix} s_{b} = \int_{- h^{-}}^{+ h^{+}} σ_{m} ζ d ζ = \int_{- h^{-}}^{+ h^{+}} [C_{mm} (ζ ε_{m}^{M} + ζ^{2} ε_{b}^{M}) + ζ C_{md} ε_{d}^{M}] d ζ \\ = B_{bm} ε_{m}^{M} + E_{bb} ε_{b}^{M} + F_{bd} ε_{d}^{M} \end{matrix},

(18)

E_{bb} = \int_{- h^{-}}^{+ h^{+}} C_{mm} ζ^{2} d ζ, F_{bd} = F_{db}^{T} = \int_{- h^{-}}^{+ h^{+}} C_{mb} ζ d ζ,

(19)

\begin{matrix} s_{d} = \int_{- h^{-}}^{+ h^{+}} σ_{d} d ζ = \int_{- h^{-}}^{+ h^{+}} [C_{dm} (ε_{m}^{M} + ζ ε_{b}^{M}) + C_{dd} ε_{d}^{M}] d ζ \\ = C_{dm} ε_{m}^{M} + F_{db} ε_{b}^{M} + H_{dd} ε_{d}^{M} \end{matrix},

(20)

H_{dd} = \int_{- h^{-}}^{+ h^{+}} C_{dd} d ζ .

(21)

The constitutive operator for transverse shear stress resultants is given in the explicit form

s_{s} = D_{ss} ε_{s} = [\begin{matrix} G α_{s} & 0 \\ 0 & G α_{s} \end{matrix}] ε_{s},

(22)

where $α_{s}$ is the well-known shear correction factor [30]. Transverse shear effects are treated as elastic, what is justifiable for thin shells with $J_{2}$ plasticity [40].

Summarizing the ultimate structure of elastic constitutive relations for the resultant quantities becomes, in matrix form,

{\begin{matrix} s_{m} \\ s_{s} \\ s_{b} \\ s_{d} \end{matrix}} = [\begin{matrix} A_{mm} & 0 & B_{mb} & C_{md} \\ 0 & D_{ss} & 0 & 0 \\ B_{bm} & 0 & E_{bb} & F_{bd} \\ C_{dm} & 0 & F_{db} & H_{dd} \end{matrix}] {\begin{matrix} ε_{m} \\ ε_{s} \\ ε_{b} \\ ε_{d} \end{matrix}} \Leftrightarrow S = C^{e} ε .

(23)

The above equations are simplified in a sense that they do not contain the local curvature of the shell, i.e. they do not contain the shifter tensor and its determinant. This obvious simplification is valid only for thin shells, i.e. $h << L$ , where $L$ is a typical characteristic dimension of the shell structure, with a small curvature $h << R_{\min}$ . Exact equations can be found for instance in Lewiński and Telega [31].

It is interesting to investigate the structure of equation (21) in view of the equations used so far in Chróścielewski et al. [21, 22, 25, 30] and Witkowski [27] (compare also with discussions in Chróścielewski and Witkowski [17, 18]). The elastic constitutive relations for drilling stiffness are assumed as:

M^{1} = α_{t} \frac{1}{6} G h^{3} κ_{1}^{M}, M^{2} = α_{t} \frac{1}{6} G h^{3} κ_{2}^{M} \Rightarrow H_{dd} = [\begin{matrix} α_{t} \frac{1}{6} G h^{3} & 0 \\ 0 & α_{t} \frac{1}{6} G h^{3} \end{matrix}],

(24)

where $α_{t}$ is clearly a material coefficient. In numerical studies on the influence of $α_{t}$ taken from the range [10⁺¹⁰, 10⁻¹⁰] it has been shown [21, 25] that this value should be lower than 1 ( $α_{t} < 1$ ). On the other hand, the explicit form of equations (21) and (12) becomes

H_{dd} = \int_{- h^{-}}^{+ h^{+}} 2 {Gl}^{2} d ζ = [\begin{matrix} 2 l^{2} G h & 0 \\ 0 & 2 l^{2} G h \end{matrix}] .

(25)

Hence from the equivalence between equations (24) and (25) yields the value of $α_{t}$ in an elastic state

α_{t} \frac{1}{6} G h^{3} \equiv 2 l^{2} G h \Rightarrow α_{t} = 12 {(\frac{l}{h})}^{2},

(26)

which clearly shows that $α_{t}$ is indeed a material coefficient of the same character as, for example, the micropolar characteristic length $l$ in equation (11) or the shear correction factor $α_{s}$ in equation (22) and must not be interpreted as a penalty parameter. This provides a strong argument to defend our theory against some misinterpretations, for example, Eberlein and Wriggers [32] and Tan and Vu-Quoc [33], where $α_{t}$ was wrongly understood as a penalty coefficient.

The next step is to assume additive decomposition of the rate-of-strain tensor into elastic and plastic parts,

\overset{\cdot}{ε} = {\overset{\cdot}{ε}}^{e} + {\overset{\cdot}{ε}}^{p}, {\overset{\cdot}{ε}}^{e} = \overset{\cdot}{ε} - {\overset{\cdot}{ε}}^{p},

(27)

which is valid only for small strain plasticity [34 –37].

We further assume $J_{2}$ plasticity (in the classical context also known as the Maxwell–Huber–Hencky–von Mises theory [38]). The second invariant $J_{2}$ , following De Borst [28], is given by

J_{2} = a_{1} s_{ij} s_{ij} + a_{2} s_{ij} s_{ji} + a_{3} m_{ij} m_{ij} / l^{2}, s_{ij} = σ_{ij} - \frac{1}{3} δ_{ij} σ_{kk},

(28)

where $m_{ij}$ are components of the 3D Cosserat couple stress tensor.

The internal hardening variable is introduced as the generalized effective strain

{\overset{\cdot}{\bar{ε}}}^{p} = {[b_{1} {\overset{\cdot}{e}}_{ij}^{p} {\overset{\cdot}{e}}_{ij}^{p} + b_{2} {\overset{\cdot}{e}}_{ij}^{p} {\overset{\cdot}{e}}_{ji}^{p} + b_{3} {\overset{\cdot}{κ}}_{ij}^{p} {\overset{\cdot}{κ}}_{ij}^{p} l^{2}]}^{1 / 2},

(29)

with ${\overset{\cdot}{e}}_{ij} = {\overset{\cdot}{ε}}_{ij} - \frac{1}{3} δ_{ij} {\overset{\cdot}{ε}}_{kk}$ as the rate of deviatoric strain and ${\overset{\cdot}{κ}}_{ij}^{p}$ the rate of 3D Cosserat microcurvatures.

Equations (28) and (29) require specification of material constants $a_{i}, b_{i}$ , $i = 1, 2, 3$ . Various choices are possible [28]. Letting $m_{ij} = 0$ in equation (28) yields the expression for the second invariant in the form $J_{2} = (a_{1} + a_{2}) s_{ij} s_{ij}$ . It is easy to observe that to recover a classical definition of the second invariant the constants $a_{1}$ and $a_{2}$ should satisfy $a_{1} + a_{2} = \frac{1}{2}$ . As for equation (29) the following constraint should be observed $b_{1} + b_{2} = \frac{2}{3}$ . In this paper we assume

a_{1} = \frac{1}{4}, a_{2} = \frac{1}{4}, a_{3} = \frac{1}{2}, b_{1} = \frac{1}{3}, b_{2} = \frac{1}{3}, b_{3} = \frac{2}{3} .

(30)

Now the yield function may be written as

f (σ, A) = \sqrt{3 J_{2}} - σ_{Y} ({\bar{ε}}^{p}),

(31)

where $σ_{Y}$ is the yield stress, while the flow rule and the hardening evolution equation are given by the classical plasticity formulae

{\overset{\cdot}{ε}}^{p} = \overset{\cdot}{γ} \frac{\partial f}{\partial σ}, \overset{\cdot}{α} = \overset{\cdot}{γ} \frac{\partial f}{\partial A} .

(32)

In equations (31) and (32)₂, $α$ denotes the set (tensor) of the hardening variables (possibly linear) and $A$ the corresponding set (tensor) of hardening thermodynamic forces.

4. Numerical integration of constitutive relations

4.1 General structure of the closest point projection

We use the closest point projection algorithm [34 –37]. The rate equations (32) are replaced with corresponding incremental values within some (pseudo)time span $[t_{n}, t_{n + 1}]$ . Backward-Euler approximations read

ε_{n + 1}^{p} - ε_{n}^{p} = Δ γ {\frac{\partial f}{\partial σ} |}_{σ = σ_{n + 1}}, α_{n + 1} - α_{n} = Δ γ {\frac{\partial f}{\partial A} |}_{A = A_{n + 1}} .

(33)

The goal is now to find $ε_{n + 1}^{e}$ , $α_{n + 1}$ , $σ_{n + 1}$ and $Δ γ$ having the values of elastic strain $ε_{n}^{e}$ and $α_{n}$ from the beginning of the step and the values of the total strain $Δ ε = ε_{n + 1} - ε_{n}$ , so that the following constraints (discrete Kuhn–Tucker) are satisfied:

f_{n + 1} = f (σ_{n + 1}, A_{n + 1}) \leq 0, Δ γ \geq 0, f_{n + 1} Δ γ = 0 .

(34)

Define the elastic trial terms

ε_{n + 1}^{p trial} = ε_{n}^{p}, ε_{n + 1}^{e trial} = ε_{n + 1} - ε_{n + 1}^{p trial}, α_{n + 1}^{trial} = α_{n},

(35)

and residuals (derived from equations (33)_1,2 and (34)₁):

r_{ε} \equiv - ε_{n + 1}^{e} + ε_{n + 1}^{e trial} - Δ γ \frac{\partial f (σ_{n + 1}, A_{n + 1})}{\partial σ} = 0,

(36)

r_{α} \equiv - α_{n + 1} + α_{n + 1}^{trial} - Δ γ \frac{\partial f (σ_{n + 1}, A_{n + 1})}{\partial A} = 0,

(37)

r_{f} \equiv f (σ_{n + 1}, A_{n + 1}) = 0 .

(38)

In the present case of linear elastic material ( $C^{e} = const$ ), equation (36) may be recast in the stress space as

r_{σ} \equiv - σ_{n + 1} + σ_{n + 1}^{trial} - Δ γ C^{e} \frac{\partial f (σ_{n + 1}, A_{n + 1})}{\partial σ} = 0 .

(39)

The closest point projection algorithm [34, 35, 37] hinges on the idea that for the associative flow rule the stress space is endowed with an energy norm

| | σ | |_{C^{e}} = {〈 σ, {(C^{e})}^{- 1} σ 〉}^{\frac{1}{2}},

(40)

that induces the structure of a metric space

d (σ_{a}, σ_{b}) = | | σ_{a} - σ_{b} | |_{C^{e}} .

(41)

Now we look for such a $σ_{n + 1}$ satisfying equation (34)₁ that minimizes $d (σ_{n + 1}, σ_{n + 1}^{trial})$ . The system of equations (37)–(39) is solved using the Newton method. Omitting the pseudotime subscript ${(.)}_{n + 1}$ , after linearization we get

r_{f}^{(k + 1)} \approx r_{f}^{(k)} + {(\frac{\partial f}{\partial σ})}^{T (k)} Δ σ^{(k)} + {(\frac{\partial f}{\partial f α})}^{T (k)} Δ α^{(k)} + {(\frac{\partial f}{\partial Δ γ})}^{(k)} Δ {(Δ γ)}^{(k)} = 0,

(42)

r_{σ}^{(k + 1)} \approx r_{σ}^{(k)} + {(\frac{\partial r_{σ}}{\partial σ})}^{T (k)} Δ σ^{(k)} + {(\frac{\partial r_{σ}}{\partial α})}^{T (k)} Δ α^{(k)} + {(\frac{\partial r_{σ}}{\partial Δ γ})}^{(k)} Δ {(Δ γ)}^{(k)} = 0,

(43)

r_{α}^{(k + 1)} \approx r_{α}^{(k)} + {(\frac{\partial r_{f α}}{\partial σ})}^{T (k)} Δ σ^{(k)} + {(\frac{\partial r_{α}}{\partial α})}^{T (k)} Δ α^{(k)} + {(\frac{\partial r_{α}}{\partial Δ γ})}^{(k)} Δ {(Δ γ)}^{(k)} = 0,

(44)

with $k = 0, 1, 2 . . .$ as the iteration counter. The corrections are found through additive accumulation

{(.)}^{k + 1} = {(.)}^{k} + Δ {(.)}^{k} .

(45)

To start the process the following values are used:

σ^{(0)} = σ_{n + 1}^{trial}, α^{(0)} = α_{n}, Δ γ^{(0)} = 0 .

(46)

Iterations are terminated when $r_{f}^{(k)} \leq TOL 1$ , $| | r_{σ}^{(k)} | | \leq TOL 2$ and $| | r_{α}^{(k)} | | \leq TOL 2$ .

4.2 Specification for Cosserat plane stress with linear isotropic hardening

In this paper we work in the plane stress space. From equation (28) with constants given in equation (30) we have

J_{2} = \frac{1}{3} (σ_{11}^{2} + σ_{22}^{2} - σ_{11} σ_{22} + \frac{3}{4} (σ_{12}^{2} + σ_{21}^{2}) + \frac{3}{2} σ_{12} σ_{21} + \frac{3}{2 l^{2}} (m_{1}^{2} + m_{2}^{2}))

(47)

and from equation (29) we obtain

{\overset{\cdot}{\bar{ε}}}^{p} = [\frac{2}{3} ({({\overset{\cdot}{e}}_{11}^{p})}^{2} + {({\overset{\cdot}{e}}_{22}^{p})}^{2}) + \frac{1}{3} {({\overset{\cdot}{e}}_{12}^{p})}^{2} + \frac{2}{3} {\overset{\cdot}{e}}_{12} {\overset{\cdot}{e}}_{21}^{p} + \frac{1}{3} {({\overset{\cdot}{e}}_{21}^{p})}^{2} {+ \frac{2}{3} ({({\overset{\cdot}{κ}}_{1}^{p})}^{2} + {({\overset{\cdot}{κ}}_{2}^{p})}^{2})]}^{1 / 2} .

(48)

Following De Borst [28] we introduce the matrix operator

P = [\begin{matrix} \frac{2}{3} & - \frac{1}{3} & 0 & 0 & 0 & 0 \\ - \frac{1}{3} & \frac{2}{3} & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{2} & \frac{1}{2} & 0 & 0 \\ 0 & 0 & \frac{1}{2} & \frac{1}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}],

(49)

which allows for writing equation (31) as

f = \sqrt{\frac{3}{2} σ^{T} P σ} - σ_{Y}, (σ_{H - M - H} = \sqrt{\frac{3}{2} σ^{T} P σ}) .

(50)

In the above notation the components of equations (42) and (43) become

\frac{\partial f}{\partial σ} = \sqrt{\frac{3}{4}} \frac{{P σ}^{T} P σ - P σ σ^{T} P}{\sqrt{J_{2}}} .

(51)

We assume linear isotropic hardening, so that

σ_{Y} ({\bar{ε}}^{p}) = σ_{Y}^{0} + H {\bar{ε}}^{p}, \overset{\cdot}{γ} = {\overset{\cdot}{\bar{ε}}}^{p}, \frac{\partial f}{\partial Δ γ} = H,

(52)

where $H$ is the constant hardening modulus and $σ_{Y}^{0}$ is the initial yield stress. Relation (52) is justifiable when the set of constants in equation (30) is introduced. Furthermore, we have

\frac{\partial r}{\partial σ} = I - Δ γ C^{e} \frac{\partial^{2} f}{\partial σ^{2}} = I - Δ γ C^{e} (\sqrt{\frac{3}{4 J_{2}}} P - \frac{\sqrt{3}}{4} {(J_{2})}^{- 3 / 2} P σ σ^{T} P),

(53)

\frac{\partial r}{\partial Δ γ} = - C^{e} \frac{\partial f}{\partial σ} = - \sqrt{\frac{3}{4}} \frac{C^{e} P σ}{\sqrt{J_{2}}},

(54)

Δ γ = \frac{{(\frac{\partial f}{\partial σ})}^{T} H \overset{\cdot}{ε}}{{(\frac{\partial f}{\partial σ})}^{T} H \frac{\partial f}{\partial σ} - \frac{\partial f}{\partial Δ γ}}, H = {[1 + Δ γ C^{e} \frac{\partial f}{\partial σ}]}^{- 1} C^{e} .

(55)

In numerical simulations we use the consistent elastoplastic modulus given by De Borst [28]

\overset{\cdot}{σ} = [H - \frac{{(\frac{\partial f}{\partial σ})}^{T} HH \frac{\partial f}{\partial σ}}{{(\frac{\partial f}{\partial σ})}^{T} H \frac{\partial f}{\partial σ} - \frac{\partial f}{\partial Δ γ}}] \overset{\cdot}{ε} = C^{ep} \overset{\cdot}{ε} .

(56)

5. Numerical examples

Two numerical examples presented in this section show the range of applicability of the presented formulation. Numerical analysis of other regular shells, within the five-parameter shell theory framework and elastoplastic material, can be found in Brank [46] and Li et al. [47].

The first example is concerned with a regular shell in the form of a cylinder, while the second one is concerned with a branched shell structure. In analyses we use the 16-node shell finite element, denoted CAMe16FI, with fully integrated matrices. This shell element (see for instance Chróścielewski et al. [21, 25]) is a displacement-rotation type and uses Lagrangian shape functions. For rotations, a specially designed interpolation scheme is developed as described for instance in Chróścielewski et al. [19, 21]. The numerical code is written in Fortran with a parallel solver [39].

5.1 Cylinder under shear load, regular shell

The first example, taken from Dujc and Brank [40], is concerned with a regular cylindrical shell clamped at both ends, subjected to a shear-like load due to the imposed translation on the upper edge, see Figure 1.

Figure 1.

Cylinder under shear load, geometry.

The radius of the cylinder is $R = 28.5$ , its height is $H = 85$ and the thickness of the shell is $h_{0} = 0.5$ . The material is assumed to be elastic perfectly plastic with Young’s modulus $E = 2.1 \times 10^{4}$ , Poisson’s ratio $v = 0.3$ , and yield stress $σ_{Y} = 24$ . Micropolar constants are assumed as $G_{c} = 0.8077 \times 10^{4}$ and $l = 10^{- 4}$ , and load $P_{ref} = 1$ $P = λ P_{ref}$ (for load we mean the sum of reactions in the z direction on the upper edge). We utilize the symmetry of the structure using 24 × 24 CAMe16FI elements for half of the shell. Through-the-thickness integration was performed using a seven-node Gauss quadrature.

The obtained results are presented in Figure 2 and are compared to the reference solution from Dujc and Brank [40], where 64 × 36 elements were used for the whole shell, and to ABAQUS results (element S8R, own calculations). All the results show excellent agreement even though different plasticity formulations and different finite elements were employed.

Figure 2.

Cylinder under shear load, results.

Next we concentrate on the loci of the limit point as a function of mesh density. In the course of non-linear convergence analyses with CAMe16 elements and ABAQUS we have obtained the results shown in Figure 3. The results converge to the load multiplier close to $λ \approx 3.45$ . To give better insight into the results we also present them in Table 1.

Figure 3.

Cylinder under shear load, convergence analysis of limit point.

Table 1.

Convergence analysis of loci of the limit point for a cylinder under shear load.

Mesh, formulation	Load multiplier λ
CAMe16FI 12 × 12	3.47544
CAMe16 FI 24 × 24	3.45442
CAMe16 FI 36 × 36	3.45115
S8R 18 × 18 ABAQUS	3.46684
S8R 36 × 36 ABAQUS	3.45065
S8R 54 × 54 ABAQUS	3.44777

In the last step we investigate the influence of micropolar material parameters. Our previous study [24] revealed that $G_{c}$ taken from $G_{c} = {10; 1; 0.1} \times G$ plays a negligible role. Therefore, we investigated only the influence of micropolar length $l$ on the plastic range, which is independent of plastic flow. For this purpose we use a CAMe16FI 24 × 24 mesh. The results are shown in Figure 4. As can be observed, with the growth of the ratio of micropolar length to the shell thickness $l / h$ the analyzed shell exhibits stiffening behavior, which is clearly marked by the equilibrium path for $l / h = 100$ . As the ratio $l / h$ becomes smaller, the responses become indistinguishable, which is shown in Table 2. Hence a conclusion can be drawn that for $l / h \leq 1$ differences become small and the quantitative character of the solution is close to other formulations.

Figure 4.

Cylinder under shear load, influence of the micropolar length l.

Table 2.

Influence of the $l / h$ ratio on loci of the limit point for a cylinder under shear load.

$l / h$	Load multiplier λ
100	8.28612
10	3.77147
1	3.47142
0.1	3.45487
0.01	3.45442

It is interesting to review the influence of $l / h$ on the overall deformation of the shell. To this end we examined two cases $l / h = 0.01$ and $l / h = 100$ under the same value of deflection $w = 10$ . The results are presented in Figure 5. In the first case, the shell deforms with a visible action of the tension field type, while in the second case, the overall shape of deformation is practically smooth. Analysis of the equivalent plastic strain (EPL, equation (29)) contours shown in Figure 6. reveals that within the area of deformation wave values of EPL tend to localize.

Figure 5.

Cylinder under shear load, influence of micropolar length l on deformation: left $l / h = 0.01$ , right $l / h = 100 .$

Figure 6.

Cylinder under shear load, influence of micropolar length l on distribution of effective plastic strain: left $l / h = 0.01$ , right $l / h = 100 .$

5.2 Channel-section clamped beam, irregular shell

The class of shell problems containing branches, irregularities and orthogonal intersections always poses a challenge for numerical analysis. This is due to the problem of the drilling degree of freedom (6th DOF). We wish to emphasize that this issue is absent in the present formulation, since the drilling DOF emerges naturally as the natural consequence of the six-parameter shell theory. There exists some demanding examples in the literature both for dynamics, for example, Chróścielewski and Witkowski [41] and Campello et al. [42] and for statics, presented here.

In this numerical example we consider a channel-section beam loaded with vertical force, fully clamped on the other end. It was presented originally in Chróścielewski et al. [25] with purely elastic material and studied later [32, 33, 43 –45]. This example amplifies the importance of the drilling DOF and associated drilling stiffness, and their influence on deformation of the order of the height of the structure. In this range, the solutions have exclusively shell character and cannot be treated on the grounds of beam theory. Hence, the proportion of the dimensions used here is important for analysis of correctness of various shell formulations.

The analyzed structure is shown in Figure 7. The total length of the cantilever is $L = 36$ , height $H = 6$ , width of flanges $B = 2$ , and uniform shell thickness $h = 0.05$ . The material is characterized by a Young’s modulus $E = 10^{7}$ , Poisson’s ratio $v = 0.3333$ , load $P_{ref} = 100$ and $P = λ P_{ref}$ , and yield stress of an elastic perfectly plastic material $σ_{Y} = 5 \times 10^{3}$ . The micropolar parameters assumed are reported for each result when necessary. For computation, 16-node CAM elements with discretization (4 + 8 + 4) × 36, i.e. the number of elements per: (flange + web + flange) × length and full integration in the surface was chosen. The shell resultants were integrated with a seven-point Gauss quadrature. In order to avoid local plastic effects under point load, the row of finite elements under the load is assumed to be purely elastic.

Figure 7.

Channel section clamped beam, geometry, shown is the 16-node element used for discretization.

Figure 8 shows the comparison of our results with reference solutions in the elastoplastic range, where our results are obtained for $G_{c} = G$ and $l / h = 10^{- 2}$ . Similar to the previous Figure 9 with the data given in Table 3. In general, the results converge to the load multiplier close to $λ \approx 0.77$ . The influence of $l / h$ on the overall behavior of the beam is shown in Figure 10 (data in Table 4). As in the previous example for a smooth shell, for $l / h \leq 1$ here the differences also become small, and the quantitative character of the solution is close to other formulations.

Figure 8.

Channel section clamped beam, results.

Figure 9.

Channel section clamped beam, convergence analysis for limit point.

Table 3.

Convergence analysis of loci of the limit point for a channel section.

Mesh, formulation	Load multiplier λ
CAMe16FI 8 × 18	0.773134
CAMe16 FI 16 × 36	0.770355
CAMe16 FI 24 × 54	0.769376
S8R 12 × 27 ABAQUS	0.773879
S8R 24 × 54 ABAQUS	0.771039
S8R 36 × 81 ABAQUS	0.769901

Figure 10.

Channel section clamped beam, influence of micropolar length l.

Table 4.

Influence of the $l / h$ ratio on loci of the limit point for a channel section.

$l / h$	Load multiplier λ
100	2.762495
10	0.931142
1	0.772526
0.1	0.770377
0.01	0.770355

Figure 11 shows the map of the equivalent plastic stress H–M–H (equation (50)₂, generalized Huber–von Mises–Hencky) on the undeformed configuration for $u_{(a)} = 3$ , $l / h = 10^{- 2}$ .

Figure 11.

Channel section clamped beam: left deformed configuration for $u_{(a)} = 3$ and right corresponding contour map of equivalent plastic stress H–M–H (Huber–von Mises–Hencky).

6. Conclusions

We have presented the complete formulation of elastoplastic constitutive relations for the non-linear six-parameter shell theory with drilling rotation. The relations are obtained by through-the-thickness integration of the asymmetric Cosserat plane stress relation, which naturally correspond to the underlying kinematics of the present shell theory.

By comparing our results to alternative formulations from the literature, we have shown the correctness and effectiveness of our formulation. In addition, since the obtained formulation naturally includes the characteristic length of material, we have evaluated its influence on the overall deformation and response of the shells. It has been shown that the values of the ratio of characteristic length to shell thickness greater than 10 are responsible for the stiffening behavior.

The analysis performed within this study shows that in our shell theory the growth of the value of characteristic length has a great influence on plastic response. The plastic zone becomes ‘smeared’ whereby localization effects are minimized. We have shown that $l / h > 1$ leads to large discrepancies between the present formulation and reference solutions. This fact allows for determining the allowable range of changes of $l / h \leq ~ 1$ . Further work should be focused on accomplishing regularization of results (independence from FEM discretization) in analysis with strain softening materials.

Footnotes

Acknowledgements

ABAQUS calculations were carried out at the Academic Computer Center in Gdańsk. The parallel solver for CAM elements was developed on the basis of HSL, a collection of Fortran codes for large-scale scientific computation ().

Funding

Stanisław Burzyński is supported by the National Science Centre of Poland (grant number DEC–2012/05/D/ST8/02298).

References

Capriz

Continua with microstructure. New York: Springer, 1989.

Eringen

. Microcontinuum field theories I: Foundations and solids. New York: Springer, 1999.

Cosserat

. Théorie des corps déformables. Paris: Librairie Scientifique, 1909.

Nowacki

. Theory of asymmetric elasticity. Oxford: Pergamon Press, 1986.

Rubin

. Cosserat theories: Shells, rods and points. Dordrecht: Kluwer Academic Publishing, 2000.

Altenbach

Eremeyev

. On generalized Cosserat-type theories of plates and shells: a short review and bibliography. Arch Appl Mech 2010; 80: 73−92.

Eremeyev

Lebedev

. Existence theorems in the linear theory of micropolar shells. Z Angew Math Mech 2011; 91: 468–476.

Altenbach

Eremeyev

. On the linear theory of micropolar plates. Z Angew Math Mech 2009; 89: 242–256.

Bîrsan

Neff

. Existence theorems in the geometrically non-linear 6-parameter theory of elastic plates. J Elast 2013; 112: 185–198.

10.

Pietraszkiewicz

Eremeyev

. On natural strain measures of the non-linear micropolar continuum. Int J Solids Struct 2009; 46: 774–787.

11.

Whitman

DeSilva

. A dynamical theory of elastic directed curves. Z Angew Math Phys 1969; 20: 200–212.

12.

Cardona

Geradin

. A beam finite element non-linear theory with finite rotations. Int J Numer Methods Eng 1988; 26: 2403–2438.

13.

Sander

. Geodesic finite elements for Cosserat rods. Int J Numer Methods Eng 2009: 82; 1645–1670.

14.

Nowacki

. Couple-stresses in the theory of thermoelasticity. In: Irreversible aspects of continuum mechanics and transfer of physical characteristics in moving fluids (ed H

Parkus

Sedov

), Vienna, Austria, 22–28 June 1966, pp. 259–278. Wien: Springer-Verlag.

15.

Cowin

. Stress functions for Cosserat elasticity. Int J Solids Struct 1970; 6: 389–398.

16.

Eremeyev

Pietraszkiewicz

. Local symmetry group in the general theory of elastic shells. J Elast 2006; 85: 125–152.

17.

Chróścielewski

Witkowski

. On some constitutive equations for micropolar plates. Z Angew Math Mech 2010; 90: 53–64.

18.

Chróścielewski

Witkowski

. FEM analysis of Cosserat plates and shells based on some constitutive relations. Z Angew Math Mech 2011; 91: 400–412.

19.

Chróścielewski

Kreja

Sabik

Witkowski

. Modeling of composite shells in 6-parameter nonlinear theory with drilling degree of freedom. Mech Adv Mater Struct 2011; 18: 403–419.

20.

Libai

Simmonds

. The nonlinear theory of elastic shells. Cambridge: Cambridge University Press, 1998.

21.

Chróścielewski

Makowski

Pietraszkiewicz

. Statics and dynamics of multifold shells: Nonlinear theory and finite element method (in Polish). Warsaw: Wydawnictwo IPPT PAN, 2004.

22.

Chróścielewski

Witkowski

. Four-node semi-EAS element in six-field nonlinear theory of shells. Int J Numer Methods Eng 2006; 68: 1137–1179.

23.

Bîrsan

Neff

. Existence of minimizers in the geometrically non-linear 6-parameter resultant shell theory with drilling rotations. Math Mech Solids 2014; 19: 376–397.

24.

Burzyński

Chróścielewski

Witkowski

. Elastoplastic material law in 6-parameter nonlinear shell theory. In: Pietraszkiewicz

Górski

(eds.) Shell structures: Theory and applications. London: CRC Press, 2014, vol. 3, pp. 377−380.

25.

Chróścielewski

Makowski

Stumpf

. Genuinely resultant shell finite elements accounting for geometric and material non-linearity. Int J Numer Methods Eng 1992; 35: 63–94.

26.

Konopińska

Pietraszkiewicz

. Exact resultant equilibrium conditions in the non-linear theory of branching and self-intersecting shells. Int J Solids Struct 2007; 44: 352–369.

27.

Witkowski

. 4-Node combined shell element with semi-EAS-ANS strain interpolations in 6-parameter shell theories with drilling degrees of freedom. Comp Mech 2009; 43: 307–319.

28.

De Borst

. Simulation of strain localization a reappraisal of the Cosserat continuum. Eng Comput 1991; 8: 317−332.

29.

Kreja

. A literature review on computational models for laminated composite and sandwich panels. Cent Eur J Eng 2011; 1: 59−80.

30.

Chróścielewski

Pietraszkiewicz

Witkowski

. On shear correction factors in the non-linear theory of elastic shells. Int J Solids Struct 2010; 47: 3537–3545.

31.

Lewiński

Telega

. Plates, laminates and shells, asymptotic analysis and homogenization. Singapore: World Scientific, 2000.

32.

Eberlein

Wriggers

. Finite element concepts for finite elastoplastic strains and isotropic stress response in shells: theoretical and computational analysis. Comp Methods Appl Mech Eng 1999; 171: 243−279.

33.

Tan

Vu-Quoc

. Efficient and accurate multilayer solid-shell element: non-linear materials at finite strain. Int J Numer Methods Eng 2005; 63: 2124–2170.

34.

de Souza Neto

Peric

Owen

DRJ

. Computational methods for plasticity: Theory and applications. Chichester: Wiley, 2008.

35.

Simo

Hughes

TJR

. Computational inelasticity. New York: Springer, 1998.

36.

Armero

. Elastoplastic and viscoplastic deformation is solids and structures In: Stein

De Borst

Hughes

TJR

(eds.) Encyclopedia of computational mechanics. Chichester: Wiley, 2004.

37.

Belytschko

Liu

Moran

Elkhodary

. Nonlinear finite elements for continua and structures. Chichester: Wiley, 2013.

38.

Zhuo

Vu-Quoc

Izzuddin

Wei

. A four-node corotational quadrilateral elastoplastic shell element using vectorial rotational variables. Int J Numer Methods Eng 2013; 95: 181−211.

39.

HSL. A collection of Fortran codes for large-scale scientific computation. http://www.hsl.rl.ac.uk/

40.

Dujc

Brank

. Stress resultant plasticity for shells revisited. Comp Methods Appl Mech Eng 2012; 247–248: 146–165.

41.

Chróścielewski

Witkowski

. Discrepancies of energy values in dynamics of three intersecting plates. Int J Numer Methods Biomed Eng 2010; 26: 1188−1202.

42.

Campello

EMB

Pimenta

Wriggers

. An exact conserving algorithm for nonlinear dynamics with rotational DOFs and general hyperelasticity. Part 2: Shells. Comp Mech 2011; 48: 95−211.

43.

Merlini

Morandini

. Computational shell mechanics by helicoidal modeling II: Shell element. J Mech Mater Struct 2011; 6: 693–728.

44.

Wiśniewski

. Finite rotation shells: Basic equations and finite elements for Reissner kinematics. Berlin: Springer, 2010.

45.

Klinkel

Gruttmann

Wagner

. A mixed shell formulation accounting for thickness strains and finite strain 3D material models. Int J Numer Methods Eng 2008; 74: 945–970.

46.

Brank

Peric

Damjanic

. On large deformations of thin elasto-plastic shells: implementation of a finite rotation model for quadrilateral shell element. Int J Numer Methods Eng 1997; 40: 689–726.

47.

Zhuo

Vu-Quoc

Izzuddin

Wei

. A four-node corotational quadrilateral elastoplastic shell element using vectorial rotational variables. Int J Numer Methods Eng 2013; 95: 181–211.