Modelling the Shear-Tension Coupling of Woven Engineering Fabrics

An approach to incorporate the coupling between the shear compliance and in-plane tension of woven engineering fabrics, in finite-element-based numerical simulations, is described. The method involves the use of multiple input curves that are selectively fed into a hypoelastic constitutive model that has been developed previously for engineering fabrics. The selection process is controlled by the current value of the in-plane strain along the two fibre directions using a simple algorithm. Model parameters are determined from actual experimental data, measured using the Biaxial Bias Extension test. An iterative process involving finite element simulations of the experimental test is used to normalise the test data for use in the code. Finally, the effectiveness of the method is evaluated and shown to provide qualitatively good predictions.


Introduction
Press forming of woven engineering fabrics can be used to create complex geometries, suitable for subsequent liquid moulding and cure for the manufacture of composite parts [1].During the press-forming process, in-plane tension is generally used to mitigate process-induced defects such as wrinkling and to some degree to control the final fibre orientation distribution across the component after forming [2][3][4].Tension is controlled through boundary conditions applied to the perimeter of the material using a blank-holder [3][4][5][6][7].The deformation kinematics of woven engineering fabrics during the forming process is dominated by trellis shear.However, the tension along tows that occurs as a result of the blank-holder load applied around the perimeter of the forming blank and due to the forming process itself can influence the shearing resistance of the woven fabric [2,8,9].As such, consideration of the shear-tension coupling, when formulating constitutive models, could possibly result in improved accuracy both in terms of shear angle and wrinkling predictions.With the exception of Lee et al. [8,9], current constitutive models for woven engineering fabrics assume no coupling between the shear resistance and the tension in the fabric despite strong experimental evidence showing that such a coupling does exist, for example, [10][11][12][13].This paper describes a method of introducing a shear-tension coupling into finite element (FE) simulation predictions [14].Experimental data measured recently using a novel shear test for woven engineering fabrics, the Biaxial Bias Extension (BBE) test, [10] is used to fit model parameters and the predictions of the model are then evaluated using two simple numerical tests.The structure of the remainder of this paper is as follows.A brief description of the FE model used in the fitting process is given, the method of implementation of the shear tension coupling in the constitutive model is described, and the iterative procedure used to fit the model to experimental results is discussed.Finally, predictions of the model are compared against experimental shear force measurements produced using the BBE test.

Finite Element Modeling Strategy
The commercial FE code Abaqus Explicit has been used throughout this investigation.The FE model uses the same combination of mutually constrained truss elements (representing the high tensile stiffness fibres) and membrane elements (representing the shear properties of the fabric) as that described in [15] (see Figure 1).

