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

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.


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 Latin American Journal of Solids and Structures, 2021, 18(6), e392 3/26 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 andCoda andCarrazedo, 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.

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 φ  which has the dimensionless space ) in which 0 i f is the initial configuration mapping, i x represents the initial domain points coordinates, j ξ are the dimensionless coordinates, i X α are initial coordinates of nodes α .In the same way, 1 i f is the current configuration mapping, i y represents the current domain points coordinates and i Y α are current nodal positions.( ) and its gradient is given by: ( ) It is important to mention that 0 A is known at all integration points (calculated only once) and 1 A 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: S is the complete second Piola-Kirchhof stress, E is the Green-Lagrange strain, 0 V is the initial volume of the analyzed body and 0 Γ is its initial surface.
Remembering that current positions of nodes are the unknowns and using the approximation defined by equation (1), one writes: As Y δ  is an arbitrary position variation, the nonlinear system of equations to be solved is: in which different shape functions are used for domain φ  and surface ϕ  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 c S 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 , i.e., the volumetric part of the complete stress will be elastic and the two deviatoric parts of the stress will be elastoplastic.

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 Latin American Journal of Solids and Structures, 2021, 18(6), e392 5/26 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): where ( ) here Â is the volumetric part of the deformation gradient A and A is its isochoric part.
As the adopted strain is the Green-Lagrange ( ) , written as a function of the right Cauchy stretch tensor t ⋅ C = A A , it is necessary to write, from equation (10), the isochoric part of the right Cauchy-Green stretch as: In which ( ) 1 Det C ≡ .The volumetric part of C is defined from equation (9) as: where 2 ( ) ( ) Multiplying equations ( 11) and ( 12) results the necessary multiplicative decomposition of the right Cauchy-Green stretch tensor as: For hyperelastic materials the specific strain energy is written using an elastic potential of the form in which 1 I and 2 I are, respectively, the first and second invariants of the isochoric part of the Cauchy-Green stretch C and J is the third invariant of its volumetric part Ĉ .
Using the concept of energy conjugate, the second Piola-Kirchhoff stress tensor is written as: or in compact form 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: ( ) Latin American Journal of Solids and Structures, 2021, 18(6), e392  6/26 in which K is the bulk modulus, G is the shear elastic modulus and 0 n > 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: ( ) where the different strain components, that results in stress "directions", are given by: In what follows it is shown that the second Piola-Kirchhoff stress components ( ) , , J S S S defined by equation ( 17) are related, respectively, to the hydrostatic and deviatoric Cauchy stress components ( ) . This is done here applying the well known relation among the second Piola-Kirchhof stress and the Cauchy stress, as follows: ( )

A S A A S S S A A S A A S A A S A (21)
To make this operation simple one starts rewriting equation ( 16) in a general form: ( ) Thus, for each component one writes: a) Volumetric The candidate to the volumetric component of the second Piola-Kirchhoff stress is written as: In which α is a scalar, see equation ( 17) and J E is given by equation ( 18).
Applying Eq. ( 21) over the second Piola-Kirchhoff stress of Eq. ( 23) results: Latin American Journal of Solids and Structures, 2021, 18(6), e392  7/26 which means that J E is the Lagrangian strain direction corresponding to the hydrostatic stress direction in the Cauchy space.Thus, J S is hydrostatic in Lagrangian sense.

