A Co-Rotational Based Anisotropic Elasto–Plastic Model for Geometrically Non-Linear Analysis of Fibre Reinforced Polymer Composites: Formulation and Finite Element Implementation

Abstract

Geometrical non-linearity is one of the aspects to be taken into account for accurate analysis of fibre reinforced polymers (FRPs), since large displacements and rotations may be observed in many of its structural applications such as in aircraft wings and wind turbine blades. In this paper, a co-rotational formulation and implementation of an invariant-based anisotropic plasticity model are presented for geometrically non-linear analysis of FRPs. The anisotropic constitutive equations are formulated in the format of isotropic tensors functions. The model assumes an anisotropic pressure-dependent yield function, and in addition to this, a non-associated plastic potential function in order to model realistic plastic deformations in FRPs. The formulation is then cast in the co-rotational framework to consider the geometrical non-linear effects in an efficient manner. The developed model is implemented in the commercial finite element (FE) software ABAQUS/Implicit via the means of the user-defined material subroutine (UMAT). The kinematics within the co-rotational frame is explained briefly while the important aspects regarding the numerical treatment and implementation are discussed in detail. Representative numerical examples at different scales are presented to demonstrate the applicability and robustness of the proposed development.

1. Introduction

Modern industry demands materials which are environmentally friendly by reducing the carbon footprint, improving safety by offering higher strengths and resistance to fatigue etc., and decreasing operational costs through virtue of fewer inspections and repairs required [1]. Recent advances in composites materials, more specifically fibre reinforced polymers (FRPs), are helping to replace traditional materials across a host of engineering applications by offering a combination of high strength to weight ratio, high stiffness, better fatigue response, reduced environmental effects, and faster manufacturing among others [2,3].

With a continuously evolving trend of shifting to composite materials, there is an ever-present need for better understanding of material behavior. Starting from simple analytical approaches to explain the material behavior, the focus gradually shifted to a more realistic and complex three-dimensional representation in the past few decades. Hence, the need to use numerical modeling came to the fore. This pronounced complexity poses a stern challenge as FRPs pose temperature, pressure, size, and rate dependencies along with the more obvious anisotropic behavior [4,5,6], and progressive failure [7,8,9].

For an accurate mechanical response prediction of FRPs along with failure behavior by means of numerical modeling techniques, the representation of anisotropy (fibre orientation) plays a fundamental role. To account for the inherently present anisotropy in the material modeling of this class of composites, generally, two main strategies are used: (i) multi-scale approach which involves in principle modeling microscopic constituents separately at corresponding scale, and (ii) macroscopic phenomenological approach which takes advantage of using extra so-called internal variables (damage, plasticity among others) to represent the characteristic non-linear material behavior under distinct loading cases. Some detailed reviews on multi-scale modeling concerning composite materials and corresponding comparisons can be found in the referenced literature [10,11,12] among many others. As the motivating thought behind numerical modeling as a virtual testing solution is efficiency along with detailed understanding of material response, the main drawback from multi-scale analysis i.e., increased computational costs goes against the soul of the objectives [13]. As a consequence, the employment of the multi-scale technique in practical engineering problems can become rather limited and impractical.

Opposed to the multi-scale approach, the anisotropic macroscopic phenomenological material modeling approach accounting for fibre orientation is promising for large engineering problems having practical real-world implications. In addition to reduced computational costs because of modeling at a single scale, only a handful number of experiments is needed for calibration and subsequent validation purposes [14,15,16]. Incorporating anisotropy into macroscopic phenomenological models can be achieved in a number of ways. One such framework is based on invariant theory [17]. In the context of this approach, the response of the material is described using scalar-valued functions through several tensorial variables such as deformation and stress tensors. To account for anisotropy, the argument list in these functions definition is extended by the so-called structural tensors which reflect the inherent symmetry of the composites material. The resulting general form of the constitutive equations is automatically invariant under coordinate transformation. For further details and a comprehensive review on the topic, the following references are useful [18,19].

In many of its structural applications (such as wind turbine blades, aircraft wings etc.), FRPs undergo large deflections and rotations, but the strains are usually within the small to moderate range because of high in-plane stiffness. Considering such behavior of FRPs, it is advantageous to use the small strain constitutive modeling which is relatively easier to handle and computationally less expensive compared to the finite strain modeling strategy (see [20]), but additionally, account for large deflections and rotations [21]. The co-rotational Lagrangian formulation provides the solution, where the idea is to decompose the motion of the body into rigid body motions i.e., deflections and rotations, and pure deformations. It has been mostly employed for beam and shell formulations for isotropic materials [22,23,24,25] as such beam and shell elements are used for applications within small deformations, but it is not limited to that. Rather, it can be employed in any finite element (FE) formulation where the basic assumption of small strains and arbitrary rotations is fulfilled as highlighted in the references [26,27,28].

The fundamental concept in the co-rotational formulation is the split of motions of a continuum body into two steps. In the first step, the rigid translations and rotations of the undeformed body are considered. The rigid translations are defined by displacements expressed in the global frame of reference. The rigid rotations, defined by an orthogonal rotation matrix, defines the orientation of local frame in the deformed configuration. In the second step, the local deformation of the body with respect to the local frame of reference is considered. This approach predates finite element methods (FEM) by over a century. Recently, the idea has successfully found application in FEM [25,28]. The pure deformation part of the displacement field, obtained by subtracting rigid body motions from the total displacement field, tends to be small when the incremental motion is sufficiently small. This argumentation is the basis for the infinitesimal magnitude of strains in the rotated frame. In a spatially discretized domain such as in FEM, this decomposition of the motions of the body is achieved by defining a local co-rotational frame for each discretized element. This local frame does not deform but rather translates and rotates with the element. The pure deformational part of the motion measured with respect to this local frame is small. Hence, the discrete gradients of the deformational displacement field in the local frame are of the order of small strains [29]. This key concept helps to simplify the updated Lagrangian formulation to the co-rotational formulation. For more details, the reader is referred to [24,30].

In this contribution, an invariant-based anisotropic elasto–plastic model is formulated and implemented within the co-rotational framework for its application in geometrically non-linear analysis of FRPs. The anisotropic constitutive equations are represented in the form of isotropic tensor functions. Accordingly, an anisotropic pressure-dependent yield surface is introduced along with a non-associative plastic potential function to account for the non-linear inelastic material behavior [31]. Employing the non-associative flow rule allows for modeling realistic plastic deformation compared to associative plasticity, especially with regard to contractility/dilatancy effects resulting in different behavior under compression and tension as it is observed in composites. The model is then cast into the co-rotational framework so that the geometrical non-linear effects (large deflections and rotations) can be included. It is to be noted that although the strains are assumed to be within small to moderate range, they are not exactly the small strains obtained using linear deformation theory [24,30]. Afterward, the computational aspects corresponding to the algorithmic treatment of the proposed model and its numerical implementation are detailed. Novel closed form expressions, necessary for a consistent FE implementation, are also derived.

For the sake of transparency, this paper focuses on the extension of the geometrically linear plasticity model presented in [15] for unidirectional (UD) FRPs, to take into account the geometrical non-linear effects due to large displacements and rotations. In comparison to the constitutive model in [15], modifications to the yield and plastic potential function definitions are proposed. These modifications allow for an easier calibration of the yield surface and plastic potential function with the experimental data. In this regard, explicit expressions for the yield surface and plastic potential parameters are provided. From the computational side, herein an explicit expression for the algorithmic consistent tangent moduli is derived.

The paper is organized into the following sections: Section 2 discusses the constitutive formulation of the invariant-based anisotropic elasto–plastic material model within the co-rotational framework in detail. Section 3 details the numerical treatment of the proposed model including the FE implementation procedure in the commercial software ABAQUS 2017. Thereafter, some numerical results are presented in Section 4 to highlight the validity and range of application of the proposed formulation. Finally, the main conclusions of the current contribution are drawn in Section 5.

2. Constitutive Formulation

