Mathematical modelling of internal geometry and deformability of woven performs

: The paper presents an approach to model the behaviour of a representative volume element (unit cell) of textile reinforcement in in-plane deformation (bi-axial tension and shear) and in compression. The model is a further development of a virtual textile concept implemented in the WiseTex software, and is based on the concept of hierarchical description of textile properties and systematic application of the principle of minimum energy to calculate the textile geometry in the relaxed and deformed state. With the internal geometry of the unit cell built, the model computes overall parameters of the deformed textile, such as fibre volume fraction, porosity etc. The internal geometry is visualised and such properties as pore structure in typical cross-sections are analysed. The load-deformation curves for compression, tension and shear are computed via the balance between change of the internal energy of the unit cell and mechanical work of the applied loads. The internal geometry description is further fed into flow modelling software, which allows computing local permeability of the deformed reinforcement, and micro-mechanical modelling to calculate homogenised local stiffness of the composite .


Introduction
Deformability of textile preforms plays a key role in quality of a composite part formed into 3D shape and processed using Liquid Transfer Moulding (LTM) technique. Ill-chosen placement of the preform, disregarding its behaviour under large strain in complex deformation may result in the preform wrinkling or even damage, deteriorating the performance of the composite part. This explains an importance of predictive modelling of deformability of textile preforms. The deformation modes of primary importance are in-plain deformation (tension and shear) and compression of the preform. Deformability of woven preforms in these modes is the subject of the present paper. Out-of-plane bending may also be considered, as it can affect internal geometry of the preform, especially for small bending radii; a model of the woven fabric bending can be found in (Lomov et al., 2000). Naturally the deformability of woven fabrics is important as well for apparel textiles, and has attracted attention of textile materials researchers. Works of Kawabata, Niva and Kawai (Kawabata et al., 1973;Kawabata et al., 1973;Kawabata et al., 1973), de Jong and Postle (de Jong et al., 1978), Hearle and Shanahan (Hearle et al., 1978) have established an approach to mathematical modelling of deformation of woven fabrics, which can be summarised by three principles. First, the model uses deformations, rather then loads, as input for in-plane deformation (for compression model the applied pressure is the input). An overall deformation pattern is imposed over the woven fabric repeat (unit cell) to change the spacing of the yarns in tension and the angle between them in shear. As formulated by Komori and Ito (Komori et al., 1991), the unit cell is subject to transformation of coordinates defined by the given deformation. The contacts between yarns stay unchanged in tension, and experience rotation (not sliding) of the contacting yarns in shear.
Second, the principle of minimum energy is applied to compute the internal geometry of the deformed fabric. With the spacing and orientation of the warp and weft given, the yarn paths are defined using one of available geometrical models (Peirce's, elastica, splines), with crimp heights and dimensions of the cross-sections of the yarns (which can change under transversal force caused by tension) as parameters. These parameters are calculated via the principle of minimum total energy, associated with the yarns tension, bending and compression. Tension of the yarns is computed using the experimental tension diagram. Elongation of the yarn is estimated by the difference between yarn length in the repeat after and before the deformation. Experimental bending and compression diagrams are used to compute bending energy and resistance to compression. All these experimental diagrams are non-linear. Application of the principle of minimum energy to fibrous assemblies must be considered as heuristics, as it can be applied rigorously to conservative systems only. This means that frictional effects are not taken into account in the solution of the minimisation problem. The internal friction between fibres in the yarns enters the calculations via non-linear bending and compression diagrams. The inter-yarn friction is absent in tension, and manifests itself in the rotation of the contacting yarns in shear.
Third, the applied loads are computed via the balance between, on one hand, the mechanical work done by the loads on deformations of the unit cell, and, on the other hand, the sum of the change of the total energy of the deformed yarns and the work of friction (if any).
We will use this scheme for all three types of deformation under consideration: compression, bi-axial tension and shear.
Returning to deformability of woven reinforcements for composites, one found quite a number of publications on modelling. The models complexity range from simple empirical models to elaborate finite element descriptions. Certain important points have been investigated, which were not covered by earlier "apparel-oriented" models.
In studies of compression the attention was given to compressibility of the reinforcement at high loads, which are characteristic for composite processing. The compression curve is broken into three regions (low, medium and high loads), with the different phenomena playing in each of them (Chen et al., 1999a;Chen et al., 2000;Chen et al., 2001b;Chen et al., 2001c). The nesting of layers of the reinforcement is taken into account (Chen and Chou 2000;Lomov et al., 2000;Kurashiki et al., 2002). Models of shear of woven reinforcements (Long 2000;Long et al., 2001a;Long et al., 2001b;Crookston et al., 2002) have to deal with very high shear angles (up to 60-70°) occurring in forming of complex 3D parts. This is dealt with by an introduction of models of lateral compression of the yarns, which come into contact when the shear angle reaches and exceeds the locking angle of the fabric "trellis". The simple, but not true-to-life concept of preserving the volume of the unit cell to calculate the change of the thickness of fabric in shear, has been advanced to more correct considerations of the yarns compression.
Recently finite element descriptions of deformability of textiles have been introduced for bi-axial tension (Boisse et al., 1999a;Boisse et al., 1999b;Launay et al., 1999;Gasser et al., 2000;Boisse et al., 2001a;Boisse et al., 2001b;Kuwazuru et al., 2001;Sakakibara et al., 2001) and compression (Kurashiki, Zako et al. 2002) of woven reinforcements. Based on the increasing power of the computers, these approaches aim to describe in detail the 3D behaviour of the fabric constituents including contact and friction and to obtain result fields at local level.
The finite element modelling encounters two difficulties. First, the geometry of a textile unit cell is very complex, and creating a solid model manually is not an easy task. The solution is provided by the use of textile geometry modeller as a preprocessor, capable to create a finite element model automatically (Kondratiev 2001;Van Genechten 2002) Second, the description of the material behaviour used in finite element model must realistically represent the actual behaviour of the fibrous assemblies -yarns. Development of such a library of material model for textiles, available in finite element packages, presents a serious challenge to researchers.
Models of textile deformability developed so far have built a solid foundation for their generalisation, encapsulating the achievements of textile material science in a modelling software tool, allowing wide variability of textile structure and yarn parameters, instrumented with visualisation features and able to transfer the models of textile geometry into specialised micro-mechanical and flow modelling software as well as into general purpose finite element packages. Such a tool can be considered as a preprocessor for calculation of homogenised propertiespermeability tensor and stiffness matrix -of deformed textile reinforcement. These properties, in their turn, are used as input to provide local parameters in modelling of Darcy flow through the deformed preform and structural finite element analysis of a 3D shaped composite part. This work is in progress in Composite Materials Group in the Department MTM, K.U.Leuven (Lomov et al., 2000;Lomov and Verpoest 2000;Lomov et al., 2001a;Lomov et al., 2001b;Belov et al., submitted). It has resulted in the development of textile modelling software WiseTex, available via the Department MTM. The present paper describes models of deformability of woven fabrics, implemented in it. Early versions of the models described here, were published in (Lomov et al., 1992;Lomov et al., 1995a;Lomov et al., 1995b;Lomov et al., 1997).

Description of the geometry of woven fabric in relaxed state
The comprehensive description of the model of the relaxed state of a woven fabric can be found elsewhere (Lomov et al., 2000;Lomov et al., 2001a;Lomov et al., 2001b). Here we state its main components used in simulation of the fabric deformation.

Weave pattern and elementary crimp intervals
A weave pattern (for one-and multilayered fabrics) is coded with matrix coding (Lomov, Gusakov et al. 2000;Lomov and Verpoest 2000;Lomov, Huysmans et al. 2001;Lomov, Huysmans et al. 2001). It allows separation of the crimped shape of the warp and weft yarns into elementary bent intervals (Figure 1), representing sections of the yarn between interlacing sites. The shape of the yarn on an elementary interval is described using a parameterised function z(x;h/p), where z and x are coordinates of the yarn middle line, h is the crimp height and p is the distance between the interval ends (spacing of the yarns). The shape z(x;h/p) for a given relative crimp height h/p is computed using the principle of minimum of bending energy of the yarn on the interval and has a form where the first term is a spline function, corresponding to the solution of the linearised minimum energy problem, and the second term represents a correction for a non-linear formulation. The function A(h/p) is calculated from the solution of the minimum energy problem and is tabulated. With this function known, the characteristic function F of the crimp interval is computed, representing the bending energy of the yarn: where B(κ) is the (measured experimentally) bending rigidity of the yarn, which depends (non-lineary) on the local curvature κ(x), or, after the integration, on an average curvature over the interval Function F(h/p) is also tabulated. With the function F known, the transversal forces acting on the interval ends can be estimated as

Compression of the yarns in the relaxed fabric
Warp and weft yarns in the relaxed fabric are compressed by the transversal forces Q [4] according to experimental diagrams, measured on "virgin" yarns where subscript "0" refer to the uncompressed state of the yarn, d 1 and d 2 are dimensions of the yarn cross-section ( Figure 1). These dimensions and crimp heights of the yarns are interconnected: where superscripts refer to the warp and weft yarns, subscripts "1" and "2" refer to two weft yarns in different layers, ΔZ is the distance between fabric layers ( Figure  1).
With crimp heights of weft yarns given, equations [4-6] provide a closed system of non-linear equations for calculation of the transversal forces Q and yarn dimensions d 1 and d 2 .

Minimum energy problem -calculation of the weft crimp heights
The weft crimp heights are found using the principle of minimum bending energy of the yarns inside the unit cell. Using [2], it is written as where subscripts i,j refer to different warp and weft yarns, k -to the elementary crimp interval of the warp/weft yarn. The minimum problem [7] is solved for the weft crimp heights, all other parameters defined inside the minimisation algorithm via solution of the system [4-6] for the given current crimp heights. It takes about 15s on Pentium II (405 MHz) PC to compute parameters for a 3D fabric with 20 yarns in the repeat, and about 0.5 s -for a plain weave fabric.

Compression
A comprehensive description of the compression model of a woven fabric can be found elsewhere (Lomov and Verpoest 2000). Here we give an outline of the algorithm and discuss two issues not covered in that paper: behaviour of z-yarns in thick 3D fabrics and nesting of the layers in compression of laminate.

Outline of the algorithm
When a fabric is compressed, the following changes in geometry take place: -Warp and weft yarns are compressed; -The less crimped yarn system increases its crimp and vice versa.
These two processes are treated in the model separately. This follows from the assumption of an even distribution of the compressive force over warp/weft intersections, because this assumption implies that force per one intersection, which compresses the yarn cross-sections and bends the yarns, is independent of any changes of warp and weft crimp or cross-section dimensions.

Compression of yarns
To compute this, the compression force per one intersection is evaluated: where F is the pressure force on fabric repeat. This value is added to all the Q ijtransversal forces acting on the intersections and computed with [4], to evaluate the dimensions of the yarns with [5]. Hence, both the compression due to yarn bendind and the compression due to external force are accounted for. The algorithm presented above is then applied to yield the compressed dimensions of the yarn cross-sections and the new values of the yarn crimp.

Change of crimp
The change of crimp in compression (increasing for warp and decreasing for weft, or vice versa) leads to a decrease of the fabric thickness. Therefore the basic mechanical equation governing this process is work of compressive force Q on change of thickness db = [9] = change of bending energy of yarns dW The compression of yarns has been accounted for before and the resulting changed cross-section dimensions are "frozen", that's why the work of yarn compression does not enter the balance [9].
Changes of the fabric thickness db and bending energy of the yarns dW depend on the change of the set of weft crimp heights {dh j We } and therefore [9] has a set of unknown variables. A reasonable assumption to cope with this difficulty is: "The crimp changes in such a way as to provide the maximum possible change of thickness". This means that if we consider the function b({h j We }) (b being the fabric thickness) then changes of crimp will follow the direction of the maximum slope (in the opposite direction): where x is computed to satisfy [8]; slashed values refer to changed crimp. Equation [9] is then written as follows: to be solved numerically for x, which should also satisfy

Behaviour of the z-yarns in 3D fabric
The algorithm described above does not preserve the length of the yarns after compression, hence introducing an error on the final fabric geometry. In the compaction of a 2D fabric this error is small and can be considered to be negligible. When a 3D fabric is compacted, its z-yarns (going through the thickness), deviate considerably from their paths, as the length of the yarn must be preserved when the thickness of the fabric is reduced. In a 3D fabric with z-yarns initially almost vertical, they will acquire S/Z shapes. If the interlacing yarns are oblique, their initial sinusoidal shape is transformed into meander.

measured (lines) and computed (dots). Inlet: Computed and observed shape of z-yarn at 300 kPa
To describe such behaviour rigorously, one would have to account for all kinds of contacts between yarns occurring during the compaction. This task may be an interesting challenge for the finite element modelling. In the present model, implemented in WiseTex, a simple geometrical approach is followed. A spline correction is added to the yarn path, with a factor chosen as to preserve the yarn length after the deformation. The result is illustrated in Figure 2. Figure 2 shows the results of measurement and simulation of compression of 3D glass fabric (yarns 4x4x33 tex, 35 yarns/cm in warp and weft, areal density 3900 g/m 2 , weave is shown in Figure 2). The fabric is proposed in (Parnas et al., 1995) as a benchmarking case for study of LTM composite processing. Compression of the fabric has been measured on the KES-F textile compression tester for the low load range and on Instron for higher loads. The input data on compressibility of the yarns and their bending rigidity has been measured on KES-F. Other fabric data was taken as specified in (Parnas, Howard et al. 1995). We see that the described algorithm provides a reasonable prediction of the compression diagram as well as the fabric internal structure after compression.

Nesting in a laminate
In composites forming processes normally stacked 2D fabric layers are compressed to form a plate or a 3D shaped laminate. It is well known that due to nesting of the fabric layers thickness of the laminate per one layer is decreasing with the number of layers (for a given pressure) (Luo et al., 1999;Chen and Chou 2000;Lomov and Verpoest 2000) (there exists an evidence of an anomalous opposite behaviour (Pearce et al., 1995)). The presented model of the fabric geometry and compression accounts for this fact.  Consider a plain woven glass fabric, made of glass rovings 480 tex, with the yarn density of 36/34 yarns/cm. The fabric has been thoroughly studied in (Lomov, Gusakov et al. 2000;Lomov and Verpoest 2000), where additional information can be found. In these papers a WiseTex model of the fabric geometry and compression has been verified against experimental data.
Consider a laminate made out of L layers of this fabric, compressed at a given pressure. The layers are randomly shifted relative to each other and nested. The nesting algorithm searches for the minimum vertical coordinate of a layer mid-plane when it is just in contact with the layer beneath it. Figure 3 shows the results of a Monte-Carlo calculation of the average thickness of the laminate (500 runs), at pressure 100 kbar (compression of the fabric taken into account). The results correspond very well with the experimental variation of the thickness with the number of layers, measured in compression tests.  The model of a laminate allows also studying the statistical characteristics of the nesting of the layers. Figure 4 depicts results of such calculations. The distribution of fibre volume fraction becomes more narrow and smooth with increase of the number of layers.
It is interesting to note the correspondence between these simulations and the variability of the permeability of the laminate studied in (Hoes et al., 2001). The distribution of permeability of the laminate has been found to have a width of one binary order of magnitude (3 times, from 90 to 280 μm 2 ). On the experimental dependency of the permeability on fibre volume fraction (Parnas et al., 1997;Belov et al., 2002) the change of the permeability by 3 times corresponds to the change of fibre volume fraction by 5%, which is exactly the width of the distributions shown in Figure 4.

Bi-and uniaxial tension algorithms
Consider first a woven fabric under bi-axial tension characterised by deformations in warp (x-axis) and weft (y-axis) directions e x = Y/Y 0 -1, e y = X/X 0 -1, where X and Y are sizes of the fabric repeat, subscript "0" designates the undeformed state. As discussed above, inside the WiseTex model the internal structure of the fabric is described based on weft crimp heights h j We and weft and warp cross-section dimensions at the intersections d ij Wa and d ji We (subscripts designate different yarns in the fabric repeat). These values change after the deformation. Tension of the yarns induces transversal forces, which compress the yarns, changing d's. The same transversal forces change the equilibrium conditions between warp and weft, which leads to a redistribution of crimp and change of crimp heights. When the mentioned values in the deformed configuration are computed, the internal geometry of the deformed fabric is built using the WiseTex algorithms as explained above. Change of length of the yarns determines their average (in the repeat) deformations, which, through the tension-deformation diagrams of the yarns allow computing tensions of the yarns. When summed up, with yarns inclinations due to the crimp accounted for, the yarn tensions are transformed into loads, caused the fabric deformations. This computational scheme has been proposed for plain weaves in 70s and also recently implemented via FEA (see the Introduction); the present implementation of it puts it into the scope of WiseTex modelling, covering a wide range of woven structures.
The key problem in the bi-axial modelling is computation of crimp heights and transversal forces in the deformed structure. Assuming that the spacing of the yarns in the fabric is changed proportionally to the change of the repeat size, we compute the x and y positions of intersections of warp and weft in the deformed structure. The configuration of the yarns in the crimp intervals between the intersections is determined by these positions and (unknown) crimp heights. Consider some values of the crimp heights. Then the WiseTex geometrical model determines positions of the ends of crimp intervals (warp/weft intersections) and bent shape of the yarns in the intervals. where Q bend is the transversal force due to the yarn bending [4], T 1,2 are the yarn tensions on two crimp intervals adjacent to the point of application of the transversal force, θ 1,2 are the angles of inclination of the yarn on these crimp intervals. We assume that the tension of the yarn can be computed based on the average deformation ε of the yarns (therefore T 1 =T 2 ): where l is the yarn length. Note that T depends on yarn length after the deformation, which in its turn depend on crimp heights and yarn dimensions.
The transversal forces compress the yarns according to an experimental law of the compression [5]. When the yarn dimensions are computed and "frozen", then crimp heights are determined using the minimum energy condition: where W bend and W tens are the bending and tension energy of the yarns. The former is computed summing up bending energies of the yarns in crimp intervals between yarn intersections [7], the latter is the sum of tension energies of all the yarns, which are computed using their (linear or non-linear) tension diagrams and yarn deformations.
Step Step 1. Set initial deformations and tensions.
Step 2. Compute dimensions of yarns and transversal forces.
Step 3. Compute length of the yarns.
Step 4. Compute deformations and tensions Step 5. Check convergence for the deformations. If not, go to Step 2. Step

Figure 6. Biaxial tension algorithm
The computations described above determine one step in the iteration process: starting from current values of the crimp heights we compute yarn lengths, yarn tensions, transversal forces, yarn compressed dimensions and then new values of the crimp heights. The full algorithm is depicted in Figure 6.
If one side of the fabric is kept free (uni-axial tension, say, along the warp), then the described algorithm has another, the outmost iteration loop, searching for X<X 0 (negative e y ) which would lead to zero loads along weft (y) direction. This allows computing Poisson coefficient for the fabric.

Experimental
Experimental bi-axial tensile properties of fibre fabrics have been tested at LMSP (UMR CNRS ESAM-ESEM Paris/Orléans). The device presented on Figure  7 is based on two deformable parallelograms (Ferron, 1992). When the system is compressed (on a classical traction-compression machine), it generates tensile deformation in a given ratio in each direction of the cross specimen located at the centre of the device. One of the lozenges has adjustable dimensions, in order to set various deformation ratio (the ratio is denoted k). The measurements of loads in both directions are made using captors, positioned close to the specimen. Strain measurements are done by optical methods or by mechanical extensometers. The optical measurement permits to check the homogeneity of the strain field within the active part of the specimen. Both methods give equivalent global results, but the optical measures give the strain field and allow checking its homogeneity (Launay et al, 2002). The bi-axial tension experiments can be done for warp/weft angles from 90° to 60°. The cross-shaped specimen is well adapted to the biaxial test of fabrics because of the very weak in-plane shear stiffness. The results of a bi-axial test on a balanced glass plain weave fabric (2.2 yarns/cm, 1220 tex) are given in comparison with calculations in Figure 8 and Figure 9. More details on the experimental procedure and experimental results on other fabrics can be found in (Boisse et al, 2001a,b).

Finite element analysis
The biaxial tensile behaviour of a fabric can also be obtained from 3D finite element analyses. The main feature is the specific mechanical behaviour of the single yarn. For the fabrics under consideration, this yarn is composed of thousand of fibres with very small sections which are very flexible and that can slide with regards to others. The behaviour of the yarn is assumed to be orthotropic. Shear modulus are very small. Young moduli in the direction perpendicular to the yarn are all very small in comparison to the modulus in the direction of the yarn. A very important point is to give at each time, the orthotropic properties in the frame based on the direction of the fibres. A hypoelastic orthotropic model is used and the rotation of the orthotropic frame (as well as the rotational objective derivative) is based on the rotation of the yarn. Numerical difficulties in the analysis are due to very weak values of some coefficients (especially shear modulus) in comparison to Young modulus in the fibre direction. In order to avoid these problems, the stabilization technique used for finite elements with reduced integration, called hourglass control, is used (Flanagan 1981).

Figure 9 Computed and measured tension force. Uniaxial tension
The crushing of the yarn is very important in bi-axial tension, because the undulation variations are directly depending of these thickness changes. Generally the models relate the thickness to the contact force for a given yarn. In finite element analysis of fabric, the local value of the young modulus in the direction perpendicular to the yarn. When the yarn is under tension, the crushing is more difficult. Consequently he transverse young modulus is assumed to be in the form. n m 3 0 33 11 E 0 , m and n are three material parameters. E ε is the transverse Young modulus of the unloaded state. It is very weak (nearly equal to zero) for the single yarns of the studied fabrics. The parameters E 0 , m, n are determine using an inverse method from the biaxial test for k=1. For this ratio crushing is maximum. The results of the finite element analysis are in a very good agreement with the experimental data (Gasser et al 2000, Boisse et al, 2001a.

Comparison with WiseTex model
The current formulation of the model does not envisage the dependency of the compression diagrams of the yarns on the applied tension. It would not be a problem to modify equations [5] accordingly, to be used in the iterations shown in Figure 6. However, this information is not normally available and standard textile KES-F compression tester does not have a tension installation. Therefore the comparison of the experimental data and results of finite element simulations can answer two questions: -Does the iterative algorithm of Figure 6, simplified vis-à-vis finite element modelling, provide results close to it (and to experiment)?
-What are the errors introduced by using compression diagrams obtained without the yarn tension?
To answer these questions, WiseTex modelling of the tension of the fabrics, described above, have been performed. The yarn tension resistance was assumed to be the sum of tension resistance of the fibres. Bending rigidity of the yarns has been measured on KES-F to provide the value of 0.5 N mm 2 . Compression law was derived from equation [13] using different constant levels of tension (ε 1 ). These diagrams are referred below as "Compression ε 1 =…%". Figure 9 show the results of the calculations (note that the ends of the curves do not represent failure). We can conclude the following:

Figure 8 and
-The calculations describe well the qualitative difference between deformation regimes of uniaxial and biaxial tension.
-For the uniaxial tension the error of calculations is small whatever compression diagram is used. The standard textile compression tester can be used to gather input data for simulation of this tension regime.
-For the biaxial tension the difference between calculations with compression laws corresponding to different tensions can be as large as 30%. A reasonable correlation is found when the compression diagram used corresponds to the highest level of strain of the yarns. For the low strain region the experimental curve corresponds better to the calculations with compression diagrams for low tension, and vice versa for higher strain.

Internal geometry of a sheared woven fabric
When a woven fabric is sheared, the orthogonal directions of warp and weft become skewed. Such a configuration is similar to a braided structure. A model for internal geometry of braids is proposed in  and is used to construct the sheared woven fabric internal geometry.
The most important feature of non-orthogonal unit cell of a sheared fabric is the presence of twist of the yarns. The twist is evident in cross-sections of nonorthogonal fabrics ( Figure 10) and can be mathematically described using the algorithms of (Lomov, Nakai et al. 2002). As shown in this paper, the twist can be considerable, attaining the values as high as 30-40°.
Consider a sheared unit cell of a woven fabric (Figure 11). For simplicity the illustrations below show a plain weave unit cell. The equations are however applied to elementary crimp intervals and the forces are summed up for the actual fabric repeat. An account is taken therefore of the differences imposed by the weave pattern, as shown in comparison with experimental data below.
Our aim is, given a value of the shear angle, γ, compute the shear force, T, in the presence of (pre)tension of the fabric. The tension is dealt with according to the algorithms of the previous section, resulting in values of the tension of yarns and transversal forces Q, associated with it [12]. The shear force T is related to a moment M at the point O and to the mechanical force A done at the shear deformation: We will take into account the following mechanisms of the yarns deformation, determining the shear resistance: Accordingly the mechanical work A, moment M and shear force T are subdivided: The lateral compression of the yarns is introduced via the transversal forces. We do not consider intra-yarn friction here, as suggested in (Harrison et al., 2002). It is felt that this factor is accounted for by the lateral c7ompression calculations, but the question needs more careful examination in future work.

Transversal forces: Lateral compression of the yarns
The transversal forces, acting on the yarns, and determining the friction between them, are caused by tension and bending of the yarns, as we have seen before [12]. During the shear, yarns are subject of lateral compression. This will certainly happen after the locking angle is reached: γ*=arcos(w/p), where w is the yarn width, p is the spacing (Long 2001a;Long et al., 2002). The lateral compression is present also before this moment, under a pressure of crimped interlaced warp and weft yarns. The lateral compression creates a pressure inside yarns, which results in an additional component of the transversal force acting on the yarns of the interlacing system. Therefore the transversal forces at yarns intersections will be n compressio tension bending where the first two terms are computed with [4] and [12]. The last term in [16] is computed based on the change of fibre volume fraction in the deformed configuration ( Figure 12). The geometric algorithm for non-orthogonal unit cell provides change of the shape of yarn cross-section and therefore the change of the fibre volume fraction. The pressure P inside the yarn and the transversal force are computed then using the experimental compression diagram of the yarn: where the superscripts "Wa" and "We" designate two intersecting yarns.

Components of the shear resistance
The components of the shear resistance are computed as follows.
Friction moment: where f is the coefficient of friction, reffective radius of the zone of friction for the normal force Q evenly distributed over a circle with radius R, d 2 is width of the intersecting warp and weft yarns.
Mechanical work of torsion: where C is the torsional rigidity of the yarn, τ is the full angle of torsion, computed by integration of rotation of vector a, determining orientations of the yarn crosssection axis, about the tangent to the yarn middle line t, over the yarn length. The torsional rigidity of the yarn can be estimated as where B is the bending rigidity of the yarn, d 1 -its thickness (Morton et al., 1970).
Mechancial work of (un)bending of the yarns: where δz is the difference between z-coordinate of the centre line of the yarn before and after the deformation (at a given coordinate along the yarn s).
Mechanical work of vertical displacement of the yarns (this displacement goes against the transversal forces Q): where l is the yarn length, and the integration is performed over the contact zones between warp and weft yarns, determined by the geometrical model (Lomov, Gusakov et al. 2000;Lomov et al., 2000). Figure 13 shows the results of simulation of the initial study of shear for a twill 1/3 glass (480 tex) fabric. The experimental curve was measured on KES-F shear tester. Bending properties of the yarns and friction coefficient was also measured on KES-F, resulting in the following values: B = 0.25 N mm 2 , f = 0.24. The results of compression tests on the yarns and yarn dimensions are given and discussed in (Lomov and Verpoest 2000). Torsional rigidity of the yarns was estimated with [20]. The comparison of the calculations and experimental values is quite satisfactory for this low shear region. Calculations for larger shear angles produce qualitatively satisfactory results, representing the main features of the shear process: fast increase of the shear force when the angle of shear increases to a locking state, increase of the fabric thickness for large shear angles due to the lateral compression of the yarns. However, quantitative comparison with shear frame experiments should involve careful control over parameters of the experiment, especially (pre)tension of the sample (Long 2001;Long, Souter et al. 2002) and will be a subject of future work. The model in the present formulation is not robust: it is sensitive to the variations of the input parameters, such as coefficient of friction or compression diagrams, which are measured with significant scatter.

Examples of calculations
To assess the sensitivity of the model to variations of parameters, consider a set of numerical experiments: simulations of shear of glass plain woven fabrics. Available literature data (Lomov and Verpoest 2000;Long 2001a;Long, Souter et al. 2002) and special measurements of parameters of glass rovings allow establishing of a set of master curves for estimation of dimensions of cross-sections, bending rigidity and compression curve of a glass roving with given linear density. The averaged curves, shown in Figure 14 and Figure 15, were used to determine input data for three types of yarns: 500, 1000 and 2000 tex. Fabrics with different spacing of the yarns (inverse to the ends/picks count) were considered, the ratio of the spacing to the yarn width being 0.6, 0.75 and 0.9. This ratio characterises the tightness of the fabric. Friction coefficient was set to be 0.2 in all the cases.

Conclusion
We have presented a family of models of internal structure of woven fabrics in relaxed and deformed state and of the load-deformation response for compression, bi-and uniaxial tension and shear. The models are based on the following common principles: − Description of the 2D and 3D weaves by a unified matrix coding methodology; − Decomposition of the crimped shape of the yarns into a set of elementary crimp intervals, with the shape of each interval corresponds to the solution of the minimum bending energy problem; − Calculation of the compressed dimensions of the yarn cross-sections by the balance between transversal forces generated by bending of the yarns and compression resistance of the yarns; − Calculation of the crimp heights of the yarns via a solution of minimum energy problem.
The models are implemented in WiseTex software package. The results of the simulations have been compared with the experimental observations and finite element simulations, resulting in a satisfactory agreement between them. However, more work has to be done to achieve robust and reliable predictions of the loaddeformation curves. A bottleneck here is creation of a database of measurements of properties of yarns in bending and compression, which can be readily used for predictive simulations of behaviour of preforms in composite forming.
The internal geometry description after deformation can further fed into flow modelling software, which allows computing local permeability of the deformed reinforcement, and micro-mechanical modelling to calculate homogenised local stiffness of the composite. This opens way to creation an integrated design software tool, combining local meso-analysis on the scale level of unit cell, homogenisation, resulting in permeability/stiffness parameters varying according to the local deformations and macro analysis of the filling/mechanical properties of the composite part.