Indentation of Sandwiches Using a Refined Zig-Zag Model and a Mesoscale Damage Model

An accurate and cost effective zig-zag plate model with variable kinematics and fixed degrees of freedom recently developed by the authors, which a priori fulfills the stress and displacement continuity requirements, is applied to study indentation of sandwiches with honeycomb/foam core. The variable elastic properties of the core during crushing are evaluated apart and once at a time through a 3D finite element analysis. Shell elements with elastic-plastic isotropic properties are adopted to discretize honeycomb cores and faces, instead solid elements with nonlinear material behavior are used for foam cores. To keep low the computational burden, the structural response is computed in closed form by the zig-zag model using these elastic properties. The damage analysis is carried out through a mesoscale model that determines the degraded properties of failed regions in a physically consistent way. A good agreement with experiments taken from the literature being shown, the present simulation procedure with a low computational effort is proven to be an efficient alternative to the ones currently employed.


Introduction
The continuous quest for lighter and stronger structures has favored the spread of sandwiches with laminated faces in the aerospace industry and in many other engineering applications. However, a full exploitation of the advantages offered by these materials (low weight, excellent strength, stiffness, energy absorption and fatigue properties) requires an accurate modeling of local damage rise and growth in service, due to the inability of sandwiches to yield and dissipate a large amount of energy via plasticity, as thoroughly discussed in the book by Daniel and Ishai [1] and in the recent papers by Davies et al. [2] and de Freitas et al. [3].
Damage of sandwiches starts by the crushing of core followed by tearing of the loaded face. A comprehensive discussion of the mechanisms of damage formation and evolution and of their modeling is given by Garnich and Akula Venkata [4] and Liu and Zheng [5]. Parameters such as topology of cells, their relative density and the thickness of the foil strongly affect this phenomenon. Instead, the damage of laminated faces begins as matrix cracking, fibres breakage and local puncture. Then, when cracks reach an interface and deflect onto it, delamination develops. Distribution, magnitude and duration of stresses, load history, stacking sequence, strengths and pre-existing defects have a significant bearing.
Great attention should be paid to the analysis of indentation (i.e. low-velocity impact), since it can produce barely visible damages with a significant influence on the structural performances and the service life of laminated and sandwich composites. Accordingly, low-velocity/ low-energy impacts are extensively investigated in the literature. As examples among many others, the recent papers [6]- [11] are cited. Warping, shearing and straining deformations of the normal and out-of-plane stresses rise, due to the much smaller moduli and strengths in the thickness direction than in the in-plane direction, which constitute the so-called layerwise effects. As they strongly influence the rise and growth of damage, they should be accurately accounted for by the models, which should maintain an affordable effort in an industrial environment. Moreover, an accurate modeling of the transverse normal stress and deformation is required for describing the core's crushing.
Indentation is often studied using three-dimensional (3-D) finite element analyses (FEA) and fracture mechanics or cohesive zone models (see, e.g. Wang et al. [12], Borg et al. [13] and Nishikawa et al. [14]). In this way, accurate results are obtained, though at the expense of a high computational effort. A detailed finite element simulation of the cellular structure of honeycomb by plate elements is often given for capturing localized phenomena like the buckling of cell walls, while foam core is simulated using solid elements.
Recently developed high-order layerwise plate (HLW) models based on a combination of global higher-order terms Universal Journal of Mechanical Engineering 2(1): 6-19, 2014 7 and local layerwise functions have been proven to be equally accurate than 3-D FEA with a much lower computational effort (see, e.g. . Frostig [15], Plagianakos and Saravanos [16], Zhen and Wanji [17], Elmalich and Rabinovitch [18] and Giunta et al. [19]). An extensive discussion of the techniques used by these models to account for the layerwise effects, extensive assessments of their structural performances and of related finite element models are given among many others in the recent papers by Matsunaga [20] and Chakrabarti et al. [21]. The crucial point is whether HLW models have a number of unknowns that depends or not on the number of physical or computational layers, as this determines accuracy and computational costs. The models with a variable number of unknowns are widespread, as they can be refined across the thickness for improving accuracy, but the so-called zig-zag models embodying continuity functions that a priori satisfy the displacements and stress contact conditions can be equally accurate and less expensive. Successful indentation studies using (HLW) models and stress-based failure criteria have been recently performed e.g. by Diaz-Diaz et al. [22] and Qiao and Yang [23] and using a zig-zag model by Icardi and Ferrero [24], which proven that accurate results can be obtained with a low computational effort.
In this paper, a special HLW model is adopted as structural model, the zig-zag model with variable kinematics across the thickness and fixed degrees of freedom (d.o.f.) recently developed by Icardi and Sola [25] being used, which is here referred as FD_VK. This model with a hierarchic representation of displacements across the thickness was developed for obtaining results as accurate as those of three-dimensional models directly from constitutive equations, employing just 5 functional d.o.f., as in the equivalent single layers models. Differently to what done in Refs. [25] and [26], here a reduced order of expansion and an increased number of computational layers are used in the FD_VK model, since numerical tests have proven this to be more efficient for sandwiches. Please notice that this choice does neither result into a considerably larger computational effort nor into a larger memory storage dimension, because irrespectively of the used number of computational layers the number of d.o.f. is fixed. In effects, as shown in [25], the model requires a processing time that is always comparable to that of equivalent single layer models.
As the FD_VK cannot properly treat the crushing of soft cores, being a homogenized macro-mechanic model, the crushing behavior is evaluated apart once at a time through a detailed 3D FEA from which the variable, apparent elastic properties of core under compressive loading are determined. Specifically, plate elements are used to discretize the face sheets and the cellular structure of honeycomb, while foam core is discretized by solid elements with non-linear material behavior. In this way, the most significant local phenomena can be described, since the properties are assumed to vary with the applied load and from the points over the contact area to the neighboring points, according to the results by the detailed finite element analysis.
The overall computational burden is kept as low as possible, since after having determined the properties of core, the analysis is carried out by the low cost FD_VK model. Here the studies are carried out in closed form because just simple test cases are considered, however structures of industrial interest can be analyzed via finite elements using the Strain Energy Updating Technique (SEUPT), as outlined in [24] and discussed forward. As customarily, the onset of damage is predicted using stress-based criteria.
Instead of guessing the residual properties after failure through appropriate multiplication factors applied to the elastic moduli, e.g. in the ply discount theory, then extending a pre-existing damage to the points where the ultimate condition is reached, in this paper the residual properties and the damage growth are described using the mesomechanic model by Ladevèze et al. [27]- [28]. It considers the effects of the damages on the microscale through a modified expression of the strain energy. This damage model was chosen because, as explained in [29], it provides greatly accurate results with a smaller computational burden than structural scale models.
Even if many authors agree that indentations can be studied using just static simulations, here a dynamic solver based on the Newmark's implicit time integration scheme is developed, which can treat general transient dynamic problems of practical interest. The indentation depth and the contact area are computed assuming the distribution of the contact force to be Hertzian and the projectile as a rigid body, as customarily. At any time step, the iterative algorithm by Palazotto et al. [30] provides the contact area, forcing the surface of the target to conform the shape of the impactor, as required by soft media. Figure 1 sums up the main steps of the procedure here adopted, while its details are provided in the following Sections.
As a general remarks it could be underlined that the present simulation methodology differs from the ones already proposed in the literature since: (i) the adopted structural model tries to more efficiently simulate the behavior of the structure, having a fixed number of unknowns variables though it allows for a variable kinematic representation in different regions across the thickness; (ii) as damage model a mesoscale damage with a physically consistent representation of the damage mechanisms is applied.
The present simulation methodology differs from that of [24] because the improved structural model here used can more accurately describe the effects of transverse shear and normal stresses and deformation involved by the mechanisms of damage. Since the structural model allows for a hierarchic representation across the thickness, it enables an accurate prediction of the stresses from constitutive equations, while the model of Ref. [24] required post-processing. Moreover, a refined damage model is used. In addition, in the present case, the solution is found in closed form, while, previously, SEUPT was used for carrying out a finite element analysis, in order to assess the technique.  Applications will be presented to sample test cases available in literature. The comparison of the present results to experiments will show the capability of the present simulation procedure to accurately predict the indentation behavior of sandwiches keeping low computational costs. This is due to the fact that the memory storage dimension are lower than for currently available models used in the impact studies and the processing time is comparable to that of smeared laminate models.

