Thin-Walled Structures A novel semi-analytical approach for predicting the buckling and postbuckling behaviour of Glare fibre metal laminates

A semi-analytical model based on the variational principle, the Rayleigh–Ritz method and the Ramberg– Osgood formula is developed to analyse the buckling and postbuckling analysis of Glare ® fibre-metal laminates (FMLs) under axial compression. The model is applied to examine the failure behaviour of standard Glare ® 4B specimens. It allows computational efficient modelling of composite plates under pure or mixed stress boundary conditions with geometric imperfections, which have been shown to cause a considerable effect on their buckling and postbuckling behaviour. It is implemented in MATLAB. Results are compared with those obtained from a 3D Finite Element (FE) explicit dynamic nonlinear analysis implemented in the Abaqus/Explicit solver. The FE model incorporates progressive damage and failure using a cohesive zone model for inter-laminar layers and continuum material damage models for constituents, considering both geometric imperfections and load eccentricity. Finally, a series of specimens under compression loading are tested for the purpose of validation. Tests are monitored using Digital Image Correlation (DIC) for the determination and visualisation of principal strains, and Acoustic Emission (AE) for detection and location of damage initiation and evolution. Excellent correlation is observed between the analytical predictions, and both FEA and experimental results.


Introduction
In aerospace applications, thin-walled structures manufactured from lightweight materials such as Glare ® fibre-metal laminates are widely used. The application of these structures which primarily fail due to buckling can be extended by allowing them to operate in the postbuckling region, where they retain a certain amount of load carrying capacity. This clearly requires an accurate prediction of both buckling and postbuckling characteristics. Since the postbuckling behaviour of Glare plates involves both material and geometric nonlinearities due to plasticity, large deflections etc, there remains an ongoing need to develop analytical models or computationally inexpensive methods which incorporate these behaviours in a computationally efficient way.
When deriving analytical or semi-analytical approaches for laminates undergoing large in-plane deflections such as those found in postbuckling failure, the deflection of the middle surface of the plate due to out-of-plane displacement needs to be considered [1]. The essential nonlinear strain-shortening relations and partial differential equations for the large deflection of thin plates were first derived by von Kármán. These have been the basis of numerous analytical solutions employed to study the postbuckling behaviour of plates and are investigated in this work to ascertain their suitability for the postbuckling analysis of fibre-metal laminate plates. The closed form derived closed-form expressions for the evaluation of the initial postbuckling stiffness of composite plates using the principle of virtual work. Diaconu and Weaver [8,9] proposed an approximate closedform solution for the postbuckling analysis of infinitely long composite plates under in-plane compression. The postbuckling problem has also been solved by minimising the potential energy given in terms of the unknown deflection variables based on the nonlinear von Kármán strain-displacement relationship [10][11][12][13]. Other researchers have focused on the effect of anisotropic coupling of composite materials in postbuckling analysis as can be found in [5,7,14].
In previous works, the Rayleigh-Ritz (RR) energy method has been shown to provide an efficient methodology for the analysis of the behaviour of variable angle tow (VAT)laminated composite plates. Alhajahmad et al. [15,16] considered variable stiffness design tailoring for the nonlinear pressure-pillowing problem of fuselage skin panels based on the (RR) energy method. Wu et al. [17] proposed an energy method integrated with Airy's stress equations to study the prebuckling and buckling analysis of VAT composite plates. They adopted a semi-analytical approach using a single variational equation to solve the compatibility equations rather than solve them separately which required more effort. This single variational equation was also applied by Bisagni and Vescovini [18] to perform the postbuckling analysis of constant stiffness composite structures containing stiffeners. As well as the fact that the compatibility and equilibrium equations in addition to the boundary conditions (including both prescribed displacements and stresses) can be treated synchronously, this single variational approach avoids the need to calculate derivative terms of the stiffness in models of VAT laminates. The (RR) method is then applied to minimise the variational method resulting in a system of nonlinear algebraic equations from which the postbuckling equilibrium paths are followed using a modified Newton-Raphson procedure. Legendre polynomials are utilised in modelling the effects of flexural-twist anisotropy which has been shown to achieve an efficient and robust simulation for the postbuckling behaviour. Hanet al. [19] implemented a buckling analysis to investigate the response of VAT composite plates with holes and geometrical imperfections subject to compression. An angle variation formula to describe the reference fibre plies path and a method for generating a series of new paths were presented. Coburn et al. [20], developed an analytical model for the buckling analysis of a new blade stiffened VAT composite laminate to allow its characteristics to be fully investigated. The prebuckling and buckling analyses, performed on a representative section of the panel, are based on a generalised (RR) procedure. Wu et al. [21], extended their previous work on VAT laminates by developing a semi-analytical model to solve the postbuckling problem of VAT plates. Oliveri and Milazzo [22] presented a numerical model based on the (RR) approach for generally restrained multi-layered variable angle tow stiffened plates in the postbuckling regime. Their model, based on the first order shear deformation criterion introduces geometric nonlinearity through von Kármán's assumptions. Stiffened panel structures are modelled as an assembly of plate-like elements and penalty techniques are used to join the elements in the assembled structure and to apply kinematic boundary conditions. Lopatin and Morozov [23], presented a solution to the buckling of an orthotropic plate subject to a uniformly distributed compressive load applied to two fully clamped edges with the other two edges free (CCFF), similar to the boundary conditions used in this work. The problem was solved using the Kantorovich and generalised Galerkin methods. Talens [24] studied the buckling and postbuckling of composite plates with elastic restraints under different combinations of in-plane loading. The implemented solutions are based on thin plate theory for mid-plane symmetric plates. The governing equations are solved using a semi-analytical formulation combining the advantages of both analytical and numerical analyses. Galerkin and Ritz formulations are adopted for a more general solution. Zhang [25] proposed a new postbuckling analysis procedure to simulate both isotropic and anisotropic plates under combined loading based on the geometric relationship between the postbuckling axial stiffnesses calculated from the actual mode shapes and those calculated under the sinusoidal assumption. This enabled comparison studies to be conducted to correct the previous conservative postbuckling analysis.
Although extensive buckling and postbuckling studies using both analytical and finite element approaches have been conducted on composite laminates, relatively few have been carried out on fibre-metal laminates (FMLs). Bi et al. [26] developed an elasto-plastic model to simulate the buckling and postbuckling behaviour of FMLs under compressive loading. Considering both the geometric nonlinearity of the structure and the elasto-plastic deformation of the metal layers, an incremental von Kármán geometric relation with initial imperfection was established. An elasto-plastic constitutive relation adopting the mixed hardening rule was also introduced for the metal layers. Nonlinear governing equations for the FML were then derived, and the whole problem solved iteratively using the finite difference method. Bikakiset al. estimated the critical buckling stress of rectangular Glare laminates using a rule of mixtures with results validated using an FE model and eigenvalue buckling analysis [27]. Kamocka and Mania [28], studied the buckling and postbuckling behaviour of a FML channel section incorporating delamination using the ANSYS software. The authors [29,30] conducted both experimental and numerical analyses to investigate the behaviour of Glare laminates containing splice and doubler features. Load eccentricity and geometric imperfections were introduced in a three-dimensional nonlinear analysis performed in the Abaqus software using the Abaqus/Explicit solver. Banat and Mania [31], discussed nonlinear buckling and failure analysis of thin-walled FML structures under axial compression. Experimental procedures included macroscopic and microscopic examinations performed using stereoscopic, optical, and scanning electron microscopes (SEM). The results of these experimental investigations were compared with FE simulations wherein Tsai-Wu, Hashin and Puck failure criteria were applied.
Wittenberg et al. [32] proposed a methodology to calculate plasticity reduction factors for Glare flat plates loaded in either compression or shear, with a plasticity correction factor formula introduced by Gerard and Becker [33]. The bifurcation buckling load of the perfect plate was determined by following the elasto-plastic pre-buckling path, using Ramberg-Osgood (RO) relations to describe the stressstrain behaviour of the aluminium layers. The inelastic buckling mode was then enforced upon the plate geometry, serving as an initial imperfection. Finally, the nonlinear plate response was calculated for different amplitudes of imperfection. Kolakowskiet al. [34], studied the elastic-plastic buckling behaviour of thin-walled fibre-metal laminate structures subjected to axial compression. Three elastic-plastic criterion was adopted for the constitutive relations for the aluminium layers -fully elastic, the J2-plastic deformation criterion and the J2flow criterion. The latter two criterion have been used to derive the nonlinear constitutive relationships between stress and strain for a singular elastic-plastic component layer alongside the application of the (RO) formula. Recently, Prasad and Sahu [35] conducted both numerical and experimental studies on the buckling behaviour of Glare plates using first order Reissner-Mindlin theory in the finite element formulation. The effects of different parameters such as aspect ratio, width to thickness ratio, fibre orientation and various boundary conditions were examined. Mania et al. [36], evaluated the buckling and postbuckling behaviour of FML columns with different cross sections experimentally. Z-section specimens with different layups were tested under compression. Results were validated using finite element and analytical models based on Koiter's asymptotic theory.
In this paper, an efficient semi-analytical model is proposed to simulate the buckling and postbuckling behaviour of 'Glare ® 4B'fibremetal laminates. The method, based on a variational principle, the (RR) method and the (RO) criterion is implemented in a self-defined MAT-LAB programme. Geometrical imperfections are incorporated. Results are validated against a 3D dynamic explicit FE model implemented using Abaqus/Explicit. The FE model incorporates the full range of damage mechanisms and plasticity of the aluminium layers. To ensure the results are representative of the real structures both geometric imperfections and load eccentricity have been included. The former are based on eigenmode analyses with amplitudes based on measured imperfections, while the latter have been introduced based on an asymmetric linear load distribution applied in the form of a nonuniform distributed displacement. Results are validated experimentally with tests monitored using Digital Image Correlation (DIC) to visualise 3D full-field deformations and strains and Acoustic Emission (AE) to detect and localise the damage initiation and evolution.

