Philippe Boisse

Simulations of Woven Composite Reinforcement Forming

Introduction

Complex preforms can be obtained by forming an initially flat textile composite reinforcement. Then the resin is injected on the preform in LCM processes (Liquid Composite Moulding) (Advani, 1994, Rudd & Long, 1997[START_REF] Parnas | Liquid Composite Molding[END_REF]. These processes and especially the RTM process (Resin Transfert Moulding) can be used to manufacture highly loaded composite part for aeronautical application (e.g. the helicopter frame and the motor blade presented in figures 1 and 2). Numerical optimization of products and production processes becomes increasingly important in the design phase of composite structures. Numerical simulations of the composite forming processes are an essential part of these optimization tools. They permit to determine the conditions of the feasibility of a process without defect (wrinkling, fracture of yarns, porosities…) but above all, they give the fibre orientations after shaping. This is mainly important because redistribution of the fibres is inevitable when double curved products are considered. The fibre orientations strongly influence the mechanical behaviour of the final part and the permeability of the reinforcement and thus the injection of the resin in the case of a liquid moulding process The first method used for analysing composite forming processes, and especially draping in woven composite reinforcements, is "kinematic models". The fibre distribution of a woven cloth is predicted on a given geometry based on a pin jointed net assumption (fishnet algorithm). The yarns of the fabric are assumed to be inextensible and the rotation between warp and weft yarns is free. The woven reinforcement is placed progressively from initial lines. Several packages are commercially available. This method is briefly described in section 2. It is fairly efficient for hand draping in classical prepreg fabrics, but the models do not account for load boundary conditions, for possible sliding of the fabric in relation to the tools, and for the mechanical behaviour of the woven reinforcement For a physical analysis of a composite forming process, the complete model must include all the equations for the mechanics, especially equilibrium, constitutive equations, and boundary conditions. These equations must be solved numerically, with some approximations. Finite Element Analysis of the composite forming process includes the tools modelling, the contact and friction between the different parts, and above all, the mechanical behaviour of the composite during forming. If these models can be numerically costly, problems of computation time are steadily reduced through improved processing capabilities. The main problem for the FE approach therefore lies in the requirement for accurate models of all the significant aspects of the forming process. The mechanical behaviour of fabrics is complex due to the intricate interactions of the yarns and fibres. It is a multi-scale problem. The macroscopic behaviour is very dependent on the interactions of yarns at the meso-scale (scale of the woven unit cell) and at the micro-scale (level of the fibres constituting yarns). Despite the many works in the field, there is no widely accepted model that accurately describes all the main aspects of a composite woven reinforcement mechanical behaviour. The approaches to model the forming of textile composite reinforcements belong to two main families that are related to the scale at which the analysis is made. The textile reinforcement is a set of yarns (or fibres). The analysis of the forming can be made considering and modelling each of these yarns (or fibres) and their interactions (contact with friction). In this case the approach is called discrete or mesoscopic. Of course the number of yarns is high and the interactions are complex. On the opposite, the continuous approaches consider a continuous medium juxtaposed with the fabric and the mechanical behaviour of which is equivalent to those of the textile reinforcement. This mechanical behaviour is complex because it concerns large strains and strong anisotropy. Furthermore, it strongly changes during the forming. The present chapter aims to present continuous and discrete approaches for composite reinforcements forming simulations. First two continuous approaches are described within a membrane assumption. The first one is based on a hyperelastic model and the second on a hypoelastic one. Then simulations of woven fabric forming based on a discrete approach are presented. Finally a semi-discrete approach which can be seen as an intermediate method between continuous and discrete ones is presented. This approach is extended to 3D interlock forming simulations. The advantages and drawback of the different approaches are discussed.

Kinematic models

"Kinematic models" also called fishnet algorithms are commonly used in industry for analysing composite forming processes, and especially draping in woven composite reinforcements [START_REF] Mark | The fitting of woven cloth to surfaces[END_REF][START_REF] Van Der Ween | Algorithms for draping fabrics on doubly curved surfaces[END_REF][START_REF] Long | A simulation of reinforcement deformation during the production of preform for liquid moulding processes[END_REF], Borouchaki et al, 2002). Several packages are commercially available. The method is based on the (strong) following assumptions : i/ The yarns are inextensible, ii/ there is no sliding at the intersection between warp and weft yarns, iii/ The rotations between warp and weft yarns are free, iv/ there is no sliding between the fabric and the tool. As shown on figure 3, the position of the node C can be determined if those of its neighbours A and C are already known. AC and BC have prescribed length. C is defined as the intersection of two geodesics coming from A and B and that intersect in C. This constitutes a small scalar problem, generally non-linear that can be solved very fast. The surface of the tool must be defined analytically of by curved elements. In order to initiate the draping as shown figure 3, it is necessary to position a first node and to fixe two initial draping directions. These directions are the symmetry axes if they exist. The result of the draping depends on these directions. This method is very fast and fairly efficient for hand draping in classical prepreg fabrics, but the models do not account for load boundary conditions, for possible sliding of the fabric in relation to the tools, and above all for the mechanical behaviour of the woven reinforcement.

Mechanical behaviour of textile composite reinforcement

The diameter of each fibre constituting the textile composite reinforcements is very small: 5 to 7 μm for carbon, 5 to 25 μm for glass, 10 to 20 μm for aramid. A yarn is made up of several thousands juxtaposed fibres (usually 3. 10 3 to 96. 10 3 ). These yarns are woven following standard weaves (plain, satin, twill) or more complex structures such as braiding or ply to ply interlock weaves (Figure 4c). An alternative consists in stitching a ply made of parallel fibres. This leads to the so-called NCF reinforcement (Non Crimp fabric) in which the fibres are not undulated (Figure 4d). The material resulting from this assembly of continuous fibres exhibits a very specific mechanical behaviour since relative motions are possible between the yarns and the fibres. The textile reinforcement pre-forming stage takes advantage of these possible motions. The forming is made on dry reinforcement (i.e. without resin) since it is performed before the injection stage. A woven fabric is intrinsically a multiscale material and, depending on the specific application of interest, one or more scales of the woven fabric have to be explored. Three scales can be distinguished (Figure 5). The macroscopic scale refers to the whole component level, with dimensions in the order of some centimetres to several meters (Figure 5a). At the mesoscopic scale, the woven reinforcement is seen as a set of yarns, respectively the warp and the weft (or fill) yarns in case of a woven fabric (Figure 5b) . Consequently, the corresponding working scale is the one of the yarn dimension, typically one to several millimetres. For periodic materials, mesoscopic models consider the smallest elementary Each yarn is made up of thousands of continuous fibres which interact (Figure 5c), and thus the interactions of the reinforcement can be analyzed at the microscopic scale. At the microscopic level, the characteristic dimension is about one to several micrometers. This is the only scale at which the material is actually continuous.