Linear truss elements
Figure 1: FE unit cell representation of textile structure modelled using mutually constrained membrane and truss elements (from Harrison et al. [16]).
The mesh is automatically generated using an in-house mesh generation code.A simple approximate homogenisation method has been used to calculate truss dimensions and mechanical properties.Using the equation where  1 is the cross-sectional area per unit length of the ends of either the warp tows of a typical glass fabric (e.g., ∼0.000086 m 2 per metre in Harrison et al. [10]) and  2 is the combined cross sectional area per unit length of the truss elements in the mesh,  1 is the tensile stiffness of typical glass tows (e.g., 30-73 GPa [16][17][18][19][20][21]) and  2 is the stiffness of the truss elements used in the FE mesh.The truss properties chosen for the truss elements here (stiffness = 6 GPa, length = 0.01237 m, circular crosssectional area 1 × 10 −6 m 2 gives an area per unit length,  2 , of 0.000082 m 2 per m) produce a sheet with a tensile response between about 5 and 13 times lower than an actual woven glass fabric, and for simplicity, the nonlinear tensile behaviour in the tows due to fabric crimp, for example, [22][23][24], is neglected.In this investigation, decreasing the tensile modulus of the truss elements in this way has been found to produce improved performance when modelling a sheartension coupling and also tends to reduce simulation times when using the explicit FE method (due to the Courant stability condition).Previous researchers have also used this technique to improve computational efficiency [24,25].If this is done, care has to be taken to ensure that this reduction in stiffness has a negligible influence on the final complex forming simulation predictions.For example, in one forming case study, Willems [24] found that reducing the tensile stiffness by factor of 20 caused a 2 ∘ of change in the resulting shear deformation predictions.In this investigation, as will be shown, the method of implementing the sheartension coupling is based on the tensile strain along the fibre directions.The latter influences the coupling behaviour and is determined by both the truss properties and the mesh density.Thus, once the shear-tension coupling is calibrated to a given mesh density, any subsequent change in mesh density has to be compensated for by an appropriate change in the truss properties (either by changing the modulus or cross section of the truss elements).Given that the main aim of this work is to examine the possibility of modelling the experimentally observed coupling between in-plane tension and shear stiffness [10], more accurate modelling of the tensile response of the fabric is deferred to future work.Ideally, this will involve correctly capturing the coupling between tensile strains in the two fibre directions due to fabric crimp.
The membrane elements provide no contribution to the tensile stiffness of the mesh and are only used to add shear resistance to the sheet.The membrane elements have an initial thickness of 0.0002 m with a Poisons ratio of 0. The shear stresses within the membrane elements are modelled using an enhanced version of the shear part of the original non-orthogonal constitutive model [15,26,27] (S-NOCM), as discussed in the following section.By replacing the tensile part of the original non-orthogonal constitutive model (T-NOCM) [15,26,27] with truss elements, the stress field within the membrane elements can be completely decoupled from the tensile stresses occurring along the fibre directions within the membrane element.The shear stress in the membrane elements can consequently be precisely controlled as a function of any of the state-dependent variables defined within the user-subroutine used to implement the constitutive model (e.g., shear angle, angular shear rate, temperature, or strain along the fibre directions).This strategy has been used recently to create a rate-dependent or viscous constitutive model for thermoplastic advanced composites [16,28,29].The original implementation of the S-NOCM VUMAT usersubroutine has been modified in order to implement a sheartension coupled version of the model, as described in the next section.

Implementation of Shear-Tension Coupling in the S-NOCM
Implementation of the shear-tension coupled S-NOCM involves linking the shear parameters in the original S-NOCM model with the tensile stresses (or equivalently the tensile strains) acting along the warp and weft fibre directions in the fabric.Like the shear angle, the tensile strains are accessible as state-dependent variables within the Abaqus user-subroutine.In this section, a method of producing the same shear-tension coupling in the numerical model as that measured in actual woven engineering fabrics is described.
The technique involves a four-stage process, as follows.

