A VARIATIONAL MODEL FOR THE WRINKLING OF CARBON FIBRE COMPOSITES DURING CONSOLIDATION

A semi-analytical variational model for the formation of wrinkles following consolidation of a carbon fibre composites over an external radius is presented. Understanding how wrinkles form in the manufacturing of carbon fiber composites during the debulk and curing processes is important because, depending on their severity, they may compromise the structural integrity of the final part. A nonlinear elastic model is described, and energy minimization of the resulting potential energy is achieved using a gradient flow approach. The critical conditions for the appearance of wrinkles and the complex post-buckling behaviour are derived and their implications to the manufacturing processes are discussed. .


INTRODUCTION
Although the formation of wrinkles in carbon fibre composites (Fig. 1 left) and folding sedimentary rocks (Fig. 1 right) appear two disparate systems, both are examples of layered systems which are strongly influenced by geometrical constraints of layers fitting together [6].From a modelling perspective much can be gained from the cross-fertilization of ideas between each field [10].In this paper we adapt an energy based gradient flow method, previously applied to the modelling folding of sedimentary rocks, to the wrinkling of carbon fibre composites during manufacture.
Multilayered structures found in the Earth's crust may arise from a sequence of sedimentary deposits, subsequently held together by pressure from additional overlying layers.Tectonic plate movement causes in-plane compression of these layers of rocks that can fold them into convoluted patterns.The resulting deformation is not only dependent on the complex mechanical and material properties of each layer, but on the highly nonlinear geometric constraints of them fitting together.Figure 1 (right) shows an exposed escarpment of folding, where the geometric constraints have been overcome by the rocks adopting straight limbs and sharp corners, in what is known as chevron folding [15].In some cases weak layers can be obliged to fit into very tight or even singular geometries [6,7,8], and as a result various patterns, named accommodation structures, may form on a localized length scale.Traditional layers of rock have been modelled as elastic or visco-elastic beams supported by a surrounding foundation [1,15].Finding constrained energy minimizers for such multilayered systems may suffer from the existence of multiple local minimisers, due to lack of convexity of the resulting energy landscape.This is particularly the case in nonlinear buckling problems, [3,2,12,6], typical of geological folding models.Previously a constrained gradient flow method has been adopted [2,13] to overcome such problems, and this is the approach we now apply to modelling wrinkling in carbon fibre composites.
Understanding how wrinkles form in the manufacturing of carbon fiber composites during the debulk and curing processes is important because, depending on their severity, they may compromise the structural integrity of the final part. Figure 1 shows small amplitude wrinkles or folds, which have developed in the corner radius of a component.
Typically CFC parts are made up by layering a series of pre-impregnated carbon fibre layers onto a tool surface.During this lay-up process the stack of plies is consolidated at moderate temperatures and pressures to remove air trapped between layers.This debulk process aims to insure correct seating onto the tool surface and to promote adhesion between laminae.However, as a laminate consolidates over even a simple geometry the plies are forced to accommodate the imposed geometry of the tool surface.For example, consider consolidation over an external radius.As the outermost ply consolidates it is forced into a tighter geometry, Fig. 3 (right).The ply is forced to accommodate the imposed geometry by shortening.For layers which slips with relative ease, the length can be accommodated by producing so called 'book ends', Fig. 3 (right).If the resistance to slip is too high, layers may form wrinkles.
In this paper we describe a simple elastic model for the potential wrinkling of composites as a laminate consolidates over an external radius.By deriving the total potential energy for the system, the gradient flow method is adopted to find global minimisers under quasi-static loading by consolidation.Minimization of the resulting potential energy leads to representation as a nonlinear fourth order differential equation.Critical conditions for the appearance of wrinkles are derived, and the snaking post-buckling behaviour is briefly discussed.
In the work presented we see a departure from a general trend of modelling advanced manufacturing processes using large complex finite element simulations (see for example the commercial software PAMFORM/ABAQUS or various academic contributions [20,19,14]).As well as having clear visual appeal, such finite element studies have been geared towards answering a host of important questions.However, with the vast range of parameters often involved in such formulations, it is sometimes unclear how their physical counterparts might accurately be represented.While we fully appreciate the value of such approaches, our objective here is to complement them with a simplified nonlinear elastic model in which the parameters influencing the localised bucking process can be closely defined and carefully monitored.For this particular study finite element approaches struggle with the lack of uniqueness of post-buckling solutions observed, however the semi-analytical model presented is specifically designed to overcome such problems.

