Anumerical study on the influence of compositewrinkle defect geometry on compressive strength

• A newmethodology is presented to rapidly create finite element models of wrinkled composites based on NDT in-

A numerical study on the influence of composite wrinkle defect geometry on compressive strength Ningbo

H I G H L I G H T S
• A new methodology is presented to rapidly create finite element models of wrinkled composites based on NDT information. • The suggested approach allowed a rigorous study to rank wrinkle-parameter importance using a simulator of NDT information. • Maximum angle in load direction was the key parameter affecting compressive strength, with wavelength being of secondary. • Increasing wrinkled region extent perpendicular to load direction also correlated with decreasing compressive strength. E-mail address: ningbo.xie@bristol.ac.uk (N. Xie). [1,2]. To add further complexity to the prediction of failure, variations can occur in the material internal meso-structure as a result of manufacturing processes. Out-of-plane misalignment of the fibre paths, generally known as 'wrinkling', is a relatively common phenomenon in thick components or curved sections [3] and is known to have a significant influence on mechanical performance. Wang et al. [4] manufactured out-of-plane fibre waviness in AS4/8552 (resin/fibre) composites by three methods and, based on their experimental work, Lemanski et al. [5] developed FE models to simulate the compressive failure process. Both experimental results and modelled predictions indicated a reduction in compressive strength of about 54% when peak misalignment angle increased to about 22°. Mukhopadhyay et al. [6] also showed the significance of wrinkles on compressive strength. They introduced wrinkles by selectively inserting or removing prepreg strips in the 90°plies of quasi-isotropic laminates. Both experimental results and numerical predictions showed that a maximum misalignment angle of 12°resulted in a 30% knock-down of compressive strength. This reduction of compressive performance related to wrinkle angles has also been investigated experimentally and numerically by Ferreira et al. [7]. Other research has identified this negative influence and summarised the dependence of the compressive-strength knockdown on amplitude [8], or the ratio of amplitude-to-wavelength [9,10]. Therefore, during the design, test, simulation and performanceassessment stages for a composite component, the potential negative influence of any wrinkles needs to be carefully considered. It is important to have a full understanding of how wrinkles influence the compressive failure of laminates, in order to understand which wrinkle parameters (severity, shape or extent) play a more important role. This information can then be used to determine the metrics that should be measured non-destructively when wrinkles occur during manufacture. A number of authors have attempted to characterise the fibre paths and variations from ideal trajectories in composite materials using nondestructive testing (NDT) methods [11][12][13][14][15][16][17][18], covering the eddy current [11,13], X-ray computed tomography (CT) [12] and ultrasonic techniques [14][15][16][17][18]. The development of this field using eddy currents is currently limited to the two dimensional (2D) response at the surface, so it is not currently capable of full 3D characterisation. Hence the latter two methods are more likely to be used in the industry. For the X-ray CT technique, its 3D data has been used to assist the creation of numerical models for estimating the influence of detected defects on structural integrity. Alghamdi et al. [19] successfully created FE models based on Xray CT data, but their models did not include any failure mechanisms. The group of Makeev et al. has made a significant contribution to this field [20][21][22]. In their most recent research [22], the in-plane and outof-plane fibre waviness was characterised based on the X-ray CT technique and was transferred into FE models with modified LaRC04 [23] failure criteria included, the predicted results showed good agreement with the experimental results. Ultrasonic characterisation has been more widely used in the component-quality evaluation stage as it has greater capacity for large-size components than X-ray CT. Freemantle et al. [24] developed FE models from ultrasonic NDT data for composites containing delaminations after suffering an impact load. Their model predictions showed good agreement with the test results for flat panels but their method did not include wrinkle characterisation or modelling of any fibre or matrix damage mechanisms. Sandhu et al. [25] characterised the fibre path for wrinkles, using 'multiple field image analysis' on ultrasonic B-scan data to determine wrinkle amplitude, wavelength and location. They proposed that this 2D information could be used with a prismatic assumption to create a 3D FE mesh through the extrusion of the 2D geometries but no results from this method were reported. Smith et al. [14][15][16][17] demonstrated inversion of ultrasonic full-waveform scans to obtain three-dimensional (3D) maps for both out-of-plane ply wrinkling and in-plane fibre waviness. Their method is the basis on which numerical models are created in this paper. In their research [26], they pointed out the measurement capability for wrinkle angle at every location could be possible up to a wrinkle angle of 15°, and the resolution for a complete map of fibre angles is dependent on scan pitch, with a minimum of around 0.2 mm. This stateof-the-art capability has not been available in industry, but is being transferred into commercial software applications. The development of such non-destructive methods raises the necessity for measurement guidance for manufacturing imperfections.

