## Abstract

This paper revisits a three-dimensional analytical approach to study internal instability in layered composites, when the behaviour of each component of the material is described by the three-dimensional equations of solid mechanics. It shows the development of a unified computational procedure for numerical realization of the three-dimensional analytical method as applied to various constitutive equations of the layers and fibres, and different loading schemes (uniaxial or biaxial loading). The paper also contains many examples of calculation of critical controlled parameters for particular composites as well as analysis of different buckling modes. The results of this method can be used as a benchmark for simplified models.

This article is part of the themed issue ‘Multiscale modelling of the structural integrity of composite materials’.

## 1. Introduction

The new composite materials, made up of advanced fibres and toughened epoxy resins, are lighter/stiffer/stronger and increase fuel efficiency in aeroplanes, compared with the aluminium currently used. In the twenty-first century carbon fibre reinforced plastics (CFRP) can and will contribute more than 50% of the structural mass of an aircraft [1]. They are extensively used in the Airbus A380 super jumbo, the first fully double-decked passenger jet; and the new Boeing 787 ‘Dreamliner’ aircraft—the primary structure (including the wing and fuselage) of the 787 are built mostly from composite materials.

The compressive strength is often a design-limiting factor for advanced layered materials: it is generally 30–40% lower than the tensile strength [2–5]. A better understanding of the compressive strength and failure mechanisms is therefore fundamental to the development of improved materials. It should be emphasized that zones of compressive stresses can appear in composite structures even under tensile loads. They could be due to the presence of holes, cut-outs and cracks, or generated by impact.

Previous experimental studies have revealed that a possible mechanism of failure initiation is fibre or layer microinstability (microbuckling) that may usually occur in regions where high stress gradients exist, for instance, near the edge of a hole or free edges. Dow & Grunfest [6] were the first to consider the microbuckling of fibres as a form of fracture of a unidirectional composite under compression. Since then, the beginning of the fracture process under compression is usually associated with the buckling of the microstructure of the material when the critical load is determined by parameters characterizing the microstructure of the composite rather than by the dimensions and shape of the specimen or structural member, i.e. with the internal instability phenomena according to Biot [7]. In this paper, we adopt the same assumption of linking the onset of fracture and the loss of stability in the internal structure of the material. The task of deriving three-dimensional analytical solutions to describe the response of layered materials was always considered as one of great importance [8]. Analytical solutions, if obtained, enable us to analyse the behaviour of a structure for the wide range of material properties, and loading schemes, without the restrictions imposed by simplified approximate methods. This paper revisits a three-dimensional analytical approach to study internal instability in layered composites, when the behaviour of each component of the material is described by the three-dimensional equations of solid mechanics. It shows the development of a unified computational procedure for numerical realization of the three-dimensional analytical method as applied to various constitutive equations of the layers and fibres, and different loading schemes (uniaxial or biaxial loading). The paper also contains many examples of calculation of critical controlled parameters for particular composites as well as analysis of different buckling modes. The results of this method can be used as a benchmark for simplified models.

## 2. On analytical models for microbuckling

Probably, the first solutions to the problem of internal instability for a layered material obtained within this approach were reported by Guz [9], where the problem for linear-elastic layers under uniaxial compression was solved. This solution was included in numerous books and the comprehensive review on the topic [10]. Later the exact solutions were derived also for more complex problems: for orthotropic, nonlinear elastic and elastic–plastic, compressible and incompressible layers including the case of large (finite) deformations (see, for example, [11–14] and the reviews [10,15,16]).

The importance and the complexity of the considered phenomena led to a large number of publications which put forward various approximate methods aimed at tackling the problems with different levels of accuracy: the early papers by Rosen [17], Schuerch [18], Sadovsky *et al.* [19] and the numerous later publications reviewed by Budiansky & Fleck [2], Soutis & Turkmen [3], Schultheisz & Waas [4], Niu & Talreja [5]. It was concluded after the detailed analyses [3,5,10,15,16,20,21], that the approximate methods are not very accurate when compared with experimental measurements and observations. For instance, one of the earlier models suggested by Rosen [17] involves considerable simplifications, modelling the reinforcement layers by the thin beam theory and the matrix as an elastic material using one-dimensional stress analysis. It makes the results of this method inaccurate even for simple cases. It was shown [10,16,22] that such approximate models can give a significant discrepancy in comparison with the three-dimensional approach used in this paper and with experimental data even for the simplest case of a composite with linear elastic compressible layers undergoing small precritical deformations and considered within the scope of geometrically linear theory. For small fibre volume fractions, the approximate model gives physically unrealistic critical strains and predicts a different mode of stability loss from that obtained by the accurate three-dimensional analysis. The three-dimensional approach used throughout this paper allows us to take into account large deformations, geometrical and physical nonlinearities and load biaxiality that the simplified methods cannot consider.

