A novel criterion for the prediction of meso-scale defects in textile preforming

In numerical simulations of reinforcement preforming, fabrics are often modelled using continuum approaches with homogenised material properties in order to reduce the computational costs. To predict the occurrence of meso-scale defects, in particular when membrane ﬁ nite elements are used for modelling, a novel wrinkling criterion is proposed here. This criterion relates the onset of reinforcement buckling to material and process parameters. It is demonstrated that the criterion correctly predicts the occurrence of meso-scale defects and is more accurate than criteria based on a constant shear locking angle or critical strain for the defect onset, which do not to take into account local processing conditions such as friction or normal pressure.


Introduction
Textile reinforcements are widely used in the manufacture of composite components because of their superior drapability compared to multi-layer preforms from unidirectional reinforcements. Non-crimp fabrics (NCF) are successfully used in industrial mid-volume production, e.g. of lightweight vehicle components. Component production requires good quality of preforms with as few defects as possible. Wrinkles in the preform are of particular concern since they can result in fibre misalignment in the finished component, which may lead to reduced strength [1]. Numerical simulations can be used to assess the process design, and evaluation of simulation results can ultimately help to avoid or reduce defects.
Several approaches to textile forming simulation have been developed over the years, offering different levels of fidelity as outlined in two reviews [2,3]. Kinematic models represent a textile reinforcement as a pin-jointed net with inextensible links, but ignore the material properties of the textile. Based purely on geometrical considerations, these models can predict fibre orientations with good accuracy for scenarios with low to medium geometric complexity. One of the main advantages of kinematic models is their computational efficiency, with running times in the order of seconds. At the other end of the spectrum of modelling approaches, forming of textile reinforcements is simulated by representing each deformable yarn separately. The yarns are discretised, typically employing finite element (FE) methods and solid elements, taking into account the mechanical properties of individual yarns. Yarn interlacing and contact interactions between yarns are modelled explicitly. However, this level of fidelity comes at the expense of a significant increase in computation time. A middle ground can be found in the continuum approach, which represents a textile reinforcement as a homogenised medium modelled employing FE methods. Here, volume (solid) elements are used for thick reinforcements, and shell or membrane elements for thin reinforcements. This approach lies in between the kinematic and discrete approaches in terms of both fidelity and computational cost; computation time is in the order of minutes or hours, depending on the complexity of the model.
In the continuum approach, the homogenised fabric behaviour is described by effective moduli in the fibre directions and a shear resistance; the behaviour is frequently characterised using picture frame shear tests [4] and modelled using a non-orthogonal constitutive model [5]. This type of constitutive model traces reorientation of the yarns during a fabric deformation and describes the behaviour of the fabrics in a sheared configuration (i.e. with the yarns non-orthogonal). Models of thin textiles based on membrane elements (or their variants) were shown to be capable of predicting the local shear angles of reinforcements in matched tool forming for cases where no macro-scale wrinkles were expected to occur [6][7][8][9][10][11]. However, due to the lack of bending stiffness, it was shown that membrane element models are not capable of realistic prediction of large wrinkles [12][13][14]. Shell element models, which include additional degrees of freedom to account for the bending stiffness, are more computationally expensive and require additional experimental data for calibration. However, they are capable of capturing realistic shapes of large wrinkles which can arise when a reinforcement is formed with little or no constraint as was shown in [15][16][17].
Whilst accurate enough to predict formed shapes and fibre orientations, the continuum approach cannot reproduce meso-scale defects such as wrinkling of individual yarns since it is based on homogenisation of meso-scale features. It is possible to employ a mixed modelling approach by adding beam elements to shells or membranes to account for the compressive behaviour of the yarns [18][19][20], but this increases the computational complexity. Alternatively, a wrinkling criterion, based on locking angle or critical in-plane strain, can be employed in the post-processing stage of simulation results. A frequently used wrinkling criterion is based on the assumption that reaching a critical shear angle results in shear locking and fabric wrinkling, although this was shown to be inaccurate in cases where a fabric experiences not only shear but a combination of shear, in-plane and bending deformations [15]. The main limitation of these approaches is that wrinkling threshold values in terms of angles or strains are postulated to be constant for all local conditions such as normal pressure or friction despite evidence suggesting otherwise.
This paper aims to explore the ability of a membrane element model to predict meso-scale wrinkling using both existing and novel wrinkling criteria. The proposed novel criterion, based on buckling theory, relates local textile deformations to the given material and process parameters in order to predict the onset of wrinkling. The wrinkling criterion is applied to the results of forming simulations, and the predictions show good agreement with experimental observations of defect locations in the preform.