G R A P H I C A L A B S T R
Composite components with out-of-plane wrinkles under compressive loading have been shown by the authors to fail due to a combination of delamination, matrix cracking and fibre failure [6]. There have been many theories proposed for the prediction of onset and propagation of each mode of failure for embedding in numerical models. The following failure-mode theories, validated by the authors [6], are those used in the present work reported in this paper. For the initiation of delamination, a quadratic damage-initiation criterion [27] was applied and then a Power-Law [28] criterion was used to govern the mixedmode interface failure. Matrix cracking is controlled by the criterion suggested by Puck et al. [29,30], with the assumption that a crack happens on a plane inclined at a specific angle. Finally for the fibre damage, the kinking model is based on Pinho's work [31,32], which suggests that the fibre failure occurs on a rotated plane, triggered by shear stress due to the initial fibre misalignment and the final region of fibre failure forms a kink band. The authors [6] combined the failure mechanisms stated above in finite-element models for wrinkled composite laminates and demonstrated that the resulting models possessed the capability to predict the whole failure process in detail.
Even though the effect of wrinkling on composite strength has been shown to be significant, there is still no consensus on which parameters, or combinations thereof, have the most significant influence on failure stress. Amplitude (A), wrinkle wavelength (L) and maximum angle (θ), defined in Fig. 1(a), have each been proposed by different authors to describe a wrinkle. Lemanski et al. [5], Mukhopadhyay et al. [6] and Sutcliffe et al. [12] chose maximum angle to define wrinkle severity, while Fedulov et al. [33] used wrinkle height as a fraction of laminate thickness. Hsiao and Daniel [10] used amplitude and wavelength to characterise wrinkles, assuming a sinusoidal wrinkle shape, and proposed a representative volume truncated at a single period of the wrinkle, with amplitude reducing linearly from the mid-plane to the sample surfaces. Later Caiazzo et al. [34] used a polynomial to describe the wrinkle shape, with amplitude reducing linearly to the sample surfaces, and used the peak height and the wrinkle extent in the load direction as the 'gross' measures of the defect size. More recently, El-Hajjar and Petersen [35] proposed a Gaussian function to characterise 'bell-curve' wrinkles, the maximum amplitude (waviness height) at one surface was designed to be diminished at the opposite one, thus forming samples with one flat surface and one concave 'bell' surface, and the wrinkleamptitude reduction through the thickness related to the distance from the maximum-amplitude surface. In the results, they chose 'waviness height' to describe the wrinkle severity or extent. The application of Gaussian functions to describe the wrinkle distribution was adopted in this paper, but the range of wrinkle metrics investigated is significantly greater.
This review of the literature shows a lack of consensus over the key wrinkle metrics for non-destructive measurement. It illustrates the need for a rigorous and controlled multi-dimensional parametric study to identify and understand the interdependencies of parameters and determine a hierarchy of wrinkle-parameter importance, which could provide guidance for the industrial quality-assurance process for components with manufacturing induced wrinkles. The objective of this paper is to use a validated modelling approach to explore the parameter space in terms of a wrinkled component's compressive strength. This could not be achieved experimentally in a systematic and rigorous way due to the very large number of highly controlled specimens that would be required and the difficulty of precisely controlling specific parameters of wrinkles during the manufacturing process; nor has it been possible in the past using numerical models because of the time and manual effort required to create each different model. This aim has now been realized by a methodology involving a novel combination of: 1) a custom MATLAB wrinkle simulator, which creates material-property maps embodying the wrinkle geometry based on a few input parameters defining wrinkle severity, shape and extent; 2) a novel transfer process automatically creating FE models based on the simulator's outputs and some validated assumptions, originally designed for NDT-based performance modelling; 3) a powerful numerical analysis using validated failure mechanisms. The numerical modelling techniques and experimental validation [6] that form the basis of the rigorous study in this paper were focused on pre-preg laminates with flat top and bottom surfaces, representing the case of an internal wrinkle in a component manufactured by 'rigid' tooling. The work presented here has thus also been limited to flat laminates manufactured using a hand layup technique.
The material-property maps can be generated in various ways including from inverted NDT data, from inverted micro-sectional analysis or, as in this case, from a custom MATLAB-based wrinkle simulator, in which the simulated maps were based on the coupons used for the validation. Section 2 defines the wrinkle topology used in this paper, and Section 3 explains this novel combination of the wrinkle simulator, the numerical analysis and the transfer process to automatically create large numbers of FE models with controlled variations of specific parameters. The validation is described using coupon-based models, controlled by the simulator and the transfer process, against the same experimental test data reported previously [6]. Section 4 contains the results of the simulations created to realize the main objective of this papera parametric study of the knock-down in compressive strength due to wrinkled load-direction fibres. In Section 5, the wrinkle orientation was rotated by 90°relative to the load direction, for further analysis. Finally, a discussion is presented in Section 6 and conclusions are given in Section 7.

