Impact Damage Analysis with Stress Continuity Constraints Fulfilment at Damaged-Undamaged Regions and at Layer Interfaces

Low-velocity impacts have a relevant importance for safety of laminated and sandwich composite structures, because they are highly susceptible to damage. In this paper, a 3D cost effective zig-zag model is developed in order to efficiently simulate such impacts. It a priori fulfills the continuity of out-of-plane stresses at layer interfaces, the continuity of stresses under in-plane variation of properties across undamaged and damaged regions and it is suitable for general boundary conditions. Its main advantage is its capability to accurately predict stresses from constitutive equations at a low cost, along with being refined across the thickness keeping fixed the number of unknowns. A modified Hertzian contact law that forces the target to conform to the shape of the impactor and the Newmark’s implicit time integration scheme are used for computing the contact force time history. The progressive damage analysis is carried out using stress-based criteria. Non-classical feature, a continuum damage mesomechanic model is employed that accounts for the effects of local failures in homogenized form by modifying the strain energy expression. Fiber, matrix and de-lamination failures are predicted using stress-based criteria, and then the modified strain energy expression is employed for computing stresses. Such modeling options enable to account for the residual properties as they are in the reality, as shown by the comparison with the damage experimentally detected. As shown by the comparison with experiments, a closed-form solution by the Ga-lerkin’s method, obtained as a series expansion of trial displacement functions, accurately simulates the contact force and the damage progressively accumulated. The results show the importance of in-plane stress continuity for obtaining accurate predictions.