Structural Model
As extensively discussed e.g. by Zhang and Yang [31], who reviewed HLW models, the displacements should be continuous across the thickness, but their derivatives should be discontinuous at the material interfaces in order to make continuous the transverse shear and normal stresses at the layer interfaces, as prescribed by the elasticity theory. In addition, the stress-free boundary conditions at the upper and lower bounding faces should be satisfied. Moreover, in order to respect the 3D solution, the same order of polynomials cannot be used for all the displacements and for any point across the thickness, thus the representation should vary with the type of loading, the boundary conditions, the thickness ratio and the constituent materials.
Trying to reach a compromise between accuracy and computational costs saving, many models have been developed that fully or partially satisfy these requirements. As mentioned in the introductory section, the most important feature by the viewpoint of accuracy and computational costs is whether the number of variables depends on the number of physical layers (see, e.g., Moreira and Rodrigues [32], Chrysochoidis and Saravanos [33], Hohe and Librescu [34], and Shariyat [35]) or it is fixed (see, e.g., Zhen and Wanji [17] and Icardi and Sola [25]). In the former case, accurate estimates are achieved at the expense of a high computational burden, while in the latter case complex algebraic manipulations are required in order to derive the expression of the unknowns as function of the d.o.f. However these operations can be carried out apart once at a time using symbolic calculus [25], thus consistently speeding-up the analysis. This aspect fosters the application of models with fixed d.o.f. to the analysis of structures of industrial complexity and of highly iterative problems, such as impacts.
In the present paper, a revised version of the zig-zag model [25] is used as structural model where a reduced order of expansion and an increased number of computational layers are used that make it more efficient, as discussed forward along with its features, merits and advantages. Sandwiches are treated as multi-layered composites, with a thick intermediate layer constituting the core, whose properties vary with the applied load, as done by all modern high-order layerwise models. They are referred to a rectangular Cartesian reference frame (x, y, z) placed on their middle plane Ω, with z normal to Ω, and x, y on Ω. The symbols (k) z + and (k) zindicate the position of the upper + and lowersurfaces of the k th layer, the quantities that belong to a generic layer k are denoted with the superscript (k) . The symbols u, v, w represent respectively the elastic displacements in the direction of x, y and z.
The through-the-thickness variation of the displacements is postulated in the following general, piecewise form [25]: Universal Journal of Mechanical Engineering 2(1): 6-19, 2014 The symbols with the superscript 0 represent the contributions that are fixed in all the layers, those with the superscript i are the terms that vary between the different layers and finally the superscript c characterizes the contributions that enable to fulfill the continuity conditions prescribed by the elasticity theory. The explicit expressions of the fixed terms are: Equations (4)-(6) just contain the in-plane displacements of the points over the middle plane u 0 (x, y), v 0 (x, y), the transverse displacement w 0 (x, y) and the shear rotations of the normal γ x 0 (x, y), γ y 0 (x, y) to these points. These quantities are the functional d.o.f.
The expression of the layers contributes is: These terms are valid within a single layer of the structure; however the number of d.o.f. is fixed, since their expressions are computed apart as function of the d.o.f using a symbolic calculus tool. Their expression are obtained imposing the constraints prescribed by the elasticity theory (i.e. the stress-free boundary conditions at the upper and lower bounding faces) and the equilibrium conditions at n p =N lay • ord_u -2 points across the thickness, N lay being the number of computational layers and ord_u the chosen expansion order. The locations of the n p points are picked arbitrarily, however they should not be placed excessively near to the interfaces in order to avoid singularity. By properly choosing the expansion order and the number of subdivision across the thickness, accurate predictions of the transverse shear and normal stresses are obtained directly from constitutive equations, even for thick structures with abruptly changing material properties.
Finally, the expressions of the terms that enable the fulfillment of the continuity constraints are: As prescribed by the elasticity theory for keeping equilibrium, the interlaminar stresses, the transverse normal stress and its gradient should be continuous at the interfaces. To this purpose, the expressions of Φ x k , Φ y k are obtained enforcing the continuity of the transverse shear stresses, while Ψ k , Ω k are obtained enforcing the continuity of the transverse normal stress and its gradient. Finally, C u k , C v k and C w k make continuous the displacements at the points across the thickness where the representation is varied. H k is the Heaviside unit step function that activates the contribution of the continuity functions from the pertinent interface k th .
Similarly to what done for the high order terms, the explicit expressions of the continuity functions are derived in closed form using a symbolic calculus tool. As shown in [25], the closed-form expressions of high order terms and continuity functions involve the functional d.o.f. and their derivatives. As a consequence, when finite elements are developed, high-order shape functions and exotic nodal d.o.f. should be adopted. This drawback can be overcome with the strain energy updating technique (SEUPT) [24], which improves the accuracy of standard C 0 shear-deformable plate elements up to the level of the FD_VK model through fast post-processing operations.