Wrinkle topology definition
There are a variety of origins for out-of-plane wrinkles during the manufacturing process. Potter et al. [3,36] categorised them into four groups, with the wrapping of prepreg during storage and transportation, friction between the tool and composite during autoclave cure and residual thermal stress being the major sources causing local fibre misalignment in flat laminates. The wrinkle topology defined here aims to reflect the range of misaligned fibres and fibre undulations that can occur in the out-of-plane direction, with some idealisations in terms of shape and distribution for ease of analysis and quantification.
Wrinkled fibres in both the 0-degree (x) and 90-degree (y) directions are considered in the paper. The wrinkle topology defined along the x direction (0-degree fibre) is expressed as a Gaussian-modulated cosine-wave profile. Fig. 1(a) illustrates the cosine profile, the wrinkle displacement (fibre misaligned undulation), d, is governed by a generic equation for a 3D wrinkle in Eq. (1). A wrinkle in the x direction is described by six parameters (noted in Fig. 1): Amplitude, A; Wavelength, L; Wrinkle Gaussian Half-Width, w 1 ; maximum wrinkle angle, θ; wrinkle-centre location, x 0 and an x-offset phase parameter, x 1 , which can convert this to a sine wave if x 1 ¼ L 4 . The maximum angle, θ, is a non-independent parameter that is determined by the other parameters. w 1 describes the wrinkle extent in the x direction. For the y (width) and z (thickness) directions, the profile of wrinkle amplitude, A, from its value A 0 at the wrinkle's centre, (x 0 , y 0 , i 0 ), can also be controlled by a Gaussian function requiring two more extent parameters ( Fig. 1(b)), the Transverse Wrinkle Half-Width, w 2 and Wrinkle Half-Height, n, as shown in the definition of A in Eq. (2), where i and n are in units of plies. An additional non-Gaussian step-reduction profile for wrinkle amplitude distribution in the z direction was required to better represent wrinkles that were embedded in the experimental validation samples, and this specific distribution is covered in Section 3, Section 4.1 and Section 4.2.

Simulator, transfer process and validation
The transfer process, first presented by the authors in [37], can be used to create finite-element models from 3D material-property maps. These can be derived either from real NDT data or from a custom MATLAB-based simulator developed for this work. The simulator constructs a specific wrinkle geometry for a coupon based on userspecified values for certain wrinkle parameters: • maximum amplitude (in the centre of the coupon), A • wavelength, L • shape (sine or cosine) • depth-dependence of amplitude (Gaussian or customised) • orientation relative to load direction • extent (Gaussian 1/e half-width) in all three dimensions, w 1 , w 2 and n (plies) as well as parameters that describe the size and in-plane pitch of the material-property maps to simulate realistic NDT-based scans. The two maps generated are 1) a ply-interface depth map defining the local depth of each ply interface (resin-rich layer between plies) and 2) an in-plane orientation map containing the local in-plane fibre-tow orientation angles.
The MATLAB-based transfer process is applied to the two materialproperty maps, to automatically transform them into a finite-element model to be run in the Abaqus/Explicit package. Mesh geometry and out-of-plane ply angles are obtained in the transfer process from the ply-interface depth map. Cohesive elements, to allow the delamination failure mode, are inserted at the interfaces between plies and these have a thickness approximately equivalent to the resin-layer thickness between plies. Fibre orientation is stored per element in a distribution table and is calculated based on a fibre-orientation vector derived from the two material-property maps. The transfer process allows selection of either a uniform fibre volume fraction (FVF) for the whole model, or a ply-thickness dependence of FVF. In this paper, all models were created with a constant FVF, since during cure the tendency is for the FVF to become as uniform as possible. For compatibility with previous work by the authors, the fibre kinking, matrix cracking and delamination failure modes were governed by the appropriate equations implemented in VUMAT [6].
The FE models developed by the transfer process were each solved in Abaqus under compression loading which is parallel to the 0-degree fibre direction. The whole process is described schematically in Fig. 2.
The Matlab-based simulator was applied to search for the best Gaussian-envelope cosine-shape fit to wrinkles in images of the edge of the experimental coupons used in [6] (Fig. 3(a)), as well as an optimised depth dependence of the wrinkle amplitude to match the unusual 90-degree ply thicknesses in these coupons. The two materialproperty maps were simulated for the amplitude, wavelength and envelope-width corresponding to this best fit. The coupons had 24 plies with lay-up [45/90/− 45/0] 3S to form a quasi-isotropic ply sequence, using the IM7/8552 [38] fibre/resin system, as explained in [6], with extra 90-degree strips of different widths being inserted in, or removed from, the 90-degree plies to form the artificial wrinkled region with a designed wrinkle severity. The wrinkles appear in the other plies by forming groups of 3 plies: [−45/0/45] (but 4 plies at the midplane: [−45/0/0/−45]), separated by 90-degree plies with the thickness variations. The wrinkled region in Fig. 3(a) shows a reducing amplitude from the central plane to the top/bottom surfaces with ratios of 1.0: 0.63: 0.39: 0.0 for the ply groups separated by 90-degree plies. The simulator uses these ratios to create the wrinkles in the materialproperty files from which the models are created by the transfer process. By comparison of the overall geometrical shape of the wrinkled region, photographed in Fig. 3(a), with the finite-element mesh in Fig. 4, where the measured maximum angles at the mid-plane (determined by the software 'ImageJ') were 11.7°(Coupon edge) and 11.4°(FE mesh) respectively, it is clear that the characteristics of the wrinkled region have been successfully captured and transferred from the experimental coupons into the FE models.
The wrinkle displacement, d, at the interface between two plies, was governed by Eq. (3) in the x (0°fibre ply) direction, and the wrinkle amplitude in the y (90°fibre ply) direction was uniform. The parameters of Eq. (3) are the same as in Eq. (1), except that A i is defined to represent the depth-dependent maximum displacement and subscript i denotes the interface index.
In the thickness (z) direction, the amplitude A i reduced stepwise from a maximum at the central interface to zero at the surfaces, following the ratios: 1.0: 0.63: 0.39: 0.0, which reflected the real distribution of the wrinkle amplitudes in the mechanical-test coupons [6].
Mukhopadhyay et al. [6] manufactured specimens with three different severity levels (maximum wrinkle angles of 6°, 10°and 12°) to study the wrinkles' influence on the compressive failure process. For each wrinkle severity level, three coupons were chosen to validate the capability of the prediction by models developed by the transfer process, two of each severity (sample numbers 1 and 8) were from edges of the 8-coupon panel and one (number 5) was chosen from the middle of the panel. Material-property maps for these samples were simulated and passed to the transfer process as indicted by the green line in the workflow shown in Fig. 2. The transfer process was then used to create finite-element models for strength predictions, using the actual specimen width of 30 mm. Comparisons of predictions by transfer-process models with the mechanical test results are shown in Fig. 5 for each coupon. Also shown in Fig. 5 are results from specimens with a reduced width of 20 mm following the same procedures as full-size models from the image geometric measurement to final model computation. In the 20 mm-width models only the most severely wrinkled central two 0°plies had the fibre-kink failure model activated. Both the reduction in model size and the confinement of the fibre-kink failure to the central plies were aimed at reducing computation time, which is of benefit in the subsequent parametric study involving a large number of models.
For both model sizes, similar failure stress levels were obtained and both predictions are in good agreement with the experiment for both types of model. Therefore 20 mm-width models with the fibre-kink failure mode only in the central two 0°plies gives acceptable predictions for compressive failure stress influenced by wrinkles, and this method was applied in the following parametric studies in Section 4 and Section 5.
On closer examination of the results it can be noted that, for both model widths, the predictions of the edge-coupon strengths (no. 1 Fig. 2. Flow diagram expressing the process from NDT data or simulated NDT data into FE model. Use of the photographic section of a wrinkle for simulation is optional. and 8) are not as good as for middle coupons (no. 5) for each wrinkleseverity case. This may originate from the local variation in the wrinkle shape within the edge regions. Further investigations will be conducted to study this issue by creating models directly from inverted NDT data on specific coupons, rather than simulated 'idealised' data. However, considering the standard deviation in the mean strength of each wrinkle severity, indicated by black error bars in Fig. 5, the edge-coupon predictions are still acceptable compared with the overall test results.

