Modelling defect formation in textiles during the double diaphragm forming process.

As the requirement for low cost, high volume production of composites increases, so does the requirement for modelling capability to help inform and optimise the composite manufacturing processes. In this paper a finite element approach is introduced to simulate the behaviour of a plain weave fabric during the double diaphragm forming process. Using a mutually constrained shell-membrane method to capture the characteristic behaviours of the textile, the approach is assessed in its ability to predict the severe deformations present within experi- mental trials. The model is validated against experimental results for different fibre orientations and for the simultaneous forming of multiple layers. The approach accurately predicts the severe wrinkles observed in the experiments and is used to help further understand the catalysts for defect formation in the double diaphragm forming process.


Introduction
Dry textile reinforced composites are gaining increased attention in high performance industries as the manufacturing processes involved offer increased production rates and reduced costs. Unlike their preimpregnated alternative, the material is formed and consolidated prior to the injection of resin, consequently there is very little cohesion between the individual fibres and yarns, with the weave architecture being the main factor that binds them together. The result is a more formable material, which lends itself well to the forming of complex shaped parts from flat preforms.
Composite components are typically made from multiple layers of fibrous material, these are laid up at different fibre orientations to optimise the structural performance of that part. The forming of these layers can be done successively or simultaneously. The more attractive of the two approaches is simultaneous forming of multiple layer as the procedure can be reduced to just a single step, improving production time significantly.
Double diaphragm forming (DDF) is a method capable of multi-layer forming which has potential for low cost, high volume applications. In the DDF process, layers of material are compressed between two flexible diaphragms and then deep drawn, under hydrostatic pressure, to shape. While the process has largely been associated with pre-impregnated preforms [1][2][3][4] they have more recently been used to produce binder-stabilised dry fabric preforms ready for infusion of liquid resin [5].
One of the main benefits of DDF is the reduction in material waste as, unlike blank holder processes, no excess material is required to constrain the preform during forming, which makes net-shaped preforms a possibility. Furthermore Krebs et al. [6] highlighted that through optimising the ply shapes, instabilities such as wrinkling and out-of-plane buckling can be minimised, as the redundant areas contributing to compressive hoop stresses are removed. The process of optimising ply shapes for DDF has been further enhanced by the development of numerical capabilities to aid in the optimisation process by Chen et al. [7]. In this work finite element analysis was used to optimise blank shape in order to minimise defects, feasible locations of darts (local cuts to the preform) were also identified to reduce fabric bridging and improve surface conformity.
The finite element modelling of textile materials during the forming process is typically performed at one of two scales. Firstly, the mesoscopic scale where the textile is considered as a construct of individual tows (hundreds or even thousands). To correctly model the behaviour of the textile at this scale the model must satisfy two key requirements; an accurate description of the meso-scale architecture and a constitutive model which describes the behaviour of the individual fibrous tows. If such a model is achieved, under the correct loading and boundary conditions, the macro-scale behaviour of the fabric is naturally present through inter-tow contacts and the tow behaviour. This is an attractive approach and one that has been actively perused with success [8][9][10][11]. However, the computational expense of these models is significant and limits their use for many larger scale applications. Continuum modelling approaches, which consider the fabric to be a continuous structure, are more efficient and commonly used. With these approaches the characteristic behaviours of the underlying structure is represented through constitutive models. Original attempts to model the material behaviour, using continuum methods, focused mainly on the anisotropy of the material behaviour and how this is modified through shear transformations [12]. Using these original approaches as a foundation, recent efforts have become more comprehensive, including a wider range of characteristic material behaviours to better describe the complexity of woven materials. Now, models are able to include specific behaviours from in-plane bending stiffness [13,14], which resists abrupt changes in fibre directions, out-of-plane bending stiffness [13,[15][16][17], key for capturing the formation, shape and size of out-of-plane wrinkles, to shear-tension coupling [18] and more recently the irreversibility of some of these behaviours [19].
While the inclusion of these behaviours and their effect on the forming of a single layer has been analysed in depth, the effect of stacking and forming multiple layers of different ply orientations simultaneously has seen relatively little attention. Recent experimental studies have shown that the simultaneous forming of multiple layers induces and/or magnifies the severity of wrinkles with the more extreme deformations occurring between plies of different orientations [20][21][22]. The dominant factor being the relative movement of the layers due to the deformations of the individual layers. A numerical study performed by Guzman-Maldonado et al. [23], which examined the interactions between adjacent layers during stamp forming, drew similar conclusions, highlighting that the deformations of contacting plies of different orientations increases wrinkle severity. Huang et al. [24] also showed, through both numerical and experimental investigations, the effect that laminates of different ply orientations has on wrinkle formation during simple bending. Stacks with alternating 0/90 • ± 45 • layers produced severe wrinkles and ply separation as opposed to laminates comprised of just ± 45 • plies, where no wrinkles were present. This was attributed to the in-extensibility of the 0/90 • in the bending direction, causing each of the plies to buckle, dragging their adjacent ± 45 • ply with them.
One of the main drivers of process modelling tools is to identify the occurrence of defects. As discussed so far, there has been substantial work to achieve this, with models becoming increasingly comprehensive to capture defects arising from a range of different mechanisms. To date, this has largely been performed in the context of blank holder forming. In this paper a numerical model is introduced which includes the key behaviours of orthogonal textiles identified to be a prerequisite of forming simulations; high fibre stiffness, low bending stiffness and nonlinear shear stiffness. This modelling approach is applied to the doublediaphragm forming of a plain weave fabric. Focussing specifically on the forming over a tetrahedron tool, the model is assessed in its ability to capture the more severe wrinkling deformations observed in experimental trials for both single layer and multi-layer forming. The possibility to tailor the blank shape to minimise deformations in the DDF process make it favourable for numerical optimisation, but this necessitates a model which is able to capture the deformations accurately. The work in this paper validates a solution for such applications as well as to generate further understanding of the catalysts for defect formation in the process.

Double diaphragm forming experiment
Double diaphragm forming is used to form planar materials into three dimensional shapes. The method creates an even tension across the material, via frictional forces, and allows for unconstrained, out-ofplane deformation. The process consists of three steps. The preform is situated between two elastomeric diaphragms. These are positioned and secured over a hollow vacuum forming box with the forming tool beneath the preform on a liftable surface. In the first step full vacuum is applied to the cavity between the two diaphragms, compacting the material. The tool is then raised so that its base is aligned to the external boundaries of the of the membranes. The final step is then to evacuate the air between the bottom diaphragm and the tool surface. The atmospheric air pressure pushes the diaphragms down, forcing the perform to form to the tool geometry. A schematic of this process is shown in Fig. 1.
The diaphragms used in this experiment were 1.5 mm thick sheets of silicone. The vacuum forming box was 1320 × 620 mm. The sheets were stretched and clamped over the vacuum forming box. A vacuum pump was attached to the top diaphragm so that the air could be evacuated from the cavity between that and the bottom diaphragm. Square fabric specimens were used, the dimension of these were 350 × 350 mm. The fabric was placed between the diaphragms and central to the vacuum forming box. Flow mesh was used to bridge the vacuum from the vacuum connection point to the corner of the fabric, the fabric was in direct contact with the diaphragms apart from at this point. The tool was aligned below so that the apex of the tetrahedron contacted the central point of the fabric. Details of the tetrahedron tool dimensions are shown in Fig. 2. When in place, a vacuum was drawn between the two diaphragms, following this the tool was lifted its full height of 82 mm into  A.J. Thompson et al. place and the air was evacuated from the forming box. Once the fabric was formed, images of the deformed surface were captured from directly above.
Three forming experiments were performed, each considering the forming of a balanced plain woven, carbon fibre fabric on to a tetrahedron tool. The plain woven fabric is constructed of 6 k tows, with 3.85 yarns/cm and an areal weight of 320 gsm. The first two tests considered a single layer of the plain weave fabric, one with a 0/90 • orientation and a second with ±45 • orientation, in respect to the tool. The final case considered two layers of the fabric with a layup of [±45 • 0/90 • ], with the ±45 • layer placed at the bottom. These were formed simultaneously to examine the effect that forming of multiple layers, at different orientations, has on the final deformed geometry.

Single layer forming trials
Frames taken from videos of both of the single layer forming experiments are presented in Fig. 3 and Fig. 4. It can be seen that the deformation behaviour of the two orientations diverge during the second step in the forming procedure -raising of the tool. The surface topology is different in the two cases. The varying surface topology then evolves into out-of-plane wrinkles following the third step, where the air is evacuated from between the bottom diaphragm and the tool. In both cases the wrinkles are orientated radially from the centre of the tool and are aligned with the fibre directions.
Images of the two single ply cases in their final formed state are presented in Fig. 5. The shear deformation in the ±45 • configuration seems to be minimal which is suggested by the net shape of the preform remaining largely unchanged, as opposed to the 0/90 • configuration which experiences significant shape change. However, the wrinkles in the ±45 • configuration are more severe.

Multi-layer forming trials
For the case where two layers of different orientations are considered ( Fig. 6), in test case 3, the wrinkle severity is amplified with wrinkles appearing along each fibre direction, forming radially to the tool. It appears that each layer tries to deform independently to accommodate the change in shape. Wrinkles form in each layer due to their own intraply behaviour, but are forced to accommodate the wrinkles that occur in the adjacent layer due to differing ply orientations. The final formed geometry is therefore a superposition of both the deformed single layer structures, but with magnified wrinkles.

Numerical investigation
A numerical model has been developed to analyse the deformation behaviour of the textile during the double diaphragm forming process. For the development of this model, it is assumed that 2D textile materials can be represented as continuous single surfaced structures using 2D finite elements. This is based on the assumption that the transverse compression behaviour of the material has no effect on its formability and can be considered rigid. Further to this it is assumed that all material behaviour is fully elastic, strain in the fibre direction is negligible, and that there is no coupling between deformation modes.

Fibre tracking algorithm
For biaxial fabrics, the material has two strong fibre directions, which are typically orthogonal in their undeformed state. At 45 • to the fibre direction the material is highly compliant and loading it at this angle causes the material to shear and the fibre directions to rotate. One of the fundamental requirements of textile models, when using a continuous approach, is the tracking of the fibre direction to ensure changes in the anisotropy of the textile are captured and the resulting stresses are computed correctly.
In this work a hypoelastic material model is used to ensure correct computation of the stresses in the frame of the fibre, this approach was first introduced by Peng et al. [12] and has since been used successfully in a number of studies [25,26]. Hypoelastic material models are used to model reversible, non-linear stress-strain behaviour and use a relationship between an objective stress derivative and the rate of deformation to capture this. In Abaqus the objective stress derivative follows the Green Naghdi framework, where the rotation frame is obtained based on polar rotation. However, for textiles, a further rotation is required so that the constitutive model can be applied in the frame of the fibre, this is due to the rigid body rotation of fibres during shear deformations.
In this approach a local 2D orthogonal coordinate system is defined by two vectors g α , where α = 1, 2 signifies the axes relative to the two fibre directions, this coordinate system is the Green-Naghdi (G-N) work frame that the element strains are received in and the stresses are to be returned in. During loading this coordinate system can undergo a rigid body rotation, the rotation tensor R, defining this rotation, can be extracted by performing a polar decomposition of the current deformation gradient F: The rotated local coordinate system g i α is then calculated using: where g o α is the vectors in their initial state. Assuming the initial fibre orientations are aligned with g o α , the actual deformed fibre directions f α are calculate using equation (3).
Strains can then be converted from the G-N work frame to the two fibre directions using the transformation matrix T where: The strain increment dε, which is provided in the G-N work frame on each iteration, is then mapped to each fibre axis using the equation: Following the transformation of the incremental strains to the individual fibre directions, the constitutive matrix C α is used to compute the incremental stress: The incremental stresses are accumulated in the frame of the fibre and then rotated back to the G-N work frame to be passed back to the solver using equations (9) and (10) respectively.
This material model has been implemented as a user material subroutine (VUMAT) for Abaqus/Explicit and is distributed for free from the Bristol Composite Institute (ACCIS) Github page (https://accis.github.io/H ypoDrape/).

Non-linear shear stiffness
The in-plane shear behaviour of textiles is unique, typically at low strain they exhibit very little resistance to in-plane shear. The fibres are able to rotate freely with the only resistance caused by tow bending and inter-tow friction. As the shear angle is increased, parallel tows come into contact and compress. When this occurs the textile experiences rapid stiffening until there is no-more intra-yarn compliance left and the tows are forced to buckle out-of-plane to accommodate for the excess material. The angle at which this occurs is commonly known as the locking angle.
The non-linear shear behaviour of the textile is derived from experimental data obtained from a picture frame shear test [27]. A polynomial regression curve is fitted to this which describes the shear stiffness as a function of the shear angle (γ). Within the hypoelastic material subroutine the shear angle is determined using equation (11), where ε f 12 are the accumulated shear strains along each fibre direction.
The constitutive tensor, C α , oriented along the fibre direction, is updated on each increment based on γ.

Out-of-plane bending behaviour
The flexural stiffness of textiles, unlike conventional solid materials, is significantly lower than the tensile stiffness of the fibre direction. This is due to the low cohesion between individual fibres and yarns which allows for large relative displacements to occur. Textiles are, therefore, highly susceptible to buckling when subject to in-plane compression. To simplify numerical models, it is common practice to disregard bending stiffness and use a membrane assumption [7,25,26,28], however, recent studies have highlighted that the bending stiffness of the material, although significantly lower than the tensile stiffness, plays a key role in capturing wrinkle onset, shape and size [29,30].
Including the low flexural stiffness, while preserving the high tensile stiffness of the fibre directions, in a model is not trivial. Membrane elements have no rotational degrees of freedom and therefore have no resistance to out-of-plane bending. Standard shell elements do have rotational degrees of freedom, however, their resistance to bending is a direct function of the applied Young's modulus. It is, therefore, not possible to separate the high tensile stiffness of the fibre direction with the much lower flexural stiffness -through material definition alone. This problem has been discussed in detail by Boisse et al. [31] in a recent review paper which concludes that a new bending theory specific for fibrous reinforcements is necessary. Steps have been made towards this with the introduction of a number of modified element definitions which have produced promising results, however, these are often difficult to implement within commercial finite element software [15,16]. Laminated shell elements, which are typically available in commercial FE packages, have also been used to decouple bending and membrane behaviours by pragmatically tailoring the properties of each integration point through the thickness [32,33]. An alternative solution is to couple the behaviour of different element types to separate the different material behaviours. Harrison [13] proposed a pentagraphic beam and membrane approach, the membrane elements provide the in-plane shear behaviour while the beam elements provide the axial stiffness and their cross-sectional width and heights are manipulated to produce the correct in-plane and out-of-plane bending behaviours. Another approach presented by Nishii et al. [34] coupled membrane and shell elements to capture the low bending stiffness and high tensile stiffness, here two shell elements are duplicated on the membrane element with offsets applied creating a sandwich structure, the offsets are introduced to increase their effect on the bending stiffness while minimising their effect on the tensile stiffness.
In the present modelling approach, shell and membrane elements are superimposed, sharing the same nodes so that they are mutually constrained, the in-plane textile behaviour is controlled by the membrane element and the out-of-plane behaviour (i.e bending) is controlled by the shell element. The hypoelastic material model is applied to both membrane and shell elements, however, the membrane elements are assigned the tensile moduli of the fibres and the non-linear shear stiffness of the fabric. The shell elements are assigned a Young's modulus, equivalent to the flexural modulus of the material, this is derived from experimentally obtained values of the textile's bending rigidity. The shear rigidity of the shell elements is set to zero and so the fabric shear is controlled solely by the mutually constrained membrane elements. By coupling the behaviour of these two element types, the high tensile stiffness of the fabric is fully separated from the significantly lower bending stiffness.
The flexural modulus along each of the fibre directions is calculated using Equation (12), according to Euler-Bernoulli assumptions, D is the bending rigidity per unit width determined experimentally from a simple cantilever bending test, E is the Young's modulus of the material and h is the thickness of the material. By rearranging this to find E, the flexural modulus can be determined and prescribed to the shell elements.
The use of the hypoelastic material model for the shell elements, as well as the membrane elements, assumes that the bending stiffness is convected along the fibre direction and rotates during any shear deformation. This results in an anisotropic bending response that evolves with any shear transformation. To verify that the flexural stiffness is projected along the fibre direction as expected, a simple cantilever bending test has been simulated, similar to that presented in Ref. [13].
A surface, 10,0 × 25 mm was created and meshed using rectangular elements. Shell and membrane elements were superimposed, and assigned the same nodes. The hypoelastic material VUMAT was assigned to both element types, for this verification test notional properties were used. For the membrane elements a modulus of 10000 MPa was assigned to both fibre directions and a constant in-plane shear modulus of 30 MPa was given. The shell elements were assigned a flexural modulus of 100 MPa for the fibre directions to represent the low bending stiffness, no inplane shear stiffness is given to the shell elements.
Both element types were given a thickness of 0.5 mm and the shell section transverse shear properties, which ABAQUS requires to run shell elements with a VUMAT, were assigned as K11 = 20, K12 = 0.0 and K22 = 20. An acceleration of 9.81 m s − 2 was applied to the model to simulate gravity, no mass scaling was applied, however, mass proportional damping of 30 m s − 2 was applied to the model to reduce any inertial effects.
Two simulations were run, one where the fibre direction in the VUMAT was specified to be 0/90 • for both shell and membrane elements, and the second simulation had the fibre direction specified as ±45 • . Fig. 7 shows the maximum principal stresses in the shell elements, highlighting the distribution of bending stress in the strip. The difference between the 0/90 • and ±45 • stress plots show that the bending stress is accommodated along the fibre directions and that bending stiffness is convected correctly. The 0/90 • fibre orientation has a stiffer bending response than the ±45 • as the fibres are aligned with the bending direction. The bending rigidity of the two models was calculated following Pierce [35], this produced values of 0.001073 Nmm for the 0/90 • and 0.000472 Nmm for the ±45 • , which equates to a flexural modulus of 103 MPa and 45.3 MPa respectfully. This verifies that the model projects the flexural stiffness along the fibre direction as expected and is independent of the higher prescribed tensile stiffness.

Verification of numerical model
In section 3.1.3, it was shown that by mutually constraining shell and membrane elements and applying the hypoelastic subroutine to both, the desired anisotropic bending response could be captured. Two further tests have been performed to verify that the tension and shear stiffness are preserved following the coupling of the two element types.
The boundary conditions of the two tests are presented in Fig. 8. Both models consisted of 5 × 5 quadrilateral membrane elements, with 5 × 5 quadrilateral shell elements superimposed and sharing the same nodes. The properties prescribed to the hypoelastic material model and applied to each element type are shown in Table 1, these are properties deduced from experimental tests on a plain weave fabric which is discussed in greater detail in the proceeding section.
For the tension test, 1% tensile strain was applied to the model along the first fibre direction. The stress-strain response due to this loading condition resulted in a modulus of 40,036 MPa being observed in the model. This is slightly above the 40,000 MPa inputted for the membrane elements, however, this is an expected result as the global response of the material is equal to the sum of the membrane and shell element stiffness. In the present case, the stiffness applied to the shell elements is negligible when compared to the value applied to the membrane elements, but when the flexural stiffness is larger, this value should be subtracted from the required tensile stiffness and the subsequent value should be applied to the membrane elements.
The shear test was set up to mimic the picture frame shear test, one of the commonly used experimental tests to measure the shear rigidity of biaxial fabrics [27]. To verify the outputs of the model, the normalised   Thompson et al. shear force as a function of shear angle is compared to the expected result. The expected result being the one that the non-linear shear stiffness was derived from. The results presented in Fig. 9 show the model to trace the mechanical behaviour well, matching the results from which the input shear curve was derived. This verifies that both the hypoelastic material subroutine performs as expected and that the coupling of shell and membrane elements does not affect its response and that the high tensile stiffness of the fibre direction, low anisotropic bending stiffness and non-linear shear modulus are all present within the model. It is also important to assess the computational efficiency of this approach and the added cost of including bending stiffness. By comparing the run-times of these models with the same models where only membrane elements were used, it has been observed that including bending stiffness, by the introduction of shell elements, increases the runtime by a factor of 1.48.

Numerical analysis of double diaphragm forming
To assess the capabilities of the modelling approach the experimental test cases are modelled to determine if the approach is able to capture the deformations present in the final formed geometries. The approach is then used to assess the parameters which influence formation of the defects in the double diaphragm process.

Model setup
The properties of the plain woven carbon fibre fabric, used in the experimental tests, and required as input for the model are presented in Table 1.
A picture frame shear test was performed on the plain woven fabric to measure the non-linear shear behaviour, and the bending rigidity was measured by performing a cantilever bending test on the material. For simplicity, the tensile stiffness along the fibre direction was assumed to be linear and was calculated under the assumption that the fabric was perfectly balanced, the crimp is negligible and the as-woven geometry had a fibre volume fraction of 0.4. A number of studies have highlighted that the tensile behaviour of fabrics does show some non-linearity [36], due to the crimp in the yarns, but this non-linearity is negligible compared to that observed in the shear behaviour and is commonly ignored in forming simulations.
For each of the simulations presented in this manuscript, the fabric was modelled using square elements aligned to the initial, orthogonal, fibre directions. A high resolution mesh was used with an element edge length of 1 mm to ensure the mesh was refined enough to capture the wrinkle shape and size observed in the experiments.
The modelling of the double diaphragm forming process follows that presented in Ref. [7]. The full diaphragms were modelled using Ogden's hyperelastic material model, which is a predefined material model available in ABAQUS. The input parameters have been experimentally determined for a silicone diaphragm in Ref. [7] and are listed in Table 2. Linear shell elements, with reduced integration were used for the diaphragms. A gap of 0.7 mm was placed between the internal surfaces of the diaphragms to allow placement of the textile. The external boundaries of the diaphragms were fixed in all degrees of freedom. The interactions of the diaphragms with the fabric were captured using the general contact algorithm available in Abaqus Explicit. The membrane elements representing the in-plane properties of the fabric were used as its contact surface and the shell elements were removed from the contact algorithm altogether, this ensured the contact was not accounted for multiple times. A Coulomb friction model was used for the diaphragm-diaphragm and diaphragm-fabric tangential interactions with a coefficient of 0.6 and 0.52 respectively, these values were deduced from experimental results in Ref. [7].
The pressure differentials were applied directly to the surfaces representing the diaphragms, the loading paths for each surface are shown in Fig. 10. The time A-B is the evacuation of the air between the diaphragms. During time B-C the tool is lifted and between time C-D air is evacuated between the bottom diaphragm and the tool, forming the material into shape. The modelling was performed under the assumption that the process is rate independent and hence the actual loading rates, i. e. vacuum rate and punch speed, were not preserved in the model. Instead, the loading rates were adjusted to minimise inertial effects that are common in explicit finite element analysis.

Single layer test cases
The simulation results for both the 0/90 • and ±45 • case are presented along with the experimental results in Figs. 11 and 12 respectively. In both cases quite severe wrinkles are present. The model  appears to capture the wrinkles as well as the net shape of the preform. Fig. 13 shows a more quantitative comparison of the boundary profiles along with the wrinkle lengths and positions. The boundary profile is typically not overly sensitive to the constitutive model used, however, when wrinkles are present they can severely distort this. The prediction for the 0/90 • configuration compares very well with the experimental results, capturing the profile as well as the wrinkle locations and lengths. A more detailed comparison of the wrinkle sites for the 0/90 • case is presented in Fig. 14. Both the shape and size of the wrinkles agree with those present in the experiment. The ±45 • predicts the sites of the wrinkles but their magnitudes are slightly under-predicted, which subsequently leads to discrepancies between the boundary profiles. Nevertheless, the approach is able to predict the presence of wrinkles and their locations, capturing the effect that ply orientation has on their formations.
It should be noted that the fabric in both the experiment and simulation is drawn on to the tool through the application of atmospheric pressure. This gives further confidence to the constitutive model as the predicted geometry is achieved with the same energy that is applied in the experiment, which suggests that the mechanical behaviour is captured with sufficient accuracy.

Multi-layer test case
To simulate the simultaneous forming of multiple layers, the individual layers were modelled separately and their interactions captured using the general contact formulation available in Abaqus. Based on the experimental findings of [22], where the tangential inter-ply interactions between a ±45 • and 0/90 • plain carbon fibre weave was studied, the fabric-fabric friction coefficient was set to 0.28. The relative orientations of each ply was set within the user material subroutine and all other parameters of the model remained the same as those used in the single ply test case.
The topology of the final deformed surfaces is compared with the experiment in Fig. 15 and the boundary profile and wrinkle sites are overlayed in Fig. 16. The simultaneous forming of multiple layers, at different orientations, has a significant effect on wrinkle formation, this is evident in both the experiment and simulation. Each ply tries to deform independently to accommodate the change in shape, forming wrinkles in the plies due to their own intra-ply behaviour but are forced to accommodate the wrinkles of the adjacent ply. The model captures this behaviour well, with the wrinkle sites and shapes being accurately predicted.

Influence of inter-ply friction
Friction plays a key role during the forming process as it affects the in-plane constraint on the material. For multi-layer forming processes, where different plies are included, this can have significant effects as it constrains the relative sliding of each ply to the deformations  A.J. Thompson et al. mechanisms of the adjacent contacting plies. A second model of the simultaneous forming of the [±45 0/90] lay-up was performed again but with no inter-ply friction (μ = 0).
In the double diaphragm process, the evacuation of the air between the two diaphragms, prior to forming, compresses the fabric, therefore, large inter-ply friction forces are likely to be present, this is evident from the comparison of the in-plane shear angles presented in Fig. 17. For the case where inter-ply friction was considered, the distribution of in-plane shear angles is severely affected when compared to the no friction case. The no friction case appears to allow the layers to deform almost independently in-plane, the wrinkle severity is also reduced however wrinkles are still present (see Fig. 18). This appears to be due to the different bending behaviour of each ply in relation to the tool which modifies the normal forces early in the process.

Loading procedure
The vacuum applied draws the material onto the tool after the tool has been lifted. This results in significant bending deformation of the fabric early in the process which matures into severe wrinkles once the fabric is drawn onto the tool. The possibility of applying the vacuum to the forming box, prior to protrusion of the tool is examined here.
The comparison of the two loading scenarios is shown in Fig. 19. The severity of the wrinkles are reduced by applying the vacuum prior to the protrusion of the tool. The application of the vacuum, earlier in the process, applies more tension to the diaphragms during forming which minimises bending early in the process and subsequently, the pinching effect observed in section 2.1 is reduced. This tension also transfers to tensioning of the fabric as the pressure draws the fabric towards the tool as it is lifted. This produces a similar effect to that observed when applying increased blank holder pressure in blank holder forming processes [37].

Conclusion
Experimental investigations and accompanying numerical analyses were performed to examine the formation of out-of-plane defects during the double diaphragm forming of a plain weave fabric. The large bending strains associated with DDF, early in the process, produce an undulated surface that, when drawn on to the tool under hydrostatic pressure, results in a pinching phenomenon producing severe wrinkles. The susceptibility of the DDF process to such phenomena is magnified by the use of atmospheric pressure to form the material to the tool. As the forming force is applied normal to the surface, any undulations become trapped and evolve into wrinkles. Layers with different fibre orientations relative to the tool result in their own wrinkle pattern. When formed together the wrinkles from each layer are superposed to create an overall greater level of wrinkling.
A mutually constrained shell-membrane modelling approach has been introduced to simulate this phenomenon. This captures the anisotropy of the low bending stiffness and high tensile stiffness as well as the non-linear shear stiffness of textiles. Applied to the double diaphragm forming of a plain woven fabric over a tetrahedron tool, the model has shown to accurately predict the wrinkle patterns present in different ply orientations as well as in the simultaneous forming of multiple plies. Capturing both the severity of the defects as well as the  boundary profile of the deformed fabric. Through a further parametric study it was found that reduced inter-ply friction and increased in-plane tension reduce wrinkle severity, similar to what has been observed in blank holder forming processes, but this does not completely remove the defects.
This study has shown that by including the identified requirements for state-of-the-art forming simulations, the DDF process can be accurately modelled and severe defects, common to the process, can be predicted. Nevertheless, future work on this subject would benefit from an examination of some of the necessary assumptions made in this work. More precisely, the assumption that the compressive stiffness of the material is negligible and that the compaction of the material does not affect other deformation modes. Although these are common assumptions to make for forming simulations, new approaches recently introduced [38] to include the compressive behaviour of textiles into forming simulations provide a promising means to examine the implications of such assumptions.

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