Modeling of the chip geometry in orbital drilling

Keywords: Orbital drilling, TiAl6V4, tool geometry, cutting forces

The orbital drilling is a complex operation. Due to the tool trajectory, which is helical, chip thickness is highly variable. This is why the cutting forces are very difficult to estimate. The aim of this study is to develop a model to predict cutting forces depending on the tool geometry and cutting conditions in order to control the final quality of the machined hole. At first, the geometry of the chip is modeled taking into account the parameters defining the trajectory and the tool. A cutting force model based on the instantaneous chip thickness is then set up. An experimental part validates the cutting force model through measures of cutting force made during orbital drilling tests.

Introduction

The orbital drilling, also called helical milling, is a process notably used in the aeronautical industry. It is very different from the axial drilling. It consists in machining a hole with a tool which has a smaller diameter, driven on a helical trajectory (Figure 1). This process is used especially because it generates low cutting forces [Lutze, 2008]. As a consequence, the burr formation at the entry and at the exit of the hole is largely reduced, and also the risk of delamination when drilling composite materials is highly reduced. For the cutting force model in orbital drilling, three approaches are possible: the numerical approach, the analytical and the semi analytical approach. The numerical approach used in the work of Lorong [Lorong et al, 2002] aims to simulate chip formation by a finite element study. Analytics approaches as Merchant's model Merchant, 1945] or Oxley 's model [Oxley, 1989] consist to describe the most physically possible the cutting phenomenon in order to calculate the cutting forces. The semi-analytics approaches, also called mechanistic approaches, are the most used one to calculate machining cutting forces and particularly when the process is complex [Fontaine, 2004], like helical milling. This approach is based on analytics approaches by simplifying all the parameters of the tool geometry or of the tool/workpiece couple by a constant experimentally defined with calibration test.

Modeling of orbital drilling

2.1.

Geometrics settings of the orbital drilling

In order to be able to develop a model of the chip geometry, the input data of the model and the references that are necessary to form the equation, need to be previously defined.

1.1.1

The defining data of the tool's geometry:

The global tool geometry (Figure 2) is generically defined using the following parameters:

-Tool outside radius: ( )

With: |

The angle defines the tool angular position considered in the machine reference. The tool center position in position « i » is called CL i , at a distance R off to HL.

Tool reference: This reference can get the location of every points of the tool « i » in the orbital reference.

( )

With: |
The angle defines the angular position of considerate point of the tool in the orbital reference

Representative function of cutting edge profile

In order to calculate the chips section, it's necessary to define the cutting edge profile. The function H k (r) is established, to give the elevation of each cutting edge point versus radius r.

The lowest tool point « k » is defined at the elevation Z = 0. This function is decomposed into several successive functions, to take into consideration the different tool geometries.

Examples, the evolution of the function H k (r) for a tooth without center cut, is represented Figure 5.

 If :
( ) d is a value arbitrarily set to ensure the non-participation of the area during the machining.

 If : ( ) [ ( )] ( )  If : ( ) √ [ ( )]