The model
We consider a stack of N plies each of initial uniform thickness h which have been laid over an external corner/tool surface, which is characterised by radius of R t being swept through 1  2 π rads, see Fig. 3.A uniform debulk pressure consolidates the laminate by a percentage γ.Each layer, numbered from the outside layer (layer 1) inwards, is described by the radius of curvature, and the small wrinkle deformations u i (x i ) away from this radius .Here x i characterises the length of arc swept out by the radius of i th layer (R i ).At this point it is useful to the note the change of variable dx i = R i R t dx t , where x t is the arc-length parameter around the fixed tool face.
The number of plies in a typical component can be large, see for example Fig. 1 where N = 40; by considering each layer individually this would lead to a system of N coupled fourth-order ordinary differential equations.To simplify the analysis at this stage we make the assumption the wrinkle deformation of the layers decays linearly to zero at the tool face the surface.This assumption can justified by viewing micrographs of wrinkles in corner sections, like those seen in Fig. 1.Consequently, the complete deformation of the laminate can be written in terms of the outer layer, described by the displacement u 1 , which is abbreviated to u, We now consider a quasi-static model of the wrinkling process, where the controlled parameter is the percentage consolidate γ.The consolidation of the laminate is modelled as a conservative system and we seek to find energy minimizing solutions for that percentage consolidation.Each ply is modelled as an inextensible linear elastic beam of length As the layer consolidates, it is forced into a small arc length, and accommodates this extra length by either slip over other plies or by forming a wrinkle.The slip between plies is resisted by two linear elastic spring of stiffness k s .Finally, we consider each layer to lie on a nonlinear Winkler type foundation, this offers resistance to wrinkle deformations.
The setup and the parameters of the model are summarized in Fig. 4. The total potential energy function for the system is now derived, which is given by the sum of the strain energy in bending, work done against the axial springs while slipping and work done into the foundation.

Bending Energy
The bending energy over a small segment of the i th layer dx i is given by where κ(x i ) and x i are the curvature and arc length of the i th layer respectively.Integrating over the domain of each x i , the total bending energy is the sum of the bending energy in each layer, therefore Since the wavelength of the buckling deformations are small, the changes in curvature of the mid-surface of each layer can be approximated by the curvature of the i th layer κ i ≈ Here primes denote differentiation with respect to x i . where (5)

Work done in slip
We now calculate the total potential energy stored in the axial springs.The strain the springs can be determined by the total shortening imposed by consolidation, given by minus the shortening achieved by the wrinkle deformation Since the springs are assumed linear elastic, the strain energy in the springs is We believe this a conservative modelling approach since the true behaviour is more closely modelled by a visco-elastic law.In this case k s represents the initial maximal resistance to slip, any strain relaxation due to visco-elastic behaviour will reduce this stiffness, and hence reduce the likelihood of wrinkle formation.This simplified elastic model provides an upper bound the true physical counterpart.

Lateral resistance to wrinkle deformations
We assume that wrinkle deformations in each layer are resited laterally by the force (f i = f (u i )) due to lateral displacements being strictly local and normal [16].We consider the nonlinear stress-strain relationship This choice of stress-strain relationship allows for a response that is initially symmetric and linear elastic, followed by an asymmetric elastic softening response away from the corner.The geometry of large deflections introduces nonlinearities, these have a de-stiffening effect [12], encouraging layers to localise and form tight Ω-shaped folds.However, for large deflections locking-up is mimicked by the stiffening symmetric cubic term, either by large deflections into the laminate or outwards against the debulk pressure.
The strain energy stored in the foundation is then given by where

Total potential energy of the system
The total potential energy of the system is then the sum of each of the contributing terms given by Here we also denote the Fréchet differential of F by

