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

Abstract

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 delamination 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 Galerkin’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.

Keywords:
Impact-induced damage; 3D variable-kinematics zig-zag model; undamaged/damaged stress continuity; interlayer stress continuity; fixed d.o.f; mesomechanic damage model

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 C⁰ 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:

$\sigma (r)=\sigma (0)\sqrt{1-\left(r^2/R_{contact}^2\right)}\text{if (}\sigma (r)=0\text{if}r>R_{contact}\text{)}$ (1)

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:

$F=K_c{\alpha }^{\upsilon }$ (2)

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.

Figure 1
Stages of the procedure to solve the problem.

The exponent υ is set in agreement with experiments in the literature. Also in the unloading phase

$$F=F_m{\left(\frac{\alpha -{\alpha }_0}{{\alpha }_m-{\alpha }_0}\right)}^q$$ (3)

the exponent q is set according to experiments. In the previous equation, F_m is the maximum of the contact force, α_m is the relative indentation depth at each loading and α₀ is the permanent indentation. In case of bounces, the following reloading law is assumed:

$F=k^{\text{'}}{}_{}{(\alpha -{\alpha }_0)}^p$ (4)

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:

$${\delta }^i={\left[K\left(d^{(i-1)}\right)\right]}_S{d}^{(i-1)}-F_c^{(i-1)}\ne 0$$ (5)

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ⁱ in order to make the structure in equilibrium with ^{_Fc i}. Displacements are applied during iterations till the impactor and the indentation radii conform:

$${\left[K\left(d^{(i-1)}\right)\right]}_T\Delta d^i={\delta }^i$$ (6)

[K(d)]_T being the tangent stiffness matrix computed using the VK-ZZ model. The updated vector of displacement d.o.f. dⁱ is computed as:

$$d^i=d^{(i-1)}+\Delta d^i$$ (7)

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

$${\overline{d}}^i=\sqrt{{\displaystyle \sum _n{\left(d_j^n\right)}^2}}$$ (8)

to compute

$$D=\frac{\left[{\overline{d}}^{i+1}-{\overline{d}}^i\right]}{{\overline{d}}^i}$$ (9)

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:

${\left(\frac{{\sigma }_{11}}{X^t}\right)}^2+\frac{1}{S_{12=13}{}^2}\left({\tau }_{12}{}^2+{\tau }_{13}{}^2\right)=1,({\sigma }_{11}>0)$ (10)

_^Xt being the tensile strength of fibres, S ₁₂₌₁₃ the in-situ shear strength of the matrix resin and σ ₁₁, τ ₁₂, τ ₁₃ the tensile and shear stresses acting on the fibres. Compressive failure is predicted by:

${\sigma }_{11}=-X^c\text{,}({\sigma }_{11}<0)$ (11)

_^Xc being the compressive strength of fibres. Matrix failure under traction is described by:

${\left(\frac{{\sigma }_{22}+{\sigma }_{33}}{Y^t}\right)}^2+\frac{1}{S_{23}{}^2}({\tau }_{23}{}^2-{\sigma }_{22}{\sigma }_{33})+{\left(\frac{{\tau }_{12}}{S_{12=13}}\right)}^2+{\left(\frac{{\tau }_{13}}{S_{12=13}}\right)}^2=1({\sigma }_{22}+{\sigma }_{33}>0)$ (12)

while under compression it is described by:

$\begin{array}{l}\frac{1}{Y^c}\left[{\left(\frac{Y^c}{2S_{23}}\right)}^2-1\right]({\sigma }_{22}+{\sigma }_{33})+\frac{{({\sigma }_{22}+{\sigma }_{33})}^2}{4S_{23}{}^2}+\frac{({\tau }_{23}{}^2-{\sigma }_{22}{\sigma }_{33})}{S_{23}{}^2}+\frac{({\tau }_{12}{}^2+{\tau }_{13}{}^2)}{S_{12=13}{}^2}=1\\ ({\sigma }_{22}+{\sigma }_{33}<0)\end{array}$ (13)

Delamination failure is predicted using the criterion developed by (Choi and Chang 1992):

$e_d^2=D_a{}^{}{\left[\frac{{\overline{\sigma }}_{yz}^n}{S_i^n}+\frac{{\overline{\sigma }}_{xz}^{n+1}}{S_i^{n+1}}+\frac{{\overline{\sigma }}_{yy}^{n+1}}{Y_i^{n+1}}\right]}^2>1$ (14)

where ^{_Yn +1}= ^{_Yn +1} _t if ${\stackrel{-}{\sigma }}_{yy}$ ≥ 0, or ^{_Yn +1} = ^{_Yn +1} c if ${\stackrel{-}{\sigma }}_{yy}$< 0, _^Da is an empirical constant that is set after consideration of the material properties, ${\stackrel{-}{\sigma }}_{ij}$ is the average stress at the interface between the n ^th ply and the n+1 ^th ply:

${\overline{\sigma }}_{ij}^{n+1}=\frac{1}{h_{n+1}}{\displaystyle \underset{t_{n-1}}{\overset{t_n}{\int }}{\sigma }_{ij}}dt$ (15)

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):

$${\left(\frac{{\sigma }_{zz}}{{\sigma }_{cu}}\right)}^n+{\left(\frac{{\sigma }_{xz}}{{\sigma }_{lu}}\right)}^n+{\left(\frac{{\sigma }_{yz}}{{\sigma }_{lu}}\right)}^n=e_{core}\text{(}{e}_{core}>1\text{)}$$ (16)

_^σ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 ₁₁/^_σv:

$${\sigma }_v=\frac{\sqrt{(12a_2+12a1+12)I_2+[4a_2{}^2+(4a_1+4)a_2+a_1{}^2]I_1{}^2}+a_1{I}_1}{2a_2+2a_1+2}$$ (17)

where I ₁ =σ ₁₁ + σ ₂₂ + σ ₃₃ and $$I_2=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.\left[{\displaystyle \sum _{i=1}^3{\sigma }_{ii}^2}-{\sigma }_{11}{\sigma }_{22}-{\sigma }_{11}{\sigma }_{33}-{\sigma }_{22}{\sigma }_{33}+3{\displaystyle \sum _{ij=\mathrm{12,13,23}}{\sigma }_{ij}{}^2}\right]$$

Parameters α₁ = k ²(d - 1) / d, d = R ^- ₁₁ / R ⁺ ₁₁, α₂ = (k ² / d) - 1 and $k=\sqrt{3}{R}_{12}/R_{11}^{+}$ are determined from experiments, while R ⁺ ₁₁, R ^- ₁₁ and 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):