κ r Therefore, to calculate the chip section machined by a tool tooth « i » at the present time (the tool identified by its tool center's point CL i ) two calculation steps are made: First, each point belonging to the envelope surface of the tool « i » is detected in from its coordinates

Tooth profile with center cut Tooth profile without center cut

( ) in function of et ( [ ], [
]), and the function applied at the tool « i » :

{ ( ) ( ) ( ) ( )
The set of the points for a constant , represent the tooth where the chip section is assessed on.

Then, to identify the surface location previously machined on each point , all the locations of the tool that interve before the location « i » on a complete orbit revolution are considered. It will then be possible to estimate the height of the machined material at the point . Each considered previous location is identified by its tool center point CL k.

The difference of altitude following between two tool locations CL i and CL k is called . Knowing that the tool describes a helical trajectory characterized by its pitch P (Figure 3), the height is defined according to the relative angular position of CL i and CL k on the trajectory characterized by (Figure 6):

( )
Then, for each tool point considered, the set of the point s of which the coordinates in the locator are such that ( ) are identified. The coordinate of the point can be calculated through the function :

{ ( ) ( ) ( ) ( ) ( )
For this the radius has to be previously determined. From Figure 6, by considering the triangle( The resulting force is then:

Experimental procedures

All tests were performed in a test piece in TiAl6V4 of 19mm thick. The cutting conditions are as follows: Vc= 30mm/min; Fza=0.005mm/tooth; Fzt= 0.04mm/tooth. The tool used is a tool with three teeth (Z=3) ,with a single tooth with center cut and of diameter D t = 9mm.

3.1.

Identification of specific cutting coefficient

To identify the specific cutting coefficient Kc, tests of slot machining were carried out on a machining center. The force measurement is performed using a Kistler dynamometer. The cutting coefficient, assumed to be constant, can be identified from the measurement of the maximum cutting force. The chip thickness is 0.04mm per tooth, the depth of cut is 0.75mm, and so the maximum chip section is 0.03mm ².

It can be deduced:

Maximum force Fxy observed during the test is 105N. This gives a specific cutting coefficient for our tool/work piece couple of: Kc = 3500N/mm ². This Kc value is a first approximation for model validation. It doesn't take into account influence of chip thickness or possible influence of speeds. And thereafter, two coefficient will be identified, one for the axial part and one for radial part, because the radial and axial parts present different cutting angles.

Force measurement in orbital drilling

Tests are also conducted in orbital drilling, in order to have a database to validate the modeling. These tests are performed on a drill bench equipped with an orbital spindle. On all tests, cutting forces are recorded using the dynamometer Kistler and sampled at 10 kHz. On the signals obtained when measuring just a taring is performed. The hole diameter is 11.1mm and it is made with the same tool which has diameter of 9mm.

Results and discussions

The Kistler 9257B dynamometer allows measuring the resultant force in the plane normal to the axis of the tool Fxy With √ (Figure 11) The force F cT obtained by modeling differs from the real effort measured Fxy (Figure 11)

The evolution of the modeled force on a tool revolution shows the successive passage of three teeth. This is explained by the fact that the radial chip is not constant on the tool revolution. In addition there is a single tooth with center cut, therefore causing disruption on the tool revolution. The force measurement Fxy shows only two peaks of effort. Either the force is compensated perfectly on the third revolution, or either the chip is not constant for all three teeth. This difference between the two graphs is the lack of taking into consideration of cutting phenomena. To improve the model, the tool will be divided into two different parts, in order to better understand how each part works. The axial and radial part of the tool is defined as (Figure 12).

The tool performance and the cutting phenomena caused by the radial part of the tool are known and already considerably studied in the literature. Modeling the cutting force of the axial part is more complex because the tool profile of the tip (Figure 12) combined with the tool trajectory is not classical. In order to better understand the mechanisms during drilling, the evolution of the tangential force F t at the contact point tool / work piece (according to Y oi ) and of the radial force F r (Figure 13) are modeled. These two efforts strongly oscillate on a tool revolution, causing a variable bending force of the latter. So a dynamic phenomenon is established and must be taken into account in the model to adjust the chip thickness of the radial cut. The Figure 13 shows that when the tooth with center cut is located in the "A" area ( ), the radial force and tangential force are positive but quite low, not causing a significant modification in radial section. When the tooth with center cut (z 1 ) is located in the "B" area (

), the radial force and tangential force are positive. Therefore this resultant force tends to push the tooth being machined towards the surface. When the tooth with center cut is located in the area "C", the radial force is positive and tangential force is negative. Therefore this resultant force tends to push the tooth (z 2 )

Tool revolution (in degrees) And when the tooth with center cut is located in the area "D", the radial and tangential force are negatives. Therefore this resultant force tends to withdraw the tooth (z3) of hole surface.

One tool revolution

Figure 13 Modeling of the radial and tangential force of the tool axial part, taking into account the three teeth.

The hypothesis is that the radial section isn't homogeneous for the three teeth, because bending forces are different to each passage of various teeth on the radial cut area. The radial chip thickness machined by the tooth z 2 is greater than that machined by the tooth z 3 . After integration of this dynamic phenomenon in the model, the force profile is similar to the force measure by its evolution and maximum values (Figure 14). A further difference is found between the model and the measure (denoted "B" in Figure 14). In these areas, modeling tends to decrease cutting force versus measured force. This difference is explained by the hypothesis made on the "Kc" at baseline: the cutting coefficient is not constant and tends to increase in areas where the chip thickness is very low. This evolution of cutting coefficient will be therefore taken into account in order to improve the model.

Conclusions

This model, with the integration of the tool geometry, allowed knowing the exact geometry of the chip. With this knowledge, the modeling of the cutting force is possible which allows a better understanding of the cut phenomena when orbital drilling with a specific tool. It is therefore possible to vary the tool geometry or the cutting conditions in the model to predict the cutting forces, and thus be able to optimize these parameters. If the link between effort and drilling defects are known (ex: bending force → decrease in diameter), it is possible to optimize these parameters to improve the quality of the hole.

  [Denkena and al., 2008] [Brinksmeier et al., 2008]. Operations such as cleaning or deburring are then considerably reduced.

-

  Figure 1 Orbital drilling's trajectory

Figure

  Figure 2 : Tool geometry

  Figure 5 : Representation of the cutting edge profile « k » without cutting tool center

Figure 12 :

 12 Figure 11: Comparison of the evolution of the cutting force over three revolution of tools (modeling at top; measure down)

Figure 14 :

 14 Figure 14: Modeling the cutting force from the corrected model.

Copyright of IDMME -Virtual Concept -Improve -Ingegraf

Acknowledgements

This work was carried out within the context of the working group Manufacturing'21 which gathers 18 French research laboratories. We would also like to thank the project called OPOSAP (Optimisiation du Perçage orbital avec Surveillance Active du Process) and above all the partners who helped us to make this research.

Copyright of IDMME -Virtual Concept -Improve -Ingegraf

The chip section on a tooth is then obtained by plotting the set of points and the associated points , for a fixed angle (representing a tooth of the tool) and for [ ] (Figure 7). This section is then split into trapezoids of thickness and width . The thickness is the normal distance to the tool profile between the profile and the previously machined surface , it represents the chip thickness. The width is the distance between points , it is set arbitrarily.

If the decomposition is sufficiently thin (small width compared to the radius of the tool), each trapezoid can be approximated and simplified by a rectangle. The chip section can then be written:

Calculation of the cutting force

The chip section in each moment being determined, the forces applied on the tool during drilling can be modeled.

In this article, only the modeling of the cutting force normal to the chip section is presented and discussed. This effort is the most important and is the main cause of tool bending.

The other two components can be modeled in the same way.

The modeling approach is semi-analytical [Merchant, 1945], which gives for the same tool/work piece couple:

Where is a specific cutting coefficient that will remain constant for the first part of the study and will be identified with calibration tests. This model is applied to each trapezoid of the chip section and then summed to obtain the total cutting force when the tooth is in position .

The evolution of the cutting force of a tooth on a tool revolution is calculated by varying from 0 to 2π (Figure 9). This cutting force ( ) must be calculated for each tooth of the tool, taking into account the differences in the tooth profile which may exist (ex: with or without center cut).

To calculate the resulting force, every effort is expressed in the orbital reference , so (Figure 10):

Cutting force evolution (in N)


Yves Davila
email: yves.davila@iut-tlse3.fr

Laurent Crouzeix

Bernard Douchin

Francis Collombet

Yves-Henri Grunevald

Identification and modelling of the in-plane reinforcement orientation variations in a CFRP laminate produced by manual lay-up

Keywords: variability, in-plane fibre orientation, FE models, meta-models

Reinforcement angle orientation has a significant effect on the mechanical properties of composite materials.

This work presents a methodology to introduce variable reinforcement angles into finite element (FE) models of composite structures. The study of reinforcement orientation variations uses meta-models to identify and control a continuous variation across the composite ply. First, the reinforcement angle is measured through image analysis techniques of the composite plies during the lay-up phase. Image analysis results show that variations in the mean ply orientations are between -0.5 and 0.5° with standard deviations ranging between 0.34 and 0.41°. An automatic post-treatment of the images determines the global and local angle variations yielding good agreements visually and numerically between the analysed images and the identified parameters.

A composite plate analysed at the end of the cooling phase is presented as a case of study. Here, the variation in residual strains induced by the variability in the reinforcement orientation are up to 28 % of the strain field of the homogeneous FE model. The proposed methodology has shown its capabilities to introduce material and geometrical variability into FE analysis of layered composite structures.

1.

Introduction

Composite structures show a rather large inherent dispersion in their mechanical properties compared to metallic structures. The scatter in the mechanical properties on composite structures has multiple origins mainly related to the geometry and manufacturing. Yet, composite structures are typically analysed through finite element (FE) models without including these variations (i.e. the calculations done assuming symmetric lay-ups, equal ply thicknesses and constant reinforcement orientation, constant volume fractions, etc.).

To introduce local variations (point to point and ply per ply) into FE models that do not only have the same order of magnitude, but exhibit similar variation trends, an original methodology has been proposed [START_REF] Davila | Étude multi-échelle du couplage matériaux-procédés pour l'identification et la modélisation des variabilités au sein d'une structure composite[END_REF]. This methodology uses meta-models that control the continuous spatial variations of the material properties at each element of the FE model as observed in real composite structures. To evaluate a broad set of scenarios, the measured property variations are not introduced into the FE models 'as is'. Instead, variables in the metamodels are described through their probability distribution functions (PDF). These PDF are obtained through the observation, measurement and analysis of the composite part before, during and after its fabrication. The present methodology was first used to assess thickness variations in a laminated composite plate [START_REF] Davila | Spatial Evolution of the Thickness Variations over a CFRP Laminated Structure[END_REF]. The evaluation of a second source of variability, variations in the reinforcement orientation, is hereby presented.

Variations in fibre orientation have an appreciable impact on the performance of a composite structure [START_REF] Chamis | Probabilistic Composite Design[END_REF]. These variations, usually perceived as fibre misalignments. Fibre misalignments are present at different scales of the composite structure and are the result of several factors during the different manufacturing stages of the structure, and are present. A global misalignment is considered when the mean orientation of the ply, or fibre bundles, differs from the intended ply orientation. Local misalignments can be produced by the periodic undulation of the prepreg fibres and dry performs [START_REF] Potter | Variability, fibre waviness and misalignment in the determination of the properties of composite materials and structures[END_REF] and by the preform manipulation during the lay-up phase of the composite structure [START_REF] Trask | Influence of process induced defects on the failure of composite T-joint specimens[END_REF].

The reinforcement misalignments and their effects in unidirectional composite parts have been subject of different studies. In these studies misalignment are determined by both indirect and direct measurement methods. Arao [START_REF] Arao | Effect of ply angle misalignment on out-ofplane deformation of symmetrical cross-ply CFRP laminates: Accuracy of the ply angle alignment[END_REF] obtained the fibre misalignments indirectly by comparing the numerical analysis to the experimental results composite plate paced into a moist environment obtaining a normal distribution of the fibre misalignments. For this study, the misalignments were found normally distributed with a mean of 0° and a standard deviation of 0.4°. However this method only describes overall fibre misalignments without producing any information on local misalignments. Yurgartis [START_REF] Yurgartis | Measurement of small angle fiber misalignments in continuous fiber composites[END_REF] proposed a method to determine fibre orientations by examining the cross section of a unidirectional and 0/90 laminates. The fibre orientation is obtained by measuring the resulting ellipsoidal form of the cross-section of the reinforcement. This method allows the determination of both in-plane and out-ofplane local fibre misalignments. The main inconvenient of this technique is that the variation in the fibre orientations must be calculated through the mean of several hundred to thousands fibres [START_REF] Barwick | Identification of fiber misalignment in continuous fiber composites[END_REF]. The measured zones are in the range of 0.015 mm² to 100 mm², the latter after the assembly of 6 624 micrographs. Automatic image treatment algorithms are introduced to reduce the number of images needed and the analysis time [START_REF] Creighton | A multiple field image analysis procedure for characterisation of fibre alignment in composites[END_REF].

The use of a fast Fourier transform (FFT) coupled with an autocorrelation function and filtering techniques permits the analysis of images taken with lower quality and resolution [START_REF] Kratmann | A novel image analysis procedure for measuring fibre misalignment in unidirectional fibre composites[END_REF]. This allows the analysis of images obtained by microcomputed tomography (µCT) [START_REF] Sutcliffe | Measurement of fibre waviness in industrial composite components[END_REF] . Other type of µCT technique, is to identify, at voxel size, an isolated fibre along its axis following its 3D reconstruction, of and obtaining its misalignments respect all major planes [START_REF] Requena | 3D-Quantification of the distribution of continuous fibres in unidirectionally reinforced composites[END_REF].

However, these techniques only allow a post-mortem characterisation of the fibre orientations with specimens that need to be cut and prepared for micrograph inspection. Taking into account factors such the image size, image resolution and the number of images needed to characterise a millimetre size specimen, it is evident that a broad study becomes time and resource intensive.

The identification of reinforcement misalignments can be also done on the prepreg before the curing phase. Skordos [START_REF] Skordos | Stochastic simulation of woven composites forming[END_REF] and by Mesogitis [START_REF] Mesogitis | Non-crip fabrics geometrical variability and its influence on composite cure[END_REF] used a mounted a digital camera placed on a coordinate measuring machine to take a set of photos of a woven prepreg. 630 images were necessary to cover a surface of 103 500 mm² (300 x 345 mm). Each image is processed using a FFT coupled with an autocorrelation function to determine the angle of the warp and weft orientation and as well as the length of each unit cell.

Potter et al. [START_REF] Potter | Variability, fibre waviness and misalignment in the determination of the properties of composite materials and structures[END_REF], [START_REF] Potter | Sources of variability in uncured aerospace grade unidirectional carbon fibre epoxy preimpregnate[END_REF] identified that the reinforcement of UD prepreg present an undulation in form of a sinusoidal wave. This effect is due that the composite prepreg is stored in rolls. The identified sin wave has a mean amplitude of 0.03 mm and a wavelength of 3 mm. The present study of the reinforcement orientation and its impact on the structural behaviour of a composite plate here presented is divided in two main phases. The first phase concerns the identification and measurement of actual fibre orientation at different composite scales (local and complete ply) by means of optical image analysis. The second phase is the introduction of a representative set of values of the previously identified trends into a FE model that assess the structural behaviour of a composite plate with variable reinforcement angles.

Materials and fabrication process description

Six composite plates are produced by manually stacking the unidirectional carbon/epoxy prepreg HexPly ® M10.1/38%/UD300/CHS, furnished by Hexcel Composites. The dimensions of the manufactured plates are 615 mm long and 300 mm wide. The 16 ply stratification is chosen to be quasi-isotropic. The lay-up sequence is [90/-45/0/(+45)2/0/-45/90]s.

The fabrication of the composite plate starts with cutting the plies from the prepreg roll into different shapes to fit different ply orientations. The plies are cut following a pattern designed to reduce material waste.

Due to the staking sequence, some plies need to be assembled from two or more shapes so that fibres run continuously in the desired angle. The lay-up is done manually by two operators. To reduce the variability induced by the operators, all composite plies were laid by the same operators using the same technique. For each composite plate, after the placement of the first eight plies the uncured preform is moved to a debulking area. This step is needed to ensure that the entrapped air between plies is removed to reduce inter laminar porosity in the cured plate. To continue the lay-up of the last eight plies, the repositioning of the uncured composite plate at the exact same position in the working zone cannot be assured. For this reason, the lay-up must be self-referenced in any moment during fabrication. The first ply (oriented at 90°) is 10 mm larger than the rest of the plies, being 625 long and 310 mm wide. Plies # 2 to # 16 are thus referenced to ply # 1.

The six plates, labelled C-11, C-12, C-13, C-21, C-22 and C-23, are cured in the autoclave using the recommended cure cycle, with a consolidation dwell at 90°C with a pressure of 2 bar for 15 minutes and a cure dwell at 120°C with a pressure of 5 bar for 60 minutes. The heating and cooling rates are 2 °C min -1 .

Optical analysis and image treatment

Experimental setup

During the lay-up process, a series of images are taken using a Canon ® EOS 550D digital camera with an EF 50MM/macro lens. The digital camera is mounted 1.7 m above the working area (cf. Figure 1). The goal of this setup is to take a single image of each complete prepreg ply. This eliminates the need to move the camera across the composite surface and eliminate errors due to image stitching during post-treatment. The working area is comprised of a reference grid of 700 x 400 mm made to place the prepreg lay-up (cf. Figure 2). This reference grid serves to three purposes. The first is to correlate the pixel size to the actual dimension of the composite plies. The second is to control the image distortion. And the third is to determine the position of punctual defects in the lay-up.

Barrel and perspective distortions on the acquired images are the result of the positioning and alignment of the digital camera relative to the grid, as well as the optics of the camera and lens. These image distortions were corrected using the commercial software DxO Optics Pro ® .

Description of the Image treatment

An in-house automatic measurement method is developed to determine the fibre orientation angles.

This procedure is based on the methodology for line detection of Shi and Wu [START_REF] Shi | Automatic Fiber Orientation Detection for Sewed Carbon Fibers[END_REF] and uses a Hough transform on the fibre bundles at the surface of the prepreg to determine the fibre bundle angle. The different steps for obtaining the reinforcement orientation in the acquired images are illustrated in Figure 3. The image is cropped to exclude the ply edges because the latter disrupt the cut-off values of pixel intensity in the image treatment (cf. Figure 3a and Figure 3b). The images are segmented into a predefined number of elements that can be treated individually (cf. Figure 3b). To enhance the lines formed by the fibre bundles a Laplacian of Gaussian (LoG) kernel is used. The LoG filter is coupled with a series of convolution masks (Robinson type) [START_REF] Redon | Automatic image analysis and morphology of fibre reinforced concrete[END_REF] to refine the filtering of the image (cf. Figure 3c). An edge detection function [START_REF] Canny | A computational approach to edge detection[END_REF] is applied to highlight the lines formed by the fibre bundles (cf. Figure 3d). Afterwards, the Hough transform is performed to extract the angle information of the treated fibre bundles (cf. Figure 3e). Each 5184 x 3456 pixel image has a resolution of 6.48 pixel/mm with a pixel aspect ratio of 1. An element size of 30 mm per side (in a grid of 20 by 10 elements, thus 200 elements in total) is chosen as good compromise between measurement accuracy and processing time. Figure 4a shows that an element of 30 mm per side contains the largest amount of Hough lines. For all element sizes, the average of the Hough lines orientation is less than 0.4° (cf. Figure 4b). However, the standard deviation of the Hough lines orientations increases drastically when the size of the element is 20 mm per side or lower. The processing time for each ply, with an element size of 30 mm per side is approximately 40 s using a 64 bits Intel ® Core TM i3 processor running at 2.40 GHz. Hough lines with angles not within ± 7° tolerance are removed as they are considered as artefacts of the image treatment. The mean angle of the Hough lines is calculated for each element (image segment) to map local orientations. Figure 5 shows orientation maps for 4 different plies in the orientations of 0°, +45°, -45° and 90° respectively. Certain elements in the plies oriented at 0° and 90°are discarded since none of the Hough lines are within the ± 7° tolerance. In average, there are 48 discarded elements for the 0° angle plies and 19 discarded elements for the 90° angle plies. There are not discarded elements for the ± 45° orientation plies. The main difference of the proposed methodology with the one proposed by Skordos [START_REF] Skordos | Stochastic simulation of woven composites forming[END_REF], is the method used for capturing the images and the post-processing techniques to obtain the reinforcement orientation. In the proposed method, the complete surface of the prepreg is taken in a single image. The image is corrected for image distortions produced by the optical lens. By having one image, a positioning device for the camera is not necessary. Another difference is in the use of an autocorrelation function coupled with a fast Fourier transform (FFT). FFT algorithms best operate on square images with an area being a power of two, thus limiting the size of the image. The proposed method here uses a Hough transform of a filtered image, the image size is not restricted.

From the orientation maps, variations in the reinforcement orientation across the ply can be noticed with the naked eye. These variations appear to be non-periodic and non-repetitive through the different plies.

An identification algorithm is developed to quantify the variation trends in the local fibre orientations.

Fibre in-plane misalignments meta-models

It is assumed that the fibres do not move during the curing process. Although, this assumption is not completely true since there is resin flow and evacuation of entrapped air during the first temperature ramp of the curing cycle where the viscosity of the resin is drastically reduced, for practical reasons, it is considered that the shift in the reinforcement position is negligible. Further studies are required to quantify the shift in the fibre position before and after curing.

A number of equations describe the changes in the reinforcement angle for an individual element taking into account a continuity in the values from one element to another (i.e. a maximum value cannot be beside a minimum value).

The reinforcement orientation at any point θk (x,y) of the composite ply can considered as the result of three different types of orientations (cf. eq. 1): the mean ply orientation 𝜃 𝑘 ̅̅̅̅ , the periodic fibre misalignments θwav and the non-periodic perturbations θpert.

𝜃 𝑘 (𝑥, 𝑦) = 𝜃 𝑘 ̅̅̅̅ + 𝜃 𝑤𝑎𝑣 (𝑥, 𝑦, 𝜃 𝑘 ̅̅̅̅ ) + 𝜃 𝑝𝑒𝑟𝑡 (𝑥, 𝑦)

The mean ply orientation 𝜃 𝑘 ̅̅̅̅ contains the theoretical ply orientation θth k (stratification) and an overall ply misalignment δθk (cf. eq. 2).

𝜃 𝑘 ̅̅̅̅ = 𝜃 𝑡ℎ 𝑘 + 𝛿𝜃 𝑘 (2)
The local continuous in-plane undulations θwav are modelled through the form of a sinusoidal wave in the direction of the ply orientation (cf. eq. 3). The mean amplitude Awav of 0.03 mm and the mean wavelength λwav of 3 mm, with extreme values between 2.2 and 4 mm, are obtained from the literature [START_REF] Potter | Variability, fibre waviness and misalignment in the determination of the properties of composite materials and structures[END_REF], [START_REF] Potter | Sources of variability in uncured aerospace grade unidirectional carbon fibre epoxy preimpregnate[END_REF].

                   wav k th k k th k wav k th k wav y x A y x           sin 2 cos 2 sin 2 , , (3)
It is observed that the in-ply fibre orientations θpert present a non-periodic continuous variations result from the prepreg manipulation during the lay-up. These local variations of the reinforcement angle can be described by the sum of N pseudo-Gaussian surfaces or perturbations (cf. eq. 4 and Figure 6). This methodology is inspired by Goshtasby and Oneill for 2D curve fittings [START_REF] Goshtasby | Curve Fitting by a Sum of Gaussians[END_REF], which was adapted for a 3D

where Xi and Yi are the coordinates of the central point of the perturbation in mm, αi and βi are respectively the length and the width of each i th pseudo-Gaussian surface in mm, and are homologues to the standard deviation in a normal distribution. Finally, Bi is the amplitude of the perturbation in degrees.

Results analysis of the optical measurements

Mean ply misalignment (ply scale)

In Table 1 are indicated the values of the difference between the mean ply orientation and its intended orientation δθk, with respect to ply # 1 for all plates. These δθk values were confirmed by manual measuring of the ply orientation [START_REF] Davila | Quantification of sources of variability in CFRP plates cured in autoclave[END_REF]. For plates C-11, C-12, C-13, C-21 and C-22, δθk values range between -0.5° and 0.5°.

For plate C-23, δθk values range from 0.4° to 1.1°. The standard deviation of the ply misalignment for all plates ranges from 0.30° to 0.40°, the overall misalignment is shifted towards a positive angle. It is assumed that this shift has its origins in the lay-up procedure with the same operators. No correlation was found between the intended ply orientation and its measured misalignment.

1 90° - - - - - - 2 -45° 0,7° 0,3° 0,0° 0,1° 0,1° 0,4° 3 0° 0,1° 0,4° -0,2° -0,1° 0,0° 0,8° 4 +45° 0,3° 0,0° -0,3° -0,1° -0,2° 0,4° 5 +45° -0,1° 0,2° -0,3° 0,1° -0,1° 0,6° 6 0° 0,1° 0,2° -0,1° 0,1° -0,1° 0,5° 7 -45° -0,1° 0,4° 0,1° 0,0° 0,3° 0,8° 8 90° -0,3° 0,1° 0,0° 0,1° 0,0° 0,5° 9 90° -0,2° -0,1° -0,2° 0,3° -0,4° 0,7° 10 -45° 0,0° 0,2° 0,3° 0,5° 0,0° 1,0° 11 0° 0,0° 0,4° 0,2° 0,3° -0,1° 1,1° 12 +45° -0,1° -0,2° 0,1° 0,0° -0,4° 0,5° 13 +45° 0,0° -0,1° -0,1° -0,1° -0,4° 0,6° 14 0° 0,4° 0,2° -0,1° 0,3° -0,2° 0,9° 15 -45° 0,4° 0,4° 0,3° 0,6° 0,0° 1,1° 16 90° 0,2° 0,1° 0,2° 0,3° -0,3° 0,8° Mean 0,1° 0,2° 0,0° 0,2° -0,1° 0,7°
The average in the global differences between the intended orientation and the measured fibre directions is inferior to 1.0° with a standard deviation (SD) of 0.35°. The measured values are lower than the values reported by Yurgartis [START_REF] Yurgartis | Measurement of small angle fiber misalignments in continuous fiber composites[END_REF] (mean angle of 3.7° and SD of 0.79°) and Olave [START_REF] Olave | Internal Geometry Variability of Two Woven Composites and Related Variability of the Stiffness[END_REF] (mean angle of 1.7° and SD of 0.69°). However, they are in good agreement with the values obtained by Skordos [START_REF] Skordos | Stochastic simulation of woven composites forming[END_REF] (mean angle of -0.1° and SD of 0.36°). The values of the mean and the SD are characteristic and unique to the studied composite part since they are function of the reinforcement nature and the actual manufacturing techniques used to conform the piece. However, the identification of variation trends are useful to construct and enlarge a database of the different property variations linked to each manufacture steps.

The value of the mean ply misalignment is lower than the initially expected value. In this case, the small deviation is partially due to the size of the ply, where a very noticeable misplacement of 10 mm at one of the corners of the 600 mm long ply is roughly translated into a 1.0° of ply misalignment. Nevertheless, global misalignments can be affected by the shape, size and cutting precision of the prepreg plies. The combination of these factors can conceal significant variations in the ply positioning up to 5° [START_REF] Torres | Valeur ajoutée d 'une encapsulation de capteurs silicium pour l'instrumentation à coeur des structures composites[END_REF]. Coefficient of linear thermal expansion the composite in the transverse direction α2 (°C -1 ) 14.7•10 -6

Figure 8 shows the strain field εxx of the upper skin. The deterministic solution strain field εxx, without taking account any variations, shows a uniform value of -83.4•10 -6 mm/mm (marked as a reference value in Figure 8), whereas the model accounting for variation in the reinforcement orientation shows compressive strains values ranging from -74.4•10 -6 up to -98.6•10 -6 mm/mm. This represents a variation up to 28 % of the reference value determined for the deterministic solution, which is significant.

Figure 1 .

 1 Figure 1. DSLR camera mounting.

Figure 2 .

 2 Figure 2. Working area (reference grid).

Figure 3 .

 3 Figure 3. Fibre orientation measurement in a composite prepreg a) ply laid over the working zone b) detail of an element image c) the image after the application of the LOG filter and convolution masks d) edge detection and e) Hugh lines along the fibre direction.