Although the fibrous reinforcement is not strictly continuous because of possible relative sliding between fibres, several mechanical behaviour models have been proposed that consider the textile reinforcement as an anisotropic continuum (Spencer, 2000;[START_REF] Dong | Processing of composites: simulations of the draping of fabrics with updated material behaviour law[END_REF][START_REF] Yu | Non-orthogonal constitutive equation for woven fabric reinforced thermoplastic composites[END_REF][START_REF] King | A continuum constitutive model for the mechanical behavior of woven fabrics[END_REF][START_REF] Peng | A continuum mechanics-based non-orthogonal constitutive model for woven composite fabrics[END_REF][START_REF] Thije | Large deformation simulation of anisotropic material using an updated Lagrangian finite element method[END_REF], Badel et al, 2009) Nevertheless there is no widely accepted model that describes accurately all of the main aspects of fabric mechanical behaviour. Actually, such a model must convey the specificities of the composition of the textile made of yarns and fibres and above all take into account the variation of the properties during the forming. These changes are very large because of the variations of fibre directions and of the local lateral compression of the yarns due to the forming.

Experimental texts

The specific mechanical behaviour of textile composite reinforcements has led to the development of specific experimental tests in order to quantify the tensile, in-plane shear and bending behaviour.

Biaxial tensile behaviour

The tensile behaviour of woven material is specific mainly because of the decrimping of tows when they are stretched. This leads to tensile behaviour non-linearities. The fabric is much softer than the tow for small axial strains. Because of the weaving, the decrimping phenomenon in warp and weft directions are interdependent and the tensile behaviour is biaxial. Some biaxial tests have been developed in order to measure these properties [START_REF] Kawabata | The Finite Deformation Theory of Plain Weave Fabrics Part I: The Biaxial Deformation Theory[END_REF][START_REF] Buet-Gautier | Experimental analysis and modeling of biaxial mechanical behavior of woven composite reinforcements[END_REF][START_REF] Carvelli | Mechanical modelling of monofilament technical textiles[END_REF], Willems et al, 2008).

In-plane shear behaviour

Two experimental tests are used to determine the in-plane shear behaviour of textile composite reinforcements: the picture-frame (Figure 7) and the bias-extension tests (Figure 8). A great literature is dedicated to those tests [START_REF] Prodromou | On the relationship between shear angle and wrinkling of textile composite preforms[END_REF][START_REF] Rozant | Drapability of dry textile fabrics for stampable thermoplastic preforms[END_REF][START_REF] Potter | Bias extension measurements on cross-plied unidirectional prepreg[END_REF][START_REF] Lebrun | Evaluation of bias-extension and picture-frame test methods for the measurement of intraply shear properties of PP/glass commingled fabrics[END_REF][START_REF] Sharma | Characterisation of material properties for draping of dry woven composite material[END_REF][START_REF] Peng | Experimental and numerical analysis on normalization of picture frame tests for composite materials[END_REF][START_REF] Harrison | Shear characterisation of viscous woven textile composites: a comparison between picture frame and bias extension experiments[END_REF][START_REF] Lomov | Picture frame test of woven composite reinforcements with a full-field strain registration[END_REF][START_REF] Launay | Experimental analysis of the influence of tensions on in plane shear behaviour of woven composite reinforcements[END_REF], Lomov et al, 2008[START_REF] Cao | Characterization of Mechanical Behavior of Woven Fabrics: Experimental Methods and Benchmark Results[END_REF] mainly because the in-plane shear is the most dominant deformation mode in woven composite forming when the manufactured part is doubly curved. The shear angle can reach 50° (and even more in some cases such as presented in section 5.6). For large values, wrinkling can occur depending on the process parameters and on the material properties. In addition, the two tests are difficult both from the experimental point of view and concerning the interpretation of the results. For these reasons an international benchmark has been launched recently. Eight laboratories of six different countries have performed picture frame and bias-extension tests on the same textile composite reinforcements [START_REF] Cao | Characterization of Mechanical Behavior of Woven Fabrics: Experimental Methods and Benchmark Results[END_REF].

The picture-frame (or trellis-frame) is made of four hinged bars. The fabric specimen is initially square and the tows are parallel to the bars. Consequently it is theoretically subjected to pure in-plane shear and the shear angle γ is function of the displacement d.

1d 2arccos 22 L 2 π ⎛⎞ γ= - + ⎜⎟ ⎝⎠ (1)
Neglecting the dissipation due to friction in the hinged bars, the in-plane shear moment on a unit shell M s is related to the force on the frame F using the equality of the power expressions. S c is the surface of the unit woven cell in the initial state.

()

c s γγ S M γ =F cos -sin 22 2L ⎛⎞ ⎜⎟ ⎝⎠ (2)
The bias-extension test is an alternative to the picture-frame test. It consists in clamping a rectangular specimen of woven fabric with warp and weft directions initially oriented at 45° with respect to the tensile load applied by a tensile machine (Figure 8). The initial length of the specimen L must be larger than twice the width ℓ. The zone C in the centre of the specimen is submitted to a pure shear γ if the yarns are assumed to be inextensible. That is a correct assumption for the type of yarns used as composite reinforcements. This inextensibility imposes that the shear angle in the zone B is half the value in the central region C. The shear angle γ is related to the specimen elongation d in equation (c). The inplane shear moment is related to the force on the frame F in equation (d).

Dd 2 arccos 2 2D

⎛⎞ ⎜⎟ ⎝⎠ π+ γ= -