$\frac{{\sigma }_{zz}}{Z^c}=1_{}{,}_{}\frac{{\sigma }_{xz}}{S^x}=1_{}{,}_{}\frac{{\sigma }_{yz}}{S^y}=1$ (18)

_^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 $\tilde{P}$ in which damage is removed and the solution of a residual problem $\stackrel{-}{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):

$\begin{array}{l}{\frac{2E_p^S}{\left|S\right|}}^{}=[\pi \overline{\epsilon }\pi {]}^t[{\overline{M}}_1({\overline{I}}_{22},{\overline{I}}_{12})][\pi \overline{\epsilon }\pi ]+{\tilde{\sigma }}_{33}[{\overline{M}}_2({\overline{I}}_{22},{\overline{I}}_{12})]{\tilde{\sigma }}_{33}+\\ \text{+}{\tilde{\sigma }}_{33}[{\overline{M}}_3({\overline{I}}_{22},{\overline{I}}_{12})][\pi \overline{\epsilon }\pi ]-\frac{(1+{\overline{I}}_{23}){\tilde{\sigma }}_{23}^2}{{\tilde{G}}_{23}}-\frac{(1+{\overline{I}}_{13}){\tilde{\sigma }}_{13}^2}{{\tilde{G}}_{23}}-\frac{(1+{\overline{I}}_{33})<{\tilde{\sigma }}_{33}^2{>}_{+}^2{}_{}}{{\tilde{E}}_3}\end{array}$ (19)

The meaning of symbols is as follows: Ī _22, Ī ₁₂ , Ī ₂₃ , Ī ₁₃ and Ī ₃₃ 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; [${\stackrel{-}{M}}_1$], [${\stackrel{-}{M}}_2$], [${\stackrel{-}{M}}_3$] are 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):

${\frac{2E_p^j}{\left|{\gamma }_j\right|}}^{}=-\frac{(1+{\overline{I}}^1){\tilde{\sigma }}_{13}}{{\tilde{k}}_1}-\frac{(1+{\overline{I}}^2){\tilde{\sigma }}_{23}}{{\tilde{k}}_2}-\frac{(1+{\overline{I}}^3){\tilde{\sigma }}_{33}}{{\tilde{k}}_3}$ (20)