Comparison of shell and membrane elements
The difference in computational time between models using shell and membrane elements is studied by simulating the forming of plates using an isotropic elastoplastic material model (Young's modulus E = 210 GPa, Poisson's ratio ν = 0.3, yield stress σ y = 110 MPa). The geometry for forming using a hemispherical punch is shown in Fig. 1. All parts of the tool, i.e. punch, blank holder and die, are modelled as rigid bodies. The force applied to the blank holder is 7000 N, and the coefficient of friction between the blank and any part of the tool is μ = 0.2. The punch velocity is set to 50 mm/min which ensures the absence of dynamic effects in the simulations. The Abaqus/Explicit solver was used for the simulations. Mesh convergence studies showed that models of at least 150 × 150 elements are required for the shell and membrane models to achieve converged solutions. The computational time for the membrane model consisting of 150 × 150 elements was 30 min, less than half the time for the shell model, 66 min, on a standard desktop computer. If required, further reduction of the computational time is possible employing mass and/or time scaling [21,22]. The final shapes predicted by both models were sufficiently similar to be considered identical for the purpose of forming simulations. However, it should be noted that this might not be the case for more complex geometries.

Textile forming model
The numerical model used for the simulation of forming of a noncrimp fabric (NCF) employs a non-orthogonal constitutive model which is pre-defined in Abaqus/Explicit under the keyword *FABRIC and can be used to trace yarn orientations during the deformation. Input data for the shear behaviour were obtained from picture frame shear tests and friction tests. In-plane moduli were calculated using manufacturer's data for carbon fibre.

In-plane shear behaviour
A biaxial carbon fibre NCF with pillar stitch pattern (Hexcel FCIM359 [ ± 45°] with stitches running at 0°) was characterised by Chen et al. [23] using picture frame shear tests. In these tests, the shear resistance as a function of the shear angle was derived from data for the force applied diagonally to corners of a square picture frame as a function of the displacement, which induces shear in the fabric specimen clamped in the frame [4]. Experimental data for the recorded force as a function of shear angle are plotted in Fig. 2. The asymmetry of the curves and the non-monotonic increase in shear resistance with increasing positive shear angle is related to the presence of stitches in the textile. In positive shear, the stitches are in tension and hence resist the shearing of the textile. At a threshold angle (here, approximately 0.5 rad), the stitches start to fail and the force reduces rapidly. For shear angles greater than approximately 0.8 rad, the shear resistance is dominated by lateral compression of the yarns in the fabric. For negative displacement, the stitches buckle easily in compression and do not contribute to the shear resistance of the textile.

Friction
Measurements of tool-fabric and fabric-fabric friction coefficients were performed according to ASTM D1894. For relative movement at constant velocity, friction coefficients for the surface pairings were calculated from the ratio of the applied tangential force to the normal force. The friction coefficients, which were assumed to be isotropic, were 0.23 and 0.36, respectively, for fabric in contact with a flat tool surface and with another fabric layer when the layers have the same orientation. The actual friction between two layers of fabric is anisotropic but it was assumed that the reorientation of two layers in this problem is not large enough for this effect to have a significant influence on results.

Elementary model verification
Non-orthogonal constitutive models are thought to be the most appropriate models for composite forming, and various implementations have been described in the literature [5,24]. Commercial implementations, such as the *FABRIC model in Abaqus/Explicit, are less versatile than some implementations described in the literature, but can be used without the need for any additional coding. As input for the model, it is suggested [19] to use the shear stress, τ, derived from picture frame shear tests, as a function of the shear angle, α, given by Here, ψ 12 is the angle between yarn directions, where F is the force applied diagonally to the picture frame, L 0 is the edge length of the picture frame, v 0 is the initial volume of the material, and d is the cross-head displacement as shown in Fig. 3. The in-plane tensile behaviour of the bi-axial fabric was assumed to be linearly elastic and the Young's modulus in both directions was computed using the rule of mixtures assuming the Young's modulus of carbon fibre to be 276 GPa and a fibre volume fraction of 0.5 in each of the directions [23].
A single finite element, where kinematic boundary conditions were applied to simulate a picture frame shear test, was used to verify the material model input. The results derived from the simulations accurately reproduced the input data, verifying that the material model is capable of dealing with the non-linear asymmetric shear curves.