(3) () ( ) FDS c Mc o s s i n M ss 2D 2 2 2D 2 γ γγ ⎛⎞⎛ ⎞ γ= - - ⎜⎟⎜ ⎟ --⎝⎠⎝ ⎠ ˋ(4)
A shear curve M s (γ) measured in the case of glass plain weave is presented in Figure 9. In the first part of the curve (i.e. for small shear angles), the in-plane shear stiffness is small. For larger shear angles this rigidity increases and becomes significant. The optical field measurements performed within the tow show that during the first part of the loading, the tows rotate in a rigid body motion (Figure 9b). When the shear angle becomes larger lateral contacts between the yarns occur [START_REF] Dumont | Field measurements for shear tests on woven reinforcements[END_REF]. The tows are progressively compressed and the shear rigidity increases significantly. This increase of shear stiffness leads to wrinkling onset. The corresponding shear angle is called locking angle (in order of 40°-45° for textile composite reinforcements [START_REF] Prodromou | On the relationship between shear angle and wrinkling of textile composite preforms[END_REF][START_REF] Cao | Characterization of Mechanical Behavior of Woven Fabrics: Experimental Methods and Benchmark Results[END_REF].

Bending behaviour

The bending stiffness of the fibrous reinforcement is very low due to the possible motion between fibres. It is often neglected and membrane approaches are often used. Nevertheless this bending stiffness is important to obtain accurate wrinkle shape. Specific tests have been developed for these textile materials, the bending stiffness of which is much smaller than the one given by the plate theory [START_REF] Kawabata | The Standardization and Analysis of Hand Evaluation[END_REF][START_REF] Lahey | Mechanical Modeling of Fabrics in Bending[END_REF][START_REF] De Bilbao | Experimental Study of Bending Behaviour of Reinforcements[END_REF]. Figure 10 shows a cantilever bending test (de Bilbao et al, 2010). The own weight of the fabric defines the bending moments for a given length of the specimen. The curvatures are deduced from the measured geometry. These measurements, made for different length of a specimen give the bending moment in function of the curvature. They are made in warp and weft directions.

Continuous approach for the simulation of textile composite forming

In the continuous approaches, a woven fabric is seen as a continuous material with a specific mechanical behaviour, including high anisotropy and the ability to exhibit very large shearing and bending deformations. Investigation at the macroscopic level is the most popular for reinforcement forming simulations, as it allows using finite elements codes with standard shell or membrane elements and does not ask the description of the internal textile material structure. Despite the large amount of work in this field (Spencer, 2000;[START_REF] Dong | Processing of composites: simulations of the draping of fabrics with updated material behaviour law[END_REF][START_REF] Yu | Non-orthogonal constitutive equation for woven fabric reinforced thermoplastic composites[END_REF][START_REF] King | A continuum constitutive model for the mechanical behavior of woven fabrics[END_REF][START_REF] Peng | A continuum mechanics-based non-orthogonal constitutive model for woven composite fabrics[END_REF][START_REF] Thije | Large deformation simulation of anisotropic material using an updated Lagrangian finite element method[END_REF], Badel et al, 2009) there is no widely accepted model that accurately describes all aspects of the mechanical behaviour of fabrics. Two continuous approaches are described below.

Hyperelastic behaviour

The formulation of a hyperelastic behaviour law lies on the proposition of a potential energy from which derives the hyperelastic constitutive model. This potential aims to reproduce the non linear mechanical behaviour of textile composite reinforcements. The proposed potential is a function of the right Cauchy Green and structural tensor invariants defined The third term s W in equation ( 5) is a function of the second mixed invariants of C.

2 12 11 22 12 1 IT r ( C L C L ) c o s II = ⋅⋅ ⋅ = θ (8)
The second Piola Kirchhoff stress tensor is derived from this potential [START_REF] Basar | An overview of semidiscretisation and time integration procedures[END_REF]:

W S2 C ∂ = ∂ (9)
And in the case of the present potential (5):

( ) (10) In order to define the form of the potential two complementary assumptions are made taking into account the specific woven fabric behaviour and its deformation modes. As assumed above, i/ The tensions in the yarns and the in-plane shear are independent. ii/ The tensions in the warp and weft directions are uncoupled. The potential has to vanish in a stress free configuration. Polynomial functions of the invariants are considered in the present work. The global form of the proposed potential energy is given by: ( )

( ) ( ) rs t j1 i1 k i1 j 2k 1 2 i0 j0 k1 11 1 WC A I 1 B I 1 CI i1 j1 k + + === =-+-+ ++ ∑∑∑ (11)
The resulting second Piola Kirchhoff tensor is: For strain-free configuration, stresses have to vanish. This condition imposes:

rs ij i0 j0 A0 ; B0 == = = ∑∑ (13)
To determine the constants A i , B j and C k , three experimental tests are necessary: two tensile tests in the warp and weft directions and one in-plane pure shear test. The details of the calculations to obtain equations ??? to ??? are given in (Aimene et al, 2010). In this paper it is also shown that the form of the potential given above gives correct results concerning the direction of the loads on the boundary of a picture frame while other forms of the potential lead to boundary loads that are not correct for a woven fabric.

The proposed hyperelastic model is implemented in a user routine VUMAT of Abaqus/Explicit and it is applied to membrane elements. The simulation of a hemispherical punch forming process is performed in the case of strongly unbalanced twill [START_REF] Daniel | Forming simulation of very unbalanced woven composite reinforcements[END_REF]. The warp rigidity is 50 N/yarn and the weft rigidity is 0.2 N/yarn. The experimental deformed shape are shown figure 12 (a) together with the results of the simulation figure 12 (c). The computed deformed shape (made using the hyperelastic model proposed above) is in correct agreement with the experimental one. Especially the strong difference of the deformation in warp and weft directions is well verified.

Hypoelastic behaviour

Hypoelastic models have been proposed for material at large strain [START_REF] Truesdell | Hypo-elasticity[END_REF], Xiao et al, 1988)

∇ = C:D σ ( 14
)
where D and C are the strain rate tensor and the constitutive tensor, respectively. ∇ σ , called the objective derivative of the stress tensor, is the time derivative for an observer who is fixed with respect to the material.

( )

TT d dt ∇ ⎛⎞ = ⋅⋅ ⋅ ⋅ ⎜⎟ ⎝⎠ QQ Q Q σσ ( 15
)
Q is the rotation from the initial orthogonal frame to the so-called rotating frame where the objective derivative is made. The most common objective derivatives are those of Green-Naghdi and Jaumann. They use the rotation of the polar decomposition of the deformation gradient tensor ⋅ F=R U, (standard in Abaqus explicit), and the corotational frame, respectively. These are routinely used for analyses of metals at finite strains (Belytschko et al, 2000). It has been shown that, in the case of a material with one fibre direction the proper objective rotational derivative is based on the rotation of the fibre (Badel et al, 2009) In the current configuration, the unit vectors in the warp and weft fibre directions are respectively:

1 1 1 = 0 0 F.f f F.f 2 2 2 = 0 0 F.f f F.f (17)
Where 00 12 (,) ee and 00 12 (,) ff are assumed to coincide initially (Figure 3). Two orthonormal frames based on the two fibre directions are ( ) The strain increment dε is given as a code's output in calculation loop from time t n to time t n+1 . (The matrix of the components of this strain increment is given in the GN frame in the case of Abaqus explicit, but it could be any other frame). The components of the strain increment in the two frames g and h are considered (α and β are indexes taking value 1 or 2): The fibre stretching strain and the shear strain are calculated for the two frames g and h.

The stress at time t n+1 in the fabric is the addition of the stresses in the two fibre frames:

( ) ( ) More detail on this approach can be found in (Badel et al, 2009[START_REF] Khan | Numerical and experimental analyses of woven composite reinforcement forming using a hypoelastic behaviour. Application to the double dome benchmark[END_REF]. This approach is used to simulate the forming double dome shape corresponding to an international benchmark [START_REF] Khan | Numerical and experimental analyses of woven composite reinforcement forming using a hypoelastic behaviour. Application to the double dome benchmark[END_REF]). An experimental device has been realised in INSA Lyon in order to perform this forming (Figure 13). The woven fabric is a commingled glass/polypropylene plain weave that has been tested in the material benchmark study conducted recently [START_REF] Cao | Characterization of Mechanical Behavior of Woven Fabrics: Experimental Methods and Benchmark Results[END_REF]. The computed and experimental geometries after forming are compared figure 14 and 15. The measured and numerical geometries and shear angles are in good agreement.

n1

Discrete approach for the composite reinforcement forming