Parametric study
As noted, the 'severity' of a wrinkle has a significant influence on compressive strength, but from the various studies in the literature it is not clear precisely which parameters for severity are most influential on the strength and should be measured non-destructively. Smith et al. [14] highlighted the range of parameters, such as wrinkle shape, number of wrinkle cycles in a given wrinkled region, percentage of wrinkled fibres, maximum deviation of fibres, etc. that may also have influential effects. It is important to understand how wrinkle severity, shape and extent influence the mechanical performance. The new modelling technique described above provides the capability to investigate this in a systematic and rigorous manner, as is described in this section.
Wrinkle shape can vary depending on the cause of the wrinkle in the production process. For the purposes of this study, only cosine-wave and sine-wave shapes are considered. For these two shapes, Eq. (4) shows that the out-of-plane wrinkle severity can be defined in the load (x) direction by an envelope parameter (w 1 ), an amplitude (A i ) and a wavelength (L). The maximum misalignment angle (θ) is linked to these parameters geometrically ( Fig. 1(a)) but it is not an independent variable. The first study, reported in Section 4.1, is an investigation of the influence of maximum amplitude, wavelength and maximum angle in the load (x) direction on compressive failure stress. A comparison of the cosine and sine shapes is then provided in Section 4.2. A study of wrinkle-extent parameters in the x, y and z directions is reported in Section 4.3. All the models created for the parametric studies had 24 0.25-mm thick plies with a stacking sequence of [45°/90°/ − 45°/ 0°] 3S and the modelled part of the coupon was 20 mm by 20 mm.