INTRODUCTION
Laminated and sandwich composites are increasingly used as primary structural components by virtue of superior specific elastic moduli, strength, energy absorption and dissipation properties.Higher speed, longer range, larger payloads, less engine power and better operating economy can be achieved using these materials.
Owing to elastic moduli and strengths in the thickness direction that are much smaller than their in-plane counterparts, warping, shearing and straining deformations of the normal and out-ofplane stresses rise that are responsible of local micro-failures and, consequently, of relevant loss of strength and stiffness.As the damaged area may increase in size under loading in service, till causing a premature failure, the basic requirement is that composites must carry the ultimate design load even when they are damaged.So, accurate numerical simulations are required to understand the effects of the damage accumulated in service.Manufacturers have found a high susceptibility to impact damage from dropped tools, vehicular collisions and to damage caused by drilling.A critical case for aircraft structures is impact from crushed stones during taxiing manoeuvres and take-off.
The readers are referred to the papers by Richardson and Wisheart (1999), Morinière et al. (2014) and by Abrate (2001Abrate ( , 2014) ) for a comprehensive review of low-velocity impact studies on laminates.Methods for predicting delamination induced by low-velocity impacts are discussed by Elder et al. (2004).A review of numerical and experimental impact studies on sandwich panels, a discussion of their damage mechanisms and of failure criteria used in the simulations are given by Chai and Zhu (2011).An overview of damage models is also given by Garnich and Akula Venkata (2009), Liu and Zheng (2009) and Berthelot (2003).The damage analysis of honeycomb sandwich panels carried out by Horrigan and Staal (2011) is also cited as a relevant contribution in the field.A comparison of simulations and experiments is given by Hongkarnjanakul et al. (2013).The effects of stacking sequence on the impact and post-impact behaviour is shown by Aktas et al. (2013), those of the impactor shape are studied by Mitrevski et al. (2005).Multiple impacts are investigated by Damanpack et al. (2013) and Chakraborty (2007).The impact behaviour of different honeycomb core cells is shown by Yamashita and Gotoh (2005).
The structural models play a central role in the simulation procedure, since they should accurately describe the behaviour of composites requiring a low computational effort, as many iterations are required to solve the impact problem.An extensive discussion of structural models used in the simulations and extensive assessments of their performances are given by Chakrabarti et al. (2011), Chen and Wu (2008), Kreja (2001), Zhang and Yang (2009), Tahani (2007) and Matsunaga (2004).
In a broad outline, three-dimensional finite element discretization (3D-FEA), refined plate models (R-M) and related plate elements are of interest for damage analysis.R-M can be further split into: Discrete-Layer models (DL) with a separate representation across each computational/physical layer, High-Order Layerwise models (HLW) based on a combination of global higher-order terms and local layerwise functions and physically-based zig-zag models (ZZ-F) with fixed degrees of freedom (d.o.f.) that a priori fulfil the displacement and stress continuity conditions at the material interfaces.
3D-FEA, DL and HLW are inadequate for impact damage analysis of structures of industrial interest, as too large number of unknowns may overwhelm the computational capacity.Studies by Palazotto et al. (2000a), Kärger et al. (2007), Diaz Diaz et al. (2007), Icardi and Ferrero (2009) and Latin American Journal of Solids and Structures 14 (2017) 1416-1442 Oñate et al. (2012aOñate et al. ( , 2012b) ) have shown that ZZ-F and even equivalent single layer models (ESL) can be can be successfully used in order to save costs.
Recent studies have shown that stress predictions as accurate as those by 3D-FEA and DL can be obtained using global-local models, as shown by Zhen andWanji (2006, 2010), Wanji et al. (2007), Zhen et al. (2009), as well as incorporating zig-zag functions based either upon kinematic assumptions, e.g.Murakami (1986), or physical assumptions, e.g.Tessler et al. (2009Tessler et al. ( , 2010) ) and Iurlaro et al. (2014).Accurate results have been obtained by Icardi andSola (2014a, 2014b) using ZZ-F of the latter type, with fixed d.o.f. and variable kinematic, with a processing time and memory storage occupation comparable to those minimal of ESL.Although these features make such models suitable for low-velocity impact studies, no applications have been presented yet in literature.
To contribute to fill the existing gap, in this paper a low velocity impact analysis is carried out using as structural model a refined version of the variable kinematic zig-zag model (VK-ZZ) with fixed d.o.f.developed by Icardi and Sola (2014a), having as functional d.o.f. the two in-plane displacement components, the transverse one and the two shear rotations at the reference middleplane.So, the characteristic feature of this model is that its number of unknowns neither increases with the number of physical or computational layers (i.e.refining the representation across the thickness), nor with the number of constraints enforced.The present refined version of the VK-ZZ model makes continuous the stresses at the interfaces of undamaged and damaged regions, where material properties suddenly vary, not just across the thickness, i.e. at layer interfaces as in the former version.
Zig-zag functions making piecewise variable displacements are incorporated in order to a priori fulfil the 3D stress contact conditions at the interfaces of physical/computational layers and of undamaged and damaged regions, as prescribed by the elasticity theory for satisfying equilibrium.The assumed displacements a priori satisfy the stress boundary conditions at the upper and lower bounding surfaces and general boundary conditions at the edges as well.The fulfilment of local equilibrium equations is enforced at selected points across the thickness, in order to define the coefficients of high-order contributions of displacements that enable a variable representation from point to point.
By virtue of an accurate representation of strain energy density, displacement, strain and stress fields obtained through a realistic representation of warping, shearing and straining deformations, the present model efficiently computes the stress field at each time step from constitutive equations, even when the constituent layers have distinctly different properties.
To speed-up computations and to overcome algebraic difficulties, closed form expressions of its coefficients and zig.zag amplitudes are obtained once for all in terms of the 5 functional d.o.f. and of their derivatives using symbolic calculus.A C 0 formulation of the model useful for developing efficient finite elements free from nodal d.o.f.derivatives can be obtained from the present VK-ZZ model applying the technique as shown in Icardi and Sola (2014c).However, in the present paper this technique is left apart.
The progressive damage analysis is carried out at each time step using stress-based criteria, extending a pre-existing damaged area at the previous iteration to the points where the ultimate condition is reached at the next one, as customarily.Instead of guessing the residual properties of the Latin American Journal of Solids and Structures 14 (2017) 1416-1442 damaged area, they are determined using the continuum damage mesomechanic model by Ladevèze et al. (2006) that accounts for the effects of hard discontinuities in homogenized form through a modified version of the strain energy.The transverse cracking rate and delamination ratios are computed using the stress-based criteria at each step, then the homogenized energy that accounts for the damage is computed and used for evaluating stresses at each point.
The modified Hertzian contact law of Icardi and Ferrero (2009) is adopted to compute the contact force for laminates.The refined contact model by Palazotto et al. (2000b) that forces the target to conform to the shape of the impactor is employed to compute the contact radius at each time step when sandwiches are analysed.
Hertzian contact is chosen being proven to be accurate by a plenty of studies (see, e.g.Yigit and Christoforou (1995) for laminates and Choi (2006) for sandwiches), nevertheless it was developed for isotropic media.A different alternative approach accounting for the heterogeneous nature of composites is that developed by Chao and Tu (1999).In this case, the contact stresses are derived via an energy variational approach, accounting for three-dimensional boundary conditions and interlaminar continuity.Another alternative approach is developed by Zhou and Stronge (2006), where the load-indentation law of sandwiches is obtained via FEA (ABAQUS/Explicit) assuming the elastic face sheets and rigid-perfectly plastic the core.Application of these more advanced contact models is left to a future study.The Newmark's implicit time integration scheme is used in the present paper for solving transient dynamic equations, despite experiments have shown that low-velocity impact studies could be carried out also in static form (see, Li et al. (2012)).One practical reason for carrying out a transient analysis is that dynamic effects can be captured progressively increasing the impactor velocity, thus a unique scheme can be used for a wide range of applications.Another reason is that the computational cost of Newmark's implicit time integration scheme is not much higher of that required by a static analysis carried out at each load step in order to capture the evolving damage.
Differently to paper by Icardi and Ferrero (2009), where non-linear strains of von Karman type are adopted, here the updated Lagrangian approach is chosen to efficiently account for geometric nonlinearity.Accordingly, the deformations at each new time step are obtained from the geometry of the structure at the previous step, not from the initial unloaded configuration.
Since core crushing is governed by the properties of the cellular structure of honeycomb, it cannot be appropriately accounted for by the present multi-layered VK-ZZ model that describes the core as a thick layer with homogenized properties.Therefore, the variable, apparent elastic moduli of the core while it collapses/buckles under transverse compressive loading are computed apart once for all through a detailed three-dimensional finite element analysis and provided to the VK-ZZ as variable properties at each load level, like in Icardi and Sola (2015).
Numerical applications presented are aimed at assessing whether the in-plane stress continuity has a relevant bearing on results.To this purpose, the results obtained in closed form adopting a truncated series expansion of trial functions for the d.o.f. with and without in-plane continuity enforcement are compared.The sample case considered are those already investigated in Icardi and Ferrero (2009) or which experimental results giving the contact force time history and the damage detected using ultrasonic cartography in Icardi and Zardo (2005) are available.The comparison of the present results with those obtained by Icardi and Ferrero (2009) will also show the progress Latin American Journal of Solids and Structures 14 (2017) 1416-1442 obtained considering a damage mesomechanic model that accounts for the effects of hard discontinuities in homogenized form.

