An alternative finite strain elastoplastic model applied to soft core sandwich panels simulation

Abstract

An alternative elastoplastic model based on the Flory’s right Cauchy-Green stretch tensor decomposition is proposed and applied to model soft cores of sandwich panels. It does not follow usual methodologies as the additive decomposition or the Kröner-Lee multiplicative decomposition of strains. It is based on an important hyperelastic relation, Flory’s decomposition, from which the total strain is separated in two isochoric and one volumetric parts. Using this decomposition, the volumetric strain energy continues to be elastic during all elastoplastic analysis and the isochoric parts are managed to produce the plastic evolution. As a consequence of Flory’s decomposition, the plastic flow direction is known and independent of the yielding surfaces. Moreover, it provides the well known deviatoric nature of plastic strains. For validation purposes, the resulting formulation is implemented using a special 3D prismatic element in a geometrical nonlinear positional FEM computational code. The achieved numerical results are compared with literature experimental data of soft core laminated structural elements.

Keywords:
Flory’s decomposition; Sandwich panels; Static core crushing; Finite strain elastoplasticity; Positional FEM

Graphical abstract

1 INTRODUCTION

The study of laminated beams, plates and shells is very extensive as it is motivated by the high importance of this type of structural arrangement in the mechanical, aeronautical, marine and civilindustries. The extent of this area is so large that the up-to-date literature reviews are divided into sub-areas. In order to find recent reviews readers are invited to consult, for example, the work of Abrate and Di Sciuva (2017) for general Equivalent-Single-Layer theories (ESL), Caliri et al. (2016) for FEM shell theories, Garg et al.(2021) to access theories including functionally graded materials, Nsengiyumva et al. (2021) for a review associated with non-destructive experimental analysis and the work of Sayyad and Ghugal (2017) for vibrations and stability of laminated beams and plates. Some current works of great interest and related to the present study can be also cited, for example, the work of Icardi and Urraci (2019) regarding sandwich panels vibration, Chen et al. (2016) related to impact analysis, Marques et al. (2018) for structural health monitoring and Naik and Sayyad (2018) and Shekari et al. (2017) for general plate theories.

The present study is dedicated to propose an elastoplastic model for the analysis of soft cores usually present in sandwich panels. When loads are applied over small areas of sandwich panels the small stiffness of face-sheets allows the core suffer high stress levels developing plastic strains that compromise the behavior of the whole panel (Koissin et al., 2004). This phenomenon is usually called core crushing and may occur in various situations as, for example, collisions, dropping objects, local buckling, overloaded connections etc. Some interesting works that studied the core crushing can be cited, for example Thomsen (1993), Soden (1996) and Olsson(2002) did important analytical developments and Lolive and Berthelot (2002) realized interesting numerical simulations.

The soft cores of sandwich panels are usually made of polymers and, thus, have a nonlinear behavior. This physical nonlinearity is present in both elastic and plastic phases and greatly limits the use of most laminate theories cited by the above mentioned literature reviews. This occurs because simplified kinematic are not able to describe the mobility of the continuum during the deformations developed by the soft core, mainly in a crushing situation. Thus, the use of solid finite elements and the improvement of existing material modeling is an important research field.

In order to coherently study the core crushing phenomenon, it is necessary to have in hand a constitutive elastoplastic model capable of considering finite strains. The finite strain theory of plasticity is a very complex subject and has been studied over several decades also resulting in a huge volume of information. A summary of the evolution of this subject in the scientific literature (linked to numerical methods) can be extracted from the work of Zhang and Montán (2019) that specifies two main approaches for modeling plasticity. The first is called multiscale and is important to improve the understanding of the microscopic behavior of materials, but its computational cost makes prohibitive its implementation in current commercial softwares (Coda et al., 2020; Abraham et al., 2002; Buehler et al., 2004). The second plasticity approach (based on continuum mechanics) has a relative low computational cost and, as a consequence, is the most used in the development of engineering softwares.

According to Brepols et al. (2014), when large strains are present, the continuous approach can be divided into two main groups. The first, called hypoelastic, uses the additive decomposition of strain rates and searches for coherent constitutive stress/strain relationships. In this area one may cite works of Atluri (1984), Hughes and Winget (1980), Kojic and Bathe (1987) and Bruhns et al. (1999) that present alternatives to solve objectivity problems that arise from the additive decomposition since the pioneering work of Argyris and Kleiber (1977). Even with these improvements, the application of hypoelastic formulations is limited to materials that present small elastic strains, which is not the case of soft cores.

The second group, called hyperelastic-based plasticity, uses the Kröner-Lee's multiplicative decomposition (Kröner, 1959; Lee, 1969) for theoretical developments. In this strategy an intermediate space appears, where only elastic strains has been developed, which is in general incompatible, usually accompanied by residual stress (Mandel, 1971). When using this kind of formulation it is necessary to operate tensors that are of difficult use and interpretation. There are tensors defined in the intermediate space and tensors dedicated to carry out pushing-over and pulling-back operations among intermediate, current and initial configurations of the studied problem. The elastic portion of this kind of formulation is modeled by hyperelastic potentials and, therefore, admits large elastic strains. Several numerical studies that use this strategy can be cited as, for example, Simo and Ortiz (1985), Weber and Anand (1990), Eterovic and Bathe (1990) and Simo (1992). The works of Zhang and Montán (2019), Brepols et al. (2014), Wallin et al. (2003) and Casey and Naghdi (1981) should be consulted for a good review on different formulations and interesting achievements in finite strain plasticity.

In the present work, with the interest of modeling sandwich panels with soft cores, an alternative and original elastoplastic framework to model continuum mechanics is presented. The proposed formulation cannot be directly classified in hypoelastic or hyperlastic-based frameworks mentioned above, as it is not based on additive or Kröner-Lee strain decompositions. The proposed constitutive model is based on the decomposition of the right Cauchy-Green stretch tensor (Ogden 1984) in its isochoric and volumetric parts, which separates the mechanical behavior of the material in such a way that volumetric strains are only elastic and isochoric strains may develop plastic flow. The differences of the proposed technique and the existent frameworks are: (i) the plastic flow direction is not governed by a plastic surface, (ii) there are two yielding surfaces and (iii) the volumetric change is controlled by a unique potential. It is important to mention that, compared to the Kröner-Lee’s decomposition strategy; the proposed formulation reduces the necessary mathematical operations to configure an elastoplastic procedure in finite strains. Moreover, due to the direct choice of isochoric strain directions, the proposed model avoids material self intersection, called by Zhang and Montán (2019) as “Jacobean issues” when solving a pure shear problem.

The proposed constitutive model is implemented in an in-house geometrically nonlinear computational code that uses the position based Finite Element Method with high order prismatic triangular base elements (Carrazedo et al., 2018;Carrazedo et al., 2020 and Coda and Carrazedo, 2017). High order elements improve the stress representation and allow the use of complete numerical integration (without locking) and prismatic elements adapt very well to plates, shells and sandwich panels geometries. The applications in sandwich panels, in which the finite plasticity occurs in the core, are compared with experimental results present in specialized literature, demonstrating the precision and applicability of the proposed model.

2 NUMERICAL PRELIMINARIES

This section briefly presents the kinematics adopted to write the used positional FEM. Those interested in more details on the development of the adopted finite element may consult Carrazedo et al. (2018), Carrazedo et al. (2020) and Coda and Carrazedo (2017), for instance. Here solids are modeled by prismatic elements with triangular base of any approximation order. Figure 1 shows an element with a cubic approximation at the base and linear in its height. The approximations of initial and current configurations are similar as they are written by shape functions $\overrightarrow{\varphi }$ which has the dimensionless space $\overrightarrow{\xi }=({\xi }_1,{\xi }_2,{\xi }_3)$ as domain, see equation (1):

$f_i^0=x_i={\varphi }_{\alpha }\left(\overrightarrow{\xi }\right)X_{\alpha i}and{f}_i^1=y_i={\varphi }_{\alpha }\left(\overrightarrow{\xi }\right)Y_{\alpha i}$ (1)

in which $f_i^0$ is the initial configuration mapping, $x_i$ represents the initial domain points coordinates, ${\xi }_j$ are the dimensionless coordinates, $X_{\alpha i}$ are initial coordinates of nodes $\alpha$. In the same way, $f_i^1$ is the current configuration mapping, $y_i$ represents the current domain points coordinates and $Y_{\alpha i}$ are current nodal positions.

Figure 1:
Deformation function described by prismatic elements and auxiliary dimensionless space