Hemisphere forming model
A model which replicates forming of a reinforcement with a hemispherical punch ( Fig. 1(a)) was created in Abaqus/Explicit using the dimensions shown in Fig. 1(b). Two square plies of the NCF with dimensions 300 mm × 300 mm were held by a blank holder applying a normal force of 1200 N. The NCF plies were placed such that the fibres were oriented along the blank edges. The blank was modelled using membrane elements M3D4 with a size of 2 mm × 2 mm. The explicit solver was run using an automatic time step selection.
Shear angles predicted by the forming simulations were in good agreement with results from corresponding experimental measurements obtained using a coordinate measuring machine and grid strain analysis. They were also in good agreement with simulation results obtained by Chen et al. [23] using the Abaqus/Explicit VFABRIC subroutine, which provides a framework for user-defined non-orthogonal constitutive models. Maps of local shear angles from the present simulations are compared with results presented by Chen et al. [23] in Fig. 4. The simulations capture the asymmetry of the final shape of the blank as well as local shear angle distributions. The results show that the pre-defined material model for textile deformation, implemented in Abaqus/Explicit and called by the *FABRIC keyword, can be reliably used instead of user-defined subroutines for the materials models, which require additional effort for implementation and validation.

Wrinkling onset
In composites forming, wrinkling is often assessed using a shear locking [25] or compressive strain concept. For example, Chen et al. [23] correlated experimental and numerical results and equated the strain at wrinkling onset to a critical value of 0.03. However, both criteria assume that the critical value is universal, independent of processing conditions, material and contact behaviour. Moreover, there is rarely a physical way to determine the critical value. For sheet metal forming, wrinkling during forming has been assessed by energy-based criteria [26,27]. A similar approach is adopted here to predict the wrinkling onset.
It is assumed that wrinkling occurs at the meso-scale level, implying that individual yarns should be considered. Instead of trying to predict a final shape of a wrinkle, which would depend on how exactly a yarn would interact with other yarns and/or the tool, this model aims to predict wrinkling onset only. The wrinkling onset of individual yarns can be related to in-plane or out-of-plane buckling modes. These two modes require different boundary conditions to account for the structure of a textile.
The out-of-plane buckling of a yarn can be modelled as a buckling of a beam of infinite length, which is pressed against a substrate by a uniformly distributed load, q, as shown in Fig. 5. This distributed load represents pressure from tool/fabric and fabric/fabric interaction as well as a pressure generated by tension induced by blank holder. In the presence of the distributed load, an external axial compressive force, P 0 , is applied to the beam as a result of (local) fabric deformation. A friction force opposes axial compression, which depends on the friction coefficient between the beam and the substrate, μ, and the distributed load, q, and is assumed to be independent of the sliding velocity. When  the compressive force is increased above a critical value, which depends on the beam stiffness and the friction behaviour, the beam buckles with a buckling length L, and the length of slippage of the beam on the substrate is equal to L s . The in-plane buckling of a yarn can also be modelled as buckling of a beam but the main difference to the out-of-plane buckling is the presence of neighbouring yarns which constrain lateral movement. This constraint can be modelled as an elastic foundation with high stiffness (to represent neighbouring yarns which can deform but do not move) in addition to the effects already included in the model for the out-ofplane buckling. Alternatively, the in-plane buckling can be modelled as immediate buckling of several neighbouring yarns in a way that they would not constrain each other. Addition of this extra constraint will increase the buckling force required to generate an out-of-plane buckling and, therefore, it is expected that the out-of-plane buckling will be initiated first.
The differential equation for the buckled part of the beam in out-ofplane mode is [28]: / , E is Young's modulus of the beam, I is the second moment of area of the beam crosssection, and P is the force in the buckled part of the beam. From the boundary conditions of zero slope at the ends of the buckles, it can be shown that = nL nL tan( /2) /2 and hence that the buckling force P for the lowest root is given by: The external axial compressive force away from the buckle, P 0 , is then a combination of the force in the buckle, P, the friction force in the sliding part of the beam, and the friction forces at the ends of the buckle: The slippage length, L s , is then found using Eq. (6) and the equation for the change in length: 27 is the shortening of the beam [28]. It can be seen that the first term in Eq. (6) is equal to the classical Euler buckling solution, which predicts a monotonic decrease in buckling load with increasing beam length, and the third term is the contribution of the frictional forces, which results in a non-monotonic behaviour, as shown in Fig. 6. An increase in the friction coefficient results in an increase in the buckling force for a particular buckling length. However, the most important difference is that in the presence of friction the buckling load has a minimum below which no buckling is possible for any length.
In contrast to the wrinkling criteria described earlier, which are based on a somewhat arbitrary single critical value, the criterion based on Eq. (6) can be linked to the material parameters and local pressure.  The value of Young's modulus, E, used in Eq. (6) was set equal to that of the reinforcement material. The local contact load predicted by the FE model was used as the normal load q for every finite element in the model. Other parameters can be found in Table 1. The upper and lower bounds on the bending rigidity, EI, in Eq. (6) can be estimated using the second moment of area of a solid beam of a rectangular cross-section, bh /12 3 , and that of N f individual fibres, × N πr /2 f f 4 , which assumes that all fibres lie on the neutral axis. A more sophisticated estimation of the second moment of inertia can be based on assuming a square fibre arrangement within a rectangle with sides equal to b and h and applying the parallel axis theorem but it would result only in a moderate increase of the second moment of inertia compared to the previous assumption. The difference between these two values is up to several orders of magnitude but it is likely that it tends to be closer to the upper value due to the inter-filament friction and other factors. It was shown by Wang [29] that using the upper bound value leads to overestimation of the buckling load. The same work used a scaling factor of 0.041 applied to the second moment of area of a solid rectangular cross-section to fit the experimental data for buckling of thin strips of textile reinforcements. In a general case, the scaling factor depends on many parameters such as fibre diameter and inter-filament friction. In the absence of experimental data on bending or buckling of tows, the scaling factor was selected equal to that in the work of Wang [29]. A sensitivity analysis of Eq. (6) showed that the predicted buckling force increases by 20% when the scaling factor is doubled.
Once buckling is initiated, the shape of the beam (yarn) changes. The exact buckling amplitude is not known since the buckled yarn is constrained by e.g. the tooling surface and the stitches in the fabric. The former obviously directs an out-of-plane wrinkle into an in-plane waviness despite the initial buckling mechanism. The latter imposes additional conditions on the buckling problem (6) as the buckling length cannot be greater than the distance between the stitches. Furthermore, the presence of stitches can result in the occurrence of a higher buckling modes. However, higher buckling modes would require much higher buckling force. For example, for a fixed length, the Euler buckling force in Eq. (5) would be about 3 times higher for the 2nd buckling mode, almost 6 times higher for the 3rd mode and 10 times for the 4th mode. The transition between out-of-plane and in-plane buckling as well as appearance of higher buckling modes will be discussed below.
In post-processing of the results, the critical axial force and strain are evaluated locally for every element by finding the minimum buckling force and a corresponding buckling length, L. This is performed by evaluating buckling forces for a range of buckling lengths and then finding the minimum. The critical strain is then compared to the local strain, and a wrinkling map is created identifying wrinkle formation at the elements in which the local strain exceeds the critical strain.
Predictions for wrinkle formation according to three criteria, the locking angle criterion, the maximum strain criterion, and the new criterion based on friction-modified Euler buckling proposed here, are shown in Fig. 7. In the locking angle criterion, wrinkling was assumed to occur at angles greater than 40°for both shear directions [30]. In the maximum strain criterion, the threshold compressive strain was set to 0.03 [23]. Fig. 7 shows that the locking angle criterion predicts the smallest ply area affected by wrinkles with the selected value of the critical angle. The wrinkling maps predicted with the maximum strain criterion and the new criterion proposed here are similar. They agree well with experimental observations as highlighted in Fig. 8. The flat parts of the preform, where the blank is compressed between the blank holder and the die, exhibit in-plane wrinkling which is predicted correctly by both criteria. The radius between the flanges and the hemisphere and the areas at the sides of the hemisphere on a diagonal between the blank corners exhibit out-of-plane wrinkling of yarns. The areas with wrinkles predicted by the maximum strain criterion have a very sharp boundary, while the areas predicted by the new criterion proposed here are fuzzy and visually smaller. In addition, the proposed criterion also predicted localised wrinkling closer to the pole of the hemisphere. The difference between these two predictions arises from the dependence of the proposed novel criterion on friction and local contact pressure. Friction can delay the onset of wrinkling, and correlates to the reduction in size of the wrinkling areas at the sides of the hemisphere.
The criterion proposed here assumed that the out-of-plane buckling mode is the first to initiate and did not include any consideration of inplane buckling mechanism. However, it can be noted that the areas of predicted wrinkling onset correlate well with areas which have both inplane and out-of-plane wrinkling. Comparison of the predicted buckling length in these areas between each other and the distance between rows of stitches, which is equal to 4.2 mm, gives additional information for future improvement of the model. The regions which correspond to outof-plane buckling had buckling lengths between 2.6 mm and 5 mm, while regions corresponding to in-plane buckling had a buckling length between 10 mm and 25 mm. It was observed that in most of the regions the in-plane waviness spans between more than one rows of stitches in a continuous manner. It is currently unclear if the observed in-plane waviness was initiated as out-of-plane buckling which then transitioned into in-plane waviness or through some other mechanism. However, the presented observations can be used to improve the presented model further.
While the improvement over the maximum strain criterion is not significant in terms of predictive accuracy, the proposed criterion eliminates the need for calibration with a somewhat abstract critical strain value and relates the onset of wrinkling to the known material and process parameters. The sensitivity of the predicted location and amount of wrinkling to the selected critical value for the locking angle and maximum strain is shown in Fig. 9. The dependence of the wrinkling predictions on the chosen critical values make it more difficult to validate and use these criteria.
Further improvements to the proposed criterion are possible, such as the introduction of a 2D buckling model. This could account for the Fig. 6. Effect of substrate friction on critical force for buckling of a beam (qualitative).   biaxial stress-strain state as well as frictional slippage, and has the potential to improve the accuracy of the predictions. Moreover, the present work shows that the onset of wrinkling can occur at a range of compressive strain values rather than at a single value. Finally, the presented criterion can be used to identify the areas where a suitable friction coefficient can be used to delay and control the wrinkling onset.
In other works, tool-tow and tool-fabric friction were clearly shown to be dependent on surface roughness [31,32] which can be adjusted using various manufacturing methods. For example, increasing the local friction coefficient in areas with high compressive strain may lead to a postponed or even entirely prevented wrinkling onset. At the same time, a controlled decrease of friction in areas with low compressive strain may reduce the required force applied to a forming punch without leading to the formation of wrinkles. Changes in friction can also lead to the change of the overall draping behaviour.