Specimen geometry
The case examined in this work is a rectangular Glare FML plate with unsupported length 'a' = 100 mm and supported width 'b' = 80 mm.ie having clamped loading edges and free transverse edges (CCFF boundary conditions) (Fig. 1).
Specimens were sized and manufactured by the industrial partner to industry standards to be representative of the behaviour of typical aircraft panels in terms of their aspect ratio and width to thickness ratio and to ensure they were manufactured to the same standard as those used in practice.

Material properties
Specimens with standard lay-up type 'Glare 4B-3/2' consisting of 3 layers of 0.4 mm thick aluminium alloy 2024-T3 sheets and 2 'layers' of unidirectional (S2-glassfibre/FM94 epoxy resin) GFRP [37] were modelled in this work (with the particular grade of Glare from which the specimens were manufactured selected by the industrial sponsor). In this grade of Glare each GFRP 'layer' consists of a 3ply sub-laminate with the layup [90 o ∕0 o ∕90 o ] with each ply having a thickness of 0.133 mm after curing. This gave a total thickness of 2 mm ((3 × 2 mm) + (2 × 3 × 0.133 mm)).
Mechanical properties were based on the experimental results (Exp) provided by Wu and Yang [38][39][40], and presented in Table 1, which are very similar to those calculated by many researchers [38,39,[41][42][43] based on a Rule of Mixtures (ROM) formulation the values used Whilst these properties are given for Glare ® 4-3/2 laminates we use them here for Glare ® 4B-3/2 since the elastic properties are reported by references such as [44], to be 'almost the same' for all Glare4 laminate types prior to yielding in the aluminium layers. Since in our semi-analytical analysis we use only elastic properties to calculate critical buckling stress and a plasticity reduction factor including the introduction of geometric imperfections, we consider the use of the material properties given for Glare ® 4-3/2 laminates is appropriate. Furthermore, our semi-analytical model estimates the postbuckling stiffness based on a modified critical buckling stress, which introduces an exponential negative postbuckling stiffness resulting in a snapthrough structural behaviour. Through the analysis of the laminate as an orthotropic plate with a combination of aluminium and GFRP layers (the GFRP layers comprising of three plies with orientations [90 o /0 o /90 o ]) representing the Glare 4B-3/2 laminate, we are able to predict postbuckling stiffness and hence structural behaviour accurately using the modified semi-analytical model.

Semi-analytical model
The semi-analytical model developed for Glare plates consist of three different analyses stages: pre-buckling, buckling and postbuckling.
In many actual structural applications, a stable postbuckling path is favoured with unstable equilibrium paths often considered as unacceptable, especially when they are followed by dynamic snap-through. Structures exhibiting unstable buckling are generally imperfection sensitive, with their critical buckling load being very sensitive to even small changes in the amplitude and shape of initial imperfections as shown in Fig. 2. While for the perfect structure we generally see an unstable bifurcation point ''B'', the imperfect structure presents a limit point ''E'' [37,45] normally at a lower load . In this case the limit load − the load which causes the perfect structure to collapse is of less importance than the bifurcation load . In the particular case shown the structure is only mildly sensitive to initial imperfections as indicated by the fact that point E is only a little below point B, which pertains to the perfect shell A limit point is defined as a maximum or minimum on the loadcrosshead displacement curve at which the structural stiffness is zero. The maximum limit point seen in Fig. 3 corresponds to the point where a positive stiffness prior to buckling becomes a negative one during postbuckling. Further along the postbuckling path a second, minimum limit point is observed. Note, whilst the unstable path which is shown by the dotted red line is the one most likely to be followed in real life-structures, this will not be the case for theoretical models where instantaneous snap-through buckling will occur as shown by the solid blue curve.
Although the buckling of a perfect structure normally corresponds to a bifurcation point on the equilibrium path, this is not necessarily the case. Unstable buckling via a limit point can also occur if the prebuckling mode shape contains symmetries imposed by the boundary conditions and which therefore cannot be changed. For this reason, bifurcation or Euler buckling is sometimes described as symmetry breaking [46]. In such cases it is necessary to introduce a small initial imperfection when simulating a nonlinear analysis to break the symmetry. Were this imperfection not to be included, the solution may continue along the fundamental equilibrium path instead of branching onto the true postbuckling path at the bifurcation point (Fig. 4). Examples of where these problems might occur are straight columns and flat plates under axial compression.
For shell structure under axial compression, an instantaneous and dramatic failure will normally occur when the critical load limit is reached, with significant out-of-plane displacement where the structure 'snaps' as shown in (Fig. 4) [46] in either a positive or a negative direction. On each of the secondary paths, sudden high levels of out-ofplane displacement will be seen with the structure releasing membrane strain energy as kinetic energy. This has the potential to have a strong effect on structural stability as well as load-carrying capacity. In the postbuckling region, after a certain level of out-of-plane displacement has occurred, the structure begins to recover its load carrying capacity and can often be loaded to a higher level than its critical load value (Fig. 4). This is dependent however on the material retaining elastic properties which is often not the case for FMLs such as Glare where deformation in the metallic layers is often permanent due to plasticity.
In this work, an efficient semi-analytical model that was developed in [21] is applied to analyse the pre-buckling, buckling, and postbuckling behaviour of Glare fibre-metal laminated plates. This method is derived from a single variational formula, which is expressed in terms of Airy's stress function and transverse deflection function as, where is the inverse of the in-plane stiffness matrix, and is the bending stiffness matrix of the plate. ( ) are the in-plane force resultants of the plate expressed in terms of Airy's stress function , ( ) are the components of curvature in terms of the out-of-plane deflection , ′ ( ) are the midplane strains due to the deflection . 1 and 2 denote general forms for the prescribed displacement and stress boundary conditions, respectively.
To solve Eq. (1) using a standard Rayleigh-Ritz method, the Airy's stress function (x, y) and the transverse deflection function ( , ) are expanded into the following series forms [21], where ( ) , ( ) are the admissible functions for ( , ) that satisfy the essential boundary conditions for the out-of-plane displacement.
( ) , ( ) are the admissible functions for Airy's stress function and satisfy the stress-free boundary conditions. 0 ( , ) is introduced as an additional function or series to satisfy the prescribed boundary stresses.
After substituting Eqs. (2) and (3) into Eq. (1) and applying the Rayleigh-Ritz method, the nonlinear post-buckling behaviour of Glare plates is modelled by a set of nonlinear algebraic equations, as given by Eq. (22) in [21]. By eliminating the nonlinear terms for the in-plane equations, the model reduces to a pre-buckling analysis procedure for the Glare plate. If only the equations associated with the bending behaviour are considered, the model is then reduced to a standard eigenvalue problem for the buckling analysis of the Glare plate.

Plasticity localisation
The metal layers of the FMLs show significant plasticity under compression which must therefore be included in the model to accurately predict the buckling and postbuckling behaviour. Fig. 5 presents the effect of this plasticity in reducing the buckling stress to a plastic critical buckling stress ( ) . A correction factor for plastic buckling to represent the decreasing gradient of the stress-strain curve is employed in this study to provide an average compressive stress at buckling. Fig. 5 illustrates the tangent modulus which may replace Young's modulus E for thin plates, especially where through-thickness stress levels are less severe. Alternatively, the secant modulus gives provides the total strain at a reference point under stress as: The Ramberg-Osgood expression [47] of a stress-strain curve gives the total strain under a particular level of engineering stress as: where ( ∕ ) represents elastic strain, ( ∕ ) represents the plastic strain and , characterise the hardening behaviour of the material. Combining Eqs. (4) and (5) at the reference stress gives:  And differentiating Eq. (6) gives the gradient = ∕ at point : The tangent modulus E T is normally obtained from the material stressstrain curve, however, = ∕2 has been used in [48] for the buckling analysis of inelastic shell plates under uniaxial compression and is introduced here to simplify the semi-analytical solution. Substituting into Eq. (7) this gives: Substituting Eq. (8) into Eq. (5) and multiplying through by ∕ leads to a normalised stress-total strain relation: where and are material properties found from fitting Eq. (9) to a stress-strain curve for any metal layers. Furthermore, Fig. 5 shows that the elastic strain under critical elastic buckling stress ( ) is = ( ) ∕ . Substituting this value for in Eq. (9) allows = ( ) to be found. Design data [49] adopts a graphical solution to ( ) , employing a plasticity reduction factor < 1: Hence, the critical plastic buckling stress is evaluated as: This elastic-plastic buckling model for plates under compression loading was originally proposed by Wittenberg and de Jonge [32], and has been used by many researchers to predict inelastic buckling behaviour for isotropic and orthotropic shell structures such as Glare laminates [50][51][52][53]. This phenomenological model is based on comprehensive testing programmes and has been validated for inelastic buckling on orthotropic shell plate structures. Here, a plasticity correction factor based on the formula derived by Gerard and Becker [33] is introduced: = 1 + . and n is the Ramberg-Osgoodhardening parameter.  The instantaneous elasto-plastic Poisson's ratio is calculated according to Stowell& Pride [54,55] as: where is the fully plastic value of Poisson's ratio and has a value of a half for isotropic materials and is the elastic value of Poisson's ratio. The reduction coefficient represents the effect of inelastic behaviour [33].
The above analysis, programmed in MATLAB, takes around 11 s to predict the buckling and postbuckling behaviour of the Glare laminate using a Dell Precision 7540 64-bit system with a 2.6 GHz processor, 1 TB of memory and a Windows 10 operating system. The same system was later used to perform the numerical analyses in ABAQUS.

Semi-analytical results
The prebuckling, buckling and postbuckling behaviour predicted by the semi-analytical model is illustrated in Fig. 6. As mentioned earlier, this model follows a snap-through path which begins by tracking a positive stiffness path up to the critical buckling load (maximum limit point) after which unstable post buckling collapse occurs. The unstable postbuckling path is predicted using an exponential relationship with a negative stiffness and a decay constant: where is the initial postbuckling stiffness and is an exponential decay constant. Shortly after this collapse the unstable postbuckling path recovers to a stable path after passing zero postbuckling stiffness at the minimum limit point (as discussed in Section 2.1). This path is however unlikely to occur in practice due to the onset of plasticity and other failure criteria not included in this model causing a more severe drop in load carrying capability and a failure to return to a stable postbuckling path. Results are therefore presented only up to the minimum limit point (Fig. 6 point 2).
To study the ability of this proposed approach to incorporate the effect of initial geometric imperfections, a sensitivity analysis was performed. Due to the efficiency of the method such a study is relatively quick to perform and provides useful insights into the effect of the imperfection on the buckling and postbuckling behaviour of the plate.
A range of different magnitudes of imperfection represented by non-dimensional amplitudes ( ∕ ) (0.025, 0.05, 0.075, 0.1, 0.125, 0.15,  0.175, 0.2); (where is the amplitude of the initial imperfection and is the thickness of the plate) were considered in the postbuckling analysis. The obtained results are shown in Fig. 7.
As expected, the results of this study show that an increase in the amplitude of the imperfection causes a reduction in both the buckling stiffness and the peak load, which reduced from 13.49 kN with an initial imperfection of 0.025 mm to 8.52 kN when the imperfection amplitude was increased to 0.2 mm. This highlights the importance of using a representative, measured initial imperfection to obtain good predictions for buckling stiffness and critical buckling load.

Finite element mesh
A three-dimensional, ply-by-ply FE model was generated in Abaqus/ CAE. The geometry and thickness of each layer were extracted from high-resolution scans of real specimens with a high-accuracy structural mesh developed to represent the internal features. The resulting mesh is shown in Fig. 8. The layers of aluminium were meshed using linear 'continuum solid' elements (C3D8R) with the GFRP-metal interfaces   [30]). Three dimensional 'continuum shell' elements (SC8R), each consisting of 8 nodes with 3 degrees of freedom (DoF) per node, were used for the GFRP plies to enable the use of Hashin theory for modelling damage in the composite plies. Each layer was meshed individually and then assembled using 'TIE 'constraints between adjacent layers connecting all nodal degrees of freedom at the interfaces. A 1 mm mesh size was used for all constituents and cohesive elements with one element through the thickness following a mesh convergence check.