In discrete modelling (also called meso-modelling in the case of textile material), the modelling does not directly concern the textile material but each fibre bundle. This one is modelled by elements simple enough to render the computation possible because it concerns the forming of the whole composite reinforcement and the number of yarns and contacts between these yarns is very large. The interactions between warp and weft directions are taken into account explicitly by considering contact behaviour and relative motions between the yarns are possible [START_REF] Pickett | Simplified and Advanced Simulation Methods for Prediction of Fabric Draping[END_REF][START_REF] Duhovic | Simulating the deformation mechanisms of knitted fabric composites[END_REF][START_REF] Boukaber | Discrete models of woven structures. Macroscopic approach[END_REF]. At the microscopic level, each fibre is satisfactorily described as a beam but this approach is time consuming. The main difficulty is the great number of contacts with friction that have to be taken into account, especially for a woven fabric. For this reason, only very small elements of the fabric have been modelled to date [START_REF] Durville | Numerical simulation of entangled materials mechanical properties[END_REF][START_REF] Miao | Mechanics of textile composites: Microgeometry[END_REF]. Nevertheless, this approach is promising because it does not necessitate any assumptions regarding the continuity of the material, the specific mechanical properties resulting at the macroscopic level naturally follow the displacements and deformations of the yarns and it provides an interesting way of taking the weaving operation into account. The fibres constituting the yarns can be modelled directly, but their very large number (3K to 48K per yarn) requires that the computations are made for a number of fibres per yarn significantly smaller than in reality. An alternative possibility is to use a continuous behaviour for each yarn (meso-modelling). This implies that the fibrous nature of the yarn is taken into account in this model especially in order to have rigidities in bending and transverse compression very small in comparison to the tensile stiffness. In any case, a compromise must be found between a fine description (which will be expensive from the computation time point of view) and a model simple enough to compute the entire forming process. Figure 16 the finite element model used for discrete simulations of forming processes (216 dof (degrees of freedom)). It is compared to another FE model of the unit cell used in (Badel et al, 2009)(Figure 16(a)) to analyze the local in plane shear of a plain weave unit cell (47214 dof). It cannot be considered (at least today) to use this last FE model to simulate the forming of a composite reinforcement that is made of several thousands of woven cells. In the simplified unit cell (Figure 16(b)) each yarn is described by few shell elements and the contact friction and possible relative displacement of the yarns are considered. The in-plane mechanical behaviour is the same as the one defined in (Badel et al, 2009). The bending stiffness is independent of the tensile one and very much reduced in comparison to the one given by plate theories. Two examples are presented in figures 17 and 18 based on a discrete modelling using the unit cell of figure 16(b) [START_REF] Gatouillat | Advantages of the meso/macro approach for the simulation of fibre composite reinforcements[END_REF] The first one is a picture frame test for which the wrinkles appear naturally in the simulation when the shear locking angle is reached. It must be noticed that the in-plane shear behaviour of the fabric is not an input data of the analysis and does not need to be known. It results at the macroscopic level of contact and friction between the yarns and lateral compression of the yarns. Figure 18 shows the results of a hemispherical forming simulation. It must be said that this study concerning forming simulation at the meso-scopic scale is beginning at INSA Lyon. If the discrete or mesoscopic modelling is a promising approach because a large part of the mechanical specificity of fabric behaviour is due to yarn and fibre interactions, and following fibre directions is simpler than for continuous models, it must be recognized that the forming simulations made with approaches that permits the relative sliding of the yarns in contact are not many.

The semi-discrete finite elements for the composite reinforcement forming

This approach can be seen as intermediate between the continuous and the discrete approaches. The textile composite reinforcement is seen as a set of a discrete number of unit woven cells submitted to membrane loadings (i.e. biaxial tension and in-plane shear) and bending (Figure 19) [START_REF] Hamila | A semi-discrete shell finite element for textile composite reinforcement forming simulation[END_REF] In any virtual displacement field η such as η = 0 on the boundary with prescribed loads, the virtual work theorem relates the internal, exterior and acceleration virtual works: function of the gradient of the virtual displacement field. T 11 and T 22 are the tensions on the unit woven cell in warp and weft directions. M 11 and M 22 are the bending moments on the woven cell respectively in warp and weft directions. M s is the in-plane shear moment. The loads on an edge of the woven unit cell (presented Figure 19) result in the tensions T 11 and T 22 in one hand and shear forces in the other hand. This shear forces on a warp and weft sections have a moment at the centre of the woven unit cell in the direction normal to the unit cell. The component of this moment is called in-plane shear moment M s . This quantity is conjugated to the in-plane shear angle γ . The internal virtual work of in plane shear is directly given from M s and the virtual shear angle (Equation 30). In the case of a textile material, the shear angle γ is a significant and clearly defined quantity and it is interesting to express the internal virtual work of in plane shear in function of this quantity.

ext int acc W() W() W() η -η = η (27)
The mechanical behaviour of the textile reinforcement defines a relation between the loads T αα , M s , M αα and the strain field. The experimental tests specific to textile composite reinforcements that have been presented in Section 4 are used to obtain T 11 , T 22 , M s , M 11 and M 22 in function of ε 11 , ε 22 , γ, χ 11 and χ 22 .

An alternative consists in virtual tests i.e. in 3D simulations of the deformation of a unit woven cell submitted to elementary loadings such as biaxial tensions or in plane shear (Badel et al, 2008, Badel et al, 2009) 7.1 Shell finite element made of woven cells A three node shell finite element have been defined from the simplified form of the principle of virtual works given in equations 28 to 31. The details of its formulation are given in [START_REF] Hamila | A semi-discrete shell finite element for textile composite reinforcement forming simulation[END_REF]. It is summarized below.

A finite element interpolation is introduced within the principle of virtual work. The displacement u and virtual displacement η of any point within an element are in the form: N is the interpolation matrix of the element under consideration and e u and e η the single column matrices of its nodal displacements and virtual displacements respectively. Equation 28 leads to:

( )

tsb int int int ext 0 + ++ -= Mu F F F F $$ (33)
M is the mass matrix, u is the single column matrices of the nodal displacements. tsb int int int ,, FFF are the single column matrices of the nodal internal forces respectively for tension, shear and bending. This dynamic equation can be solved using an explicit scheme (central differences):

( ) i1 i 1 i t i s i b i De x t i n ti n t i n t +- =- - - uM FF F F $$ (34) ( ) i1 2 i1 2 i1 i i1 1 tt 2 +- - + =+ Δ + Δ uu u $$ $ $ (35) i1 2 i1 i i t + + = +Δ uu u $ (36)
There is no system to solve since M D is a diagonal matrix calculated from M [START_REF] Zienkiewicz | The finite element method[END_REF]. This explicit scheme requires the time step to be small enough to insure the stability of the scheme (Belytschko, 1983). It is effective for many dynamic applications and also in material forming. For the sake of numerical efficiency, the speed can be larger than the real one under the condition that the dynamic effects do not modify the solution.

The 3 node shell finite element M 1 M 2 M 3 made up of ncelle woven cells is considered (Figure 20). The vectors k 1 =AM 2 and k 2 =BM 3 respectively in warp and weft directions are defined.

() M = γη = ∑ η F (39)
The internal in-plane shear force components are calculated from the in-plane shear moment:

( ) ( ) se s int ij ij Fn c e l l e B M γ = γ (40)
In order to avoid supplementary degrees of freedom and consequently for numerical efficiency, the bending stiffness is taken into account within an approach without rotational degree of freedom [START_REF] Onate | Rotation-free triangular plate and shell elements[END_REF][START_REF] Sabourin | Detailed formulation of the rotation-free triangular element "S3" for general purpose shell analysis[END_REF]. In theses approaches the curvatures of the element are computed from the positions and displacements of the nodes of the neighbouring elements (Figure 20). The internal virtual work of bending on the element defines the element nodal bending internal forces be int F :

ncelle pp p pp peT be 11 22 i n t p1 ()M L ()M L 11 22 = χη + χη = ∑ η F (41)
The internal bending force components are calculated from the bending moments M 1 and M 2 :

( )