Wrinkle severity -angle, amplitude and wavelength (fixed wrinkle volume)
In this section, the displacement, d, of ply interfaces in the wrinkled region is defined along the x direction by the product of a cosine wave, with wavelength L, a Gaussian envelope with a 1/e half-width, w 1 and an amplitude A i (see Eq. (3)). The volume of the wrinkle envelope is fixed in these studies since the Wrinkle Gaussian Half-Width (w 1 ) remains unchanged at w 1 = 4 mm while the wrinkle extends uniformly across the width of the coupon (y direction)see Fig. 6, and the zdirection variation is fixed as described for each study.
The maximum wrinkle angle (θ), based on a cosine shape, can be calculated taking the arc-tangent of the maximum value of the derivative of Eq. (3). This derivative is given in Eq. (4).
It is clear that maximum angle is dependent on w 1 , A i and L, and is not an independent variable. Therefore, for a single interface, only the three parameters of w 1 , A i and L were changed independently. For a fixed amplitude (A i ) and Wrinkle Gaussian Half-Width (w 1 ), wrinkle shapes will change with different wavelengths (L). When the ratio L/ w 1 is less than approximately 1.5, the wrinkle shape tends to be repeated cosine cycles ( Fig. 6(a)) but it gradually changes into a single Gaussian-governed general shape when the ratio becomes larger than 4 ( Fig. 6(c)). Three examples are shown in Fig. 6.
In terms of the thickness direction (z-direction), two studies were undertaken, with the aforementioned two kinds of amplitudereduction ratios, A i , from the central interface to the top and bottom surfaces of the specimens, as shown in Fig. 7. These two types of models used an amplitude variation following a Gaussian reduction and a reduction in stepped blocks of four interfaces, which more closely follows the experimental geometry. In the Gaussian thickness profile, local amplitude A i of ply-interface i reduces continuously from central ply, i mid having the maximum amplitude A i mid , towards the upper and lower   Fig. 8 shows examples of these two amplitude profiles.
For all models, the wrinkle Gaussian Half-width (w 1 ) was fixed at 4 mm, while the wavelength (L) was changed from 6 mm to 16 mm and mid-ply amplitude (A i mid ) was changed from 0.2 mm to 0.6 mm, with the maximum misalignment angle (θ) being calculated in each case. Fig. 9 illustrates the dependence of modelled failure-stress knockdown (colour scale) and maximum angle (vertical axis) on the amplitude and wavelength for both types of model. The percentage knockdown in failure stress is defined as the modelled pristine value (647.0 MPa) minus the wrinkle predicted failure stress, and then divided by the pristine value. For both types of model, the knock-down in failure stress follows the same trend, being related to the maximum angle, rather than wavelength or mid-ply amplitude alone. This indicates that maximum angle is the major influence on failure stress, while wavelength and amplitude have an indirect influence due to their joint involvement in governing maximum angle. The failure stress knock-down, in both model types, varied by N 40% when the maximum angle varied from 4°to 28°, showing the importance of careful nondestructive angle measurement when assessing wrinkled-component  performance. Moreover, when comparing the distributions of failurestress knock-down between the two types of amplitude throughthickness profile ( Fig. 9(a) and Fig. 9(b)), it is clear that the distributions are quite similar, which suggests that the through-thickness profile has only a slight influence on failure stress, not as significant as the maximum angle. The same data is shown as contour plots in Fig. 10(a) and (b), from which three different regions can be identified (for w 1 = 4 mm), as follows.
i. Wavelengths N 7 mm (L/w 1 ≥ 2) with maximum angles N8°T he failure-stress knock-down (magenta solid lines) follows the trend of the maximum angle (black dashed lines) very closely when the wavelength is N 7 mm and maximum angle is N8°.
ii. Wavelengths b7 mm (L/w 1 b 2) However, for shorter wavelengths, the wavelength parameter appears to become more important. A possible explanation is that, below a wavelength of approximately 8 mm (L/w 1 b 2), there are multiple cycles in the wrinkle, as shown in Fig. 6(a), and this could be influential in introducing a wavelength dependence for smaller wavelengths. By contrast, when the wavelength increases to N12 mm (L/w 1 N 3) the wrinkle tends to be a single Gaussian-governed shape as shown in Fig. 6(c).
iii. Maximum angles b8°A t small maximum angles, there is evidence of a periodic effect, possibly linked to the ratio of wavelength to Gaussian half-width. A further study of low-angle wrinkles would be required to investigate this further.