Another approach, which is commonly used, is based on the investigation of fibre kinking. From the early literature on compressive fracture, it was easy to get the impression that fibre instability (microbuckling) and kinking are competing mechanisms. However, it is now accepted that a kink band is an outcome of the microbuckling failure of actual fibres, as observed experimentally [23,24]. Fibre microbuckling occurs first, followed by propagation of this local damage to form a kink band. A comprehensive comparative analysis of the Rosen model, Argon–Budiansky (kinking) model, and Batdorf–Ko model was presented by Soutis & Turkmen [3]. Studies of the kinking phenomenon were also reviewed by Budiansky & Fleck [2]. Soutis & Turkmen [3] showed that the existing kinking analyses are able to account for some, but not all, of the experimental observations. They correctly predict that shear strength and fibre imperfections are important parameters affecting the compressive strength of the composite. However, within this model it is not possible to say exactly how the strength will vary with fibre content; and the value of misalignment is chosen arbitrarily. This model requires knowledge of the shear strength properties, the initial fibre misalignment and, most importantly, the kink-band orientation angle, which is a post-failure geometric parameter. The analysis of this approach is outside the scope of the present paper.

## 3. Developing solutions for three-dimensional microbuckling problems

Let us briefly consider the statement of the problem of internal instability (microbuckling) for layered materials. The detailed formulations for particular types of layers were given previously, for example, in [15,20,21].

The material consists of alternating layers with thicknesses 2*h*_{r} and 2*h*_{m} (figure 1). Two different loading schemes are studied: the uniaxial compression and the biaxial compression in the plane of the layers. The solution of the problem is sought for four modes of stability loss (e.g. [16,20]). Using the equations of the three-dimensional stability theory [25], the following eigen-value problem must be solved:

The stability equations for each layer are
3.1where *t*_{ij} is the non-symmetrical Piola–Kirchhoff stress tensor (nominal stress tensor), which has the following form for compressible solids:
3.2and for incompressible solids
3.3The components of the tensors *κ*_{ijαβ} and *ω*_{ijαβ} depend on the properties of the layers and the loads. The most general expressions for *κ*_{ijαβ} and *ω*_{ijαβ} can be found in [25]:
3.4and
3.5where and are the quantities which characterize the axial and shear stiffnesses. The quantity characterizing the precritical state (the stress component or the strain component ) is the controlled parameter in respect to which the eigen-value problem should be solved.

To complete the problem statement, the boundary conditions should be defined for each interface. The layer interfaces could consist of zones of perfectly connected (bonded) layers and defects such as cracks or delaminations. In this study, we consider composites with perfectly bonded layers or ‘perfectly lubricated’ (sliding without friction) interfaces. For the perfectly bonded layers, we have the continuity conditions for the stresses and displacements 3.6and for ‘the perfectly lubricated layers’ [26,27], only the continuity of the normal components is retained at the surface, with boundary conditions for perturbations of stresses and displacements in the form of 3.7

Solutions of equation (3.1) (i.e. perturbations of stresses and displacements) for each of the layers can be expressed through the functions *X* and *Ψ*, which are the solutions of the following equations [25]:
3.8The parameter *ξ*_{j} depends on the components of the tensor *κ*_{ijαβ} (or *ω*_{ijαβ}) and, therefore, on the properties of the layers and on the loads. It was proved in [11,25] that for elastic compressible and elastic incompressible layers
3.9and for elastic–plastic incompressible layers
3.10The characteristic determinants associated with the four modes of stability loss were derived earlier, for example, in [11,12,15,16] for various constitutive equations of the layers, different loading schemes (uniaxial or biaxial loading) and different precritical conditions (large or small precritical deformations). The four modes of stability loss include all possible periodic modes with periods which are equal to one or two periods of the internal structure. Similarly, characteristic equations can be derived for other modes of stability loss. The described method can give the solutions for modes with periods which are equal to 3, 4, 5, … periods of the internal structure. Other modes with periods which are not multiples of the period of the internal structure can also be examined. The solution for them would be based either on the Floquet theorem for ordinary differential equations with periodic coefficients using the formalism of Brillouin [28], or on reducing the problem to an infinite set of equations with the consequent solution by a numerical method [10]. However, the modes with the larger periods in transverse direction are usually not of practical interest, as shown in [10,20]. In this paper, we attempt to present the characteristic determinants in a unified form in order to facilitate a uniform computational procedure for solving them:

— for perfectly bonded layers 3.11

— for perfectly lubricated layers 3.12

The expressions for *β*_{ij} of the determinant for different types of materials and for different loading schemes are given in [1,10,13,14,20–22] and other publications mentioned there.

## 4. Results and discussion

To facilitate the analysis of characteristic determinants (3.11) and (3.12), a software package with a graphical user-friendly interface was developed using Matlab v. 7.6.0. The software contains the database of material properties for typical layered composites and the library of components of the tensors *κ*_{ijαβ} and *ω*_{ijαβ} given by equations (3.4) and (3.5). The fully automated numerical procedure consists of the following steps. First, the characteristic determinants, equations (3.11) or (3.12), are computed depending on the user’s choice of loading schemes (uniaxial or biaxial loading), initial conditions (large or small precritical deformations) and interfacial properties (perfectly bonded and perfectly lubricated layers). Then the results are analysed, and the critical controlled parameters of the internal instability (including the critical wavelength) are searched for. This analysis is conducted for all four considered modes of stability loss. At the final stage, the modes are compared and the critical mode is found. Some of the results for the cases of perfectly bonded and perfectly lubricated layers are presented below for two different types of layers—compressible linear elastic and incompressible elastic–plastic.

Let us consider the following layered composite (figure 1): the reinforcement behaves as a linear-elastic isotropic compressible material, and the matrix response is elastic–plastic incompressible described by the following relationship for equivalent stress and strain :
4.1where *k*_{m} and *A*_{m} are material constants for the elastic–plastic matrix. The constitutive equation (4.1) is typical for metal matrix composites [12,29,30]. For the reinforcement layer, we have
4.2where *E* is Young’s modulus and *ν* is Poisson’s ratio.

The components of tensors *ω*_{ijαβ} for such materials are given in [25,31] and the components of tensor *κ*_{ijαβ} in [20,25] for different types of loading. Following the procedure described in the previous section, i.e. substituting the expressions for *κ*_{ijαβ} and *ω*_{ijαβ} into the characteristic equations (3.11) and (3.12), the characteristic equation can be specified for the considered type of composite material, see [12,15] for more details. For all modes, we will have transcendental equations in terms of two variables: (applied strain) and *α*_{r} (normalized half-wavelength). Solving the characteristic equations for different modes of stability loss, the dependences are obtained (*N*=1,2,3,4 is the number of the mode). A minimum of the corresponding dependence is the critical value for the particular mode—*ε*^{(N)}_{cr}. The critical strain of internal instability for the considered layered material is the minimal of these four values (*ε*^{pl}_{cr} in the case of perfectly lubricated layers, and *ε*^{pb}_{cr} in the case of perfectly bonded layers):
4.3The computed values of critical strain for metal matrix elastic–plastic composites under biaxial and uniaxial loadings are presented in figures 2–5. The results show how the bonds between the layers and the material properties affect the solution for the first two modes of stability loss. Figures 2 and 3 correspond to the case of perfectly bonded layers and figures 4 and 5 to the case of perfectly lubricated layers.

In practical cases, the assumption of perfect bonding between neighbouring layers along the entire interface does not always correspond to reality due to different imperfections of interfacial adhesion. When considering a composite with such defects it is sometimes difficult to identify a set of the defects and its influence on the internal instability. Hence, we suggest the following bounds for the controlled parameter. The critical strain, *ε*_{cr}, for a composite with weakened interfacial adhesion must be larger than the critical strain *ε*^{pl}_{cr} for the same structure with perfectly lubricated layers, but smaller than the critical strain *ε*^{pb}_{cr} for the structure with perfectly bonded layers. Thus, we obtain the following bounds for the critical strain:
4.4These bounds are particularly suitable for the so-called ‘defects with connected edges’ (figure 6) also considered in [20,22]. These model defects are an idealization for the case when a change in the nature of the interlaminar contact occurs, when an interaction of the layers is implemented so that infinitesimal sliding is allowed, but still there are no gaps between the layers (e.g. molecular chains in some kinds of glue connection, figure 6). For defects with connected edges, the continuity at the interface is retained for normal components only, with boundary conditions for the perturbations of stresses and displacements in the form of equation (3.7). Composites can contain not only interlaminar, but also various sorts of intralaminar defects. The effect of intralaminar damage can be accounted for by considering layers with reduced stiffness properties following, for example, Kashtalyan & Soutis [32] or employing the methods developed by Kashtalyan & Rushchitsky [33] or Winiarski & Guz [34].