Energy minimization using a gradient flow approach
The gradient flow method solves the partial differential equation Starting with random data, the long term behaviour of this system converges to stationary solutions of E which are assumed to be constrained local minima.By repeating this process for a 'large' number of different initial positions the global minima is selected.
We now solve the gradient follow problem using a standard finite element approximation, over the domain x t ∈ [− , ], coupled with simply supported boundary conditions.By defining the perturbed function as u ε = u + εϕ for ε ∈ R and a suitable test function ϕ, the weak form of ( 8) is given by This is equivalent to Since this weak formulation has differential of order two and no higher we approximate u by piecewise cubic functions upon a uniform mesh ξ = 2 /M , where M ∈ N is the number of node points.Therefore u is approximated by where φ k and ψ k form a basis of piecewise cubic functions on ξ with simply supported boundary conditions [18].By defining the vector the finite element equivalent of (2.5) is the first order system in U , Here A, B and C are sparse 2M × 2M matrices given by and with similar values to A ij for other ranges of i and j.Here primes denote differentiation with respect to x.We note that most elements of A, B and C are zero.The components D i of D, to be found by a suitable quadrature rule, are given by and The first order system (11) can then be solved using a typical first order solver in Matlab, for example ode45 [17].

Solutions profiles, wrinkling bifurcation and snaking behaviour
The individual solutions profiles, for example Fig. 5, and the overall system behaviour will be sensitive to choice of modelling parameters.Here we discuss the general behaviour of the system, and derive the relationship between these parameters for which the system has wrinkling type solutions; we then discuss the post-buckling behaviour of the system for a set fixed parameters.
Figure 6 shows a plot of percentage consolidation γ against load in the axial springs P .Starting from zero consolidation the solution path follows the fundamental path (u ≡ 0) shown by the red dot-dash line.For this case the outer layer accommodates its extra length by slipping over the other layers a total distance λ, and the axial springs exert a reactive compressive load k s λ.Solutions follow this path until the load in the axial springs reaches the critical buckling load of the system, and bifurcates into an unstable post-buckling path, shown by the blue lines.
The location of the initial bifurcation point can be derived from the Euler-Lagrange equation associated with minimizing solutions of E close the fundamental path (u ≡ 0), given by where we recall that By seeking solutions of the form u = exp(ξx) the characteristic equation of ( 15) is given by At the parameter values the eigenvalues ζ coalesce to form two imaginary pairs, and correspond to the linear buckling load of a strut.Figure 7 shows a plot of k c s against γ and demonstrates that for moderate consolidation the resistance to slip must be low if wrinkles are not to form.How low would require careful parameterisation of such a model, which is discussed briefly in the following section.
We now consider the more complicated post-buckling behaviour, which is more clearly indicated by plotting L 2 (u ) a measure of the total size of the wrinkle deformation, against the load in the axial springs P , Fig. 8.
The plot is made up of a collection of discontinuous curves, where at various points the global minimiser jumps between different localised solution profiles.The solid blue lines in Fig. 8 represent solutions with an odd number of localised humps, whereas the dash line  gives the even counterparts.Physically, the initial de-stiffening of the foundation results in a subcritical Hamiltonian-Hopf bifurcation into a single hump homoclinic package, away from an initially periodic solutions.The subsequent re-stiffening effect brought about by the symmetric cubic term of f restricts the amplification of the central hump.Snaking behaviour ensues, where the global minimisers winds its way either side of Maxwell load P M [12,5,4].With sufficient external disturbance the Maxwell load provides a limiting load capacity under dead loading, indicating that this system is particularly sensitive to imperfections.