b) First isochoric
The candidate to the first isochoric component of the second Piola-Kirchhoff stress is written as: where β is a scalar and 1 E 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: ( ) And, as one achieves: ( ) which means that 1 E is the first isochoric Lagrangian strain direction and 1 S 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: with γ and 2 E 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: Using equation ( 27) results: From Eq. (31) one calculates ( ) meaning that 2 E is the second isochoric Lagrangian strain direction and 2 S 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 1 S or 2 S (without volumetric components) reach a critical value.

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.

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: in which y σ is the reference yielding stress and dev σ is the elastic deviatoric part of the Cauchy stress tensor with 2 J being its second invariant.In order to develop calculations (choosing the shear stress as a reference) one defines 2 y σ τ = , resulting: As shown in section 3, before yielding takes place, the deviatoric elastic stress can be represented by the Lagrangian components ( 1 S and 2 S ) 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: is assumed, following the elastic relation (17).Other way to write (36) is:  and Structures, 2021, 18(6), e392 9/26 in which the isochoric "directions" 1 E and 2 E are present.In order to simplify operational inferences, one rewrites equation ( 38) in an unified form, with E representing both 1 E or 2 E : 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 1 x ).For small strains and the proposed test one has, in which i λ represents stretch in direction i .Remembering that the material is assumed isotropic, for this uniaxial test one writes 2 3 Making some algebraic simple steps one achieves from equations ( 19) and ( 20): and Using equation ( 17), the following particular relation is achieved: Using equation ( 21) and remembering that for small strains one has 1 1 λ ≅ results: verifying the proposition for small strains.For large strains a simple calibration of τ can be made.

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.a -Plastic evolution
The local plastic evolution is established from the observation that the general strain "direction" E is isochoric, i.e.: : where / : E E E represents the unitary deviatoric directions of plastic strain evolution p ∆E , and ∆λ 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 p ∆S 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 p p p ∆ = + S S S cannot be used because, for large strains, a past value of p ∆S 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:

4
: and the evolution is imposed on the scalar plastic strain ζ λ as: ( ) in which ac represents cumulated and ζ is the sign of plastic evolution, to be defined after the introduction of kinematic hardening, see equation (50).As p S is proportional to E , hypothesis (ii) and hypothesis (iii) are fulfilled.After the determination of the plastic multiplier ∆λ , 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 0 J > 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 ( c H ) is considered constant at iterations of the numerical solution process, but depends upon the scalar plastic strain ζ λ for each load level.Defining the ratio between the kinematic hardening and the shear elastic modulus as: the internal variable of the kinematic hardening becomes: Latin American Journal of Solids and Structures, 2021, 18(6), e392 11/26 ( ) and the signal of the scalar plastic strain evolution becomes ( ) The internal variable q is multiplied by the isochoric direction and the elastic modulus to generate the back stress tensor as

4
: 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 κ associated with isotropic hardening is calculated by a cumulative process of the plastic multiplier ∆λ , see equation ( 47).The isotropic hardening ( i H ) is considered constant during iterations, but depends upon the cumulated plastic multiplier.Similarly to the kinematic hardening, one defines: The plastic multiplier λ is cumulated by and the evolution of the internal variable κ is given by: ( )

Preparation of the objective function
Adopting the symbol t S for both isochoric Piola-Kirchhoff stresses 1 S or 2 S of equation ( 17) and considering previous definitions, the stress to be introduced in the yielding function (objective function) is defined by: or 1 4 4 4 4 : : : : Thus, one calculates: It is important to remember that ( / 4)(1 ) :  5) and ( 7); and that ( )

S
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 κ ), equation ( 39) is upgraded to the current yielding function or objective function, as: The objective function f given in equation ( 57) has the usual meaning, i.e, for 0 f < 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)).

Calculating the plastic multiplier ∆λ
Plastic evolution occurs when a trial objective function violates the plastic criterion (51), that is: 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: As hardenings are considered constant at iterations, equation ( 59) is easily solved and the plastic multiplier is the smaller positive value of: For cyclic applications, when plastic stress ( p S ) exists and : C C 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.

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 Resulting the proposed model algorithm tangent constitutive tensor as: It is important to mention that the volumetric part of the complete constitutive tensor is elastic.