Crushing Behavior
The soft material constituting the core of sandwiches plays an important role in determining the energy absorption of the structure. The typical behavior is characterized by an initial linear phase at low strains, followed by a collapse phase, which finally turns into densification, as shown in many works taken from literature (see, e.g. Mohmmed et al. [36], Aminanda et al. [37], Aktay et al. [38], Flores-Johnson et al. [39] and Ivañez et al. [40]).
This behavior is strongly affected by shear stress, because the ratio of this stress to other stresses determines the energy absorption of the core. Other parameters that play an important role in determining the core crushing behavior are the ratio of skin thickness to side length and the core density.
As far as simulation techniques are concerned, crushing is generally described through a detailed finite element simulation of the cellular structure of the honeycomb (Aminanda et al. [37] and Aktay et al. [38]), or using solid elements for foam core (Mohmmed et al. [36] and Ivañez et al. [40]). In details, honeycomb's cells are simulated using shell elements coupled with a self-contact algorithm that prevents from interpenetration between the folds in the cell walls. Instead, a proper modeling of foam core calls for the employment of suited homogenized material properties that are determined from experiments.
Of course, the refinement of the mesh used to discretize the soft core has an important role in describing the material behavior. In fact, a coarse discretization determines a higher initial load peak and higher load in the densification regime. On the other hand, very refined mesh could be unsuited if the element size is not larger than the thickness of the honeycomb foil.
In order to keep the computational burden as low as possible, the soft material constituting the core can be treated as homogeneous material, determining the properties according to Gibson and Ashby [41]. This approach was followed in previous applications of the FD_VK model (see, e.g., [25] and [26]), since therein only the global response of the structure was required. Here the attention is focused on indentation, which instead needs to capture with the due accuracy all the local phenomena, such as the buckling of cell walls, in order to accurately describe the onset of damage and its growth. Accordingly, the apparent elastic moduli of the core while it collapses/buckles under transverse compressive loading are here determined with a detailed finite element analysis, which is performed apart once at a time. The evaluation of the elastic moduli is done in any sub-region in which the area around the contact point is discretized. Using the results by this analysis, the FD_VK model can then determine stress and displacement fields taking into consideration all the localized phenomena. In fact, though the core is treated as a homogeneous material by the FD_VK model, the effects of crushing are accounted for locally assuming elastic properties variable with position. The updated Lagrangian methodology is used, which evaluates the deformations at each new time step from the geometry of the structure at the previous step, not from the initial unloaded configuration. In this way, the nonlinear effects, considered when computing the apparent, variable elastic properties of the core, are accounted for during the structural analysis, even if the FD_VK model is based on infinitesimal strain-displacement relations. This methodology is chosen because it has an improved numerical efficiency than the standard Lagrangian approach, as claimed by Bathe [42].
With the above described simulation procedure, computational costs are reduced, since the detailed finite element analysis is limited to the preliminary phase, but, at the same time, accuracy is preserved, because the accuracy of the FD_VK model is similar to that of 3D models with a much lower computational effort.
During the preliminary finite element analysis, the honeycomb structure is described employing shell elements with elastic-plastic isotropic material, paying attention to adopt a self-contact algorithm to prevent from interpenetration of the folds in the cell walls. The cell walls bonded together are modeled using double shell elements. It is well known that compared to metallic core, Nomex TM honeycomb can present a relatively brittle behavior, since fractures develop after reaching a critical load. Accordingly, in metallic honeycomb the collapse mechanism is plastic and starts after buckling, while Nomex TM honeycomb undergoes buckling of cells, followed by debonding failure at the cell interfaces and fracture. In order to capture this behavior, failure criteria have been considered during the finite element analysis. Foam core is instead discretized using solid elements with nonlinear material behavior, the properties being derived from experiments.
The finite element analysis is started below the lowest failure load predicted by the criteria of Lee and Tsotsis [43], Petras and Sutcliffe [44], Besant et al. [45], the rule suggested by Evonik for Rohacell foam, or the analytical model by Koissin et al. [46] for predicting wrinkling.
To assess accuracy, the Divinycell foam cores analyzed by Flores-Johnson et al. [39] are considered. Experiments were carried out in order to detect the crushing behavior of two cubes with 25 mm edges in Divinycell H100 and H130 foams undergoing compressive load. Figure 2(a) reports the strain-stress behavior for foam cores with different density, under uniaxial compression test. It could be noticed that the numerical results by the present model behave in accordance with the experimental ones.
As a further assessment it is considered the experimental force-displacement curve detected by Aktay et al. [38] for an aluminum 5052 honeycomb with a density of 27 kg/m 3 . The honeycomb has 6 cells in the width direction and 5 in the ribbon one, while its thickness is 50 mm. The cell side is 13,5 mm, while the cell wall thickness is 0,07 mm. Figure 2(b) shows the comparison between experimental and numerical results. Also in this case the simulation procedure provides satisfactory results.

