Preforming of non-crimp fabrics with distributed magnetic clamping and Bayesian optimisation

A novel preforming process was developed for non-crimp fabric (NCF) materials that generated in-plane tension through discontinuous blank boundary conditions. The method employed magnetic clamps and was designed to be both flexible and scalable, with clear routes to industrialisation. The capability of the process was explored in physical trials for a hemispherical and a cubic geometry. Characterisation of a biaxial veiled NCF showed the veil had a dominant effect on the bending mechanics. Subsequently a macroscale finite element model was developed to include an efficient bending idealisation and non-orthogonal in-plane material behaviour. Finally, global process optimisation of the preforming process was demonstrated. The optimisation approach used Gaussian process modelling with a periodic kernel to estimate the wrinkle size for untested clamping arrangements and then deployed Bayesian optimisation to find the optimal configuration. Results indicated that distributed magnetic clamping was effective and amenable to surrogate modelling.


Introduction
Dry fabric composites offer accelerated cycle times through efficient processing, out-of-autoclave manufacture, longer shelf life, and reduced dwell times. In comparison with prepregs, fabrics are also much easier to work with and manage on the shop floor. Non-crimp fabric (NCF) composites are now widely adopted across the aerospace and automotive sectors due to their reduced mass, as a result of keeping the fibres inplane, when compared with woven materials. Manufacture and production processes, however, remain inflexible, 1,2 resulting in added barriers in developing open, multifunctional systems toward the goals of industry 4.0. Even within highly engineered settings, such as the aerospace sector, composite manufacture is predominantly guided by a practical understanding of forming processes. Academic research into NCF forming is often outpaced by pressures on industry. Leading manufacturers are continually developing new protocols and manufacturing techniques, frequently with proprietary intellectual property rights. Research has also suffered from a lack of standardisation in what is a relatively new field. This adversely affects knowledge transfer and minimises the development of novel predictive modelling tools that can leverage diverse data sets. Herein, a unique adaptable manufacturing process is demonstrated, which allows tailored blank holder conditions using minimal hardware. The process employs magnetic clamps to provide scalable and flexible control of local mechanics through material tensioning. Bayesian optimisation is subsequently employed in combination with a predictive numerical model to effectively control the process.
Thus far, forming processes for NCFs have largely been borrowed from their crimped woven counterparts. These processes, including techniques such as diaphragm forming (DF), are beset with variability with little opportunity to accurately control boundary conditions (BCs) during preforming. They are also difficult to model, requiring an understanding of dynamic vacuum rates and fabric pinning, exacerbating variation in the quality of resultant parts. 3 Increasing demand for parts of additional geometrical complexity indicates a need for bespoke forming operations. External to this, material development is also ongoing, meaning additional lead time is required to fully characterise any given material. NCFs are highly engineered, and hence many variations of this material category exist. Within NCFs, biaxial fabrics are most common (Figure 1). They contain uniaxial fibre plies in two directions which are usually initially orthotropic. Additions such as toughening veils 4 are now common in the aerospace sector as well as a plethora of other modifiable parameters via the stitch. 5 Asymmetric shear of some NCFs 6-8 adds complexity to the material mechanics relation to geometry and has been shown to cause significant wrinkling. 6,9 This asymmetry is the result of resistance to shear when one of the bias directions is aligned with a stitch direction. Asymmetric shear, however, can be reduced through choosing a fabric with balanced in-plane shear properties, such as biaxial NCFs which do not have a stitch path in the bias direction. A second method is through boundary control, which is a convenient and impactful process control parameter. Simple demonstrator geometries are typically used in forming studies to isolate desired effects. For example, a hemispherical tool, due to its implicit axisymmetry, facilitates evaluation of intraply shear. This is an active area of research, with the principle of differential shearing explored through spring-clamping, 10,11 overstitching, 12 resin printing and pressure pads. [13][14][15][16] These methods provide feasible routes to process control, but still carry barriers that prevent adoption in industry. These include high tensile forces and damage to the material, whilst requiring high accuracy in placement. Similarly, many methods require large mechanical assemblies and secondary processes, all of which increase production times. Continuous blank holders, which prevent out-of-plane bending can be advantageous in reducing unwanted out-of-plane wrinkles, however, continuous blank holders can also prevent beneficial out-of-plane deformation, particularly in relation to positive curvature. By damping this bending mechanism, material can remain inside the part and increase in-plane wrinkle severity. Discontinuous variable clamping, that allows more out-ofplane deformation in discretised locations, combines several of the benefits of the aforementioned techniques and few of their challenges. Specifically, the use of distributed magnetic clamping is explored in this paper. Distributed magnets allow for variable clamping forces at discrete locations, and can be applied to: the periphery or interior of the NCF sheet, against a supporting plate or against the tool itself. Magnetic consolidation has shown promise in reducing void content in vacuum assisted resin transfer moulding and demonstrates that permanent magnets could be easily utilised in production. 17 Large mechanical assemblies can be avoided and, thanks to the development of increasingly cheap suitable alloys, magnetic clamps are available to all levels of part production. The ease of application of magnetic clamps makes them appealing for studying boundary control, localised deformation, rate control, differential boundary forces, and multi-step protocols.
Simulation is a key tool towards understanding material behaviour. The commercial finite element analysis (FEA) package Abaqus is frequently used, making use of both membrane models 10,12 and shell models. 18 Membrane models have the advantage of computational efficiency, but at the detriment of mechanical validity, particularly in forming protocols that include bending, whereas typical shell models overestimate the compressive stiffness. Novel developments to these models include the use of submodelling, 19 coupled elements, 20 laminated shell elements 21,22 and more recently, coupled membrane and bending idealised shell models. 3,[23][24][25] Despite these advances, however, accurate computational models require substantial runtimes, meaning most models cannot be directly applied towards the global optimisation of industrial processes. Similar considerations apply to physical testing. For instance, there is a large reliance in research on digital image correlation (DIC). However, scaling this into industrial settings is often impractical. Stereo vision systems require line-of-sight access to the operation, whilst also being expensive, labour intensive and difficult to implement on a production line. Other methods, such as photogrammetry, are operationally similar, yet offer more flexibility. 26,27 By focusing on flexible and scalable sensing techniques, it is possible to merge data sets from both research and industry, including from both physical experiments and simulations. In turn, increased volume and diversity of data fuel the deployment of Machine Learning (ML) techniques to optimise manufacture.
ML has already been shown to be powerful in the context of composite materials. Chen et al. [10][11][12] recently demonstrated that genetic algorithms and artificial neural networks can successfully optimise component in-plane properties for a chosen geometry. Optimisation in this manner however, is often expensive to evaluate for new geometries. Hence others have explored alternative techniques which reduce computational requirements such as Bayesian optimisation and Gaussian process modelling, especially when coupled with convolutional neural networks (CNNs). 16,23,25 Zimmerling et al. 16 highlight that significant computational savings can be made in feature and process level optimisation through optimised offline CNNs. These are, however, still limited by the size of training data sets, with features sufficiently separated. 23 Future approaches utilising feature and process level optimisation in a single package consequently seem likely. Bayesian optimisation, [28][29][30][31] when coupled with Gaussian process modelling, 32 allows efficient optimisation that balances exploration and exploitation of the search space whilst outperforming stochastic gradient based methods in uncertainty quantification. 33 It appears that to improve the development of ML strategies, standardisation of some parameters is needed, in much the same way as was done for bias-extension characterisation. 34 This would allow much more inclusive training data sets that could be extrapolated to wider forming scenarios. An extensive review of ML for production processes is available in Ref. 35. In this paper, an experimental method is proposed that can impose tailored BCs, which activate in-plane deformation modes, and allows for the subsequent evaluation of defects. The process is based on open-mould stamp forming, following similar deformation mechanics to that of DF. Neodymium magnets are used to clamp the material onto a forming bed, allowing interply slip whilst providing greater flexibility in their placement. The authors previously gave a summary of the experimental procedure in Ref. 26. That study showed that tailored boundary control in the bias directions prevented wrinkling. The same preliminary study highlighted that the distance between the magnetic clamps and the tool had measurable effects, in lieu of a changing magnetic force. Further research, undertaken in Ref. 36, revealed distributed clamping was also applicable to forming operations with multiple plies. In the following section, this experimental procedure is summarised for completeness and extended to quantitatively measure defects and to increase the geometrical complexity of the tool.
The magnetic clamping process was first demonstrated over a hemispherical tool, and second, over a more complex cubic geometry including features designed to induce large deformations. Separately, in-plane and out-of-plane properties of a veiled biaxial NCF ( Figure 1) are characterised for the first time to facilitate creation of a parametric Abaqus model. Bayesian optimisation is then employed in conjunction with the FEA model to indicate the placement of magnetic clamps that minimise wrinkling.