Stage
One.This involves simulating the BBE test; details of the actual experiments can be found in [10].A BBE test sample with dimensions 210 × 210 mm and a clamping length of 70 mm is modelled (see Figure 2) using mutually constrained truss and membrane structural elements (572 truss and 264 membrane elements) as shown in Figure 1.
The typical computation time for each simulation was about 10 minutes using a Dell OptiPlex 760 Intel (R) Core(TM)2 Duo CPU E7500@2.93 GHz and 3.25 GB of RAM running Abaqus Explicit v 6.9.Faster simulation speeds could have been obtained using the symmetries of the test (e.g., by simulating just a quarter of the test), though this would   require modification to the automatic mesh generator to create triangular elements along the centrelines of the specimen, also since future work will involve exploring the influence of fibre orientation variability on test results (e.g., see [30]), this would negate the existing symmetries of the test, and so full specimen simulations have been conducted.These have been conducted in two steps.Step one involves application of a constant transverse load,    equal to the loads used in [10] (5, 37, 50, 75, and 100 N).The superscript  is the experiment number ( = 1 to 5) with each experiment using a different transverse load ( = 1 corresponds to 5 N,  = 2 corresponds to 37 N, etc).The transverse load is applied to nodes at the edge of the central section of the right and left sides of the blank (Region C in [10] see Figure 2).Step two involves applying a displacement controlled boundary condition on the upper and lower centrally located nodesets at the middle of the top and bottom side lengths of the blank (corresponding to the edge of Region C in [10]) see Figure 2. The corresponding experimental shear forces versus shear angle curves,    (), measured on a plain weave glass engineering fabric were used as input curves in the standard S-NOCM to conduct these preliminary simulations, and for simplicity, , the shear angle in Region A, is taken from one of the central elements of Region A (see Figure 2).This approximation assumes the shear angle across all elements in Region A is uniform.In practice, a variation in the shear angle within each of the regions A, B, and C exists.It will be shown later that the size of this variation is small and depends on the shear angle and the size of the transverse load applied to the specimens.   () are initially approximated from the axial load,   (), [10] using (2).In Stage 4 of the fitting process, this estimate is improved using a simple normalisation procedure Note that to determine   , contributions to the measured total axial force,   , from the reaction force,   , which is caused by application of the transverse clamping load,   , must first be removed before applying (2).The method of doing this for experimental results is described in [10].To do this for the numerical results, the following equations are used: where   is determined from   using where V  and V  are the vertical and horizontal velocities of the nodes at the upper, bottom, right, and left node sets, respectively.

Stage
Two.This involves determining the average tensile strains, , along the warp and weft fibre directions,   and  weft , as a function of the shear angle for  = 1 to 5. The tensile strains are given as state-dependent variables within the VUMAT user-subroutine and have been verified to be the same as the tensile strains occurring along the truss elements Advances in Materials Science and Engineering bounding the corresponding membrane element.The average tensile strain across the entire specimen along the two fibre directions is determined as a function of the shear angle in Region A, by taking an average of  warp and  weft from a selection of elements across both Regions A and B. The average fibre tensile strain is determined for each value of the transverse loads,    , as a function of the shear angle, and a polynomial curve is fitted to the data from each of the five simulations,    (), the coefficients of which are stored for later reference by the enhanced S-NOCM code during the course of the simulations (the p subscript indicates this is a fitted polynomial function).Thus, each shear force input curve    () has a corresponding average fibre strain curve    ().