be 12 int celle 1km 1 2km 2 22 km 12 LL Fn B b M B b M kk ⎛⎞ ⎜⎟ =+ ⎜⎟ ⎝⎠ (42)
The finite element presented above is used to simulate the hemispherical forming of the very unbalanced fabric presented in section 5.1. The experimental forming by a [START_REF] Daniel | Forming simulation of very unbalanced woven composite reinforcements[END_REF]. A 6 kg ring was used as blank-holder avoiding reinforcement wrinkling in the hemispherical zone (Figure 21a). The experimental shape obtained after forming is shown in Figure 21b. In warp direction (with the strongest rigidity) large fabric sliding is observed relatively to the die. On the contrary, in the weft direction (weak direction) no edge movement is depicted and the yarns are subjected to large stretch deformations. The computed shape after forming is shown in Figure 21c. It is in good agreement with experiments. Especially the extension ratio at the centre of the hemisphere (l weft /l warp = 1.8) is correctly computed. The shape of the many wrinkles in the flat part of the textile is also properly simulated.

Extension to 3D interlock textile reinforcements

When the thickness of a composite part is large, the use of these laminated composites is restricted by manufacturing problems and their low resistance to delamination cracking. To overcome these difficulties composites with 3D fibre architecture called ply to ply interlock fabric have been proposed [START_REF] Tong | 3D Fibre reinforced polymer composites[END_REF]. This material is not fully 3D since there is no third yarn set in the transverse direction but the properties through the thickness are much improved. Above all, the possible delaminations of the 2D laminated composites are overcome. (Figure 22). The semi-discrete approach has been extended to 3D hexahedral finite elements for interlock forming simulations; These elements are made of yarns as shown

Conclusion

Different approaches for composite forming simulations have been presented in this chapter. They are continuous discrete or semi-discrete element. These different approaches are based on the strong multi-scale nature of the textile reinforcements. The discrete approach is attractive and promising. The very specific mechanical behaviour of the textile material due to the contacts and friction between the yarns and to the change of direction is explicitly taken into account. If some sliding occurs between warp and weft yarns, they can be simulated. This is not possible by the continuous approaches that consider the textile material as a continuum. This is an important point because it can be necessary to prevent such a sliding in a process. Nevertheless, the main difficulty of the discrete approach is the necessary compromise that must be done between the accuracy of the model of the unit woven cell and the total number of degrees of freedom. The modelling of the unit cell must be accurate enough to obtain a correct macroscopic mechanical behaviour, but the number of degrees of freedom of each cell must remain small in order to compute a forming process for which there will be thousand of woven cells. There are a lot of improvements to achieve in the meso-modelling of different textile reinforcements. The continuous increase of the computer power is a strong argument in favour of this approach. The continuous approach is the most commonly used in composite reinforcement forming today. The main advantage is to use standard shell or membrane finite element. The only mechanical behaviour has to be specified in order to take the very particular behaviour of textile materials into account. Many models exist, but none of them is clearly admitted. The semi-discrete approach aims to avoid the use of stress tensors and directly define the loading on a woven unit cell by the warp and weft tensions and by in-plane shear and bending moments. These quantities are simply defined on a woven unit cell and above all they are directly measured by standard tests on composite reinforcements (biaxial tension, picture frame, bias extension and bending tests).

Advani SG. Flow and rheology in polymeric composites manufacturing. Elsevier, 1994. Aimène Y, Vidal-Sallé E, Hagège B, Sidoroff F, Boisse P, A hyperelastic approach for composite reinforcement large deformation analysis, Journal of Composite Materials Vol. 44, No. 1/2010, 5-26

Fig. 1 .

 1 Fig. 1. Preform/RTM parts in NH90 (Dumont et al. 2008) (Courtesy of Eurocopter, EADS Group)

Fig. 3 .

 3 Fig. 3. Fishnet algorithm: Calculation of the position of point C knowing those of A and B

Fig. 4 .

 4 Fig. 4. Textile composite reinforcements (a) plain weave, (b) twill weave (c), interlock, (d) NCF

Fig. 5 .

 5 shows a biaxial tensile device using a cross shape specimen[START_REF] Buet-Gautier | Experimental analysis and modeling of biaxial mechanical behavior of woven composite reinforcements[END_REF]. The measurements of tensions in warp and weft directions T 1 (ε 11 , ε 22 ) and T 2 (ε 11 , ε 22 ) are shown Figure5bfor different ration between warp and weft strains. It has been experimentally shown that the influence of the shear angle on the tensile behaviour is usually weak and can be neglected[START_REF] Buet-Gautier | Experimental analysis and modeling of biaxial mechanical behavior of woven composite reinforcements[END_REF].

  Fig. 6. (a) Biaxial tensile test on cross-shaped specimen (b) Load versus strain for carbon twill weave. k = ε warp / ε weft[START_REF] Buet-Gautier | Experimental analysis and modeling of biaxial mechanical behavior of woven composite reinforcements[END_REF]

Fig. 7 .Fig. 8 .

 78 Fig. 7. Picture frame test. (a) Experimental device. (b) Kinematics of the test

  Fig. 9. (a) Shear curve of glass plain weave, (b) relative displacement field inside a yarn

Fig. 10 .λ

 10 Fig. 10. Cantilever bending test from the fibre directions. This potential is based on the assumption that tensile and shear strain energies are uncoupled. It is the sum of three terms.

Fig. 11 .

 11 Fig. 11. Fibres axes and GN axes after deformation. Initially these axes are superimposed.

Fig. 12 .

 12 Fig. 12. Deformed shape in the case of an unbalanced fabric (experimental shape (a); Result of the simulation (c)) A membrane assumption is used. The Green-Naghdi's frame (GN) is the default work basis of ABAQUS/Explicit. Its unit vectors 12 (,) ee in the current configuration are updated from the initial orientation axes, 00 12 (,) ee using the proper rotation R : 00 1122 = ⋅= ⋅ eR e eR e (16)

Fig. 13 .Fig. 15 .

 1315 Fig. 13. Double dome forming: experimental device

  Fig. 16. Meso-modelling of a unit cell of a plain weave. (a) FE model for the analysis of the behaviour of the unit cell. 47214 Dof. (b) FE model for simulations of the whole composite reinforcement forming. 216 Dof.

Fig. 17 .

 17 Fig. 17. Simulation of a picture frame test using the unit cell model of figure 16(b)

Fig. 20 .Fig. 21 .

 2021 Fig. 20. Three node finite element made of unit woven cells,

  Figure 23. The simulation of a thick twisted plate made of interlock textile reinforcement is shown figure 24. The formulation of these finite elements for interlock forming simulations is given in (De Luycker et al, 2009).

Fig. 22 .

 22 Fig. 22. Example of complex layer interlock weave

Fig. 24 .

 24 Fig. 24. Forming of a 3D interlock twisted plate

  are the stiffness in warp and weft fibre directions respectively and G the inplane shear stiffness of the fabric (They are not constant, especially G depends strongly on the in plane shear). Following the scheme of Hughes and Winget the stresses are then integrated on the time increment from time t n to time t n+1[START_REF] Hughes | Finite rotation effects in numerical integration of rate constitutive equations arising in large deformation analysis[END_REF]):

	g 11 dε =⋅⋅ gd ε g 11	g 12 dε =⋅⋅ gd ε g 12			(19)
	h 22 dε =⋅⋅ hd ε h 22	h 12 dε =⋅⋅ hd ε h 12			(20)
	From these strain components the axial stress component and shear stress components are
	calculated in each frame g and h:
	gg 11 dE d σ =ε 11 g	gg 12 12 dG d σ =ε			(21)
	hh 22 dE d σ =ε	h 22	hh 12 12 dG d σ =ε			(22)
	E g and E ( )( ) n1 gg 11 11 + σ= σ	n	n1 / 2 11 g + + σ d	( ) ( ) n1 gg 12 12 + σ= σ	n	n1 / 2 12 g + + σ d	(23)
	( )( ) n1 hh 11 11 + σ= σ	n	n1 / 2 11 h + + σ d	( ) ( ) n1 hh 12 12 + σ= σ	n	n1 / 2 12 h + + σ d