CONTACT FORCE MODELLING
Since the literature shows that the contact pressure profile is less important than the load and the area over which it acts, the contact stress is assumed to be distributed according to the Hertzian law, like in homogeneous isotropic materials as: In the former equation, contact R is the radius of the contact area, while ( ) r s , (0) s are the Hertzian stress intensity at a distance r from the center and at the centre, respectively.The impactor is assumed to be spherical, as no details regarding the real shape of the contact area and the real distribution of the force acting on it are available in the literature for other impactor shapes.The analysis of laminates is carried out assuming contact R as a fixed variable, while for sandwiches it is computed at each step by forcing the target to conform to the shape of the impactor.

Laminates
In the loading phase, the contact force F is related to the indentation depth α through the contact stiffness Kc in standard form as: A 3D Finite element analysis, here referred as PFEA-1 (see Figure 1) is performed apart once for all as a virtual material test (VMT) for determining Kc and Rcontact.A non-linear model (RA-DIOSS -MAT25) is used to account for the failures occurring in the target structure, as described in section 3.
The exponent υ is set in agreement with experiments in the literature.Also in the unloading phase the exponent q is set according to experiments.In the previous equation, Fm is the maximum of the contact force, αm is the relative indentation depth at each loading and α0 is the permanent indentation.In case of bounces, the following reloading law is assumed: Again, the contact stiffness k is computed by VMT accounting for contact stiffness degradation due to the damage arose in the loading phase.The exponent p is still set according to experiments.Structures 14 (2017) 1416-1442 The literature shows that υ= p= 1.5 and q= 2.5 is very often the most appropriate choice for laminates.

Sandwiches
A constant R contact being inappropriate for simulating sandwich panels with lightweight soft core, the iterative algorithm by Palazotto et al. (2000b) is applied to compute the contact radius, as it forces the impact area to conform to the shape of the impactor at any load increment.The problem is solved using a modified version of the Newton-Raphson method: ( ) because the solution depends both on the current configuration, the previous history and local deformations can be quite large.The secant stiffness matrix [K(d)]S is computed using the VK-ZZ model of section 4. The symbols d (i-1) and δ i represent the converged solution at the previous load increment Fc (i-1) and the residual force, respectively.The displacement amplitude vector d (i-1) is updated by Δd i in order to make the structure in equilibrium with Fc i .Displacements are applied during iterations till the impactor and the indentation radii conform: ( ) being the tangent stiffness matrix computed using the VK-ZZ model.The updated vector of displacement d.o.f.d i is computed as: Latin American Journal of Solids and Structures 14 (2017) 1416-1442 The process is repeated till the convergence tolerance is reached.The convergence tolerance is confronted with the percentage of variation of the solution, from the current to the previous iteration using the norm of displacements ( ) to compute The iterative process is stopped when D ≤ 1 %.Once a converged solution is obtained, the failure analysis is carried out.

DAMAGE AND POST-DAMAGE MODELLING
Stress-based criteria with a separate description of various failure modes are used to estimate where the constituent materials are failed.

Failure Criteria
The 3-D criterion developed by Hashin and Rotem (1973) with in-situ strengths is used to predict the fibre and matrix failures.Tensile failure is ruled by: ( ) while under compression it is described by: Latin American Journal of Solids and Structures 14 (2017) 1416-1442 Delamination failure is predicted using the criterion developed by Choi and Chang (1992): where set after consideration of the material properties, ij s is the average stress at the interface between the n th ply and the n+1 th ply: The symbol 'i' stands for in situ, while 't' and 'c' stand for traction and compression, respectively.This criterion is chosen because numerical tests in literature have shown its accuracy for laminates undergoing low velocity impact loading.
The failure of honeycomb core under compression and transverse shear is predicted using the criterion by Besant et al. (2001): ( 1) cu s and lu t being the core strengths in compression and transverse shear.It is assumed n = 1.5, since numerical tests in literature have shown that results are quite accurate.However, results do not considerably vary changing the exponent from 1 to 2. The failure of foam core is estimated according the rule suggested by Evonik for Rohacell core (see Rohacell (2011) and Li et al. (2000)), using as safety factor 11 / v R s : where are determined from experiments, while 11 R + , 11 Rand 12 R represent the strengths of the foam under traction compression and shear, respectively.Core crushing failure under out-of-plane normal and shear stresses is also predicted using the criterion by Lee and Tsotsis (2000): Latin American Journal of Solids and Structures 14 (2017) 1416-1442 c Z , x S , y S being the compressive yield strength and the out-of-plane shear strengths, respectively.As both criteria ( 16) and ( 18) refer to honeycomb failure, the failed region is computed as the envelope of failures predicted by each of these two criteria.Strengths and material properties assumed in section 5 are taken from material databases and from experiments published in the literature.