Method
Distributed clamping process and experimental framework All tests were conducted with a Tenax biaxial NCF (Teijen IMS65 E23 24K). This material consists of 0°and 90°u nidirectional (UD) fibres with a polyamide copolymer (CoPA) tricot-pillar stitch. Below each UD ply, chopped strand toughening veils (ET205) are present. A powdered binding agent (EPIKOTEÔ Resin 05311-15 gm À2 ) was also applied to the outer veil by the manufacturer. Figure 1 shows a schematic of the material composition and constituent plies [90v0v], with both front and back views. The front is distinguished as the 'fibre side', with stitches orthogonal to the 90°fibre, the back is defined as the 'veil side', faced with a veil and tricot stitches ±45 + to the 0 + fibre. The fibre elastic modulus, as stated by the manufacturer was 290 GPa, with the measured fibre volume fraction of the NCF, V f = 0.33, owing to the fibre architecture in the unconsolidated form.
Material sample preparation followed a rigorous process to reduce variability and induced pre-shear. This was achieved through masking the periphery of the cut lines with adhesive tape and using a rotary cutter to cut square samples with side length of 400 mm. Each specimen was sparingly sprayed with a suspension of hydrated magnesium silicate (Talc, Mg 3 Si 4 O 10 [OH] 2 ) in isopropanol (C 3 H 8 O)in a ratio of 1:4to facilitate image analysis. The alcohol evaporated leaving a fine white powder, which did not affect structural behaviour, contrary to paint. This approach, coupled with the effects of the veil, reduced reflections from carbon fibres whilst also providing trackable features.
The axisymmetric hemispherical geometry implicitly amplified the influence of material anisotropy on forming mechanics and was used for one set of investigations. Using distributed magnetic clamping, a forming strategy was established and extended to a cubic geometry. Through these trials the impact of distributed clamping on the forming behaviour was assessed.
The distributed clamping test-rig developed for the hemisphere (radius = 60 mm) geometry is presented in Figure 2(a)-(c). As the hemispherical tool allowed evaluation of the magnetic clamping mechanism without the influence of local geometrical variability, developing the process with further complexity was desired. Typical aerospace parts usually have many planar sections, particularly in panels, spars and wings, and these parts have more features such as radii and tapers. A cubic geometry was therefore designed to replicate some of these complexities, see Figure 3. This part was a natural follow-up to the hemisphere, with comparable tool dimensions and only positive Gaussian curvature. Further, it was not axisymmetric but had four transverse planes of symmetry.
The rig was fixed to an Instron tensile test machine via mounting holes. A mild steel forming bed (400 × 500 × 3 mm) was attached to the supporting frame with a circular exclusion (radius = 61 mm). An insert with a 10 mm radius was then bonded to the internal exclusion. This part was vital for reducing stress concentrations that would lead to fibre or stitch damage and for efficiently transferring tension between the magnets and tool. The hemispherical and cubic tools were machined with a smooth surface finish from aluminium and steel billets respectively. The hemispherical and cubic tool included a 20 mm tangential extension from the equator to ensure no edge effect in forming, Figure 3. The NCF blank was placed on top (veil side up) of the forming bed and clamped in place using magnets at the specified test locations. NCFs tested were orientated with the material x-direction aligned with the positioner X-direction. This was achieved through use of a magnet positioner with pre-defined discrete occlusions, laser-cut from 3 mm MDF (medium density fibreboard), this was placed for positioning the magnets and removed before initiating Magnet positioner and coordinate system where four experimental test configurations are defined by θ, which represent the four polar coordinates the magnetic clamps take: ½r,θ, ½r,θþ90,½r,θþ180,½r,θþ270. The material coordinate system is shown in red and (d) shows the cubic geometry used in this study, which was machined from a mild steel billet.
the forming procedure. These corresponded to 16 angles at 22:5 + increments and four radii at r = 90, 120, 150 and 180 mm, which were expressed in polar coordinates [r,θ]. For the purpose of this study, only one radius was used, r = 150 mm. The influence of the radial position of the magnetic clamps was explored in Ref. 26 as mentioned above. Given the material symmetry, four magnetic clamps were the minimum required to promote orthotropic inplane shear deformation. We defined this arrangement with the coordinates ½150 mm, θ þ i90 + , where i ¼ 0, 1, 2, 3. We performed four separate tests where θ ¼ 0 + , 22:5 + , 45 + , 67:5 + . The material architecture suggests the simplest way to promote in-plane shear deformation is when θ ¼ 45 + , as illustrated in Figure 4. This works similarly for any geometry with four transverse planes of symmetry. Tricot stitch architecture is distinct from chain stitches, 37 with connections between each stitch line. The stitch architecture is shown in Figure 4, both in natural, and deformed configurations. In the experiments, 20 mm diameter circular N42 neodymium magnets with a thickness of 5 mm were used. The vertical and shear pull, as rated by the manufacturer, are 7.30 kg and 1.46 kg respectively when in contact with a mild steel surface, however as no test standard or specification is stated, calculating the compressive force from these values is not valid. Our own testing found a higher compressive force of 90 N. This was achieved by bonding a magnet to an Instron tensile test machine platen (non-ferrous), placing this fixture on top of the steel forming bed and subsequently raising the Instron fixture to measure the reaction force.
The test procedure was as follows: First the test rig was attached to an Instron 3369 tensile testing machine via the mounting plate and brackets. The mould tool was then fixed to a 1 kN load cell via a mounting bracket and pin. The NCF blank was placed centrally on the forming bed and orientated to the required position. Following this, the magnet positioner was employed to place the magnets at the test positions and removed once the magnets were all placed. The tool was first placed in contact with the NCF sheet and then lowered at a constant rate of 10 mm/ min, moving a total of 61 mm along the negative Z direction, see Figure 3. As the experiment progresses the horizontal offset between tool and the radial insert decreases. During the final stages, this results in compaction of the fabric between tool and radial insert. This effectively clamps the material, promoting increased tensile forces that draw the material closer to the tool surface, allowing small wrinkles to re-orientate. On completion, both sides of the NCF were photographed using a Canon EOS 2000D SLR camera, ensuring full coverage of the specimens with 50-100 photos. Analysis through photogrammetry (Autodesk's ReCap) produced a 3D reconstruction of the surfaces in the form of point clouds. These were subsequently scaled with the use of calibration markings on the radial insert, see Figure 5.

Extracting wrinkle morphology
Whilst photogrammetry is less accurate than 3D DIC, it is advantageous in other ways. Photogrammetry allows optimisation of image acquisition through its non-fixed set-up and therefore more data can be collected around particular features of concern. Readily available cameras permit reconstructed 3D models with high quality texture mapping, aiding analysis. We previously demonstrated the experimental capabilities of photogrammetry in Ref. 26, which has since been corroborated by Harrison and Gonzalez Camacho, 27 in a comparative study with white light scanning. A parallel study, outside the main scope of this work, showed photogrammetry achieved sub-mm precision over similar NCF preforms. 38 Photogrammetry's fast deployment is more aligned with industrial capability and can be used at any scale. High quality strain information is less accessible to photogrammetry, particularly as this method only captures the final material morphology, but accurate . The hemisphere and cubic geometry used in this study, which were machined from an aluminium and mild steel billet respectively. Radial insert schematics, with a 10 mm radius, were bonded to the forming bed and are shown taking a cross section through the Xx-axis ( Figure 2). measurement of wrinkle size can be achieved through numerical analysis of the point cloud. The process used follows Figure 5, where the top surface is shown to indicate the plausibility of image analysis to understand the magnet trajectories. The workflow is also amenable to automation and shows high quality scanned surfaces with texture mapping. As hemispheres can be described analytically, it is trivial to calculate the distance of the formed part from the origin of the hemisphere. This was measured at each coordinate in the point cloud, α. The distance between α and the tool surface, for every node was then calculated using the tool radius and original material thickness (1 mm) to give, δ.
The calculation of this metric was modified to deal with geometries that cannot be described analytically. Specifically, in those instances, a triangular mesh representing the tool surface, β, was used in MATLAB to obtain normal vectors, n ! at each node. The intersection of the normal, n ! , with the scanned or simulated final shape, α, was used to calculate δ. 39,40 This technique failed to capture the complexity of wrinkles that included multiple folds, or in distinguishing between poor tool/fabric conformity and wrinkling, but it produced a simple estimate of defects. Further development, utilising local curvatures could further improve defect classification. We further leveraged the metric, δ, by calculating δ max and δ 75th (75th percentile of δ) for a given instance of the test, allowing for further analysis of the distribution of results.

Biaxial non-crimp fabric material characterisation
In-plane material characterisation. Bias extension tests with an aspect ratio of 2 were conducted according to the method described in Ref. 41 and 42 on an Instron 3369 tensile testing machine. The fibrous nature of the NCF did not allow successful DIC strain information to be gathered in this instance, hence axial strain information was based on output from the Instron, which included the displacement of the sample edge. Using the axial strain data from the Instron . Symmetrical shear mechanisms present in the tricot-pillar stitched NCF used in this study. When θ ¼ 45°(inset) shear modes are generated by the magnetic clamps. When θ ¼ 0°no in-plane shear mode is promoted. Dark grey is 90°fibre tows, light grey 0°fi bre tows, green is the chain stitch on the front, red dashed line is the tricot stitch viewed through the material (back). The key to the right highlights stitch architecture and shows where the unit cell of the tricot stitch is connected at points (1) and (2). therefore produced an overestimate of the shear angles at high strains. This was quantified using ImageJ, (Fiji), 43 to extract the shear angle and deviation from the theoretical value, Figure 6. However, as this deviation is also attributed to edge effects, such as stitch breakdown and slip, the kinematic model was deemed suitable as these would break the pure shear assumption. While slip is not modelled, we can hypothesise that, in the physical forming trials, some amount of slip will occur between the shear angles of 30°and 60°. Whilst this study did not aim to develop a state-of-the-art material model, macroscale modelling results presented later highlight the efficacy of this simplification. The resulting normalised relationship between shear stress and shear angle was fitted via MATLAB's fitting toolbox as a rational polynomial. To allow the development of simplified hyperelastic behaviour at high strain, a second exponential fit was added to the shear stress formulation. The shear stress, illustrated in Figure 6, therefore had the following formulation (2) with coefficients Out-of-plane material characterisation. Unlike continuous blank holder restricted processes, 14 where M is the weight per unit area and C is the bending length. Rectangular strips (30 × 250 mm) were used for warp and weft orientations, with a larger width (100 × 250 mm) used for the bias orientation. 46 Six test repetitions were performed to eliminate the influence of experimental variability. Strip displacement was measured using a steel rule. Whilst this method could not evaluate non-linear bending, 21 it did offer an approximate representation for the secondary deformation mechanism in distributed clamping.
The bending experiments showed that the veil and stitch architecture influenced the out-of-plane mechanics. This can be evaluated from Figure 7 orientation 2, whereby having the 0°fibres on the top surface gave a bending stiffness nearly half the magnitude to orientation 3. This is not only attributed to displacing the 0°fibres 0.5 mm above the fulcrum, but also to the addition of the veil and binder. When the veil sandwiched a 90°fibre layer on the compressive bending region [0v90v], it acted as a stiffener in preventing rotations of the tows. When the veil was removed from this location [v0v90], the 90°tows were free to rotate, significantly reducing the bending stiffness. Further, Figure 1(d) shows binder particulates appear to have partially cured and fused with the fibre sizing. This would increase the compressive stiffness on the mesoscale, allowing load to transfer between the 90°and 45°fibres when in compression (see Figure 7, orientations 4-5). Without the binder or veil, compressive resistance comes from the stitch. This may have important implications for macroscale homogenisation strategies that aim to capture the full effect of the anisotropy.

Finite element model
A dynamic explicit analysis was chosen in Abaqus for all FEA models. To model the NCF efficiently, we assume an elastic material with a zero Poisson's ratio in the direction of the fibres and we ignore tow slip that can often emerge at large shear. These assumptions enable a simplified homogenisation of the material properties. By definition, membrane models cannot capture the bending behaviour. Thus, a coupled membrane and shell element approach was pursued in this study similar to Ref. 3. Using a shell model in solitude, would not allow the decoupling of the in-plane and out-of-plane behaviours. To achieve this in Abaqus, shell elements and membrane elements were defined, whilst sharing superimposed nodes. This is distinct from coupling or tying and allowed efficient computation of the bending response without affecting the constitutive in-plane behaviour. 3,16,[23][24][25] First, we discuss the membrane material model, followed by the shell material model.
Non-orthogonality emerges as the material shears; this necessitates updating the material properties with respect to the new fibre orientations. Towards that end, a FORTRAN VFABRIC non-linear material subroutine was chosen. VFABRIC calculates and updates the incremental fabric stresses for the given incremental fabric strains. The constitutive equation for the material, following our assumption about the Poisson's ratio, takes the form shown in equation (4).
Where the directional subscripts 1 and 2 represent the fibre directions, which initially align with the material coordinate system direction x and y for the undeformed sheet. Here, dσ is the in-plane incremental stresses, dε, the in-plane incremental direct strains, dγ, incremental in-plane shear angle and G 12 , the shear modulus. The Young's modulus applied to the membrane elements in the fibre direction was calculated from the given fibre properties after accounting for the fibre volume fraction. To allow the shell model to work in parallel with the in-plane membrane response, 3 the elastic modulus was satisfied according to equation (5).
Where the superscripts, T, M and S represent the total, membrane and shell quantities respectively. This accounts for the intrinsic in-plane stiffness of the shell element, by subtracting it from the total in-plane stiffness of both elements. Using VFABRIC, the direct incremental stresses are calculated as shown in equation (4), but the shear stress, σ 12 , is updated by making use of equations (1) and (2) for the current shear strain, γ 12 , without using incremental stresses and strains or the need to calculate an equivalent non-linear shear modulus G 12 . This model has been shown to successfully represent orthotropic fabric materials, such as woven glass fibre. 12,47 In this instance, this results in an overestimation in shear stress, σ 12 , (∼30°), similar to the behaviour observed in Ref. 6. This may be due to edge effects and slip in the bias extension experiments. The resulting membrane material model was used in conjunction with the 4 mm (M3D4) membrane elements. The bending material model used S4 shell elements defined using the same nodes as the M3D4 elements in Abaqus. A single orthotropic lamina was used to define the material properties. The bending properties observed in the cantilever tests above are represented by an equivalent, average Young's modulus along the fibre directions. Specifically, the mean of the orientations where the long fibres were orthogonal to the fulcrum (see Figure 5, cases 1-4) was used. This dictates the bending model cannot capture the through thickness anisotropy and instead the model homogenises the warp and weft bending stiffnesses in an isotropic bending model. A mean bending stiffness, G, was calculated as 3.007 × 10 À3 N m for a bending length of 164 mm -see equation (3). Then, Bernoulli-Euler beam theory was deployed to estimate an equivalent in-plane Young's modulus. This was calculated by E S ¼ 12 G t 3 and produced an equivalent modulus of 36.09 MPa with t = 1 mm for the NCF. The shell model also required transverse shear stiffness properties, which were set at 15.04 MPa, following Abaqus' documentation on shell section behaviour. 48 With these properties, a cantilever simulation was performed, similar to Ref. 3, where a bending angle of 41:5 + was found. Further simulations were undertaken with the fibres in the bias direction (Figure 4 cases 5 and 6) for validity. These gave a 3.9% underestimate in bending angle, 39:9 + for an 127 mm bending length, which was the experimentally derived mean of bending length for the bias direction samples. This is to be expected as the cantilever experiment neglects the contributions from slip, shear and fraying at the edges, which would be exacerbated in the bias orientation.
A dynamic explicit analysis again was chosen to model the distributed magnetic clamping experimental procedure. This allowed for better contact stabilisation in comparison with a standard analysis and allowed the use of the Abaqus VFABRIC material subroutine. In pursuit of a model that could easily be adapted, a parametrised Python script was written that allows the user to modify many of the key variables. The magnets were modelled as analytical rigid parts with a mass equal to 35 g, a radius of 10 mm, a 0.5 mm fillet and a height of 5 mm. The penalty contact method was chosen, following a Coulomb friction model, with the magnets free to slide relative to the forming bed, whilst allowing the material to also slip. A coefficient of friction equal to 0.2 49 was used for all surface interactions, except for the interaction between the magnets and the veiled upper surface of the NCFthere 0.22 was used. The higher coefficient represents the higher frictional response of the magnet compacting the upper material surface. The model has two steps: the first initialises the clamping forces to promote numerical stability in the initial increment; the second step is the lowering of the tool, which emulates the tool stroke in the physical experiment for either geometry. With the natural time scale preserved, a small stable time increment could ensue due to the relatively small elastic modulus used for the shell elements impacting the wave speed propagation, thus fixed mass scaling was applied at 10 5 . All simulations were terminated after 361 s, indicating successful draw. The simulations were run on the University of Bath's high performance computing (HPC) service Balena, utilising 1 Intel Ivybridge node (7 cpus with 2.6 GHz).

Bayesian optimisation
In this section, a brief introduction to Bayesian optimisation is presented. A more detailed background to Bayesian optimisation can be found in Ref. 31. This method is used in this paper to find the optimal magnet position. In a general Bayesian optimisation framework, the aim is to minimise a so-called black-box function, denoted by f ðθÞ, of varying inputs θ. A black-box function refers to the situation in which the process for obtaining the function output at a given input is unknown or difficult to model. In this paper f ðθÞ denotes the material conformity at magnet angle θ, in the absence of any noise. Bayesian optimisation is able to minimise such functions by maintaining a statistical model representing a belief about the function based on data δ 1 ,…,δ n from past experiments, corresponding to noisy versions of the function f ðθÞ at tested inputs θ 1 ,…,θ n . The model is updated sequentially with the acquisition of new data, and the next test input is proposed based on improved knowledge of the minimiser, Figure 8. The model used in this paper is a Gaussian process model, described in more detail later. The data also provides information about the overall shape of the Gaussian process surrogate, which reduces the uncertainty in the model, but at the same time, as the simulations are computationally expensive, resources should be spent towards locating the minimum accurately. Of course in practice, this method is not bespoke to only computational evaluation of the objective function, but could also be used with physical experiments in a parametric trial. The learning of the function is known as exploration and the locating of the minimum as exploitation. A decision about which input to test next is made in a way that balances exploration and exploitation. At the early stage, when not enough observations are available, exploration is more desirable. The statistical model carries some uncertainty about the value of the function at untested inputs. As more data is obtained, the uncertainty is reduced, leaving the focus to become exploitation. Figure 9 illustrates the algorithm used in this study. The application of this method for the optimal positioning of the magnets is presented in the results section.
Gaussian process surrogate. A statistical model is a relationship between inputs θ and their corresponding random outputs, δðθÞ. In this paper, we assume δðθÞ ¼ f ðθÞ þ eðθÞ, where f ðθÞ denotes a structured component and eðθÞ denotes a random perturbation around the value f ðθÞ that is assumed to be normally distributed with mean 0 and standard deviation τ 1 . A Gaussian process model 32 is used to represent the function f ðθÞ. The Gaussian process model is a class of statistical models that contains functions with normally distributed outputs. A Gaussian process model is characterised by a mean and a covariance function (also called the kernel); the former corresponds to an average value in the absence of any data, and the latter models the correlation between different inputs. We denote the mean by μ and the kernel by Kðθ,θ 0 Þ, where θ and θ 0 are two inputs. The kernel can be used to impose a certain behaviour to f ðθÞ. For example, the choice imposes periodicity to the function with period ω. The notation jθ À θ 0 j means the distance between the inputs, τ 2 is the standard deviation of the normally distributed output, and f is a scaling parameter which determines how much the function is allowed to fluctuate, and generally, the further θ and θ 0 are, the less information about f ðθ 0 Þ is provided by f ðθÞ. For the statistical model used in this paper, we set the noise standard deviation to τ 1 ¼ 0:5. Due to the inherent periodicity at 90°due to the magnet positions and the symmetrical shear properties leading to material symmetry at 45°, the following kernel was chosen where K p ðθ,θ 0 Þ is the periodic kernel given in equation (6). The Gaussian process standard deviation is set to τ 2 ¼ 2, the period to ω ¼ 90, and the scale parameter to f ¼ 0:5. A sensitivity analysis of τ 1 , τ 2 and f revealed the global optima had little sensitivity to these parameters. The Gaussian process mean is set to μ ¼ 4 for the hemispheric geometry and to μ ¼ 10 for the cubic geometry.
Bayesian optimisation process. Suppose n experiments at inputs θ 1 ,…,θ n gave correspoding outputs δ 1 ,…,δ n . The Gaussian process model was used to provide a prediction for an arbitrary input θ. As mentioned at the beginning of this section, the objective was to propose the next experiment/input θ nþ1 that would help locate the minimum of the function. This is determined by evaluating what is known as an acquisition function at untested inputs and choosing the input with the highest value. Suppose that the current minimum value after the n experiments isδ n . Then when a new test at θ is conducted, the acquisition function calculates how much the minimum value is expected to decrease after seeing the new outcome δðθÞ. Of course, the tested input θ must be decided before the experiment, i.e., before seeing δðθÞ, which is why it can be regarded as the expected decrease, or more commonly, the expected improvement. Mathematically, the expected improvement (EI) at a test input θ is given by Here E stands for "expected value", i.e., an average over all possible outcomes of the term in the square brackets over the uncertain output δðθÞ for fixed θ. Since the outcome is unknown, all possible outcomes are considered and then averaged. The term maxf0,δ n À δðθÞg is the improvement over the current minimum value, assuming that δðθÞ is observed. If it turns out that δðθÞ is smaller than the current minimum (δðθÞ <δ n ), then we have improved by an amount δ n À δðθÞ. However, if δðθÞ turns out to be larger than or equal to the current minimum (δðθÞ ≥δ n ), then there is no gain, so the improvement is 0. This proceeds sequentially for n ¼ 1; 2,…, each time computing equation (8) for all the available experiments, and choosing the one that corresponds to the largest value, i.e., θ nþ1 ¼ argmax EIðθÞ. This is illustrated in Figure 8. In each case it was proposed to test close to the existing minimum, taking into account the function uncertainty. Here, we implemented this optimisation approach through a Python script. This was flexible to which statistical input parameter we chose.
In this work, the magnetic clamping process is treated as a black-box function, with the input θ corresponding to the positioning of the magnet, θ 2 ½0 + ,90 + , and the output corresponding to that choice of θ being the wrinkle metric δ max ðθÞ. This metric was chosen as it was less susceptible to the exclusion of slip in the FE model. We considered angles at 2:5 + steps, i.e., we choose among θ 2 f2:5 + ,5 + ,7:5 + , …,90 + g, with the first experiment being at θ 1 ¼ 2:5 + . The optimisation procedure was implemented via a python script that follows the flowchart in Figure 9.

Application of the distributed clamping process
For the hemisphere geometry ( Figure 10 and Table 1), large wrinkles were mitigated when θ ¼ 45 + . In the fibre mode example (when θ ¼ 0 + Þ large macroscale wrinkles were present in low tensile regions. These low-tensile regions do not allow shear deformation as the material is drawn in, causing the initiation of compression-folding which prevents conformity. In the mixed-mode experiments, when θ ¼ 22:5 + or θ ¼ 67:5 + , twist of the fibre tows was observed. Wrinkles however, were again largely mitigated in comparison with the 0°case, although both show evidence of tow slip. In the case of the cubic geometry, Figure 11 and Table 1, when θ ¼ 45 + , forming over the complex geometry highlighted the plausibility of the forming technique, with wrinkles above 2 mm reserved to the bottom of the formed part indicating compression at the radii. The 75th percentile of deformation shows the majority still conformed well. When examining the material qualitatively, wrinkles were reduced over the main section of the geometry, with a wrinkle-free region over 55 mm of stroke depth (St). Fibre tow compaction was also present, Figure 11(b), due to the large strain required to form the corner. This is due to the tricot stitch illustrated in Figure 4. When θ ¼ 0 + wrinkles were largely symmetric over each corner. This represents two areas of material which fold back in on themselves as the shear mode is not activated. When θ ¼ 22:5 + or θ ¼ 67:5 + , twist was again observed, with less impact on the fibre direction. This is possibly due to a decoupling or reduced dominance of in-plane tension, with the contact point from each magnet not coincident, as it is in the hemisphere (pole). Wrinkling differs, with major macroscale folds now asymmetric on each corner. Again, both mixed-mode experiments showed more evidence of fibre slip, as is to be expected from tension-shear coupling.
It is evident from the experiments that distributed magnetic clamping had a strong effect on forming for both hemisphere ( Figure 10) and cubic geometries ( Figure 11). In-plane radial tension was generated between magnetic clamps and points of contact on the tool. Tangential forces around each magnet were also present, however these seemed to be less dominant as they induce bending, which relieves the corresponding stresses. Magnetic clamping was also shown to control the material draw-through when activating the shear mode, (θ ¼ 45 + Þ, for both the hemisphere and the more complex cubic geometry. Discontinuous blank clamping, unlike continuous blank holder clamping, 6 allowed the material to deform with out-of-plane bending (Figures 10 and 11). Hence reducing the available energy that could induce mesoscale wrinkling. This may indicate an important method of facilitating optimum shear in a similar way to manual folding and tension shearing in hand lay-up. 50 FEA models that incorporate bending behaviour, could therefore be necessary to fully unlock the forming potential of NCFs in this manner. Figures 10(c) and 11(c) show that in the mixed-mode experiments, θ ¼ 22:5 + , fibres slip to accommodate the differential force induced by the unbalanced boundary condition. Whilst this could lead to defects, particularly post-cure, this does indicate that magnets could be utilised to induce steered fibre tows and non-standard deformation.

Abaqus process simulation
Axial force values, F a , from the experimental studies were compared to the virtual process model; Figure 12 shows the comparison when θ ¼ 45°. Here, it is seen that the general   Figure 13. The results indicate that the macroscale forming mechanics have been successfully reproduced. In the optimal cases, when θ ¼ 45 + , the hemisphere FE results show excellent agreement with the experiment, however the cubic case shows some wrinkles are present on the faces which were not in the experiment. This may be a result of not including slip in the model, instead the model must accommodate the curvature of the corner by pulling more material over the faces. Comparisons when θ ¼ 22:5 + and θ ¼ 67:5 + again reveal that the major deformation mechanics have been  reproduced. The 's' shape present in these experiments was not as clear in FE, although again this is mainly driven by slip in the physical experiments. The direction of the folds and the locations of material thickening in both geometries do remain good experimental matches. When θ ¼ 0 + , for the hemisphere, the locations of the wrinkles approximately matched, however the wrinkle pattern differs. The FE model here appears to highlight the symmetrical shear property that was assigned to the macroscale model that prevents the model capturing these slight material asymmetries. In the cubic example, θ ¼ 0 + , the key deformation has been reproduced, albeit the wrinkle pattern on the corner shows differences in folding. Figure 13 does show these differences in macroscale folding, which is possibly due to not modelling the through thickness bending anisotropy.
The in-plane shear stress results for the hemispherical tool configurations when θ ¼ 45 + are shown for every 15 mm of stroke depth in Figure 15. It is evident that magnets activated the shear deformation modes, as explained in Figure 4. Here, shear stress initiating from the  change in curvature at the forming bed radius increased (by St ¼ 30 mmÞ as the tension induced by the magnets increased. Subsequently, in-plane shear modes were promoted as the tension between the magnet and hemisphere pole continued to increase. At the same time, out-of-plane bending outside of the magnetic tensile region was initiated, allowing the material to accommodate the shear deformation. This facilitated the material to stretch and develop positive Gaussian curvature. When θ ¼ 45 + , the direct fibre modes 26 have relatively little tension, hence these conform to the tool with direct fibre stresses, σ 11 and σ 22 . Mild wrinkling was observed at locations with no wrinkling in the experiment. We hypothesise that slight deviations represent some of the homogenisation strategy for the material model, particularly the exclusion of tow slip and the coarseness of the mesh imposed by computational limitations. Here compaction was not modelled efficiently and has been shown to be significant in similar heavy, thick stitched NCFs. 51 This cannot be represented by standard plane-stress elements. Further, from the out-of-plane bending behaviour, we did not model the throughthickness anisotropy. These factors however did not alter the major characteristics of the deformation, particularly in the membrane-driven deformation area around the tool geometry. In the physical experiment, the compaction of the fabric between the radial insert and the tool is thought to drive the final important stage in inducing peak in-plane tension and shear. In the simulation this was also present, but to a lesser extent due to a reduction in compressibility. This may have led to an underestimation of the effect of compression at this interface, reducing the in-plane tension generated across the material, and reducing conformity to the tool. As can be seen in the cubic case, when θ ¼ 45 + , this higher tensile period could have resulted in a better experimental match.

Application of Bayesian optimisation to the distributed clamping process
The optimisation procedure for the two geometries took 5.25 days (hemisphere) and 9 days (cubic) to complete nine simulations. The discrete rigid description of the cubic tool explains the increased computational demands. In Figures  16 and 17 it is noted that after five simulated experiments, both the cubic and hemisphere geometry results already exploit the area around θ = 45°, indicating this as the optima. After nine tests, both the hemisphere and cubic geometries minima is clear, with relatively little uncertainty around this point. In more complicated geometries, the optimisation procedure may have to be run for further iterations, here a stop criterion could be beneficial. This could declare an optimal result when selected θ values are consecutively within a predefined tolerance. The Gaussian process approach highlighted its intrinsic ability to balance exploration and exploitation, and through each stage it demonstrated the ability of the algorithm to select the magnet location that offered the most probable expected improvement. Figures 16 and 17 also illustrate δ max for the experimental results, plotted for comparison after nine tests. There appears to be phase changes in response throughout the domain, particularly for in the cubic case. These could indicate where wrinkles changed between a larger amplitude wrinkle and multiple smaller wrinkles. This eludes to the need for a metric or set of metrics that would be able to encompass many more statistics to fully characterise the forming behaviour and resulting deformation. Metrics such as wrinkle area and curvature have previously been demonstrated in Ref. 37. Statistical inference models could therefore allow combinations of metrics in an automated post-processing application. Whilst the optimisation had a feasible run-time, this could be further reduced to improve the utility and allow efficient global optimisation models based on process segmentation and image analysis.

Discussion
Magnetic control was simple and flexible in application. By decreasing the variability in the manufacturing process, coupled with increasing the flexibility of boundary control, more translatable data sets can be built. Bayesian optimisation was presented here to highlight the ease of optimisation when using distributed clamping. This technique opens the door to more control with variables such as magnet position radius, number of magnets, magnet angle and magnet size. While these could be seen as increasing the search domain, feature segmentation could further simplify complex geometries, reducing the sampling burden and increasing diversity of data.
Whilst tools here were closed form and had at least four transverse planes of symmetry, magnets could readily be deployed to control shear deformation in parts such as aircraft secondary structures. Here, magnets could be used Figure 16. Bayesian optimisation stages for the hemispherical tool using FE simulation results, δ max , (red) and comparison with experimental δ max (green) after nine tests. The mean of the distribution is represented by the thick blue line (current estimate) and the uncertainty denotes the 95% confidence interval for the value of the function.
to control material-draw in around particular features by deploying specific magnet configurations. Magnetic constraints could also be used with matched tooling forming processes. Further, a variety of magnet shapes such as bar shaped magnets could allow control of large sections of fabric in drape forming.
While other process control measures exist, few offer the same ease of application. Here, a flexible process has been described. This will allow manufacturers and researchers alike to develop comprehensive training data banks, fuelling predictive modelling capabilities of the future.

Concluding remarks
Magnetic distributed clamping offered an efficient process control measure for forming positive Gaussian curvature parts. Providing an effective means of boundary control, with high levels of repeatability, in-plane magnetic constraints were effectively modelled and amenable to optimisation. Feature segmentation and feasible offline optimisation would however require further characterisation of the material mechanics. This study has therefore resulted in the following remarks.

Distributed magnetic clamping successfully induced
in-plane tension to alter the deformation mechanics and highlighted the advantages of magnets for application in fabric forming. Promotion of out-ofplane bending in the material, because of discontinuous blank clamping, enabled a reduction in defects when the shear mode was activated. A complex cubic geometry with industrial applicability was successfully shown through experiment to have good material conformity. 2. The non-fixed non-stereo photogrammetry technique was successful at capturing the surface morphology of the material and allows further flexibility and evaluation of otherwise inaccessible regions and could be used in industrial settings. 3. Characterisation of the veiled biaxial NCF revealed complicated bending mechanics. The veil was shown to act as a stiffener in bending, with highly Figure 17. Bayesian optimisation stages for the cubic tool using FE simulation results, δ max , (red) and comparison with experimental δ max (green) after nine tests. The mean of the distribution is represented by the thick blue line (current estimate) and the uncertainty denotes the 95% confidence interval for the value of the function.
anisotropic bending properties. However, the efficient macroscale FEA strategy used in this study, with coupled idealised shell and membrane elements, remained a viable approach. 4. Gaussian process modelling of the amount of deformation, using a periodic kernel, and the use of Bayesian optimisation to find the optimal angle, were demonstrated. Results indicated that the magnetic distributed clamping process was amenable to surrogate modelling. However further statistical inference between metrics may be needed to reduce unexpected noise which results in local minima in more complex geometries. 5. This novel magnetic forming method allowed more flexible, scalable, and industrially relevant processing. The method would work with many existing processes and facilitates integration and extension of data sets.
Future research looks to explore the distributed magnetic clamping method with multiple plies, broadening industrial relevance, and to demonstrate the reduction in cycle times that could be attained. Magnetic clamps have been shown to offer broad possibilities that warrant future investigation including more complex geometries, fibre angle offsetting, magnetic clamping on the tool rather than the boundary, multi-stage variable clamping, active magnetic shear actuation and diverse magnet sizes. Such a process in combination with powerful non-parametric models may facilitate the creation of large databanks, leading to dramatically improved flexibility in design and manufacturing.