Damage Modeling
The damage in laminates begins as diffuse fibre-matrix debonding, whose accumulation determines transverse matrix microcracking. Their pattern is periodic, at least locally, and it is described in micromechanics by the cracking rate ρ=L/H at the tip of cracks, L being the distance between two adjacent cracks and H being the length of the crack across the thickness. These microcracks determine at their tip a stress concentration, whose consequence is a degradation of the resin film between adjacent plies (i.e. the interface). The next step of the damage evolution is the growth of local cracks perpendicular to matrix microcracks (i.e. local delaminations) that, once coalesced under out-of-plane macroloading, turn into macroscopic interlaminar debonding. The pattern of local delamination cracks is also periodic, and it is quantified at each crack tip by the local delamination ratio τ=e/H, where e is the length of the microcrack. It could be noticed that, according to experiments, the range of practical interest for ρ is [0; 0,7] and for τ is [0; 0,4], because for higher levels of degradation materials are fully damaged.
This micromechanical behavior can be accurately described using fracture mechanic models, but at the expense of a high computational burden. To overcome this drawback, here the damage mesomodel (MSD) by Ladevèze et al. [27]- [28] is chosen, considering that this kind of models has an accuracy comparable to that of the above mentioned models, but an improved efficiency, as claimed by McQuigg [29].
The MSD easily describes the effects of damage, adopting an intermediate scale between the macroscale of the structure and the microscale of the damage mechanism. This method replaces the discretely damaged medium with a continuous homogeneous medium, equivalent from an energy standpoint and incorporating into the expression of the strain energy damage indicators. They represent the homogenized result of damage micromodels and they are obtained through micro-meso relations, which are derived from damage variables and associated forces and which have an intrinsic meaning, being independent of the stacking sequence. The MSD assumes laminates as a stack of homogeneous plies and interfaces, which are considered as mesoconstituents, whose stiffness loss due to the damage is quantified by the variables ρ and τ. The MSD makes the assumption of linear elasticity, therefore it expresses all the involved quantities as the sum of the solution of a problem P in which damage is removed and the solution of a residual problem P where a residual stress is applied, correcting the undamaged solution around each damaged area. Accordingly, the homogenized potential energy density of the generic ply S is expressed as [ 3 ] are operators that depend on the material properties, |S| is the deformation, while <> + represents the positive part operator. This relation is "quasi exact", since it disregards the coupling effects, because they are always very small in the practical cases. It features an equivalent state of damage on the mesoscale, which is nearly intrinsic for a given state of microdegradation. As far as the generic interface j is concerned, the potential energy density can be written as in [11]: 1 k  1 , k  2 and k  3 are the elastic stiffness coefficients of the interface, I 1 , I 2 and I 3 the three damage indicators and γ j the deformation.
Equations (13) and (14) determine the relationship between the damage indicators and the microdegradation variables ρ and τ of the elementary cell. In fact, they are obtained making the potential energy stored in the plies and in the interfaces the same than in the micromodel. It could be noticed that these equations relate the new elastic properties of the homogenized model to the damage indicators, since the coefficients of the stresses in these relations represent the elastic properties of the homogenized model equivalent to the discretely damaged real model. Therefore a continuum damage model is obtained that is quasi-equivalent from an energy standpoint to the damage micromodel.
Using brick elements to discretize each elementary cell, it is possible to numerically solve the residual problems for determining the damage indicators of (13) and (14). Then, the FD_VK model performs the analysis employing the modified elastic properties given by the mesoscale model, in order to take into consideration the effects of microdamage. Instead of considering the materials healthy as in the unloaded state at the initial iteration, the failure criteria are applied at any time step considering the damage as it is in the reality, since its effects computed at previous iterations are incorporated as modified properties.
In conclusion, the damage growth is described using the initial elastic properties and the stresses that result by this physically based modified expression of the strain energy, instead of simulating damage by guessing appropriate parameters. The stress-based failure criteria used to predict the onset of damage are reported in the Appendix.