Stage Three.
This involves implementing the sheartension coupling in the VUMAT user-subroutine.To do this, code has been added within the original VUMAT usersubroutine for the S-NOCM to compare the value of  in each membrane element at each time increment against the values of    () using the shear angle within the element (also given as a state dependent variable in the VUMAT user subroutine).Depending on the value of , the code assigns the appropriate shear force curve    () to the element using the algorithm given in flow chart Figure 3.The shear stress within the element is then determined using the S-NOCM.
Thus, the shear force input curve is now a function of both the shear angle and the fibre strain within the membrane element.The process is illustrated in Figure 4, which shows actual shear force data measured in experiments and values of the average tensile strain along the fibre directions predicted in the FE simulations of the BBE test (see Figure 4).The process of assigning the appropriate shear force versus shear angle curve is described and illustrated in Figure 4 using a specific example.Note that in Figure 4, only data corresponding to tranvserse loads of 5, 50, and 100 N are shown in order to simplify the figure.
Consider an element that has a shear angle of 45 ∘ at time, .The average tensile strain, , inside the element is determined, and in this case, the value is 0.03.An orange point indicates the (, ) coordinate in Figure 4.The algorithm in the flow chart shown in Figure 3 is run to determine where the average tensile strain in the element, , lies in relation to the average tensile strain versus shear angle polynomial curves,    () (plotted as black lines in Figure 4).Once the appropriate polynomial is identified and assigned to the element (the assignment is indicated by a blue arrow in Figure 4, in this case,  = 3 for the 50 N transverse load), the corresponding shear force versus shear angle curve    () (plotted as red lines in Figure 4) is also assigned to the element, indicated by a red arrow in Figure 4.    () is used to determine the shear stiffness of the membrane element using the S-NOCM, as has previously been described in detail in [15].
At this point, it is possible to compare the results of the enhanced S-NOCM against the experimental input data, as shown in Figure 5. Here, experimental data from [10] are plotted as thin continuous lines with error bars (a different colour for each transverse load), and numerical predictions are plotted as thick continuous lines (the same colour as the corresponding experimental curve).Agreement between numerical prediction and experimental input curve is quite poor at this stage as the experimental shear force input curves supplied to the code are not yet normalised.The issue of normalisation of shear test results for advanced composites in bias-extension tests is well known, and some author have used gauge sections [31], while others have considered the energy or force contributions from the entire specimen, including Regions A, B, and C [32][33][34][35][36]. In these cases, the precise normalisation procedure depends on the test method (uniaxial or biaxial), the specimen geometry, and the material's response during shear (rate dependent, rate independent, or showing a coupling between tensile strain and shear resistance).A theoretical method to normalise BBE test results for materials with a strong shear-tension coupling was described in detail in [37].The method requires custom software to retrieve the underlying normalised data via an automated iterative process.Future work will involve use of this theory for accurate and fast normalisation.For now, a simpler approximate normalisation technique is described in the final stage, Stage 4, of the fitting process.The aim of normalisation procedures is to find the shear response of the fabric per unit length or per unit area, in order to determine the parameters governing the shear behaviour in the material's constitutive model.Normalisation is relatively simple for the picture frame test [33], where the entire specimen undergoes homogenous deformation but is more complex for bias-extension tests.Here, the specimen undergoes different deformations in different regions (e.g., see Figure 2), and this has to be taken into account when interpreting test results.