Mesoscale Damage Model
In order to account for the effects of hard discontinuities in homogenized form, the mesoscale damage model (DML) by Ladevèze et al. (2006) is used that provides a modified strain energy expression from which stresses are computed.With this model, no factors are guessed to modify the initially health properties in the failed regions as the damage is simulated as it is in the reality inside energy functionals.
DML replaces the discretely damaged medium containing matrix cracks, delaminations and fibre failures with a continuous homogeneous medium, equivalent from an energy standpoint, whose strain energy expression incorporates damage indicators computed as the homogenized result of damage micromodels.Micro-meso relations, representing the homogenized result of micromodels, are derived from damage variables and associated forces.Displacement, strain and stress fields on the microscale are denoted as U m , ε m and σ m , respectively.These quantities are expressed 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 around the damaged area.The homogenized potential energy density of the generic ply is expressed as Ladevèze et al. (2006): The meaning of symbols is as follows: 22 12 23 13 , , , I I I I and 33 I are the damage indicators, each defined as the integral of the strain energy of the elementary cell for each basic residual problem under the five possible elementary loads in the directions 22, 12, 23, 13, 33; operators that depend on the material properties, S is the deformation, .+ < > represents the positive part operator.The homogenized potential energy density of the generic interface is expressed as Ladevèze et al. (2006): which are damage variables that links the damage state to the loss of stiffness of the damage mesoconstituent, L being the distance between two adjacent cracks, H being the length of the crack across the thickness and e being the length of the microcrack.Experiments show that the range of practical interest for r is [0; 0.7] and for t is [0; 0.4].
In this paper, the damage indicators are determined numerically through 3-D FEA VMT (PFEA-3).The elementary cell is discretized by 3-D elements once for all outside the time-loop cycle, considering discrete values of fibre failure and matrix cracking rates and delamination ratios of the elementary cell obtained at various load levels applying criteria (10) to ( 14).
The modified elastic properties by the mesoscale model are provided to the VK-ZZ model that is used to carry out the analysis.The progressive failure analysis is performed at each time step applying criteria ( 10)-( 18) under the forces provided by the contact model ( 1)-( 9) and the modified elastic coefficients of the mesoscale model ( 19)-( 20).The damage computed at the previous step by the criteria (10)-( 18) is extended at the next step to the points where the ultimate stress is reached.Instead of carrying out the analysis at an infinite number of points, a discrete representation is used that identifies small square damaged cells of the domain over which the ultimate stress is reached or just passed (up to an assumed magnitude tolerance, see Fig. 9), in a fashion similar to that of finite elements.Accordingly, the domain is subdivided into fictitious small square cells, whose damage state computed at the central point is assumed constant over the area of the cell.A layer loop computes the stresses across the thickness at the central point of each cell, which is used to carry out the damage analysis.Obviously, if the dimension of cell decreases, the accuracy and computational effort rise.

THE VK-ZZ STRUCTURAL MODEL
The VK-ZZ model, as developed by Icardi and Sola (2014b), is chosen as structural model.The displacement field is postulated to piecewise vary across the thickness as the superposition of continuous functions and layerwise functions having appropriate discontinuous derivatives at the layer interfaces.
Here x, y represent the in-plane coordinates (the mid-plane is assumed as reference plane) and z the thickness coordinate (-h/2 <z< h/2), while u, v and w denote the elastic displacement components in the directions x, y and z, respectively.
The through-the-thickness variation of elastic displacements is assumed as follows: The explicit expressions of each contribution are given forward.In order to determine the structural response through a closed form approach, they are expressed as truncated series expansions of unknown amplitudes computed by solving the transient dynamic problem and assumed in-plane functions.

Linear Contribution
It provides the basic linear contribution across the thickness and contains the five functional d.o.f.
( , , ) ( , ) x x y y U x y z u x y z x y w x y V x y z v x y z x y w x y W x y z w x y u 0 (x, y), v 0 (x, y), w 0 (x, y), γx 0 (x, y) and γy 0 (x, y) respectively represent the two in-plane elastic displacement components of the points laying on the reference middle-plane, the transverse component and the two shear rotations at these points.

High-Order Contribution
They are aimed at refining the model across the thickness in the regions with step stress gradients, because they can vary the representation from point to point across the thickness: ( , , ) ...
The expressions of 1 x A … zn A are determined by enforcing the stress-boundary conditions at the upper and lower faces of laminated and sandwich plates: and the fulfilment of local differential equilibrium equations at discrete points chosen across the thickness: So the expansion order can be increased adding more points.Please notice that the computations of Ax1 … Azn don't imply adding new d.o.f.Symbols bi represent volume loading (e.g., inertia terms), while as customarily, a comma is used to indicate differentiation (for instance σzz,z represents the stress gradient of σzz across the thickness).