The deformation function $\overrightarrow{f}$ is represented by the following composition:

$\overrightarrow{f}={\overrightarrow{f}}^1\circ {\left({\overrightarrow{f}}^0\right)}^{-1}$ (2)

and its gradient is given by:

$A=A^1\cdot {\left(A^0\right)}^{-1}$ (3)

with

$A_{ij}^0=\frac{\partial x_i}{\partial {\xi }_j}and{A}_{ij}^1=\frac{\partial y_i}{\partial {\xi }_j}$ (4)

It is important to mention that $A^0$ is known at all integration points (calculated only once) and $A^1$ is known as a trial and enables numerical calculations of all necessary values to assemble the iterative solution of geometrically nonlinear problems, see Coda and Carrazedo (2017), Siqueira and Coda (2016) and Pascon and Coda (2013) for instance.

The weak form of the static equilibrium to be solved is:

${\displaystyle {\int }_{V_0}^{}{S}^c:\delta Ed{V}_0}-{\displaystyle {\int }_{V_0}^{}{\overrightarrow{b}}_0^{}\cdot \delta \overrightarrow{y}d{V}_0}-{\displaystyle {\int }_{{\Gamma }_0}^{}{\overrightarrow{p}}_0\cdot \delta \overrightarrow{y}d{\Gamma }_0}=0$ (5)

in which ${\overrightarrow{b}}_0$ is the body force, ${\overrightarrow{p}}_0$ is the surface force, $\delta \overrightarrow{y}$ is an arbitrary position variation, $S^c$ is the complete second Piola-Kirchhof stress, $E$ is the Green-Lagrange strain, $V_0$ is the initial volume of the analyzed body and ${\Gamma }_0$ is its initial surface.

Remembering that current positions of nodes are the unknowns and using the approximation defined by equation (1), one writes:

$\delta E=\frac{\partial E}{\partial \overrightarrow{Y}}\cdot \delta \overrightarrow{Y}$ (6)

As $\delta \overrightarrow{Y}$ is an arbitrary position variation, the nonlinear system of equations to be solved is:

$\left({\displaystyle {\int }_{V_0}^{}{S}^c:\frac{\partial E}{\partial \overrightarrow{Y}}d{V}_0}-{\displaystyle {\int }_{V_0}^{}\overrightarrow{\varphi }\otimes \overrightarrow{\varphi }d{V}_0\cdot {\overrightarrow{B}}^0}-{\displaystyle {\int }_{{\Gamma }_0}^{}\overrightarrow{\phi }\otimes \overrightarrow{\phi }d{\Gamma }_0\cdot \overrightarrow{P}}\right)=\overrightarrow{0}$ (7)

in which different shape functions are used for domain $\overrightarrow{\varphi }$ and surface $\overrightarrow{\phi }$ approximations as the number of nodes is not the same.

This equation was solved several times by the positional FEM using elastic or elastoplastic constitutive models for small strains, see Carrazedo et al. (2018), Carrazedo et al. (2020) and Coda and Carrazedo (2017) {40-41], for example. In the present work the complete second Piola-Kirchhoff stress tensor $S^c$ will be achieved by an alternative finite strain elastoplastic constitutive model in order to solve sandwich panes with elastoplastic soft cores. It is interesting to advance, see section 4.4.2, that the complete stress will be written in the following form $S^c={(S^{vol})}^e+{(S^1)}^{ep}+{(S^2)}^{ep}$, i.e., the volumetric part of the complete stress will be elastic and the two deviatoric parts of the stress will be elastoplastic.

3 ISOCHORIC AND VOLUMETRIC DECOMPOSITION - FLORY’S DECOMPOSITION

In this section, the decomposition of the right Cauchy-Green stretch tensor in its isochoric and volumetric parts is presented, not only to proceed with the usual definitions of hyperelastic constitutive models, but to support the elastoplastic model originally presented here. It is also shown, through simple operations, a well-known property of nonlinear continuum mechanics for which the isochoric portions of the Green-Lagrange strain are directly associated with the Cauchy deviatoric stresses, and the volumetric portion of the Green-Lagrange strain is associated with the Cauchy hydrostatic stress.

It is important to remember that the applied decomposition is not the Kröner-Lee that splits strains into plastic and elastic parts. The multiplicative decomposition applied on the deformation gradient $A$ (equation (3)) is defined as (Flory, 1961):

$A=\widehat{A}\cdot \overline{A}$ (8)

where

$\widehat{A}=J^{1/3}I\Rightarrow Det\left(\widehat{A}\right)=J$ (9)

$\overline{A}=J^{-1/3}A\Rightarrow Det\left(\overline{A}\right)=1$ (10)

here $\widehat{A}$ is the volumetric part of the deformation gradient $A$ and $\overline{A}$ is its isochoric part.

As the adopted strain is the Green-Lagrange$E=\left(C-I\right)/2$, written as a function of the right Cauchy stretch tensor $C=A^t\cdot A$, it is necessary to write, from equation (10), the isochoric part of the right Cauchy-Green stretch as:

$\overline{C}=J^{-2/3}{A}^t\cdot A$ (11)

In which $Det(\overline{C})\equiv 1$. The volumetric part of $C$ is defined from equation (9) as:

$\widehat{C}={\widehat{A}}^t\cdot \widehat{A}=J^{2/3}I$ (12)

where $Det(\widehat{C})=J^2=Det(C)$.

Multiplying equations (11) and (12) results the necessary multiplicative decomposition of the right Cauchy-Green stretch tensor as:

$C=\widehat{C}\cdot \overline{C}=\overline{C}\cdot \widehat{C}$ (13)

For hyperelastic materials the specific strain energy is written using an elastic potential of the form

$\psi ={\psi }^{vol}(J)+{\psi }^{iso1}({\overline{I}}_1)+{\psi }^{iso2}({\overline{I}}_2),orincompactform\psi ={\psi }^J+{\psi }^1+{\psi }^2$ (14)

in which ${\overline{I}}_1$ and ${\overline{I}}_2$ are, respectively, the first and second invariants of the isochoric part of the Cauchy-Green stretch$\overline{C}$ and $J$ is the third invariant of its volumetric part $\widehat{C}$.

Using the concept of energy conjugate, the second Piola-Kirchhoff stress tensor is written as:

$S=\frac{\partial \psi }{\partial E}=\frac{\partial {\psi }^J}{\partial E}+\frac{\partial {\psi }^1}{\partial E}+\frac{\partial {\psi }^2}{\partial E}orincompactformS=S^J+S^1+S^2$ (15)

In this work, during the elastic phase a constitutive law based on the Rivlin and Saunders (1951) model with a volumetric control similar to the one proposed by Hartmann and Neff (2003) is adopted, as follows:

${\psi }^J=\frac{K}{8n^2}\left(J^{2n}+J^{-2n}-2\right),{\psi }_{}^1=\frac{G}{4}\left({\overline{I}}_1-3\right),{\psi }_{}^2=\frac{G}{4}\left({\overline{I}}_2-3\right)$ (16)

in which $K$is the bulk modulus, $G$ is the shear elastic modulus and $n>0$ is a constant that can help the control of the volumetric stiffness. Other hyperelastic models may be used if desired, respecting equation (14), however all of them should be coincident for small elastic strains.

Using the potentials defined in equation (16) and performing the operations defined in equation (15), for elastic phase one writes:

$S^J=\left\{\frac{K}{4n}\left(J^{2n-1}-J^{-(2n+1)}\right)\right\}{\mathfrak{E}}^J,S^1=\frac{G}{4}{\mathfrak{E}}^{1}and{S}^2=\frac{G}{4}{\mathfrak{E}}^2$ (17)

where the different strain components, that results in stress "directions", are given by:

${\mathfrak{E}}^J=\frac{\partial J}{\partial E}=J{C}^{-1}$ (18)

${\mathfrak{E}}^1=\frac{\partial {\overline{I}}_1}{\partial E}=-\frac{2}{3}{J}^{-2/3}Tr(C)C^{-1}+2J^{-2/3}I$ (19)

${\mathfrak{E}}^2=\frac{\partial {\overline{I}}_2}{\partial E}=2J^{-4/3}\left(-\frac{2}{3}{C}^{-1}{I}_2+\left\{Tr(C)I-C^t\right\}\right)$ (20)

in which $I_2$ is the second invariant of $C$.

In what follows it is shown that the second Piola-Kirchhoff stress components $\left(S^J,S^1,S^2\right)$ defined by equation (17) are related, respectively, to the hydrostatic and deviatoric Cauchy stress components$\left({\sigma }_{}^{vol},{\sigma }^{iso1},{\sigma }^{iso2}\right)$. This is done here applying the well known relation among the second Piola-Kirchhof stress and the Cauchy stress, as follows:

$\sigma =\frac{1}{J}A\cdot S\cdot A^t=\frac{1}{J}A\cdot \left(S^J+S^1+S^2\right)\cdot A^t=\frac{1}{J}A\cdot S^J\cdot A^t+\frac{1}{J}A\cdot S^1\cdot A^t+\frac{1}{J}A\cdot S^2\cdot A^t$ (21)

To make this operation simple one starts rewriting equation (16) in a general form:

${\psi }^J={\psi }^J(J),{\psi }_{}^1={\psi }_{}^1({\overline{I}}_1),{\psi }_{}^2={\psi }_{}^2({\overline{I}}_2)$ (22)

Thus, for each component one writes:

a) Volumetric

The candidate to the volumetric component of the second Piola-Kirchhoff stress is written as:

$S_{}^J=\frac{\partial {\psi }_{}^J}{\partial J}\frac{\partial J}{\partial E}=\alpha {\mathfrak{E}}^J$ (23)

In which $\alpha$ is a scalar, see equation (17) and ${\mathfrak{E}}^J$ is given by equation (18).

Applying Eq. (21) over the second Piola-Kirchhoff stress of Eq. (23) results:

${\sigma }_{}^{vol}=\frac{1}{J}A\cdot S_{}^J\cdot A^t=\frac{1}{J}A\cdot \alpha J{C}^{-1}\cdot A^t=\alpha I={\sigma }^{hid}$ (24)

which means that ${\mathfrak{E}}^J$ is the Lagrangian strain direction corresponding to the hydrostatic stress direction in the Cauchy space. Thus, $S^J$ is hydrostatic in Lagrangian sense.

b) First isochoric

The candidate to the first isochoric component of the second Piola-Kirchhoff stress is written as:

$S_{}^1=\frac{\partial {\psi }_{}^1}{\partial {\overline{I}}_1}\frac{\partial {\overline{I}}_1}{\partial E}=\beta {\mathfrak{E}}^1$ (25)

where $\beta$ is a scalar and ${\mathfrak{E}}^1$ is a symmetric tensor of the same order of the Green strain, see equation (19). Considering Eq. (19) and applying the transformation (21) over Eq. (25), results:

${\sigma }^1=\beta \left\{2J^{-5/3}\left(A\cdot A^t-\frac{tr\left(A^t\cdot A\right)}{3}I\right)\right\}$ (26)

And, as

$tr\left(A^t\cdot A\right)=A:A^t=A^t:A=tr(A\cdot A^t)$ (27)

one achieves:

${\sigma }^{iso1}=\beta \left\{2J^{-5/3}\left(A\cdot A^t-\frac{tr\left(A\cdot A^t\right)}{3}I\right)\right\}={\sigma }^{dev}$ (28)

which means that ${\mathfrak{E}}^1$ is the first isochoric Lagrangian strain direction and $S^1$ is the first isochoric second Piola-Kirchhoff stress.

c) Second isochoric

The candidate to second isochoric component of the second Piola-Kirchhoff stress is written as:

$S^2=\frac{\partial {\psi }^2}{\partial {\overline{I}}_2}\frac{\partial {\overline{I}}_2}{\partial E}=\gamma {\mathfrak{E}}^2$ (29)

with $\gamma$ and ${\mathfrak{E}}^2$being respectively a scalar and a symmetric tensor of the same order as the Green-Lagrange strain tensor, see equation (20). Using Eq. (21) over Eq. (29) and considering Eq. (20), one writes:

${\sigma }^{iso2}=\gamma 2J^{-7/3}\left[\left(tr(C)\left(A\cdot A^t\right)-\left(A\cdot A^t\right)\cdot \left(A\cdot A^t\right)\right)-\left(\frac{2}{3}{I}_2\right)I\right]$ (30)

Using equation (27) results:

$tr(A\cdot A^t)tr(A\cdot A^t)-tr((A\cdot A^t)\cdot (A\cdot A^t))=2I_2$ (31)

From Eq. (31) one calculates $Tr\left({\sigma }^{iso2}\right)$ as

$Tr\left({\sigma }^{iso2}\right)=2\gamma J^{-7/3}\left(2I_2-2I_2\right)=0$ (32)

thus:

${\sigma }^{iso2}={\sigma }^{dev},$ (33)

meaning that ${\mathfrak{E}}^2$ is the second isochoric Lagrangian strain direction and $S^2$ is the second isochoric component of the second Piola-Kirchhoff stress.

Thus, from the compact form of equation (15) one understands that an isotropic ductile material enters in plastic flow when the deviatoric elastic stressess$S^1$ or $S^2$ (without volumetric components) reach a critical value.

4 PROPOSED CONSTITUTIVE MODEL

The finite strain elastoplastic model proposed here fulfill three hypothesis stated as follows: (i) The volumetric changes are exclusively elastic, (ii) At any instant the total plastic strain tensor is deviatoric and (iii) The multiplicative decomposition of the Cauchy-Green stretch tensor (in volumetric and isochoric parts) guaranties independent evolutions of hydrostatic and deviatoric stresses even when plastic flow occurs.

4.1 Elastic limit - before occurrence of plastic flow

The usual von- Mises condition that guaranties the elastic behavior of an isotropic ductile material is given by the following expression:

$J_2<\frac{{\sigma }_y^2}{3}or\frac{1}{2}{\sigma }^{dev}:{\sigma }^{dev}<\frac{{\sigma }_y^2}{3}$ (34)

in which ${\sigma }_y$ is the reference yielding stress and ${\sigma }^{dev}$ is the elastic deviatoric part of the Cauchy stress tensor with $J_2$ being its second invariant. In order to develop calculations (choosing the shear stress as a reference) one defines ${\sigma }_y=2\overline{\tau }$, resulting:

$\sqrt{\frac{3}{2}{\sigma }^{dev}:{\sigma }^{dev}}-2\overline{\tau }<0$ (35)

As shown in section 3, before yielding takes place, the deviatoric elastic stress can be represented by the Lagrangian components ($S^1$ and $S^2$) established in equations (25) and (29) or (17). As there are two deviatoric (isochoric) strain directions (that will be used to define plastic flow direction, see item 4.2.1.a) in the proposed model it is necessary to split the von-Mises criterion to consider each Lagrangiandeviatoric stress components. It is summarized by imposing the following Lagrangian conditions:

$\sqrt{\frac{3}{2}{S}^1:S^1}-{\overline{\tau }}_1<0and\sqrt{\frac{3}{2}{S}^2:S^2}-{\overline{\tau }}_2<0$ (36)

in which ${\overline{\tau }}^1={\overline{\tau }}^2=\overline{\tau }$ is assumed, following the elastic relation (17). Other way to write (36) is:

$\frac{3}{2}{S}^1:S^1-{\overline{\tau }}^2<0and\frac{3}{2}{S}^2:S^2-{\overline{\tau }}^2<0$ (37)

Using equation (17) one writes equation (37) as a function of the deviatoric part of the hyperelastic constitutive model, i.e.:

$\frac{3}{2}\frac{G^2}{16}{\mathfrak{E}}^1:{\mathfrak{E}}^1-{\overline{\tau }}^2<0and\frac{3}{2}\frac{G^2}{16}{\mathfrak{E}}^2:{\mathfrak{E}}^2-{\overline{\tau }}^2<0$ (38)

in which the isochoric “directions” ${\mathfrak{E}}^1$ and ${\mathfrak{E}}^2$ are present. In order to simplify operational inferences, one rewrites equation (38) in an unified form, with $\mathfrak{E}$ representing both ${\mathfrak{E}}^1$ or ${\mathfrak{E}}^2$:

$\frac{3}{2}\frac{G^2}{16}\mathfrak{E}:\mathfrak{E}-{\overline{\tau }}^2<0$ (39)

The initial elastic limiting surfaces are defined imposing equality in equation (38), or schematically in equation (39). Here appears a difference of the proposed formulation and others, i.e., there are two yielding surfaces, not one.

In order to verify the validity of substitutingequation (35) by equation (36), one proceeds with a small strain uniaxial test (for example, following direction $x_1$). For small strains and the proposed test one has,

$J={\lambda }_1{\lambda }_2{\lambda }_3\cong 1$ (40)

in which ${\lambda }_i$represents stretch in direction $i$. Remembering that the material is assumed isotropic, for this uniaxial test one writes ${\lambda }_2={\lambda }_3=\lambda$or even ${\lambda }^2=1/{\lambda }_1$