This section presents the constitutive formulation of the anisotropic invariant-based model for FRPs. It is to be noted that the constitutive model proposed here is a modification of the one presented in [15]. These modifications include a new form of the yield and plastic potential functions. Nevertheless, for the sake of clarity and completeness, the constitutive formulation is provided in detail.

It should be noted that the constitutive equations are formulated with respect to the co-rotational frame.

2.1. Transversely Isotropic Free-Energy Definition

From the modeling standpoint, the anisotropic mechanical response admits a tensor-based representation through the definition of a second order structural tensor $\mathbf{A}$ in the rotated frame. The structural tensor represents the anisotropic material inherent structure and is defined as:

$$\mathbf{A}:=\mathbf{a}\otimes \mathbf{a},$$

(1)

where $\mathbf{a}$ identifies the fibre orientation vector in the rotated frame.

Based on the hypothesis of the flow theory of plasticity, the total strain tensor $\mathit{\epsilon }$ is additively decomposed into elastic ${\mathit{\epsilon }}^e$ and plastic ${\mathit{\epsilon }}^p$ counterparts as follows:

$$\mathit{\epsilon }={\mathit{\epsilon }}^e+{\mathit{\epsilon }}^p.$$

(2)

For the constitute formulation, the existence of a Helmholtz free-energy function, $\mathsf{\Psi }\left({\mathit{\epsilon }}^e,\mathbf{A},\mathit{\varpi }\right)$ is assumed. This free-energy is a function of the elastic strain ${\mathit{\epsilon }}^e$, the structural tensor $\mathbf{A}$, and the internal variable set $\mathit{\varpi }$ that accounts for the inelastic material response along the deformation process:

$$\mathsf{\Psi }\left({\mathit{\epsilon }}^e,\mathbf{A},\mathit{\varpi }\right)=\frac{1}{2}{\mathit{\epsilon }}^e:{\mathbb{C}}^e:{\mathit{\epsilon }}^e+{\mathsf{\Psi }}^{\mathrm{hard}}\left(\mathbf{A},\mathit{\varpi }\right),$$

(3)

where ${\mathbb{C}}^e$ is the constitutive elastic tensor and ${\mathsf{\Psi }}^{\mathrm{hard}}\left(\mathbf{A},\mathit{\varpi }\right)$ is the hardening part of the free-energy function due to plastic effects.

Having on hand the free energy function definition, the constitutive stress tensor $\mathit{\sigma }$ is obtained as the first derivative of the free energy function with respect to the elastic strain tensor, while the elastic constitutive operator ${\mathbb{C}}^e$ is defined as the second derivative of the free energy with respect to elastic strain tensor:

$$\mathit{\sigma }:=\frac{\partial \mathsf{\Psi }}{\partial {\mathit{\epsilon }}^e}={\mathbb{C}}^e:{\mathit{\epsilon }}^e,$$

(4)

For transversely isotropic materials, the constitutive transversely isotropic elasticity tensor is represented as follows:

$${\mathbb{C}}^e:=\frac{{\partial }^2\mathsf{\Psi }}{\partial {\mathit{\epsilon }}^e\partial {\mathit{\epsilon }}^e}=\lambda \mathbf{1}\otimes \mathbf{1}+2{\mu }_T\mathbb{I}+\alpha (\mathbf{1}\otimes \mathbf{A}+\mathbf{A}\otimes \mathbf{1})+2({\mu }_L-{\mu }_T){\mathbb{I}}_A+\beta \mathbf{A}\otimes \mathbf{A},$$

(5)

where $\mathbb{I}$ refers to the fourth-order identity tensor, whereas ${\mathbb{I}}_A=A_{im}{\mathbb{I}}_{jmkl}+A_{jm}{\mathbb{I}}_{mikl}$, and $\lambda$, $\alpha$, $\beta$, ${\mu }_T$ and ${\mu }_L$ are the elastic constants. Their definition and relationship to the engineering constants are given in [4].

2.2. Thermodynamics Considerations

The constitutive equations are restricted by the second-law of thermodynamics in the form of the Clausius–Duhem inequality. Under the assumption of isothermal deformations, this inequality reads the following for internal energy dissipation ${\mathcal{D}}_{\mathrm{int}}$:

$${\mathcal{D}}_{\mathrm{int}}=\mathit{\sigma }:\dot{\mathit{\epsilon }}-\dot{\mathsf{\Psi }}\ge 0.$$

(6)

Recalling the previous definitions, the restriction over the internal dissipation reads:

$${\mathcal{D}}_{\mathrm{int}}=\mathit{\sigma }:{\dot{\mathit{\epsilon }}}^p+\mathsf{\Gamma }*\dot{\mathit{\varpi }}\ge 0.$$

(7)

where $\mathsf{\Gamma }$ denotes the so-called hardening force and * stands for any arbitrary product.

2.3. Yield Function

The elastic domain $\mathbb{E}$, assuming the maximum dissipation principle, is defined as:

$$\mathbb{E}=\left\{(\mathit{\varpi },{\overline{\epsilon }}^p)\right|\mathcal{F}(\mathit{\sigma },\mathbf{A},{\overline{\epsilon }}^p)\le 0\},$$

(8)

where ${\overline{\epsilon }}^p$ identifies the equivalent plastic strain. The definition of the equivalent plastic strain in the present formulation is given by:

$${\overline{\epsilon }}^p:=\sqrt{\frac{1}{2}}\parallel {\mathit{\epsilon }}^p\parallel .$$

(9)

The construction of a transversely isotropic yield surface $\mathcal{F}(\mathit{\sigma },\mathbf{A},{\overline{\epsilon }}^p)$, which accounts for the pressure-dependency and plastic-inextensibility in FRPs along the fibre direction yields:

$$\mathcal{F}(\mathit{\sigma },\mathbf{A},{\overline{\epsilon }}^p)={\zeta }_1{I}_1+{\zeta }_2{I}_2+{\zeta }_3{I}_3+{\zeta }_4{I}_3^2-1\le 0,$$

(10)

$I_i$ ($i=1,3$) are the stress invariants which symbolize the integrity basis of an isotropic tensor function representing a transversely isotropic response:

$$I_1=\frac{1}{2}{\left(\mathrm{tr}\left[{\mathit{\sigma }}^{\mathrm{pind}}\right]\right)}^2-\mathrm{tr}\left[\mathbf{A}{\left({\mathit{\sigma }}^{pind}\right)}^2\right];I_2=\mathrm{tr}\left[\mathbf{A}{\left({\mathit{\sigma }}^{\mathrm{pind}}\right)}^2\right];I_3=\mathrm{tr}\left[\mathit{\sigma }\right]-\mathrm{tr}\left[\mathbf{A}\left(\mathit{\sigma }\right)\right],$$

(11)

where ${\mathit{\sigma }}^{\mathrm{pind}}$ is the plasticity-inducing stress:

$${\mathit{\sigma }}^{\mathrm{pind}}:=\mathit{\sigma }-\frac{1}{2}(\mathrm{tr}\left[\mathit{\sigma }\right]-\mathbf{a}\mathit{\sigma }\mathbf{a})\mathbf{1}+\frac{1}{2}(\mathrm{tr}\left[\mathit{\sigma }\right]-3\mathbf{a}\mathit{\sigma }\mathbf{a})\mathbf{A},$$

(12)

Here, ${\zeta }_i\left({\overline{\epsilon }}^p\right)$ ($i=1,4$) refers to four yield parameters which together with their corresponding invariants represent different loading states.

A compact representation of the yield function takes the form:

$$\mathcal{F}(\mathit{\sigma },\mathbf{A},{\overline{\mathit{\epsilon }}}^p)=\frac{1}{2}\mathit{\sigma }:\mathbb{K}:\mathit{\sigma }+\mathbf{L}:\mathit{\sigma }-1\le 0,$$

(13)

where