Piecewise Contribution
These terms are aimed at a priori fulfilling the continuity of displacements and interlaminar stresses at the material interfaces through appropriate discontinuous derivatives across the thickness: Hk being the Heaviside unit step function (Di Sciuva (1985)) , i.e.Hk=0 for z < zk, while Hk=1 for z ≥zk.The contributions to in-plane displacements Φx k , Φy k are similar to those postulated by Di Sciuva (1985) in order to fulfil the continuity of σxz and σyz at the interfaces between layers.The contribution Ψ k to the transverse displacement is aimed at satisfying the continuity of σzz, while Ω k satisfies the continuity of the gradient σzz,z (see, Icardi and Sola (2014a)), as prescribed by the elasticity theory (it directly follows from local equilibrium equations ( 25) and from conditions ( 24)).
The contributions Cu k , Cv k and Cw k restore the continuity of displacements where the representation is changed across the thickness through (23).In this way, all stress and displacement continuity conditions are fulfilled.Summing up, the continuity functions are computed by enforcing the following set of stress continuity conditions: Algebraic difficulties are overcome using symbolic calculus (MATLAB® symbolic software package).As a result, expressions for Φx k , Φy k , Ψ k , Ω k , Cu k , Cv k and Cw k are computed once for all for any lay-up that speeds up computations.SEUPT technique, as shown in Icardi and Sola (2014c) can be applied for obtaining a C 0 formulation that overcomes the presence of unwise derivatives of the d.o.f. in the expressions of continuity functions, so to develop efficient finite elements.

Variable Piecewise In-Plane Representation
These contributions are aimed at making continuous the stresses under in-plane variation of properties, namely crossing the interface between undamaged and damaged regions.They constitute the new contribution brought by the present paper.
New zig-zag functions of the following form are added, in order to restore the continuity of stresses under a discrete in-plane variation of properties: Latin American Journal of Solids and Structures 14 (2017) 1416-1442 The continuity functions j k u x q make continuous xx s in the x direction, while their counterparts q do the same in y.
q are aimed at making continuous yy s in the y direction, while j k v x q make it continuous in x.In a similar way, å makes continuous the stresses across the thickness.Whether desired, the continuity of the gradients of order n of these stresses could be enforced considering power of order (n-1).Instead the in-plane continuity of σzz and of its gradient should not be enforced when in-plane discontinuity occurs.
A former application of the idea discussed in this section was presented in Icardi and Sola (2015) for studying bonded joints (suddenly changing lay-up).In that case, the material variation was allowed only in one direction.Therefore (28) constitute a generalization of Icardi and Sola (2015).
Also the continuity functions in (28) are computed using symbolic calculus through enforcement of the continuity of the chosen stresses.It could be noticed that just a fraction of second is required on a laptop computer, once their closed form expressions are obtained for the first time through symbolic calculus in a couple of minutes.
As previously shown by Icardi andSola (2014a, 2014b) the apparently cumbersome and time consuming zig-zag model ( 21)-( 28) requires an overall computational effort for the analysis that is still comparable to that of the basic model as used by Icardi and Ferrero (2009) and by Icardi and Sola (2016) with stresses computed by integrating equilibrium equations.
A comparison of simulations with and without contributions (28) will be presented in order to show the effect brought by the in-plane continuity of stresses.A case previously studied in Icardi and Zardo (2005) with a less refined version of the simulation procedure will be retaken in order to show the progress obtained by refining the structural model and using the DML for the damage analysis.

Core Crushing Simulation
With the VK-ZZ model, sandwiches are described in homogenized form as multilayered structures.In order to account for the core crushing mechanism, the technique developed by Icardi and Sola (2015) that accounts for the buckling of cell wall of honeycomb core is adopted.
A detailed finite element analysis (PFEA-2) is carried apart once for all in order to compute the apparent elastic moduli of the core while it collapses/buckles under transverse compressive loading.They are determined as the tangent moduli at each load level from the average ratio of stresses and strains, then they are provided to the VK-ZZ model for the analysis.
In PFEA-2, a detailed representation of the honeycomb structure as it is in the reality is employed.A rather refined meshing is used to discretize the cellular structure of core, because the prediction of micro-buckling phenomena requires a dimension of shell elements comparable to the wall thickness.An elastic-plastic isotropic material and a self-contact algorithm are used to prevent from interpenetration between the folds in the cell walls.The cell walls bonded together are modelled using double shell elements.Criteria ( 16) and ( 18) reported previously are used to have an indication of the load range at which carrying out the analysis PFEA-2.Foam core is discretized by solid elements with nonlinear material behaviour derived from experiments and materials databases.
In Icardi and Sola (2015), this technique was shown to provide accurate predictions for a variety of sample cases dealing with sandwich indentation whose experimental results are published in the literature.
As PFEA-2 is carried apart once for all in order to determine the apparent elastic moduli of the core at each magnitude of transverse loading, the simulation of the core crushing mechanism while it collapses/buckles under the contact load is accounted for by the VK-ZZ model with a low computational effort at each time step and quite accurately, as it will be shown by the numerical results.