Concluding remarks
The work presented in this paper shows an adaptation of models previously used to successful understand some of the aspects of folding rocks, to the a study on the wrinkling of carbon fibre composites.Although a simple model is presented it demonstrates complex non-linear behaviour, which encapsulates some of the key aspects of the localised wrinkles observed in practice, for example Fig 1 .Such folding is fundamentally non-linear, since the wrinkle deformations are sufficiently large.In this case closed-formed solutions to such nonlinear problems cannot be expected, and robust numerical methods are the only means to attain such convoluted solutions.The gradient flow method outlined in this paper is particularly powerful, as it offers a stable numerical scheme by which the equilibrium equations are solved alongside calculating the potential energy of solutions.This is vital when evaluating on such an undulating energy landscape.
Although the current model is simple, and future extensions are envisaged, it provides some suggestions to manufacturing on the important parameters to produce wrinkle free parts.This is particularly highlighted in the critical condition for the formation of wrinkles (17).A few observations are: Lower k s : Wrinkles do not form if k s is sufficiently small.Strategies to encourage interlayer slip will be important.These may include providing optimum heat to the laminate during consolidation to decrease the viscosity of the resin.

Decrease γ:
The critical value of k s is inversely proportional to γ. Therefore the chance of forming wrinkles can be reduced by consolidating the laminate at regular intervals, and therefore reducing γ at each debulk step.
Increase R t : The chance of forming wrinkles decreases as the R t increases.By rearranging equation ( 17) a design rule could be obtained for the smallest allowable value of R t .
How regular debulk actually needs to be, the exact heating temperature during consolidation or the limiting value of R t all require careful parameterisation of the model.The next step of this study is to carry out a series of experiments, to obtain the required modelling values.This will particularly focus of the understanding the behaviour of inter-layer slip, parameterised by the elastic stiffness k s in this model.
The gradient flow approach here is used to model a quasi-static system.However, this method could easily be adapted to model dissipative systems; and the work could be readily extended to consider visco-elastic material behaviour, which may be necessary if good experimental comparisons are to be achieved.However the approach here does have certain modelling advantages over more complex rheologies, Modelling has sometimes failed in the past not because of the formulation, nor the resulting governing equations, but because of poor choice of solution [6,12].Extra criteria for selection are sometimes neither necessary or useful.A vital extension to this work will be to consider the consolidation of 3D parts.The directional nature of different plies through out a laminate means that in directions orthogonal to the fibre directions these plies may offer no or little bending stiffness, yet may be heavily constrained by how the consolidation takes place in the perpendicular direction.These interactions while consolidate over a 3D geometry are complex, and we believe they are the key factors in the wrinkling of carbon fibre composites during manufacture.

Figure 1 .
Figure 1.Left: A micrograph of small amplitude wrinkles in a typical corner radius.Right: Chevron folds developed in thinly bedded shales and sandstone at Millook Haven, Cornwall, UK.The geometric constraints of layers fitting together have force layers to adopt straight limbs punctures by sharp corners.

Figure 2 .
Figure 2. (Left) shows the geometry of the tool face, with N layers of initial thickness h laid on top.(Right) Demonstrates the book end effect created as the the laminate consolidates without wrinkles.Here λ is the amount of slip in the outer layer which consolidate from an initial radius of curvature R 1 (0) to R 1 (γ), during a percentage consolidation γ.

Figure 3 .
Figure 3. Setup of the model, simplified to a single layer of adjusted stiffness B, laid around the corner of radius R(γ), and on a nonlinear Winkler foundation, describe by the local normal force f (u).Linear elastic axial springs of stiffness k s resist slip as the laminate consolidates by a uniform percentage γ.

Figure 4 .
Figure 4.An example of a typical wrinkle deformation found by solving the gradient flow method described in this paper.The parameters for this solutions sample is given by B = 1, k f = 1, K f = 0.3 h = 0.8, M = 20 and R t = 40.

Figure 5 .
Figure 5. Plot of percentage consolidation γ against the load in the axial springs P , for the following parameter values: B = 1, k f = 1, K f = 0.3 h = 0.8, M = 20 and R t = 40.The red (dash-dot) line shows the initial fundamental path, and the blue lines (solid and dash) shows a snaking post-buckling path of wrinkling solutions.

Figure 7 .
Figure 7. Plots of L 2 (u ) = t 0 u 2 dx, a measure of the size of the wrinkle against load, for the parameter value B = 1,k f = 1,K f = 0.3 h = 0.8, M = 20 and R t = 40.The plot is made up of a collection of discontinuous curves, where at various points the global minimiser jumps between different localised solutions profiles.The solid line shows solutions with an odd number of localised humps, where the broken line shows those with an even number.