h

  is the index of the direction (i=1 to 3), j is the index of the node (j=1 to 3). B 1ij and B 2ij are strain interpolation components. They are constant over the element because the interpolation functions in equation 32 are linear in the case of the 3 node triangle. The internal virtual work of in-plane shear on the element defines the element nodal tensile internal forces se

	The internal virtual work of tension on the element defines the element nodal tensile
	internal forces te int F :
	ncelle p1 = ∑	pp ()T L p 11 11 εη	p + ε η ()T L p peT te 22 i n t 22 = η F	(37)
	The internal tensile force components are calculated from the tensions T 1 and T 2 :
	( ) te int Fn c e l l e B T 1ij 1 ij ⎛⎞ 12 2ij 2 22 12 LL B T kk ⎜⎟ =+ ⎝⎠ ⎜⎟	(38)
	int F :
		ncelle	pp	eT se
					si n t
		p1

i


Benjamin Legros

Reservation, a tool to reduce the balking effect and the probability of delay

Keywords: queue, balking, reservation, threshold policy

We investigate a threshold reservation policy implemented within a single customers' class. From an explicit performance analysis, we prove that the potential of reservation is into the reduction of the balking effect together with a higher servers' utilization. However, it also may result in a higher expected waiting time and more abandonment. We conclude that reservation can be efficiently implemented in large systems, under high workload situations, with a low waiting aversion and a low impatience.

Introduction

In a system where all resources are occupied, the service of new arriving customers has to be delayed. The waiting aversion leads some customers to reconsider asking for service in case of a non-immediate service.

This phenomenon is referred to as the balking effect. Balking can be significant in situations where customers know that they can easily find a similar service elsewhere, eventually without wait. For instance in a street with many restaurants, a significant proportion of customers may refuse waiting for an available table if they can find another one elsewhere. Another example is call centers; [START_REF] Brown | Statistical analysis of a telephone call center: A queueing-science perspective[END_REF] show that the balking phenomenon results in a significant loss of customers.

Since the loss of customers is undesirable, it may be useful to develop strategies which allows the service provider to maintain an amount of resource availability if preemption is not possible. In a context where a unique group of servers has to serve urgent and non-urgent customers, it has been shown that a reservation strategy can significantly improve the service level of urgent customers. The idea of a reservation strategy is to force some servers to remain available for new arriving urgent customers even if some non-urgent customers are waiting. This strategy is sometimes implemented in health care systems (see [START_REF] Dobson | Reserving capacity for urgent patients in primary care[END_REF]) and in call centers with inbound and outbound calls (see [START_REF] Koole | Call Center Optimization[END_REF]).

An important stream of literature is devoted to the analysis of such reservation strategies with two customer's classes. These strategies are referred to as blending policies. Some papers focus on performance evaluation, and others address the analysis of blending policies or staffing decisions. [START_REF] Deslauriers | Markov chain models of a telephone call center with call blending[END_REF] develop various continuous Markov chain models for a call center with inbound and outbound calls. The authors consider a threshold policy and characterize the rate of outbound calls and the waiting time distribution of inbound calls. Other papers address the analysis of reservation policies. [START_REF] Gans | A call-routing problem with service-level constraints[END_REF] and [START_REF] Bhulai | A queueing model for call blending in call centers[END_REF] prove that a threshold policy on the number of idle agents is optimal to maximize the outbound throughput under a service level constraint on the inbound waiting time. Similar results are also found in [START_REF] Legros | Adaptive threshold policies for multi-channel call centers[END_REF], for a non-stationary model where inbound calls arrive according to a non-homogeneous Poisson process or in [START_REF] Legros | Optimal scheduling in call centers with a callback option[END_REF] in a setting with a callback option.

[15] consider a large call center with a reservation model and propose a logarithmic safety staffing rule, combined with a threshold control policy to ensure that agents' utilization is always close to one with always idle agents present.

In a broaden perspective, reservation is the idea that some resources should remain idle when some work is available. For instance in a context with heterogeneous servers and the objective to minimize the time spent in the system, it is optimal to consider non-workconserving policies where slow servers are used only if the queue length exceeds some given thresholds. A large literature investigates this question through the so-called "slow-server" problem (e.g., see [START_REF] Krishnamoorthi | On poisson queue with two heterogeneous servers[END_REF], [START_REF] Larsen | Control of a heterogeneous two-server exponential queueing system[END_REF], [START_REF] Rykov | Monotone control of queueing systems with heterogeneous servers[END_REF]). Other illustrations are overflow policies where agents are initially reserved for only one customer type but can treat another customer type in case of high congestion (e.g., see [START_REF] Barth | Analysis of two-level support systems with time-dependent overflow-a banking application[END_REF], [START_REF] Koole | Optimization of overflow policies in call centers[END_REF]). When there is a switching time or a switching cost to change the job type in service, it may also be preferred for a server to remain idle and ready to serve the same job type instead of serving the other job type. For this purpose, cycling reservation strategies are studied (e.g., see [START_REF] Tran-Gia | A discrete-time analysis of the cyclic reservation multiple access protocol[END_REF], [START_REF] Boxma | Cyclic reservation schemes for efficient operation of multiplequeue single-server systems[END_REF]).

The particularity of our article compared to the existing ones on the same subject is to consider a threshold reservation strategy with a unique class of customers who may balk. The choice for a reservation threshold policy follows from the knowledge that this type of policies outperforms randomized policies since they take into account the system state. In addition, threshold policies are deterministic; for each state there is a unique action. This makes these policies easy to implement in practice. Finally, threshold policies are controlled by only one parameter which also makes them possible to analyze and implement.

We choose to adapt the threshold policy developed in the call blending references for a unique class of customers. We consider a multi-server single queue with infinite capacity and s identical, parallel servers.

The arrival process of customers is Poisson with rate λ. Service times are independent and exponentially distributed with rate µ. At a customer's arrival, if at least one server is available then this customer is directly served, otherwise with probability r he/she accepts waiting and is routed to a first-come-first-served queue. The reservation policy is of threshold type. We define a threshold c on the number of servers. After a service completion, if the number of busy servers is higher than or equal to s -c then no customer is routed to service. Otherwise, if the queue is non-empty then the first customer in line is routed to service. In other words, c servers are reserved for new arriving customers (0 ≤ c < s).