Wrinkle shape (cosine to sine shape in a Gaussian envelope)
There are many options for shape investigations, but the change from a cosine shape to a sine shape was chosen for this study in order to move the maximum angle to the centre of the wrinkle, where the Gaussian envelope always has a value of unity (Fig. 11). This was in preparation for the wrinkle-extent study reported in Section 4.3. The experimental validation of the model used a cosine shape with maximum displacement in the wrinkle centre. The validated model can then be used to investigate whether a change in shape is likely to make a significant difference to the compression strength.
A series of models were created using the simulator and the transfer process to investigate if a sine shape gives similar failure-stress levels to the cosine, provided maximum angle remains the same. Seven maximum angles were chosen for comparison, using the cosine shape as the comparison baseline. While the main criterion was to keep maximum angle the same as for the baseline (cosine) shape, due to interactions between amplitude, wavelength and angle, there are three different methods of achieving this: 1) Keep angle exactly the same and change amplitude and wavelength 2) Keep angle and amplitude exactly the same and change wavelength 3) Keep angle and wavelength exactly the same and change amplitude Fig. 11 shows the shapes of the three sine curves and one baseline cosine curve (amplitude is 0.4 mm and wavelength is 6 mm), with a maximum angle of 20°, as an example. Both Stepped models and Gaussian models were studied and Fig. 12 shows the results for (a) stepped models and (b) Gaussian models. It can be seen that the difference in failure stress level when moving from cosine-shape to sine-shape is not significant if the maximum angle remains the same, for both Stepped and Gaussian models. According to the study in Section 4.1, angle and wavelength are the most important factors. Case 3) here, fixing these parameters gives similar predictions to the cosine baseline and so was used for the study in Section 4.3. Overall this provides sufficient evidence that the shape difference between sine and cosine is insignificant and that the sine shape is suitable for the following wrinkle-extent parameter studies.

Investigation of wrinkle extent
In the case of the cosine curve, wrinkle extent is controlled by the Gaussian Half-Width (w 1 ), which also affects the maximum angle. In this section, a sine curve is used to remove this relationship as the maximum angle of the sine curve depends only on amplitude and wavelength. The model-creation process and details are the same here as in  , was used to determine the wrinkle displacement, d.
where (x 0 , y 0 , i 0 ) is the centre of the wrinkle, at which the displacement is A 0 , and i is in units of the number of plies. This also allows a variation of the wrinkle parameters across the width of the specimen. As aforementioned, the amplitude envelope in the load (x), width (y) and thickness (z) directions are defined as Gaussian functions, so the shape of the wrinkled region described by any isoamplitude surface is thus an ellipsoid. As indicated in Fig. 13(c), the particular case of a 1/e isoamplitude surface is an ellipsoid where the semi-axes along the x, y and z directions are respectively: the Wrinkle Gaussian Half-Width (w 1 ), the Transverse Wrinkle Half-Width (w 2 ) and the Wrinkle Half-Height (n), where n is in units of the number of plies.
An amplitude, A 0 , of 0.35 mm and a wavelength, L, of 6.0 mm were fixed in the Gaussian models for the study of wrinkle extent in this section, resulting in a maximum angle of 20°for all cases. The rationale behind this choice was that the amplitude should be small enough to avoid unrealistically large changes in ply thickness, while the maximum angle needs to be large enough that the modelled knock-down effects are significant and could vary both up and down to determine the effect of the wrinkle extent. The use of Gaussian envelopes provides a continuous reduction of amplitude through the thickness, which is more realistic for naturally occurring manufacturing defects than the stepped method. It is straightforward to control the wrinkle extent in all three directions by changing w 1 , w 2 or n (see Fig. 13(c)).
The area of the coupon cross section within the 1/e isoamplitude surface, Area_wrinkle, encompasses all wrinkle-envelope locations greater than A 0 /e (i.e. 37% of A 0 ). Having explored different isoamplitude values to define the extent of the wrinkle, 1/e was chosen for convenience on the basis that different isoamplitude levels have no unexpected influence on the relationship between failure stress and cross-sectional wrinkled area, as shown in Fig. 14. The distributions of six different polynomial-fitted curves have similar trends and positions to expectations given that the area within the isoamplitude (see Eq. (7)) changes predictably.
The aim of this section is to investigate whether n, w 2 or the Area (%) is most significant in governing the knock-down of failure stress, since these parameters decide the wrinkle extent. Numerous models were created spanning a wide range of extents in all three dimensions. All models were simulated with a sine-wave shape along the load direction. The wrinkle Gaussian Half-Width (w 1 ) was varied between 2 mm and 10 mm. The Transverse Wrinkle Half-Width (w 2 ) was varied in the range 3 mm to 20 mm, and the wrinkle half-height (n) was varied between 1 and 5 ply thicknesses.
The first comparison shows the effects of cross-sectional dimensions on failure stress. Fig. 15(a) presents the effect on percentage knockdown of failure stress by changing w 2 and n when w 1 is fixed at 4 mm. It is clear that both w 2 and n have an influence on the trend of knockdown of failure stress since the colour map is not uniform in either the w 2 or n axes. However, it does not follow the cross-sectional area contours in Fig. 15(b), indicating a more complex relationship between w 2 , n and failure-stress knock-down. Hence, at this stage it can be concluded that both w 2 and n have an effect on the failure stress. The combination as a percentage of cross-sectional area also has some influence, the variation of which from 2% to 40%, caused N30% change in the failure stress knock-down.
Considering the poor fit of cross-sectional area and failure-stress contours in Fig. 15(b), a parameter, a, combining w 2 and n is suggested in Eq. (8) to better fit the dependence of failure stress on w 2 and n.
where N stands for the total number of plies, y width is the width of the coupon in the y direction, and parameters C and D adjust where the 'corners' occur in the contours of parameter a. C is most likely to depend on the aspect ratio of the coupon: y width /z thickness , although this has yet to be investigated. Similar plots to Fig. 15(b) but replacing cross-sectional area (%) with parameter a, are presented for C = 4 and D = 1 in Fig. 16. The magnitudes of C and D are chosen from several trials to obtain the best-fit. From Fig. 16(a), it is clear that, when plotting against parameter a, the contours of a follow those of failure stress better than cross-sectional area (%). While comparing Fig. 15(b) and Fig. 16(a), it is obvious that the two plots have similar colour-scale distributions, which means both the area and the smallest linear dimension in the cross section have effects on failure stress.
The n − w 2 space in Fig. 16(a) divides into two zones, shown in Fig. 16(b). The interface between the two zones is determined from the transition points of each contour curvei.e. where each contour transitions from horizontal to vertical. Two zones are thus formedone where the knock-down is independent of w 2 because the contours are approximately horizontal, and the other independent of n where the contours are approximately vertical. Thus, in the top left corner, covered by yellow shading, the more vertical contours indicate that the knockdown is dominated by w 2 , while in the green-shaded zone at the bottom right, the knock-down is controlled by n, as the contours are more horizontal. The angle of the diagonal line between these two zones is dependent on C and its location depends on D. At present a physical mechanism for this dependence has not been determined, but one theory is that the dominant parameter is the one that governs the proportion of straight fibres. When one extent parameter is large, the wrinkle affects the whole coupon in that corresponding dimension, so the other parameter dictates the proportion of fibres that are not wrinkled.
The final investigation focused on wrinkle extent in the load direction (x), which is controlled by w 1 . Fig. 17 shows plots of polynomial fitting curves based on five different w 1 values. At each w 1 value, the variances of w 2 and n are in the same range and maximum wrinkle angle is constant at 20°. When the value of w 1 is small, the curve is at a higher position, with a larger knock-down, while as w 1 becomes larger, the knock-down value reduces. This implies that small w 1 causes higher stress concentrations in the wrinkled region, which influences the failure stress negatively.