Figure 4 .

 4 Figure 4. Hough lines results comparison at different element sizes a) Number of lines and b) mean measured orientation. Error bars indicate ± 1 standard deviation.

Figure 5 .

 5 Figure 5. Orientation maps from plate C-11 showing four different ply directions with a) 0° b) +45° c) -45° and d) 90°.

Figure 6 .

 6 Figure 6. Schematic description of the pseudo-Gaussian surface representing the local perturbation in fibre orientations.

Figure 8 .

 8 Figure 8. Strain field εxx of the upper skin of a 16-ply 600 x 300 mm composite plate subject to a ΔT = -100°C simulating the cooling phase.

Table 1 . Difference between the intended ply orientation and measured orientation for all plies and plates referenced to ply # 1.

 1

	Ply #	𝜃 𝑘 ̅̅̅̅	C-11	C-12	C-13	δθk	C-21	C-22	C-23

Table 2 . Mechanical properties of the M10.1/CHS composite ply.

 2

	Property	Magnitude
	Young's modulus of the composite in the longitudinal direction E1 (GPa)	131.4
	Young's modulus of the composite in the transverse direction E2 (GPa)	9.3
	Shear modulus of the composite G12 (GPa)	4.3
	Poisson's ratio of the carbon fibre ν12	0.34
	Coefficient of linear thermal expansion of the composite in the longitudinal direction α1 (°C -1 )	-0.5•10 -6