$$\mathbb{K}:={\zeta }_1{\mathbb{P}}^{\mathrm{pind}}+({\zeta }_2-{\zeta }_1){\mathbb{P}}_{\mathbf{A}}^{\mathrm{pind}}+2{\zeta }_4(\mathbf{1}-\mathbf{A})\otimes (\mathbf{1}-\mathbf{A});\mathbf{L}:={\zeta }_3\left(\mathbf{1}-\mathbf{A}\right),$$

(14)

where the operators ${\mathbb{P}}^{\mathrm{pind}}$ and ${\mathbb{P}}_{\mathbf{A}}^{ind}$ are defined as:

$${\mathbb{P}}^{\mathrm{pind}}=\mathbb{I}-\frac{1}{2}\left(\mathbf{1}\otimes \mathbf{1}\right)+\frac{1}{2}\left(\mathbf{A}\otimes \mathbf{1}+\mathbf{1}\otimes \mathbf{A}\right)-\frac{3}{2}\left(\mathbf{A}\otimes \mathbf{A}\right);{\mathbb{P}}_{\mathbf{A}}^{\mathrm{pind}}:={\mathbb{P}}_{\mathbf{A}ijkl}^{\mathrm{pind}}={\mathbf{A}}_{im}{\mathbb{P}}_{mjkl}^{\mathrm{pind}}+{\mathbf{A}}_{mj}{\mathbb{P}}_{imkl}^{\mathrm{pind}}.$$

(15)

In comparison to the six-parameter yield surface definition in [15], herein a four-parameter yield surface is proposed. The herein proposal allows for an easier calibration of the yield surface and reduces the experimental effort. Nevertheless, the six-parameter yield function definition regards a better description of biaxial stress states which is crucial for accurate modeling of FRPs undergoing high hydrostatic pressures. This is achieved in [15] via the case differentiation concerning the invariant $I_3$ based on its sign.

2.4. Plastic Potential Function

To predict realistic plastic deformations, a non-associative flow rule is assumed. The construction of a non-associative transversely isotropic plastic potential function $\mathcal{G}(\mathit{\sigma },\mathbf{A})$ yields:

$$\mathcal{G}(\mathit{\sigma },\mathbf{A})={\varsigma }_1{I}_1+{\varsigma }_2{I}_2+{\varsigma }_3{I}_3^2-1,$$

(16)

where ${\varsigma }_i$ ($i=1,3$) denotes the plastic potential parameters. A condensed expression of the plastic flow potential is given by:

$$\mathcal{G}(\mathit{\sigma },\mathbf{A})=\frac{1}{2}\mathit{\sigma }:\mathbb{M}:\mathit{\sigma }-1\le 0,$$

(17)

where the fourth-order tensor $\mathbb{M}$ is expressed as:

$$\mathbb{M}:={\varsigma }_1{\mathbb{P}}^{\mathrm{pind}}+\left({\varsigma }_2-{\varsigma }_1\right){\mathbb{P}}_{\mathbf{A}}^{\mathrm{pind}}+2{\varsigma }_3\left(\mathbf{1}-\mathbf{A}\right)\otimes \left(\mathbf{1}-\mathbf{A}\right).$$

(18)

2.5. Evolution Equations

The evolution equations of the internal variables (${\mathit{\epsilon }}^p$ and $\mathit{\varpi }$) are expressed as follows:

$${\dot{\mathit{\epsilon }}}^p=\dot{\gamma }\frac{\partial \mathcal{G}(\mathit{\sigma },\mathbf{A},{\overline{\epsilon }}^p)}{\partial \mathit{\sigma }}=\dot{\gamma }{\mathbf{n}}_{\mathcal{G}}=\dot{\gamma }\mathbb{M}:\mathit{\sigma }\mathrm{with}{\mathbf{n}}_{\mathcal{G}}=\mathbb{M}:\mathit{\sigma },$$

(19)

$$\dot{\mathit{\varpi }}=\dot{\gamma }\frac{\partial \mathcal{G}(\mathit{\sigma },\mathbf{A},{\overline{\epsilon }}^p)}{\partial \mathsf{\Gamma }},$$

(20)

where $\gamma$ represents the so-called plastic multiplier.

As customary, the Kuhn–Tucker loading/unloading conditions are defined by:

$$\dot{\gamma }\ge 0;\mathcal{F}(\mathit{\sigma },\mathbf{A},{\overline{\epsilon }}^p)\le 0;\dot{\gamma }\mathcal{F}(\mathit{\sigma },\mathbf{A},{\overline{\epsilon }}^p)=0,$$

(21)

and the consistency condition as:

$$\dot{\gamma }\dot{\mathcal{F}}(\mathit{\sigma },\mathbf{A},{\overline{\epsilon }}^p)=0.$$

(22)

2.6. Parameter Identification

In addition to the elastic material constants, the yield function parameters ${\zeta }_i$ ($i=1,4$) and the plastic potential parameters ${\varsigma }_i$ ($i=1,3$) are to be determined.

The parameters ${\zeta }_i$ ($i=1,4$) control the size and shape of the elastic domain $\mathbb{E}$ as a function of the equivalent plastic strain variable ${\overline{\epsilon }}^p$. For each parameter, the relation ${\zeta }_i\left({\overline{\epsilon }}^p\right)$ is determined from an independent experiment, thus a total of four different experiments is required for calibration. For instance, the following four experiments can be employed for calibration: (i) in-plane shear test, (ii) transverse shear test, (iii) uniaxial transverse tension test, and (iv) uniaxial transverse compression test. The corresponding yield stress states are denoted as ${\sigma }_{is}^y,{\sigma }_{ts}^y,{\sigma }_{tt}^y$, and ${\sigma }_{tc}^y$, respectively. Similar to the procedure in [15], the four parameters ${\zeta }_i({\sigma }_{is}^y,{\sigma }_{ts}^y,{\sigma }_{tt}^y,{\sigma }_{tc}^y)$ ($i=1,4$) can then be obtained by entering the stress states from each experiments above in Equation (10) and setting the yield function state to yielding i.e., $\mathcal{F}=0$. Accordingly, the coefficients ${\zeta }_i$ ($i=1,4$) are explicitly given in the following.

From the in-plane shear test the first coefficient ${\zeta }_1$ is expressed as:

$${\zeta }_1=\frac{1}{{{\sigma }_{ts}^y}^2},$$

(23)

and from the transverse shear test the second coefficient ${\zeta }_2$ is given by:

$${\zeta }_2=\frac{1}{{{\sigma }_{is}^y}^2}.$$

(24)

The third coefficient ${\zeta }_3$ controls the tension-compression yield asymmetry and therefore is expressed in terms of the uniaxial transverse tension and uniaxial transverse compression tests as:

$${\zeta }_3=-\frac{1}{{\sigma }_{tc}^y}+\frac{1}{{\sigma }_{tt}^y},$$

(25)

Lastly, the coefficient ${\zeta }_4$ is associated with transverse loading, hence is expressed as:

$${\zeta }_4=-\frac{1}{4{{\sigma }_{ts}^y}^2}+\frac{1}{{\sigma }_{tc}^y{\sigma }_{tt}^y}.$$

(26)

To comply with the maximum dissipation principle, the convexity of the yield surface must be insured. This imposes the following restrictions to the relations ${\zeta }_i\left({\overline{\epsilon }}^p\right)$ ($i=1,4$) which must hold for any ${\overline{\epsilon }}^p$:

$${\sigma }_{tt}^y\le \frac{4{{\sigma }_{ts}^y}^2}{{\sigma }_{tc}^y}.$$

(27)

Similary, the parameters ${\varsigma }_i$ ($i=1,3$) control the size and shape of the plastic potential surface. However, one of these parameters is a scaling parameter and can be set to any value since the size of the plastic potential has no inherent physical meaning. Accordingly, there are only two remaining parameters to be determined and to associate with experimental data. In the present case, ${\varsigma }_1$ is arbitrarily set to unity.

