nova página do texto(beta)
Inglês (pdf)
Artigo em XML
Tradução automática
Citado por SciELO 
Citado por Google 
Similares em
SciELO 
Similares em Google 
Permalink
1 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-of-plane 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 2001, 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 (Oñate et al. 2012a, 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 and Wanji 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. 2009, 2010) and (Iurlaro et al. 2014). Accurate results have been obtained by (Icardi and Sola 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 middle-plane. 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 C0 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 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 obtained considering a damage mesomechanic model that accounts for the effects of hard discontinuities in homogenized form.
2 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, Rcontact is the radius of the contact area, while σ(r), σ(0)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 Rcontact 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.
2.1 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 (RADIOSS -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. The literature shows that υ= p= 1.5 and q= 2.5 is very often the most appropriate choice for laminates.
2.2 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 Δdi in order to make the structure in equilibrium with Fc i. Displacements are applied during iterations till the impactor and the indentation radii conform:
[K(d)] T being the tangent stiffness matrix computed using the VK-ZZ model. The updated vector of displacement d.o.f. di is computed as:
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.
3 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.
3.1 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:
Xt being the tensile strength of fibres, S 12=13 the in-situ shear strength of the matrix resin and σ 11, τ 12, τ 13 the tensile and shear stresses acting on the fibres. Compressive failure is predicted by:
Xc being the compressive strength of fibres. Matrix failure under traction is described by:
while under compression it is described by:
Delamination failure is predicted using the criterion developed by (Choi and Chang 1992):
where
Yn
+1=
Yn
+1
t
if
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):
σcμ and τ lu 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 R 11/ σv :
where I
1
=σ
11
+ σ
22 + σ
33 and
Parameters α1 = k
2(d - 1) / d, d = R
-
11 / R
+
11, α2 = (k
2 / d) - 1 and
Core crushing failure under out-of-plane normal and shear stresses is also predicted using the criterion by (Lee and Tsotsis 2000):
Zc , Sx , Sy 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.
3.2 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
Um
,
εm
and
σm
, respectively. These quantities are expressed as the sum of the solution of a problem
The meaning of symbols is as follows: Ī
22,
Ī
12
, Ī
23
, Ī
13 and Ī
33 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; [
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.
4 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.
4.1 Linear Contribution
It provides the basic linear contribution across the thickness and contains the five functional d.o.f.
u0 (x, y), v0 (x, y), w0 (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.
4.2 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 Ax 1… Azn 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).
4.3 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 C0 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.
4.4 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:
The continuity functions ujθkx make continuous σ xx in the x direction, while their counterparts ujθky do the same in y.
Functions
vjθky
are aimed at making continuous σ
yy
in the y direction, while
vjθkx
make it continuous in x. In a similar way,
uj⋋kx
and
vj⋋ky
make continuous σ
xy
moving in x and y;
wjηkx
and
wjηky
make continuous σ
xz
, while,
wj
Υ
kx
and
wj
Υ
ky
make continuous σ
yz
. Each contribution
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 and Sola 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.
4.5 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.
4.6 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 u0, v0, w0, γx 0 and γy 0 expressed as:
for clamped edges, or as:
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 eigen-frequencies, with just a reduction of the processing time less than 15%. Hence, the use of a simplified matrix is left to a future study
5 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.
5.1 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).
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)).
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).
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 of materials with dissimilar properties. Figure 3 reports the in-plane variation of the shear stress σxy across the material interface normalized as follows:
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 u0, v0, w0, γ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 uj⋋kx and vj⋋ky 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.
5.2 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 105 Hz is used.
The results for this case by the present simulation presented in Figures 4 and 5 constitute a contemporaneous assessment of the contact force simulation, of the VK-ZZ model used to represent the 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 4 Comparison between simulated and experimental contact force ((Sola 2016)) for a laminated [45/0/-45/90]2s plate impacted at 1.5 J.
Figure 5 Comparison between simulated and experimental contact force ((Sola 2016)) for a laminated [45/0/-45/90]2s plate impacted at 2.5 J.
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.
5.3 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.
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 by (Icardi and Ferrero 2009).
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/m3, 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.
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.
In the next section, the damage predicted by the present simulations is compared to that detected via ultrasonic cartography.
5.4 Importance of In-Plane Continuity
A composite panel with I stiffeners having a length of 800 mm (Lx), a width of 330 mm (Ly) and a skin thickness (t) of 3 mm formerly tested by (Icardi and Zardo 2005) is now considered. The panel is clamped at the short sides and is kept free at the long sides. It is impacted at the centre with impact energy of 40 J (mass of the impactor 5.45 kg, velocity 3.83 m/s). The skin layers are 0.25 mm thick, they are stacked with a [45°/-45°/0°/0°/45°/-45°/-45°/45°/0°/0°/-45°/45] sequence and have the following mechanical properties: E1= 130 GPa, E2= 8 GPa, G12= 5 GPa, G13= 5 GPa, G23= 2.5 GPa, ν12= 0.3, ρ= 1557 kg/m3. The strengths are: St11= 1.67 GPa, St22= 0.06 GPa, Sc11= 1.08 GPa, Sc22= 0.17 GPa, S12= S13= S23= 0.07 GPa, where Stii represents the strength to traction in directions i, Scii represents the strength to compression in directions i while Sij represents the shear strength. The impactor is made of steel (E= 210 GPa, ν= 0.3) and has a nose sphere with a diameter of 25.4 mm.
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.
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.
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.
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.
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).
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 (light-grey 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.
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.
Table 1 Delaminated area [mm2] at the interfaces for a laminate with I stiffeners by the present procedure and by results of experiment in (Icardi and Zardo 2005) and by simpler models (FSDPT and (Icardi and Zardo 2005)).
| Physical interface | Experimental by (Icardi and Zardo 2005) | Present | Present without continuity | Basic Model | Model by (Icardi and Zardo 2005) |
|---|---|---|---|---|---|
| 1st | 960 | 950 | 930 | 1423 | 950 |
| 2nd | 790 | 758 | 740 | 1096 | 258 |
| 3rd | 430 | 400 | 380 | 590 | 376 |
| 4th | 310 | 250 | 210 | 490 | 143 |
| 5th | 160 | 115 | 107 | 289 | 114 |
| 6th | 135 | 102 | 98 | 190 | 108 |
| 7th | 95 | 75 | 62 | 125 | 75 |
| 8th | 50 | 44 | 38 | 100 | 43 |
6 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 interfaces and at the interfaces of undamaged and damaged regions (where material properties suddenly vary), as prescribed by the elasticity theory for keeping equilibrium.
It gives an accurate representations of strain energy density, displacement, and strain and stress fields at each time step with a low computational effort from constitutive equations.
Five unknown variables describe the displacement fields, nevertheless the representation can be refined across the thickness as desired. The continuity of out-of-plane stresses at the layer interfaces and across damaged/undamaged interfaces is a priori fulfilled with an appropriate choice of zig-zag amplitudes. A modified version of the Hertzian contact law is used for solving the contact problem. It forces the target to conform to shape of the impactor when soft core sandwiches are analysed. The Newmark’s implicit time integration scheme is used for solving transient dynamic equations.
Stress-based criteria with a separate representation of failure modes are employed to identify the damaged regions at each time step under the contact force. The Hashin and Rotem’s criterion with in-situ strengths predicts fibre and matrix’s failures. Delaminated regions at the interfaces are predicted by the Choi-Chang’s criterion. The failure of honeycomb is predicted separately by the criteria by Besant et al. and Lee and Tsotsis, computing the failed region as the envelope of failures predicted by these two criteria. The failure of foam core is estimated according to the rule suggested by Evonik for Rohacell core.
The progressive damage analysis is carried out accounting for the effects of hard discontinuities (i.e. matrix cracks, fibre failures and delaminations) in homogenized form by a continuum damage mesomechanic model. A modified version of the strain energy is obtained, from which stresses are computed. The damage computed at the previous step is extended at the next step to the points where the ultimate stress is reached. Instead of guessing the residual properties, they are computed in a physically consistent way.
Small square cells of the domain are identified, where the damage of layers is assumed uniform and equal to that computed at the central point. The damage indicators of the mesoscale model are determined via 3-D FEA at various load levels (computations carried apart once for all). 3-D FEA is also used to predict the variable apparent elastic moduli of the core while it collapses/buckles under transverse compressive loading. In this way, stresses are accurately computed across the core, even if a homogenized structural model is used.
The impact problem is solved in analytic form within the framework of the Galerkin’s method assuming appropriate displacement trial functions for representing the in-plane variation of the functional d.o.f.
The application to a two-material wedge problem shows the present structural model able of capturing continuous stresses even under a sudden variation of properties like across undamaged and damaged regions, with a discontinuous gradient, as in the exact elasticity solution.
The simulated contact force time-history is shown to be in a close agreement with experiments for all laminated and sandwich composite sample plates considered and for different level of impact energy. This result shows the reliability of contact force simulations and of the progressive damage model adopted. The comparison with previous analyses reveals the improved predictive capability of the structural models and of the progressive failure simulation here adopted, while the computational time is not considerably increased.
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.