The Contact Problem
As explained in the introductory section, indentation represents for structures a serious design concern, since it can determine the growth of barely visible damage, which may produce a large reduction of their mechanical performance, or even failure.
In the present paper, the indentation depth and the contact area are evaluated assuming the distribution of the contact force to be Hertzian, as customarily. The impactor is described as a rigid body and the effective stiffness of the target structure is simulated by the structural model. The iterative algorithm developed by Palazotto et al. [30] is particularized to the FD_VK model, with the aim of forcing the surface of the target to conform the shape of the impactor, as required by soft structures like sandwiches. At any load increment, an iterative procedure that computes the contact radius (equivalent radius for non-spherical indenter) is repeated till the impactor and the indentation radii are in agreement, then the failure analysis is performed.
At the first instant, the contact force F and the vertical displacementare assumed to be zero. After a small time interval Δt, during which the assumption of no damage is postulated, the contact force becomes F=ΔF and, according to the Hertzian law, is distributed over the contact area. As a consequence of the applied load F=ΔF, the impactor moves on at a distance influenced by the effective plate stiffness that depends on the plate stiffness and the contact stiffness. The contact radius determines the applied pressure corresponding to the load and it is computed at every load step using an iterative algorithm that forces the impacted top surface to conform, in the least-squares sense,to the shape of the impactor, applying appropriate constraints. The iterations are repeated till the impactor and the indentation radii are in agreement, then the failure analysis is performed with the damage model described in the previous section. The process is repeated for all the following time steps till the end of the contact time. Since the evaluation of the contact radius at any time step requires very few iterations, this technique gives an accurate description of the contact force and its distribution over the variable contact area not requiring a much larger effort than distributing the contact force over a fixed surface as customarily.
As the solution depends on the current configuration, the previous history and, in addition, deformation can be locally quite large, a modified version of the Newton-Raphson method is adopted to solve the contact problem. The analysis is linearized over any load increment taking into consideration the tangent [K(q)] T and the secant [K(q)] S stiffness matrices. The residual force R i is computed using the secant stiffness matrix, the load at the next iteration F (i) and the solution at the previous iteration q 0 (i-1) . The tangent stiffness matrix is used to evaluate the updated solution that makes the structure in equilibrium from the residual force balance.