${\stackrel{-}{k}}_1$, ${\stackrel{-}{k}}_2$ and ${\stackrel{-}{k}}_3$ being the elastic stiffness coefficients of the interface, Ī ¹, Ī ² and Ī ³ the three damage indicators and γ_j the deformation. The coefficients of stresses in equations (19) and (20) represent the elastic properties of the homogenized model equivalent to the discretely damaged real model. The damage is quantified in terms of cracking rate ρ = L/H and delamination ratio τ = e/h, 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 ρ is [0; 0.7] and for τ 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.

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:

$\begin{array}{l}\begin{array}{l}u(x,y,z)=U^0\left(x,y,z\right)+U^i\left(x,y,z\right)+U^c\left(x,y,z\right)\hfill \end{array}+U^{c\_ip}\left(x,y,z\right)\\ \begin{array}{l}v(x,y,z)=V^0\left(x,y,z\right)+V^i\left(x,y,z\right)+V^c\left(x,y,z\right)\hfill \end{array}+V^{c\_ip}\left(x,y,z\right)\\ \begin{array}{l}w(x,y,z)=W^0\left(x,y,z\right)+W^i\left(x,y,z\right)+W^c\left(x,y,z\right)+W^{c\_ip}\left(x,y,z\right)\hfill \end{array}\end{array}$ (21)

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.

$$\begin{array}{l}{U}^0(x,y,z)=u^0(x,y)+z\left[{\gamma }_x^0(x,y)-w^0{}_{,x}(x,y)\right]\\ V^0(x,y,z)=v^0(x,y)+z\left[{\gamma }_y^0(x,y)-w^0{}_{,y}(x,y)\right]\\ W^0(x,y,z)=w^0(x,y)\end{array}$$ (22)