Solution of the Structural Model
Solution is found in closed form at each time step within the framework of the Galerkin's method, assuming the functional degrees of freedom u 0 , v 0 , w 0 , γx 0 and γy 0 expressed as: for clamped edges, or as: Latin American Journal of Solids and Structures 14 (2017) 1416-1442 for edges that are freely movable in the in-plane directions.At clamped edges, the mechanical boundary conditions are satisfied by an appropriate choice of terms ( 23), Γ being the contour of the plate.Qx, Qy are assumed to be uniformly distributed over Γ when the edges are sufficiently far from the impact point located at the centre of the plate.This assumption is valid for large length side-to-thickness ratios.
The dynamic governing equations in matrix form are obtained from ( 29) or (29'), then they are solved using Newmark's time integration scheme.Once unknown amplitudes Amn, …, Emn are computed at each time step, stresses are calculated and the damage analysis is carried out.
The expressions of linear, secant and tangent stiffness matrices of the VK-ZZ model are deduced from ( 29) or (29'), once they are substituted in ( 22) to (28) according to Zienkiewicz and Taylor (2000).The consistent mass matrix of the VK-ZZ is computed considering the inertia contributions within governing functionals, like the dynamic version of the principle of virtual work.The consistent mass matrix is used in the computations, because numerical tests using the mass matrix derived just considering displacements ( 22) have shown errors up to 10% in the first five eigenfrequencies, with just a reduction of the processing time less than 15%.Hence, the use of a simplified matrix is left to a future study

NUMERICAL RESULTS AND DISCUSSION
The contact force time history for laminates and sandwiches and the damage predicted by simulations are compared to experiments and to the results of previously developed zig-zag models and finite elements.The results of simulations are presented for different impact energies in order to assess: (i) which advantages are brought by the refined VK-ZZ model by the viewpoint of progressive damage analysis; (ii) whether consideration of the effects of hard discontinuities (matrix cracking, broken fibres and delamination) in homogenized form improves accuracy; (iii) whether consideration of in-plane stress continuity has a significant bearing for the quality of results and, finally; (iv) whether the impact damage simulation can be successfully carried out in closed form.
As no damage detections are available for sandwiches with laminated faces, except for a stitched one (see Potluri et al. (2003)), the comparison with experiments is limited to the contact force time history.Modelling of stitching (Icardi and Sola (2013)) could be applied, but as it represents an uncertainty source that could prevent from the assessments of (i)-(iv), it is left to a future study.

In-Plane Stress Continuity
In order to assess the capability of the VK-ZZ structural model to predict continuous stresses under sudden in-plane variation of properties (Figure 2), first the sample case of a two-material wedge is considered (Figure 3).
The goal is to compare the in-plane variation of the shear stress σxy across the material discontinuity to exact elasticity solutions by Hein and Erdogan (1971) and Yamada and Okumura (1983).Structures 14 (2017) 1416-1442 In this sample case, the variation of solutions in the transverse direction with respect to the plane x, y of Figure 3 has no relevance.Here just the in-plane shear stress is considered, being the one with the highest singularity strength at the corner interface.A Young's modulus of 7.3 GPa and a Poisson's ratio of 0.3 are assumed for the deformable sector, while a modulus larger by a factor 100 is assumed for simulating the rigid sector.The exact elasticity solution shows a strong stress concentration at the edge corners, which could produce a catastrophic failure in service when dissimilar materials are bonded or welted together.Hence the stress analysis should be carried out with the maximal accuracy at the interface Latin American Journal of Solids and Structures 14 (2017) 1416-1442 of materials with dissimilar properties.Figure 3 reports the in-plane variation of the shear stress σxy across the material interface normalized as follows:

Latin American Journal of Solids and
The present solution is obtained assuming a subdivision of the domain of the two half plates into 16x16 subdomains whose side length is progressively reduced close to the corner interface.Within each subdomain, the variation of the functional d.o.f. of the VK-ZZ model u 0 , v 0 , w 0 , γx 0 and γy 0 is assumed to be linear (product of Lagrange's polynomials in x and y).Hence, the continuity of the in-plane shear stress is enforced at the sides of subdomains by computing the appropriate expressions of the continuity functions j k u x l and j k v y l appearing in Eq. ( 28).
In this sample case, it is sufficient to enforce the continuity of σxy across the material discontinuity line.For multi-layered materials, i.e. with different properties across the thickness direction z, the stress continuity should be simultaneously enforced at each interface through Eqs. ( 27) and ( 28) and at each point in x, y.Thus, this continuity should be enforced at the interfaces of undamaged and damaged regions and for all stress components excepted σzz.
The results of Figure 3 show the capability of the VK-ZZ model to accurately reproduce the variation of the shear stress σxy across the interface even in this case.The present structural model is shown able to capture a continuous stress with a discontinuous gradient like in the exact solution.
The strong stress concentration occurring at the interface of material with dissimilar properties shows the importance of enforcing the stress continuity under in-plane material discontinuity, as in impact damaged structures, to accurately evaluate stresses.The effective importance of in-plane stress continuity is further assessed in section 5.4 in a specific case for which the impact-induced damage is detected by ultrasonic cartography.
In the next sections 5.2 and 5.3 the contact forces for laminates and sandwiches predicted at different levels of energy are compared to experiments.

Contact Force of Laminates
First, the sample case of a 230x127 laminated plate with a [45°/0°/-45°/90°]2s lamination scheme and with all the edges clamped is presented.The plate is subjected to the impact of a steel hemispherical nose of 1.15 kg at different levels of energy (1.5 J and 2.5 J).Each of its 16 layers is made of AS4/3502 Graphite-Epoxy and its total thickness is 2.48 ± 0.0254 mm.
The acceleration time-history is detected through an accelerometer mounted on the impactor nose and connected to a data acquisition system.The contact force time-history is obtained multiplying the impactor mass by the measured acceleration.Since the arm carrying the impactor nose follows a circular trajectory, attention has been paid to assure that the impact force is perpendicular to the impacted panel, as assumed in numerical simulations.A sampling frequency of 10 5 Hz is used.
The results for this case by the present simulation presented in Figures 4 and 5    The results of Figures 4 and 5 show the numerical simulation to be appropriate for this sample case, its results being in good agreement with experiments.Note that just 350 s are required to carry out the analysis on a laptop computer, thus the efficiency of the present simulation is proven.The results with and without enforcement of the target to conform the shape of the impactor are shown to be equivalent for this case, as the iterative algorithm of section 2.2 is unnecessary when laminates are analysed.
In the subsequent sections, further tests are made considering other samples, with the aim of showing if accurate predictions can be always obtained by the present simulation technique.

Contact Force of Sandwiches
The experimental contact force time history detected by Kärger et al. (2007) is compared to that predicted by the present simulation, with and without enforcing the panel to conform to the shape of the impactor and by the previous model by Icardi and Ferrero (2009).The aim is to show that soft media requires the application of algorithm of section 2.2, as well as to show that accuracy is improved with respect to the previous structural model by Icardi and Ferrero (2009).
The thin top face sheet (0.633 mm) is made of Cytec 977-2/ HTA.The thick bottom face sheet (2.7 mm) is made of CFRP fabric plies and UD tapes, with a [(0°2/45°/90°/-45°)s]s lay-up.The core, which is 28 mm thick, is made of NOMEX 4.8-48 honeycomb.This sandwich panel is completely supported at the edges and it is impacted by a steel sphere with energy of 8 J, having a diameter of 25.4 mm carried by an impactor with a mass of 1.10 kg.The results of simulations reported in Figure 6 show that a quite accurate prediction of the contact force is obtained at each instant of time only if: (i) the iterative algorithm of Section 2.2 that computes the contact radius at each step is considered; (ii) the detailed FEA representation of the core behaviour is carried out through PFEA-2, in order to account for the buckling of cell walls under transverse compressive loading, otherwise totally wrong results are obtained.
The results also show the superior predictive capability of the present model with respect to the former model by Icardi and Ferrero (2009), which already considered the core crushing mechanism.Hence the present VK-ZZ model has a significant bearing on accuracy of results.
As in this case the algorithm that accounts for the core crushing behaviour is activated, the computations at each time step are a little bit slower than for the laminates previously considered.However, just 550 s are required to compute the contact force time history, hence the simulation procedure is still proven to be accurate and cost-effective.
Whether accurate predictions of the contact force of sandwiches can be obtained at different levels of energy is discussed next.
Attention is now focused on the sandwich panel with PVC foam core (Divinycell H250) and laminated carbon fabric/epoxy faces (AGP370-5H/ 3501-6S) studied by Schubel et al. (2005), for impact energies ranging from 7.75 to 108 J.A sphere with a radius of 12.7 cm mounted on an impactor with a mass of 6.22 kg is left to fall onto the panel.Increasing the heights of the mass drop, the above mentioned impact energies are obtained, with velocities ranging from 1.6 to 5 m/s.The sandwich plate has sides of 27.9 cm, thickness of 2.82 cm and it is simply supported.The core is 25.4 mm thick, while the laminated face sheets are obtained stacking four plies of prepreg for a total thickness of 1.37 mm.According to Schubel et al. (2005), the mechanical properties of the faces are: ρ= 1600 kg/m 3 , E1= 77 GPa, E2= 75 GPa, G12= 6.5 GPa, G13= 5.1 GPa, G23= 4.1 GPa, υ12 = 0.07.
Since at low energies the damage may occur with a different mechanism than at high energies, consideration of different impact energies is required to assess whether the range of applicability of the present model is quite large to cover practical applications.However, only cases for which strain-rates effects are unimportant will be considered Figure 7 shows the comparison between the experiments of Schubel et al. (2005) and the contact force time histories predicted by the present simulations at the three energy levels of 7.75, 41.1 and 108 J.The results show that the present simulation obtains accurate contact force time histories in all three cases.No increasing errors are seen increasing the impact energy.Thus, it is reputed that the present modelling approach can be successfully used at least for energies up to a hundred of J.A further refinement is required in order to more closely represent the jumps in the force and the oscillating behaviour detected during experiments.
In the next section, the damage predicted by the present simulations is compared to that detected via ultrasonic cartography.
For this case the results by various zig-zag models by Icardi and co-workers are available for comparisons, in addition to experimental results.The comparison with the former models will permit to understand the improvements achieved over the years.
Latin American Journal of Solids and Structures 14 (2017) 1416-1442 The comparison between the contact force predicted by the present simulations to that experimentally detected of Figure 8 shows that a better predictive capability is achieved with respect to previous simulations by Icardi and Zardo (2005), thanks to the improved structural and damage models here used.Note that while in Icardi and Zardo (2005) the residual properties of the failed regions are guessed within the framework of the ply-discount theory, here they are computed as outlined in Section 3.2.
The centre of each arbitrarily chosen sub-region for damage analysis corresponds to a point pi where the damage criteria of Section 3 are applied in order to determine whether or not materials are failed.If in pi one of the damage criteria predicts failure, the corresponding sub-region is considered damaged and therefore it is coloured in grey in Figure 9. Stresses in pi are determined assuming the modified expression of the strain energy resulting by the damage mesoscale model (3.2) with the damage indicators computed by PFEA-3.
The damage predicted with and without enforcement of the in-plane stress continuity is compared to that detected via ultrasonic cartography in Figure 9.This figure shows the considerable practical importance of contributions (28) in the VK-ZZ model, which thus don't merely represent a theoretical achievement.An elapsed time of 490 s being required to complete the analysis, consideration of in-plane continuity does not adversely affect computational costs.
Figure 9 shows that the overlapped damaged area measured through C-SCAN (dashed line) by Icardi and Zardo (2005), is much better represented considering in-plane stress continuity (lightgrey regions) than disregarding it (dark-grey regions).Thus, in-plane stress continuity is shown to have a significant bearing in improving the accuracy of results However, a correct prediction of the overlap damage is obtained in any case.
Latin American Journal of Solids and Structures 14 (2017) 1416-1442 Table 1 shows the extension of the delaminated area at each interface, as measured in Icardi and Zardo (2005) and as computed by the present simulation procedure and by previous ones.The results show that the present refined structural model achieves at each interface a better accuracy.It is reminded that an analytical model is used in this paper, instead of the finite element of Icardi and Zardo (2005).However, square sub-regions similar to finite elements are considered for the damage analysis.The results of Figures 8 and 9 show that the present closed-form approach does not result in an accuracy loss with respect to finite element modelling, though its computational cost is much lower.So use of finite elements can be limited to cases with a complex shape.Icardi and Zardo (2005) and by simpler models (FSDPT and Icardi and Zardo (2005)).

CONCLUDING REMARKS
A refined, 3D zig-zag structural model has been developed for simulating low-velocity impacts on laminated and sandwich composite plates.It a priori fulfils the continuity of stresses at the layer Latin American Journal of Solids and Structures 14 (2017) 1416-1442 The comparison of simulations with and without enforcement of in-plane stress continuity shows the non-negligible effects of such continuity on accuracy of results.
The capability of the present impact simulation is corroborated by the comparison with the damage detected through ultrasonic cartography.

Figure 1 :
Figure 1: Stages of the procedure to solve the problem.

l
make continuous xy s moving in x and y;

Figure 2 :
Figure 2: Illustration of the interface where in-plane continuity is imposed.

Figure 3 :
Figure 3: In-plane variation of the normalized in-plane shear stress for the two material wedge by the present model (inset: exact solution by Hein and Erdogan (1971)).
constitute a contemporaneous assessment of the contact force simulation, of the VK-ZZ model used to represent Latin American Journal of Solids and Structures 14 (2017) 1416-1442the structural response and of the damage model, as they are simultaneously employed.So there is no possibility to separately ascertain these three contributions by the comparison to experiments.

Figure 6 :
Figure 6: Contact force for a sandwich plate impacted at 8 J by experiment by Kärger et al. (2007),the present model and the previous model byIcardi and Ferrero (2009).

Figure 7 :
Figure 7: Comparisons between experimental results by Schubel et al. (2005) and numerical contact forces for a sandwich plate with foam core considering an impact energy of a) 7,75 J, b) 41,1 J and c) 108 J.

Figure 8 :
Figure 8: Comparison between the experimental contact force in Icardi and Zardo (2005) and that computed with the present procedure for a laminate with I stiffeners.

Figure 9 :
Figure 9: Overlap induced damage for a laminate with I stiffeners by: C-SCAN in Icardi and Zardo (2005) (dashed line), present procedure without considering in-plane continuity (dark-grey square sub-regions) or considering it (light-grey square sub-regions).

Table 1 :
Delaminated area [mm2]at the interfaces for a laminate with I stiffeners by the present procedure and by results of experiment in