Numerical Results
The accuracy of the structural model has been extensively assessed in [25]. In order to show its capability to treat clamped edges, this case being usually considered in the indentation studies, the analysis of a square sandwich panel under transverse distributed loading is premised. After, sample cases of indentation of sandwich panels will be discussed.
As a general remarks it could be underlined that in all the numerical applications the governing equations are solved using the Rayleigh-Ritz's method, assuming trigonometric expansions for the d.o.f., since all the sample cases considered have simple geometry, loading and boundary conditions.
The mass matrix is evaluated without introducing simplifying assumptions, because at certain frequencies the acceleration of the face sheets can be different from that of core, thus the higher order inertia terms do not always provide minor contributions. However, these effects are not relevant in the quasi-static sample cases considered in the numerical applications. In order to demonstrate the advantages of the FD_VK model, results will be compared to those by a formulation equivalent to the layerwise models of Refs. [18] and [19], that is here referred as VD_VK. In this model, the unknowns defined in (4)-(6) are assumed within each computational layer and not for the whole laminate as in FD_VK, while the expansion is limited to the third order in the core and to the first order in the faces. Please note that all the computational times next reported will refer to analyses performed on a laptop computer with dual-core CPU 2.20 GHz, 64 bit operating system and 4 GB RAM.

Assessment of the Structural Model
Clamped edges generally cannot be properly treated by the models having the mid-plane displacements and the shear rotations as functional d.o.f. In fact, these models are unable to feature non-vanishing shear forces at the edges as in the reality, since all the d.o.f. must vanish. On the contrary, the present model does not fail because a non-vanishing transverse shear force can be enforced as a constraint meanwhile the higher-order coefficients of displacements (7)-(9) are computed. So with the present model, the vanishing of the d.o.f. at the clamped edges does not necessarily mean a vanishing shear force, as the elastic solution for clamped edges can be enforced for computing coefficients of displacements.
The sandwich square plate, clamped on all four sides, analysed by Ramtekkar et al. [47] with a 3D mixed finite element model is considered. The material constituting the face sheets has the following mechanical properties: E 1 =90,121 GPa, E 2 = E 3 =7,303 GPa, G 12 = G 13 =4,134 GPa, G 23 =2,68 GPa, υ 12 = υ 13 =0,28, υ 23 =0,34. The mechanical properties of the material constituting the core are: E 1 = E 2 = 0 GPa, E 3 =2,076 GPa, G 12 =0 GPa, G 13 =0,358 GPa, G 23 =0,122 GPa, υ 12 =υ 13 =υ 23 =0. The stacking sequence is [0°/90°/core/0°/90°], the ratio between the thickness of the core and that of the faces is 8, while the length to thickness ratio S is 5. Solution to this case is found in closed form with sin cos  Please note that these trial functions are valid considering a reference system at the centre of the plate. This representation will also be adopted forward for studying indentation. The stresses reported in Figure 3, which are compared to the reference results by Ramtekkar et al. [47], are obtained by the present model in the following normalized form: The symbol p 0 being the intensity of the uniformly distributed load. These results are computed considering a third order expansion of the transverse displacement and a second order expansion of the in-plane displacements. The comparison of the present results with the reference results shows that the present model does not fail in the prediction of deflection and stresses with clamped edges.
Hereafter, numerical results obtained using the simulation procedure previously described are compared to experiments dealing with indentation of sandwiches taken from literature. In all the sample cases, except the one by Castaniè et al. [49], just a second order expansion for the in-plane displacements and a third order expansion for the transverse displacement will be considered, but each physical layer is subdivided into two computational layers, this choice being the best compromise between accuracy and containment of computational cost, as appeared form numerical tests.
Forward, several test cases of indentation of sandwich panels will be solved using the simulation procedure outlined above and the result will be compared to experiments taken from the literature.