The objective of this article is to evaluate the performance of this system and determine the effect of reservation on the proportion of lost customers and on the waiting time. The idea is to show the motivations and the risks in implementing such a policy. In Section 2, we present the Markov chain associated to this system and develop a method to obtain the explicit expressions of the stationary probabilities. Section 3 is devoted to the impact of reservation on balking. We prove that reservation can reduce the proportion of lost customers. This might be unexpected by service providers. By forcing some servers not to work, the overall productivity may increase. The improvement can be significant in large systems, under high workload situations and with a low customers' waiting aversion. In Section 4, we develop a method to compute the Laplace Stieltjes Transform (LST) of the waiting time distribution of an arriving customer. This allows us to obtain explicit expressions of the first and second moments of the waiting time. As expected, although reservation helps to reduce the probability of delay, it deteriorates the expected waiting time due to the existence of situations where a server may remain idle while a customer is waiting. In addition, in Section 5 we include abandonment in the model. From a numerical analysis, we show that reservation may be counterproductive if customers' impatience is high.

Stationary probabilities

In this section, we derive explicitly the stationary probabilities. The system is modeled using a two dimensional continuous-time Markov chain. We denote by (x, y) a state of the system for 0 ≤ x ≤ s and y ≥ 0, where x represents the number of busy servers and y represents the number of customers in the queue. The

state space S is S = {(0, 0), (0, 1), • • • , (0, s)} ∪ {s -c, s -c + 1, • • • , s} × N * .
We denote by p x,y the stationary probability to be in state (x, y) and by a the ratio λ/µ. We next describe the 4 possible transitions in the Markov chain.

1. An arrival with rate λ while a server is available (0 ≤ x < s), which changes the state to (x + 1, y).

The number of busy server is increased by one.

2. An arrival with rate rλ while all servers are busy (x = s), which changes the state to (x, y + 1). The number of customers in the queue is increased by 1.

3. A service completion with rate (s -c)µ while the queue is not empty (y > 0) and the number of busy servers is equal to s -c (x = s -c), which changes the state to (x, y -1). The number of customers in the queue is reduced by 1.

4.

A service completion with rate min(s, x)µ while the queue is empty (y = 0) or the queue is not empty but the number of busy servers is strictly higher than s -c (x > s -c, y > 0), which changes the state to (x -1, y). The number of busy servers is reduced by 1.

The Markov chain is depicted in Figure 1 of Section 1 of the online supplement. In Theorem 1, we give the stationary probabilities and the stability condition. The proof of Theorem 1 is given in Section 2 of the online supplement.

Theorem 1 Under the stability condition r

(s-c-1)! s! a c+1 < 1, we have p x,0 = a x x! p 0,0 , for 0 ≤ x ≤ s -c, (1)
p x,0 = a x x! 1 + r s-x-1 ∑ k=0 (x+k)! s! a s-k-x 1 + r c-1 ∑ k=0 (s-c+k)! s! a c-k p 0,0 , for s -c ≤ x ≤ s, ( 2
)
p x,y = a x   1 -a s-(x+c+1) [ 1 -r (s-c-1)! s! a c+1 ] s-x-1 ∑ k=0 (x+k)!a -k c ∑ k=0 (s-c-1+k)!a -k   rx! c ∑ k=0 (s-c-1+k)! s! a c+1-k     r c ∑ k=0 (s-c-1+k)! s! a c+1-k 1 + r c ∑ k=1 (s-c-1+k)! s! a c+1-k     y+1 p 0,0 , (3)
for s -c ≤ x ≤ s, y > 0, with

p 0,0 =     s-c-1 ∑ x=0 a x x! + s ∑ x=s-c a x x! 1 -r (s-c-1)! s! a c+1     -1 . ( 4
)

Impact on Balking

We evaluate in Proposition 1 the impact of the reservation threshold on the proportion of customers who balk, P .

Proposition 1 P is decreasing and convex in c.

The consequence of this result may be unexpected by service providers. By forcing some servers not to work when some work is available, these servers will at the end work more. The reason is due to the reduction of the flow of customers who balk. This result is stated in Corollary 1 and follows Proposition 1 using the relation: Servers' utilization = λ(1-P ) sµ .

Corollary 1 Servers' utilization is increasing and concave in c.

Proof. From Theorem 1, we deduce the explicit expression of

P from P = (1 -r) ∞ ∑ y=0
p s,y . This leads to

P = (1 -r) a s s! ( 1 -r (s-c-1)! s! a c+1 ) s-c-1 ∑ x=0 a x x! + s ∑ x=s-c a x x! . ( 5
)
The proportion of customers who balk, P , can be rewritten as

P = (1 -r) [ s ∑ x=0 s! x!a s-x -r s-c-1 ∑ x=0 (s -c -1)! x!a s-x-c-1 ] -1 = (1 -r) [ B(s, a) -1 -rB(s -c -1, a) -1 ] -1 ,
where B(t, a) =

( t ∑ x=0 t! x!a t-x ) -1
is the Erlang loss function for t ∈ N.

For analytic purposes, it will be convenient to extend the domain of the Erlang loss function to noninteger values. Following [START_REF] Jagers | Convexity of functions which are generalizations of the Erlang loss function and the Erlang delay function[END_REF], the continuous extension of the classic Erlang loss function is defined for any t ∈ [0, +∞) and a > 0 by

B(t, a) = (∫ +∞ 0 ae -ax (1 + x) t dx ) -1
As observed by [START_REF] Jagers | Convexity of functions which are generalizations of the Erlang loss function and the Erlang delay function[END_REF], the above integral expression coincides with the initial definition of B(t, a) for integer values of t.

This allows us to obtain the monotonicity properties by derivation. We have

∂P ∂c = -(1 -r)r ∫ +∞ 0 ae -ax ln(1 + x)(1 + x) s-c-1 dx [B(s, a) -1 -rB(s -c -1, a) -1 ] 2 < 0.
So, P is decreasing in c. We also have

∂ 2 P ∂c 2 = (1 -r)r ∂ 2 B(s-c-1,a) -1 ∂c 2 [ B(s, a) -1 -rB(s -c -1, a) -1 ] + 2r ( ∂B(s-c-1,a) -1 ∂c ) 2 [B(s, a) -1 -rB(s -c -1, a) -1 ] 3 . Since ∂ 2 B(s-c-1,a) -1 ∂c 2 = ∫ +∞ 0 ae -ax (ln(1 + x)) 2 (1 + x) s-c-1 dx > 0, we have ∂ 2 P ∂c 2 > 0 and P is convex in c. 2
In Figure 1, we represent the proportion of customers who balk as a function of the reservation threshold for different values of the system parameters. For each curve, we only evaluate P when the system is stable. At the end of each curve, for the highest possible value of the reservation threshold, we compute the relative improvement obtained in comparison with a situation without reservation, where the relative improvement is defined by relative improvement = P (c=cmax)-P (c=0)

P (c=0)
. This computation is done to show the maximal potential in applying a reservation policy. In addition to illustrate Proposition 1, we observe that the highest improvement is obtained for the highest value of r. In other words, the fewer customers are tempted to balk, the more balking can be reduced by a reservation strategy. This makes sense intuitively; for low values of r most customers balk, so only a small proportion of customers accept to wait. This means that the chance to find an available server is high at customer's arrival even without reservation. So, reservation has only a limited effect in this case. The opposite holds when r is high.

We also observe that we get higher improvement in larger systems. Recall that balking is due to situations where all servers are busy at a new customer's arrival. The idea of reservation is to create idle capacities to absorb the new arriving customers. From a situation where all servers are busy, the first service completion which allows the system to have a first idle server occurs after an expected duration of 1/(sµ). This duration reduces with the system size. So, large systems have a stronger ability to recreate idle capacity which in turn leads to fewer customers who balk.