The results of computation for layered composites with elastic–plastic matrix are shown in figures 7 and 8. For the first and second modes of stability loss, the critical strain remains constant and the difference between the results for perfectly bonded and perfectly lubricated layers does not change, while the ratio of the layer thicknesses is lower than a certain value (0.027 for the case of figure 7). Then with the increase of the difference between layer thicknesses, the bounds for critical strain become narrower. The bounds for critical strain are shown in figure 8 as a function of *k*_{m}. With the increase of the coefficient *k*_{m}, the distance between the upper and the lower curves significantly decreases for the first mode of stability loss and remains almost the same for the second mode.

The computed bounds are presented for model generic material systems. They appear to give a reasonable estimation for the critical controlled parameters and may be considered as the first approximation on the way to the exact solution of the problem of stability in compression along interfacial defects. According to Guz & Soutis [13,14], Menshykova *et al.* [15] and Guz *et al.* [16], the method works better for layered materials with small to medium volume fraction of stiffer (reinforcing) layers. The applicability of the method to practical materials, e.g. composite materials used in aerospace, automotive and other industries, or layered rocks, should be discussed separately for each class of such materials. Further studies are required to compare the results with experimental observations and measurements.

## 5. Concluding remarks

The purpose of this work was to outline a procedure which allows us to develop three-dimensional analytical solutions and investigate the internal instability of different types of layered composites consisting of compressible elastic or/and incompressible elastic–plastic constituents. The analysis of different loading schemes and precritical conditions was carried out using the developed procedures. We have intentionally chosen to present the computations for model generic material systems only. Some comparisons with available experimental data were discussed earlier [10,22–24]. The applicability of the method to practical materials, e.g. composite materials used in aerospace, automotive and other industries, or layered rocks, should be discussed separately for each class of such materials. It would depend on many factors, such as the ability of the equations of Newtonian solid mechanics to fully capture the influence of fine microstructure, various types of defects usually present in real-life materials, the importance of considering more complex loading schemes, etc. In order to take such factors into account, some simplifying assumptions may be required when developing a robust solution. Then the presented analytical solution obtained within the three-dimensional theory of stability (albeit for a very particular model configuration with a particular loading scheme) can be used as a benchmark for those simplified methods.

Zhuk *et al.* [35] gave an example of one possible application of the model presented in this paper. Carbon fibre composite materials are sensitive to open holes, defects and low-velocity impact that can cause barely visible damage (BVID) that can significantly reduce their stiffness and strength properties. To develop structures which are more damage resistant and tolerant, it is necessary to understand how the damage is caused and how it can affect residual performance. A typical aircraft structure such as a fuselage shell or a wing surface usually consists of a skin reinforced with stiffeners. An analytical formula, based on three-dimensional stability theory, was presented by Zhuk *et al.* [35] for calculating the unnotched compressive strength of a multidirectional composite plate. Then the maximum stress failure criterion was employed to estimate the critical load of a stiffened panel with an equivalent open hole loaded in compression. In the range of the model applicability, critical loads predicted by the model were very close to the measured data from [36].

## Authors' contributions

I.A.G. and C.S. conceived of the study, designed the study, coordinated the study and drafted the manuscript; I.A.G. and M.M. designed and carried out the numerical analysis; M.M. conceived and designed the software tools and helped drafting the manuscript. All authors gave final approval for publication.

## Competing interests

We have no competing interests.

## Funding

Financial support of part of this research by The Royal Society, The Royal Society of Edinburgh, The Royal Academy of Engineering and The Carnegie Trust for the Universities of Scotland is gratefully acknowledged.

## Footnotes

One contribution of 22 to a Theo Murphy meeting issue ‘Multiscale modelling of the structural integrity of composite materials’.

- Accepted February 22, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.