Sandwiches with Honeycomb Core
As a first test case, it is considered a sandwich with brass faces of 0,1 mm thickness and Nomex TM honeycomb core (HRH78, 1/8, 3) of 15 mm thickness. Castaniè et al. [49] presented experimental and numerical results of quasi-static indentation tests carried out on 100x100mm fully supported sandwiches. The indentations tests were carried out using spherical indenters with radii of 21,75, 30,125 and 57,25 mm. The numerical results by Castaniè et al. [49] were obtained using a nonlinear finite element model. An elastic modulus of 103100 MPa and a yield stress of 433 MPa were assumed for the brass faces. The honeycomb core was simulated with nonlinear spring, placed at the same location as the honeycomb cell angles and the faces with Mindlin plate element. From Figure 4 it could be noticed that the results by the present simulation procedure are in accordance with the numerical and experimental ones by Castaniè et al. [49]. It could be underlined that the finite element analysis calls for 20 elements across the thickness of the honeycomb, while, as far as the mesoscale model is concerned, the reduction of the mechanical properties of the core is influenced by the indenter radius. Specifically a 30% reduction is necessary, for the smallest indenter, a 33% reduction is required by the indenter with 30,125 mm radius and in the last case a 37% reduction is needed. The computational time is no longer affected by the indenter dimension, since in all the cases the analysis takes about 39 s. The same analyses performed with the VD_VK model require a computational time of 150 s. Castaniè at al. [49] showed also the experimental and numerical results for another sandwich plate with honeycomb core. As in the previous case, the plate is fully supported, it is square with a 100 mm side and the Nomex TM honeycomb core (HRH78, 1/8, 3) is 15 mm thick. The plate has brass faces of 1 mm thickness and according to Castaniè et al. [49], the elastic modulus of the brass is 70400 MPa, its yield stress is 104 MPa. Once again, the indentations are carried out using spherical indenters with radii of 21,75, 30,125 and 57,25 mm. Also in this case, the numerical results by the current simulation procedure are in good agreement with those presented by Castaniè et al. [49], as shown by Figure 5. During the finite element analysis 20 elements across the thickness of the honeycomb are again required, and, also in this case the indenter size has a bearing in determining the results by the mesoscale model. In fact the first indenter calls for a 27% reduction of the mechanical of the core, while 29% and 35% reductions are required in the other cases. As far as computational burden is concerned, nearly 43 s are required in the three cases.
As already mentioned, the numerical results of Figure 4 and 5 by the FD_VK model are obtained using a representation, which is different from that adopted in the following cases. Specifically in-plane displacements and transverse displacement of the faces are represented using a second order expansion, while a third order expansion of the transverse displacement and a second order expansion of in-plane displacements are adopted across the core. A single computational layer is considered for each physical layer. Sandwiches with laminated faces are now considered. The attention is focused on the cases analysed by Williamson and Lagace [50] that presented experimental indentation curves for a rigidly supported sandwich plate indented by a 25,4 mm diameter hemispherical nose. The faces are 3,25 mm thick and they are realized with AS4/3501-6 carbon/epoxy prepregs, whose mechanical properties are E 11 = 142 GPa, E 22 = 9,8 GPa, G 12 = 7,1 GPa, υ 12 =0,3. The 25,4 mm thick core is in HRH 10 1/8-3,0 Nomex TM honeycomb. The comparison between experimental and numerical results is reported in Figure 6, where also the numerical results by Hoo Fatt and Park [51] are considered. They were obtained using the principle of minimum potential energy. Differently to what done in this paper, Hoo Fatt and Park [51] did not consider the contact problem to be Hertzian, as they supposed that the top face sheet is resting on a deformable plastic foundation (an elastic, perfectly plastic core hypothesis). Comparison between experimental and numerical force-indentation curves for a rigidly supported Nomex TM honeycomb sandwich plate indented by a 25,4 mm diameter hemispherical nose.
The finite element analysis needs 30 elements across the thickness of the core, while the mesoscale model calls for a 36% reduction of the mechanical properties of the core and a 41% reduction of those of the faces. As far as computational times are concerned, the analysis takes 50 s.
The results by the present simulation procedure are shown in Figure 7, where they are compared with those by Williamson and Lagace [50] and by Hoo Fatt and Park [51]. In this case the finite element analysis calls for 15 elements across the thickness, while the mesoscale model provides a 38% reduction of the mechanical properties of the core and a 42% reduction of those of the faces. The presence of multilayered faces determines an increase of the computational times, which for this case is 55 s. Also in this case the FD_VK model is advantageous since the analysis by the VD_VK model takes 321 s. Comparison between experimental and numerical force-indentation curves for a rigidly supported Nomex TM honeycomb sandwich plate with multilayered faces indented by a 25,4 mm diameter hemispherical nose.