Parametric study for change of wrinkle orientation
The parametric study presented above is for the case of the loading being in the same plane as the main sine-curve wrinkle direction, but there are many manufacturing scenarios where wrinkles can run at 90°to thissee Fig. 18. Hence, this section seeks to identify the key parameters of wrinkles in this orientation that govern the compressive strength. If the wrinkle's cross-sectional profile (in the y_z plane in Fig. 18) is constant along the whole component in the load direction (the x direction in Fig. 18), all the load-bearing fibres are still straight and so there will be no knock-down in failure stress. However, this is rarely the case -a wrinkle generally has a finite extent in the load (x) direction. It has not been possible to consider all shapes governing the load-direction extent but a simple Gaussian shape has been modelled here, to determine whether there is a need for further investigation. The same process and wrinkle-definition scenarios are used as previously, but the wrinkle is rotated by 90°in the x_y plane as shown in Fig. 18. The ply out-of-plane deviation can be governed by a similar equation (Eq. (9)) to Eq. (6), but swapping the directions of w 1 and w 2 . Then maximum angle in the load direction is determined by the combined effect of A and w 2 . Here L was fixed as well as the crosssection parameters w 1 and n; A was changed from 0.2 mm to 0.6 mm and w 2 in the range 1 mm to 5 mm.
A similar modelling process was carried out and the distributions of failure-stress knock-down (%) with maximum angle as colour map and contour map are presented in Fig. 19(a) and (b). Both maps indicate that the distribution of failure-stress knock-down follows the trend of maximum angle in the load direction, which proves that this is again the main factor influencing the failure stress, even when the load is aligned transverse to the direction use for the study in Section 4.1.