Acknowledgements

The authors would like to acknowledge the Consejo Nacional de Ciencia y Tecnologia (CONACyT) of Mexico for providing the funding for the PhD thesis work of Yves Davila.

It has been shown that the measuring method can be used to accurately determine the mean misalignment of the complete ply. A more detailed analysis is required to determine the perturbation misalignments contained within each ply.

Local misalignments (local scale within a given ply)

The representative values of amplitude, position and shape of the pseudo-Gaussian perturbations that conform the local fibre misalignments, as described in equation 4, are obtained by means of an in-house detection algorithm. The algorithm is set to find values of the 5 variables, Bi, Xi, Yi, αi, and βi, that minimize the mean square error between the measured misalignments and the generated surface from equation 4.

The algorithm operates in the normalised to zero-mean angle. It is possible that the placement of the centre of any i th perturbation can be outside the edges of the ply, thus this point can be located between -60 and 660 mm in the x-direction and between -60 and 360 mm in the y-direction. Afterwards, the size of the perturbations λi and γi are considered between 30 and 210 mm. Finally, during the analysis of the treated images, it was found that in any case the maximum difference between any point and the mean ply orientation was within the rage of ± 1°. Thus, it is assumed that the amplitude Bi which the fibres can deviate between -1° and +1°. It is chosen to use N = 12 pseudo-Gaussian surfaces to form local perturbations in a trade-off between the accuracy, the calculation time and the capability of convergence of the detection algorithm. The computation time for each ply goes from 60 s to 80 s using the same CPU described in section 3.