As mentioned above, the motive behind adopting a non-associative plasticity scheme is to model realistic plastic deformation behavior as compared to associative plasticity. Accordingly, the parameters ${\varsigma }_i$ ($i=2,3$) are used to enforce certain plastic Poisson’s ratios ${\nu }_{23}^p={\epsilon }_{22}^p/{\epsilon }_{33}^p$ and plastic distortion behavior through the relation ${\mu }_{12}^p={\epsilon }_{12}^p/{\epsilon }_{23}^p$:

$${\varsigma }_1=1,$$

(28)

$${\varsigma }_2={\mu }_{12}^p,$$

(29)

$${\varsigma }_3=\frac{-1+{\nu }_{23}^p}{4(1+{\nu }_{23}^p)}.$$

(30)

Similarly, for the plastic potential function $\mathcal{G}$, the following must hold:

$${\mu }_{12}^p\ge 0\wedge -\frac{-1+{\nu }_{23}^p}{4(1+{\nu }_{23}^p)}\ge 0.$$

(31)

In contrast to the time-consuming iterative procedure presented in [15] for the determination of the plastic potential parameters, herein explicit expressions for the parameters are provided.

3. Numerical Treatment

In this section, the numerical treatment of the constitutive model proposed in Section 2 is discussed.

The construction of a numerical scheme for the solution of the initial boundary value problem (IBVP) associated with the current elasto–plastic model involved two main aspects [32]. The first concerned the local (at the Gauss point in FE context) integration of the evolution equations. The second regarded the employment of the result stemming from the previous step in the constitutive block of the weak formulation of the balance of linear momentum, which was discretized in space by means of FEM and solved by means of a standard incremental-iterative Newton–Raphson scheme.

It should be noted that all quantities presented in this section are computed in the rotated frame ${\mathcal{B}}_{n+1}^{\mathrm{rot}}$.

3.1. Numerical Integration: General Return Mapping Algorithm

For a prescribed motion of an arbitrary body, let us consider the time interval $[t_n,t_{n+1}^{\left(i\right)}]$, with $t\in {\mathbb{R}}_{+}$, where $t_n$ identifies the previous converged time step and $t_{n+1}^{\left(i\right)}$ denotes the current prospective time step at the global Newton–Raphson iteration i. The strain rate within the time step were given by:

$$\dot{\mathit{\epsilon }}=\frac{{\mathit{\epsilon }}_{n+1}-{\mathit{\epsilon }}_n}{\mathsf{\Delta }t};\mathrm{with}\mathsf{\Delta }t=t_{n+1}-t_n.$$

(32)

To simplify the notation, the superscript i is omitted.

The internal variables ${\mathit{\epsilon }}_n^p$, ${\overline{\epsilon }}_n^p$ and ${\mathit{\varpi }}_n$, and the prospective total strain ${\mathit{\epsilon }}_{n+1}$ are assumed to be available. Then, the elasto–plastic constitutive boundary value problem at the material (Gauss) point level is stated as follows:

Given: ${\mathit{\epsilon }}_n^p$, ${\overline{\epsilon }}_n^p$, ${\mathit{\varpi }}_n$, and ${\mathit{\epsilon }}_{n+1}$,

Find: ${\mathit{\epsilon }}_{n+1}^p$, ${\overline{\epsilon }}_{n+1}^p$, and ${\mathit{\varpi }}_{n+1}$ at the end of the time interval $\left[t_n,t_{n+1}\right]$,

Such that:

$${\dot{\mathit{\epsilon }}}^e=\dot{\mathit{\epsilon }}-\dot{\gamma }{\mathbf{n}}_{\mathcal{G}};{\dot{\overline{\epsilon }}}^p=\dot{\gamma }\sqrt{\frac{1}{2}}\parallel {\mathbf{n}}_{\mathcal{G}}\parallel ,$$

(33)

with

$$\dot{\gamma }\ge 0;\mathcal{F}(\mathit{\sigma },\mathbf{A},{\overline{\epsilon }}^p)\le 0;\dot{\gamma }\mathcal{F}(\mathit{\sigma },\mathbf{A},{\overline{\epsilon }}^p)=0.$$

(34)

The central point for the local integration of the model is the adoption of the backward-Euler (fully implicit, first-order accurate and unconditionally stable) integration scheme. Accordingly, the discrete version of the rate expressions given in Equations (32) and (33) within the interval $[t_n,t_{n+1}]$ are obtained as follows:

$${\mathit{\epsilon }}_{n+1}^e={\mathit{\epsilon }}_n^e+\mathsf{\Delta }\mathit{\epsilon }-{\gamma }_{n+1}{\mathbf{n}}_{\mathcal{G},n+1};{\overline{\epsilon }}_{n+1}^p={\overline{\epsilon }}_n^p+{\gamma }_{n+1}\sqrt{\frac{1}{2}}\parallel {\mathbf{n}}_{\mathcal{G},n+1}\parallel ,$$

(35)

with

$${\gamma }_{n+1}\ge 0;\mathcal{F}({\mathit{\sigma }}_{n+1},\mathbf{A},{\overline{\epsilon }}_{n+1}^p)\le 0;{\gamma }_{n+1}\mathcal{F}({\mathit{\sigma }}_{n+1},\mathbf{A},{\overline{\epsilon }}_{n+1}^p)=0,$$

(36)

where $\mathsf{\Delta }\mathit{\epsilon }={\mathit{\epsilon }}_{n+1}-{\mathit{\epsilon }}_n$.

Next, the classical two-step predictor-corrector procedure [33] is applied. The first step concerns the computation of the predictor elastic trial step as follows:

$${\mathit{\epsilon }}_{n+1}^{e,\mathrm{trial}}={\mathit{\epsilon }}_n^e+\mathsf{\Delta }\mathit{\epsilon }\mathrm{and}{\overline{\epsilon }}_{n+1}^{p,\mathrm{trial}}={\overline{\epsilon }}_n^p,$$

(37)

$${\mathit{\sigma }}_{n+1}^{\mathrm{trial}}={\mathbb{C}}_e:{\mathit{\epsilon }}_{n+1}^{\mathrm{e},\mathrm{trial}}.$$

(38)

The corresponding trial yield function is given by:

$$\mathcal{F}({\mathit{\sigma }}_{n+1}^{\mathrm{trial}},\mathbf{A},{\overline{\epsilon }}_n^p)=\frac{1}{2}{\mathit{\sigma }}_{n+1}^{\mathrm{trial}}:{\mathbb{K}}^{\mathrm{trial}}:{\mathit{\sigma }}_{n+1}^{\mathrm{trial}}+{\mathbf{L}}^{\mathrm{trial}}:{\mathit{\sigma }}_{n+1}^{\mathrm{trial}}-1,$$

(39)

where the operators ${\mathbb{K}}^{\mathrm{trial}}$ and ${\mathbf{L}}^{\mathrm{trial}}$ are function of the trial equivalent plastic strain ${\overline{\epsilon }}_{n+1}^{p,\mathrm{trial}}$.

As customary, if the elastic trial state lies within the elastic domain i.e., $\mathcal{F}({\mathit{\sigma }}_{n+1}^{\mathrm{trial}},\mathbf{A},{\overline{\epsilon }}_n^p)<0$, then the solution is elastic with ${\gamma }_{n+1}=0$ and the trial step is accepted as the correct solution. Otherwise, the solution is plastic with ${\gamma }_{n+1}>0$ and is obtained via the plastic corrector step fulfilling the constraint:

$${\mathcal{F}}_{n+1}({\mathit{\sigma }}_{n+1},\mathbf{A},{\overline{\epsilon }}_{n+1}^p)\stackrel{!}{=}0.$$

(40)

Based on this, the computation of the plastic multiplier ${\gamma }_{n+1}$ follows the procedure outlined in Algorithm 1.

Algorithm 1 Plastic corrector step: algorithmic computation of the plastic multiplier and update of the internal variables.