Consequences on the waiting time

We provide here a method to compute the Laplace Stieltjes transform (LST) of the waiting time distribution of an arriving customer in order to evaluate the consequences of reservation on the waiting time.

Recall that an arriving customer waits only if all servers are busy. We focus then on an arriving customer in state (s, y), signifying that there are s busy servers and y other customers in the queue (y ≥ 0). The state of the system is denoted by (x, y) where x is the number of busy servers and y is the number of customers in front of the considered customer.

The state transition diagram for the finite two-dimensional birth-and-death process modeling the behavior of a tagged customer in the queue is shown in Figure 2. This process has (c + 1) × (y + 1) transient states and one absorbing state called "Service". Let W x,y (t) be the unconditional LST of the density function for

ݏμ ሺݏ െ 1ሻμ ሺݏ െ 2ሻμ ሺݏ െ ܿ 1ሻμ λ λ λ λ ሺݏ െ ܿሻμ ሺ,ݏ 0ሻ ሺݏ െ 1,0ሻ ሺݏ െ 2,0ሻ ሺݏ െ ܿ, 0ሻ ݏμ ሺݏ െ 1ሻμ ሺݏ െ 2ሻμ ሺݏ െ ܿ 1ሻμ λ λ λ λ ሺݏ െ ܿሻμ ሺ,ݏ ݕሻ ሺݏ െ 1, ݕሻ ሺݏ െ 2, ݕሻ ሺݏ െ ܿ, ݕሻ ݏμ ሺݏ െ 1ሻμ ሺݏ െ 2ሻμ ሺݏ െ ܿ 1ሻμ λ λ λ λ ሺݏ െ ܿሻμ ሺ,ݏ ݕ െ 1ሻ ሺݏ െ 1, ݕ െ 1ሻ ሺݏ െ 2, ݕ െ 1ሻ ሺݏ െ ܿ, ݕ െ 1ሻ ^Ğƌ|ŝĐĞ Figure 2: Transition diagram
the passage time from state (x, y) to state "Service". In Theorem 2, we provide a method which allows us to obtain the expression of W s,y (t) after a finite number of computations. The proof is given in Section 3 of the online supplement.

Theorem 2 For y ≥ 0, we have

W s,y (t) = 1 f c+1 (t) ( f c (t) f c+1 (t) ) y , ( 6
)
where the functions f k (t) for 0 ≤ k ≤ c + 1 are given by f 0 (t) = 1, f 1 (t) = 1 + t sµ and

f k (t) = ( 1 + λ + t (s -(k -1))µ ) f k-1 (t) - λ (s -(k -1))µ f k-2 (t), (7)
for 2 ≤ k ≤ c + 1.
In Corollary 2, we give the expected value and the variance of the waiting time for an arriving customer in state (s, y) for y ≥ 0; E(W y ) and V (W y ). The proof is given in Section 3 of the online supplement.

Corollary 2

We have

E(W y ) = 1 µ ( (y + 1) c ∑ k=0 k ∑ i=0 a i (s -k -1)! (s -k + i)! -y c-1 ∑ k=0 k ∑ i=0 a i (s -k -1)! (s -k + i)! ) , V (W y ) = 1 µ 2   2(y + 1) ( c ∑ k=0 k ∑ i=0 a i (s -k -1)! (s -k + i)! ) 2 -(2y + 1) ( c ∑ k=0 k ∑ i=0 a i (s -k -1)! (s -k + i)! ) ( c-1 ∑ k=0 k ∑ i=0 a i (s -k -1)! (s -k + i)! ) +y   c ∑ k=1 k-1 ∑ i=1 i-1 ∑ j=0 j ∑ m=0 a m+k-i+1 (s -k)!(s -j -1)! (s -1)!(s -j + m)!   -(y + 1)   c+1 ∑ k=1 k-1 ∑ i=1 i-1 ∑ j=0 j ∑ m=0 a m+k-i+1 (s -k)!(s -j -1)! (s -1)!(s -j + m)!    
From Corollary 2, one can compute the first and the second moments of the waiting time for an arbitrary customer who do not balk by considering the system at arrival instants using the stationary probabilities in Theorem 1. For the expected waiting time, E(W ), we get

E(W ) = r µ c ∑ k=0 (s-c-1+k)! s! a c-k s ∑ x=s-c a x x! - ( 1 -r (s-c-1)! s! a c+1 ) s ∑ x=s-c s-x-1 ∑ k=0 (x+k)! s!x! a s-1-k ( 1 -r (s-c-1)! s! a c+1 ) ( s-1 ∑ x=0 a x x! -r (s-c-1)! s! a c+1 s-c-2 ∑ x=0 a x x! ) . ( 8
)
Another way to obtain E(W ) is to apply Little's law by first computing the expected number of customers in the queue and next by dividing this quantity by λ(1 -P ). The complexity of E(W ) does not allow us to prove the first and second order monotonicity results in c. Yet, from an extensive numerical study, it seems clear the E(W ) is increasing and convex in c (see Figure 3).

Reservation produces two phenomena which may degrade the waiting performance measures. First, as shown in Section 3, reservation reduces the proportion of customers who balk and therefore increases the workload. This is somehow positive since this also allows the service provider to serve more customers. Yet, it is clear that a higher workload will result in longer waiting times for customers who wait. Second, with reservation one may encounter situations with one idle server and a customer waiting in the queue, this may also have a negative effect on the waiting performance measures. To isolate this second phenomenon, we present in Figure 4 different waiting time performance measures in the case r = 1 (no balking); the The first observation in Figure 4(a) is that the probability of delay reduces with reservation. This is particularly apparent in high workload situations. This is already proven in Proposition 1 since the probability of delay is proportional with the proportion of customers who balk; P = (1 -r)P D . The second is that the expected value and the standard deviation of the waiting time increase with reservation (Figures 4(b) and 4(c)). As mentioned above, the reason is related to the loss of work-conservation. It is however interesting to note that the variability of the waiting time measured by the coefficient of variability reduces with reservation in high workload situations. This is mainly due to the high proportion of customers who do not wait.

Benefits of reservation with abandonment?

In Section 4 of the main document, we have shown the degradation of the waiting time with reservation.

Given that high waiting times may induce a high abandonment, one may wonder if reservation instead of allowing to serve more customers would instead lead to more customers' loss. We therefore propose here to include abandonment in the modeling. For this purpose, we consider an exponential abandonment with parameter β. The complexity added by abandonment does not allow us to obtain simple closed-form expressions for the stationary probabilities. However, in Section 4 of the online supplement, we propose a numerical method to derive them.

Figure 1 :

 1 Figure 1: P as a function of c (µ = 1, λ = sµ)

Figure 3 :Figure 4 :

 34 Figure 3: E(W ) as a function of c (µ = 1, s = 40, r = 50%)

Figure 5 :

 5 Figure 5: Impact of abandonment (s 5, λ = 8, µ = 1, r = 0.5)

Acknowledgment. The author wants to express his gratitude to the reviewer and the associated editor for their useful comments, that significantly improved this paper. The author would like also to thank Sébastien Thorel from INTERACTIV GROUP for our helpful discussions.