The shape of the generated surface through the present algorithm using equation 4 is in good agreement with the experimental data by both visual and quantitative analysis. Figure 7 shows the misalignment maps of the reinforcement orientation of a -45° ply for the values as measured by the optical treatment (cf. Figure 7 a) compared to a reconstruction of the same ply (cf. Figure 7b). The measurement error has a mean of 0° with values ranging from -0.2° to 0.2° (cf. Figure 7c). For all treated plies the maximum mean error found is under 0.02° and the maximum standard deviation value is 0.2°. Furthermore, the error scatter in each image is normally distributed, hypothesis confirmed by a Kolmogorov-Smirnov goodness-of-fit test with a confidence of 95 %. This confirms that the differences between the measured image and the reconstructed parameters are mainly due to noise in the measurements. The objective of this methodology is not to introduce the measured misalignments directly into a FE model, but to introduce a computer-generated set of misalignments representative of the observed misalignments. To be representative, the virtual misalignments must have the same order of magnitude and follow the identified variation trends. The image analysis results are used as the basis to obtain the different meta-model coefficients needed to generate the corresponding FE models with continuous variations in the reinforcement angle, using equations 1 through 4.

6.

Finite element modelling of the residual strains of a composite plate with varying reinforcement orientation.

In a previous study [START_REF] Davila | Spatial Evolution of the Thickness Variations over a CFRP Laminated Structure[END_REF], a FE model was used to determine the variations in residual strains at the end of the cooling phase during manufacturing a composite plate taking into account variations in the ply thickness.

The example here presented serves to determine the influence in the strain field due to local variations in the reinforcement orientation. The analysis is based on the assumption that the composite is stress-free at the end of the curing phase (120°C for the M10.1/CHS composite material). Due to the mismatch of the thermal properties of the constituent materials, a temperature variation induces residual strains and stresses [START_REF] Mulle | Assessment of cure-residual strains through the thickness of carbon-epoxy laminates using FBGs Part II: Technological specimen[END_REF].

The FE model uses a shell element type in the LMS Samcef ® (Siemens, formerly Samtech, Belgium) code. The shell composite element allows local control reinforcement orientation for each element and each ply. The complete 600 x 300 mm composite plate is analysed using a 2D Mindlin composite shell (element type T028/T029). The element size is 1 mm per side in order to model reinforcement undulation corresponding to a minimum wavelength of 2 mm. The applied temperature change is ΔT = -100 °C, simulating the cooling step in the polymerisation cycle.

The material properties are presented for a 0.300 mm thick unidirectional ply of M10.1/CHS with a 56 % fibre volume ratio. The typical values for this composite material are shown in Table 2. In this example, only the variation of the reinforcement orientation is taken into account.

Conclusions

A methodology to introduce the variation of the reinforcement orientation in FE models of composite structures was presented. This methodology includes the generation of meta-models that describe the orientation variation in the different scales of the composite, for local misalignments and the average ply misalignment. The values that were introduced into the meta-models were obtained from a real composite layup by means of automatic image post-treatment.

Image analysis results show that variations in the mean ply orientations are found between -0.5 and 0.5° and the standard deviations range between 0.34° and 0.41°. No correlation was found between the overall ply misalignment and the intended orientation of such ply. The identification algorithm that generates the set of parameters introduced into the meta-models yield extremely good results in terms of calculation time and the misalignments measurements accuracy validated both visually and quantitatively. Finally for the FE model results show variations up to 28 % in the strain field of a composite plate subjected to cooling due only to variations in fibre orientation.

An aspect still unexplored is the interaction between two or more sources of variability. In this paper only the in-plane variation of the reinforcement angle are assessed. Further work is required to introduce the effects of out-of-plane waviness, which can be coupled with a corresponding ply thickness variation. It is important to consider that the sources of variability are highly dependent of the process and material semiproducts (i.e. raw materials, composite prepreg, etc.). Hence, great efforts should be conducted to identify and quantify the trends and magnitudes of the different sources of variability that are specific to each application.


Yves Davila
email: yves.davila@iut-tlse3.fr

Laurent Crouzeix

Bernard Douchin

Francis Collombet

Yves-Henri Grunevald

Spatial Evolution of the Thickness Variations over a CFRP Laminated Structure

Keywords: Variability, Composite structure, Thickness, Discrete Fourier transform, Finite element analysis (FEA), Autoclave, Unidirectional prepreg

Introduction