Stage
Four.This involves a simple normalisation procedure aimed at normalising the experimental input curves (which have to be supplied as shear force per unit length of fabric).By correctly normalising the experimental biaxial bias-extension curves, the numerical simulations should produce approximately the same shear force versus shear angle predictions as those observed in experiments.To do this, an approximate procedure is used here by the following simple iterative method.(i) The input shear force versus shear angle curves are divided by the predicted shear force versus shear angle curves to produce a ratio (also a function of the shear angle).(ii) Polynomial functions,    (), are fitted to each ratio curve (iii) Input curves are multiplied by the ratio curves to produce a next generation of input curves (iv) The process is repeated until reasonable agreement between numerical BBE test predictions and experimental results is obtained.Normally around three iterations are required.This is a simple method designed only to examine the possibility of introducing a shear-tension coupling in the model.Future work will involve employing the more rigorous normalisation developed in Harrison [37].Figure 6 shows the comparison between the original experimental results and the final predicted shear force versus shear angle curves after conducting this normalisation process.The horizontal error bars given on the numerical results indicate the variation in shear angle across Region A, calculated using the standard deviation of the shear angle of all elements in Region A. The vertical error bars on the experimental results indicate the variation in the measured force, calculated using the standard deviation of 3 tests.Thus, the full length of each error bar represents two standard deviations.The agreement between numerical predictions and experimental data is clearly improved compared to Figure 5.
To test the effectiveness of the modelling approach, two final BBE simulations are conducted, this time using transverse loads increasing linearly in time from 5 N to 100 N rather than using constant transverse loads.In Figures 7(a) and 7(b), the grey curves are experimental results originally reported in [10], and the black curves are the numerical predictions following the approximate normalisation process described in Stage 4, when applying constant transverse loads of 5, 37, 50, 75, and 100 N (the same information is shown in Figure 6).The blue curves in Figures 7(a) and 7(b) are the results predicted by the coupled S-NOCM when increasing transverse loads are applied over the course of the test.In Figures 7(c) and 7(d), the applied transverse load is plotted against  rather than against time, creating slightly nonlinear profiles.In Figures 7(a transverse load is held constant at 5 N for the first 60% percent of the total simulation time, then increased linearly to 100 N for a further 80% percent of the total simulation time, and then held constant at 100 N until the end of the simulation.As expected, the axial force predictions of the enhanced shear-tension coupled S-NOCM, made using increasing transverse loads, move across the normalised numerical predictions generated using constant transverse loads (the black curves).The different transverse loads versus shear angle profiles,   (), shown in Figures 7(c) and 7(d), produce different axial force predictions, as can be seen by comparing Figures 7(a) and 7(b).The result in Figure 7(a) is close to that which might be expected from the woven glass fabric used in the experimental investigation [10].However, while the result of Figure 7(b) appears correct until around 30 ∘ , an unrealistic softening is apparent above this shear angle.Thus, at this point the predictions of the model have been found to be qualitatively correct under simple loading conditions though they can show unexpected behaviour under more complex loading.Possible explanations for the unexpected predictions could be related to the following.
(i) The first is the choice of elements used to create the average strain curves,    ().The resulting predictions have been found to be sensitive to this choice, and future work may involve using a more refined mesh to model the BBE test and use a larger selection of elements to examine this sensitivity.
(ii) The second is the normalisation technique used in this work.The very simple normalisation procedure used here takes no account of the shear-tension coupling in the fabric, and a more rigorous method was recently proposed in [37].Future work will aim to employ this method to improve accuracy and reduce the uncertainty in the shape of the input curves passed to the S-NOCM.
(iii) The third is the method of calculating the stress increment at each time step.A tangent stiffness matrix has been used to determine this stress increment; that is, The linearisation process is known to reduce the sensitivity of the technique of using multiple input curves to control the shear compliance of the membrane elements, a point discussed in detail in [16].Nevertheless, the linearised increment was used in this first attempt to model to the shear-tension coupling, as the method has the advantage of being particularly robust.Future work will focus on improving the sensitivity of the approach, using the methods described in [16].
Despite the irregularities in the predictions of the sheartension coupled model under certain in-plane loading conditions, it is clear that the technique proposed here produces a shear-tension coupling similar to that seen in actual experiments.Future work will focus on improving the accuracy of the method, though the model predictions are considered to be sufficiently accurate at this stage to begin to examine the question of whether or not and also under which conditions, the influence of a shear-tension coupling on the shear angle and wrinkling predictions of complex forming simulations, is important.

Conclusion
A method of modelling the coupling between shear compliance and in-plane tension in woven engineering fabrics has been demonstrated.The method is similar to that used previously to create rate-dependent "viscous" behaviour using a hypoelastic model [16] though here the average in-plane strain along the two tow directions, rather than the angular shear rate, is used to control the selection of the shear force versus shear angle curve for use in the non-orthogonal constitutive model (used to relate the shear force and shear stress) [8,9].A simple normalisation procedure has been proposed.The sensitivity of the modelling approach is assessed and found to give reasonable results, clearly showing a coupling between shear compliance and in-plane strains in the fibre directions.Future work will involve refining the modelling and normalisation process in order to improve the accuracy of the predictions and could also involve reimplementing the technique using fibre stress rather than strain to control input shear curve selection.The shear-tension coupled model will be used to evaluate the importance of a shear-tension coupling on the predictions of complex forming simulations.

Figure 2 :
Figure 2: The BBE FE model.Force boundary conditions are applied to the right and left centrally located node sets, and vertical displacement boundary conditions are applied to the upper and lower centrally located node sets.The colour legend indicates the shear angle.The three different deformations occurring in Regions A, B and, C of the test specimen are clearly visible.The shear angle in Region A is taken from the highlighted element.

5 𝑠Figure 3 :Figure 4 :Figure 5 :
Figure 3: The flow chart of shear-tension coupling algorithm which runs for each membrane element at every time increment during a simulation.

Figure 6 :
Figure 6: Comparison between the experimental and the predicted results using the coupled S-NOCM and normalized   −  input curves from the BBE 3 : 1 simulations.