Compute ${\mathit{\epsilon }}_{n+1}^{e,\mathrm{trial}}={\mathit{\epsilon }}_{n+1}^e+{\gamma }_{n+1}{\mathit{n}}_{\mathcal{G},n+1}$. Substitute ${\mathit{n}}_{\mathcal{G},n+1}=\mathbb{M}:{\mathit{\sigma }}_{n+1}$→${\mathit{\epsilon }}_{n+1}^{e,\mathrm{trial}}={\mathit{\epsilon }}_{n+1}^e+{\gamma }_{n+1}\mathbb{M}:{\mathit{\sigma }}_{n+1}$. Compute ${\mathbb{C}}^e:{\mathit{\epsilon }}_{n+1}^{e,\mathrm{trial}}={\mathbb{C}}^e:{\mathit{\epsilon }}_{n+1}^e+{\gamma }_{n+1}{\mathbb{C}}^e:\mathbb{M}:{\mathit{\sigma }}_{n+1}$. Identify ${\mathit{\sigma }}_{n+1}^{\mathrm{trial}}={\mathit{\sigma }}_{n+1}+{\gamma }_{n+1}{\mathbb{C}}^e:\mathbb{M}:{\mathit{\sigma }}_{n+1}$. Compute ${\mathit{\sigma }}_{n+1}={\left[\mathbb{I}+{\gamma }_{n+1}{\mathbb{C}}^e:\mathbb{M}\right]}^{-1}:{\sigma }_{n+1}^{\mathrm{trial}}=\mathbb{H}:{\sigma }_{n+1}^{\mathrm{trial}}$; with $\mathbb{H}={\left[\mathbb{I}+{\gamma }_{n+1}{\mathbb{C}}^e:\mathbb{M}\right]}^{-1}$. Solve the equation to determine the consistency parameter ${\gamma }_{n+1}$ via local iterative process (local Newton–Raphson index denoted by the superscript k). (a) Set $k=0$ and the initial values $({\mathit{\sigma }}_{n+1}^{(k=0)}={\mathit{\sigma }}_{n+1}^{\mathrm{trial}},{\overline{\epsilon }}_{n+1}^{p,(k=0)}={\overline{\epsilon }}_n^p,{\gamma }_{n+1}^{(k=0)}=0)$.  (b) Compute ${\mathcal{F}}^{\left(k\right)}({\mathit{\sigma }}_{n+1},\mathbf{A},{\overline{\epsilon }}_{n+1}^p)$.  (c) IF ${\mathcal{F}}^{\left(k\right)}({\mathit{\sigma }}_{n+1},\mathbf{A},{\overline{\epsilon }}_{n+1}^p)\le \mathrm{TOL}$ GOTO 7, ELSE  (d) Set residual for local Newton-Rapshon iteration ${\mathcal{R}}_{n+1}^{\left(k\right)}={\mathcal{F}}^{\left(k\right)}({\mathit{\sigma }}_{n+1},\mathbf{A},{\overline{\epsilon }}_{n+1}^p)$.  (e) Perform linearization of ${\mathcal{R}}_{n+1}^{\left(k\right)}$: $\widehat{L}\left[{\mathcal{R}}_{n+1}^{\left(k\right)}\right]\simeq {\mathcal{R}}_{n+1}^{\left(k\right)}+\mathsf{\Delta }{\gamma }^{\left(k\right)}\left[\frac{\partial {\mathcal{F}}_{n+1}^{\left(k\right)}}{\partial {\mathit{\sigma }}_{n+1}^{\left(k\right)}}:\frac{\partial {\mathit{\sigma }}_{n+1}^{\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}+\frac{\partial {\mathcal{F}}_{n+1}^{\left(k\right)}}{\partial {\mathbb{K}}_{n+1}^{\left(k\right)}}····\frac{\partial {\mathbb{K}}_{n+1}^{\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}+\frac{\partial {\mathcal{F}}_{n+1}^{\left(k\right)}}{\partial {\mathbf{L}}_{n+1}^{\left(k\right)}}:\frac{\partial {\mathbf{L}}_{n+1}^{\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}\right]=0$.  (f) Compute ${\displaystyle \mathsf{\Delta }{\gamma }^{\left(k\right)}=\frac{-{\mathcal{R}}_{n+1}^{\left(k\right)}}{\frac{\partial {\mathcal{F}}_{n+1}^{\left(k\right)}}{\partial {\mathit{\sigma }}_{n+1}^{\left(k\right)}}:\frac{\partial {\mathit{\sigma }}_{n+1}^{\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}+\frac{\partial {\mathcal{F}}_{n+1}^{\left(k\right)}}{\partial {\mathbb{K}}_{n+1}^{\left(k\right)}}····\frac{\partial {\mathbb{K}}_{n+1}^{\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}+\frac{\partial {\mathcal{F}}_{n+1}^{\left(k\right)}}{\partial {\mathbf{L}}_{n+1}^{\left(k\right)}}:\frac{\partial {\mathbf{L}}_{n+1}^{\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}}}$.  (g) Correct ${\gamma }_{n+1}^{(k+1)}={\gamma }_{n+1}^{\left(k\right)}+\mathsf{\Delta }{\gamma }^{\left(k\right)}$.  (h) $k\leftarrow k+1$ GOTO (b)   Update the internal variables ${\mathit{\sigma }}_{n+1}={\mathit{\sigma }}_{n+1}^{\left(k\right)}$, ${\overline{\epsilon }}_{n+1}^p={\overline{\epsilon }}_{n+1}^{p,\left(k\right)}$, ${\mathit{\epsilon }}_{n+1}^p={\mathit{\epsilon }}_{n+1}^{p,\left(k\right)}$. Compute algorithmic tangent operator, see Section 3.2.

The expressions required for the computation of Algorithm 1 are provided in the following. The first term within the denominator of the linearization $\frac{\partial {\mathcal{F}}_{n+1}^{\left(k\right)}}{\partial {\mathit{\sigma }}_{n+1}^{\left(k\right)}}$ takes the form:

$$\frac{\partial {\mathcal{F}}_{n+1}^{\left(k\right)}}{\partial {\mathit{\sigma }}_{n+1}^{\left(k\right)}}={\mathbb{K}}_{n+1}^{\left(k\right)}:{\mathit{\sigma }}_{n+1}^{\left(k\right)}+{\mathbf{L}}_{n+1}^{\left(k\right)},$$

(41)

where

$${\mathit{\sigma }}_{n+1}^{\left(k\right)}={\mathbb{H}}_{n+1}^{\left(k\right)}:{\mathit{\sigma }}_{n+1}^{\mathrm{trial}},$$

(42)

and

$${\overline{\epsilon }}_{n+1}^{p,\left(k\right)}={\overline{\epsilon }}_n^{p,\left(k\right)}+{\gamma }_{n+1}^{\left(k\right)}\sqrt{\frac{1}{2}}\parallel {\mathbb{M}}_{n+1}:{\mathit{\sigma }}_{n+1}^{\left(k\right)}\parallel .$$

(43)

The second term $\frac{\partial {\mathit{\sigma }}_{n+1}^{\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}$ is expressed as:

$$\frac{\partial {\mathit{\sigma }}_{n+1}^{\left(k\right)}}{\partial {\gamma }_{n+1}}=-{\mathbb{H}}_{n+1}^{\left(k\right)}:\left[\left({\mathbb{C}}_e:{\mathbb{M}}_{n+1}\right):{\mathit{\sigma }}_{n+1}^{\left(k\right)}\right].$$

(44)

The third and fifth terms take the form, respectively:

$$\frac{\partial {\mathcal{F}}_{n+1}^{\left(k\right)}}{\partial {\mathbb{K}}_{n+1}^{\left(k\right)}}=\frac{1}{2}\left[{\mathit{\sigma }}_{n+1}^{\left(k\right)}\otimes {\mathit{\sigma }}_{n+1}^{\left(k\right)}\right];\frac{\partial {\mathcal{F}}_{n+1}^{\left(k\right)}}{\partial {\mathbf{L}}_{n+1}^{\left(k\right)}}={\mathit{\sigma }}_{n+1}^{\left(k\right)}.$$

(45)

The fourth term reads:

$$\frac{\partial {\mathbb{K}}_{n+1}^{\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}=\frac{\partial {\mathbb{K}}_{n+1}^{\left(k\right)}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}\frac{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}.$$

(46)

The term $\frac{\partial {\mathbb{K}}_{n+1}^{\left(k\right)}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}$ in Equation (46) is expressed as:

$$\frac{\partial {\mathbb{K}}_{n+1}^{\left(k\right)}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}=\sum _{i=1,2,4}\frac{\partial {\mathbb{K}}_{n+1}^{\left(k\right)}}{\partial {\zeta }_{\mathrm{i}}^{\left(k\right)}}\frac{\partial {\zeta }_{\mathrm{i}}^{\left(k\right)}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}=\left[{\mathbb{P}}^{\mathrm{ind}}-{\mathbb{P}}_{\mathbf{A}}^{\mathrm{ind}}\right]\frac{\partial {\zeta }_1^{\left(k\right)}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}+{\mathbf{P}}_{\mathbf{A}}^{\mathrm{ind}}\frac{\partial {\zeta }_2^{\left(k\right)}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}+2\left(\mathbf{1}-\mathbf{A}\right)\otimes \left(\mathbf{1}-\mathbf{A}\right)\frac{\partial {\zeta }_4^{\left(k\right)}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}},$$

(47)

where

$$\frac{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}=\sqrt{\frac{1}{2}}\parallel {\mathbb{M}}_{n+1}:{\mathit{\sigma }}_{n+1}^{\left(k\right)}\parallel +{\gamma }_{n+1}^{\left(k\right)}\sqrt{\frac{1}{2}}\left[\frac{\left[{\mathbb{M}}_{n+1}:{\mathit{\sigma }}_{n+1}^{\left(k\right)}\right]:{\mathbb{M}}_{n+1}}{\parallel {\mathbb{M}}_{n+1}:{\mathit{\sigma }}_{n+1}^{\left(k\right)}\parallel }\right]:\frac{\partial {\mathit{\sigma }}_{n+1}^{\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}.$$

(48)

Lastly, the sixth terms $\frac{\partial {\mathbf{L}}_{n+1}^{\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}$ takes the form:

$$\frac{\partial {\mathbf{L}}_{n+1}^{\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}=\frac{\partial {\mathbf{L}}_{n+1}^{\left(k\right)}\left({\zeta }_3^{\left(k\right)}\right)}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}\frac{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}{\partial {\gamma }_{n+1}^{\left(k\right)}}.$$

(49)

where

$$\frac{\partial {\mathbf{L}}_{n+1}^{\left(k\right)}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}=\frac{\partial {\mathbf{L}}_{n+1}^{\left(k\right)}}{\partial {\zeta }_3^{\left(k\right)}}\frac{\partial {\zeta }_3^{\left(k\right)}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}=\left(\mathbf{1}-\mathbf{A}\right)\frac{\partial {\zeta }_3^{\left(k\right)}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^{p,\left(k\right)}}.$$

(50)

3.2. Algorithmic Consistent Tangent Moduli

For the solution of the non-linear FE equations (discretized weak form of the balance of linear momentum) on a global level, the incremental-iterative Newton–Raphson scheme is used [32]. Therein, in order to obtain a quadratic convergence, the computation of the algorithmic consistent tangent moduli is required, i.e., consistent with the chosen algorithmic time integration scheme.

The form, $\mathrm{d}{\mathit{\sigma }}_{n+1}={\mathbb{C}}_{n+1}^{\mathrm{ep}}:\mathrm{d}{\mathit{\epsilon }}_{n+1}$ describes the sensitivity of the stress with respect to an infinitesimal increment in the strain at time $t_{n+1}$ When the local integration algorithm described has converged is looked for.

The starting point to derive the algorithmic consistent tangent moduli is forming an expression for the infinitesimal increment of the total stress at time $t_{n+1}$. Using the relation ${\mathit{\sigma }}_{n+1}=\mathbb{H}:{\sigma }_{n+1}^{\mathrm{trial}}$ in Algorithm 1, the increment of the total stress reads:

$$\mathrm{d}{\mathit{\sigma }}_{n+1}={\mathbb{H}}_{n+1}:\left[{\mathbb{C}}^{\mathrm{e}}:\mathrm{d}{\mathit{\epsilon }}_{n+1}-\mathrm{d}{\gamma }_{n+1}\left[{\mathbb{C}}^e:{\mathbb{M}}_{n+1}\right]:{\mathit{\sigma }}_{n+1}\right],$$

(51)

Next, an explicit expression for the differential of the plastic multiplier $\mathrm{d}{\gamma }_{n+1}$ is to be obtained. This is achieved through the consistency condition given in Equation (22). In case of plastic loading i.e., ${\gamma }_{n+1}\ge 0$, ${\dot{\mathcal{F}}}_{n+1}=0$ and therefore $\mathrm{d}{\mathcal{F}}_{n+1}=0$. Accordingly from the condition $\mathrm{d}{\mathcal{F}}_{n+1}=0$, $\mathrm{d}{\gamma }_{n+1}$ is obtained as:

$$\mathrm{d}{\gamma }_{n+1}=-\frac{\frac{\partial {\mathcal{F}}_{n+1}^{*}}{\partial {\mathit{\epsilon }}_{n+1}}:\mathrm{d}{\mathit{\epsilon }}_{n+1}}{\frac{\partial {\mathcal{F}}_{n+1}^{*}}{\partial {\gamma }_{n+1}}},$$

(52)

where the term $\frac{\partial {\mathcal{F}}_{n+1}^{*}}{\partial {\mathit{\epsilon }}_{n+1}}$ reads:

$$\frac{\partial {\mathcal{F}}_{n+1}^{*}}{\partial {\mathit{\epsilon }}_{n+1}}=\frac{\partial {\mathcal{F}}_{n+1}}{\partial {\mathit{\sigma }}_{n+1}}:\left[{\mathbb{H}}_{n+1}:{\mathbb{C}}^e\right]+\frac{\partial {\mathcal{F}}_{n+1}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^p}\sqrt{\frac{1}{2}}\frac{{\gamma }_{n+1}}{∥{\mathbb{M}}_{n+1}:{\mathit{\sigma }}_{n+1}∥}\left[\left[{\mathbb{M}}_{n+1}:{\mathit{\sigma }}_{n+1}\right]:{\mathbb{M}}_{n+1}\right]:\left[{\mathbb{H}}_{n+1}:{\mathbb{C}}^e\right],$$

(53)

where

$$\frac{\partial {\mathcal{F}}_{n+1}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^p}=\frac{\partial {\mathcal{F}}_{n+1}}{\partial {\mathbb{K}}_{n+1}}····\frac{\partial {\mathbb{K}}_{n+1}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^p}+\frac{\partial {\mathcal{F}}_{n+1}}{\partial {\mathbf{L}}_{n+1}}:\frac{\partial {\mathbf{L}}_{n+1}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^p}.$$

(54)

The term $\frac{\partial {\mathcal{F}}_{n+1}^{*}}{\partial {\gamma }_{n+1}}$ takes the form:

$$\begin{array}{c}\frac{\partial {\mathcal{F}}_{n+1}^{*}}{\partial {\gamma }_{n+1}}=-\frac{\partial {\mathcal{F}}_{n+1}}{\partial {\mathit{\sigma }}_{n+1}}:\left[{\mathbb{H}}_{n+1}:\left[\left[{\mathbb{C}}^e:{\mathbb{M}}_{n+1}\right]:{\mathit{\sigma }}_{n+1}\right]\right]\\ +\frac{\partial {\mathcal{F}}_{n+1}}{\partial {\overline{\mathit{\epsilon }}}_{n+1}^p}\sqrt{\frac{1}{2}}\left[∥{\mathbb{M}}_{n+1}:{\mathit{\sigma }}_{n+1}∥-\frac{{\gamma }_{n+1}}{∥{\mathbb{M}}_{n+1}:{\mathit{\sigma }}_{n+1}∥}\left[\left[{\mathbb{M}}_{n+1}:{\mathit{\sigma }}_{n+1}\right]:{\mathbb{M}}_{n+1}\right]:\left[{\mathbb{H}}_{n+1}:\left[\left[{\mathbb{C}}^e:{\mathbb{M}}_{n+1}\right]:{\mathit{\sigma }}_{n+1}\right]\right]\right].\end{array}$$

(55)

Finally, by substituting the expression for $\mathrm{d}{\gamma }_{n+1}$ in Equation (51), the algorithmic consistent tangent moduli ${\mathbb{C}}_{n+1}^{ep}$ is given by:

$${\mathbb{C}}_{n+1}^{ep}=\frac{\partial {\mathit{\sigma }}_{n+1}}{\partial {\mathit{\epsilon }}_{n+1}}={\mathbb{H}}_{n+1}:\left[{\mathbb{C}}^e+\left[\left[{\mathbb{C}}^e:{\mathbb{M}}_{n+1}\right]:{\mathit{\sigma }}_{n+1}\right]\otimes \frac{\frac{\partial {\mathcal{F}}_{n+1}^{*}}{\partial {\mathit{\epsilon }}_{n+1}}}{\frac{\partial {\mathcal{F}}_{n+1}^{*}}{\partial {\gamma }_{n+1}}}\right].$$

(56)

3.3. FE Implementation in ABAQUS

Herein, the numerical implementation of the model in the general purpose FE code ABAQUS/Implicit via the user-defined subroutine UMAT is described.

During the global computation, the subroutine UMAT was called at all material calculation points of elements for which the material definition includes a user-defined material behavior. The subroutine must update the stress ($\mathit{\sigma }$) and solution-dependent state (internal) variables (${\mathit{\epsilon }}^p$ and ${\overline{\epsilon }}^p$) to their values at the end of the increment for which it is called and also provide the material Jacobian matrix (${\mathbb{C}}^{ep}$), see [34].

The incremental strain ($\mathsf{\Delta }\mathit{\epsilon }$) and the total strain (${\mathit{\epsilon }}_{n+1}$) in the rotated frame were passed in by the UMAT and their components are rotated to account for rigid body motion in the increment before UMAT was called.

The stress at the beginning of the increment (${\mathit{\sigma }}_n$) is also passed in. The stress is already rotated to account for rigid body motion in the increment and must be updated in the routine to be the stress at the end of the increment (${\mathit{\sigma }}_{n+1}$). For this reason, only the co-rotational part of the stress integration should be computed in UMAT as described above.

One major concern is the solution-dependent state variables. These variables are also passed in as the values at the beginning of the increment (${\mathit{\epsilon }}_n^p$ and ${\overline{\epsilon }}_n^p$). However, the vector-valued or tensor-valued internal variables (e.g., ${\mathit{\epsilon }}_n^p$) must be rotated to account for rigid body motion of the material in the increment. For this purpose, the rotation increment tensor (the increment of rigid body rotation of the element local co-rotational coordinate system) is also passed in so that the passed in vector- or tensor-valued internal variables are rotated appropriately in the UMAT subroutine (see [24] for the computation of the rotation increment tensor). Thereafter, the state variables must be updated based on the constitutive behavior to their values at the end of the increment (${\mathit{\epsilon }}_{n+1}^p$ and ${\overline{\epsilon }}_{n+1}^p$).

4. Representative Applications

The previously described formulation is implemented into ABAQUS/Implicit by means of the user-defined subroutine UMAT. In reference [35], the model is calibrated for carbon fibre reinforced polymer (CFRP) IM7/8552 carbon/epoxy using test data from experiments on UD laminates. Furthermore, the performance of the elasto–plastic model is verified and validated via the FE simulation of the characterization tests performed in reference [36].

In the following, two numerical examples at two different scales are presented in order to demonstrate the applicability and capability of the proposed development in the context of geometrical non-linear analysis of composites. The examples discussed in the sequel are: (i) micro-buckling of UD composites subjected to compressive loading, and (ii) structural application involving laminated composites cylinder with free edges subjected to a point load. In these examples, mesh and time step convergence studies are carried out to ensure the validity of the results. In the time step convergence study, in each study, the maximum step size is controlled.

4.1. Micro-Buckling

To assess the capabilities of the current model at micro-scale, the failure under axial compression of unidirectional glass fiber reinforced polymers (GFRP) E-Glass/MY750 glass/epoxy ply with 60% fibre volume fraction is considered. Fibres naturally show a sinusoidal misalignment in continuous unidirectional fibre reinforced polymers [37] resulting in geometrical non-linearities. These geometrical non-linearities needed to be considered in the modeling, along with obvious material non-linearities, as this defines the accurate prediction of the UD compressive behavior. The schematic representation of the fibre waviness in the model is shown in Figure 1 along with its boundary and loading conditions. The model is a 3D homogenized representation of a layer of 15 glass fibres from a unidirectional ply to show the effect of misaligned fibres on the failure under compression, termed as kinking or micro-buckling in literature [38]. Overall fibre lengths of 500 $\mathsf{\mu }$m are modeled, whereas width and thickness of the model come naturally from the fibre volume fraction and the number of fibres considered and are 93.75 $\mathsf{\mu }$m and 6.25 $\mathsf{\mu }$m, respectively. It should be noted that for the prediction of the different competing mechanisms leading to final kinking failure under compression, the micro-mechanical approaches with separate fibre-matrix modeling are useful [38]. However, the global stress-strain response can accurately be obtained through the current approach, with the advantages of significantly easier modeling and higher degree of computational efficiency.

The geometrical non-linearity of the fibres is introduced as in-plane sinusoidal angular misalignment in the model following [39] to initiate a kink band. The sinus waviness is over a length of 85 $\mathsf{\mu }$m in the central region of the model with variable amplitudes, starting with an amplitude at one end and decreasing smoothly to an amplitude of 0 at the other end of the region bounded axially by $x_1\le x\le x_2$ in global x-direction as plotted in Figure 1. The fibre misalignment function is given below:

$$y=\left\{\begin{array}{c}(i-1)hx<x_1\hfill \\ (i-1)h+\lambda (1-\frac{i}{N})(1-cos\frac{\pi }{l}x)x_1\le x\le x_2\hfill \\ (i-1)h+2\lambda (1-\frac{i}{N})x>x_2,\hfill \end{array}\right.$$

(57)

where N is the number of fibres in thickness of the layer, h refers to the distance between the center of adjacent fibres, l denotes the half wavelength, $\lambda$ is the maximum value of amplitude, and $x_1$ and $x_2$ are the starting and ending positions of the waviness region, respectively.

A 3D finite element analysis (FEA) is performed to highlight the necessity of accounting for geometrical non-linearities at micro-level in the simulation of compressive failure of FRPs and to show the gained advantage of reduced computational costs through a homogenized modeling approach. The FE discretization consists of 9600 second-order, structured topology (3D 20-node brick elements—C3D20R).

The left face of the model is bounded in-plane i.e., global x- and y-axis, and the bottom left edge is bounded out-of-plane i.e., global z-axis. The right face of the model is coupled with a reference node through kinematic coupling, and axial force load is applied in negative x-direction. Since the kinking failure of unidirectional FRPs show a snap-back behavior, the riks method is used to capture the equilibrium path beyond limit points.

The material data needed for the model calibration are taken from reference [40]. The elastic material properties are reported in

Table 1. GFRP E-glass/MY750: elastic properties.

${\mathit{E}}_{11}$ (MPa)	${\mathit{E}}_{22}$ (MPa)	${\mathit{G}}_{12}$ (MPa)	${\mathit{\nu }}_{12}$	${\mathit{\nu }}_{23}$
$\mathrm{55,000}$	$\mathrm{45,600}$	$\mathrm{16,200}$	$0.0987$	$0.40$

. Beside the elastic material constants, utilizing Equations (23)–(26), the yield function parameters ${\zeta }_i$ ($i=1,4$) that characterize the onset of yielding are listed in

Table 2. GFRP E-Glass/MY750: yielding parameters ${\zeta }_i$ at the onset of yielding.

${\mathit{\zeta }}_1$	${\mathit{\zeta }}_2$	${\mathit{\zeta }}_3$	${\mathit{\zeta }}_4$
$0.00261641$	$0.00189036$	$0.0112808$	$0.000163349$

. Furthermore, the plastic potential function parameters ${\varsigma }_i$ ($i=1,3$) are provided in

Table 3. GFRP E-Glass/MY750: plastic potential parameters ${\varsigma }_i$.

${\mathit{\varsigma }}_1$	${\mathit{\varsigma }}_2$	${\mathit{\varsigma }}_3$
$1.0$	$1.0$	$-0.1071428$

. These values are determined based on the plastic Poisson’s ratio ${\nu }_{23}^p=0.4$ and plastic distortion ratio ${\mu }_{12}^p=1.0$. Due to the lack of experimental data concerning the transverse shear, reasonable assumptions were made for transverse shear behavior.

The results in Figure 2a show the axial compression response curve for the geometrically linear and non-linear cases. In the plot, the axial stresses are calculated by taking the ratio of the applied incremental load with the initial cross-sectional area. Whereas, the strains are calculated by the ratio of the axial end shortening to the initial micro-model length. Under the applied compressive load, the shear stress concentrates at the misalignment region resulting in shear yielding and a sudden drop in load carrying capacity because of the instability which is seen as snap-back in the equilibrium path. This point of instability corresponds to the peak load. The shear localization, in turn, rotates the already misaligned region and forms the so-called kink band. The kink band formation represented by the equivalent plastic strain is depicted in Figure 2b. For the E-Glass/MY750 material, the calculated compressive strength through geometrically non-linear analysis is 860 MPa whereas the measured strength according to reference [40] is 800 MPa. On the other hand, the geometrical linear analysis with the same parameters shows an unrealistically high strength value of 1800 MPa. Considering the stochastic nature of compressive strength and limited experimental data available, it can be concluded that the current formulation is able to predict the compressive behavior reasonably well. Another thing to note is the highly reduced numerical size of the problem along with simpler modeling due to the homogenized material representation.

Using the micro-mechanical modeling approach where fibres and matrix are modeled separately, the detailed mechanism of the compressive failure mode can be investigated and observed, see reference [38]. Employing both the micro-mechanical and the current homogenized approach shows the same qualitative global response. However, the current approach is much more numerically efficient as compared to the micro-mechanical approach. For example, for the same model dimensions, the micro-mechanical approach in reference [38] required 20 times more elements for FE discretization. It should be noted that a direct quantitative comparison of the results obtained employing the presented approach with the micro-mechanical approach in reference [38] is not possible here since the materials investigated are different.

4.2. Laminated Composites Cylinder under Point Loads

Herein, elasto–plastic co-rotational framework-based geometrical non-linear analysis of a cross-ply ${[0/90]}_{\mathrm{s}}$ IM7/8551-7 carbon/epoxy laminated cylinder with free edges subjected to two opposite point loads is presented. The geometric description of the cylinder, FE mesh, boundary conditions and loading are depicted in Figure 3. The dimensions of the cylinder are: (i) length L = 5000 mm, (ii) mid-surface radius R = 2470 mm, and (iii) thickness t = 60 mm.

The elastic and plastic material properties needed for model calibration are given in

Table 4. Carbon fibre reinforced polymer (CFRP) IM7/8551-7: elastic properties.

${\mathit{E}}_{11}$ (MPa)	c${\mathit{E}}_{22}$ (MPa)	${\mathit{G}}_{12}$ (MPa)	${\mathit{\nu }}_{12}$	${\mathit{\nu }}_{23}$
$\mathrm{165,000}$	8400	5600	$0.0173$	$0.50$

,

Table 5. CFRP IM7/8551-7: yielding parameters ${\zeta }_i$ at the onset of yielding.

${\mathit{\zeta }}_1$	${\mathit{\zeta }}_2$	${\mathit{\zeta }}_3$	${\mathit{\zeta }}_4$
$0.00176541$	$0.00127551$	$0.00926641$	$0.000110219$

and

Table 6. CFRP IM7/8551-7: plastic potential parameters ${\varsigma }_i$.

${\mathit{\varsigma }}_1$	${\mathit{\varsigma }}_2$	${\mathit{\varsigma }}_3$
$1.0$	$1.0$	$-0.08333333$

. Herein, the plastic Poisson’s ratio ${\nu }_{23}^p=0.5$ and plastic distortion ratio ${\mu }_{12}^p=1.0$. The 20-node quadratic brick element type C3D20R is used. After mesh convergence study, 31,600 elements are generated.

The load level and laminate stacking sequence are selected so that the strains remain small. The occurrence of material failure is checked by the invariant-based pressure-dependent quadratic asymmetric failure criteria (IQC) proposed in references [7,14].

The deformed configuration of the analyzed cylinder is shown in Figure 4a. The load-displacement diagrams at point A (directly under the load) is depicted in Figure 4b.

To point out the significance of including the geometrical non-linear effects, a geometrical linear elasto–plastic analysis is performed and the load-displacement diagram at point A is added to Figure 4. By comparing the load-displacement diagrams in Figure 4 obtained from the geometrical non-linear and geometrical linear analysis, a significant difference in the response of the structure under the same load level is observed. At point A, under the applied load, geometrically non-linear analysis resulted in a deflection of about is 0.83 m whereas the geometrical linear analysis under the same load level resulted in a deflection value of about 1.35 m.

In Figure 5 the fibre orientation (represented by the normal to the nominal fibre orientation) change in the outer ply predicted by the geometrical non-linear analysis is depicted. In this graph, a significant variation of the fibre direction throughout the process is estimated. This fact stems from the large displacements and rotations experienced by the cylinder, which can notably affect the performance in service and cannot be captured using a geometrically linear model. This becomes evident in the current example and highlights the necessity of triggering the evolution of the fibre orientation along the deformation process. This issue can be only performed using a geometrically non-linear setting.

5. Conclusions

This paper was focused on the co-rotational formulation of an invariant-based anisotropic elasto–plastic model including detailed aspects of its numerical treatment and implementation in the finite element (FE) framework for geometrically non-linear analysis of fibre reinforced polymers (FRPs).

The proposed plasticity formulation assumed a pressure-dependent yield surface and a non-associate flow rule to capture realistic evolution of the inelastic behavior. In comparison to the yield function definition in [15], herein a new definition of the yield function that eases the calibration procedure and reduces the experimental effort was proposed. Hence, explicit expressions for the determination of the model parameters were provided.

On the computational side, the full computational algorithm of the proposed model was developed. Locally, the integration of the model evolution equations was given. Therein, explicit expressions necessary for the algorithmic computation of the model variables were provided. Globally, the consistent algorithmic tangent moduli was derived. Moreover, the important aspects of the model implementation in the general purpose FE code ABAQUS/Implicit were discussed.

Finally, two numerical examples at two different scales were presented pointing out the relevance of including the geometrical non-linear effects in the finite element analysis of FRPs. One key aspect was the possibility to allow for finite fibre rotation concurrently with the deformation process and thus the change of the material orientation.

The development of realistic models for complex materials usually requires a combination of more than one basic dissipative phenomena. Quasi-brittle materials, like FRPs, show damage and plasticity at the same time. Therefore, coupling the proposed plasticity anisotropic formulation with damage in order to describe the interaction between these processes represents the upcoming research focus. Among the different available options for damage modeling, those associated with kinematic enrichment of the FE mesh represent an appealing candidate.