Polymer matrix composite materials present a higher variability on their mechanical properties compared to classic engineering materials, such as metals. Such variabilities can be traced from the mesoscopic through structural scales. However the great number of variables make of this task a great challenge [START_REF] Davila | Étude multi-échelle du couplage matériaux-procédés pour l'identification et la modélisation des variabilités au sein d'une structure composite[END_REF]. One of these sources of variability is the structure thickness. For design purposes, the thickness of a composite structure is given by the sum of the ply thicknesses forming the laminate, or alternatively, the ply thickness is obtained by dividing the total structure thickness by the number of plies contained in the stratification. The thickness of the structure is considered constant and equal for all plies, provided that there is no change of geometry in the cross-section of the composite structure [START_REF] Barkanov | Optimal weight design of laminated composite panels with different stiffeners under buckling loads[END_REF]. In reality, each ply in the stratification does not present a constant thickness mainly due to manufacturing conditions and particularly the resin flow effect before polymerisation in autoclave. Figure 1 exhibits, for different scales, the variations pf ply thickness across a 16-ply unidirectional carbon/epoxy laminate. The effects of the ply thickness variations have been studied in the literature as localized defects due to manufacturing conditions, in the form of wrinkles and warpings [START_REF] Barkanov | Optimal weight design of laminated composite panels with different stiffeners under buckling loads[END_REF][START_REF] Li | Numerical and Experimental Study on the Effect of Lay-Up Type and Structural Elements on Thickness Uniformity of L-Shaped Laminates[END_REF], and on singular zones such as L-shaped stringers [START_REF] Lightfoot | A new mechanism for the formation of ply wrinkles due to shear between plies[END_REF]. Conversely, there are few publications that address the representation of thickness variations spread over a laminated structure. Collombet et al. [START_REF] Collombet | Contribution of Embedded Optical Fiber with Bragg Grating in Composite Structures for Tests-Simulations Dialogue[END_REF] reported a mean ply thickness of 0.26 mm with a standard deviation of 0.002 mm, with extreme values between 0.26 mm and 0.27 mm, in a 28-ply M21/T700 carbon/epoxy laminate. However, the determination of the ply thickness was performed through the measurement of the laminate thickness. A further inspection of the cross section of the specimens revealed that the ply thicknesses deviated up to 21 % from the mean ply thickness, with extreme values of 0.25 mm and 0.31 mm [START_REF] Collombet | Contribution of Embedded Optical Fiber with Bragg Grating in Composite Structures for Tests-Simulations Dialogue[END_REF] It is well-known that the composite structure cannot be dissociated from its manufacturing conditions. The volume fractions of constituent materials (fibre and matrix) as well as porosities are directly linked to the stratification of the composite, the nature of the raw materials and the curing cycle [START_REF] Collombet | Contribution of Embedded Optical Fiber with Bragg Grating in Composite Structures for Tests-Simulations Dialogue[END_REF][START_REF] Olivier | Comparison between longitudinal tensile characteristics of thin and thick thermoset composite laminates: influence of curing conditions[END_REF]. These changes in the volume fraction can be a source of warpages for thin laminates (global on the scale of the plate and "local" on the scale of the plies). Daniel et al. [START_REF] Li | Numerical and Experimental Study of the Bleeder Flow in Autoclave Process[END_REF][START_REF] Hsiao | Effect of fiber waviness on stiffness and strength reduction of unidirectional composites under compressive loading[END_REF] have studied the effects of the fibre waviness on thick unidirectional stratifications (up to 150 plies) generated by local and periodic perturbations in the stratification. The proposed method calls for a modification in the compliance matrix by adding an out-of-plane component to the fibre orientation. Nonetheless the problem is treated as a local effect related to defects present in the material.

A classic computational analysis of a composite structure is done by employing a unique value of ply thickness for all the elements [START_REF] Chun | Effects of material and geometric nonlinearities on the tensile and compressive behavior of composite materials with fiber waviness[END_REF][START_REF] Lee | Determination of robustness for a stiffened composite structure using stochastic analysis[END_REF]. Conversely in the literature, no methodology was found that accounts for the local changes spread over the composite structure applied to multidirectional laminates fabricated with UD prepregs.

The aim of this study is to include the "continuous" ply thickness variation in a finite element (FE) model considering actual variations of thicknesses in the composite laminate. This is essential in order to make more realistic property predictions.

To retrieve the required information needed by the numerical model, the cross-section of a 16-ply composite plate is analysed to obtain the real profiles of the ply thicknesses spread over the plate. A discrete Fourier transform (DFT) of the measured ply thickness profiles provides the model parameters.

These parameters are analysed to be then introduced into a mathematical model to produce a digital representation of the ply variations. This mathematical model can be modified by the introduction of random parameters that modify the shape of digital lay-ups for their use in different simulations. The set of digital lay-ups containing the variations of the plate and ply thicknesses is then used to calculate the material property of the composite ply that will be afterwards introduced into a FE model.

Material

A set of composite plates was fabricated using a unidirectional carbon/epoxy prepreg HexPly ® M10.1/38%/UD300/CHS supplied by Hexcel Composites. The dimensions of the manufactured plates are 600 mm long and 300 mm wide. The 16-ply stratification was chosen to be quasi-isotropic with the lay-up sequence [90/-45/0/(+45)2/0/-45/90]s. The plates were cured in autoclave using the recommended cure cycle adding a consolidation dwell at 90 °C for 15 minutes. The cure dwell was maintained at 120 °C for 60 minutes. The heating and cooling rates were set to 2 °C min -1 . The consolidation pressure at the first dwell and the curing pressure in the second dwell were 2 and 5 bar respectively.

Thickness measurements

Overall plate thickness

The overall thickness of the cured plates is measured using a 3D coordinate-measuring machine (CMM) with a mechanical probe and a resolution of 0.01mm. With a mesh size of 30 x 30 mm, between 190 and 220 measurement points are needed to cover the complete plate. Figure 2 illustrates the thickness mapping reconstructed using the measurements. The 5.45 mm maximum thickness value is located at the centre of the plate while the 4.68 mm minimum thickness value is closer to the edges and corners of the plate. The difference between the maximum and minimum values is 0.77 mm. This difference is approximately 2.5 times the thickness of the cured ply, considering a nominal ply thickness of 0.32 mm.

Olave [START_REF] Olave | Internal Geometry Variability of Two Woven Composites and Related Variability of the Stiffness[END_REF] reported variations in the order of 0.80 to 1.00 mm for a plate of 300 x 300 mm dimensions fabricated with a HexPly ® M10.1/T700 2/2 twill woven textile. Notice that both types of plates, the ones studied by Olave and the ones presented here, were fabricated under similar manufacturing conditions.

Ply thickness profile

To determine the profiles of the ply thickness, a 45 x 45 mm specimen was removed from one of the manufactured plates. It exhibits a 4.66 mm mean plate thickness with 4.84 mm and 4.53 mm, maximum and minimum values respectively. Being a section of the composite plate, the difference between the maximum and minimum values is 0.31 mm (compared to 0.77 mm for the complete composite plate).

The specimen was placed on mounting resin to prepare its cross-section for micrographic inspection. A series of micrographs of the cross-section was taken with a Dino-Lite ® optical microscope with a 161 pixels/mm resolution. For determining the spatial evolution of the ply thicknesses, the micrographs were assembled into a single image containing the full length of the specimen to avoid duplicated points on the ply interfaces, as each micrograph only covers approximately 7 mm (of the 45 mm) of the length of the specimen. An automatic determination of the ply interfaces by image analysis is preferred. However, the low contrast and fuzziness at the ply interfaces makes it difficult to analyse the images with good precision. Hence, the plies were then delimited by manually selecting the points along the ply interface. A grid on top of the image provides a reference in order to have a rough control of the number of points and their spacing. Approximately 220 points were selected for each interface. The selected points were interpolated into a regular mesh consisting of 1024 points along the x-axis. Figure 3 shows the ply interface selection for the interfaces # 2, # 3 and # 4 (delimiting plies # 2 and # 3). A summary of the basic statistical measurements for all the plies and the composite plate thicknesses is shown in Table 1. Due to the lay-up sequence, the pairs of plies # 4/# 5, # 8/# 9 and # 12/# 13 share the same orientation, thus their interfaces were not discernible and they were counted as single plies. Nonetheless, for the present analysis their mean value, as well as the extreme values, are divided by 2 and their standard deviations are obtained thus by √𝜎 ² 2 ⁄ . 4. While for the ply thickness profile, using the case of ply #2 as example (cf. Figure 5), the CV comprises a difference at both ends of the specimen of 0.02 mm and, more importantly, large variations up to ± 29 % of the mean ply thickness (cf. Table 1). It is noteworthy that the difference between the slope of the linear regression (m) of the ply profiles and the mean thickness (zero slope) is not negligible (cf. Figure 5). All the plies profiles and the overall plate profile also exhibit negative slopes (m < 0). As expected, the sum of the slopes of the ply profiles is equal to the slope of the plate profile itself. Indeed, the variation of the overall plate thickness is strongly coupled with the ply morphology.

From the thickness profiles, it can be seen that there is an apparent repetition of the thickness variation along the length of the specimen. For this reason, the ply thickness profile can be analysed by means of signal treatment techniques. Indeed, the profiles can be considered as a sum of N sine waves of different frequencies and amplitudes. This consideration leads towards a digital representation of the ply profiles that can be generated by a set of frequencies representative of the cross-section measurements.

Modelling the ply thickness variations

Analysis of the thickness profiles

The collection of thickness profiles are analysed in order to determine a set of parameters representative of the material and manufacturing conditions. For this analysis the thickness profiles are converted from the spatial domain to the frequency domain using a discrete Fourier transform (DFT).

The DFT is applied to the thickness profiles normalised to zero-mean zk as shown in equation 1.

      m k k p m k m k t x t x t x z     (1)