Discussion
In the future, when NDT information from real samples is available and used to create FE models of as-manufactured wrinkle defects, there will be no need for assumptions in terms of shape or distribution of wrinkles. For the purposes of this study, to determine which wrinkle parameters are structurally significant, it has been necessary to bound the scope and assume flat external-surface laminates, cosine or sinewave wrinkle shapes and Gaussian profiles. The limitations of this study are listed below, followed by a comparison of the results with other research reported in the literature.
1) The modelling techniques used in this paper build on the authors' previous work [6]. Hence, the current study was limited to flat prepreg laminates. The modelling technique validated in [6] may be not suitable for other manufacturing techniques, and has not yet been validated for wrinkles that appear at the surface. Thus caution should be exercised when applying the results in this paper to components that are beyond its current scope.
2) The wrinkle topologies defined and studied in this paper are idealised for the purposes of providing trend information, to aid understanding of the dependencies of mechanical compressive strength on various wrinkle parameters. There are many more morphologies of wrinkles in industrial components that have not been studied here, such as non-symmetrical or with random shapes that are difficult to depict by a simple mathematical equation. The interdependencies of wrinkle parameters summarised from the sine and cosine shapes can however still provide useful insight on the overall sensitivities. 3) Although the experimental coupons were manufactured with artificially induced wrinkles and were represented by stepped models, the small differences in strength between these and Gaussianprofiled models shown in Fig. 9 and Fig. 10 suggests that there is little dependence of mechanical performance on this profile.  Due to the lack of either a consensus on wrinkle metrics or an approved standard to evaluate wrinkles, it is difficult to thoroughly compare the influence of wrinkle-parameters from data presented in the literature. Seon et al. [22] have similarly investigated the effect of imperfections (e.g. wrinkles and waviness) on mechanical performance using a method of automatically transferring the NDT data into structural FE models. Most their research focused on other loading scenarios, like fatigue or bending. For the static in-plane compression case, the wrinkled samples in Lemanski et al. [5], had similarities to those used for this study [6], but the wrinkles were introduced differently: two central plies were dropped to form a half-cycle wrinkled region (the exact equation to simulate the wrinkle geometry was not explained). Based on the discontinuous-plies direction (0°or 90°) and the ply quantities, there were four cases: L2, L10, T2 and T10 (L referred to 0°and T for 90°). The wrinkle distribution was kept constant in their work. The maximum angle for L2 and T2 was 8°, while this angle for the other two cases was 30°. In both their experimental and numerical modelling results, the failure load reduction was be around 55% for either L or T cases, when the angle was changed from 8°to 30°, but it is likely that the dropped plies will have an influence on the failure, alongside that of the wrinkled plies. Also for uniformly distributed wrinkles along the width, research conducted by Elhajjar and Shams [8] identified that compressive strength reduction would be~30% when the wrinkle angle changed from 20°to 30°, in samples with the wrinkles introduced using metallic rods. In Section 4.1 of this paper, where wrinkles also were uniform across the coupon width, the failure strength reduction was about 40%, when the maximum angle ranged from 4°to 28°, and this reduction would be around 10% if angle changed from 20°to 28°. The differences observed in compression strength knock-down can be accounted for by the different wrinkle-introduction methods.
As stated above, the material-property maps used to create FE models were from a simulation process and governed by Eqs. (1) and (2). In future work, direct NDT-inversion material maps will be applied to assist the model creation and control the wrinkle topology. This will allow the suggested methodology in this paper to be applied to more general wrinkle geometries. Within the stated limitations, the study presented here is sufficient for initial guidance on wrinkle metrics for flat laminates under in-plane compressive loads.

Conclusions
To assist the performance evaluation of composite components containing wrinkles, a systematic and rigorous parametric study was undertaken, to compare the dependence of compressive failure stress on various wrinkle parameters. A novel methodology, combining simulation and numerical analysis by an automated transfer process, was implemented to achieve this aim, by creating a large number of numerical models with controlled wrinkle geometries. The investigations showed a primary dependence on the maximum wrinkle angle in the load direction, which should be the main parameter to be measured nondestructively, as well as secondary dependencies on the wrinkle wavelength and wrinkle extent in all three dimensions. Quantitative results applying to the modelled case of the IM7/8552 [38] fibre/resin system with 24 plies in the sequence [45/90/ − 45/0] 3s , are summarised as follows.
1) In terms of the influence of potential wrinkle-severity parameters of amplitude, wavelength and angle in the load (x) direction, for a fixed wrinkle volume and shape, the maximum wrinkle angle is more important than amplitude and wavelength, and caused a predicted 40% knock-down of failure stress when the wrinkle angle increased from  4°to 28°when wrinkles were uniform across the coupon width, with an additional dependence on wavelength when the wavelength is small (b 8 mm in this case).
2) Changing wrinkle shape from a cosine to a sine wave, within a Gaussian envelope, has little influence on compressive failure stress. 3) For wrinkle-extent parameters: • A very localised wrinkle in the load direction, where there is less than one cycle of wrinkle, concentrates the stress and enhances the knockdown in compressive failure stress. • In the plane perpendicular to the load direction, knock-down in compressive failure stress increased by 30% as the 1 e isoamplitude surface cross-sectional area of the wrinkled region increased from 2% to 40% and when the maximum angle remained at 20°. Between the two dimensions of the wrinkle in cross-section, the smaller as a proportion of the coupon size seems to have a dominant role, possibly because it governs the number of non-wrinkled fibres remaining in the coupon. 4) When the wrinkle orientation is rotated by 90°in the x_y plane, the maximum angle to the load direction is still the main factor determining the compressive failure stress for a Gaussian profile in all directions.
Based on the study above, when wrinkles exist in composites, the recommended NDT wrinkle metrics are: 1) Maximum wrinkle angle, 2) Wrinkle wavelength, 3) Extent of wrinkled regioni.e. equivalent 1 e Gaussian half-width in all three dimensions if possible. Of these, the maximum angle is the most critical measurement.
In this work, the only wrinkle shapes considered were cosine and sine wrinkles in one direction, weighted with Gaussian profiles in both in-plane directions. The range of real manufacturing wrinkle shapes is much wider than this and it is recommended that a further study of load-direction shape should be performed, for obtaining a more accurate evaluation of component performance with wrinkles.