Conclusions
Forming of a NCF reinforcement was simulated employing a FE method with membrane elements and a non-orthogonal constitutive material model. Using a commercial solver, Abaqus/Explicit, and a standard material model reduced the required coding and validation associated with a custom material model. Abaqus/Explicit also provided flexible control over the solver options and contact formulation. The use of membrane elements reduced the computational time by a factor of two when compared with the same model using shell elements. Input data for the material model were obtained from picture frame shear tests and friction measurements. This ensured good correlation between simulations of hemisphere forming and the corresponding experimental observations. Formation of meso-scale defects, in particular wrinkling of yarns, was first assessed using two existing criteria based on locking shear angle and compressive strain. While the first criterion underestimated the size of ply areas affected by the appearance of wrinkles, the second criterion required tuning for particular cases as the material and processing conditions can vary over the component. A novel criterion based on modified Euler buckling was developed to encompass several process parameters such as the coefficient of friction and the local normal pressure. An additional term in the critical buckling load, related to frictional slippage, increased the buckling resistance of the yarns. The novel criterion also provides a relationship between the material parameters and the formation of wrinkles and which can be used to study sensitivity of wrinkling to these parameters.
A case study of hemisphere forming showed that the presented methodology of coupling macro-scale continuum forming simulations with a criterion to detect meso-scale defects allows experimental observations to be predicted with good accuracy and at low computational cost. Fig. 9. Predicted onset of wrinkling using locking angle of (a) 35°, (b) 40°and (c) 45°, and using the maximum strain of (d) −0.025, (e) −0.030 and (f) −0.035. Red colour represents areas affected by wrinkles; blue colour represents areas not affected by wrinkles.