where zk m is the zero-mean signal from the measured k th ply thickness profile tk m and t̅ k m is the mean of the measured thickness for each k th ply. Since DFT peaks at very low frequencies of the spectrum can be accounted for in the overall plate thickness variation, a function δtp k is considered. For this analysis, δtp k is the equation obtained from the linear regression of each thickness profiles.

The higher frequencies in the spectrum are mostly attributed to the manual selection of the ply interfaces. A low-pass Butterworth filter is applied to the zero-mean signal zk m to remove the frequencies greater than 1 mm -1 present in the thickness profile to remove the uncertainty of the point selection at the ply interface. It should be noted that these frequencies cannot be modelled in a finite element model.

A total of 14 DFT spectra were analysed, 13 for the ply profiles (considering three double plies) and one for the plate profile. An example of the resulting DFT spectrum for the ply # 2 is shown in Figure 6.

Measurement window

The peaks of each spectrum are selected within a measurement window imposed by required finite element (FE) model mesh fineness which means it is impossible to model physics finer than the mesh itself. In this study, the measurement window is chosen to go from 0 to 0.50 mm -1 in frequency, and from 0.0035 mm upwards in amplitude (cf. Figure 6). The analysis of the collected data was difficult due to the coupling of at least three variables, frequency, amplitude and phase shift, without mentioning that each profile spectrum contains a different number of peaks.

Therefore, to introduce the ply thickness variations into the FE model, a single set of N = 10 sine waves for all K plies (K = 16 for the present study) is created using an in-house pseudo-random generator of the zero-mean signal zk:

        N n n n n n k x f A p x z 1 2 sin 2   (2)
To have a precise control of the overall plate profile, the coupling of the ply thickness profiles is done by assigning a specific phase shift ϕn to each n th sine wave and each k th ply. As noted, the complete dataset of the selected DFT peaks is divided into ten frequency bands of 0.05 mm -1 each from 0 to 0.50 mm -1 . A single frequency fn and amplitude An is obtained for each n th frequency band. After a goodness-of-fit test for each n th band it was observed that the frequency values were better represented by a uniform distribution (with lower bound a and upper bound b limits) while the amplitude is better described by a normal distribution (with mean x ̅ and standard deviation σ). It is noteworthy that not all the frequency bands contain DFT peaks for all plies. For this reason, the parameter pn is drawn from a binomial distribution with values being either 0 or 1. The probability to generate a value equal to 1 is the probability of appearance within the n th frequency band. The model parameters, which are obtained from the analysis of the different DFT spectra, are shown in Table 2.

Mathematical representation of ply thickness variations

The overall shape of the composite plate is controlled through the contribution of the thickness variation of the plate δtp k . Nevertheless, the exact distribution of the thickness variations across the ply is difficult to know. However, the overall change of the composite plate δtp is known and it is evenly spread over all the composite plies K, thus δtp k can be estimated as:

K t t p p k    (3)
A 3D sine wave was chosen to be representative of the overall plate shape, δtp(x,y) is thus:

                                  c p c p p p Y y X x A y x t y x      2 cos 2 cos 4 , ( 4
)
where Ap is the difference between the maximum and minimum thickness value, λp x and λp y are the plate length and width respectively, and Xc and Yc are the coordinates of the maximum thickness zone.

Figure 7 shows the resulting thickness mapping on a 600 x 300 mm composite plate which is in good agreement, qualitatively and quantitatively, with the reality (cf. Figure 2). The digital lay-up can be formed combining equations 1, 2 and 3 using the pseudo-random generator explained in section 4.1 and employing the parameters established in Table 2. The proposed model of each ply thickness adapted for a 3D representation is thus:

                  N n n k n n k n n n p k k y f x f A p K y x t t y x t 1 , , 2 sin 2 sin 4 , ,      (5)
where tk(x,y) is the generated thickness in both x and y-directions of the k th ply, t k is the mean thickness of the k th ply, An is the amplitude of the corresponding frequency fn for the n th peak, and np takes the value 1 or 0 depending on the probability of appearance of such n th peak. As previously explained, δtp(x,y) is the desired plate profile divided by the number of plies, K.

To ensure a control on the coupling of the different ply thicknesses profiles, a pseudo random phase shift k,n is used for each k th ply and n th frequency peak. The overall plate profile (i.e. the sum of the K ply profiles tk (x)) can maintained constant if the k th sum of the ply shifts must be equal for all n th sinusoids (from 1 to N):

         K k N k K k k K k k 1 , 1 2 , 1 1 , ...    (6)
The phase shift ranges from -π to π for each n th sinusoid:

+ (Xc, Yc)       n k n k K , , 2 (7)
where n k ,  is a vector containing the ply numbers 1 to K randomly shifted for each n th frequency peak.

Ply thickness generation discussion

The main goal of the digital thickness generator is not to reproduce the exact form of the composite ply spread over the plate, but to generate a geometry having the same characteristics, in terms of variability, as the actual ply thickness profiles. Even though the results associated with the proposed methodology (DFT spectra and range of values) are indeed dependent on the geometry, material and manufacturing conditions, the methodology itself does not depend on any of these conditions.

The variations found in the composite exhibit similar patterns but also exhibit differences compared with any point of the plates and with different plates manufactured with the same conditions.

For this reason, the inclusion of random parameters into the mathematical representation allows the generation of new sets of virtual thicknesses that share the same characteristics in terms of variations than the studied plate.

As observed in the cross-section of the composite material, the real stratification exhibits variations in the ply thickness (cf. Figure 8a). In the classic modelling of a composite stratification, the thickness of the different k th is considered as constant and it does not take in account the overall change in the stratification thickness (cf. Figure 8b).

Any effort made in order to completely model the behaviour of the thickness variation will take considerable computing and analytical resources since there are four highly interdependent variables that contribute to the thickness profile. The model must account for the number of the different frequencies and the spectral spacing between DFT peaks, the associated amplitude for each frequency and phase shifts. All these variables also depend on the ply position in the stratification. Different solutions to the problem, such as using the frequency values "as is" from one single ply, or randomly selecting the frequencies and their associated amplitudes from all the identified peaks were discarded, since the selected points come from a very specific zone in the plate that may or may not describe completely the behaviour of the thickness variation along the full composite plate. Indeed, the data collection cannot be performed in all types of composite structures and for different composite plates due to its destructive nature, as well as the time and resources that this task demands. Therefore, simplifications must be made in order to obtain accurate data to be introduced into the model.

The variability of the ply thickness, as well as the coupling between the different plies to obtain a constant plate thickness, can be controlled by using pseudo-random algorithms to vary the parameters involved in the thickness profile generation (cf. Figure 8c).

Determination of the material variations according with the geometrical variations

Accounting for only thickness variations across the composite plate is not enough to determine their effect on local properties. The variation of volume fraction could be crucial since much of the thickness variation is due to differences in the amount of resin. This means that variations in properties could be less than the ones expected based on the thicknesses alone.

To introduce the local properties of the material linked to both ply thickness and the fibre volume fraction, the mass per unit area of the reinforcement ρA f is assumed constant over all the manufactured plates (i.e. for a given plate and from one plate to another). A micromechanical model of the fibre distribution along the ply could be implemented [START_REF] Huang | Effects of Fiber Arrangement on Mechanical Behavior of Unidirectional Composites[END_REF]. However, at the moment, it is difficult to determine with good precision this property due to the rearrangement of fibre bundles during the curing phase. Also the computation time would increase significantly. Assuming that the quantity of fibres cannot drastically vary along the reinforcement direction (in a continuous fibrous reinforcement), a rough and quick estimation of the fibre volume fraction Vf k (x,y) at any point in the stratification can be thus made without the need of a micromechanical model of the fibre distribution inside the ply:

    y x t t y x V k f k f , ,  (8)
where tk(x,y) is the local thickness in the xy-coordinates and in the k th ply, and an tf is the equivalent thickness of the reinforcement defined by:

f A f f t    (9)
where ρA f is the mass per unit area of the reinforcement and ρf the density of the reinforcement. Since both parameters are considered constant, the equivalent thickness of the reinforcement tf is constant. The ρA f value is calibrated to 305 g•m -2 from observations of the M10.1/CHS manufactured plates.