Material properties
An extensive literature review was implemented on the mechanical properties of the Glare ® material constituents and interfaces. Elastic properties for the metal layers represented by aluminium alloy 2024-T3 were taken from [37] and are summarised in Table 2, with plastic properties obtained from [56] being presented in Table 3. The mechanical properties and damage characteristics for the S2-glass fibre/FM94epoxy GFRP material were given in [57] and are shown in Table 4. Interfaces were simulated using Cohesive Zone Model (CZM) with bilinear traction-separation relations 'defined independently for modes I and II (further details are discussed in Section 3.2.6). The cohesive properties for GFRP-metal interfaces were extracted from [58] and are presented in Table 5.
Finally, the mode I and mode II cohesive stiffnesses and , respectively, were calculated based on the elastic properties of bulk FM94 resin [60], assuming that the interfacial stiffness is dominated by the deformation of a 10 μm thick epoxy resin layer (more details of the cohesive properties are provided in [30]).  Table 4 Mechanical properties for S2-glass fibre/FM94 epoxy prepreg material [57].

Loading and boundary conditions
Boundary conditions were applied at the top and bottom edges of the specimen to simulate the boundary conditions of the compression rig (Fig. 9). The bottom edge was constrained in all directions to represent the built-in end condition, while the top edge was free to move axially but restricted in all other directions. Both sides of the specimens were left free in all directions. Following this, a compressive Table 5 Cohesive zone properties for the GFRP-metal interfaces [30,[60][61][62].  axial load under velocity control was defined at the top edge of the model.

Geometric imperfections
In order to represent 'as-built' structures geometric imperfections were introduced. Since their exact form was not known they were modelled in the form of the first buckling mode shape (Fig. 10) with an amplitude scaled to give a maximum representation of the deviations measured in the specimens (due to the manufacturing process).

Load eccentricity
In addition to geometric imperfections introduced into the FE model, load eccentricity caused by misalignments in clamping the specimen leading to an asymmetrical load distribution was also considered. This was achieved by implementing two load steps in Abaqus/Explicit analysis as shown in Fig. 9. In the first step, the load eccentricity was represented in the form of an asymmetric linearly distributed crosshead displacement applied on the top edge. In the second step, a uniform compression load was imposed up to the final failure.

Cohesive zone model
One of the greatest difficulties of modelling the behaviour of FMLs is the complexity of the damage mechanisms compared to the modelling of metallic or composite laminates. In order to accurately model the initiation and propagation of damage in the specimens, a range of damage and fracture criteria were therefore introduced. The mixedmode bi-linear traction-separation cohesive zone model (CZM) was used to model delamination initiation and evolution in the GFRP-metal interfaces. This model uses quadratical nominal stress formula to simulate the damage initiation. Evolution of damage is then predicted based on the strain energy release rates for both mode I and II. The interfacial mechanical properties required to define the CZM are provided in Table 5.

Continuum damage and plasticity
Plastic deformation of the aluminium sheets was represented using the stress-strain properties of aluminium given in Table 3. A ductile damage theory was implemented to model damage in the aluminium layers [67]. This theory presumes the equivalent plastic strain at the onset of damage as a function of strain rate and stress triaxiality. The elastic mechanical properties used for the aluminium are shown in Table 2. The Hashin damage criterion [68] was used to model damage in the GFRP layers. The model available in Abaqus for continuum shell elements is a 2D version of the original Hashin criteria with four damage variables -for both compressive and tensile failures along and transverse to the fibre direction (for simplicity the former is defined as 'fibre failure' and the latter is termed as 'matrix failure'). The damage initiation and evolution variables for each damage mode of the composite material follows a bilinear traction-separation relation more details of which are described in [68].

Eigenvalue analyses
Linear eigenvalue-buckling analyses were implemented to provide estimates of the compression load and to obtain the eigenmode shapes required for modelling geometric imperfections. Abaqus/Standard (2019) was used for the eigenvalue analysis [67]. Since the Lanczos solver is generally faster for a system with multi degrees of freedom it was utilised in this study. Scaled eigenmodes based on the amplitude of imperfections measured in the test specimens were used to model imperfections, providing a conservative technique when the exact geometry of the imperfections is unknown [67].

Explicit dynamic analysis
Buckling and postbuckling behaviour were performed in Abaqus/ Explicit (version 2019) solver [67] incorporating the different factors that representing various forms of nonlinearity, including plasticity, damage, and crack propagation.

Sensitivity analysis for FE mesh
To examine mesh sensitivity in terms of its effect on predicting buckling behaviour, the Glare specimen described in Section 3.2.2 was studied. An analysis was implemented to compare axial load versus in-plane displacement curves for Glare ® 4B specimen models with different mesh densities. The FE model incorporated three layers of aluminium and two layers of GFRP representing the Glare ® 4B-3/2 standard Glare specimen. Cohesive layers with a constant thickness of 0.01 mm were inserted at each metal-fibre interface. Further details in relation to material properties and element types can be found in Section 3.2.2. Different element lengths (2 mm, 1.5 mm, 1 mm and 0.5 mm) were chosen to mesh the continuum and cohesive layers. Fig. 11 shows load versus in-plane displacement curves for the four different mesh sizes under compressive loading for the Glare ® 4B model. The buckling and postbuckling results in Fig. 11 show identical behaviour in the purely elastic region until buckling occurs. Following this, in the postbuckling region from approximately 7.5 kN axial load onwards, for mesh sizes 2 mm and 1.5 mm the load starts to be overestimated in comparison to the results obtained using a size of 1 mm. For both 1 mm and 0.5 mm mesh sizes load versus in-plane displacement curves show similar behaviour, therefore a mesh size 1 mm was chosen as the optimum value.

Finite element results
The buckling and postbuckling behaviours predicted using FE analysis will be compared with those from the semi-analytical approach and experimental work in Sections 3.1 and 4, respectively. At this stage however, it is worth considering the insight they provide.
In terms of damage initiation Fig. 12 highlights damage due to fibre tension ( ) in the middle section of the plate and fibre compression ( ) nearthe upper and lower grips, where local curvatures were high, while only partial matrix tension ( ) and matrix compression ( ) damage initiations are predicted in the top and bottom ends of the free edges of the plate.
Plastic strain contours at a crosshead displacement x = 1.175 mm in the postbuckling region are presented in Fig. 13. The results show  The total time taken to complete the full FE analysis of the buckling and postbuckling of the specimens was over 40 h on the same computer described in Section 3.1.1.

Specimen design
As discussed in Section 2.2 specimens measured 140 × 80 mm with a 100 × 80 mm unsupported section between clamping areas and a total thickness of 2 mm as shown in Fig. 15. Each specimen had a layup corresponding to 'Glare ® 4B'. They were manufactured by Airbus Germany GmbH from aluminium alloy 2024-T3 sheets with 0.4 mm thickness and Hexcel S2-glass fibre/FM94 epoxy unidirectional prepreg.

Test setup
A specially designed rig was used with a Dartec ® servo-hydraulic testing machine (equipped with a 500 kN load cell) as shown in Fig. 16.
Tests were conducted under displacement control with loading rate of 0.001 mm s −1 . (Rig design is discussed in further detail in [29]). Whilst only two specimens were able to be supplied by the industrial partner, the results for example those seen in Fig. 18 confirm the high level of consistency in both the specimens and the testing procedure providing confidence that the results obtained are representative.

Digital Image Correlation (DIC)
Specimens were monitored using a Dantec™ Dynamics Q-400 system to record three dimensional (3D) measurements of full-field displacements and strains to investigate the buckling and postbuckling behaviour. Specimens were prepared by applying a sprayed speckle pattern to white primed surfaces to enhance contrast. Two 2/3-inch greyscale CCD Limess™ sensors, each with a resolution of 1600 × 1200 pixels, were used with two lenses of focal length 28 mm to achieve a working distance of 300 mm to 600 mm. A HiLis™ monochromatic LED system was used as a high-quality light source. DIC images were recorded manually every 0.5 kN step up to the end of the test and postprocessed using the ISTRA™ 4D software with subset size 17 pixels and spatial resolution 0.2 mm.

Acoustic Emission (AE)
Four (Nano-30) acoustic emission sensors manufactured by Mistras™ Group were used to localise and detect the initiation and evolution of damage. These sensors have a frequency response range between 125-750 kHz and diameter of 8 mm. They were bonded to the specimen using multi-purpose silicone sealant as shown in Fig. 17.
The sensors were connected to a Mistras™ Group PCI2 acquisition system with a 45 dB threshold (as recommended by previous studies including Büyüköztürk and Taşdemir [69], Pearson [70] and Al-jumaili [71] for a wideband differential transducer) with a sampling rate of 5 MHz, recording over 1.2 ms to enable capture of full waveforms, using pre-amplifiers with a 40 dB gain and a built-in band pass filter with a frequency range of 20-1200 kHz.
Event locations were calculated using a bespoke location algorithm called 'Delta-T Mapping' [72]. This technique was developed to improve location in comparison with commonly adopted Time of Arrival (ToA) techniques which are based on triangulation and assume constant wave velocity in different directions, and direct paths between acoustic source and sensors. For composite structures, material anisotropy causes great variation in wave velocity with respect to in-plane direction. The Delta-T technique overcomes this limitation by mapping the structure and then using these maps in the location of any AE-generating event (impact, crack evolution). The sensitivity and accuracy of the method make it particularly useful for the detection and localisation of internal damage in FML structures, which contain anisotropic composite laminae [73]. A more comprehensive discussion of the Delta-T algorithm can be found in [72].

Experimental results
Specimens buckled with one half wavelength along the axial direction. The full field of the out-of-plane displacement is shown in Fig. 18, in which similar results are obtained for each specimen. Axial load versus in-plane displacement curves for both specimens are presented in Fig. 19. Excellent agreement is observed between the two specimens in terms of pre-buckling stiffness with a slightly lower ultimate compressive strength of 10.60 kN for specimen (2) when compared with specimen (1) which failed at 11.18 kN. Postbuckling failure was consistent in the tested specimens although specimen (2) presents slightly lower postbuckling stiffness compared with specimen (1).
AE results for the two specimens show similar behaviour with those for specimen (1) shown in Fig. 20. A large jump in cumulative energy is observed at the maximum compressive load (11.18 kN) followed by a steadily increase in energy throughout the postbuckling region due to increasing levels of damage such as matrix cracking. The AE location results show damage is concentrated in areas subjected to the highest curvature, with intensity increasing as the test progresses. A higher level of events towards the left-hand side of the specimen correspond  to the slightly asymmetric buckling mode seen in Fig. 18 resulting in higher levels of displacement towards the left-hand edge.
The location of the damage correlates well with that predicted using the Hashin criteria with a high level of activity in centre as highlighted by the damage initiation criteria and towards the upper and lower grips, where local curvatures were high (Fig. 12). As a consequence, it confirms the capability and accuracy of the Hashin damage criterion to predict damage for Glare structures in both buckling and postbuckling regimes. Fig. 21 compares the load-displacement path predicted by the semi-analytical approach with that from the FE analysis and the experimental results.

Results and discussions (validation of the semi-analytical and FE models)
An excellent correlation is seen between the analytical results and experiments in terms of the buckling behaviour, ultimate compressive strength, snap-through and postbuckling behaviour of the tested specimens. For the FE model, specimens initially buckled at around 7.5 kN which is in excellent agreement with experimental results. The ultimate compressive strength of 10.54 kN from the FE model also shows excellent agreement with specimen (2) in particular although it slightly underestimates the strength of specimen (1). Whilst the FE model predicts a slightly higher pre-buckling stiffness than that found experimentally this is believed to be as a result of using shell elements when modelling damage in the composite layers which neglect the through thickness stresses [74], leading to an overestimation of their stiffnesses. A greater difference however is seen in the results for the postbuckling region where the model again overestimated the experimentally determined stiffness. This might be as a result of unstable behaviour in the test structures leading to collapse, not represented in the FE model where calculations are based on the central-difference integration method giving a more stable postbuckling analysis. In contrast to the FEA analysis described in Section 3.2 which took more than 40 h to produce similar fidelity results, a semi-analytical approach enables us to simulate the optimum buckling and postbuckling design of Glare laminate structures in only a few seconds.
In terms of failure modes, Fig. 22 shows the out-of-plane displacement behaviour from the FE models alongside experimental results from the DIC system. The contour plots present the out-of-plane deformation at initial buckling, peak load and postbuckling regimes, with corresponding in-plane displacements. The plate is seen to buckle with a single half wavelength in the loading direction as expected for a plate with built-in ends and free longitudinal edges subjected to axial compression. Deformations occur in the middle of the specimen as would be expected, however, the buckle initiates and becomes more pronounced towards the left-hand side of the specimen due to the load eccentricity found in the experiments. Good agreement is noted between experimental and FE results in capturing the mode shape, although the FE model overestimates the out-of-plane displacement slightly.

Conclusions
A semi-analytical approach has been developed to simulate the buckling and postbuckling behaviour of fibre-metal laminate specimens (Glare ® 4B-3/2) subjected to compressive loading. Results have been validated against both an FE model incorporating geometric imperfections, load eccentricity and the full range of damage mechanisms and experimental results.
Good agreement was observed between axial loads, out-of-plane displacements and event location using AE which has been shown to be particularly sensitive to early damage activity in FML structures.
The FE model developed has been demonstrated to provide an accurate prediction of FML behaviour based on measured geometrical imperfections and load eccentricities. A combination of damage models was able to determine a lack of delamination, as well as plastic deformation of the metallic sheets in areas of high curvature. Damage in composite layers was predicted and found to be due mostly to fibre compression and in-plane shear. These results are validated by the AE data from the experiments.
Finally, the semi-analytical model results predicted using the Rayleigh-Ritz method and the Ramberg-Osgood formula showed an excellent prediction of the overall buckling and postbuckling behaviour of FML structures using 'as measured' values of initial geometrical imperfection. This was particularly apparent in terms of pre-buckling stiffness, buckling load, ultimate compressive strength and initial postbuckling behaviour. In addition to its excellent ability to predict this behaviour the semi-analytical analysis was found to be computationally much less expensive compared to FEM. While the FE analysis takes more than 40 h to produce buckling and postbuckling failure analysis results, the same analysis can be implemented in only 11 s using the semi-analytical model, which enables us to obtain optimum buckling and postbuckling designs for real structures using a semi-analytical approach in hours, whilst a similar FE analysis would take weeks.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.