Making some algebraic simple steps one achieves from equations (19) and (20):

${\mathfrak{E}}^1=\frac{2}{3}\left[\begin{array}{ccc}2(1-{\lambda }_1^{-3})& 0& 0\\ 0& 1-{\lambda }_1^3& 0\\ 0& 0& 1-{\lambda }_1^3\end{array}\right]$ (41)

and

${\mathfrak{E}}^2=\frac{2}{3{\lambda }_1^{}}\left[\begin{array}{ccc}2(1-{\lambda }_1^{-3})& 0& 0\\ 0& 1-{\lambda }_1^3& 0\\ 0& 0& 1-{\lambda }_1^3\end{array}\right]$ (42)

Using equation (17), the following particular relation is achieved:

$S^1={\lambda }_1{S}^2$ (43)

Using equation (21) and remembering that for small strains one has${\lambda }_1\cong 1$results:

$S_1=S_2={\sigma }_1^{dev}={\sigma }_2^{dev}={\sigma }^{dev}/2$ (44)

Thus one achieves

$\sqrt{S_2:S_2}=\sqrt{S_1:S_1}=\frac{1}{2}\sqrt{{\sigma }^{dev}:{\sigma }^{dev}}$ (45)

verifying the proposition for small strains. For large strains a simple calibration of $\overline{\tau }$ can be made.

4.2 Proposed elastoplasticity

Equation (39) needs to be improved to consider previous plastic flow and the possibility of plastic evolution. Thus, some classic definitions are adapted here to the elastoplastic model to be proposed in this study.

4.2.1 Definitions

4.2.1.a - Plastic evolution

The local plastic evolution is established from the observation that the general strain “direction” $\mathfrak{E}$ is isochoric, i.e.:

$\Delta {\mathfrak{E}}^p=\Delta \lambda \frac{\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}or\Delta S^p=\frac{G}{4}\Delta \lambda \frac{\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}$ (46)

where $\mathfrak{E}/\sqrt{\mathfrak{E}:\mathfrak{E}}$ represents the unitary deviatoric directions of plastic strain evolution $\Delta {\mathfrak{E}}^p$, and $\Delta \lambda$ is the variation (or evolution) of the plastic multiplier. As the deviatoric elastic directions, see equation (17), are the same defined for plastic flow, one defines the evolution of an internal variable $\Delta S^p$called here, by similarity with equation (17), the plastic stress change. The statement of equation (46) is useful in the definition of the tangent algorithm constitutive tensor, but has local validity, that is, cannot be applied for the plastic stress (above defined) cumulative process.

In order to be clear, an accumulation of the type $S^p=S^p+\Delta S^p$ cannot be used because, for large strains, a past value of $\Delta S^p$ stop being isochoric (at present time) and, therefore, the use equation (21) with the up-to-date deformation gradient ($A$) implies loss objectivity for past values.

As it is not possible to make a simple sum to state the accumulated plastic stress, here one writes:

$S^p=\frac{G}{4}{\lambda }^{\zeta }\frac{\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}$ (47)

and the evolution is imposed on the scalar plastic strain ${\lambda }^{\zeta }$ as:

${\lambda }^{\zeta }={\left({\lambda }^{\zeta }\right)}^{ac}+\zeta \Delta \lambda$ (48)

in which $ac$ represents cumulated and $\zeta$ is the sign of plastic evolution, to be defined after the introduction of kinematic hardening, see equation (50). As $S^p$ is proportional to $\mathfrak{E}$, hypothesis (ii) and hypothesis (iii) are fulfilled. After the determination of the plastic multiplier $\Delta \lambda$, a simple evolution of the plastic stress tensor is given by equation (47) in association with equation (48). The reader understands that in this item an important difference of the proposed formulation and the classical frameworks appear, i.e., the plastic flow direction did not depend on a plastic surface, but follows phenomenological isochoric kinematic directions. This choice prevents material self intersection, i.e., guaranties that $J>0$ as the volumetric changes are controlled by a consistent hyperelastic potential, see equation (16).

4.2.1.b - Kinematic hardening

For convenience, the kinematic hardening ($H^c$) is considered constant at iterations of the numerical solution process, but depends upon the scalar plastic strain${\lambda }^{\zeta }$ for each load level. Defining the ratio between the kinematic hardening and the shear elastic modulus as:

$\beta ({\lambda }^{\zeta })=H^c({\lambda }^{\zeta })/G$ (49)

the internal variable of the kinematic hardening becomes:

$q=q^{ac}+\beta ({\lambda }^{\zeta })\zeta \Delta \lambda$ (50)

and the signal of the scalar plastic strain evolution becomes $\zeta =Sign\left(\sqrt{\mathfrak{E}:\mathfrak{E}}-{\lambda }^{\zeta }-q\right)$.

The internal variable $q$ is multiplied by the isochoric direction and the elastic modulus to generate the back stress tensor as

$\chi =q\frac{G}{4}\frac{\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}$ (51)

The back stress tensor is present in equation (54) and simplified by the double contraction performed in equation (55).

4.2.1.c - Isotropic hardening

The internal variable $\kappa$associated with isotropic hardening is calculated by a cumulative process of the plastic multiplier $\Delta \lambda$, see equation (47). The isotropic hardening ($H^i$) is considered constant during iterations, but depends upon the cumulated plastic multiplier. Similarly to the kinematic hardening, one defines:

$\eta (\lambda )=H^i(\lambda )/G$ (52)

The plastic multiplier $\lambda$is cumulated by $\lambda ={\lambda }^{ac}+\Delta \lambda$ and the evolution of the internal variable $\kappa$ is given by:

$\kappa ={\kappa }^{ac}+\eta (\lambda )\Delta \lambda$ (53)

4.2.2 Preparation of the objective function

Adopting the symbol $S^t$ for both isochoric Piola-Kirchhoff stresses $S^1$ or $S^2$ of equation (17) and considering previous definitions, the stress to be introduced in the yielding function (objective function) is defined by:

$\begin{array}{l}S=S^t-S^p-\chi \\ or\\ S=\frac{G}{4}\mathfrak{E}-\frac{G}{4}{\lambda }^{\zeta }\frac{\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}-\frac{G}{4}q\frac{\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}=\frac{G}{4}\left(1-\frac{{\lambda }^{\zeta }}{\sqrt{\mathfrak{E}:\mathfrak{E}}}-\frac{q}{\sqrt{\mathfrak{E}:\mathfrak{E}}}\right)\mathfrak{E}\end{array}$ (54)

Thus, one calculates:

$S:S=\frac{G^2}{16}{\left(1-\frac{{\lambda }^{\zeta }}{\sqrt{\mathfrak{E}:\mathfrak{E}}}-\frac{q}{\sqrt{\mathfrak{E}:\mathfrak{E}}}\right)}^2\mathfrak{E}:\mathfrak{E}$ (55)

It is important to remember that

$S^{ep}=S^t-S^p=(G/4)(1-{\lambda }^{\zeta })\frac{\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}$ (56)

represents the elastoplastic deviatoric parts (${(S^1)}^{ep}$ and ${(S^2)}^{ep}$) of the complete stress $S^c={(S^{vol})}^e+{(S^1)}^{ep}+{(S^2)}^{ep}$ defined after equations (5) and (7); and that ${(S^{vol})}^e$ is the elastic volumetric part given in the first of equations (17).

Taking into account translation (kinematic hardening $q$) and the change of size of the plastic surface (isotropic hardening $\kappa$), equation (39) is upgraded to the current yielding function or objective function, as:

$f=\frac{3}{2}\frac{G^2}{16}{\left(1-\frac{{\lambda }^{\zeta }}{\sqrt{\mathfrak{E}:\mathfrak{E}}}-\frac{q}{\sqrt{\mathfrak{E}:\mathfrak{E}}}\right)}^2\mathfrak{E}:\mathfrak{E}-{\left(\sqrt{\frac{3}{2}}\frac{G}{4}\kappa +\overline{\tau }\right)}^2<0$ (57)

The objective function $f$ given in equation (57) has the usual meaning, i.e, for $f<0$ the regime is elastic and there is no plastic evolution. As mentioned before, see equations (47) and (48), the plastic evolution follows isochoric direction and any volume change is controlled by the volumetric term of equation (17), thus elastic (hypothesis (i)).

4.2.3 Calculating the plastic multiplier $\Delta \lambda$

Plastic evolution occurs when a trial objective function violates the plastic criterion (51), that is:

$f^{tr}=\frac{3}{2}\frac{G^2}{16}{\left(1-\frac{{\lambda }^{ac\zeta }}{\sqrt{{\mathfrak{E}}^{tr}:{\mathfrak{E}}^{tr}}}-\frac{q^{ac}}{\sqrt{{\mathfrak{E}}^{tr}:{\mathfrak{E}}^{tr}}}\right)}^2{\mathfrak{E}}^{tr}:{\mathfrak{E}}^{tr}-{\left(\sqrt{\frac{3}{2}}\frac{G}{4}{\kappa }^{ac}+\overline{\tau }\right)}^2\ge 0$ (58)

in which the superscript $tr$ represents elastic trial without plastic evolution. The usual procedures for plastic analysis (permanence of stress levels on the plastic surface) are summarized here as simple operational, i.e., imposing the nullity of equation (58) allowing the evolution of the plastic multiplier at the same time, as follows:

$\frac{3}{2}\frac{G^2}{16}{\left(1-\frac{({\lambda }^{ac\zeta }+q^{ac})}{\sqrt{{\mathfrak{E}}^{tr}:{\mathfrak{E}}^{tr}}}-(1+\beta )\zeta \frac{\Delta \lambda }{\sqrt{{\mathfrak{E}}^{tr}:{\mathfrak{E}}^{tr}}}\right)}^2{\mathfrak{E}}^{tr}:{\mathfrak{E}}^{tr}={\left(\sqrt{\frac{3}{2}}\frac{G}{4}({\kappa }^{ac}+\eta \Delta \lambda )+\overline{\tau }\right)}^2$ (59)

As hardenings are considered constant at iterations, equation (59) is easily solved and the plastic multiplier is the smaller positive value of:

$\Delta \lambda =\left[\left(\sqrt{{\mathfrak{E}}^{tr}:{\mathfrak{E}}^{tr}}-({\lambda }^{ac\zeta }+q^{ac})-{\kappa }^{ac}\right)-\overline{\tau }/\left(\sqrt{\frac{3}{2}}\frac{G}{4}\right)\right]/\left(\zeta (1+\beta )+\eta \right)$ (60)

$\Delta \lambda =\left[\left(\sqrt{{\mathfrak{E}}^{tr}:{\mathfrak{E}}^{tr}}-({\lambda }^{ac\zeta }+q^{ac})+{\kappa }^{ac}\right)+\overline{\tau }/\left(\sqrt{\frac{3}{2}}\frac{G}{4}\right)\right]/\left(\zeta (1+\beta )-\eta \right)$ (61)

For cyclic applications, when plastic stress ($S^p$) exists and $\sqrt{ℭ:ℭ}$ is small, equation (47) may become ill conditioned. Thus, in this situation it is necessary to use an alternative objective function to calculate the plastic evolution, see item 4.3.

4.2.4 Isochoric elastoplastic tangent constitutive tensor

Equation (46) is valid at an instant (not used to assemble evolution), thus it is possible to write equation (56) as:

$S^{ep}=S^t-{\left(S^p\right)}^{ac}-\Delta S^p$ (62)

Observing that ${\left(S^p\right)}^{ac}$ does not vary and $\Delta S^p$ is a finite value, one writes:

$\delta S^{ep}=\frac{\partial S^t}{\partial E}:\delta E-\frac{\partial \Delta S^p}{\partial E}:\delta E=\frac{G}{4}\frac{\partial \mathfrak{E}}{\partial E}:\delta E-\frac{G}{4}\Delta \lambda \frac{\partial \mathfrak{E}}{\partial E}:\delta E$ (63)

or

$\delta S^{ep}=\frac{G}{4}(1-\Delta \lambda )\frac{{\partial }^{2}I}{\partial E\otimes \partial E}:\delta E$ (64)

Resulting the proposed model algorithm tangent constitutive tensor as:

$ℭ=\frac{G}{4}(1-\Delta \lambda )\frac{{\partial }^{2}I}{\partial E\otimes \partial E}$ (65)

It is important to mention that the volumetric part of the complete constitutive tensor is elastic.

4.3 Small strain range - plastic evolution

When small strains occur with $S^p\ne 0$ (cyclic loads) a small strain procedure is necessary due to an ill conditioning of equation (47), i.e., when the quantity $\mathfrak{E}/\sqrt{\mathfrak{E}:\mathfrak{E}}$ passes through its nullity, it loses the ability of representing the total plastic stress value $S^p$, given in equation (47), but when $S^p$ is actually null the problem does not exist.

Here it is also done using Flory’s decomposition, thus the isochoric and volumetric strains follow the same separation used for large strains. Using the unified notation for $S^1$ and $S^2$ the total quantity to be verified in the yielding criterion is:

$S=S^t-S^p-\chi$ (66)

in which $\chi$ is the back-stress tensor. The pure elastic stress $S^t=(G/4)\mathfrak{E}$ and the accumulated plastic stress $S^p$ come from a previous existence of the large strain procedure, so are isochoric. From these tensors one writes the plastic direction as:

$\vartheta =\frac{S^t-S^p}{\sqrt{(S^t-S^p):(S^t-S^p)}}=\frac{S^t-S^p}{\Vert S^t-S^p\Vert }$ (67)

resulting:

$\Delta S^p=\frac{G}{4}\Delta \lambda \vartheta and\Delta \chi ={\overline{H}}^c\frac{G}{4}\Delta \lambda \vartheta$ (68)

Thus, the objective function becomes:

$f^{tr}={\left(\sqrt{\frac{3}{2}\left(S^{tr}:S^{tr}\right)}\right)}^2-{\left(\sqrt{\frac{3}{2}}{\overline{\kappa }}^{ac}+\overline{\tau }\right)}^2\le 0$ (69)

Evolutions are defined by:

$\{\begin{array}{l}\overline{\kappa }={\overline{\kappa }}^{ac}+{\overline{H}}^i\Delta \lambda \\ \overline{\lambda }={\overline{\lambda }}^{ac}+\Delta \lambda \\ \chi =\chi +\Delta \chi \\ S^p=S^p+\Delta S^p\end{array}$ (70)

in which $\Delta \lambda$ is the evolution of plastic multiplier. When $f^{tr}\ge 0$ one imposes equality in equation (69) together with plastic multiplier evolution, resulting:

${\left(\sqrt{\frac{3}{2}\left((S^e-S^p-\chi -\vartheta \Delta {\lambda }^{*}-{\overline{H}}^c\vartheta \Delta {\lambda }^{*}):(S^e-S^p-\chi -\vartheta \Delta {\lambda }^{*}-{\overline{H}}^c\vartheta \Delta {\lambda }^{*})\right)}\right)}^2={\left[\left(\frac{G}{4}\sqrt{\frac{3}{2}}{\kappa }^{ac}+\overline{\tau }\right)+\sqrt{\frac{3}{2}}{\overline{H}}^i\Delta {\lambda }^{*}\right]}^2$ (71)

in which one uses $\Delta {\lambda }^{*}=\Delta \lambda (G/4)$, ${\overline{H}}^c=\beta ({\lambda }^{\zeta })=H^c({\lambda }^{\zeta })/G$ and ${\overline{H}}^i=\eta (\lambda )=H^i(\lambda )/G$ to write it short. After some algebraic manipulations equation (71) turns into:

$f^{tr}+\left(\frac{3}{2}{\varpi }^2-\frac{3}{2}{\left({\overline{H}}^i\right)}^2\right){\left(\Delta {\lambda }^{*}\right)}^2-2\left(\frac{3}{2}\varpi \left(S^t:\vartheta \right)+(\frac{G}{4}\sqrt{\frac{3}{2}}{\kappa }_{ac}+\overline{\tau })\sqrt{\frac{3}{2}}{\overline{H}}^i\right)\Delta {\lambda }^{*}=0$ (72)

in which $\varpi =1+{\overline{H}}^c$. As hardening is considered constant during iterations, one solves equation (72) as:

$\Delta {\lambda }^{*}=\frac{-b\pm \sqrt{b^2-4a{f}^t}}{2a}or\Delta \lambda =\frac{-b\pm \sqrt{b^2-4a{f}^{tr}}}{2a}\frac{4}{G}$ (73)

with

$a=\frac{3}{2}\left({\varpi }^2-{\left({\overline{H}}^i\right)}^2\right)$ (74)

and

$b=-2\left(\frac{3}{2}\varpi \left(S^t:\vartheta \right)+(\sqrt{\frac{3}{2}}{\overline{\kappa }}_{ac}+\overline{\tau })\sqrt{\frac{3}{2}}{\overline{H}}^i\right)$ (75)

being the smaller positive value of $\Delta \lambda$ the searched solution. The same algorithm constitutive elastoplastic tensor of finite strains, given by equation (65), is used for small strains.

4.4 Transition between models

Small and large strain variables are the same, but for small strains the tensor values $S^p$ and $\chi$ should be stored, while for large strains the correspondent scalar values ${\lambda }^{\zeta }$ and $q$ are stored. When running the large strain part of the model the tensor variables are directly known by equations (47) and (51). However, when running the small strain part of the model, one should calculate the scalar values by multiplying equations (47) and (51) by the isochoric strain tensor $\mathfrak{E}$, as follows

$\mathfrak{E}:S^p=\frac{G}{4}{\lambda }^{\zeta }\frac{\mathfrak{E}:\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}and\mathfrak{E}:\chi =\frac{G}{4}q\frac{\mathfrak{E}:\mathfrak{E}}{\sqrt{\mathfrak{E}:\mathfrak{E}}}$ (76)

resulting

${\lambda }^{\zeta }=\frac{4}{G}\frac{S^p:\mathfrak{E}}{\mathfrak{E}:\mathfrak{E}}andq=\frac{4}{G}\frac{\chi :\mathfrak{E}}{\mathfrak{E}:\mathfrak{E}}$ (77)

Finally the transition is made from large strains to small strains when $\sqrt{\mathfrak{E}:\mathfrak{E}}<Tol$ and from small strains to large strains when $\sqrt{\mathfrak{E}:\mathfrak{E}}>Tol$. It is worth mentioning that the process always start using the large strain procedure, because in this situation $S^p=0$.

5 NUMERICAL EXAMPLES

Three examples are presented in this section. The first one intends to show that the proposed model is able to represent finite strain plasticity including cyclic loads, isotropic and kinematic hardenings. The next two examples are used to model sandwich panels with core crushing and use experimental results to validate the proposed model in real situations.

5.1 Cyclic axial test

This example consists in a unitary cub allowed to slide at faces $x=0$, $y=0$ and $z=0$, see Figure 2a. The cub is stretched (tension and compression) in $z$ direction. As the specimen develops large strains - keeping rectangular shape - it is interesting to compare the achieved second Piola-Kirchhoff stress (related to the initial area) with the Cauchy stress values (related to the current area). The difference between them indicates the influence of large strains in the measured physical quantities.

The adopted dimensionless physical properties are: Bulk modulus $K=80$, shear elastic modulus $G_1=G_2=24$, yielding stress ${\overline{\tau }}_1={\overline{\tau }}_2=0.4$ and different isochoric and kinematic hardenings for each case. A simple discretization with only two linear elements is employed, see Figure 2b.

Figure 2
Sliding conditions and simple discretization

Figure 3 presents the behavior of the second Piola-Kirchhoff stress ($S_{33}$) and the Cauchy stress (${\sigma }_{33}$) for a moderate stretch interval $0.9\le \lambda \le 1.1$ considering $H_1^i=H_2^i=0.5$. As one can see in Figures 3, for example, the model reproduces the one dimensional problem, as the beginning of the plastic flow is at ${\overline{\sigma }}_y=0.8$. For this small range of strain the difference between Piola-Kirchhoff and Cauchy stresses is small and the plastic surface grows as expected, i.e. linearly.

Figure 3
Isotropic hardening for moderate stretching (tension) interval.

Figure 4
Kinematic hardening for moderate stretching

Figure 4 presents the Piola-Kirchhoff ($S_{33}$) and Cauchy (${\sigma }_{33}$) stresses for the same moderate stretching $0.9\le \lambda \le 1.1$, but considering kinematic hardening $H_1^c=H_2^c=0.5$. The closure of the stress, see Figure 4, confirms the reproduction of the Mansing effect when kinematic hardening is adopted, i.e, the energy dissipation is due to plastic flow only.

Figure 5 presents the specimen behavior for the case $H_1^c=H_2^c=0.5$, with a larger stretch interval $0.4\le \lambda \le 1.6$. As it is expected, the hardening behavior is not linear when large strains appear it occurs because the volumetric part of the model is highly activated. Particularly, in compression, even with a small hardening, the stress values should present a fast grow to avoid material self intersection. This excellent answer is only possible due to complete separation of strains into isochoric and volumetric parts. One can observe that, even for this level of strain, the stress cycle closes perfectly. Moreover, due to the large change of area, the Piola-Kirchhoff and Cauchy stresses are very different. Table 1 presents important values to draw Figure 5 in order to make possible the reproduction of the proposed model for this cyclic load and large strain range.

Figure 5
Kinematic hardening for large strain - stress cycle

Table 1
Cauchy and Piola-Kirchhoff stresses for large cyclic strain reange

$\lambda -1$	Cauchy	Piola-Kirchhoff
0	0	0
0.006	0.39092	0.38922
0.012	0.77829	0.77156
0.018	0.82088	0.80914
0.3	1.27354	0.98491
0.6	1.5682	0.98665
0.3	-0.48094	-0.36921
-0.3	-1.01703	-1.44679
-0.6	-2.14791	-5.24741
-0.3	-0.10103	-0.14427
0	0.78114	0.78321
0.3	1.27344	0.98482

Figure 6
Different hardening values for monotonic stretching for finite strain.

Figure 6 presents the monotonic stress behavior for a large stretch range $0\le \lambda \le 3.75$. Three different values of hardening $H_1^{}=H_2^{}=0.15$, $H_1^{}=H_2^{}=0.5$ and $H_1^{}=H_2^{}=1.0$ are adopted, remembering that for the monotonic behavior it is irrelevant which kind of hardening is adopted.

One can see in Figure 6 that, when tension is applied, the Cauchy stress also grows (in a smaller rate than in compression) to avoid self intersection. Table 2 presents important values to draw Figure 6 in order to make possible the reproduction of the proposed model.

Table 2
Cauchy and Piola-Kirchhoff stresses for large strain tension

$\lambda -1$	Cauchy stress			Piola-kirchhoff stress
	$H^i=0.5$	$H^i=0.15$	$H^i=1.0$	$H^i=0.5$	$H^i=0.15$	$H^i=1.0$
0	0	0	0	0	0	0
0.006	0.39092	0.39092	0.39092	0.38922	0.38922	0.38922
0.012	0.77829	0.77829	0.77829	0.77156	0.77156	0.77156
0.018	0.82833	0.8178	0.84232	0.81652	0.80609	0.83035
0.024	0.84941	0.82744	0.87855	0.83246	0.81085	0.86113
0.03	0.8702	0.83685	0.91444	0.84794	0.81533	0.89121
0.3	1.63109	1.12579	2.29755	1.26331	0.8701	1.78451
0.6	2.27389	1.27579	3.58186	1.43487	0.80167	2.27286
0.9	2.90743	1.39874	4.86043	1.54919	0.74051	2.61167
1.2	3.58968	1.5569	6.22157	1.65667	0.71234	2.90452
1.5	4.35145	1.75256	7.68687	1.7731	0.70624	3.17861
1.8	5.18723	1.98433	9.24748	1.89411	0.71469	3.43792
2.1	6.10922	2.24969	10.9318	2.02339	0.7327	3.69959
2.4	7.08672	2.54649	12.6999	2.14939	0.75718	3.95097
2.7	8.13219	2.9	14.5544	2.27724	0.78659	4.19723

In order to show the extension of the applied strain, Figure 7 presents snapshots with Cauchy stress values for cases $H_1^{}=H_2^{}=0.5$ within $0.4\le \lambda \le 1.6$ and $H_1^{}=H_2^{}=0.15$ within $0\le \lambda \le 3.75$. In order to identify the initial configuration, Figure 7a shows a position corresponding to $\lambda =1.006$.

Figure 7
Snapshots including Cauchy stress ${\sigma }_{33}$ for 2 cases

5.2 Three-point bending test:

In this example, the three-point bending test performed by Salami et al. (2015) is numerically simulated. It is a sandwich beam with a length of $300mm$ and a width of $50mm$, the distance between supports is $250mm$. As shown in Figure 8a, the polyurethane core has a height of $30mm$ and the spring steel faces CK75 have a height of $1mm$.

The simulation is carried out using a mesh with $3200$ nodes and $336$ prismatic finite elements with cubic approximation (length and height of the beam) and linear approximation along its width. Taking advantage of symmetry, this discretization is used for half the problem, as shown in Figure 8b. Considering that the contact regions (actuator / upper face and supports / lower face) do not suffer degeneration and that in the experimental response there is no softening, a simple control of force is performed. A load increment of $\Delta F=25.0N$ is equally applied in 21 steps for half structure. The load is distributed over an area of $7mmx50mm$ as depicted in Figure 8a. Face-sheets are made of spring steel CK75 and its mechanical properties were experimentally evaluated by Salami et al. (2015), resulting in $K=74.074GPa$ and $G=95.238GPa$. Although no plastic strain takes place in face-sheets for the experiment force range and no information are given by Salami et al. (2015), it is assumed an isotropic hardening of $H_1^i{}^1=H_2^i=0.95GPa$ and an elastic limit of ${\overline{\tau }}_1={\overline{\tau }}_2=\mathrm{386,5}MPa$.

Figure 8
Three point bending experiment and used discretization

The mechanical behavior of the adopted polyurethane is extracted directly from the experimental graph (compression stress versus strain) achieved by Salami et al. (2015). The numerical results are compared to the experimental one in Figure 9 using $K=0.096GPa$, $G=3.10x{10}^{-3}GPa$ and $H_1^c=H_1^c=1.0x{10}^{-4}GPa$. As one can see the proposed model naturally reproduces the stiffness of the core material at the end of the curve.

Figure 9
Polyurethane behavior in compression, experimental values due to Salami et al. (2015).

To build the numerical curve, a one-dimensional test is performed using a cube of unitary dimensions, the same used in example 5.1, see Figure 2. It is important to note that the small delay in material response (presented by Salami et al., 2015) is not included in Figure 9.

Figure 10 shows the applied force (half beam) as a function of the vertical displacement (at the load application line $x=0$) for different situations: (a) Experimental, (b) numerical model with z constrained displacement (plane strain situation), (c) numerical model with free z displacements and (d) same as case (c) with $7mm$ caps at the free ends of the beam. As can be seen, the results are very close for all cases, revealing: (i) that the developed constitutive model adequately represents the plastic evolution of polyurethane and its combination with face-sheets, (ii) the stiffness loss of the sandwich beam is associated to the core elastoplastic evolution and (iii) the introduction of the caps at the free end of the beam, despite eliminating the localized yielding, does not significantly interfere in the overall behavior of the beam. Regarding cases (c) and (d), comparing Figures 11 and 12, one may see de difference in plastic strain near the end of the beam at the final load.

Figure 10
Force versus displacement at the central sandwich beam

Figure 11 shows some snapshots (associated with load levels) of the tested beam (case (c)), presenting the polyurethane plastic strain values ${\lambda }^{\zeta }$. The same is done in Figure 12 for case (d). As one can see, differences in plastic strains appears at the final experimental load level $F=500N$ and are localized at the low corner of the right side of Figures 11 and 12.

Figure 11
Plastic strain level ${\lambda }^{\zeta }$for polyurethane core, case (c)

Figure 12
Plastic strain level ${\lambda }^{\zeta }$for polyurethane core, case (d)

Figure 13 presents the normal stress ${\sigma }_y$ inside the core for case (d) and applied force $F=425N$. As one can see a stress bulb is present under the region of applied load, as well as near the support.

Figure 13
Normal stress ($y$ direction) inside the core, case (d) - $GPa$

There is no yielding of face-sheets for the experimental load interval. In order to further analyze the proposed experimental test, Figure 14a shows the plastic strain ${\lambda }^{\zeta }$ of face-sheets (zoom close to the middle of the beam) for a applied load of $F=600N$. Figure 14b shows the total displacement of the core following $z$ direction for $F=600N$. Case (c) is considered for both Figures 14a and 14b.

Figure 14
Plastic strain ${\lambda }^{\zeta }$for face-sheets and displacement in z direction for $F=600N$

The adopted tolerance in positions for this example is $1.0x{10}^{-5}$.

5.3 Core crushing experiment:

In this example, the experimental study carried out by Koissin et al. (2004) is simulated by the computational code developed here using the proposed elastoplastic model. It is the coupling of a single face-sheet $2.4mm$ thick by $50mm$ wide with a core $50mm$ thick by $50mm$ wide. The length of the specimen is $280mm$ and it is resting on a smooth sliding surface. A cylindrical mechanical actuator with$R=12.5mm$ is used to impose the vertical displacement of the face-sheet as shown in Figure 15.

Figure 15
Schematic representation of half specimen proposed by Koissin et al. (2004)

The same material used by Koissin et al. (2004)(GRFP) is adopted to build up the face-sheet with the experimental elastic modulus $E=16.5GPa$ and Poisson ratio $\nu =0.25$. From these values results $K=2.75GPa$ and $G=6.6GPa$. Rohacell WF51 rigid foam, also defined by Koissin et al. (2004), is used in the core with $E=85MPa$, $\nu =0.42$ and $\overline{\sigma }=0.9MPa$, resulting in the following properties used in the proposed model $K=177MPa$ and $G=29.93MPa$.

Koissin et al. (2004) did not inform the core post-critical behavior, only makes mention that it has been modeled as crushable foam, thus $\overline{\tau }=0.4MPa$ and $H^1=H^2=-0.485MPa$ (softening at small strains) has been adopted. Following the numerical tensile and compression tests proposed in example 5.1, the core mechanical behavior adopted for the numerical analysis is shown in Figure 16. No additional adaptation is performed in the model parameters to show that, even with softening, after a certain strain level the constitutive model provides the necessary stiffening to prevent the volumetric degeneration of the material. The higher the Poisson's ratio the earlier this phenomenon occurs.

Figure 16
Longitudinal behavior of the adopted core material

Considering symmetry, the numerical simulation is performed using the in house computer code and the employed discretization is shown in Figure 17. It has $2398$ nodes and $248$ finite elements of cubic approximation at the base and linear in depth (width of the specimen). All actuator nodes are prevented to move along horizontal directions. The actuator vertical movement is imposed by means of node $A$, see Figure 17, for which the actuator force is measured. The actuator is stiffer than the specimen with the following mechanical properties: $K=4GPa$ and $G=8.99GPa$. The adopted penalization constant (contact) is $0.5kN$.

Figure 17
Employed discretization ($xy$ plane)

In Figure 18, the forces of experimental and numerical cases are compared as a function of the actuator displacement. Considering the complexity of the problem, the numerical result conforms very well with the experimental one, revealing the capacity of the proposed constitutive model to accurately represent parts of sandwich panels subjected to near concentrate forces. It is also observed that at the end of the analysis the numerical stiffness of the set increases, it occurs because the core constitutive model entered its ascending part, always present to avoid material self-intersection. This stiffening at the end of the curve could be avoided if more information regarding the mechanical behavior of the core, for example, reduction of bulk modulus for higher deformations, was available.

Figure 18
Load versus actuator displacement

An important detail of the achieved results is the detection of the depression present in the graph of Figure 18. It occurs in the region of $0.85mm$ imposed displacement. The reproduction of this depression surprised even the author. Next to the $2mm$ displacement level, there is another depression detected by the numerical model, however this phenomenon is not evident in the experimental results presented by Koissin et al. (2004). The causes of these depressions are identified as abrupt realignments of the core material at certain deformation stages. It should be noted that the maximum number of iterations allowed for the whole analysis is 31 and the adopted tolerance (in position) is $5.0x{10}^{-5}$.

Figure 19 presents the final configuration snapshot indicating vertical and horizontal displacements, as well as, the developed plastic strain ${\lambda }_1^{\zeta }$, related to ${\overline{I}}_1$

Figure 19
Snapshot for final displacement (8mm)

6 CONCLUSIONS

In this study a new elastoplastic model is proposed and successfully implemented in the positional FEM. The model is developed to perform applications in sandwich panels with soft core, mainly representing the core crushing. The model is explored theoretically, revealing its ability to simulate plasticity in large strains, including cyclic loads and preventing the volumetric degeneration of the material. Experimental tests are reproduced with good precision by the proposed numerical technique. It is concluded that the constitutive model is suitable for the intended applications with a simple implementation. Future developments include its generalization to orthotropic materials and its application to dynamic problems, where impacts on small regions of sandwich panels can cause major core damage.

7 ACKNOWLEDGEMENTS

This research has been supported by the São Paulo Research Foundation, Brazil - Grant #2020/05393-4.

References

• Abrate S., Di Sciuva M. (2017) Equivalent single layer theories for composite and sandwich structures: A review, Composite Structures 179: 482-494
• Abraham, F.F., Walkup, R., Gao, H., Duchaineau, M., Diaz De La Rubia, T., Seager, M. (2002) Simulating materials failure by using up to one billion atoms and the world's fastest computer: Brittle Fracture. Proc. Natl. Acad. Sci. 99:577-578
• Argyris J.H., Kleiber M. (1977) Incremental formulation in nonlinear mechanics and large strain elasto-plasticity - natural approach. Part 1. Comput. Methods Appl. Mech. Eng. 11 (2):215-247
• Atluri S.N. (1984) On constitutive relations at finite strain: hypo-elasticity and elasto-plasticity with isotropic or kinematic hardening. Comput. Methods Appl. Mech. Eng. 43(2):137-171
• Brepols T., Vladimirov I.N., Reese S. (2014) Numerical comparison of isotropic hypo- and hyperelastic-based plasticity models with application to industrial forming processes. International Journal of Plasticity 2014; 63:18-48
• Bruhns O.T., Xiao H., Meyers A. (1999) Self-consistent Eulerian rate type elasto-plasticity models based upon the logarithmic stress rate. Int. J. Plasticity 15(5):479-520
• Buehler M.J., Hartmaier A., Gao H., Duchaineau M., Abraham F.F. (2004) Atomic plasticity: description and analysis of a one-billion atom simulation of ductile materials failure. Comput. Methods Appl. Mech. Engrg. 193:5257-5282
• Caliri Jr. M.F., Ferreira A.J.M., Tita V. (2016) A review on plate and shell theories for laminated and sandwichstructures highlighting the Finite Element Method, Composite Structures 156: 63-77
• Carrazedo R., Paccola R.R., Coda H.B. (2018) Active face prismatic positional finite element for linear and geometrically nonlinear analysis of honeycomb sandwich plates and shells, Composite Structures‏ 200:849-863
• Carrazedo R., Paccola R.R., Coda H.B. (2020) Vibration and stress analysis of orthotropic laminated panels by active face prismatic finite element., Composite Structures e:112254
• Casey J., Naghdi P.M. (1981) A correct definition of elastic and plastic deformation and its computational significance, J. Appl. Mech. 48(4):983-984
• Chen L., Du, B., Zhang J., et al. (2016) Numerical Study on the Projectile Impact Resistance of Multi-Layer Sandwich Panels with Cellular Cores, Latin American journal of Solids and Structures 13 (15):‏ 2576-2595
• Coda H.B., Carrazedo R. (2017) Triangular based prismatic finite element for the analysis of orthotropic laminated beams, plates and shells. Composite Structures 168:234-246
• Coda H.B., Sanches R.A.K., Paccola R.R. (2020) Alternative multiscale material and structures modeling by the finite-element method, Engineering with Computers, http//DOI: 10.1007/s00366-020-01148-y.
  » http//DOI: 10.1007/s00366-020-01148-y.
• Eterovic A.L, Bathe K.J. (1990) A hyperelastic-based large strain elasto-plastic constitutive formulation with combined isotropic-kinematic hardening using the logarithmic stress and strain measures. Int. J. Numer. Methods Eng. 30:1099-1114
• Flory, P.J. Thermodynamic relations for high elastic materials. Transactions of the Faraday Society, v. 57, p. 829-838, 1961.
• Garg A., Belarbi M-O., Chalak H.D. (2021) A. Chakrabarti. A review of the analysis of sandwich FGM structures, Composite Structures258: 113427
• Hartmann S., Neff P. (2003) Polyconvexity of generalized polynomial-type hyperelastic strain energy functions for near-incompressibility. International Journal of Solids and Structures 40:2767-2791
• Hughes T.J.R., Winget J. (1980) Finite rotation effects in numerical integration of rate constitutive equations arising in large-deformation analysis. Int. J. Numer. Methods Eng. 15(12):1862-1867
• Icardi U., Urraci A. (2019) Free Vibration of flexible soft-core sandwiches according to layerwise theories differently accounting for the transverse normal deformability, Latin American journal of Solids and Structures ‏ 16 (8): e230,‏ 2019
• Koissin V., Shipsha A., Rizov V. (2004) The inelastic quasi-static response of sandwich structures to local loading. Composite Structures 64 (2):129-138
• Kojic M., Bathe K.J. (1987) Studies of finite element procedures - stress solution of a closed elastic strain path with stretching and shearing using the updated LagrangianJaumann formulation. Comput. Struct. 26(1-2):175-179
• Kröner E. (1959) Allgemeine kontinuumstheorie der versetzungen und eigenspannungen. Arch. Ration. Mech. Anal. 4: 273.
• Lee E.H. (1969) Elastic-plastic deformations at finite strains. Journal of Applied Mechanics (ASME) 36:1- 6
• Lolive E., Berthelot J.M. (2002) Non-linear behaviour of foam cores andsandwich materials. Part II: indentation and three-point bending. J Sandwich Struct & Mater 4 (4):297-352.
• Mandel J. Plasticite Classique et Viscoplasticity. In: CISM Course 97, Springer-Verlag, Udine, 1971.
• Marques D., Flor F.R., de Medeiros R., et al. (2018) Structural Health Monitoring of Sandwich Structures Based on Dynamic Analysis, Latin American journal of Solids and Structures 15 (11): e58
• Naik N.S., Sayyad A.S. (2018) 2D analysis of laminated composite and sandwich plates using a new fifth-order plate theory, Latin American journal of Solids and Structures 15 (9): e114
• Nsengiyumva W., Zhong S., Lin J., Zhang Q., Zhong J., Huang Y. (2021) Advances, limitations and prospects of nondestructive testing and evaluation of thick composites and sandwich structures: A state-of-the-art review, Composite Structures 256: 112951
• Ogden RW. Nonlinear Elastic deformation. Ellis Horwood, England, 1984.
• Olsson R. (2002) Engineering method for prediction of impact response and damage in sandwich panels. J Sandwich Struct & Mater 4(1):3-29
• Pascon, J. P.; Coda, H. B. (2013) A shell finite element formulation to analyze highly deformable rubber-like materials, Latin American journal of Solids and Structures ‏ 10 (6): 1177-1209,
• Rivlin R., Saunders D. (1951) Large elastic deformations of isotropic materials VII. Experiments on the deformation of rubber, Philos.Trans. R. Soc. London Ser. A 243:251-288
• Salami S.J., Sadighi M., Shakeri M. (2015) Improved high order analysis of sandwich beams by considering a bilinear elasto-plastic behavior of core: An analytical and experimental investigation. International Journal of Mechanical Sciences 93:270-289
• Sayyad A.S., Ghugal Y.M. (2017) Bending, buckling and free vibration of laminated composite and sandwich beams: A critical review of literature, Composite Structures 171,: 486-504.
• Shekari A., Ghasemi F.A., Malekzadehfard K. (2017) Free Damped Vibration of Rotating Truncated Conical Sandwich Shells Using an Improved High-Order Theory, Latin American journal of Solids and Structures 14 (12):‏ 2291-2323,
• Simo J., Ortiz M. (1985) A unified approach to finite deformation elastoplastic analysis based on the use of hyperelastic constitutive equations. Comput. Methods Appl. Mech. Eng. 49:221-245
• Simo J.C. (1992) Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory. Comput. Methods Appl. Mech. Eng. 99:61-112
• Siqueira T.M., Coda H.B. (2016) Development of Sliding Connections for Structural Analysis by a Total Lagrangian FEM Formulation, Latin American journal of Solids and Structures 13 (11):‏ 2059-2087
• Soden P.D. (1996) Indentation of composite sandwich beams. J Strain Anal for Eng. 31 (5):353-360.
• Thomsen O.T. (1993) Analysis of local bending effects in sandwich plates with orthotropic face layers subjected to localized loads. Composite Structures 25:511-520.
• Wallin M., Ristinmaa M., Ottosen N.S. (2003) Kinematic hardening in large strain plasticity, European Journal of Mechanics A/Solids 22:341-356
• Weber G., Anand L. (1990) Finite deformation constitutive equations and a time integration procedure for isotropic, hyperelastic-viscoplastic solids. Comput. Methods Appl. Mech. Eng. 79:173-202
• Zhang M., Montán F.J. (2019) A simple formulation for large-strain cyclic hyperelasto-plasticity using elastic correctors: Theory and algorithmic implementation. International Journal of Plasticity 113:185-217

Edited by

Editor:

Marco L. Bittencourt.

Publication Dates

• Publication in this collection
  27 Sept 2021
• Date of issue
  2021

History

• Received
  20 May 2021
• Reviewed
  22 July 2021
• Accepted
  02 Aug 2021
• Published
  02 Aug 2021