Finally, the volume fraction of the matrix Vm k (x,y) is given as follows: [START_REF] Chun | Effects of material and geometric nonlinearities on the tensile and compressive behavior of composite materials with fiber waviness[END_REF] change is applied to calculate the residual strains at the end of the cooling phase during manufacturing.

      y x V y x V y x V p f m k k , , 1 ,    
The analysis is typically based on the assumption that the composite is stress free at the curing temperature. Due to the mismatch of the thermal properties of the constituent materials, a temperature variation induces residual strains and stresses [START_REF] Mulle | Assessment of cure-residual strains through the thickness of carbon-epoxy laminates using FBGs Part II: Technological specimen[END_REF]. The complete 600 x 300 mm composite plate is analysed using a 2D Mindlin composite shell (element type T028/T029). The element size is set to be 1 mm per side in order to model a ply thickness undulation corresponding to a minimum wavelength of 2 mm (cf. Figure 6). Three cases are analysed, the first one considers a constant plate and ply thicknesses. The second case uses variable laminate profiles generated by the algorithm explained in section 4.2 with constant material properties. For the third case, the ply and plate geometry is the same as case 2 but the material properties are calculated for each element and each ply as per section 5. The mean ply thickness for each plate is 0.30 mm, thus the plate thickness is 4.84 mm (cf. Figure 1). For all 3 cases, no porosity is assumed Vp(x,y) = 0.

The deterministic solution without taking account of the thickness variations (case 1) shows as expected a uniform strain field εxx having a value of -83•10 -6 mm/mm marked as a reference value in Figure 9. Figure 9 shows the strain field εxx of the upper skin for case 2 (Figure 9 top) and case 3 (Figure 9 bottom). The model without accounting for the material variations (cf. Figure 9 top) over estimates the compressive strains within a range of -80•10 -6 to -90•10 -6 mm/mm. It is also evident that the strain field forms uniform bands across the plate width. Conversely, the plates with variable thicknesses and variable materials (cf. Figure 9 bottom) exhibits a strain field with values ranging from -65•10 -6 up to -85•10 -6 mm/mm. As shown, an attenuation of the residual strains is noticed, however the range of values Δεxx passed from 10•10 -6 to 20•10 -6 mm/mm, signifying an increase of an order of magnitude of the internal variations. Also, the strain field distribution is completely different, with the most significant strains at the centre of the plate and a relaxation of such residual strains towards the corners of the plate.

Conclusions

The first input parameter towards the creation of a finite element model, affected by material variability, is the geometry variation of the ply thickness over the plate. Indeed, variations of this property impact the volume fractions of the constituent materials thus modifying the local material properties of the composite structure.

The thickness measurement of the composite plate studied reveals variations approximately

Strain concentration spots

Plate with thickness variations and constant material properties Plate with thickness and material property variations mm/mm mm/mm 2.5 times the thickness of a cured ply within the same composite plate. The observations made of the cross-section of the composite plate reveal continuous and periodic variations of the profiles and a non-linear coupling between the different ply thickness profiles and the resulting plate profile. This coupling is reflected by a low coefficient of variation (CV) of the plate thickness being less than ± 2 % of the mean plate thickness, while, the ply thicknesses exhibit a CV of 7 % of the mean ply thickness with maximal deviations up to ± 29 %.

A mathematical model of the ply variations is proposed to generate a digital laminate having similar dispersion on the ply thickness and having a strict control on the plate thickness shape. The frequency and amplitude values necessary to generate the digital ply are identified by the DFT analysis of the measured ply and plate thickness profiles. This modelling method produces ply thickness profiles having variations in accordance with the real ply and plate profiles.

The capabilities of the model are shown by computing the residual strains of the composite plate during the cooling phase of the autoclave processing. It is shown that accounting for geometrical and material variations generates a different strain field than for a classic composite plate without variations.

Figure 1 .

 1 Figure 1. Micrograph showing the ply thickness variations tk(x,y) spread over a layered composite plate (16ply HexPly® M10.1/38%/UD300/CHS) of approximately 4.8 mm thick t̅ p.

Figure 2 .

 2 Figure 2. Example of measured tp(x,y) thickness distribution of the HexPly ® M10.1/38%/UD300/CHS plate.

Figure 3 .

 3 Figure 3. Profiles for the interfaces # 2 # 3 and # 4 delimiting plies # 2 and # 3 with the cross markings showing the selected points along the ply interface, and the continuous lines showing the interpolated interfaces.

Figure 4 .Figure 5 .

 45 Figure 4. Profile of the measured plate thickness and the least squares linear fit (dashed line) and the mean thickness (dashed horizontal line).

Figure 6 .

 6 Figure 6. DFT spectrum of the thickness profile of ply # 2 with selected frequency peaks (square markings).

Figure 7 .

 7 Figure 7. Example of a numerical overall plate thickness generated using eq. 4.

Figure 8 .

 8 Figure 8. Representation of the cross section of the 16 ply stratification (plies 4/5, 8/9 and 12/13 shown as independent plies) with, a) reconstruction from the actual ply thickness profiles (after filtering), b) constant ply thickness accounting the contribution of the plate thickness variation in each ply and c) digital stratification generated by the proposed model.

Figure 9 .

 9 Figure 9. Strain field εxx of the upper skin of a 16-ply 600 x 300 mm composite plate subject to a ΔT = -100°C simulating the cooling process of a composite (top) no material variation assumed (bottom) material vary according with thickness.

Table 1 . Statistical summary for the ply and plate thickness measurements.

 1 Figure 4 and 5 show the measured thickness profile, mean value and linear regression for the composite plate and the # 2 ply respectively. The coefficient of variation (CV) of the composite plate is 2 %. In comparison, the CVs of the composite plies fall within the range of 4 % to 9 %. It should be noted that the CV of the plate comprises a large difference in thickness at both ends of the specimen length (0.31 mm), shown by the linear regression in Figure

	Ply #	Mean thickness	Standard deviation	CV	Maximum	Minimum
		(mm)	(mm)	(%)	(mm)	(mm)
	1	0.28	0.018	6	0.34	0.23
	2	0.31	0.020	6	0.36	0.25
	3	0.29	0.021	7	0.35	0.21
	4 and 5	0.29	0.021	7	0.35	0.25
	6	0.31	0.026	8	0.39	0.24
	7	0.29	0.026	9	0.37	0.21
	8 and 9	0.28	0.021	7	0.31	0.25
	10	0.30	0.022	7	0.35	0.25
	11	0.30	0.022	7	0.36	0.23
	12 and 13	0.29	0.019	7	0.33	0.26
	14	0.30	0.024	8	0.37	0.22

Table 2 . Parameters of the distribution for the Frequency (uniform) and Amplitude (normal) associated with the probability of appearance of the n th peak within a frequency band.

 2

	n th peak	Frequency band limits	Probability of appearance	Frequency (mm -1 )	Amplitude (mm)
	#.	(mm -1 )	(%)	a	b	x̄	σ
	1	0.00 to 0.05	100	0.020	0.049	0.0081	0.0033
	2	0.05 to 0.10	100	0.054	0.099	0.0088	0.0036
	3	0.10 to 0.15	80	0.106	0.146	0.0095	0.0026

Acknowledgements

The authors would like to acknowledge CONACyT of Mexico for providing Yves Davila the financing for his PhD program. The authors would also give special thanks to Peter Davies for his invaluable comments and remarks.

where Vp(x,y) is the local volume fraction of porosity.

Subsequently, using the computed volume fractions and the constituent materials properties expressed in Table 3, the elastic properties are calculated by means of the classic law of mixtures. The proposed modelling is only intended for obtaining a quick estimation of the material properties linked to the variation of the ply thickness. For this reason, it does not intend to model the lamina morphology (i.e. fibre bundles distribution, out-of-plane misalignments, porosity distribution, etc.). Nevertheless, the overall method does not exclude a more complete modelling of the lamina micromechanics depending on the available computational resources.

6.

Finite element modelling of the residual strains of a composite plate with variable thicknesses.

The finite element (FE) model is required to be robust and light-weight. It uses a shell element type in the LMS Samcef ® (Siemens, formerly Samtech, Belgium) code. The shell composite element allows local control of the ply thickness for each element and each ply.

To illustrate of the capability to introduce local thickness variations into a FE model, a composite plate subjected to a temperature change of ΔT = -100 °C is presented as an example. This temperature