u⁰ (x, y), v⁰ (x, y), w⁰ (x, y), γ_x ⁰ (x, y) and γ_y ⁰ (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:

$$\begin{array}{l}{U}^i(x,y,z)=A_{x1}z+A_{x2}{z}^2+A_{x3}{z}^3+A_{x4}{z}^4+\mathrm{...}+A_{xn}{z}^n\\ V^i(x,y,z)=A_{y1}z+A_{y2}{z}^2+A_{y3}{z}^3+A_{y4}{z}^4+\mathrm{...}+A_{yn}{z}^n\\ W^i(x,y,z)=A_{z1}z+A_{z2}{z}^2+A_{z3}{z}^3+A_{z4}{z}^4+\mathrm{...}+A_{zn}{z}^n\end{array}$$ (23)

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:

$$\begin{array}{l}\begin{array}{cccc}{\sigma }_{xz}\left(z=\pm \raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)=0& {\sigma }_{yz}\left(z=\pm \raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)=0& & \end{array}\\ \begin{array}{cccc}{\sigma }_{zz}\left(z=\raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)=p^0\left(z=\raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)& {\sigma }_{zz}\left(z=-\raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)=p^0\left(z=-\raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)& {\sigma }_{zz,z}\left(z=\pm \raisebox{1ex}{$h$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)=0& \end{array}\end{array}$$ (24)

and the fulfilment of local differential equilibrium equations at discrete points chosen across the thickness:

$$\begin{array}{l}{\sigma }_{xx,x}+{\sigma }_{xy,y}+{\sigma }_{xz,z}=b_x\\ {\sigma }_{xy,x}+{\sigma }_{yy,y}+{\sigma }_{yz,z}=b_y\\ {\sigma }_{xz,x}+{\sigma }_{yz,y}+{\sigma }_{zz,z}=b_z\end{array}$$ (25)

So the expansion order can be increased adding more points. Please notice that the computations of A_x1 … A_zn don’t imply adding new d.o.f. Symbols b_i 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:

$\begin{array}{l}{U}^c(x,y)={\displaystyle \sum _{k=1}^{n_i}{\Phi }_x^k(x,y)(z-z_k)H_k+{\displaystyle \sum _{k=1}^{n_i}{C}_u^k(x,y)H_k}}\\ V^c(x,y)={\displaystyle \sum _{k=1}^{n_i}{\Phi }_y^k(x,y)(z-z_k)H_k+{\displaystyle \sum _{k=1}^{n_i}{C}_v^k(x,y)H_k}}\\ W^c(x,y)={\displaystyle \sum _{k=1}^{n_i}{\Psi }^k(x,y)(z-z_k)H_k}+{\displaystyle \sum _{k=1}^{n_i}{\Omega }^k(x,y){(z-z_k)}^2{H}_k+}{\displaystyle \sum _{k=1}^{n_i}{C}_w^k(x,y)H_k}\end{array}$ (26)

H_k being the Heaviside unit step function ((Di Sciuva 1985)), i.e. H_k=0 for z < z_k, while H_k=1 for z ≥z_k. 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 C_u ^k, C_v ^k and C_w ^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:

$$\begin{array}{l}\begin{array}{cccc}{\sigma }_{xz}{|}_{{}^{{}^{(k)}}{z}^{+}}={\sigma }_{xz}{|}_{{}^{{}^{(k)}}{z}^{-}}& {\sigma }_{yz}{|}_{{}^{{}^{(k)}}{z}^{+}}={\sigma }_{yz}{|}_{{}^{{}^{(k)}}{z}^{-}}& {\sigma }_z{|}_{{}^{{}^{(k)}}{z}^{+}}={\sigma }_z{|}_{{}^{{}^{(k)}}{z}^{-}}& {\sigma }_{z,z}{|}_{{}^{{}^{(k)}}{z}^{+}}={\sigma }_{z,z}{|}_{{}^{{}^{(k)}}{z}^{-}}\end{array}\\ \begin{array}{ccc}{u|}_{{}^{\left(k\right)}z{}^{+}}={u|}_{{}^{\left(k\right)}z{}^{-}}& {v|}_{{}^{\left(k\right)}z{}^{+}}={v|}_{{}^{\left(k\right)}z{}^{-}}& {w|}_{{}^{\left(k\right)}z{}^{+}}={w|}_{{}^{\left(k\right)}z{}^{-}}\end{array}\end{array}$$ (27)

Algebraic difficulties are overcome using symbolic calculus (MATLAB® symbolic software package). As a result, expressions for Φ_x ^k, Φ_y ^k, Ψ^k, Ω^k, C_u ^k, C_v ^k and C_w ^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⁰ 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:

$\begin{array}{l}{U}^{c\_ip}=\sum _{j=1}^s\sum _{k=1}^{n_l}{}_u{}^j{\theta }_x^k\left(y\right)\left(x-x_k\right)H_k+\sum _{j=1}^s\sum _{k=1}^{n_l}{}_u{}^j{\theta }_y^k\left(x\right)H_k+\sum _{j=1}^s\sum _{k=1}^{n_l}{}_u{}^j{\lambda }_x^k\left(x\right)\left(y-y_k\right)H_k\\ V^{c\_ip}=\sum _{j=1}^s\sum _{k=1}^{n_l}{}_v{}^j{\theta }_x^k\left(y\right)H_k+\sum _{j=1}^s\sum _{k=1}^{n_l}{}_v{}^j{\theta }_y^k\left(x\right)(y-y_k)H_k+\sum _{j=1}^s\sum _{k=1}^{n_l}{}_v{}^j{\lambda }_y^k\left(y\right)(x-x_k)H_k\\ W^{c\_ip}=\sum _{j=1}^s\sum _{k=1}^{n_l}{}_w{}^j{\eta }_x^k\left(y\right)\left(x-x_k\right)H_k+\sum _{j=1}^s\sum _{k=1}^{n_l}{}_w{}^j{\eta }_y^k\left(x\right)H_k+\sum _{j=1}^s\sum _{k=1}^{n_l}{}_w{}^j{\mathrm{{\rm Y}}}_x^k\left(y\right)H_k+\\ +\sum _{j=1}^s\sum _{k=1}^{n_l}{}_w{}^j{\mathrm{{\rm Y}}}_y^k\left(x\right)\left(\mathrm{y}-y_k\right)H_k\end{array}$ (28)

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 ${\displaystyle \sum _{j=1}^s}$ makes continuous the stresses in either x or y directions, while $${\displaystyle \sum _{k=1}^{n_l}}$$ 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 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 u⁰, v⁰, w⁰, γ_x ⁰ and γ_y ⁰ expressed as:

$\begin{array}{l}{u}^0\left(x,y,t\right)=\sum _{m=\mathrm{1,2},3}^M\sum _{n=\mathrm{1,2},3}^N{A}_{mn}\left(t\right)\sin \left(\frac{m\pi }{L_x}x\right)\sin \left(\frac{n\pi }{L_y}y\right);v^0\left(x,y,t\right)=\sum _{m=\mathrm{1,2},3}^M\sum _{n=\mathrm{1,2},3}^N{B}_{mn}\left(t\right)\sin \left(\frac{m\pi }{L_x}x\right)\sin \left(\frac{n\pi }{L_y}y\right);\\ w^0\left(x,y,t\right)=\sum _{m=2,4}^M\sum _{n=2,4}^N{C}_{mn}\left(t\right)\left[\cos \left(\frac{m\pi }{L_x}x\right)-(-1{)}^{\raisebox{1ex}{$m$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right]\left[\cos \left(\frac{n\pi }{L_y}y\right)-(-1{)}^{\raisebox{1ex}{$n$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right];\\ {\gamma }_x^0\left(x,y,t\right)=\sum _{m=\mathrm{1,2},3}^M\sum _{n=\mathrm{1,2},3}^N{D}_{mn}\left(t\right)\sin \left(\frac{m\pi }{L_x}x\right)\sin \left(\frac{n\pi }{L_y}y\right);{\gamma }_y^0\left(x,y,t\right)=\sum _{m=\mathrm{1,2},3}^M\sum _{n=\mathrm{1,2},3}^N{E}_{mn}\left(t\right)\sin \left(\frac{m\pi }{L_x}x\right)\sin \left(\frac{n\pi }{L_y}y\right);\end{array}$ (29)

for clamped edges, or as:

$\begin{array}{l}{u}^0\left(x,y,t\right)=\sum _{m=\mathrm{1,2},3}^M\sum _{n=\mathrm{1,2},3}^N{A}_{mn}\left(t\right)\sin \left(\frac{m\pi }{L_x}x\right)\mathrm{co}\mathrm{s}\left(\frac{n\pi }{L_y}y\right);v^0\left(x,y,t\right)=\sum _{m=\mathrm{1,2},3}^M\sum _{n=\mathrm{1,2},3}^N{B}_{mn}\left(t\right)\mathrm{co}\mathrm{s}\left(\frac{m\pi }{L_x}x\right)\sin \left(\frac{n\pi }{L_y}y\right);\\ w^0\left(x,y,t\right)=\sum _{m=2,4}^M\sum _{n=2,4}^N{C}_{mn}\left(t\right)\left[\cos \left(\frac{m\pi }{L_x}x\right)-(-1{)}^{\raisebox{1ex}{$m$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right]\left[\cos \left(\frac{n\pi }{L_y}y\right)-(-1{)}^{\raisebox{1ex}{$n$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right];\\ {\gamma }_x^0\left(x,y,t\right)=\sum _{m=\mathrm{1,2},3}^M\sum _{n=\mathrm{1,2},3}^N{D}_{mn}\left(t\right)\sin \left(\frac{m\pi }{L_x}x\right)\mathrm{co}\mathrm{s}\left(\frac{n\pi }{L_y}y\right);{\gamma }_y^0\left(x,y,t\right)=\sum _{m=\mathrm{1,2},3}^M\sum _{n=\mathrm{1,2},3}^N{E}_{mn}\left(t\right)\mathrm{co}\mathrm{s}\left(\frac{m\pi }{L_x}x\right)\sin \left(\frac{n\pi }{L_y}y\right);\end{array}$ (29’)

for edges that are freely movable in the in-plane directions.

At clamped edges, the mechanical boundary conditions

$Q_x={\displaystyle {\int }_{\Gamma }{\sigma }_{xz}}dzdy\text{;}{Q}_y={\displaystyle {\int }_{\Gamma }{\sigma }_{yz}}dzdx$ (30)

are satisfied by an appropriate choice of terms (23), Γ being the contour of the plate. Q_x, Q_y 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 A_mn, …, E_mn 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 2
Illustration of the interface where in-plane continuity is imposed.

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:

$$\overline{{\sigma }_{xy}}=\frac{{\sigma }_{xy}}{{\sigma }_{xy}{}_{MAX}}$$ (31)

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⁰, v⁰, w⁰, γ_x ⁰ and γ_y ⁰ 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 10⁵ 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°₂/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/m³, E₁= 77 GPa, E₂= 75 GPa, G₁₂= 6.5 GPa, G₁₃= 5.1 GPa, G₂₃= 4.1 GPa, υ₁₂ = 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 (L_x), a width of 330 mm (L_y) 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: E₁= 130 GPa, E₂= 8 GPa, G₁₂= 5 GPa, G₁₃= 5 GPa, G₂₃= 2.5 GPa, ν₁₂= 0.3, ρ= 1557 kg/m³. The strengths are: S_t11= 1.67 GPa, S_t22= 0.06 GPa, S_c11= 1.08 GPa, S_c22= 0.17 GPa, S₁₂= S₁₃= S₂₃= 0.07 GPa, where S_tii represents the strength to traction in directions i, S_cii represents the strength to compression in directions i while S_ij 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 p_i where the damage criteria of Section 3 are applied in order to determine whether or not materials are failed. If in p_i 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 p_i 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.

References

• Abrate, S. (2001). Modeling of impacts on composite structures. Composite Structures; 51:129-38.
• Abrate, S. (2014). Impact on composite plates in contact with water. Procedia Engineering; 88: 2 - 9.
• Aktas, A., Aktas, M., Turan, F. (2013). The effect of stacking sequence on the impact and post-impact behaviour of woven/knit glass/epoxy hybrid composites. Composite Structures; 103: 119-135
• Berthelot, J.M. (2003). Transverse Cracking and Delamination in Cross-Ply Glass-Fiber and Carbon-Fiber Reinforced Plastic Laminates: Static and Fatigue Loading. Appl. Mech. Reviews; 56: 111-147.
• Besant, T., Davies, G.A.O., Hitchings, D. (2001). Finite element modelling of low velocity impact of composite sandwich panels, Composites Part A; 32: 1189-1196.
• Chai, G.B., Zhu, S. (2011). A review of low-velocity impact on sandwich structures. Journal of Materials: Design & Applications;225(4):207-30.
• Chakrabarti, A., Chalak, H.D., Iqbal, M.A., Sheikh, A.H. (2011). A new FE model based on higher order zig-zag theory for the analysis of laminated sandwich beam soft core. Composite Structures; 93: 271-279.
• Chakraborty, D. (2007). Delamination of Laminated Fiber Reinforced Plastic Composites Under Multiple Cylindrical Impact. Materials & Design; 28: 1142-1153.
• Chao, C.C., Tu, C.Y. (1999). Three-dimensional contact dynamics of laminated plates: Part 1. Normal impact. Composites Part B; 30: 9-22.
• Chen, W.J., Wu, Z. (2008). A selective review on recent development of displacement-based laminated plate theories. Recent Pat. Mech. Eng.; 1: 29-44.
• Choi, H.Y., Chang, F.K. (1992). A model for predicting damage in graphite/epoxy laminated composites resulting from low velocity point impact. J. Comput. Math.; 26: 2134-69.
• Damanpack, A.R., Shakeri, M., Aghdam, M.M. (2013). A new finite element model for low-velocity impact analysis of sandwich beams subjected to multiple projectiles. Composite Structures; 104: 21-33
• Di Sciuva. M. (1985). Development of an Anisotropic, Multilayered, Shear-Deformable Rectangular Plate Element.” Composite Structures; 21: 789-796.
• Diaz Diaz, A., Caron, J.J., Ehrlacher, A. (2007). Analytical Determination of the Modes I, II and III Energy Release Rates in a Delaminated Laminate and Validation of a Delamination Criterion. Composite Structures; 78: 424-432.
• Elder, D.J., Thomson, R.S., Nguyen, M.Q., Scott, L. (2004). Review of delamination predictive methods for low speed impact of composite laminates. Composite Structures; 66: 677-683.
• Garnich, M.R., Akula Venkata, M.K. (2009). Review of degradation models for progressive failure analysis of fiber reinforced polymer composites. Appl. Mech. Reviews; 62: 1-35.
• Hashin, Z., Rotem, A. (1973). A Fatigue Criterion for Fiber-reinforced Materials. J. of Compos. Mater.; 7: 448-464.
• Hein, V.L., Erdogan, F. (1971). Stress singularities in a two-material wedge. Int J of Fract Mech; 7: 317-330.
• Hongkarnjanakul, N., Bouvet, C., Rivallant, S. (2013). Validation of low velocity impact modelling on different stacking sequences of CFRP laminates and influence of fibre failure. Composite Structures; 106: 549-559.
• Horrigan, D.P.W., Staal, R.A. (2011). Predicting failure loads of impact damaged honeycomb sandwich panels- A refined model. J of Sandwich Structures and Materials; 13: 111-133.
• I.H. Choi. (2006). Contact force history analysis of composite sandwich plates subjected to low-velocity impact. Composite Structures 75 582-586
• Icardi, U., Ferrero, L. (2009). Impact Analysis of Sandwich Composites Based on a Refined Plate Element with Strain Energy Updating. Composite Structures; 89: 35-51.
• Icardi, U., Sola, F. (2013). Recovering critical stresses in sandwiches using through-the-thickness reinforcement. Composites Part B; 54: 269-277.
• Icardi, U., Sola, F. (2014a). Development of an efficient zig-zag model with variable representation of displacements across the thickness. J. of Eng. Mech.; 140: 531-541.
• Icardi, U., Sola, F. (2014b). Analysis of bonded joints with laminated adherends by a variable kinematics layerwise model. Int. J. of Adhesion & Adhesives; 50: 244-254.
• Icardi, U., Sola, F. (2014c). C0 Fixed Degrees of Freedom Zigzag Model with Variable In-Plane and Out-of-Plane Kinematics and Quadrilateral Plate Element. Journal Aerosp. Eng. In press, DOI: 10.1061/(ASCE)AS.1943-5525.0000472
  » https://doi.org/10.1061/(ASCE)AS.1943-5525.0000472
• Icardi, U., Sola, F. (2015). Indentation of sandwiches using a plate model with variable kinematics and fixed degrees of freedom. Thin-Walled Structures; 86: 24-34
• Icardi, U., Zardo, G. (2005). C0 plate element for delamination damage analysis, based on a zig-zag model and strain energy updating. International Journal of Impact Engineering; 31: 579-606
• Icardi, U., Sola, F. (2016). Assessment of recent zig-zag theories for laminated and sandwich structures. Composites Part B 97, 26-52
• Iurlaro, L., Gherlone, M., Di Sciuva, M. (2014). A mixed cubic zigzag model for multilayered composite and sandwich plates including transverse normal deformability. In Proc. 11th World Conference on Computational Methods, Barcelona, Spain 2014
• Kärger, L., Baaran, J., Teßmer, J. (2007). Rapid Simulation of Impacts on Composite Sandwich Panels Inducing Barely Visible Damage. Composite Structures; 79: 527-534.
• Kreja, I. (2001). A literature review on computational models for laminated composite and sandwich panels. Central European Journal of Engineering; 1: 59 - 80.
• Ladevèze, P., Lubineau, G., Marsal, D. (2006). Towards a bridge between the micro- and mesomechanics of delamination for laminated composites. Compos. Sci. & Tech.; 66: 698-712.
• Lee, S.M., Tsotsis, T.K. (2000). Indentation failure behaviour of honeycomb sandwich panels. Compos. Sci. & Tech.; 60: 1147-1159.
• Li, Q.M., Mines, R.A.W., Birch, R.S. (2000). The crush behaviour of Rohacell-51WF structural foam. Int. J. Solids Struct.; 37: 6321- 6341.
• Li, Y., Xuefeng, A., Xiaosu, Y. (2012). Comparison with Low-Velocity Impact and Quasi-static Indentation Testing of Foam Core Sandwich Composites. International Journal of Applied Physics and Mathematics; 2: 58-62
• Liu, P.F., Zheng, J.Y. (2009). Review on methodologies of progressive failure analysis of composite laminates. In: Koppel A, Oja J. Continuum mechanics, Nova Science Publishers, New York.
• Matsunaga, H. (2004). A comparison between 2-D single-layer and 3-D layerwise theories for computing interlaminar stresses of laminated composite and sandwich plates subjected to thermal loadings. Composite Structures; 64: 161-177.
• Mitrevski, T., Marshall, I.H., Thomson, R., Jones, R., Whittingham, B. (2005). The effect of impactor shape on the impact response of composite laminates. Composite Structures; 67: 139-148.
• Morinière, F.D., Alderliesten, R.C., Benedictus, R. (2014). Modelling of impact damage and dynamics in fibre-metal laminates-A review. International Journal of Impact Engineering; 67: 27-38.
• Murakami, H. (1986). Laminated Composite Plate Theory With Improved In- Plane Responses. ASME Appl. Mech.; 53:661-666.
• Oñate, E., Eijo, A., Oller, S. (2012a). A numerical model of delamination in composite laminated beams using the LRZ beam element based on the refined zigzag theory. Composite Structures; 104: 270-280
• Oñate, E., Eijo, A., Oller, S. (2012b). Delamination in laminated plates using the 4-noded quadrilateral QLRZ plate element based on the refined zigzag theory. Composite Structures; 108: 456-471
• Palazotto, A.N., Herup, E.J., Gummadi, L.N.B. (2000a). Finite Element Analysis of Low- Velocity Impact on Composite Sandwich Plates. Composite Structures; 49: 209-227.
• Palazotto, A.N., Herup, E.J., Gummadi, L.N.B. (2000b). Finite element analysis of low-velocity impact on composite sandwich plate. Composite Structures; 49: 209-227.
• Potluri, P., Kusak, E., Reddy, T.Y. (2003). Novel stitch-bonded sandwich composite structures. Composite Structures; 59: 251-259
• Richardson, M.O.W., Wisheart, M.J. (1999). Review of low-velocity impact properties of composite materials. Composites Part A; 27: 1123-1131.
• ROHACELL, W.F. (2011). Product Information. Evonik Röhm GmbH.
• Schubel, P.M., Luo, J.J., Daniel, I.M. (2005). Low velocity impact behavior of composite sandwich panels. Composites Part A; 36:1389-96.
• Sola, F. (2016), Simulation of local effects, energy absorption and failure mechanism in multilayeres and sandwich structures. PhD Dissertation Politecnico di Torino. Tutor: Ugo Icardi
• Tahani, M. (2007). Analysis of laminated composite beams using layerwise displacement theories. Composite Structures; 79: 535-547.
• Tessler, A., Di Sciuva, M., Gherlone, M. (2009). A Refined Zigzag Beam Theory for Composite and Sandwich Beams. J. Compos. Mat.; 43: 1051-1081.
• Tessler, A., Di Sciuva, M., Gherlone, M. (2010). A Consistent Refinement of First-Order Shear-Deformation Theory for Laminated Composite and Sandwich Plates Using Improved Zigzag Kinematics. J. Mech. Mat. Str.; 5: 341-367
• Wanji, C., Cheung, Y.K., Zhen, W. (2007). Augmented higher order global-local theory and refined triangular element for laminated composite plates. Composite Structures; 81: 341-352
• Yamada, Y., Okumura, H. (1983). Finite element analysis of stress and strain singularity eigenstate. In: Atluri, S.N., Gallagher, R.H., Zienkiewicz, O.C. Inhomogeneous Media or Composite Materials, Hybrid and Mixed Finite Element Methods, Ed. John Wiley & Sons Ltd, New York, 1983.
• Yamashita, M., Gotoh, M. (2005). Impact behavior of honeycomb structures with various cell specifications-numerical simulation and experiment. International Journal of Impact Engineering; 32: 618-630
• Yigit, A.S., Christoforou, A.P. (1995). On the impact between a rigid sphere and a thin composite laminate supported by a rigid substrate. Composite Structures; 30: 169-177
• Zhang, Y., Yang, C. (2009). Recent Developments in Finite Element Analysis for Laminated Composite Plates. Composite Structures; 88: 147-157.
• Zhen, W., Wanji, C., Xiaohui, R. (2009). Refined global-local higher-order theory for angle-ply laminated plates under thermo-mechanical loads and finite element model. Composite Structures; 88: 643-658
• Zhen, W., Wanji, C. (2006). An efficient higher-order theory and finite element for laminated plates subjected to thermal loading. Composite Structures; 73: 99-109
• Zhen, W., Wanji, C. (2010). A C0-type higher-order theory for bending analysis of laminated composite and sandwich plates. Composite Structures; 92: 653-661
• Zhou, D.W., Stronge, W.J. (2006). Low velocity impact denting of HSSA lightweight sandwich panel. International Journal of Mechanical Sciences; 48: 1031-1045.
• Zienkiewicz, O.C., Taylor, R.L. (2000). The finite element method, Butterworth-Heinemann, Oxford.

Publication Dates

• Publication in this collection
  16 Aug 2017

History

• Received
  28 Mar 2017
• Accepted
  23 May 2017