Small strain range -plastic evolution
When small strains occur with 0 p ≠ S (cyclic loads) a small strain procedure is necessary due to an ill conditioning of equation ( 47), i.e., when the quantity / : E E E passes through its nullity, it loses the ability of representing the total plastic value p S , given in equation ( 47), but when p S 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 1 S and 2 S the total quantity to be verified in the yielding criterion is: in which χ is the back-stress tensor.The pure elastic stress E and the accumulated plastic stress p S come from a previous existence of the large strain procedure, so are isochoric.From these tensors one writes the plastic direction as:

S S S S S S S S S S
(67) resulting: Thus, the objective function becomes: ( ) Evolutions are defined by: in which λ ∆ is the evolution of plastic multiplier.When 0 tr f ≥ one imposes equality in equation ( 69) together with plastic multiplier evolution, resulting: ( ) in which one uses to write it short.After some algebraic manipulations equation ( 71) turns into: in which . As hardening is considered constant during iterations, one solves equation ( 72) as: ( ) being the smaller positive value of ∆λ the searched solution.The same algorithm constitutive elastoplastic tensor of finite strains, given by equation ( 65), is used for small strains.

Transition between models
Small and large strain variables are the same, but for small strains the tensor values p S and χ should be stored, while for large strains the correspondent scalar values ζ λ 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 Finally the transition is made from large strains to small strains when

:
Tol < E E and from small strains to large strains when

:
Tol > E E . It is worth mentioning that the process always start using the large strain procedure, because in this situation 0 p = S .

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.

Cyclic axial test
This example consists in a unitary cub allowed to slide at faces  .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.  .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 1 2 0.5 , with a larger stretch interval 0.4 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.
Latin American Journal of Solids and Structures, 2021, 18(6 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.In order to show the extension of the applied strain, Figure 7 presents snapshots with Cauchy stress values for cases

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 .
Latin American Journal of Solids and Structures, 2021, 18(6), e392 19/26 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 25.0 F N ∆ = is equally applied in 21 steps for half structure.The load is distributed over an area of 7 50 mm x mm 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 .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 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  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 0 x = ) 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.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 ζ λ of face-sheets (zoom close to the middle of the beam) for a applied load of 600 F N = . Figure 14b shows the total displacement of the core following z direction for 600 F N = . Case (c) is considered for both Figures 14a and 14b

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 12.5 R mm = is used to impose the vertical displacement of the face-sheet as shown in Figure 15.(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.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 selfintersection.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.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 Latin American Journal of Solids and Structures, 2021, 18(6), e392 24/26 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

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.

Figure 1 :
Figure 1: Deformation function described by prismatic elements and auxiliary dimensionless space The deformation function f  is represented by the following composition: Using equation (17) one writes equation (37) as a function of the deviatoric part of the hyperelastic constitutive model, i.e.:LatinAmerican Journal of Solids see Figure2a.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 80 K = , shear elastic modulus 1 and kinematic hardenings for each case.A simple discretization with only two linear elements is employed, see Figure2b.

Figure 5 -
Figure 5 -Kinematic hardening for large strain -stress cycle

Figure 8 -
Figure 8 -Three point bending experiment and used discretization

Figure 10 -
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 ζ λ .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 500 F N = and are localized at the low corner of the right side of Figures 11 and 12. .

Figure 11 -Figure 12 -Figure 13 -
Figure 11 -Plastic strain level ζ λ for polyurethane core, case (c) Figure 14 -Plastic strain ζ λ for face-sheets and displacement in z direction for 600 F N =

Figure 15 -
Figure15-Schematic representation of half specimen proposed byKoissin et al. (2004) The same material used byKoissin et al. (2004)(GRFP) is adopted to build up the face-sheet with the experimental elastic modulus 16.5 E GPa = and Poisson ratio

Figure 16 -
Figure 16 -Longitudinal behavior of the adopted core material

Figure 18 -
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

Figure 19 IFigure 19 -
Figure 19 presents the final configuration snapshot indicating vertical and horizontal displacements, as well as, the developed plastic strain 1 ζ λ , related to 1 I

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