Foamed Core Sandwiches
As final example, the sandwich panel with foam core analysed by Potluri et al. [52] is considered. The panel is square with a 300 mm side and it is clamped on all the sides. Its faces are in woven glass fabric of weight 600 g/m 2 , while the core is in Divinycell H130. The face thickness is 1 mm and the core is 25 mm thick. Indentation is carried out with a 25 mm flat-faced indenter.
The results presented in Figure 8 are in good agreement with the experiments by Potluri et al. [52]. Please note that in this case 30 solid elements across the thickness have been considered in order to get an evaluation of the crushing behaviour of the Divinycell core. The mesoscale model calls for a 41% reduction of the mechanical properties of the upper face and a 37% reduction of those of the core. The computational time required by this analysis is 45 s. As a general conclusion it could be noticed that in all the sample cases the simulation procedure presented provides results in good agreement with experimental ones. However, at the same time, the computational times reported for each case are very low, since the most intricate analysis takes only 55 s. This achievement is reached primarily thanks to the efficiency of the FD_VK, but also the MSD model and the choice of carrying apart the evaluation of the core crushing behaviour play an important role in this context. These features make the present simulation procedure suited for the analysis of indentation of sandwiches.

Concluding Remarks
Indentation of sandwiches with honeycomb or foam core has been investigated using a recently developed zig-zag model with variable kinematics, together with mesoscale damage model. The paper is aimed at developing a modelling approach able to reduce the computational burden with respect to the typical procedure. Sandwiches are treated as multi-layered composites, whose properties vary with the applied load. The core crushing behaviour is determined apart once at a time using finite element analysis. Once the properties of the core have been determined, the analysis is carried out in closed form using the adaptive model by the authors. It has a variable kinematics, but only five degrees of freedom, in order to keep as low as possible the computational effort. The damage formation is determined using stress-based criteria, as customarily, while the damage growth is accurately and efficiently taken into consideration employing a mesoscale model. The present simulation procedure provides results in good agreement with experiments, in all the cases examined. At the same time the computational effort is kept also very low because the evaluation of the variable elastic properties of the core is carried out apart once at a time. For instance, the analysis of an intricate sandwich panels with 32 layers takes only 55 s, while the same analysis performed with a model with variable d.o.f. requires 321 s. This difference is due to the fact that the structural model here presented has a fixed number of d.o.f., while the customarily adopted models have a number of d.o.f. that varies with the number of layers. To conclude, the present simulation procedure is proven to be suited for analysis of indentation of sandwich structures. 1 1 , ( 0  The onset of delamination is predicted using the Choi-Chang's criterion, since as shown by many researchers it gives accurate estimations. Delamination is predicted to start if where Y n+1 = Y t n+1 if σ yy ≥0, or Y n+1 = Y c n+1 if σ yy <0, D a is an empirical constant that is defined after consideration of the material properties, σ ij , which is the average stress at the interface between the nth ply and the n+1th ply, is computed as follows: The subscript 'i' stands for in-situ, while 't' and 'c' stand for traction and compression, respectively. This criterion does not take into consideration the transverse interlaminar stress σ 33 . However numerical tests in literature have shown that, in the majority of cases, this component has not a significant bearing.
The criterion by Besant et al. [45] is chosen to predict the failure of sandwiches with honeycomb core under compression and transverse shear. Failure occurs when: σ cu and τ lu are the core strengths in compression and transverse shear. In this paper, it is chosen n= 1,5, because this value best fitted the experimental results. Nevertheless numerical tests in literature have shown that varying the exponent from 1 to 2 does not determine any remarkable effects on the results of sandwiches with laminated faces.
The rule suggested by Evonik for Rohacell core is used to predict the failure of foam cores: Using as safety factor R 11 /σ v . The empiric parameters, a 1 =k 2 (d-1)/d, d= R -11 / R + 11 , a 1 =(k 2 /d)-1 and k= 3 R 12 /R + 11 are determined from experiments, R + 11 , Rand R 12 are the strengths of the foam under traction compression and shear, respectively.
The criteria by Lee and Tsotsis [43] and Petras and Sutcliffe [44] are adopted to estimate failure indentation due to core crushing under out-of-plane normal and shear stresses. The former predicts indentation failure to occur when one of these inequalities is verified: