Véronique Michaud

Malin Åkermo

Joël Breard

Luisa Silva

Keywords: scale validation test case

for the inspiring discussions around my work. I wish to express my gratitude to Ana Carolina Quintanilha for the help in the experimental work, and to all the colleagues in GeM that shared some part of this journey with me.

I want to thank my family. Facing a global pandemic

Industrial context and motivation

Being able to predict the outcome of a manufacturing process still remains one major challenge in the industry of composite materials. Resin Transfer Moulding (RTM) is a family of Liquid Composite Moulding (LCM) processes used to produce composite structural parts in which a fibrous preform is compressed inside a rigid mould and impregnated with a polymeric resin (Fig. 1.1), which is later solidified to obtain a finished part [START_REF] Advani | Process modeling in composites manufacturing[END_REF]. When producing RTM parts, engineers are interested in predicting the filling pattern and anticipate defects that may arise due to an insufficient infiltration. This is typically achieved using numerical simulations to model the macroscopic flow inside the fibrous material using Darcy's law [START_REF] Michaud | A review of non-saturated resin flow in liquid composite moulding processes[END_REF][START_REF] Liu | Modelling and simulation of resin transfer moulding (RTM)-gate control, venting and dry spot prediction[END_REF]. However, simulations rely on assigning the proper permeability data to the material. This upscaled tensorial quantity expresses the compliance of a textile to be infiltrated by a fluid and depends on the material microstructure. In the industrial context, permeability is estimated mostly by experimental characterization on a flat sample material [START_REF] Parnas | Permeability characterization. Part 1: A proposed standard reference fabric for permeability[END_REF][START_REF] Kuentzer | Permeability characterization of dual scale fibrous porous media[END_REF]. However, for the same layup and fiber volume fraction, different values of permeability can be registered from part to part [START_REF] Arbter | Experimental determination of the permeability of textiles: a benchmark exercise[END_REF]. Variability arises mostly from the mesoscale architecture of the part, which originates from the preforming stage [START_REF] Bodaghi | On the variability of permeability induced by reinforcement distortions and dual scale flow in liquid composite moulding: A review[END_REF], but also from the intrinsic variability of textiles [START_REF] Gommer | Quantification of mesoscale variability and geometrical reconstruction of a textile[END_REF]. In recent years, significant efforts have been devoted to the numerical characterization of permeability by simulating the flow inside a digital Representative Volume Element (RVE) of the material [START_REF] Verleye | Permeability of textile reinforcements: Simulation, influence of shear and validation[END_REF][START_REF] Verleye | Permeability prediction for the meso-macro coupling in the simulation of the impregnation stage of Resin Transfer Moulding[END_REF][START_REF] Wong | Comparisons of novel and efficient approaches for permeability prediction based on the fabric architecture[END_REF]. The challenge for this approach is in the acquisition of the mesoscale geometric data to create a realistic digital material [START_REF] Lomov | Compression of woven reinforcements: a mathematical model[END_REF]. The geometry of an RVE can be measured in great detail using µ-CT, but the high cost and the small size of sample that can be studied mean that this technique is of little interest for the industry context [START_REF] Desplentere | Micro-CT characterization of variability in 3D textile architecture[END_REF]. Furthermore, due to the design of composite parts, the permeability of a preform is usually not homogeneous throughout the part (ply drop-offs, change in fiber orientations and volume fraction), hence cannot be represented by a single value extracted at one location. [START_REF] Facciotto | Modeling of anisotropic dual scale flow in RTM using the finite elements method[END_REF] Objective The objective of this thesis is to develop a novel workflow for RTM process simulation that relies on real material data, bypassing the need for traditional permeability characterization with liquid flows. The research work can be divided conceptually in two parts: geometrical data acquisition and process simulation. In the first part (chapters 2, 3) we introduce two new non-destructive methods to acquire mesoscopic geometrical data from quasi-UD textiles. The measurement is cost-effective and not limited in size, because based on using an affordable pressure-sensitive film [START_REF] Liggins | The practical application of Fuji Prescale pressure-sensitive film[END_REF] that can be laid inside the mould. In the second part (chapters 4, 5, 6) the acquired data is processed to develop an efficient model of flow simulation, based on the principle of geometry skeletonization and model reduction. The outcome is a dual-scale flow model that can be used to simulate directly the filling of a multi-layer preform, at the meso-scale, based on real material data and without the need of any permeability characterization.

Although both parts can be treated independently, they complement each other in a seamless workflow. The scientific challenges related to each task and their place in the thesis are presented in the next section.

Scientific challenges

The challenge in the characterisation of fine geometrical features of a textile preform comes from the intrinsic variability of highly deformable dry fabrics. The meso-scale structure of a fibrous preform is the result of the complex interaction between multiple layers of textile which are stacked and compacted inside a mould [START_REF] Bodaghi | On the variability of permeability induced by reinforcement distortions and dual scale flow in liquid composite moulding: A review[END_REF]. Traditional techniques to characterize the meso-scale geometry of a fibrous preform are available. Photographic inspection techniques such as the one developed by Gommer et al. [START_REF] Gommer | Quantification of mesoscale variability and geometrical reconstruction of a textile[END_REF] are effective, but limited to the visible surface layer. µ-CT acquisition is very accurate (Fig. 1.2) but the size of investigation is limited to a few centimeters. Furthermore the high cost and time-consuming operation make this option more suited for research than for industry.

The meso-scale configuration of the material can also be affected by the intrinsic variability that textiles exhibit before they are even stacked in a mould. This depends on the degree of reproducibility of textile manufacturing processes, which in turn comes from the 1.2. Scientific challenges Figure 1.2: Examples of quasi-UD textiles. a) section view [START_REF] Hemmer | Unloading during the infusion process: direct measurement of the dual-scale fibrous microstructure evolution with X-ray computed tomography[END_REF]; b) µ-CT scan [START_REF] Dharmalingam | Evolution of single carbon and glass fibrous tow crosssections in dry and lubricated states during compaction perpendicular to the fibers[END_REF] of a glass quasi-UD textile; c) schematic view of quasi-UD tows [START_REF] Binetruy | The interactions between flows occurring inside and outside fabric tows during RTM[END_REF] highlighting inter and intra yarn pores non-cohesive nature of textiles. Variability can affect not only the shape of a yarn crosssection, but also the spacing between their paths [START_REF] Gommer | Quantification of mesoscale variability and geometrical reconstruction of a textile[END_REF].

The aim of chapter 2 is to use the an affordable pressure-sensitive film to study the yarn variability in a single layer of quasi-UD textile. The pressure information obtained from the compression of the textile is used to interpret the initial thickness and width of each yarn.

The use of pressure information to interpret geometrical features in a textile is novel, but the single-layer study should be seen as a first attempt of application of the pressure sensor. In RTM composite processing, the interest in studying one single layer textile is low since all preforms are always formed by multi-layer stacks of textiles. Extending the methodology to the multi-layer case, the detailed information valid for single yarns is not obtainable anymore. The mechanical approach is abandoned, however the periodic structure of the material opens up the possibility to other types of analysis.

In chapter 3 the analysis is extended to multi-layer stacks. A pressure-sensitive film is once again used to interpret the meso-scale structure of a preform under compression. By combining Moiré theory of interference [19] and image spectral analysis, the precise orientation and spatial distribution of yarns are measured. The outcome is a three-dimensional digital model of the multi-layer preform.

Knowing the geometrical information of a preform with mesoscopic detail allows one to build realistic digital models of the material [START_REF] Sherburn | Geometric and mechanical modelling of textiles[END_REF][START_REF] Verpoest | Virtual textile composites software WiseTex: Integration with micro-mechanical, permeability and structural analysis[END_REF]. These are in turn used as input data to compute flow simulations with the purpose of material characterization [START_REF] Verleye | Permeability of textile reinforcements: Simulation, influence of shear and validation[END_REF]. The challenge associated with numerical simulation of flow in porous media is in the computational complexity. When the simulation is carried out at the meso-scale, the problem to solve can be computationally demanding. As a consequence, flow simulations are usually limited to small sample domains (RVEs). Properties extracted from one localised position however are hardly representative of the entire part, since the meso-structure of the material could be different simply a few lengths away. Instead of using local information to compute an upscaled permeability value, in this thesis we propose to directly simulate the 1.2. Scientific challenges mesoscopic flow of resin in the preform. We also look for a computationally efficient model to simulate the flow over a large domain, in line with the scale of measurement enabled by the pressure-sensitive film, which is orders of magnitude larger than a traditional RVE. Therefore, some simplification to the flow model must be introduced. The network of interyarn pores, where the resin flows, is composed by thin channels in which the gap is much smaller than the characteristic length of the flow. This allows one to reduce the computational cost of the problem by adopting the lubrication approximation, essentially reducing the complexity of the problem by one dimension. This assumption is well established in soils [START_REF] Zimmerman | Lubrication theory analysis of the permeability of rough-walled fractures[END_REF][START_REF] Ranjith | Applicability of the 'cubic law'for non-Darcian fracture flow[END_REF][START_REF] Dippenaar | On the cubic law and variably saturated flow through discrete open rough-walled discontinuities[END_REF] and has been applied successfully to textile materials [START_REF] Wong | Comparisons of novel and efficient approaches for permeability prediction based on the fabric architecture[END_REF][START_REF] Barthélémy | Effective permeability of media with a dense network of long and micro fractures[END_REF][START_REF] Levasseur | A micro-macro approach of permeability evolution in rocks excavation damaged zones[END_REF]. However, in order to use lubrication theory for RTM simulations, a compatible geometrical space is needed.

In chapters 4 and 5 we introduce the Single-Scale Skeleton model: from the threedimensional pore space of textiles, a medial skeleton mesh is extracted. Reynolds equation of lubrication is solved numerically on the skeleton mesh, for both saturated and unsaturated flow. This model can be used to numerically compute the permeability of a preform using only the geometric data coming from the pressure print. More importantly, the model allows one to carry out direct filling simulations, at the meso-scale, completely bypassing the need for an upscaled permeability. The computational efficiency of the model allows one to simulate an RTM infusion process over a domain of large size. Local defects such as racetracks can also be included directly. A series of test cases are provided and the potential of the skeleton model is discussed.

The first version of the skeleton model only accounts for single-scale flow in the pore space of inter-yarn channels. However, the flow of resin impregnating a textile preform is affected by two characteristic lengths: the size of the pores in-between the yarns (the inter-yarn channels) and the size of the pores within the yarns (the gaps between individual fibers) [START_REF] Binetruy | Tow impregnation model and void formation mechanisms during RTM[END_REF]. Consequently, the flow of resin during an RTM process is dual-scale and models exist to account for it [START_REF] Bodaghi | On the variability of permeability induced by reinforcement distortions and dual scale flow in liquid composite moulding: A review[END_REF][START_REF] Li | Derivation of modified Reynolds equation-a porous media model[END_REF][START_REF] Syerko | A numerical approach to design dual-scale porosity composite reinforcements with enhanced permeability[END_REF]. Accounting for the presence of the permeable yarns is important [START_REF] Tahir | Effect of dual scale porosity on the overall permeability of fibrous structures[END_REF] because local defects such as dry spot formation are likely to be found there [START_REF] Simacek | A numerical model to predict fiber tow saturation during liquid composite molding[END_REF].

Building on the single-scale model, in chapter 6 we propose a Dual-Scale Skeleton (DSS) model. The permeable yarns are added to the domain available for the resin flow in the form of a second skeleton mesh, interconnected with the original channels mesh. A Reynolds-Darcy model is proposed to simulate the dual-scale flow. The model is applied to a series of test cases. The DSS model can simulate the transient filling of resin on a hybrid skeleton mesh representing the meso-scale geometry of both the channels and the yarns.

The outcome of chapter 6 fulfils the objective of proposing a seamless workflow: from the affordable acquisition of geometric data to the numerical simulation of a dual-scale RTM filling process, bypassing any measurement of permeability.

Chapters structure

In chapter 7 we draw the conclusions, highlighting the potential and the downsides of the novel methodology we propose, leading to the future prospects of development.

Chapters structure

Chapters 2 and 3 have been previously published in the form of two journal articles [START_REF] Bancora | Inverse methodology as applied to reconstruct local textile features from measured pressure field[END_REF][START_REF] Bancora | Characterization of mesoscale geometrical features of a preform using spectral Moiré analysis on pressure print[END_REF]. The main body of the chapters has been directly extracted from the published versions and Introduction/Conclusions sections have been modified to blend with the flow of the thesis.

Chapters 4, 5 and 6, on the other hand, contain unpublished material and have been written expressly for this document.

The different source of the material explains some potential disparity in style and length between the first and the second parts of the thesis.

Introduction

In this chapter we begin the study of geometrical features in textiles using a novel technique. The variability in yarn dimensions of a single layer of quasi-UD textile is measured by analysing its pressure field resulting from its through-thickness compression. The data acquisition is carried out using a pressure-sensitive film. This technique is a cost-effective, non-destructive tool to conduct the analysis.

The detection of the original width and size for each yarn on the print is gathered by adopting an hyperelastic constitutive law to determine the inverse transformation. The information extracted from the pressure field analysis can be used to generate a digital representation of the original textile.

Although the geometrical study of one single layer of textile has little direct application, this novel methodology is a first attempt to use a compaction pressure field to interpret geometrical features of a textile. The work presented in this chapter represents a first milestone towards the extension to characterizing multi-layer preforms, introduced later in the thesis.

Materials

Outline of the chapter

The methodology presented in this chapter relies on the use of a pressure-sensitive film. Therefore, before detailing the technique to acquire geometric data, the Prescale pressure sensor is introduced and its calibration is carried out. This is followed by the constitutive model used to link the measured stress to the yarn deformation. Next, the characterization of the needed material parameters is performed.

Following that, the methodology for detection of the yarn mesoscale dimensions is detailed and validated on single yarns.

Finally, a test case involving a portion of one textile layer is carried out. The outcome of this step is the generation of a three-dimensional digital model of the original material.

Materials

Prescale pressure film

The measurement of the boundary pressure is achieved using a flexible pressure sensor. Technical developments of recent years spawned many examples of similar sensors based on current-resistive cells [START_REF] Chen | Progress in achieving high-performance piezoresistive and capacitive flexible pressure sensors: A review[END_REF]. These sensors are typically limited in both size and resolution due to added manufacturing costs and applications range from biological measurements [START_REF] Mannsfeld | Highly sensitive flexible pressure sensors with microstructured rubber dielectric layers[END_REF][START_REF] Lee | A transparent bending-insensitive pressure sensor[END_REF] to process monitoring [START_REF] Walbran | Measurements of normal stress distributions experienced by rigid liquid composite moulding tools[END_REF]. In this work the purely mechanical Prescale sensor by Fujifilm is used. This is a pressure-sensitive film capable of registering a distributed pressure applied normally on its surface, yielding a quantification of the field in the form of a color intensity map. Since the measurement is carried out by the rupture of color-releasing microcapsules, this system is cheaper than the aforementioned electrical sensors. Some examples of pressure prints from various textures obtained using a Prescale film are shown in Fig. 2.1. The Prescale sensor has been extensively used in biomechanics to study contact surfaces, in articular joints [START_REF] Liggins | The practical application of Fuji Prescale pressure-sensitive film[END_REF][START_REF] Macdonald | A biomechanical analysis of joint contact forces in the posterior cruciate deficient knee[END_REF] or other mechanical applications [START_REF] Aritomi | A simple method of measuring the footsole pressure of normal subjects using prescale pressure-detecting sheets[END_REF][START_REF] Kitafusa | Application of "prescale" as an aid to clinical diagnosis in orthodontics[END_REF]. Other fields of application include studies on wheel/road [START_REF] Backx | Tyre/road contact measurement using pressure sensitive films[END_REF] or wheel/rail [START_REF] Dörner | On the accuracy of the pressure measurement film in Hertzian contact situations similar to wheel-rail contact applications[END_REF] contact. This particular sensor works well with the chosen material, as it can be laid under the textile when a mould is closed during the compaction stage.

In addition to the mentioned cost effectiveness, another interesting advantage of this film is that the measurement has no limitations in size, unlike electrical sensors or direct geometry measurement techniques such as µ-CT which yield extremely accurate results but limited to a small material sample [START_REF] Desplentere | Micro-CT characterization of variability in 3D textile architecture[END_REF].

To the best of the authors' knowledge this is the first time a pressure-sensitive film is used as an input in an inverse identification method. Several gradations of Prescale films are available that address various ranges of pressure. In this thesis, the LLLW (0.2 -0.6 MPa) and the LLW (0.5 -2.5 MPa) types are used, depending on the range of interest.

Pressure film calibration

The intensity of the color displayed on the film after a pressure is applied depends on the density of ink dots on the surface. Each ink dot is released when a microcapsule breaks, the quantity of broken capsules depending directly on the pressure applied. The manufacturer provides a color scale that can be directly used to link the color intensity on the paper to a value of pressure associated. The relation between pressure and intensity is claimed to be linear. However, since the calibration procedure is straightforward and has been detailed before in literature [START_REF] Liggins | Calibration and manipulation of data from Fuji pressure-sensitive film[END_REF], the calibration of the Prescale sensor was carried out autonomously for the Prescale types LLW and LLLW.

The methodology consists in applying a range of known loads onto a patch of pressure film, over a flat contact area which is exactly known, so that the applied pressure can be calculated. At each application, the marked area is scanned and the average color intensity is extracted, as to create a map between pressure and intensity.

All compression tests for the calibration of the film were performed on an Instron 5584 machine. The detailed procedure is described below:

-In order to control exactly the area of application of the load, an aluminum compression punch is machined in-place. This punch has the shape of a cylinder of radius R = 20mm. Great care is placed in making its top and bottom faces perfectly flat, so that the load applied on the pressure film is uniform and the pressure color print does not present any gradient of intensity. The contact area between the punch and the pressure film is A = R 2 π and the pressure for a given load can be obtained as:

p = F/A.
-The marked print is digitized using a digital scanner Epson Workforce DS-5500. In order for all measurements to be consistent, the image exposure and the acquisition resolution (600 PPI) are set manually and never changed.

-Following already documented procedures [START_REF] Liggins | Image-averaging in the analysis of data from pressuresensitive film[END_REF], the average intensity of the marked area is measured using the ImageJ software [START_REF] Schneider | NIH Image to ImageJ: 25 years of image analysis[END_REF].

The procedure is repeated for imposed compressive load values between 60 and 3000 N, and carried out for both the LLW and the LLLW types of Prescale films. Values of pressure and intensity are plotted in Fig. 2.2 for the LLW and the LLLW cases respectively. The map shows indeed a linear trend, therefore a linear fit can be adopted. When the intensity (I) is expressed in a normalized grayscale value [0-1], coefficients for Eq. 2.1 were characterized for a linear map of the form:

p(I) = C 1 I + C 2 (2.1)
The coefficients were identified as follows:

-LLLW film: C 1 = 2.65MPa, C 2 = 0.05MPa -LLW film: C 1 = 5.40MPa, C 2 = 0.34MPa

Quasi-UD textile

The entirety of this research was developed on quasi-UD textiles. These were chosen over woven fabrics on account of the lower morphological complexity they present. Another reason behind the choice is that UD materials, more than woven, are likely to present long preferential channels for the resin to flow due to their intrinsic anisotropy. This geometrical characteristic of the inter-yarn channels will be of interest in chapters 4-6. Two different quasi-UD textiles are used in the course of this thesis: in this chapter, the material used is the UDT400 glass fabric. The UDT400 glass fiber textile is an unbalanced quasi-UD textile (A w = 423g/m 2 ) from Chomarat. The fabric material is shown in Fig. 2.3.

Yarn dimensions are measured later in the chapter. The quasi-UD textile presents a very thin weft yarn, which thickness is one order of magnitude lower than the warp. The presence of the weft yarns will mostly be neglected.

Inverse methods for geometry characterization

Inverse methods for geometry characterization

Inverse formulations in mechanics allow to compute unknown properties of the study domain starting from available response or external boundary conditions. Many kinds of such methods exist: perhaps the most common example in mechanics is the identification of material parameters [START_REF] Avril | General framework for the identification of constitutive parameters from full-field measurements in linear elasticity[END_REF][START_REF] Ma | An inverse approach for extracting elastic-plastic properties of thin films from small scale sharp indentation[END_REF] or boundary conditions [START_REF] Dennis | A finite element formulation for the detection of boundary conditions in elasticity and heat conduction[END_REF].

One possible application of inverse formulation is the detection of the unknown initial shape of a domain, when some information is available such as its deformed configuration or the boundary conditions applied. An extensive work on such application was presented in [START_REF] Yamada | Finite element procedure of initial shape determination for hyperelasticity[END_REF], where a full inverse finite element method was deployed to determine the initial shape of a rubber domain given some final conditions. Another example [START_REF] Lu | Inverse elastostatic stress analysis in predeformed biological structures: demonstration using abdominal aortic aneurysms[END_REF] is the reconstruction of the initial shape of an aortic aneurysm from the final configuration and the pressure experienced by the tissue. One big challenge in inverse methods is that they often rely on knowing some field variables from the final configuration: for example, in a mechanical deformation problem we might need to know the Cauchy stress experienced by the deformed domain, or the displacement field on some boundaries. An extensive theoretical treatment of the inverse application of boundary values can be found in [START_REF] Bonnet | Inverse problems in elasticity[END_REF]. The actual values however are not always easy to obtain due to physical shortcomings such as difficult sensors placement or limited measurable domain.

Methodology

To reconstruct the initial geometry of the textile, an inverse methodology is envisioned: the initial shape of a deformed domain can be identified from the knowledge of the final geometry and boundary pressure. The method is applied to the reconstruction of a textile composed of many yarns, which is subject to a planar compaction state between two par- allel rigid plates. The final thickness of the domain is known, while the initial geometry (thickness, width and path) of each yarn are unknown and subject to individual variability from one yarn to another. The method is now illustrated for one generic yarn, and later will be extended to the entire textile. The following assumptions are made:

-the length of the yarn is assumed to be much larger than the cross-sectional dimensions, therefore a 2D plane strain formulation can be justified. The yarn is oriented with its major axis (width) and minor axis (thickness) as illustrated in Fig. 2.4.

-the yarn cross section geometry of a yarn is described as a power ellipse as detailed in [START_REF] Sherburn | Geometric and mechanical modelling of textiles[END_REF], with a measured shape parameter n = 0.3. Therefore identifying the crosssection geometry reduces to finding the value of thickness (h) and width (w).

-The geometry is symmetric with respect to planes passing through the center O with normals e 1 and e 2 .

-The yarn is modeled as a continuum.

Based on this premise, the stress state on the domain is only known at the boundaries, where the magnitude of the normal distributed traction vector σe 2 can be measured (on the deformed domain) by the Prescale pressure sensor. By Cauchy's theorem the admissible stress state will imply:

-e 2 σ(x) e 2 = p film ∀x ∈ Γ film (2.2)
where x is the position vector of a point in the deformed configuration, p film is the boundary pressure measured by the sensor and Γ film is the boundary subregion in contact with the pressure film. By describing the deformation behaviour by an appropriate constitutive law the stress state can be related to the deformation state.

Yarn constitutive law

Previous efforts have described the mechanical behavior of yarns under compaction. Modeling approaches can be broadly classified into micro-mechanical or continuum descriptions of the yarn. In our study, continuum approach to describe the yarn mechanical behavior is adopted.

This choice is justified by the fact that, although composed of many filaments or fibers, a yarn often behaves as a continuum body due to the closely packed small diameter filaments that form the yarn bundle and the presence of sizing which maintains its cohesion. However, it will be highly non-linear mechanical behavior, since a yarn deformation under a load is a motion of the fibers within. Dixit [START_REF] Dixit | Finite element compression modelling of 2x2 twill woven fabric textile composite[END_REF] modeled the yarn as an elastic material, transversely isotropic in directions normal to the yarn orientation and described the evolution of the elastic modulus as a function of fiber volume fraction. Recently Dharmalingam [START_REF] Dharmalingam | Evolution of single carbon and glass fibrous tow crosssections in dry and lubricated states during compaction perpendicular to the fibers[END_REF], Lectez [START_REF] Lectez | Three-dimensional mechanical properties of dry carbon fiber tows subjected to cyclic compressive loading[END_REF] and Hemmer [START_REF] Hemmer | Influence of the lateral confinement on the transverse mechanical behavior of tows and quasi-unidirectional fabrics: Experimental and modeling investigations of dry through-thickness compaction[END_REF] studied the evolution of yarn cross sections under compaction and described the material behavior with a compressible hyperelastic constitutive law. They used a variation on the Ogden [START_REF] Ogden | Non-linear elastic deformations[END_REF] strain energy density function known as hyperfoam material, that was proposed in 2000 by Jemiolo and Turtletaub [56]. This model is adopted in the present work.

The strain energy density function of the hyperfoam material is defined in terms of principal stretches λ i as:

W (λ 1 , λ 2 , λ 3 ) = 2µ α 2 λ α 1 + λ α 2 + λ α 3 -3 + 1 β (J -αβ -1) (2.3)
where the principal stretches λ i = |dx|/|dX| are the eigenvalues of the deformation gradient F , X is the position vector of a point in the initial configuration, J(F) is the determinant of the deformation gradient and α, β, µ are material parameters that need to be identified. The orthonormal basis is placed so that e 3 is aligned along the direction of the fibers and e 2 is oriented in the thickness direction out of plane. Furthermore, inextensibility of the material in the fibers direction is assumed by setting λ 3 = 1. According to classical hyperelastic material theory, the eigenvalues of the Cauchy stress tensor are derived from Eq. 2.3 as:

σ i = λ i J ∂W ∂λ i (2.4)
where Voigt notation was adopted for the stresses:

σ = [σ 11 , σ 22 , σ 33 , σ 23 , σ 13 , σ 12 ] (2.5)
As stated in the assumptions, we simplify the problem according to plane strain theory. The principal Cauchy stresses in the cross-sectional plane are:

σ 1 = 2µ αλ 2 λ α-1 1 -λ -αβ 2 λ -αβ-1 1 (2.6) σ 2 = 2µ αλ 1 λ α-1 2 -λ -αβ 1 λ -αβ-1 2 (2.7)
The yarn is free to laterally expand under compression in the e 2 direction, so we adopt the assumption that σ 1 = 0 at all locations of the domain [START_REF] Hemmer | Compaction behaviour of tow at meso-scale: experimental investigations of morphological evolution at dry, lubricated and saturated states[END_REF]. This yields from Eq. 2.6 a 2.4. Methodology

λ 1 = λ -β 1+β 2 (2.8)
and Eq. 2.7 can be simplified to:

σ 2 = 2µ α λ α+αβ+1 1+β 2 -λ -αβ-1 1+β

2

(2.9)

The coefficients α, β, µ need to be characterized for this specific material in order to exercise the model.

Material parameters identification

A set of compression tests were performed in order to identify the parameters of the model. The experimental apparatus used is shown in Fig. 2.5. It consists of an Instron 5584 machine, on which a ball joint was mounted to ensure parallelism between the top and bottom plates. A small punch-die compression assembly was used to better control the placement of the yarn. For each test, a single yarn of 5 cm of length was placed between the plates of the device and compacted at very low deformation velocity of 0.5 mm/min. Values of force, displacement and time were recorded from the machine. The cross-section transformation during the compaction was recorded by a camera placed in front of the apparatus as shown in Fig. 2.6. Measurements were taken from 8 compression repetitions. The evolution of width w and thickness h was extracted by image analysis from the bounding box of the cross sections, using the software ImageJ [START_REF] Schneider | NIH Image to ImageJ: 25 years of image analysis[END_REF]. Consequently, the values of stretch in the principal directions are obtained as λ 1 = h/h 0 , λ 2 = w/w 0 where h 0 and w 0 are respectively the thickness and width of the initial configuration (stress free). The β parameter is calculated first by fitting Eq. 2.8 to the recorded experimental measurements of λ 1 , λ 2 . The β parameter is used to determine the yarn width according to the vertical compression λ 2 to obtain a Cauchy stress curve from the measured force. An assumption of homogeneous stress state is made. It is reasonable as the geometry of the yarn is almost rectangular. Although this simplification could lead to some inaccuracy in the material characterization, relaxing this assumption would require knowing the complete strain field in the cross section, which is not practical to achieve. The raw data collected is processed and the machine compliance is eliminated to yield a stress-deformation curve for 8 yarns tested. Measurements between all samples followed a coherent trend, with some differences due to the intrinsic variability that characterizes fiber yarns [START_REF] Desplentere | Micro-CT characterization of variability in 3D textile architecture[END_REF]. From this experimental data the material parameters α, µ are calibrated by mean square error minimization. Experimental and fitted curves are plotted in Fig. 2.7. The values found for the model coefficients are: α = 44.8178, β = 0.1857, µ = 0.8341MPa.

Inverse formulation

The objective of the problem is to determine the initial yarn cross section geometry in terms of width (w 0 ) and thickness (h 0 ) of a yarn under compression, given a final thickness h and a measured boundary pressure p film . The final thickness h is an arbitrary imposed value (h < h 0 ) and corresponds to the final gap between the compression plates. In the case of composite materials, it is a design parameter determined to obtain the desired fiber volume fraction. [START_REF] Advani | Process modeling in composites manufacturing[END_REF]. Adopting the hyperfoam material law allows one to use Eq. 2.2 to formulate the problem stated below in Eq. 2.10:

λ i = argmin λ i |e 2 σ(x, λ i ) e 2 -p film | ∀x ∈ Γ film (2.10)
One simple formulation would be to use Eq. 2.10 to determine λ 2 for a given boundary pressure p film that could be measured by the mentioned pressure sensitive film. Being λ i the eigenvalues of the deformation gradient F (λ i ), the reference configuration could be determined with the pull-back operation:

X = F -1 x (2.11)
which in the chosen basis and for our simple transformation yields:

h 0 = λ -1 2 h (2.12)
w 0 = λ -1 1 w = λ β 1+β 2 w (2.13)
This approach however would only work if the stress state was perfectly homogeneous and for a simple rectangular geometry. When a yarn is compressed instead the measured stress p meas exhibits a gradient of intensity, with a maximum value on the centerline, as shown in Fig. 2.11c. Furthermore, since the final width w is unknown, Eq. 2.13 is not applicable right away. For this reason we look for some further condition to link w 0 to known quantities. The cross section geometries of a population of 8 yarns in the initial state were measured using a confocal scanner (profilometer). As already stated in the assumptions, these cross sections fit into power ellipses [START_REF] Sherburn | Geometric and mechanical modelling of textiles[END_REF]. Width (w 0 ) and thickness (h 0 ) for each sample were extracted as the bounding box dimensions, and the area A 0 of each cross section is computed numerically from the profile, as listed in Table (2.1). From this sample space the average area A = 1.08mm 2 and the standard deviation σ = 0.049mm 2 (≈ 5%) are calculated. The very low deviation allows us to assume that for this particular material the initial cross section area remains constant in the reference state even though the yarn initial dimensions are affected by variability. Adopting this assumption the width w 0 of a yarn in the reference state can be directly expressed as a function of a given thickness h 0 :

w 0 (h 0 ) = arg min w 0 |A(h 0 , w 0 ) -A| (2.14)
Using this assumption of constant initial cross section area, a finite element solution of the stress field is adopted to compute the inhomogeneous left-hand side stress term of Eq. 2.10. The compaction of the yarn is simulated iteratively in ABAQUS using an explicit solver. An initial yarn cross section geometry is generated numerically, with an initial guess for thickness h * 0 and width w * 0 that satisfies Eq. 2.14 for a fixed value of A. The yarn is assigned the hyperfoam material law with the identified parameters α, β, µ in the plane strain framework. A rigid 1D element which represents the mould is used to impose the final thickness h on the domain. From the finite element solution, the σ 2 Cauchy stress is computed, which as expected is quite homogeneous along the e 2 direction and inhomogeneous in the e 1 direction, as shown in Fig. 2.9a. From the computed stress field, the normal pressure p num (h * 0 ) at the rigid element boundary Γ film is extracted. Finally, using the stress value p film acquired from the Prescale pressure-sensitive film, the error is calculated as:

ξ = |p num (h * 0 ) -p film | (2.15)
The initial thickness h 0 of the yarn is found as: and initial width w 0 is determined using Eq. 2.14 under the assumption of constant initial cross-section area A. This optimization is performed iteratively by assigning an arbitrary initial value h * 0 , performing the analysis and saving the resulting p num (h * 0 ) to a file before starting over with a new value of h * 0 . The steps of this numerical optimization are summarized in Fig. 2.8. The interest of saving all values of p num (h * 0 ) for all possible values of h * 0 is to create a database that can be readily used to quickly evaluate Eq. 2.15 for a given h when multiple yarns are being reconstructed, each with their local value of p film . An example of p num (h * 0 ) database for different yarn reference dimensions and a given final h is shown in Fig. 2.9b.

h 0 = argmin h * 0 ξ(h * 0 ) (2.

Cross-section validation

The validation of the initial geometry reconstruction is performed by comparison with experimental measurements. A set of yarns were compressed up to an imposed thickness h and the whole transformation was recorded from a camera placed in front of the cross section as in Fig. 2.10(left). A Prescale film was also placed under the yarn, and the pressure profile Fig. 2.10(right) was used to identify the yarn geometry following the procedure just described. The geometry calculated from the pressure information using the inverse method was compared with the real yarn geometry measured from the video information for error quantification.

The results and a quantification of error are reported in the table 2.2, where w 0 and h 0 are the reference dimensions. The superscripts "meas" and "calc" indicate the measured and calculated values respectively. The imposed final thickness is h = 0.25mm. The percentage relative errors in the calculated quantities with respect to the measured ones are reported in parentheses.

The results of this section show that:

-it is possible to reconstruct with accuracy a yarn geometry in its initial configuration

Full textile reconstruction

The interest of this work is to detail and showcase how stress boundary conditions can be applied for inverse method reconstruction of more complex domains and geometries.

To showcase the potential of the technique presented in Section 3 for a single yarn, the reconstruction will be extended to a full textile which contains several yarns.

Cross sections

To address this scenario, the method is simply repeated for each yarn, since by laying the entire textile sample over a layer of Prescale film and compressing it it is possible to know the pressure field p i film experienced by every i th yarn individually. A UDT400 textile was placed on the pressure film and compressed up to an imposed final thickness h. From the resulting pressure print an area of 98×98mm, corresponding to a collection of 33 weft yarns and 27 warp yarns, was extracted to be analyzed. This frame can be seen in Fig. 2.12b.

For each individual yarn, multiple intensity measures are sampled from the pressure print along its length and in between crossover points, as shown in Fig. 2.11b. The measures are averaged length-wise to compute the curve of pressure p i film relative to the i-th yarn as shown in Fig. 2.11c. Assigning a length-averaged pressure profile and therefore a constant cross-section identification to each yarn simplifies the geometry identification process. Each pressure curve is used in Eq. 2.16 to find the corresponding yarn reference dimensions (h 0 , w 0 ) under the assumption of known initial cross-section area A = 1.08mm 2 . It should be noted that this method can allow for the cross-section to vary along the yarn length by simply repeating the identification at every sampling instead of length-averaging the pressure print. However this was not taken into account as it would be material specific and the number of identifications to carry out increases especially on a large domain.

Yarn paths

When a full textile is being reconstructed it is not sufficient to determine the yarn cross section geometries, but yarn paths are needed as well as this is important for investigating dual scale flow when these textiles are impregnated with resin [START_REF] Slade | Investigation of unsaturated flow in woven, braided and stitched fiber mats during mold-filling in resin transfer molding[END_REF]. These could easily be measured by imaging techniques, taking a picture of the textile to have realistic yarn paths as was done by Gommer et al. [START_REF] Gommer | Quantification of mesoscale variability and geometrical reconstruction of a textile[END_REF]. However since the pressure field measured by the Prescale film contains spatial information as well, the yarn paths can be directly extracted from that measurement. The complexity of this step depends on the complexity of the textile architecture. In this study, since the textile is woven material, the regions where the weft and warp yarns cross over are clearly visible as spikes in pressure since the thickness is locally higher. The yarn paths are easily measured by interpolating a spline through these regions. At the current level of implementation, the yarn paths definition is performed manually by using the ImageJ software to select and output the nodal positions C(x, y) of crossover regions on the domain as shown in Fig. 2.11a. To these in-plane coordinates the z information is added by shifting the nodes out of plane following the weaving pattern to account for the thickness of the small weft yarn, which was assumed to be non-deformable. These nodes C(x, y, z) are used for the definition of each yarn path when the textile geometry is generated.

Geometry generation

Following the method proposed, the information of initial cross section and yarn path for each individual yarn has been obtained. In order to generate the 3D geometry, the textile geometry generator TexGen [START_REF] Lin | Modelling and simulating textile structures using TexGen[END_REF] is used through its Python API.

The original textile, compaction pressure map and reconstructed geometry can be seen in Fig. 2.12.

It is clearly visible that the yarn paths and geometries are not idealized, but exhibit instead the variability in cross-section and path that was captured from the pressure print.

Conclusions

In this chapter we have presented a novel methodology to identify the meso-scale geometrical variability in a textile. The method is based on the acquisition of boundary values of pressure experienced by the material in the process, which was achieved using a pressure sensitive film. The hyperfoam material model was used to relate the pressure to the deformation of the yarns. Validation carried out yarn-wise showed a good match between the reconstructed and measured geometry, proving that the sensitivity of the film is satisfactory even on small scales. An example of generation of a digital model of a textile was provided, by using the pressure film to extract the yarn paths as well as the cross sectional geometries of each yarn. Although only the planar scenario was investigated, the same methodology could be potentially extended to single-curvature shapes as well (but not double-curvature shapes).

This method has some potential as a non-destructive measurement technique for single layers of textiles. However it is important to state that interpreting the geometrical features of a single layer of textile has very little application interest. Furthermore, the methodology presents some unfavourable aspects: the treatment of the data yarn-wise is time consuming and requires the characterization of certain material parameters. Finally, the high geometrical detail measured (yarn-wise mesoscale geometry variability) hardly plays any relevant role in the upscaled properties of a multi-layer textile preform.

With these considerations in mind, we move on to the next chapter: the multi-layer scenario will be tackled and the mechanical analysis will be abandoned in favour of a more appropriate methodology.

Chapter 3

Multi-layer geometrical study

The material presented in this chapter has been published as a journal article: Bancora, S., Binetruy, C., Advani, S., & Comas-Cardona, S. (2021). Characterization of mesoscale geometrical features of a preform using spectral Moiré analysis on pressure print. Composites Part A: Applied Science and Manufacturing, 106608.

The body of the chapter has been maintained as the published text, only minor additions have been made.

Introduction

The objective of this chapter is to extend the acquisition of geometrical data to multi-layer stacks. In the light of proposing a method which is both cost-effective and applicable to industrial preforms, the Prescale pressure film is once again used as the measurement instrument. However, as interactions between layers become more complex, the fine mechanical analysis used in chapter 2 must be abandoned in favour of some other technique of investigation. A novel methodology of characterization is proposed: we take advantage of the periodic structure of textile layers to apply Moiré theory of interference to their interactions. In this theory, quasi-UD textile layers are seen as gratings and the patterns they exhibit on the pressure field as Moiré fringes. Through the power of image spectral analysis, the precise orientation and spatial distribution of yarns for each layer can be extracted from the pressure print.

To the best of the author's knowledge, this is the first time that Moiré theory is applied to the reconstruction of the mesoscale structure of stack of fabrics through the analysis of a pressure field. This technique has application potential because the information extracted can be used to generate realistic digital models of a multi-layer stack.

Outline of the chapter

This chapter relies on the theory of spectral Moiré analysis. Before applying the methodology to the task at hand, a brief overview of the theory is provided.

Next, we present a series of test cases in which stacks of multiple textile layers are compacted over a Prescale pressure film. The pressure field is analyzed using spectral Moiré theory to detect orientation and spatial distribution of each layer.

Finally, we provide an example of exploitation of the geometrical information acquired by generating a digital reconstruction of one of the layups analysed.

Materials

The textile used in this chapter is the quasi-UD carbon fiber textile (C-Weave UT710P 50K T13 from Chomarat). The material is shown in Fig. 3.1a. This fabric is morphologically very similar to the UD400 used in chapter 2, with the only key difference that the size of the yarns is fairly bigger. This shift in material is justified by technical limitations which are related to the larger size of the domain of investigation:

-a larger compression plate is used in this part, resulting in the possibility of studying larger domains but also in higher possibility of misalignment between compression plates;

-the larger area to acquire numerically by scanning limits the resolution of detail that can be appreciated:

A larger, thicker material compensates both these criticalities.

Similarly to what was done in chapter 2, the 3D surface measurement of one layer of uncompressed textile was acquired using a confocal microscope to provide an accurate dimensional reference (Fig. 3.1b). Average measurements from a sample size of 5 measurements yield the following values. Tow thickness: h meas = 1.3 mm, tow width: w meas = 4.80 mm and tow center-lines period: T meas = 5.12 mm. The geometry of the weft yarn, which is one order of magnitude smaller than the carbon tows, is neglected completely to limit the complexity of the problem.

For the measurement of the pressure field, a Prescale pressure film type LLLW (0.2 -0.6 MPa) is used.

Theoretical background: Spectral Moiré analysis

Moiré effect

Moiré patterns can be observed when two or more periodic structures are superimposed. A geometrical treatment of the patterns that are generated from superposition of periodic gratings can be found in [START_REF] Post | High sensitivity moiré: experimental analysis for mechanics and materials[END_REF].

Example of classical periodic structures which are concerned by the Moiré phenomenon are gratings, grids and dot-screens (Fig. 3.2).

In the present study we focus on gratings-like structures, as they resemble closely the yarn arrays of textiles commonly found in composite layups. We define a raised cosinusoidal grating by the function:

g(x, y) = 1 2 cos[2πf (x cos θ + y sin θ)] + 1 2 (3.1)
where f represents the frequency per unit length and θ defines the normal orientation of the grating. This spatial function assumes values between 0 and 1, therefore can easily be represented as a grayscale image bounded by its extrema (0 is black, 1 is white). In Figs.

3.3a, 3.4a two such gratings g 1 and g 2 are represented, over a unit domain, with spatial frequency f = 20 and orientations θ 1 = 0 • and θ 2 = 20 • respectively. Moiré effects arise when gratings are superimposed, which can be carried out in two ways: multiplicative superposition or additive superposition. In this work we focus only on multiplicative superposition, which for two gratings g 1 and g 2 is expressed as: The image G(x, y) resulting from this operation is shown in Fig. 3.5a. One can notice visually that a series of low frequency periodic bands arise from the interaction of the base arrays. We will refer to these bands as moiré fringes.

G(x, y) = g 1 (x, y) • g 2 (x, y) (3.2)

Spectral Moiré analysis

Fourier theory of spectral analysis can be used to characterize Moiré structures. In the frequency domain, spatial patterns are easy to interpret as they can be decomposed into a series of cosinusoidal functions. An extensive treatment of the topic can be found in [19]. Theoretical material on the subject can also be found in [START_REF] Amidror | The role of Fourier theory and of modulation in the prediction of visible moiré effects[END_REF][START_REF] Patorski | Fourier domain interpretation of real and pseudo-moiré phenomena[END_REF]. The procedure of analysis of Moiré patterns using spectral analysis is now detailed.

The first step consists in performing the 2D Fast Fourier Transform (FFT) on the greyscale image G(x, y). The resulting spectrum is shifted to center the zero-frequency. The result of the operation is the spectrum F G shown in Fig. 3.5b, where the central portion of the image which contains the frequency information has been magnified for the sake of clarity. The spectrum contains a series of impulses, each characterized by a position vector and an amplitude. Each impulse is paired with an identical twin impulse which is symmetrical from the origin. The position vector f of an impulse is expressed in polar coordinates by its distance f from the origin and the angle θ the vector makes with the u axis. These values match the frequency and orientation of the periodic grating g associated to that impulse. Every impulse is also characterized by an amplitude. Position and amplitude of each impulse are dictated by the convolution theorem [START_REF] Bracewell | The Fourier transform and its applications[END_REF] which states that the Fourier transform of the product of functions is the convolution of the Fourier transforms of the individual functions.

F G (u, v) = [F 1 * F 2 ](u, v) (3.3)
The convolution of the spectra F 1 and F 2 can be geometrically defined by placing over each of the impulses of F 1 a centered copy of F 2 or vice-versa [START_REF] Rosenfeld | Digital picture processing[END_REF]. As a consequence, we can identify two groups of impulses in F G . First we can locate the impulses that correspond to the generating layers, that we will name first order (O1) impulses. Their position is described by the corresponding frequency vectors:

f 1 at polar coordinates (20, 0 • ): identifies the grating g 1 , see Fig.

f 2 at polar coordinates (20, 20 • ): identifies the grating g 2 , see Fig. 3.4

We can also locate two new pairs of impulses which frequencies were not present in the original gratings. As a consequence of the convolution theorem the location of these impulses is obtained as the sum and difference of the frequency vectors of the base impulses being convoluted. Their amplitude is the product of the amplitudes of the impulses participating in the convolution. We will name these second order (O2) impulses:

-a set of low frequency impulses, at location f 1 -f 2 and f 2 -f 1 -a set of high frequency impulses, at location f 1 + f 2 and -f 1 -f 2

Index notation

A convenient index notation exists that allows to identify every impulse by expressing it as a linear combination of all the n generating impulses, and collecting the integer coefficients in the form (k 1 , ..., k n ). In Table 3.1 are listed all the impulses present in the spectrum F G of G, as well as their geometric coordinates and their indices. We define the O2 impulse (1,-1) and its twin (-1,1) as low Moiré because they generate the low frequency structures known as Moiré fringes that are clearly visible on G. On the other hand, the O2 impulse (1,1) and its twin (-1,-1) are high Moiré because they contribute to G with some high frequency patterns which are hardly distinguishable at visual observation [START_REF] Amidror | The role of Fourier theory and of modulation in the prediction of visible moiré effects[END_REF]. Finally the central impulse, traditionally named DC, has a frequency value equal to 0, undefined θ and represents the average intensity of G.

f 0 0f 1 + 0f 2 (0,0) - (0,0) (0, 0 • ) 0.25 DC f 1 1f 1 + 0f 2 (1,0) O1 (20,0) (20, 0 • ) 0.125 g 1 -f 1 -1f 1 + 0f 2 (-1,0) O1 (-20,0) (20, -180 • ) 0.125 g 1 f 2 0f 1 + 1f 2 (0,1) O1 (19,-7) (20, 20 • ) 0.125 g 2 -f 2 0f 1 -1f 2 (0,-1) O1 (-19,7) (20, -160 • ) 0.125 g 2 f 1 -f 1f 1 -1f 2 (1,-1) O2 (1,7) (7.07, -81.8 • ) 0.0625 low Moiré -f 1 + f -1f 1 + 1f 2 (-1,1) O2 (-1,-7) (7.07, 98.2 • ) 0.0625 low Moiré f 1 + f 1f 1 + 1f 2 (1,1) O2 (39,-7) (39.6, 10 • ) 0.0625 high Moiré -f 1 -f -1f 1 -1f 2 (-1,-

Methodology

Spectral Moiré analysis can provide information about the structure of multiple fabric layers by examining its compaction pressure field imprinted on the Prescale film. To do so, we take advantage of the similarity between layers of quasi-UD textile and periodic gratings. The underlying assumption is that the local thickness of one layer of textile roughly resembles a two-dimensional cosinusoidal field (grating) of the form:

r(x, y) = h 2 (cos[2πf (x cos θ + y sin θ)] + 1) (3.4)
where h is the maximum tow thickness, θ is the orientation of the grating and f is its spatial frequency (f = T -1 where T is the tow center-lines period). When multiple layers of textile are stacked and compressed, Moiré patterns appear due to the interaction between overlapping tows. To investigate this behavior a series of tests were carried out. Multiple stacked plies of textile were compressed using an Instron 8805 hydraulic press, as to replicate the compaction phase in preform manufacturing. One Prescale LLLW film was placed under the preform to capture the pressure field exerted on the bottom plate. Since the range of pressure sensitivity of the film is limited, for all tests in this chapter the compression is driven by imposing an appropriate force equal to 30 kN. The compression is carried out at a very low speed of 0.5 mm/min to reduce viscoelastic effects. Force/displacement data is recorded, from which the final thickness of the compressed preform can be obtained. Initial thickness (h 0 , measured at 0.1 kN force) and final thickness (h f , measured at 30 kN force) of all specimens are provided in Table 3.2.

For each test presented, the compression is carried out over a 30 × 30 cm plate and textile samples are 47 × 47 cm in size. The Prescale film used to capture the pressure field is placed at the center of this area and has dimensions of 21 × 27 cm. Each pressuremarked film is scanned and digitized using an Epson Workforce DS-5500 unit. To obtain consistent measurements, exposure and acquisition resolution (600 PPI) are manually set. The methodology is first illustrated for a layup of two angled plies, then extended to more layers. A complete list of all test cases is provided in Table 3.2.

Two Ply Stack

Two layers of C-Weave textile are stacked on top of each other at different orientations. The first layer (bottom) is laid down along an orientation which is assumed as reference (0 • ). The second layer is placed on top of the first one at some angle which is intentionally not measured, in order to create an unbiased generic layup. Care was taken to ensure the textiles were not sheared during manipulation. An LLLW Prescale film is placed under the stack to record the pressure field exerted on the bottom plate. The two-layer stack is compressed over the pressure film by the imposed force. The force-stretch compression curve is shown in Fig. 3.6b, from which we can also measure the final thickness (λ = h(F )/h 0 , where h 0 is the initial specimen thickness measured at a force of 0.2 kN). The pressure field retrieved after the compression is shown in Fig. 3.6a. From a visual inspection of the image, one can notice the presence of bands in the pressure field, crossing from left to right with a rather uniform period.

These are the Moiré fringes that were introduced in Section 3.3.1. The similarity between the pattern in Fig. 3.6a and the one in Section 3.3.2, Fig. 3.5a obtained from ideal cosinusoidal functions can be easily identified.

The Moiré patterns are then analysed using spectral analysis. The grayscale pressure image is first inverted and normalized, so that 1 corresponds to white and 0 corresponds to black. Black margins are added as needed to make the image square for performance optimization. No additional filtering is performed on the image. The two-dimensional Fast Fourier Transform (FFT) of the image is generated numerically and the spectrum is shifted to center the zero-frequency impulse (DC). The resulting spectrum in frequency domain F (u, v) is plotted in Fig. 3.7, where only the magnification of the central region is shown. The spectrum is populated by precisely located impulses. For the purpose of analysis, a graphical scheme of the same region collecting the main impulses is also provided on the right side of Fig. 3.7. The amplitude values are normalized, which results in the amplitude of the DC impulse being unity. To interpret the physical meaning of each impulse we will adopt the classical index notation used in Moiré spectral analysis, which is detailed in section 3.3.2. The impulses present in the spectrum are listed in Table 3.3. The Cartesian u, v coordinates of impulses in the spectrum are expressed as pixel positions in the shifted frequency domain with the origin in the center. Polar coordinates f, θ are expressed in mm -1 and degrees respectively (clockwise rotation is positive). For each impulse in the table a twin impulse symmetrical from the origin is present, which is omitted in the list.

As detailed in section 3.3.2, we classify impulses under two groups: first order (O1) and second order (O2) impulses. O1 impulses directly represent in frequency and orientation the generating cosinusoidal gratings that approximate the textile layers. Measuring the polar components (f, θ) of these impulses provides sufficient information to determine precisely the average orientation and spacing of the tow center-lines from the corresponding layer:

-the twin impulses (1,0) and (-1,0) identify a cosinusoidal grating oriented at θ 1,0 = 0 • and of period T 1,0 = f -1 1,0 = 0.1952 -1 mm = 5.12 mm. This grating corresponds to the morphology of the first layer that was placed as the reference orientation θ = 0 • and had a period of T meas = 5.12 mm. -the impulses (0,1) and (0,-1) identify a cosinusoidal grating oriented at θ 0,1 = 12.16 • and of period T 0,1 = f -1 0,1 = 0.1958 -1 mm = 5.10 mm. This grating corresponds to the second layer, which differs from the first only in orientation.

Position (f, θ) Amplitude Source (0, 0) f 0 - (0,0) (0, 0 • ) 1 DC (1, 0) f 1 O1 (
All the remaining impulses in the spectrum belong to the group of O2 impulses. These originate from the convolution of O1 impulses [START_REF] Bracewell | The Fourier transform and its applications[END_REF]. They are responsible for the band structures that one can perceive in the pressure field, i.e. the actual Moiré fringes that arise from the interference between layers.

-the twin impulses (1,-1) and (-1,1) identify some low frequency fringes (low Moiré) that run in the direction transverse to the tows. These bands are well within the circle of visibility [START_REF] Amidror | The role of Fourier theory and of modulation in the prediction of visible moiré effects[END_REF] and are visually evident in Fig. 3.6a.

-the impulses (1,1) and (-1,-1) define some higher frequency interference (high Moiré) between the layer gratings. They contribute to some extent to the definition of the pressure field, but at the same time their effect is less evident as one's sensitivity to high frequency detail decreases significantly [START_REF] Amidror | The role of Fourier theory and of modulation in the prediction of visible moiré effects[END_REF].

To visually clarify their physical meaning, the individual pairs of O1 and O2 impulses present in the spectrum F have been isolated in Fig 3.8, where the corresponding cosinusoidal gratings from Eq.3.4 are plotted.

Validation

The validity of the spectrum interpretation is compared to the direct visual measurement of the layer orientations. A picture of each layer being laid down is taken using a digital camera with 30 mm full-frame equivalent lens mounted in a fixed position from above. Due to space constraints the camera lens is slightly off-centered and tilted from the plate's vertical axis: perspective correction is applied to the pictures to recover the orthogonal view of the layers. The picture of the layers are shown in Fig. 3.9. The measurement of the layer orientations is performed manually by measuring the direction of the tows using the software ImageJ [START_REF] Schneider | NIH Image to ImageJ: 25 years of image analysis[END_REF]. Angle values are offset so that the bottom layer is the reference θ vis 1,0 = 0 • . The orientation of the second layer measured manually was θ vis 0,1 = 11.9 • , against the value measured in section 3.4.1 of θ FFT 0,1 = 12.16 • . The value measured from the pictures shows very good agreement with the one computed from the analysis of the pressure field spectrum. This proves that the physical interpretation of the impulses is satisfactory. However the values measured from the pictures hardly represent a benchmark for error estimation: in fact, it is reasonable to say that the uncertainty in the measurement of angles from the pictures (lens distortion, manual measurement precision) exceeds the uncertainty of measurement from analysing the position of an impulse. High sensitivity to small angle differences is a characteristic of Moiré analysis [START_REF] Post | High sensitivity moiré: experimental analysis for mechanics and materials[END_REF]. Therefore the visual validation is intended more to be a confirmation of the physical interpretation of the spectrum rather than an error quantification.

Discussion on Two Ply Stack case

Some observations regarding this case study can be elaborated. The first observation concerns the first order (O1) impulses. Since the layers in the stack are made of the same textile, all O1 impulses have the same frequency component f and only differ in orientation θ. In this case we can define a theoretical orientation circle, centered in the origin and with radius equal to f , on which all impulses from the layup layers fall, their position depending only on their orientation. As a consequence any O1 impulse not falling on the orientation circle becomes an indication of some uncontrolled distortion in the tows center-lines spacing of the corresponding layer that might have occurred during the compaction phase. The orientation circle then becomes a tool for quantification of the layup deviation from the baseline configuration.

The second observation concerns second order (O2) impulses. According to the convolution theorem [START_REF] Bracewell | The Fourier transform and its applications[END_REF], the presence of O2 impulses suggests that the pressure field exhibits characteristics of multiplicative interference rather than additive. An ideal additive field would not generate second-order impulses. (Moiré fringes could still be visible, but would rather be called pseudo-Moiré [START_REF] Amidror | The role of Fourier theory and of modulation in the prediction of visible moire effects[END_REF]). However the pressure field captured by the Prescale sensor is not a mathematical construct, but rather the result of a non-linear physical in-3.4. Methodology teraction between materials. The sensitivity bounds of the Prescale sensor also introduce some thresholding on the field. Therefore it is reasonable to assume that the pressure field captured contains traits of both multiplicative and additive interference. For this reason it would be more appropriate to name this behavior pseudo-multiplicative interference. This claim is further supported by the amplitude of all the impulses. Under the assumption of purely multiplicative interference of cosinusoidal gratings, all O1 impulses should have the same amplitude. However the measured amplitudes of O1 impulses are generally lower the farther away we move from the pressure sensor. This likely happens because the contact interaction between plies is highly non-linear, and affected by local variability of the material. On the other hand, the sequential decay in O1 impulses is useful because it reveals the ordering of plies in the stack.

Fabric Stacks with Three and Four Plies

In this section the same technique is applied to more configurations with higher number of plies. This is important because in real industrial applications fabric preforms are constructed with varying number of plies to build up to the thickness of the composite. The application of the method is demonstrated here with a Three Ply Stack (3PS) and a Four Ply Stack (4PS) test cases. In each case the first ply is placed at a reference orientation (0 • ), then all subsequent layers are placed at some arbitrary, not pre-determined angles with respect to the reference layer. That means that for an N -layers case, N -unique orientations are present in the stack. For the 3PS and 4PS cases, the pressure fields, power spectra and table of impulses are provided.

-Three Ply Stack (3PS)

Compression curve and pressure field are shown in Fig. 3.10. The FFT spectrum is shown in Fig. 3.10b Impulses are listed Table 3.4.

-Four Ply Stack (4PS) Compression curve and pressure field are shown in Fig. 3.11. The FFT spectrum is shown in Fig. 3.11b. Impulses are listed in Table 3.5

It can be observed that as the number of layers increase, impulses of higher order are also present in the spectrum (O2, O3, ... ON for N layers). They arise from the convolution of lower order impulses, hence their amplitude is much lower. Their presence is barely detectable in the FFT image, however their contribution to the pressure field is negligible.

The fact that they are present confirms the multiplicative character of the pattern, however the information they contain is redundant and will not be considered any further. The interpretation of the impulses remains unchanged: O1 impulses describe the orientation and spatial frequency of the base layers, (O1+) are generated from their convolution and are responsible for the patterns on the print. To keep the analysis synthetic:

-O1+ impulses with low amplitude are ignored -O1+ impulses far outside of the orientation circle are ignored 3.5. It can be concluded that the spectral analysis of the pressure field allows one to fully detect the structure of multi-layer layups containing up to 4 distinct orientations.

Discussion on Three and Four layer Cases

It is well known that when gratings are at nearly orthogonal orientations the sensitivity of Moiré analysis is at its lowest [19]. The Four Ply Stack case shows that the technique is viable even when the angles are close to 90 • . As already observed in the 2PS case, there is a progressive decrease in amplitude for all impulses when more orientations are introduced farther away from the pressure sensor. Hypothetical reasons behind this effect were discussed in section 3.4.1. On the other hand, the decrease in amplitude for all O1+ impulses can indeed be explained by the convolution theorem. For every additional N th layer a set of additional ON impulses are generated the amplitude of which decays as it is the product of the amplitudes of all ONimpulses that contribute in the convolution. This trend is qualitatively confirmed by the amplitude of all O1+ impulses in the three and four layer case. Given that the amplitude of the impulses in the spectrum directly correlates to the pattern intensity on the pressure film, a physical maximum number N of layers can be reached so that higher order impulses (ON +) would not be detectable in the spectrum, as their contribution to the pressure field would fall under the sensitivity of the sensor. For example the LLLW Prescale sensor used in the experiments has a sensitivity range of 0. (1, 0, -1, 0)

Position (f [mm -1 ], θ) Amplitude Source (0, 0, 0) f 0 - (0,0) (0, 0 • ) 1 DC (1, 0, 0) f 1 O1 (52,0) (0.1954, 0 • ) 0.37 layer 1 (0, 1, 0) f 2 O1 (52,-5) (0.1963, 5.5 • ) 0.26 layer 2 (0, 0, 1) f 3 O1 (47,-21) (0.1934, 24.0 • ) 0.27 layer 3 (1, -1, 0) f 1 -f 2 O2 (0,6) (0.0225, -90 • ) 0.16 low Moiré (0, 1, -1) f 2 -f 3 O2 (4,15) (0.0583, -75.05 • ) 0.18 low Moiré (1, 0, -1) f 1 -f 3 O2 (4,21) (0.0808, -79.35 • ) 0.14 low Moiré (1, 1, 0) f 1 + f 2 O2 (103,-5) (0.3875, 2.77 • ) 0.10 high Moiré (1, 0, 1) f 1 + f 3 O2 (99,-20) (0.3795, 11.42 • ) 0.10 high Moiré (0, 1, 1) f 2 + f 3 O2 (
(u, v) Position (f [mm -1 ], θ) Amplitude Source (0, 0, 0, 0) f 0 - (0,0) (0, 0 • ) 1 DC (1, 0, 0, 0) f 1 O1 (52,0) (0.1956, 0 • ) 0.37 layer 1 (0, 1, 0, 0) f 3 O1 (
f 1 -f 3 O2 (0,4) (0.0150, -90 • ) 0.10 low Moiré (0, -1, 1, 0) f 2 -f 3 O2 (2,12) (0.0458, -80.53 • ) 0.10 low Moiré (1, -1, 0, 0) f 1 -f 2 O2 (2,16) (0.0607, -82.87 • ) 0.11 low Moiré (1, 0, 0, -1) f 1 -f 4 O2 (40,50) (0.2409, -51.34 • ) 0.06 low Moiré (0, 0, 1, -1) f 3 -f 4 O2 (39,46) (0.2269, -49.71 • ) 0.04 low Moiré (0, 1, 0, -1) f 2 -f 4 O2 (37,34) (0.1890, -42.58 • ) 0.08 low Moiré
and normalised amplitude, we can estimate a theoretical minimum detectable amplitude, under the assumption of purely multiplicative interference, equal to 0.0625. This would allow one to barely detect O3 impulses while O4 impulses would be undetectable. Actual sensitivity could be even lower. The maximum order N of impulses that can be detected on the pressure field could theoretically be increased by expanding the sensitivity range of the sensor. This property of the convolution theorem could be interpreted as unfavorable.

However there is a way to turn this same property into an useful feature when addressing multi-layer layups, and this property is studied in the next section.

An application case study: misalignment detection in periodic layups

We propose here an example of scenario where one property of the convolution theorem is exploited to the benefit of the technique. In the layup of LCM preforms it is common practice to follow stacking sequences which are periodically repeated, for example:

[0 • , 45 • , 0 • , 45 • ] or [0 • , 30 • , 60 • , 0 • , 30 • , 60 •
] and so on. The precise placement of each ply according to the designated orientation plays a major role in the final properties of the preform such as the local permeability and fiber volume fraction. While tools exist to assist the precise placement of every layer, a certain degree of misalignment from the prescribed values can occur. This is especially true in complex shape molds, where the distortion of the textiles being laid down can hinder the precise alignment. Obtaining a precise quantification of the deviation from the theoretical layup is not an easy task, especially when the misalignment originates during the compaction of the plies inside a closed mold.

An application case study: misalignment detection in periodic layups

Table 3.6: 4PPm: impulses in the spectrum (Fig. 3.13a)

Index Order Position (u, v) Position (f, θ) Amplitude Source (0, 0, 0) DC (0,0) (0, 0 • ) 1 DC (1, 0, 0) O1 (50÷52,0÷1) (0.1934, 0 • ) 0.22 layers 1,3 (0, 1, 0) O1 ( 49

Four Ply Periodic Stack

The methodology presented in this article can be applied to this scenario to evaluate the alignment between plies. To demonstrate this, the compaction of a Four Ply Periodic (4PP) layup, containing 4 plies in sequence [0 • , θ 2 , 0 • , θ 2 ], was repeated twice under different conditions. In the first scenario (4PPm) the stacking was performed allowing some slight misalignment between the layers pairwise in the sequence. In the second scenario (4PPa) great care was observed in aligning the layers two by two as precisely as possible. To keep the study as generic as possible, the angle θ for layers 2 and 4 was not pre-determined. This means that layer 2 was arbitrarily placed over layer 1, then layer 4 was placed to match precisely the orientation of layer 2. In test case 4PPm, the slightly misaligned layers generate the fairly structured pressure field shown in Fig. 3.12a. One can notice visually the presence of two distinct low-frequency Moiré fringes on the image. After computing the FFT, shown in Fig. 3.13a, it is possible to detect three distinct pairs of O1 impulses that fall on the orientation circle. That is due to the fact that layers 2 and 4 are misaligned and are individually detectable as separate impulses. Furthermore the (1,0,0) impulse appears blurry and its position is not well defined but spans several pixels. Values and detected orientations are provided in Table 3.6. On the other hand in test case 4PPa, when the plies are pairwise aligned with precision, only one set of Moiré fringes is present in the pressure field in Fig. 3.12b. The spectrum in Fig. 3.13b appears similar to the Two Ply Stack case addressed in section 3.4.1. Only two O1 impulses are detectable: (1,0) and (0,1) mark the position of layers 1,3 and 2,4 respectively. Only two pairs of O2 impulses are detected. Their amplitudes are higher than in test case 4PPm. Period and orientation of all layers can be accurately determined. Values are provided in Table 3.7. The multi-layer periodic layup scenario allows one to take advantage of a property of the convolution theorem: under the assumption of pseudo-multiplicative interference, the amplitudes of two impulses are summed when the convolution of the spectra places them in the same position in the frequency domain. This applies to impulses of all orders. This condition is met only when two layers are well aligned (they behave as one single layer), as in test case 4PPa, but not in test case 4PPm. An important consequence is that all impulses in the misaligned scenario exhibit a lower amplitude than in the well-aligned

Nine Ply Periodic Stack

To confirm this finding, the same test was repeated for a Nine Ply Periodic (9PP) test case, with 9 layers stacked in periodic sequence

[0 • , θ 2 , θ 3 , 0 • , θ 2 , θ 3 , 0 • , θ 2 , θ 3 ].
In test case 9PPm the layers were aligned roughly, while in test case 9PPa the layup was realised with extreme care, in order to effectively have only 3 unique orientations. Pressure fields for these two scenarios are shown in Figs. 3.14a, 3.14b. The corresponding FFT spectra are shown in Figs. 3.15a, 3.15b. In the spectrum of test case 9PPm a total of 5 O1 blurry impulses can be detected, with the impulses (0,0,0,1,0) and (0,0,0,0,1) much weaker in amplitude than the others. This implies that at least 5 distinct orientations were present in the stack (instead of the nominal 3). It is reasonable to believe that even more unique orientations might have been present in the layup, but the sensitivity of the technique fails to detect more than 5 distinct O1 impulses, especially because they are closely grouped together (the misalignment is slight). On the other hand in the spectrum of test case 9PPa only 3 O1 impulses can be detected, corresponding to the orientations [0, θ 1 , θ 2 ]. Their position is well defined and their amplitude is higher than the O1 impulses from test case 9PPm, denoting a good alignment of the matching orientations. The 9PPa scenario also presents a lower number of O2 impulses than 9PPm.

Conclusions on Periodic Layups

The periodic layup test cases 4PP and 9PP demonstrate that the analysis of the pressure field can provide information about orientations and spatial frequency of textile layers even in presence of a high number of plies. The results lead to two important conclusions:

-it is not the number of layers, but rather the number of unique orientations that defines the limit of the technique.

An application case study: misalignment detection in periodic layups

Table 3.8: 9PPm: impulses in the spectrum (Fig. 3.15a) Table 3.9: 9PPa: impulses in the spectrum (Fig. 3.15b) -for a given periodic layup, it is possible to assess the misalignment of layers by looking at the number and amplitude of the impulses present in the FFT of the pressure field. For example, a threshold value could be set on the amplitude of O2 impulses (which decay faster than O1 in case of misalignment) to assess if the layup is acceptable or not.

Index Order Position (u, v) Position (f, θ) Amplitude Source (0,0,0,0,0) - (0,0) (0,0 • ) 1 DC (1,0,0,0,0) O1 (
Index Order Position (u, v) Position (f, θ) Amplitude Source (0, 0, 0) - (0,0) (0,0 • ) 1 DC (1, 0, 0) O1 (51,0) (0.1919,0 • ) 0.44 layers 1,4,7 (0, 1, 0) O1 ( 48

Digital preform reconstruction

The acquisition and analysis of the compaction pressure field allows one to precisely measure certain quantities of the preform layup, notably the layer orientations and tow centerlines distance. This information is valuable for several applications, most notably the characterization of mechanical or flow related properties. One possibility of exploiting this information consists in reconstructing a realistic, three-dimensional digital geometry of the layup, of which we will give here an example.

We provide an example, using TexGen to generate the geometry by scripting via its Python API. We demonstrate the process using the 4 layer case as a representative case, but the same could be accomplished for any other case. The layup digital geometry is generated using the data measured from the real material:

-Layer orientations and tow centerline spacing are assigned the values obtained from the Four Ply Stack (4PS).

-Tow cross-sectional geometry is approximated by an ellipse. The tow width is set to the quantity measured in section 3.4.1 (w = 4.9mm). Tow thickness of the compacted preform is assumed to be homogeneous across the layers, and is obtained by dividing the final thickness (Table 3.2) by the number of layers. For this test case, h = h f /4 = 2.67/4 = 0.67mm).

-Size of the domain is set to the same values of the pressure field print (20.5 × 26 cm).

The geometry generated is shown in Fig. 3.16. Layers have been clipped in visualisation to display the layers below. One can see that the tow orientations and distribution are very accurately matching the ones of the real material, while on the other hand the tow cross-sectional shapes are only approximated. Once the information about accurate layer orientation and tow centerlines position is secured, an additional step would be to refine the geometry of tow cross sections. To this extent, mechanical models exist that compute the evolution of tow cross sections under compression [START_REF] Bancora | Inverse methodology as applied to reconstruct local textile features from measured pressure field[END_REF][START_REF] Grishanov | The simulation of the geometry of two-component yarns. Part I: The mechanics of strand compression: Simulating yarn cross-section shape[END_REF][START_REF] Dharmalingam | Evolution of single carbon and glass fibrous tow crosssections in dry and lubricated states during compaction perpendicular to the fibers[END_REF]. Their implementation however deserves a completely separate treatment and will not be dealt with here. The numerical mesh of a multi-layer layup such as the one in Fig. 3.16 is the typical domain input for further computations, such characterization of flow properties of the textile layup. This kind of analysis will be at the center of the next chapters.

Conclusions

In this chapter we presented a methodology to determine the orientation and alignment of a multi-layer layup configuration based on the analysis of its compaction pressure field. Tow orientations and center-line periods of layups in compacted state were measured applying spectral Moiré analysis for the first time to the compaction pressure field. The theory of Moiré notoriously ensures high precision in the detection of angles between layers when the misalignment is small. However we have demonstrated that, even in presence of high angles (4PS and 9PP cases), the method yields accurate spatial description. The use of a cost effective pressure sensitive film results in the absence of size limitations in the area to study. This represents an advantage over other mesoscale direct measurement techniques such as µ-CT.

This novel methodology opens the way to several future developments. The first regarding the material: quasi-UD textiles were used, but Moiré theory is able to interpret also interaction of periodic grids, and contributions exist in literature that apply the theory to non-periodic structures as well [START_REF] Amidror | The Theory of the Moiré Phenomenon[END_REF]. It is the author's belief that the present methodology could be extended to woven fabric preforms and detection of average tow path distortions. Another future work might be to add a mechanical model: at the current state, the mechanical information in the pressure field is not really used. It might be worth trying to complement this analysis with a mechanical model, using the actual local values of pressure from the print, reprising elements from chapter 1 [START_REF] Bancora | Inverse methodology as applied to reconstruct local textile features from measured pressure field[END_REF]. Coupling the Moiré analysis of the patterns with a mechanical model for the yarn deformation [START_REF] Grishanov | The simulation of the geometry of two-component yarns. Part I: The mechanics of strand compression: Simulating yarn cross-section shape[END_REF][START_REF] Dharmalingam | Evolution of single carbon and glass fibrous tow crosssections in dry and lubricated states during compaction perpendicular to the fibers[END_REF] could provide some insight on the evolution of yarn cross-sections as well, hence improving even more the output digital reconstructions (section 3.6). These future works will be reprised in chapter 7.

The multi-layer methodology shown in this chapter, along with the single-layer study, conclude the first part of this thesis dedicated to the acquisition of geometric data. We can consider the digital textiles generated from such data, such as the one shown in section 3.6, as the output of this first part. Starting from the next chapter, the focus will shift to exploitation of these data for flow simulations.

Chapter 4

Single-Scale Skeleton model: formulation

Introduction

In the last chapter we presented a novel technique to extract information about the structure of multi-layer heterogeneous preforms. Most notably, the two pieces of information that can be obtained are the orientation of the yarns and the spacing between them, which is a consequence of their spatial distribution. As stated, the benefit of the spectral Moiré analysis is that the domain size is no longer a limitation. In order to exploit the information obtained by Moiré analysis for material characterisation, if the size of the domain is a factor, a fitting tool of analysis is needed. We can consider the three-dimensional preform models presented in chapter 3 as the starting point for the next study. Depending on the stacking configuration, the periodicity length of such preforms could be orders of magnitude larger than the yarn dimensions. Because of the fine scale of the channel cross-sections, coupled with the large size of the domain, a full three-dimensional Stokes flow becomes a computationally intensive problem to solve [START_REF] Ballard | Stress analysis of 3D textile composites using high performance computing: new insights and challenges[END_REF]. Going back to the premises laid down in chapter 1, we expressed the need for a tool capable of solving a flow analysis over a potentially large domain. Such a tool would serve as a fitting complement to the spectral Moiré analysis. This tool is represented by the Skeleton model, which is introduced in this chapter.

Outline of the chapter

This chapter is to be intended as strongly connected to the next chapter: here the theoretical treatment and validation is laid, then in chapter 5 the application to composite materials is presented. We take advantage of the typical aspect ratio of intra-yarn channels to reduce the complexity of the flow problem. Starting from a digital model of a porous domain we will generate a Medial Skeleton mesh, which simplifies and replaces the original geometry. Then by adopting Reynolds lubrication approximation we model a steady-state flow prob-lem on the skeleton domain. The problem is solved numerically and several quantities are computed. The solution obtained is compared with a reference solution obtained by solving a steady-state Stokes flow on the original domain. Both 2D (2D Stokes -1D Reynolds) and 3D (3D Stokes -2D Reynolds) cases will be addressed. While the 3D scenario is the most appealing for real application, the 2D case is a good starting point to illustrate the methodology and investigate its limitations.

Predictive methods for permeability characterisation

We will review briefly the non-experimental methods available for characterising the flow properties of a porous medium. The quantity that is typically investigated is permeability, both in its full tensorial form or simply in a given direction of interest. In order to go from the knowledge of a porous geometry to the characterisation of its permeability, we can identify two approaches for predictive procedures. The first group is represented by analytical solutions, which can directly express the flow properties of a porous domain given its geometrical parameters. The advantage of these methods is that a simple analytical formula can directly yield a quantification of the properties of the geometry, often to a very good degree. This is the case of Gebart's expression for permeability of cylinder arrays in different packing configurations [START_REF] Gebart | Permeability of unidirectional reinforcements for RTM[END_REF], or the Kozeny-Carman equation [START_REF] Kozeny | Uber kapillare leitung der wasser in boden[END_REF][START_REF] Carman | Fluid flow through granular beds[END_REF]. The downside of these models, however, is that they are either only applicable to very simple domains, with certain idealised geometrical characteristics, or they depend on empirical parameters that are unknown. The second group of characterisation procedures consist in solving (numerically) a flow problem on a porous domain, so that material properties can be computed from the flow solution. Typically these methods involve generating a digital representation of the porous geometry under study, on which a flow problem is modeled and numerically solved. All numerical methods provide a solution which is based on the arbitrary geometry provided, which is not required to observe restrictions in its complexity. This is why numerical flow methods are today the most viable option to perform the digital characterisation of any given porous material. The downside associated with these methods comes from the numerical solution of a (potentially) complex problem.

Analytical methods

Kozeny-Carman equation

The earliest constitutive model that linked permeability to volume fraction of a porous domain was developed by Kozeny [START_REF] Kozeny | Uber kapillare leitung der wasser in boden[END_REF] and modified by Carman [START_REF] Carman | Fluid flow through granular beds[END_REF], into what is generally known today as the Kozeny-Carman equation [START_REF] Scheidegger | The physics of flow through porous media[END_REF]. Initially developed for soils, the version of the equation used with composite materials is:

K = α (1 -V f ) 3 V 2 f (4.1)
where V f is the area fraction of the domain occupied by the cylinders and α is an empirical parameter depending on the fibrous structure of the domain, and linked to the (unknown) flow tortuosity. The presence of an empirical parameter requires one to obtain this value by other means such as conducting experiments: more authors have attempted to derive variations on Eq. 4.1 with the aim to relax the need for empirical parameters. This was achievable only by assuming some regularity and idealisation on the cell model geometry. Among others, Bruschke and Advani [START_REF] Bruschke | Flow of generalized Newtonian fluids across a periodic array of cylinders[END_REF] used lubrication theory to express permeability across an array of cylinders:

K = 1 3 (1 -L 2 ) 2 L 3   3L tan -1 1+L 1-L √ 1 -L 2 + L 2 2 + 1   -1 R 2 (4.2)
where L 2 = 4V f /π and R is the radius of the cylinder. Expression 4.2 has the benefit of not needing any empirical parameter to be fitted experimentally.

Gebart model

Gebart [START_REF] Gebart | Permeability of unidirectional reinforcements for RTM[END_REF] proposed an analytical formula to express the permeability of a domain populated by an array of cylinders. The way this approximation is derived assumes that most of the viscous dissipation happens at the narrowest points of the channel restrictions. The permeability in the longitudinal and transverse direction is expressed as:

K = 8R 2 A (1 -V f ) 3 V 2 f (4.3) K ⊥ = B V f max V f -1 5/2 R 2 (4.4)
where V f max is the maximum volume fraction possible, obtained when the cylinders touch, and A, B are geometrical parameter listed in Table 4.1. The parameters A, B are only known for simple packing arrangements, and become once again empirical parameters if the packing of cylinders is arbitrary. In particular, for the square and hexagonal packing scenarios: These expressions, although limited in applicability to real composite tow geometries, can be used as a term of comparison for the model we propose further on. Gebart's expressions 4.3-4.4 will be used again in Section 4.6.4.

Numerical methods

Finite Elements method

Finite element methods (FEM or FE) are general purpose methods. Although they apply to an endless variety of physical scenarios, their working principle is always the same: a system of PDEs of interest are solved numerically by discretising the continuous space into finite elements. In the case of fluid flow, the equations of interest are Navier-Stokes (or Stokes) and continuity. The differential equations are expressed in the weak form, which is an integral form of the same that relaxes the continuity requirements of the derivatives. The characteristic of FEM is that every variable, including the solution variables, can be expressed as a linear combination of basis functions and values of the variable at the nodes. One strong advantage of FEM is that a great number of solvers exist for the numerical solution of most equations. In the field of composite processing, this method can be used to compute the permeability or the transient filling of a part. For example there are many commerial codes available that use FEM to simulate the RTM process such as PAM-RTM by ESI [START_REF] Trochu | Numerical analysis of the resin transfer molding process by the finite element method[END_REF] or LIMS developed by University of Delaware [START_REF] Bruschke | A finite element/control volume approach to mold filling in anisotropic porous media[END_REF][START_REF] Simacek | Desirable features in mold filling simulations for liquid composite molding processes[END_REF]. In this thesis, FEM is used in chapters 4, 5 and 6.

Finite Volume methods

Finite Volume (FV) methods are extensively used to solve CFD problems. The reason is, the formulation of FV methods is suited to solve conservation equations, such as mass or momentum, expressed in their integral form. The domain is discretised in cells, not unlike with FEM, over which the divergence theorem is applied to transform volume integrals of quantities into surface integrals of fluxes across the boundaries of the cells. This operation, unlike when using shape functions, does not introduce approximation in the conservation of quantities, making FV methods valuable for filling simulations. Commercial softwares based on FV method which are used in composite materials simulations are GEO-DICT (Math2Market, Germany) and OpenFOAM (OpenCFD, ESI).

Meshless methods

Meshless methods stand out in a category by themselves because unlike all previous methods they are not based on a discretisation of space into elements. Two particle-based methods are worth mentioning: Lattice-Boltzmann method and Smoothed Particle Hydrodynamics. In the Lattice-Boltzmann method (LBM) the macroscopic behaviour of a fluid is retrieved by tracking the evolution of particles of fluid (or rather, their distribution). This rather advanced model is based on Boltzmann's equation for gases, which is upscaled to a solution of Navier-Stokes. Without going into the details of the method, which would require an understanding of particle dynamics, we remark that certain characteristics of LBM make this model suitable for mesoscale flow simulations in porous media, both for generic purpose [START_REF] Guo | Lattice Boltzmann model for incompressible flows through porous media[END_REF][START_REF] Pan | An evaluation of lattice Boltzmann schemes for porous medium flow simulation[END_REF] and for the estimation of permeability of composite textiles [START_REF] Belov | Modelling of permeability of textile reinforcements: lattice Boltzmann method[END_REF].

A comprehensive treatment of LBM application to CFD can be found in [START_REF] Succi | The lattice Boltzmann equation: for fluid dynamics and beyond[END_REF][START_REF] Kruger | The lattice Boltzmann method[END_REF].

In Smoothed Particle Hydrodynamics (SPH) the flow of a fluid is represented as a collection of interacting particles, each acting on the neighbours by a smoothed function of influence called kernel function. This numerical method, introduced in the 1990s [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF], has been successfully applied to mesoscale flow simulations in porous media [START_REF] Jiang | Mesoscale SPH modeling of fluid flow in isotropic porous media[END_REF] In the field of composite materials, SPH was used to simulate the resin infusion of a part in [START_REF] Su | Numerical study of liquid composite molding using a smoothed particle hydrodynamics method[END_REF][START_REF] Yashiro | Particle simulation of dual-scale flow in resin transfer molding for process analysis[END_REF] or even flow-induced deformations of a preform in RTM processes [START_REF] Comas-Cardona | A generic mixed FE-SPH method to address hydro-mechanical coupling in liquid composite moulding processes[END_REF]

Pore Network Model

In the Pore Network Model (PNM), a porous medium is represented as a network of pores and throats. This model is ideal for highly unstructured porous media such as soils, where it is possible to clearly identify openings in the channels network (pores) and narrow channels connecting them (throats). The advantage of PNM is that the number of unknowns to be resolved, which are linked to the number of pores, is much smaller than the number of DOFs a detailed mesh would require. Applications of PNM to flow in porous media have been countless, we refer to [START_REF] Xiong | Review of pore network modelling of porous media: Experimental characterisations, network constructions and applications to reactive transport[END_REF] for a comprehensive review. The critical feature in PNM resides in the proper generation of the pores and their connections [START_REF] Blunt | Simulation and theory of two-phase flow in porous media[END_REF][START_REF] Oren | Pore-scale network modelling of waterflood residual oil recovery by immiscible gas flooding[END_REF]. Although the applicability of the model is limited to the characteristic morphology of porous networks such as the ones found in soils, PNM has been used to model the flow in textiles [START_REF] Delerue | Pore network modeling of permeability for textile reinforcements[END_REF].

Other methods

Other numerical methods that are used to solve flow problems in porous media are the Boundary Element Method (BEM [START_REF] Um | A study on the mold filling process in resin transfer molding[END_REF][START_REF] Yoo | Numerical simulation of the resin transfer mold filling process using the boundary element method[END_REF]), the Finite Elements with Level Set method [START_REF] Soukane | Application of the level set method to the simulation of resin transfer molding[END_REF] and mixed Finite Element / Finite Difference method [START_REF] Hieber | A finite-element/finite-difference simulation of the injectionmolding filling process[END_REF]. These methods have been used to perform both saturated and unsaturated flow simulations.

Model equations

We will review now the equations which will be used in the remainder of the chapter to describe the flow behaviour. Two equations will be taken into account: Stokes and Reynolds equations. Stokes equation is used to model a creeping flow within any geometry. When the flow takes place in thin gaps and some assumptions are adopted, Stokes equation is simplified into Reynolds equation within the theory of lubrication approximation.

Stokes flow

The flow of a Newtonian fluid inside a domain is governed by Navier-Stokes equations. In their general form, assuming a constant density of the fluid, they read as:

ρ Du Dt = f + µ∇ 2 u -∇p (4.5)
where u indicates the velocity field, f is a force vector acting on the fluid volume and p is the pressure field in the fluid. The material derivative term Du Dt can be explicitly expressed as:

Du Dt = ∂u ∂t + ∇u • u (4.6)
The left-hand side terms in Eq. 4.5 is related to the inertial effects, and the terms on the right-hand side are the body forces, viscous effects and pressure effects respectively. Dimensional analysis [START_REF] Szeri | Fluid film lubrication[END_REF] shows that the contribution of inertial effects is closely related to the Reynolds number, which is defined as:

Re = ρuL µ (4.7)
where u is the magnitude of the characteristic velocity, and L is the flow characteristic length. In composite processes the resin flows is generally characterized by low velocity and high viscosity. In such flows, known as creeping flows, Re << 1 and the inertial term in Eq. 4.6 becomes negligible compared to the other terms. These simplifications lead to what is known as the Stokes momentum equation for fluids:

∂u ∂t = f + µ∇ 2 u -∇p (4.8)
and with some further assumptions such as -the flow is steady-state, ∂u ∂t = 0 -the fluid is Newtonian. Viscosity µ is constant and independent of time.

-source terms or body forces f are not present.

∇ • u = 0 (4.9)
results in the steady-state formulation of Stokes equations for fluids:

µ∇ 2 u -∇p = 0 (4.10) ∇ • u = 0 (4.11)
Given some boundary conditions, Eqs. 4.10 are typically solved numerically for any arbitrary domain using the tools described in Section 4.2.2. Besides conditions on inlet/outlet fluxes or source terms, one typical boundary condition that is usually adopted in single-scale Stokes flow is the no-slip condition. If we denote the boundaries as Γ w , this condition is expressed as: u = 0 on Γ w (4.12)

Viscous power One representative scalar quantity that can be useful for global quantification is the viscous power dissipated by the fluid. For a steady-state incompressible viscous flow, we can write an energy balance equation in the form:

Γ u • (σ • n) dS = Ω σ : D dV (4.13)
where the left-hand side represents the work done by boundary stresses and the right-hand side represents the internal viscous dissipation. The second-order tensorD indicates the symmetric part of the velocity gradient, expressed as:

D = 1 2 (∇u + ∇ T u) (4.14)
and σ is the stress in the fluid expressed as:

σ = -pI + 2µD (4. 15
)
where I is the identity matrix. Injecting this constitutive relation in the right-hand side of Eq.4.13 allows one to express the viscous power density as:

σ : D = σ ij D ij = -pD ij δ ij + 2µD ij D ij = -pD kk + 2µD ij D ij (4.16)
and since from 4.14:

D kk = 1 2 (u k,k + u k,k ) = 0 (continuity equation) (4.17)
we can express the viscous power dissipation as:

Ψ = Ω σ : D dV = 2µ Ω D : D dV (4.18)

Lubrication approximation

Lubrication theory leads to a much simplified form of Stokes equation which can describe the flow of a fluid film inside slender channels or fractures. The theory is valid for slow, creeping flows inside thin gaps or fluid films and was developed by Reynolds in 1886 [START_REF] Reynolds | IV. On the theory of lubrication and its application to Mr. Beauchamp tower's experiments, including an experimental determination of the viscosity of olive oil[END_REF].

It has been historically employed to model flow in bearings [START_REF] Dowson | A generalized Reynolds equation for fluid-film lubrication[END_REF][START_REF] Bair | High-pressure rheology of lubricants and limitations of the Reynolds equation[END_REF] and in porous rocks containing a network of fractures [START_REF] Zimmerman | Lubrication theory analysis of the permeability of rough-walled fractures[END_REF][START_REF] Ranjith | Applicability of the 'cubic law'for non-Darcian fracture flow[END_REF][START_REF] Dippenaar | On the cubic law and variably saturated flow through discrete open rough-walled discontinuities[END_REF]. Lubrication approximation can also be applied to creeping flow in textiles. The typical cross-sectional shape of inter-yarn pores is that of thin and elongated channels, in which the characteristic length of flow is much higher than the thickness. In Fig. 4.2 we show an example from the µ-CT scan of a textile, in which the inter-yarn channels are highlighted. Lubrication theory has been used as the micro-scale flow model for permeability upscaling in porous media [START_REF] Barthélémy | Effective permeability of media with a dense network of long and micro fractures[END_REF][START_REF] Levasseur | A micro-macro approach of permeability evolution in rocks excavation damaged zones[END_REF]. The model and its definition of hydraulic aperture was also used by Wong et al. in [START_REF] Wong | Comparisons of novel and efficient approaches for permeability prediction based on the fabric architecture[END_REF] to compute the upscaled permeability of textiles using a "stream surface" method. The theory can describe the steady-state flow of a fluid in a much simplified form of Stokes equations, provided that the following assumptions are adopted:

-the thickness h of the channel is small compared to the characteristic length of the flow: h/L 1 -no-slip boundary conditions at the channel walls -constant viscosity µ When these conditions are met, the pressure gradient in the transverse direction of the channel can be neglected:

∇P • e 2 = 0 (4.19)
where e 2 is the direction transverse to the channel gap.

Starting from the Stokes momentum equation 4.10 we can add the non-slip condition on velocity for a channel of gap h as: u(y = 0) = u(y = h) = 0, where u(y) denotes the component of velocity in the channel direction. Eq. 4.10 can be integrated along the gap dimension y, satisfying Eq. 4.9 and the no-slip condition, to obtain the well-known velocity profile:

u(y) = u(y) e 1 = y(y -h) 2µ ∇p e 1 (4.20)
where h denotes the height of the channel. The mass flow rate as a function of the gap thickness h is obtained directly from the velocity field by integration:

Q = ρ w 0 h 0 (u(y) • e 1 ) dydz = -w ρh 3 12µ ∇p (4.21)
where ρ is the fluid density. For simplicity in the derivation we will adopt a unit-width channel (w = 1). The gap-averaged component of velocity of the fluid in the channel direction is obtained from Eq. 4.21 as:

u = Q ρh = - h 2 12µ ∇p (4.22)
Finally we can apply mass conservation to Eq. 4.21 to obtain the differential form of the continuity equation describing flow in a thin gap:

∇ • (ρ h u ) = ∇ • ρh 3 12µ ∇p = 0 (4.23)
This equation is also known as Reynolds equation, and in the field of geoscience is sometimes referred to as the "local cubic law (LCL) [START_REF] Szeri | Fluid film lubrication[END_REF]".

Viscous power

The power dissipated by viscous force is expressed for Reynolds equation as:

Ψ = Ω u • ∇p dx = Ω u • ∇p dx (4.24)
This quantity is equivalent to the one expressed in Eq. 4.18

Validity of lubrication approximation in varying gap channels

The lubrication approximation theory, from which Eq.4.23 is derived, is based on the assumption of a slowly varying gap in the channel. The actual meaning of this statement has been largely investigated, however the conclusions that have been drawn on the rigidity of this assumption have sometimes been not univocal. While studying flow between periodically corrugated gaps, Zimmerman et al. [START_REF] Zimmerman | Lubrication theory analysis of the permeability of rough-walled fractures[END_REF] compared Stokes and Reynolds solutions for a series of channel profiles. They found that the error introduced is less than 10% if the channel gap to wavelength ratio is lower than 0.2, validating the idea that the smooth variation of the channel gap is what most affects the approximation. On the other hand Brown [START_REF] Brown | Fluid flow through rock joints: the effect of surface roughness[END_REF] found that that same ratio should be lower than 0.02 for the same result to hold. Mourzenko et al. [START_REF] Mourzenko | Permeability of a single fracture; validity of the Reynolds equation[END_REF] observed that, for a given set of fractal aperture realisations, the error introduced by lubrication approximation is not monotonic with respect to the gap aperture. Furthermore in the same study they showed that a best approximation of the Stokes results is obtained when the distance b d between the medial streamline and the channel walls, rather than the vertical separation between them b, is used as a value for gap size in Eq. 4.21, as shown in Fig. 4.3. Finally, if surface roughness and sudden discontinuities in the gap aperture are found, the validity of the lubrication approximation is hindered, as observed in [START_REF] Dippenaar | On the cubic law and variably saturated flow through discrete open rough-walled discontinuities[END_REF]. To address this issue, researchers proposed to adopt a spatially averaged definition of the hydraulic aperture [START_REF] De Kraker | Development of a texture averaged Reynolds equation[END_REF]. However, this aspect will not pose an issue in the application of the lubrication theory to textiles, where the channel geometries are generally smooth.

Medial Skeleton

By neglecting the pressure gradient in the transverse direction, Reynolds equation allows one to reduce the dimensionality of a flow problem. The gap-averaged velocity from Eq. 4.22 is constant in the cross-section leading to a Poiseuille flow. As a consequence, the direction of the flow is only determined by the direction of the hyperplane separating the channel walls. In the case of a 2D channel this hyperplane is represented by a medial line running between the walls. If the gap is three-dimensional, the hyperplane becomes a surface located at equal distance from the boundaries. We can generalise this definition to any arbitrary 2D or 3D geometry to express the concept of topological Medial Skeleton.

The Medial Skeleton of a given geometry is defined as the "locus of centers of maximal spheres bi-tangent to the shape boundaries" [START_REF] Blum | A transformation for extracting new descriptors of shape[END_REF]. The equidistant condition is defined locally, meaning that for every point belonging to the Medial Skeleton there exist at least 2 points belonging to the boundaries of the geometry which lie at the same distance.

The geometrical nature of the skeleton is dependent on the topological dimension of the original geometry: the topology of the skeleton domain has always codimension 1 with respect to the topology of the starting domain. The Medial Skeleton of a two-dimensional set of boundaries is represented by a set of curvilinear lines that occupy the locus of points equidistant from the boundaries. Such geometrical construct is known as "Medial Axis Skeleton". On the other hand, the Medial Skeleton of a three-dimensional set of boundary surfaces is built as a set of surfaces described by the points equidistant from the boundaries. Such geometrical object is known as "Medial Surface Skeleton". In the remainder of the chapter, the term "Medial Skeleton" is generically used to indicate both Medial Axis and Medial Surface Skeletons, the distinction falling on the obvious context. Whenever clarity is demanded, the proper labels "Medial Axis" and "Medial Surface" are used.

Background

Medial skeletons of shapes were historically introduced by Blum in 1967 [START_REF] Blum | A transformation for extracting new descriptors of shape[END_REF]. In this early work, he proposed several procedures to define a Medial Axis Function (MAF) of a shape (Fig. 4.4). Today a number of numerical algorithms have been developed for the extraction of medial skeletons from complex shapes. A comprehensive review of skeletonization algorithms was realised by Saha et al. [START_REF] Saha | A survey on skeletonization algorithms and their applications[END_REF]. Loosely following the classification they proposed, we can identify the main methods used to generate a Medial Skeleton mesh:

Voronoi diagram

The medial skeleton is obtained by performing operations on the elements of the Voronoi diagram of the original shape. Voronoi diagrams were formally developed by Georges Voronoï in 1908 [START_REF] Voronoi | Nouvelles applications des paramètres continus à la théorie des formes quadratiques. Premier mémoire. Sur quelques propriétés des formes quadratiques positives parfaites[END_REF][START_REF] Voronoi | Nouvelles applications des paramètres continus à la théorie des formes quadratiques. Deuxième mémoire. Recherches sur les parallélloèdres primitifs[END_REF] and have since been used in many fields. The extraction of the Medial Axis Skeleton using Voronoi diagram from a 2D shape is quite straightforward [START_REF] Brandt | Continuous skeleton computation by Voronoi diagram[END_REF][START_REF] Ogniewicz | Voronoi skeletons: theory and applications[END_REF]. The extension to Medial Surface Skeleton generation from 3D shapes is slightly more complex [START_REF] Amenta | The power crust, unions of balls, and the medial axis transform[END_REF][START_REF] Attali | Computing and simplifying 2D and 3D continuous skeletons[END_REF][START_REF] Sherbrooke | An algorithm for the medial axis transform of 3D polyhedral solids[END_REF]. The Voronoi diagram method is simple to implement and very robust for complex shapes. Furthermore, one can benefit from open-source libraries such as QHull [START_REF] Barber | The Quickhull algorithm for convex hulls[END_REF] and CGAL [START_REF]CGAL User and Reference Manual[END_REF] to generate Voronoi diagrams starting from a set of input points. The Voronoi diagram is the method that will be used to generate all the Medial Skeleton meshes, 2D and 3D, in the remainder of the thesis. Detailed procedures for the generation of Skeleton meshes using the Voronoi method are provided for both 2D and 3D in Section 4.4.3. The Voronoi method has some known downsides: the tendency of the skeleton being perturbed by small boundary features [START_REF] Ogniewicz | Voronoi skeletons: theory and applications[END_REF] and the bad scalability in efficiency with large sets of input points [START_REF] Jalba | Surface and curve skeletonization of large 3D models on the GPU[END_REF].

Grassfire method

The grassfire algorithms originate from an original concept by Blum [START_REF] Blum | A transformation for extracting new descriptors of shape[END_REF]: the original shape should be imagined as made of a mass of grass with the outer layer set on fire. The erosion of the mass would advance until the fronts of the fire would eventually meet at what is the medial skeleton of the shape. This early concept has been implemented numerically both for 2D [START_REF] Lam | Thinning methodologies-a comprehensive survey[END_REF] and 3D [START_REF] Tsao | A parallel thinning algorithm for 3-D pictures[END_REF]. Although this methodology of Medial Skeleton generation was not explored in the present work, it is that author's belief that the grassfire algorithm could be an alternative to the Voronoi method. In that case, however, additional care would need to be adopted to retrieve the distance map from the boundaries (which is easy to compute using Voronoi method).

Other methods A number of alternative methods, not further explored in the context of this work, do exist for the generation of complex skeleton meshes from 2D or 3D geometries. We mention the sequential Betti-number based algorithm [START_REF] Morgenthaler | Three-dimensional simple points: sereial erosion, parallel thinning and skeletonization[END_REF], the extended shrinking algorithm [START_REF] Jiang | Pore network extraction for fractured porous media[END_REF], and distance-transform (DT) based algorithms [START_REF] Borgefors | On digital distance transforms in three dimensions[END_REF].

Medial Skeleton and porous media

Flow models in porous media based on medial skeletons have been developed before. Riasi et al. [START_REF] Riasi | A feasibility study of the pore topology method (PTM), a medial surface-based approach to multi-phase flow simulation in porous media[END_REF] used a 3D-voxels representation (Pore Topology Method) of the medial surface skeleton of a porous domain to compute its macroscopic permeability. The permeability of rock fractures was estimated using Reynolds equation and the medial skeleton in [START_REF] Li | A pore-skeleton-based method for calculating permeability and capillary pressure[END_REF][START_REF] Jiang | Extraction of fractures from 3D rock images and network modelling of multi-phase flow in fracture-pore systems[END_REF], for 2D and 3D skeletons. The Medial Axis Skeleton is also sometimes used to identify the throats connectivity in Pore Network Modeling [START_REF] Jiang | Pore network extraction for fractured porous media[END_REF]. Aufort et al. [START_REF] Aufort | A new 3D shape-dependant skeletonization method. Application to porous media[END_REF] were able to generate very complex surface skeleton meshes from digital samples of sandstone. Random porous media structures were reduced to skeleton representations using Voronoi diagram method in [START_REF] Xiao | Geometry models of porous media based on Voronoi tessellations and their porosity-permeability relations[END_REF].

Numerical generation of Skeleton meshes

Skeleton domains are generated as discrete meshes from input geometries by using a series of numerical tools, following sequential steps. As mentioned, the Voronoi diagram method is used in the present work. The open-source library QHull [START_REF] Barber | The Quickhull algorithm for convex hulls[END_REF] is used to compute Voronoi diagrams of a set of input points through its Python API. The procedure for generating the skeleton domains is now detailed, both for the 2D and the 3D cases. The sequence of operations is fairly similar, with some additional steps in the 3D case.

Identification of Boundaries

Let Ω be a bounded geometrical domain of dimension d = 2, 3. This domain represents the full geometry on which a given flow problem will be modeled. In its discrete definition, Ω is a numerical mesh. The first step in the generation of its skeleton representation is the identification of all boundaries Γ w that correspond to solid walls of the domain, i.e. boundaries where the no-slip condition is enforced. The subset of all element nodes P w ∈ Γ w becomes the input data for the generation of the skeleton mesh. We shall refer to this set of points as generating nodes. The dimensionality of this set of points depends on the dimensionality of the input geometry, as P w ∈ R d .

Voronoi diagram

The Voronoi diagram is the mathematical tool for the identification of the points in space which will build the Medial Skeleton. Without going deep into the mathematical formalisms, a brief review of the theory is provided here.

Let a set of points P w = {p 1 , p 2 , ..., p n } ∈ R d with n ≥ 2 be given as generating nodes. d denotes the geometrical dimension of the initial geometry. For each point p i it is possible to identify a region in space such that all points belonging to this region are closer to p i than any other point of the set. We can define this region of space as:

V i = V (p i ) = {x | x -p i ≤ x -p j , i = j} (4.25)
where x ∈ R d is any generic position vector and • denotes the Euclidean distance. Each region V i is called the Voronoi region associated to the generator point p i , and the set given by is identified as the Voronoi diagram of the set of generator nodes P w . The topological dimension of a Voronoi region is identical to the geometrical dimension d of the domain.

V = {V 1 , V 2 , ..., V n } (4.
If the original geometry is defined in a two-dimensional space, each Voronoi region will delimit an area of space (Fig. 4.5a). Conversely if the original geometry is defined in a three-dimensional space, each Voronoi region will delimit a volume of space (Fig. 4.5b). We define two regions as adjacent if there exists at least one point x ∈ R d that belongs to both regions. Then for every pair of adjacent regions {V i , V j , i = j} we can define the locus of points which are equidistant from both generating points (p i , p j ). This locus of points, which is a hyperplane of dimension (d -1), identifies a boundary ∂V i,j between the regions V i , V j .

∂V i,j = {x | x -p i = x -p j , i = j} (4.27)
or, more succinctly:

∂V i,j = V i ∩ V j , i = j (4.28)
We will define the boundary ∂V i,j as a Voronoi boundary. Notice that:

∂V i,j = ∂V j,i (4.29)
By extension the set of all Voronoi boundaries is: ). The former shall be referred to as Voronoi ridge, the latter as Voronoi plane. In the case of d = 3, the definition of boundary hyperplanes can be applied recursively to the newly found Voronoi planes: these in turn are delimited once more by lines that we can consider Voronoi ridges themselves. Furthermore any triplets or more of adjacent Voronoi regions (that share one ridge pairwise) identify one intersection point v ∈ R d , that is:

∂V = {∂V 1,2 , ∂V 1,3 , ∂V 2,3 , ..., ∂V i,j } (4.
v i,j,k,...,n = V i ∩ V j ∩ V k ... ∩ V n , i = j = k = n (4.31)
We will refer to this intersection point as a Voronoi vertex and we can identify the set of all vertices as:

V = {v 1,2,3 , v 1,2,4 , ..., v i,j,k } ∈ R d (4.32)
Notice that as a consequence of its definition, every Voronoi vertex v i,j,k,...,n is equidistant from 3 or more associated generating points. We will refer to these as parent nodes to the Voronoi vertex. For example in Fig. 4.5a the vertex v 1,2,6 has the parent point p 1 , p 2 , p 6 .

Using the notion of parent nodes, it is possible to compute the Euclidean distance g between each Voronoi vertex and the associated parent nodes.

g = v i,j,k -p i (4.33)
The set of Voronoi vertices, along with the associated ridges, planes and regions, is the relevant data obtained from the Voronoi diagram. Before proceeding, it is necessary to reorganise all the data into classical ordered tables of nodes and connectivity. All Voronoi vertices are collected as nodal table of spatial coordinates:

V = {v 1 , v 2 , ..., v n } ∈ R d (4. 34
)
where n is the number of Voronoi vertices and d is the geometrical dimension of the original domain. We equally collect all distance values g and generate an array of local channel gap values h to be associated with all Voronoi vertices:

h = {2g 1 , 2g 2 , ..., 2g n } ∈ R (4.35)
Using the indices of the parent nodes of each Voronoi vertex, we can generate a series of connectivity tables. These shall be called Ck/0, where k and 0 are integer numbers referring to the topological dimension of the geometrical entities related by the table, with 0 referring to vertices. The connectivity tables vary according to the geometrical dimension of the problem.

For the case d = 2:

-C2/0: associates each Voronoi region with its Voronoi vertices. Elements: planar polygonal regions of varying number of vertices -C1/0: associates each Voronoi ridge with its Voronoi vertices. Elements: lines with 2 Voronoi vertices

Mesh construction: 2D skeleton

The Medial Axis skeleton mesh is constructed by one-dimensional line elements. These can be readily extracted directly from the connectivity table C1/0. Therefore the Medial Axis Skeleton generated from a two-dimensional geometry is simply the compound collection of all Voronoi ridges of its corresponding Voronoi diagram. The nodes of the skeleton are the Voronoi vertices. However, in order to only retain the portion of the skeleton that exists in the Ω 1 domain where the flow is defined, a filtering step is necessary to remove all unwanted regions of the skeleton. The filtering is achieved by associating a subdomain flag with each node in the skeleton. For sake of explanation we can imagine to assign the label 'c' (channel) to all nodes located inside the Ω 1 subdomain, 's' (solid) to all others. Once the nodal flags are identified, the filtering is performed by simply retaining all elements with both nodes flagged as 'c', and discarding all others. Finally the value of the local gap thickness h calculated from 4.35 is assigned to each node in the skeleton. The output is a line elements mesh Ω, representing the medial axis skeleton of domain Ω, with associated scalar field h. An example of Medial Axis Skeleton mesh is shown in Fig. 4.18b.

Mesh construction: 3D skeleton

The medial surface skeleton mesh is built as a compound mesh of two-dimensional elements. Unlike the simpler case of the medial axis skeleton, these surface elements cannot be directly extracted from any of the connectivity tables computed from the Voronoi diagram.

After the generation of the Voronoi connectivity tables, an additional step is necessary: the C2/0 table is composed of convex planar polygons of varying number of vertices. In order to generate a conventional mesh, the polygons must be subdivided in triangular elements. This is achieved by performing Delaunay triangulation on each polygon in the set. A new connectivity table C 2/0 is created, populated only by triangle elements. Finally the medial surface skeleton of a three-dimensional geometry can be generated as the compound surface of all triangulated Voronoi planes of its Voronoi diagram. Once again, domain filtering is performed in similar fashion to the 2D case: a flag 'c' (channel) is assigned to all nodes located inside the Ω 1 subdomain, 's' (solid) to all others. Only triangular elements with all nodes flagged as 'c' are retained, all others are discarded. The value of the local gap thickness h from Eq. 4.35 is assigned to each node in the skeleton. The output is a triangular

Methodology

In the remainder of this chapter, a series of test cases will be set up with the aim to evaluate the accuracy of the Medial Skeleton model. For each test case, a steady-state Stokes flow problem is modeled and solved numerically using the finite elements method. This is adopted as the reference solution. Consequently, the Medial Skeleton of the same domain is extracted and the same flow problem is modeled on the dimensionally-reduced mesh by Reynolds equation. The problem is solved once again by finite elements method, then the skeleton solution is compared with the benchmarked Stokes solution. The analysis is first carried out over two-dimensional domains: these cases are useful to isolate and study specific features of the skeleton model. In the second half of the chapter, the method is extended to three-dimensional domains: these cases have application potential. A list of the test cases present in the chapter is provided in Table 4.2.

Boundary conditions

A specific nomenclature will be adopted to denote domains and boundaries. Unless otherwise specified, this nomenclature will apply to all test cases and will only be detailed here, in order to avoid unnecessary repetitions. The computational domain is indicated as Ω.

Inlet
p = p out on Γ out Γ in u • n dS = Q in on Γ in u = 0 on Γ w (4.36)
-Reynolds:

p = p out on Γout Γin û • n dS = Q in on Γin (4.37)
It is worth noting that the no-slip condition at the solid boundaries (u = 0 on Γ w ) is already embedded in Reynolds equation by derivation. Conditions of Eqs. 4.36-4.37 describe a classical flow scenario: along the outlet boundary a known value of pressure p out is imposed, while at the inlet boundary a volumetric flow rate Q in is prescribed. By choosing this set of boundary conditions, the inlet pressure is unknown and becomes part of the solution.

Numerical solution

We will present the procedures that were followed to obtain a numerical solution of Stokes and Reynolds problems in all the test cases that follow. Before doing so, we introduce the FEM library that was used to achieve the goal.

FEniCS library

In order to solve numerically Eqs. 4.9, 4.10 and 4.23 a finite element solver must be used. For all the test cases in this chapter, the FEniCS finite element library [START_REF] Alnaes | The FEniCS project version 1.5[END_REF][START_REF] Logg | Automated solution of differential equations by the finite element method[END_REF] was used. This open-source library can be used through a convenient Python API, and presents many 4.5. Methodology benefits but also some drawbacks. The main advantage of using the FEniCS library is that only the weak formulation of the problem must be coded in the script, without the need to compute assembly elements. Another point in support is that FEniCS is well established in the research community. The main drawback of FEniCS is performance: being mostly an interpreted library, the speed and memory management of FEniCS cannot be compared to compiled codes such as Deal.II [START_REF] Arndt | The deal.II Library, Version 9.3[END_REF] or FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF]. This drawback however was of secondary importance for the present research, where the focus is not on the speed of the solution but rather on comparison between different models that are solved by the same solver.

Workstation specifications

All test cases presented in this chapter were computed on a machine running Ubuntu 18.04, equipped with an Intel Core i5-6500 CPU 3.20GHz, 16 GB of memory and an Nvidia Quadro P600 graphics card. CPU time results provided must be interpreted with respect to these specifications.

Stokes equation

To compute the solution of the Stokes flow problem we seek the unknown functions (u, p) that satisfy Eqs. 4.9-4.10 and boundary conditions 4.36. These are solved using the Finite Elements method (FEM). To do so, Eqs. 4.9-4.10 must be expressed in their weak form. Multiplying Eqs. 4.9-4.10 by arbitrary test functions v, w we obtain:

Ω µ ∇u : ∇v + Ω p ∇ • v = ∂Ω µ g v (4.38) Ω ∇ • u w = 0 (4.39)
This is the mixed formulation (unknowns p, u) classically adopted for Stokes equations. Eq. 4.38 can be solved for pressure and velocity by assigning the appropriate boundary conditions. The procedure to obtain the numerical solution of Eq. 4.38 will be the same in all the following test cases, therefore it will only be detailed here to avoid unnecessary repetitions. This applies to both 2D and 3D cases.

-The GMSH mesh file (*.msh) of the domain is imported in FEniCS using a custom written library.

-Domain Ω and boundaries Γ are identified.

-Linear elements are converted into Taylor-Hood elements of order (2,1) that satisfy LBB condition (needed for stability of mixed Eq. 4.38).

Methodology

-Boundary conditions from Eq. 4.36 are imposed. In particular, the inlet flow rate condition is enforced by means of Lagrange multiplier by adding to Eqs. 4.38:

λ * Γ in ((u • n) -q) dS + λ Γ in ((v • n) -q) dS (4. 40
)
where q = Q in /||Γ in || -Eqs. 4.38 are solved using solvers available in FEniCS. Depending on the test case, either a GMRES iterative solver or a MUMPS direct solver were used. Both these solvers can handle non positive-definite matrices that arise from the application of Eq. 4.40. The direct solver was used for most cases, except in the presence of large meshes (SR3/A,B) when memory constraints forced the choice on the solver.

-Solution validation: the value of the expression Γ in u • n dS is calculated from the solution to ensure it amounts to Q in up to a low tolerance (1E-5).

-The mesh with solution fields p, u is exported (*.xdmf format). CPU solution time (CPU t) is recorded.

-The total dissipated power Ψ is computed (Eq. 4.18).

Reynolds equation

To compute the solution of the Reynolds lubrication approximation flow problem we seek the unknown functions p that satisfies Eq. 4.23 and boundary conditions 4.37. Once again the Finite Elements method (FEM) is used. The weak formulation of the Eq. 4.23 by Galerkin method, after multiplying by arbitrary test function v and integrating over the Skeleton domain Ω, is:

Ω ρh 3 12µ ∇p • ∇v dx = Γ hgv ds (4.41)
where g = ρh 2 12µ ∇p • n. Eqs. 4.41 and 4.37 are solved using FEniCS library. This elliptical PDE is a much simpler problem than Stokes. Once again we detail only once the procedure to obtain the numerical solution of Eq. 4.41, which will be the same for all test cases.

-The GMSH Skeleton mesh file (*.msh) is imported in FEniCS.

-Domain Ω and boundaries Γin , Γout are identified.

-Boundary conditions from Eq. 4.37 are imposed. In particular, the inlet flow rate condition is assigned in the weak formulation by: g = Q in /|| Γin h|| -Eq. 4.41 is solved using a direct LU solver available in FEniCS. The low resource requirements for all the Skeleton mesh solutions did not require any particular care to the solver algorithm. -Solution validation: the value of the expression Γin û • n dS is calculated from the solution to ensure it amounts to Q in up to a low tolerance (1E-5).

-The mesh with solution fields p, û is exported (*.xdmf format). CPU solution time (CPU t) is recorded.

-The total dissipated power Ψ is computed (Eq. 4.24).

Comparison of Solutions

To assess the error introduced by the Skeleton Model when compared to the Stokes solution, we look for comparable quantities. Fig. 4.6 illustrates the method that will be followed below. Due to the topological difference between the meshes of the domains, direct comparison of solution fields p, u and calculation of standard error quantities becomes impossible.

In order to make a quantitative comparison between the models, we must look at equivalent global quantities that can be computed from both models. The quantity that will be used as a measure of comparison is the total viscous power dissipation. The expression of this quantity was detailed in Eqs. 4.18 and 4.24. To differentiate the cases, we maintain the notations Ψ and Ψ to denote the viscous power from Stokes and Reynolds equations, respectively. The Percent Relative Error ξ of the Skeleton Model solution compared to the Stokes solution is:

ξ = |Ψ -Ψ| Ψ × 100 (4.42)

Two-dimensional test cases 4.6.1 SR2/A: Flow in constant channel

The test case SR2/A, represented in Fig. 4.7, depicts the steady-state flow of a viscous fluid in a constant cross section channel. When the flow is laminar it is known as Poiseuille flow.

The flow is one-directional, and the velocity profile of the fluid is described analytically by Eq. 4.20. The existence of an analytical solution of this problem allows one to use this case to validate the numerical solvers. The geometrical and model parameters are provided in Table 4.4

L [m] H [m] µ [Pa s] Q in [m 3 /s] p out [Pa] 1 0.15 1 1 0

Numerical solution: Stokes

The domain is modeled using Gmsh. A series of meshes with decreasingly smaller size of elements are generated, in order to perform a convergence study of the solution. This step is carried out to determine an approximate size of elements to be adopted for future cases in order to have a good accuracy of the solution. Examples of a solution and a mesh are shown in Figs. 4.8 and 4.9 respectively. The solution p, u is computed for all the refined A convergence study on the root mean square error (RMSE) between the numerical and analytical velocity profiles for an increasing number of elements is plotted in Fig. 4.10b. Convergence is achieved for N e > 10000. This density of elements will be approximately

e 1 e 2 Γ w Γ w Γ in Γ out Ω u(0) = 0 u(H) = 0 u(y)e 1

Numerical solution: Skeleton

The Medial Axis Skeleton mesh is extracted from the most refined 2D mesh. The skeleton mesh is shown in Fig. 4.11. In this case, the gap size h is constant throughout the mesh. Eq.4.23 is solved numerically on the Skeleton mesh following procedure detailed in Section 4.5.2. The result is shown in Fig. 4.12. Results are reported in Table 4.5.

Results discussion

From the values reported in Table 4.5, one can notice that the error introduced by lubrication approximation is very low. This is expected because, in the case of a constant cross The aim of this test case is to test the validity of lubrication approximation in the presence of a varying gap size. We depict a scenario of a flow through a channel that contains a venturi-like restriction at half of its length. A parametric study is conducted: a series of simulations are carried out in which the wavelength λ of the restriction is kept constant, while the channel gaps at the beginning (h 0 ) and at the throat (h 1 ) are varied in order to study their effect on the lubrication approximation. The aspect ratios between h 0 and λ are shown in Fig. 4.14. This analysis is a follow-up on the studies from [START_REF] Zimmerman | Lubrication theory analysis of the permeability of rough-walled fractures[END_REF][START_REF] Mourzenko | Permeability of a single fracture; validity of the Reynolds equation[END_REF] on the validity of Reynolds equation already discussed in section 4.3.2.

SR2/B: Gradually varying gap size

h 0 h 1 Γ w Γ in Γ out Ω λ

Numerical solution: Stokes

A series of meshes for all the different variations of h 0 and h 1 are created in Gmsh. For all the cases, the mesh refinement is small enough to be comparable to the finest mesh of test case SR2/A (at least 10 elements across the smallest h 1 throat), to reasonably remove solution bias due to elements size. On each mesh, Stokes flow problem is solved following the procedure detailed in Section 4.5.

Comparison of the Results

The results from all simulations, for all the values of h 0 and h 1 tested, are reported in Table 4.6. The evolution of the percentage relative error with respect to the h 1 /h 0 ratio is plotted in Fig. 4.17. The results show that the error decreases rapidly as soon as the ratio between 4.6. Two-dimensional test cases inlet gap (h 0 ) and throat (h 1 ) is not extreme. For all scenarios, any h 1 /h 0 > 0.05 case yields an error less than 5%. This suggests, in the worst performing case (h 0 /λ = 1), that a channel gap is allowed to change its size at most by a factor of 10 per unit length. This condition is very reasonable if inter-yarn channels are considered. An interesting remark that emerges from the parametric study is the influence of the h 0 /λ ratio. When the inlet gap size is small relative to the domain length (h 0 /λ = [0.15; 0.3]), the error is monotonically decreasing with respect to the increase of the throat gap (h 1 /h 0 ). This is reasonable, as it would reflect for a lower angle of restriction (smooth varying gap assumption) and the values are in line with the ones observed in [START_REF] Zimmerman | Lubrication theory analysis of the permeability of rough-walled fractures[END_REF]. However, when the inlet gap/domain length ratio increases (h 0 /λ = [0.6; 1]), we observe a non-monotonic evolution of the error introduced by the lubrication approximation: the error quickly decreases as the h 1 /h 0 ratio becomes less extreme, but starts increasing again, following a somewhat linear trend, with Figure 4.17: Error evolution for varying gap channel increasing h 1 /h 0 . That means that there are configurations when a gradual gap variation introduces more error than a steep one. A similar effect was observed by Mourzenko et al. in [START_REF] Mourzenko | Permeability of a single fracture; validity of the Reynolds equation[END_REF], when they observed a non-monotonic behaviour of the Reynolds/Stokes error with respect to similar geometrical parameters. The interpretation of this phenomenon is not trivial, and still open to discussion, but we will try to provide some justification. As stated in [START_REF] Gebart | Permeability of unidirectional reinforcements for RTM[END_REF], the majority of viscous dissipation and pressure loss happen in the region of the maximal channel narrowing. As the h 1 /h 0 ratio increases, this region gets wider and wider, and so does the importance of the transverse velocity. As a consequence, there are two competing factors that affect the lubrication approximation accuracy: smooth gap variation and high narrowing dissipation. The former describes the smooth gap variation assumption formulated in [START_REF] Zimmerman | Lubrication theory analysis of the permeability of rough-walled fractures[END_REF][START_REF] Mourzenko | Permeability of a single fracture; validity of the Reynolds equation[END_REF] and observed in all scenarios as an error decreasing factor. The latter affects h 0 /λ > 0.6 cases, when it is not possible to identify a clear throat in the channel, where the contribution to the viscous dissipation is maximum, but instead the transverse velocity field component is present nearly everywhere and cannot be neglected. This is clearly a simple explanation, based on the interpretation of the numerical results and on the observations from [START_REF] Zimmerman | Lubrication theory analysis of the permeability of rough-walled fractures[END_REF][START_REF] Gebart | Permeability of unidirectional reinforcements for RTM[END_REF][START_REF] Mourzenko | Permeability of a single fracture; validity of the Reynolds equation[END_REF]. A more in-depth study of lubrication theory should be conducted to formulate a well defined quantification of its limitations, identify the influencing factors and provide justifications for the disagreements found in literature which is not the focus of this thesis. The sufficient information to be extracted from this test case is: the error introduced by the lubrication approximation, in presence of a smooth varying gap channel, is mostly below 10%, which becomes less than 5% if the gap variation is progressive.

There is another consideration that can be drawn from this study. By keeping approximately the same element size, the number of elements needed to mesh the 2D geometry is strongly affected by the geometrical parameters h 0 and h 1 , as can be seen from Table 4.6. On the other hand, in all the Medial Skeleton meshes extracted from these cases, the number of elements remains reasonably similar. The reason is the number of elements in the Skeleton mesh is only a function of the number of boundary nodes along Γ w , and is unaffected by the number of elements in Ω. This feature affects the CPU time needed to compute both solutions (Table 4.6), and is a key point in favour of the Skeleton model.

SR2/C: Steeply varying gap size

The aim of this test case is to check the validity of lubrication approximation when applied to a more arbitrary geometry. In Fig. 4.18a we show the two-dimensional domain used as test case. The fluid domain, denoted as Ω, is delimited by rigid wall boundaries and a rigid circular inclusion is placed at its center. Geometrical and model parameters are summarised in Table 4.7. This test case contains two critical elements which are of interest: abrupt change in gap size h and diverging/converging branches in the Skeleton mesh.

Stokes solution

The mesh of the domain is created in Gmsh. A sufficient number of elements are used to remove solution bias due to elements size. Result fields p, u are shown in Fig. 4.19 (c and d). Results are summarised in Table 4.8.

Reynolds solution

The Medial Axis Skeleton mesh is generated from the 2D mesh. In Fig. 4.18b the mesh is plotted with the channel gap field h. Once can notice the strong change, as well as the branching of the path. Result fields p, u are shown in Fig. 4.20. Results are summarised in Table 4.8.

Discussion of Results

Results show that the error introduced by the Skeleton model is less than 1% even in presence of steep variations in gap size and multiple branches.

SR2/D: Gebart's study

The objective of this test case is to test the application of the Medial Skeleton model to solve a classical problem in the domain of composite materials. Among the analytical expressions that link the geometrical parameters of a regular periodic cell to its permeability, we have

L [m] H [m] R [m] V f µ [Pa s] Q in [m 3 /s] p out [Pa]

Numerical solution: Skeleton

The medial axis skeleton mesh for each configuration is extracted from the full 2D domain, and Reynolds equation is solved following the procedure detailed in Section 4.5.

Permeability calculation

To compare the Stokes, Skeleton and Gebart models we will compute the permeability of the homogenized domain yielded by all three methods. The expressions proposed by Gebart [START_REF] Gebart | Permeability of unidirectional reinforcements for RTM[END_REF] for the transverse permeability of cylinder arrays have been presented in Section 4.2.1. We recall them here, this time in their square and hexagonal packing forms:

K ⊥quad = 16 9π √ 2 V f max V f -1 5/2 R 2 (4.43)
K ⊥hex = 16 9π √ 6 V f max V f -1 5/2 R 2 (4.44)
where R is the radius of the cylinders, V f max is the maximum volume fraction possible for a square packing configuration, obtained when the cylinders touch (V quad

f max = π/4, V hex f max = π/(2 √ 3
)). Using expressions 4.43, 4.44 the Gebart permeability for all the scenarios solved numerically can be calculated. Values are listed in Table 4.10.

Using instead the numerical solution obtained from the Stokes and Reynolds problems, the equivalent permeability of the cell in a given direction can be calculated as:

K = µQ in L A∆p (4.45)
where Q in is the volumetric flow rate going through the domain, µ = 1 indicates the fluid viscosity, L = 1 is the domain length and A is the total cross-sectional area of the domain in consideration. By using respectively values of pressure drop from Stokes and Reynolds solutions, we can compute the values of permeability yielded by solving the flow problem on the full 2D or on the Skeleton meshes. Values are reported in Table 4.10.

Discussion of Results

As expected the Gebart permeability values show good agreement with the reference solution, with the error decreasing with the increase of V f . This is coherent with the assumption underlying Gebart's derivation, as the pressure dissipation is driven by the channel narrowing.

The Medial Skeleton solution shows generally an even better agreement with the benchmark solution. This is justified by the fact that, unlike in Gebart's derivation, the whole domain geometry is taken into account in the computation of the pressure loss. Nonetheless, we can appreciate how the result is indeed dominated by the effect of the channel restriction, which is also the only region where the lubrication approximation is exact. For this reason, even the skeleton solution becomes more accurate for higher V f (narrower gap). In Figs. 4.25, 4.26 we plot a comparison between the values of permeability K obtained using the different methods. The full 2D Stokes solution is indicated as reference and is used to compute the percentage relative error with respect to the Gebart expression or the skeleton solution. The results are summarised in Table 4. [START_REF] Verleye | Permeability prediction for the meso-macro coupling in the simulation of the impregnation stage of Resin Transfer Moulding[END_REF].

As in all parametric test cases of the chapter, the computation time of the Skeleton solution is independent of the geometrical parameters of the domain, while the full Stokes solution is directly affected by it.

Finally we can remark that in the case of the hexagonal packing, all solutions show very good agreement and the error remains below 10 % in both Gebart and skeleton solutions. This can be explained by the fact that the hexagonal packing arrangement creates lower permeability configurations than the square packing arrangement for the same volume fraction. This results in a larger portion of the domain being dominated by the channel narrowings (in the hexagonal packing unit cell, there are two constrictions as opposed to only one in the square packing cell). Consequently, both the Gebart assumption and the lubrication approximation become less strong.

Stokes

Skeleton Gebart We can conclude that the Medial Skeleton model can be used to solve scenarios such as quantifying the permeability of a domain containing periodic cylinder arrays. The results show good agreement with the full solution, and proves to be equivalent for higher volume fraction and better for lower volume fractions than Gebart's solution.

Packing V f K [m 2 ] CPU t [s] K [m 2 ] CPU t [s] ξ K [m 2 ] ξ Square 0.

SR2/E-H: Gallery of test cases

To conclude the section, we present in the following page a gallery of test cases with generic geometries that can be addressed by the medial axis model. Each test case is presented with a short description, an illustration of the flow geometry, a plot of the solution and a table of parameters and results. To avoid having a large number of figures, only the pressure field p is shown. In order to provide a quick perception of the solution accuracy, we adopt a visual device: Stokes and Skeleton solutions are superimposed on the same plot, and their meshes are warped along the e 3 direction, normal to the mesh plane, by a quantity proportional to the local value of p. This gives a quick visual understanding of the local gradients.

No in-depth analysis is provided for these cases, as their features can be retraced to the test cases SR2/A-D.

Test case Geometry

Solution Results Description SR2E -L = 4 -H = 0.2 -µ = 1 -Q in = 1 -p out = 0 -ξ = 4.0 %
Long channel containing a long obstacle, where the dissipation is concentrated.

Features: strong gap variation

SR2F -L = 2 -H = 0.4 -µ = 1 -Q in = 1 -p out = 0 -ξ = 6.6 %
Channel containing multiple obstacles with irregular morphology.

Features: strong gap variation, complex geometry

SR2G -L = 3 -H = 1 -µ = 1 -Q in = 1 -p out = 0 -ξ = 6.0 %
Complex domain with many obstacles of different shapes, changes in flow direction and high local h/λ ratios.

Features: strong gap variation, complex geometry, high h/λ ratios

SR2H -L = 2 -H = 0.6 -µ = 1 -Q in = 1 -p out = 0 -ξ = 3.2 %
Hexagonal packing of ellipses, a transverse flow scenario using geometries typical of yarn cross sections.

Features: textile geometry, branching paths 4.7. Three-dimensional test cases

Conclusions

In this section a series of test cases were presented where the lubrication approximation theory was applied to the Medial Skeleton model and tested. The first set of test cases covered the two-dimensional scenarios. This analysis was useful to isolate and study critical geometrical features of domains that can affect the applicability of the Medial Skeleton model. It was found that for most of the scenarios solved, the Skeleton model yields a close solution to the reference Stokes solution in terms of pressure and dissipated power. Throughout all the test cases, the error rarely exceeded 5% and never 10%, with the exception of some extreme cases from Test Case SR2B, which are not common in the geometries of inter-yarn channels commonly found in textiles. Test cases SR2D showed that the Medial Skeleton model can accurately estimate permeability of periodic cells encountered in fabric architectures in the composite literature.

The study has confirmed some of the observations from literature, while still diverging from some others, confirming the difficulty in finding a univocal rule to quantify the approximation introduced by lubrication theory. The strength of the methodology relies in the possibility of generating a numerical mesh of a Medial Skeleton for any arbitrary geometry. The savings in computational time from adopting the Medial Skeleton Model instead of Stokes have been quantified. Although the comparison is strongly biased by the efficiency of the FEniCS solvers, which differ from Stokes to Reynolds equations, the quantification proves that the solution times differ by two orders of magnitude at the very least.

While the two-dimensional Medial Skeleton model is excellent to illustrate the method and analyse the influence of the parameters involved, it is of little interest for real applications. In the second half of the chapter, the method will be extended to the three dimensional cases.

Three-dimensional test cases

The objective of this section is to extend the Medial Skeleton model to three-dimensional geometries and assess its validity. As introduced in Section 4.4, the Medial Skeleton of a three-dimensional geometry is a surface. We adopted the name "Medial Surface Skeleton" to identify this version, although for the remainder of the chapter the name "Medial Skeleton" will often be used without risk of ambiguity. The study will follow the same methodology as in the previous Section (Fig. 4.6): a series of problems will be carried out, in which the solution of three-dimensional Stokes equations (Eq. 4.10) will be compared to the solution of two-dimensional Reynolds equations (Eq. 4.21) solved on the Skeleton mesh. The procedure to generate a Medial Surface skeleton mesh from a three-dimensional set of boundaries was detailed in Section 4.4. As in the previous part of the work, we will adopt the same nomenclature convention to identify domains and boundaries. However, in the three-dimensional case, the topological dimension of all boundaries is increased by one. The updated summary of boundary definitions is reported in Table 4.11.

Stokes solution

The domain is modeled and meshed using tetrahedral elements. For this particular test case, the stretched aspect ratio in the width dimension results in a very high number of 3D elements. Boundary conditions as in Eq. 4.36 are assigned. In order to represent the same flow as in the two-dimensional counterpart, one additional boundary condition becomes necessary, i.e. symmetry of the velocity field:

u • n = 0 on Γ sym (4.46)
where n denotes the outwards normal of the boundary facet. This condition will be true, in all future test cases, whenever boundaries marked as Γ sym are present. The large amount of memory needed by the library to obtain the solution (>8 GB) is a justifiable limitation to the possibilities of three-dimensional calculations over large domains. Solution is shown in Fig. 4.28. Results are summarised in Table 4.13.

Skeleton solution

The Medial Surface Skeleton mesh is generated from the three-dimensional domain. The Reynolds problem is solved following the procedure detailed in Section 4.5.2. Solution is shown in Fig. 4.29. Results are summarised in Table 4. [START_REF] Desplentere | Micro-CT characterization of variability in 3D textile architecture[END_REF].

Comparison of the Results

The relative percentage error is less than 5% when the viscous dissipation Ψ is considered. This result is not surprising, since the solution is the three-dimensional extension of problem SR2-C. Nonetheless the comparison shows that the Medial Skeleton model is a viable alternative to the full three-dimensional solution even in presence of sudden gap size changes. The trade off for the error due to the approximation is justified by the saving in computational time and resources. This is especially true with 3D problems.

SR3/B: Doubly slanted plates

In Fig. 4.30a we show an irregular geometry bounded by two solid plates, of which the top one is not planar. The plates identify a channel whose thickness varies in both e 1 and e 2 directions. This test case was explored with the aim of testing this unique feature of the Medial Surface Skeleton: the change in gap height h in two directions. This scenario is critical because the lubrication approximation might not hold if the velocity component normal to the skeleton becomes non negligible. We conduct a parametric study: the same flow problem is repeated for different h 1 /h 2 aspect ratios, ranging from 2 to 4.

L [m] W [m] h 1 [m] h 2 [m] µ [Pa s] Q in [m 3 /s] p out [Pa]

Stokes solution

The domain is modeled using Gmsh.

Skeleton solution

One Medial Skeleton mesh is generated from each 3D mesh of the parametric study. Boundary conditions as in 4.37 are assigned. The solution p, u for h 1 /h 2 = 4 is shown in Fig. 4.32.

Discussion of Results

The first observation we can draw is that the relative percentage error ξ committed by the skeleton model is always less that 3%, for all values of h 1 /h 2 . This value is especially low if the big difference in CPU time between the Stokes and the Skeleton solutions is considered. We remark that, as in all previous parametric test cases, the number of elements in the 3D mesh is not constant, while the number of elements in the Skeleton mesh remains approximately the same. The velocity field obtained on the skeleton mesh shows a good match with the three-dimensional counterpart for all h 1 /h 2 = 4 ratios. In Fig. 4.33 we plot the streamlines calculated from a flat slice of the Stokes velocity field (left) and from the full Skeleton velocity field (right), both extracted from the case h 1 /h 2 = 2. There is no point in comparing the velocity values because u = u, but one can appreciate the qualitative similarity of the paths. This test case proves that the Medial Skeleton is suited for application to three-dimensional problems with arbitrary gap variation in more directions.

SR3/C: Inter yarns channel

This test case is aimed at applying the Medial Skeleton modeling to the typical geometrical shape of a channel in-between the yarns. The space for the resin to flow between two adjacent yarns is typically hourglass-shaped, as can be seen for example in Fig. 4.2. We replicate a typical hourglass-shaped channel as shown in Figs. [START_REF] Parnas | Permeability characterization. Part 1: A proposed standard reference fabric for permeability[END_REF].34a and 4.34b. Geometrical and model parameters are provided in Table 4.16. This geometry is an oversimplification of real channel geometries, which are much more irregular and hardly constant throughout their length due to material variability.

L [m] W [m] h 1 [m] h 2 [m] µ [Pa s] Q in [m 3 /s] p out [Pa]

Discussion of Results

The relative percentage error of the Skeleton Model is less than 10 %. This value is fairly low, rendering the Medial Skeleton a viable model to describe flow in inter-yarn channels.

It should be remarked that the geometry is over-simplified when compared to real material channels. More realistic geometries will be presented in chapter 5.

Gallery of test cases

As with the two-dimensional scenarios, in the next page we collect some generic purpose test cases with the aim to further compare the Medial Skeleton model to the Stokes solution. These cases are generated on arbitrary geometries with different degrees of complexity, with similarity to variable gap channel flow scenarios. Each case is presented with a short description, an illustration of the flow geometry, a plot of the solution and a table of parameters and results. Only the velocity field v is represented, with a superposition of the Skeleton solution over a section view of the 3D solution. Errors observed remain below 10%, suggesting that the Skeleton model approximation is acceptable to address generic cases of three-dimensional channel flow.

SR3D -L = 2 -H = √ 3/2 -W = 1 -V f = 0.73 -µ = 1 -Q in = 1 -p out = 0 -ξ = 1.3 %
Longitudinal flow between array of cylinders in hexagonal packing (V f = 0.73). Application of the Skeleton model to a classical scenario as tackled by Gebart in [START_REF] Gebart | Permeability of unidirectional reinforcements for RTM[END_REF].

Features: classical case, longitudinal flow SR3E -L = 2 -H = 0.4 -W = 0.6 -µ = 1 -Q in = 1 -p out = 0 -ξ = 5.0 %
A complex geometry channel with no planar symmetries: the flow around a cylindrical and a spherical obstacles.

Features:

complex geometry, threedimensional velocity field 4.8. Conclusions

Conclusions

In this chapter we introduced the Single-Scale Skeleton model for both two-dimensional and three-dimensional flow problems. When a channel-flow-like scenario can be identified, the lubrication approximation assumption can be adopted to reduce the dimensionality of the flow by neglecting the gradient of pressure in the transverse direction. The Reynolds equation of lubrication can be numerically solved on a Medial Skeleton mesh, generated from any arbitrary geometry. The potential of this tool has been demonstrated both qualitatively and quantitatively by comparing the Skeleton model accuracy against the full-dimensional Stokes solution. The error quantification was performed by comparing the values of total viscous power dissipation from the concurrent models. It was found that, for all various scenarios tested, the error introduced by the Medial Skeleton model remains mostly below 10%. The computation time, on the other hand, is 2 to 4 orders of magnitude smaller.

We conclude that the Skeleton model is a viable tool to solve large scale flow problems at a reasonable approximation level and at a low computational cost when compared to full solutions. In the next chapter, the Single-Scale Skeleton model will be applied to composite material scenarios, specifically to simulate the flow through textile preforms typical of RTM processes.

Chapter 5

Single-Scale Skeleton model: application to textiles

Introduction

The test cases presented in Chapter 4 were used to demonstrate the applications of the medial skeleton model and to quantify its approximation as opposed to solving for the full flow solution. In this chapter the methodology will be applied to the main challenge at hand: modeling of resin flow in textile layups in RTM processing. The flow is still assumed to be single-scale, meaning the yarns are considered solid and the resin flows only inside the network of inter-yarn channels. A dual-scale extension of the model will be presented in Chapter 4. The potential of the medial surface skeleton for this application is two-fold. First, the characteristic cross-sectional geometry of the inter-yarn channels resembles the typical gap-varying channels presented in Chapter 4, which yielded accurate predictions when modeled by the medial skeleton approach. Second, the mesoscale 3D computation on an RVE containing a stack of layers generally requires large computational resources, due to the small size of elements required to get good accuracy of the field variables. Such a cost is greatly reduced on the skeleton mesh.

Outline of the Chapter

The objective of this chapter is to apply the Skeleton model to textile preforms, specifically to address the following questions:

-solve the saturated flow problem and compute the full permeability tensor of a digital multi-layered RVE -solve the unsaturated (transient) flow problem and simulate the single-scale flow within a textile layup at the meso-scale 93

Permeability measurement

This chapter is organized as follows. First the background and the methodology are described. Then, several case studies are addressed and discussed.

Permeability measurement

A large part of this chapter is dedicated to estimating permeability of fabric architectures encountered in composite processes using the Skeleton model. By conducting a flow simulation on textile architectures, we can compute an upscaled permeability of the full domain. Better yet, by following some well-established procedures, we can determine the full permeability tensor of any given fabric layup.

When the permeability tensor of an anisotropic material is unknown, several techniques exist that allow one to determine all its components. These procedures, either experimental or simulated numerically, typically consist in inducing a flow in the material and measuring the resin velocity or pressure along pre-determined directions.

The type of flow, boundary conditions and directions of measurement vary according to several possible experimental approaches. We can distinguish constant pressure measurements and constant flow rate measurements [START_REF] Weitzenböck | Measurement of principal permeability with the channel flow experiment[END_REF]. In the case of constant pressure experiments a pressure difference is imposed between the inlet and the outlet vents of a mould and the flow front progression or flow rate is measured over time to obtain the unsaturated (transient) permeability. Alternatively, the flow rate is monitored over time to measure the saturated (steady state) permeability. In the case of constant flow rate, the difference between the inlet and outlet pressure is measured.

On a flow type basis, we can distinguish one-dimensional channel flow, two-dimensional radial flow and three-dimensional hemispherical flow [START_REF] Advani | Process modeling in composites manufacturing[END_REF].

The different measurement possibilities enabled by different experiment types have prompted researchers to make some unification effort in the measurement protocols. Most notably, three international permeability benchmark exercises [START_REF] Arbter | Experimental determination of the permeability of textiles: a benchmark exercise[END_REF][START_REF] Vernet | Experimental determination of the permeability of engineering textiles: Benchmark II[END_REF][START_REF] May | In-plane permeability characterization of engineering textiles based on radial flow experiments: A benchmark exercise[END_REF] which revealed the need of an ISO measurement standard, currently under development.

Channel flow measurement

The channel flow measurement, sometimes referred to as the unidirectional flow measurement, is the simplest experimental protocol. The experimental setup consists of a flat rectangular mould inside which the material is placed and compressed to the desired volume fraction (Fig. 5.1). A unidirectional flow of fluid of known viscosity is induced by either imposing constant flow rate or constant pressure boundary conditions. The flow front position is tracked by either a camera or pressure transducers to measure the average velocity of the fluid. The measurement can be repeated along different known directions to determine the full permeability tensor of the material. One critical factor that needs to be managed is the risk of race-tracking effects along the edges of the mold [START_REF] Gauvin | Permeability measurement and flow simulation through fiber reinforcement[END_REF]. The channel flow measurement will be simulated numerically throughout the remainder of the chapter Figure 5.1: Schematic of channel flow permeability measurement apparatus [START_REF] Gauvin | Permeability measurement and flow simulation through fiber reinforcement[END_REF] to characterize the permeability of textiles using the Skeleton model. Therefore it is worth detailing the procedure that will be used for the characterization of permeability.

Channel flow permeability model

In Fig. 5.2 we represent the typical top view schematic of a channel flow experiment. A rectangular porous domain of length L and width W is delimited by inlet and outlet boundaries. If a pressure difference ∆p = p in -p out is imposed between the inlet and outlet regions, a flow of resin is induced that obeys Darcy's law:

u = - K µ ∇p (5.1)
where u is the phase averaged resin velocity. In the channel flow measurement it is assumed that the velocity across the width of the domain is null, that is: u • e 2 = 0. By expanding Eq. 5.1, we obtain:

u 2 = 0 = - K 12 µ ∂p ∂x 1 - K 22 µ ∂p ∂x 2 → ∂p ∂x 2 = - K 12 K 22 ∂p ∂x 1 (5.2) u 1 = - K 11 µ ∂p ∂x 1 - K 12 µ ∂p ∂x 2 = 1 µ K 11 K 22 -K 2 12 K 22 ∂p ∂x 1 (5.3)
When the resin flow front has reached a position x 1 , the average pressure gradient between the inlet and the flow front is computed under the assumption of linearity, as:

∂p ∂x 1 = ∆p x 1 (5.4)
For simplicity, we will consider that x 1 = L, meaning the resin has reached the outlet region of the mould and the material is saturated. By measuring the flow rate of the resin at that time, it is possible to compute the averaged velocity in direction e 1 :

u 1 = Q A = 1 A A u 1 dA (5.5)
where A is the cross-sectional area of the porous material.

Alternatively, if we have access to the full velocity field, as is the case with numerical simulations, the same quantity can be expressed as:

u 1 = 1 V Ω u 1 dV (5.6)
For a channel flow experiment, due to mass conservation and lateral solid walls boundary condition it is trivial to see that Eqs. 5.6 and 5.5 yield the same quantity.

Then we can express the averaged velocity of the resin as:

u 1 = K µ ∆p L (5.7)
If the permeability tensor of the material is assumed as homogeneous, then by comparison between Eq. 5.7 and Eq. 5.3 the following expression must be valid:

K = K 11 K 22 -K 2 12 K 22 (5.8)
We will refer to this quantity K as effective permeability in the direction of the mould (e 1 in this scenario).

From Eq. 5.7 we can express the effective permeability:

K = µ u 1 L ∆p (5.9)
The effective permeability quantity expressed in Eq. 5.9 can be measured easily from a channel flow experiment. However, one measurement alone is not sufficient to determine all the permeability tensor values K 11 , K 12 , K 22 . Three different measurements along three different directions are needed [START_REF] Weitzenböck | A unified approach to determine principal permeability of fibrous porous media[END_REF].

For that purpose, Eq. 5.8 needs to be made parametric with respect to the angle of measurement θ. This is easily accomplished by operating a change of basis to the permeability tensor. We can express a new permeability tensor K in an orthonormal basis e 1 , e 2 rotated by an angle θ around e 3 from the original basis as:

K [e 1 ,e 2 ] = R K [e 1 ,e 2 ] R T (5. 10
)
where R is the rotation matrix: R = cos θ -sin θ sin θ cos θ (5.11) Finally, we can re-write Eq. 5.8 of effective permeability as function of the angle θ:

K(θ) = K 11 K 22 -K 12 2 K 22 (5.15)
By performing three separate measurements of K(θ) along three distinct orientations at an angle θ with respect to the direction e 1 , it is possible to determine the values of K 11 , K 12 , K 22 of the permeability tensor expressed in the original basis e 1 , e 2 . This is accomplished on a practical standpoint by simply rotating the material inside the mould by known angles θ 1 , θ 2 , θ 3 and performing the same measurement three times.

Once the full tensor of permeability is determined in the [e 1 , e 2 ] basis, it can be diagonalised to find the principal permeability values (its eigenvalues) and the principal permeability directions (its eigenvectors).

Special cases may require less measurements:

-if the principal directions are known a-priori, then only two measurements of effective permeability along these directions directly yield the principal permeabilities -if the material is known to be isotropic, then one single measurement in any direction will suffice

Radial flow measurement

In the radial flow measurement, the porous material is placed in a transparent mould and compressed up to the desired volume fraction (Fig. 5.3). A radial flow of a fluid of known viscosity is induced by placing the inlet gate at the center of the mold. Either constant flow rate or constant pressure boundary conditions can be used. The flow front is tracked by camera recording or pressure transducers. The radial flow measurement presents some advantages over the channel flow method. First, race tracking effects are absent because the fluid does not reach the edges of the mould. Secondly, the entire permeability tensor is determined by one single experiment because the anisotropy in the material will affect the elliptical shape of the advancing flow front. On the other hand, the correct determination of average velocity and pressure gradients is not straightforward, because it is not possible to isolate flow directions and make an assumption of linear pressure gradient. As shown in [START_REF] Weitzenböck | Radial flow permeability measurement. Part A: Theory[END_REF], the inlet gate shape and radius also affect the measurement. In this chapter, radial permeability measurement is applied to one single test case presented in Section 5.6.4. We provide a brief overview of the methodology. The theory and procedure of measuring permeability from a radial flow experiment is thoroughly detailed by Weitzenbock in [START_REF] Weitzenböck | Radial flow permeability measurement. Part A: Theory[END_REF][START_REF] Weitzenböck | Radial flow permeability measurement. Part B: Application[END_REF]. Additional studies of radial flow permeability measurement can be found in [START_REF] Gauvin | Permeability measurement and flow simulation through fiber reinforcement[END_REF].

Radial flow permeability model

If the material is isotropic, the flow front expands with a circular shape centered in the inlet gate. We will consider here the scenario where a fixed inlet pressure is imposed. After expressing Darcy's law in polar coordinates and solving the continuity equation, the permeability is given by the expression:

K = r 2 f [2 ln (r f /r o ) -1] + r 2 o 1 t µφ 4∆p (5.16)
where r f (t) is the transient position of the flow front, t is time, r 0 is the inlet radius (r f (0) = r 0 ) and φ is the porosity of the fabric. By recording the flow front radius over time, one can plot the value of permeability given by Eq. 5.16 until the values converge, some distance away from the inlet. When the material is not isotropic, the flow front will develop into an elliptical shape. The orientation of the ellipse axes marks the principal direction of permeability. Furthermore, if we denote the flow front positions measured along the two ellipse axes as x f and y f , the ratio of principal permeabilities is equal to the ratio of the ellipse axes squared:

K 1 K 2 = x f 2 y f 2 (5.17)
The expressions of permeability in the anisotropic case become very complex [START_REF] Weitzenböck | Radial flow permeability measurement. Part A: Theory[END_REF]. However it is possible to derive much simplified expressions of principal permeabilities K 1 and K 2 if the measurements of the flow front position are performed along axes of the ellipse:

K 1 = x 2 f 2 ln x f /r o -1 + r 2 o 1 t µφ 4∆P
(5.18)

K 2 = y 2 f 2 ln y f /r o -1 + r 2 o 1 t µφ 4∆P (5.19)
The values of K 1 and K 2 from Eqs. 5.18-5.19 can be computed for the transient flow front positions x f and y f at different times t until values converge. Expressions 5.18-5.19 will be used to compute the permeability tensor of an anisotropic stack in Section 5.6.4. We remark that the use of r 0 as the inlet radius for the anisotropic case is a simplification to the original expression, because distinct values x 0 and y 0 would be required. These values are difficult to measure and can affect the accuracy of the calculation. The sensitivity to the inlet radius is a downside of the radial flow permeability measurement [START_REF] Weitzenböck | Radial flow permeability measurement. Part A: Theory[END_REF].

Methodology: saturated flow

The methodology that will be followed is largely similar to the one detailed in Chapter 4. A textile geometry is generated, the Medial Surface Skeleton mesh is extracted from it, and the Reynolds flow equation is solved numerically on it. There are however some differences from the methodology followed in Chapter 4 as related to 1. Textile geometry generation 2. Boundary conditions 3. Validation and will be discussed in this section.

Textile geometry generation

The digital geometry of textiles to be used for the test cases is generated using TexGen [START_REF] Sherburn | Geometric and mechanical modelling of textiles[END_REF], and not Gmsh. This change of tool is natural, because TexGen was created specifically for the creation of digital textile geometries. The TexGen library offers control over the position, path and cross section of individual yarns and is implemented in the study using its Python API scripting language. In particular, a custom-built library will be used that allows one to declare as input data:

-the number of layers -the orientation and vertical position of each individual layer -the yarn cross section shape, dimensions, and center-line spacings -the RVE size -the output mesh elements size By modifying this data, many possible layups can be created. For all the following test cases, we will be assigning the same yarn cross-sectional geometry. The dimensional parameters are provided in Table 5.1 and a schematic detailing the quantities is shown in Fig. 5.4. These are the same geometrical parameters used in test case 4PS in chapter 3. They correspond to four layers of C-Weave quasi-UD textile, compressed to λ 3 = 0.44 of their original thickness. Once again the thin warp yarns are neglected. Other parameters that affect the realism of These geometric modeling aspects should be further developed and included in the study to increase the level of realism of the generated geometries. However, measuring the variability is a more complex task, one that could be achieved only by sophisticated mechanical modeling, µ-CT scanning, or adopting the methodology of Chapter 1 which is so far limited to a single layer. All these efforts are out of scope in the current study, the latter will rely on simplified, idealistic layup geometries. An example of TexGen generated layup geometry is shown in Fig. 5.5a. The set of all yarn surfaces constitutes the Γ w boundary.

On the other hand, the extraction of the Medial Surface Skeleton mesh follows the same procedure as in chapter 4. The set of generating points is the collections of all nodes P w ∈ Γ w . The underlying assumption, embedded in the Reynolds equation, is the no-slip condition on the surface of the yarns:

u = 0 on Γ w
This condition, typical of single-scale models, will be partially relaxed in Chapter 6. The Skeleton mesh is created following the procedure detailed in Section 4.4.

Boundary conditions

The nomenclature used to label the regions of the skeleton meshes is identical to the one used in chapter 4. The Skeleton mesh domain is denoted as Ω. Inlet and outlet boundaries of the mesh are denoted as Γin and Γout , respectively. The boundary conditions that will be used in this chapter are slightly different from the ones used in chapter 4. The inlet pressure, p in , and outlet pressure, p out , are imposed as Dirichlet boundary conditions on the corresponding boundary edges. That is expressed as: p = p in on Γin p = p out on Γout (5.20) This adaptation is made to mimic typical boundary conditions applied in RTM processes, where it is common to impose a pressure difference between inlet and outlet, while the

Workstation specifications

All test cases presented in this Chapter were run on the same workstation, whose characteristics are detailed in chapter 4, Section 4.5.2.

Test cases: saturated flow

RT/A: UD aligned layers

As a first case study, we use TexGen to generate a unidirectional layup composed of 4 aligned layers of the same C-Weave UD textile used in chapter 3. Geometry parameters are provided in Table 5.2. A small gap of g = 0.05 mm is created vertically between the layer planes, in order to generate a non-interrupted skeleton mesh surrounding the yarns (yarns are not in contact). Although being idealised, the similarity between this geometry and the µ-CT scan of a similar quasi-UD material from Fig. 4.2 is appreciable. The purpose of this analysis is to apply the Medial Skeleton model to compute the permeability tensor of the material. In this case, the principal directions of permeability of the layup are known because all yarns are aligned along the same direction. The highest permeability (K 1 ) is in the direction along the yarns while the lowest permeability (K 2 ) will be in the direction transverse to the yarns. These correspond respectively to the directions e 1 (θ = 0 • ) and e 2 (θ = 90 • ) of Fig. 5.6. Therefore only 2 flow solutions will be necessary to fully determine the permeability tensor.

Numerical solution: Skeleton model

The set of points P w belonging to the yarn surfaces Γ w are extracted from the input geometry. The Medial Surface Skeleton is generated from the input points and is shown in Fig. 5.7. The skeleton mesh appears as a regular grid of thin surface channels, with a locally The medial skeleton mesh Ω is imported in the FEniCS library and Reynolds flow problem is solved following procedure detailed in Section 4.5.2, this time with boundary conditions from 5.20. In order to characterise the anisotropy of the material, the procedure is repeated twice: once for longitudinal flow (θ = 0 • ) and once for transverse flow (θ = 90 • ).

The solutions p, û are shown in Figs. 5.8 and 5.9 for longitudinal and transverse flow, respectively. Test case parameters and solution quantities are provided in Table 5.3.

Permeability determination

The velocity field computed from the pressure solution seems qualitatively consistent. In the longitudinal flow case, the resin travels faster in vertical slits between the arrays of yarns, and in particular at the junction between the slits and the interfaces between layers, where the gap h reaches its maximum value (see Fig. 5.7b). On the other hand in the transverse flow case, the resin shows maximum velocity in the narrowings between layers as expected (we remark that the values of velocity are always very low because ∆p = 1 Pa). From the velocity field û we can calculate the averaged velocity from Eq. 5.6 and then the permeability in the flow direction, for both longitudinal and transverse cases, using Eq. 5.9. The values are indicated in Table 5.3. The two values of permeability that we compute correspond to the principal permeability values: K 1 = 2.75E-09 m 2 , K 2 = 5.44E-10 m 2 . As expected, the textile exhibits a much lower permeability in the transverse direction.

θ Ψ [W] N e CPU t [s] u [mm/s] K [m 2 ] 0 • 7.5E-

RT/B: UD nested layers

This test case reworks elements from RT/A, with the difference that the layers of UD yarns are now stacked in a nested configuration rather than being exactly superposed. To achieve nested configuration, layers of textile are shifted pairwise by a length of T /2, T being the center-lines period, and vertically shifted. Geometrical parameters are provided in Table 5.4. In the nested configuration, the packing factor is higher (0.82) and the channels morphology and gap field h is different. The principal directions of permeability values K 1 and K 2 are known because they correspond to the direction longitudinal (θ = 0 • ) and transverse (θ = 90 • ) to the yarns respectively. Therefore we will follow exactly the same procedure as in test case RT/A.

Permeability determination

The average velocity in both longitudinal and transverse directions is calculated from the solution field and the corresponding permeability values are computed using Eq. 5.9. Values are reported in Table 5.5. Being the textile constructed along the principal domain axis, these values directly correspond to principal permeabilities: : K 2 = 4.36E-10 m 2 , K 2 = 3.83E-11 m 2 . We remark that, as should be expected, both values are lower than in the un-nested scenario RT/A.

RT/C: UD generic layup

In this test case we will use the Skeleton model to determine the permeability tensor of a generic layup, of which the principal direction of permeability are unknown. For the generation of the textile, we will be using the specifics from test case 4PS from chapter 3. The layup consists in 4 layers of C-Weave textile, which were stacked randomly. This is a "real" layup because the orientations were interpreted using Spectral Moiré analysis on the pressure print. Using the direction e 1 as reference, the layers were oriented following the sequence from bottom to top: are oriented at a rather small angle from the e 1 axis, with the exception of layer number 4 (θ 4 = 76.5 • ). This is also reflected in the principal direction of permeability K 1 , which is slanted by 8.5 • from the principal direction of test case RT/A.

Conclusions

The Medial Skeleton model has been successfully applied to generic textile layups through a series of test cases. The flexibility of the numerical tools developed for the generation of the Skeleton mesh allows this methodology to be independent of the complexity of the domain. The efficiency in terms of CPU time of the numerical solution allows this technique to be applicable to domains of considerable size, such as for the case RT/C.

N layers

Methodology: unsaturated flow

The permeability tensor of a material is usually treated as the input data to run filling simulations of industrial parts. However, by definition permeability is an averaged quantity. In practice, it means that the local variation of properties of the material and the effect of local features such as microstructural defects cannot be taken into account when a permeability tensor is assigned to the whole part. Furthermore, modeling the presence of local features inside RVEs in the estimation of permeability should be avoided, because the purpose of upscaling is to compute quantities which are representative of the entire part, not affected by the presence of local features. This is why local features are not addressed in the first half of the chapter. However, this is not true when dealing with unsaturated flow. Although their influence over the macro-scale permeability might be small, small local changes such as defects can affect and distort the flow pattern while the resin is advancing through the part. This behaviour can potentially lead to air entrapment and dry spots [START_REF] Patel | Modeling of void formation and removal in liquid composite molding. Part I: Wettability analysis[END_REF][START_REF] Bickerton | Investigation of draping and its effects on the mold filling process during manufacturing of a compound curved composite part[END_REF][START_REF] Siddig | A new methodology for race-tracking detection and criticality in resin transfer molding process using pressure sensors[END_REF]. In the second part of this chapter, we will attempt to use the Skeleton model to simulate numerically the unsaturated single-scale flow of a resin filling a preform. The finite elements/ control volume method (FE/CV) is the most widely used to compute filling simulations. The Skeleton model in its current formulation will be implemented into an existing FE/CV solver. This will allow to perform a direct filling simulation of a preform at the meso-scale, including the presence of potential local features, greatly expanding the potential of the model. Before moving to the actual implementation, we provide a brief overview of the FE/CV method and the solver that will be adopted. the operations: materials parameters assignment, boundary conditions definition, etc. However in this work, only the solver component of LIMS will be used, which can be accessed by scripting.

K 11 [m 2 ] K 12 [m 2 ] K 22 [m 2 ] K 1 [m 2 ] K 2 [m 2 ] p 1 p 2 θ princ 1.17E-
It is important to notice that LIMS is designed to solve the macro-scale Darcy flow, and not the Skeleton model. Therefore the implementation of the Skeleton model will need to be adapted to the capabilities of LIMS. where K e is the permeability tensor and µ is the viscosity. The permeability tensor K e is constant on the element, and LIMS allows one to assign an element-wise different permeability tensor to the mesh. By comparing the weak forms of the Darcy Eq. 5.28 solved by LIMS with Reynolds Eq. 4.41 solved by FEniCS, it is clear that the equations can be reduced to the same form by assigning to each element the permeability tensor:

Implementation of the Skeleton model

K e = h 2 e 12 e i ⊗ e i = h 2 e 12 0 0 h 2 e 12 (5.29)
By observing these modifications, the Skeleton model can be successfully implemented in LIMS. This implementation was validated: a Medial Skeleton model test case is set up and solved both using the established FEniCS approach from chapter 4 and the new LIMS implementation. The purpose is to prove that the same pressure field solution is obtained using both tools. This validation is detailed in Appendix B. Results substantiate the correct implementation of the Skeleton model in LIMS.

LIMS numerical procedure

The procedure to solve the transient flow simulation on a Skeleton mesh using LIMS is broken down in three steps:

-Skeleton mesh importation

4 [0 • , 0 • , 0 • , 0 • ] 20 20 2.6 0.73 1.0E+03

Workstation specifications

All test cases presented in this chapter were run on the same workstation detailed in chapter 4, Section 4.5.2. The LIMS solver was run on a Windows operating system installed on the same machine.

Test Cases: unsaturated flow

RT/AT

We propose again the same textile from the test case RT/A: four aligned layers of C-Weave textile, stacked in un-nested configuration. Geometrical parameters are reported in Table 5.10. The procedure detailed in Section 5.5.4 is followed.

Result discussion

Fill factor solution is presented in Fig. 5.17b, where the skeleton being filled is shown among the yarns, and in Fig. 5.18 as a series of time frames. The flow pattern is not straight, as the gap on the side of the yarns is larger than the interface gap between layers. The Skeleton simulation can capture this locality. After a time of t = 6s the flow front reaches the end of the mesh in the lateral gaps, but not yet in the interface gaps. A pressure difference of only 1.0E+03 Pa was imposed so that the filling would not be too fast, because the length of the domain is only 2 centimeters. Although being a meso-scale simulation, the efficiency of the Skeleton model allows to complete the simulation is less than one hour of CPU time. Values are provided in Table 5.9.

RT/BT

In this test case the same layup from test case RT/C is used. The objective is to demonstrate how the direct meso-scale simulation on the Skeleton mesh can capture the local anisotropy effect of a complex layup sequence without homogenising the permeability of the layers. A medial skeleton mesh is generated from the same layup as in test case RT/C. The layup is oriented in the θ = 0 • configuration as in Fig. 5.14a. Geometrical parameters are given in Table 5.11. This is an interesting case because in test case RT/C the traditional permeability tensor estimation was carried out following the channel flow procedure. However, when the material is highly heterogeneous due to the many orientations in the stack, the actual meaning of computing a macroscopic permeability tensor is questionable. The meso-scale simulation afforded by the Skeleton model bypasses completely this necessity, as the threedimensional connectivity of the complex channel network is directly taken into account.

Discussion of Results

Looking at the flow front progression in Figs. 5.19a and 5.19b neighbour layers. The advancing flow front develops an angle (α) due to the orientation of layers. This is in agreement to what observed in [START_REF] Bancora | Effective permeability averaging scheme to address in-plane anisotropy effects in multi-layered preforms[END_REF]. The angle is tracked (manually) in Fig. 5. 19.

A somewhat unexpected result, one can notice that the change in angle α is not monotonic during the advancement. This effect could not be captured by the macro-scale simulations of [START_REF] Bancora | Effective permeability averaging scheme to address in-plane anisotropy effects in multi-layered preforms[END_REF]. We suspect this effect is due to the alternating periods of crossover regions each layer experiences from neighbouring layers (i.e. the average flow front angle oscillates periodically). However, this behaviour might also be the result of an edge effect bias; further investigation would be necessary to identify the cause with certainty.

Finally we remark that the efficiency of the Skeleton solution allows to compute a mesoscale filling simulation over a domain orders of magnitude larger than the channels gap characteristic size. Values of number of elements and CPU time are reported in Table 5.9.

RT/CT: Racetrack defect

One critical factor that may affect the outcome of an RTM process is the presence of local defects in the part. Typical features that may be qualified as defects are race-tracks, uneven spacing between yarns, twists or distortions in the yarn paths and local changes in the compaction/fiber volume fraction. Local defects of different sorts can be represented by the Skeleton model, provided that their presence is accurately included in the textile geometrical model. In this section we will only look at one type of defect: the presence of a race-tracking channel in the layup.

Racetrack A well-known defect in LCM process is racetracking. This defect occurs when, somewhere in a layup, a gap is created either between yarns or between the part and the mould walls. This larger gap is a race-track, as it provides a preferential channel for the resin to flow. The consequence on the process is that the resin flow pattern will be distorted in an uncontrolled way by the presence of the defect. Depending on the magnitude of the defect, this may lead to quality issues.

The position of the racetrack and its orientation with respect to the filling direction is also a factor that affects the influence of the defect. The Skeleton model can be used to study the influence of racetracking. We setup a layup made by 4 layers of C-Weave textile, which are stacked following the sequence: [0 • , 45 • , 0 • , 45 • ]. The angles are measured relative to the injection direction. Using this layup, we recreate three scenarios: i) in the first scenario, no defects are present; ii) in the second scenario, a racetrack defect is introduced in the second layer (at 45 • from the flow direction); iii) in the third scenario, a racetrack defect is introduced in the third layer (at 0 • from the flow direction). For each case, the racetrack defect is created by widening one single inter-gap channel by a quantity gap rt = 1mm in the input geometry. Geometrical values about the domain are provided in Table 5.12. The Medial Skeleton mesh is generated and the filling simulation is computed for each scenario.

N layers

Orientations

L [mm] W [mm] H [mm] gap rt [mm] Packing factor p in [Pa] 4 [0 • , 45 • , 0 • , 45 • ] 45 30 2.7 1 0.73 2.0E+03

Discussion of results

Result values are reported in Table 5.9. Filling patterns from simulations are displayed in Fig. 5.21. Snapshots from the three different scenarios (absent, 45 • and 0 • racetrack) are taken at the same time steps. The snapshot comparison shows the effect of the racetrack on the flow front. We can see that the orientation of the racetrack plays a role in the alteration it causes. In the case b), when the racetrack is slanted at 45 • from the main flow direction, the effect is barely appreciable, and only results in some minor distortion of the flow front with respect to the base configuration. On the other hand in case c), when the racetrack runs in the same direction as the flow direction, a major distortion in the flow front is produced, with the flow front angle being disrupted and the resin reaching the end of the part sooner in the proximity of the racetrack. It should be noted that in scenarios b) and c) only one single racetrack channel is present in the entire layup composed by 4 layers. That means that the effect of the racetrack is somehow damped by the influence of the other 3 defect-free adjacent layers. Looking at the fill time values from Table 5.9, one can see that the racetrack has an effect on the time the flow front takes to reach the end of the part: in scenario c) the fill time is 20% lower than in scenario a). We can conclude that the Skeleton model could be used as a tool to directly integrate the presence of local meso-scale features in simulations of scale larger than the defect, in order to study their effect on surrounding regions of the part. It should be stressed that the 118 5.6. Test Cases: unsaturated flow accuracy of the prediction given by the Skeleton model is only as good as the assumptions of the skeleton model itself and the input geometry used to generate the skeleton mesh. At the present time, the reliability of the Skeleton model in simulating realistic scenarios has not been validated by any experimental comparison.

RT/DT: Radial flow

Although channel flow experiments were exclusively addressed so far, different injection conditions can easily be implemented in the Skeleton model. An example of radial flow injection is provided. For this study we use a layup of 4 layers in stacking sequence:

[0 • , 45 • , 0 • , 45 • ].
To be able to observe the radial flow front effect developing, the layup size is increased to a size of 100mm × 100mm. Layup and process parameters are provided in Table 5.13. The medial skeleton mesh is generated and the inlet pressure BC is prescribed to a small group of nodes in the center, marked as Γin (see Fig.

Conclusions

The Single-Scale Skeleton flow model has been successfully applied to multi-layer textile scenarios. Examples have been provided to show how the model can be used for both saturated flow (permeability estimation) and unsaturated flow (filling simulations to predict flow patterns). In particular, the computational efficiency that the Skeleton model offers over full-dimensional Stokes flow allows one to tackle meso-scale filling simulations over much larger domains than typical RVEs. Furthermore, the presence of local features such as racetracks can be directly included in the Skeleton flow simply by modeling their presence in the input geometry. The presence of these local features has been shown to have an effect on the flow patterns. Because of the prohibitive technical requirements associated, it was not possible at the current state to carry out a rigorous validation of the Skeleton model accuracy. Such a task would require to solve a Stokes problem in 3D on equally large domains, which is a challenging effort that comes within the competence of high performance computing. Such a validation should be carried out eventually. However, the series of validations performed in chapter 4, which focused on representative geometrical features, allow one to retain some confidence in the level of accuracy of the Skeleton model. This is supported by the fact that all simulations have shown qualitatively consistent results.

t [s] x f [m] y f [m] K 1 [m 2 ] K 2 [m 2 ]
The computational cost of solving Reynolds equation on the Skeleton mesh is reasonable already at the current level of implementation. However, it should be stressed that no effort was made on the optimisation of the computation efficiency. It is the author's belief that the employment of other solvers, especially compiled ones, as well as parallelisation, which was never used, might drastically improve the CPU time efficiency of the method.

This analysis concludes the study of the Single-Scale Skeleton model. However, in the context of composite materials, single-scale flow is only a partial solution. The presence of permeable tows, which were considered as solid in this chapter, can have an important effect on the resin flow behaviour. This is especially true in the case of unsaturated flow. To address such scenarios, in the next chapter we will resume and extend the Skeleton model to its dual-scale version.

Chapter 6

Dual-Scale Skeleton model

Introduction

Chapter 5 focused on the application of the Single-Scale Skeleton model to textiles to simulate unsaturated flow. In this chapter the Skeleton model is extended to its Dual-Scale version in which the assumption that the yarns are solid is relaxed. To model this, the porous geometry of the yarns is added to the channels network, as an additional domain available for the resin to flow and impregnate. To do so, a new type of hybrid skeleton mesh is introduced and a series of assumptions to solve the dual-scale equations are adopted. The model takes advantage of the efficiency of the medial skeleton to approximate the 3D flow problem at a lower computational cost. Similarly to what was previously done with the single-scale version, the dual-scale model is implemented in an existing commercial solver to compute filling simulations of textiles at the mesoscale.

This chapter represents the final step of the continuous workflow that started in chapter 3 with the acquisition of geometrical data from layups and continued into chapters 4, 5 with the single-scale skeleton model. The long-term objective is to create an efficient tool capable of handling dual-scale, unsaturated flow simulations at the scale of the yarns and over a domain size of a whole RTM part. While the current state of development is still far from this target, this chapter is to be interpreted as a first successful step towards this challenging objective. [START_REF] Facciotto | Modeling of anisotropic dual scale flow in RTM using the finite elements method[END_REF] simulation of the same scenario.

Outline of the chapter

Finally, the DSS model is implemented into the LIMS solver and the potential of the model is showcased through a series of representative test cases.

Dual-scale flow

The flow of resin impregnating a dual-scale porous fabric is affected by two characteristic lengths: the size of the pores in-between the yarns (the inter-yarn channels) and the size of the pores within the yarns (the gaps between individual fibers) as shown in Fig. 6.1. These characteristic lengths differ usually by one or more orders of magnitude. Consequently, the flow of resin in the medium is a dual-scale flow, which occurs at two scales: the mesoscopic flow between the yarns, and the microscopic flow within the yarns. The dual-scale flow in composite processing has been observed extensively. The effects on the resin flow front progression have been observed experimentally in [START_REF] Binetruy | Tow impregnation model and void formation mechanisms during RTM[END_REF], studied in [START_REF] Michaud | A review of non-saturated resin flow in liquid composite moulding processes[END_REF][START_REF] Bodaghi | On the variability of permeability induced by reinforcement distortions and dual scale flow in liquid composite moulding: A review[END_REF][START_REF] Binetruy | The influence of fiber wetting in resin transfer molding: scale effects[END_REF]. Dual-scale flow induced by capillarity effects has been observed in [START_REF] Verrey | Dynamic capillary effects in liquid composite moulding with non-crimp fabrics[END_REF][START_REF] Salvatori | Permeability and capillary effects in a channel-wise non-crimp fabric[END_REF].

There is a large body of literature devoted to the dual-scale flow in fabrics, both on modeling and experimental aspects. Although many models have been proposed, where some of which are presented here, this still continues to be an active area of research.

Most models agree on identifying three key factors that need to be handled: i) flow in the channels in between the fiber yarns (mesoscale region); ii) flow in the porous region within the fiber yarns (microscale region); iii) conditions at the interface between the two regions.

Channels region

The flow in the channels region is usually modeled as a creeping flow in a gap obeying Stokes equation. The behaviour is the same as discussed in chapter 4, with the difference that the no-slip condition is no-longer valid at the interface between the channel and porous regions. The region close to the geometric interface is actually a boundary layer marking the transition from the channel region to the porous region. In case of thin channel region, lubrication theory has been successfully applied even in presence of dual-scale flow [START_REF] Li | Derivation of modified Reynolds equation-a porous media model[END_REF].

Porous region

Although the flow of a fluid in a porous medium is usually described in terms of its phase averaged velocity by Darcy's law, this model fails to properly describe the behaviour of the fluid near the interface. This happens because the second-order viscous term present in Stokes equation has no corresponding terms in Darcy's law, leading to an inconsistency in the shear stress of the fluid at the interface [START_REF] Tan | Finite element implementation of stress-jump and stresscontinuity conditions at porous-medium, clear-fluid interface[END_REF]. To address this limitation, two approaches have been explored. The first consists in prescribing some slip-velocity boundary condition at the interface between the channel and the porous regions [START_REF] Beavers | Boundary conditions at a naturally permeable wall[END_REF]. This is known as Beavers-Joseph boundary condition and has been extensively studied in literature ( [START_REF] Saffman | On the boundary condition at the surface of a porous medium[END_REF][START_REF] Mikelic | On the interface boundary condition of Beavers, Joseph, and Saffman[END_REF][START_REF] Cao | Coupled Stokes-Darcy model with Beavers-Joseph interface boundary condition[END_REF]). The main drawback comes from the presence of a semi-empirical parameter (α BJ ), which depends on the material properties and needs to be determined experimentally or by some other means [START_REF] Valdés-Parada | Velocity and stress jump conditions between a porous medium and a fluid[END_REF].

The second approach uses the Brinkman equation [START_REF] Brinkman | A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles[END_REF], instead of Darcy, to model the fluid behaviour in the porous layer:

-∇p + µ ∇ 2 u -µK -1 • u = 0 (6.1)
where µ and µ are the Stokes and Darcy viscosities respectively, and K is the porous medium permeability (scalar when isotropic, otherwise second-order tensor). This is a well established model, that extends Darcy's equation with the addition of a second-order volume-averaged viscous term similar to the corresponding term in Stokes equation, that can handle the presence of the interface boundary layer. Far away from the interface, the shear stress in the fluid is negligible with respect to the viscous drag and Eq. 6.1 degenerates to Darcy's equation. It should be noted that the viscosity values in channel and porous domains are different terms, i.e. µ = µ . The Darcy viscosity µ depends on the characteristics of the porous medium, and their relation has been studied in [START_REF] Martys | Computer simulation study of the effective viscosity in Brinkman's equation[END_REF][START_REF] Kaviany | Principles of heat transfer in porous media[END_REF]. Sometimes µ = µ is adopted for simplicity. Coupling Stokes and Brinkman equations to describe dual-scale flow in fabrics made of fiber tows has been used extensively [START_REF] Syerko | A numerical approach to design dual-scale porosity composite reinforcements with enhanced permeability[END_REF][START_REF] Tan | Finite element implementation of stress-jump and stresscontinuity conditions at porous-medium, clear-fluid interface[END_REF].

Interface region When Brinkman equation is used together with Stokes equation to model the dual scale flow in a channel/porous domain, the stress continuity at the interface is achieved. However a stress-jump condition at the interface was developed by Ochoa-Tapia and Whitaker to better match experimental values when the porosity of the porous medium is high [START_REF] Ochoa-Tapia | Momentum transfer at the boundary between a porous medium and a homogeneous fluid-I. Theoretical development[END_REF][START_REF] Ochoa-Tapia | Momentum transfer at the boundary between a porous medium and a homogeneous fluid-II. Comparison with experiment[END_REF][START_REF] Valdés-Parada | Jump momentum boundary condition at a fluid-porous dividing surface: derivation of the closure problem[END_REF]. This is to be used with a modified version of Brinkman equation:

-φ∇p + µ∇ 2 u - µ φ K u = 0 (6.2)
which differs from Eq. 6.1 by µ = µ/φ, with φ being the porosity of the porous domain. Li [START_REF] Li | Derivation of modified Reynolds equation-a porous media model[END_REF] applied the stress-jump boundary condition in conjunction with the lubrication theory (Reynolds equation instead of Stokes) to address scenarios when the channel region is a thin layer. Although the methodology will be different, lubrication theory with stress jump at the interface will be applied to the model we propose in this chapter.

Macroscale numerical models

The scales at which Stokes-Brinkman equations are defined are traditionally not suitable for computation of flow over large domains, such as filling simulations in RTM processes. At the macro-scale, the geometrical features of the porous media are captured in the permeability term and interface boundary layer effects are not modeled. The macroscopic mass balance in non-deformable dual-scale porous media can be expressed as:

∇ • u = q (6.3)
where q is a sink term describing the resin entering the tows [START_REF] Pillai | Governing equations for unsaturated flow through woven fiber mats. Part 1. Isothermal flows[END_REF]. Eq. 6.3 states that some of the fluid flowing in the inter-yarn gaps (channel domain) escapes into the intra-yarn gaps (yarn domain). Simacek et al. [START_REF] Simacek | A numerical model to predict fiber tow saturation during liquid composite molding[END_REF] used arrays of 1D elements added to the nodes of the macroscale mesh to model the absorption of fluid by the fiber tows. In [START_REF] Imbert | Efficient dual-scale flow and thermo-chemo-rheological coupling simulation during on-line mixing resin transfer molding process[END_REF] and [START_REF] Babeau | Modeling of heat transfer and unsaturated flow in woven fiber reinforcements during direct injection-pultrusion process of thermoplastic composites[END_REF] the sink term in Eq. 6.3 is modeled as function of the microscale properties of the material.

Although in the present work all dual-scale modeling and numerical solution will be carried out at the meso-scale exploiting the skeleton mesh, the addition of "sink" elements remains at the core of the new DSS model that is introduced later.

Hybrid Skeleton Mesh

The Skeleton mesh, both in its 2D and 3D versions, has proven to be effective to model the channel network of a complex domain. Having been used to describe a single-scale flow, all the Medial Skeleton meshes shown so far have only represented the mesoscale network of channels running between the yarns. In the intent of developing a dual-scale model that uses the Skeleton mesh as a geometry descriptor, this definition of Skeleton mesh needs to be updated. We have defined in Section 6.2 the existence of three separate features that need modeling attention: the channels subdomain, the yarns subdomain and their interface.

We propose here a novel Hybrid Skeleton mesh that is based on the natural separation of the aforementioned elements.

Channels, yarns and their interface will be treated as three separate subdomains: each of them will be represented by an associated mesh. More precisely, the channels and the yarns subdomains will each be represented by a Medial Skeleton mesh (2D elements) that fully describes their meso-scale geometry. A third mesh, made of 1D elements, representing the interface between them will be used to spatially connect the channels and yarns subdomains. Each of these are now described and the procedure to obtain the Hybrid Skeleton mesh from a dual-scale domain is presented below.

Input mesh

The input mesh has no specific requirements in terms of geometry or element type. The only (trivial) condition that needs to be satisfied is that at least two separate subdomains can be identified in the mesh. These correspond to the channels subdomain and the yarns subdomain. We will refer to these subdomains in the input mesh as Ω 1 (channels) and Ω 2 (yarns) respectively. An example of input mesh is shown in Fig. 6.2a , with labels identifying these subdomains. Although the model was originally formulated in 2D, in the remainder of the chapter we will deal exclusively with 3D geometries, as it is the most complex scenario and has application potential.

Mesh generation

The Hybrid Skeleton mesh is composed of two distinct subdomains connected by an interface of one-dimensional elements. The physical subdomains Ω 1 (channels) and Ω 2 (yarns) are reduced to their skeleton versions Ω1 (channels), Ω2 (yarns) and linked through Ω3 (interface). An example of geometry with labeled domains is shown in Fig. 6.2.

Domains Ω1 and Ω2 :

The surface mesh Ω1 represents the channels domain. This is exactly the same type of skeleton mesh that has been used for all the single-scale test cases presented this far. Each element is equipped with the scalar value h e which describes the local gap height. The mesh Ω2 describes the yarns domain. This is also a Medial Surface skeleton mesh entirely analogous to the channels counterpart, Ω1 . Each element stores a scalar value h e that represents the local yarn thickness in the position occupied by that element. The nonconstant thickness of fiber tows of lenticular or elliptical cross-section is then captured. The skeleton meshes Ω1 and Ω2 are generated at the same time. The procedure is mostly similar to the one detailed in Section 4.4.3, the only difference being in the filtering step. Voronoi vertices are labeled as 'c' or 's' based on whether they fall respectively inside a channel or a yarn. After this step, Voronoi planes with all vertices flagged as 'c' are added to the C2/0 tables of the Ω1 mesh, while planes with all vertices flagged as 's' are added to the C2/0 tables of the Ω2 mesh. Delaunay triangulation follows for both subdomains, and the C'2/0 connectivity tables of triangular elements are created (see chapter 4, Section 4.4.3).

Interface domain Ω3 :

The generation of the one-dimensional elements that represent the channel/tow interface requires an additional procedure. After labeling all Voronoi vertices as either 'c' or 's', all Voronoi ridges containing at least one 'c' vertex and one 's' vertex are identified. This can be accomplished easily by filtering the C1/0 connectivity table of the entire Voronoi diagram. This filtering operation identifies a set of 1D lines that cross the yarn boundaries, with one vertex belonging to an element in Ω1 and one vertex belonging to some other element in Ω2 . Using the connectivity tables from the Voronoi diagram and some simple geometrical operations, each interface element is equipped with 2 scalar quantities that are essential for the dual-scale model:

r e : the radius of the associated region of influence. This is an approximation of the volume represented by the interface element, computed as the distance between the element center and the generating node of the associated Voronoi region. Using r e we can also compute the cross-sectional area associated to the element as A e = r 2 e π.

δ e : the fraction of the element length which is immersed inside the yarn . This quantity will be useful is section 6.4.

Domains and boundaries definition

The same notation scheme introduced for the single scale skeleton model is used for the dual-scale model. As already stated, Skeleton mesh subdomains are referred as Ω. Boundaries are labeled as Γ according to the same scheme used in chapters 4 and 5.

Dual-Scale Skeleton model

In the formulation of the dual-scale model, the stress continuity state guaranteed by the Brinkman equation is abandoned. The core idea is the following: the channel and yarn domains are modeled quite independently using their classical governing equations, respectively Stokes and Darcy. The interface elements are then used to ensure interaction between them. The adoption of interface elements to relate separate 2D meshes was inspired by the work in [START_REF] Chebil | Efficient numerical simulation method for three dimensional resin flow in laminated preform during liquid composite molding processes[END_REF], although in the present model, the equations and the numerical methods used are completely different. The starting point is the always valid continuity equation, which in its differential form and for constant density is:

∇ • u = f (6.4)
where f indicates a point-wise source term. Introducing the test function v we can transform Eq. 6.4 in its weak form to take advantage of the additive property of integrals:

Ω ∇ • u v dx = Ω f v dx (6.5
)

Ω 1 ∇ • u v dx + Ω 2 ∇ • u v dx = Ω 1 f v dx + Ω 2 f v dx (6.6)
The terms in Eq. 6.6 can be isolated and treated separately based on a separation of domains. We will deal with each one individually.

Channels domain

The terms in Eq. 6.6, valid in the channels domain Ω 1 , are expressed as:

Ω 1 ∇ • u v dx = Ω 1 f v dx (6.7)
Lubrication theory is adopted in the channels, similarly to what was done in [START_REF] Li | Derivation of modified Reynolds equation-a porous media model[END_REF].

We can then substitute the point-wise velocity field u with the gap-averaged Reynolds velocity u . This operation is possible by applying the volume averaging theorem of the divergence operator [START_REF] Pillai | Governing equations for unsaturated flow through woven fiber mats. Part 1. Isothermal flows[END_REF], which states:

∇ • u = ∇ • u + 1 V Γw n ⊥ • udx (6.8)
We remark the second term in Eq. 6.8 is null under the assumption of Poiseuille flow. Reynolds volume-averaged velocity is expressed as:

u = - 1 µ h 2 12 ∇p (6.9)
In the weak formulation of Eq. 6.7 we will consider the viscosity µ constant everywhere in the domain. The constitutive law from Eq. 6.9 is inserted into Eq. 6.7 to obtain:

1 µ Ω 1 ∇ • h 2 12 ∇p v dx = Ω 1 f v dx (6.10)
Eq. 6.10 is defined over the Medial Surface Skeleton mesh Ω1 . The value of the gaps h e is assumed as constant element-wise. This allows one to express Eq. 6.10 in FE discrete formulation, which after integration by parts yields: We can observe that the second term on the LHS of Eq. 6.11 has non-zero value only on domain boundaries where Neumann conditions are prescribed. To simplify notation we rewrite Eq. 6.11 as: Adding an interface flux term In order to account for the interaction between domains, a non-conventional source term is added to the local formulation Eq. 6.12 to allow some net outflow (or inflow) of mass in the elements. This term will be related to quantities modeled later on the interface domain Ω3 . Fig. 6.3 shows a typical triangular element from the mesh Ω1 . We can identify three subsets ∂ Ωi of the triangle surface Ωe , each associated to the corresponding local i node, such that: We represent the volumetric flow Q i entering the Ωe domain through a subset ∂ Ωi as:

i j k q i q j q k ∂ Ωi e ∂ Ωj
Q i = ∂ Ωi q i v dx (6.14)
where q i is the flux of resin through node i. The quantity q will have to be properly defined later. The total net flux in/out of the element is obtained by the sum of all contributions of fluxes associated with all the nodes of the element: This term can be seen as some extra surface flux entering the element through its surface. Adding this term to the expression in Eq. 6.12 we can express the weak form of the problem with the added interface flux contribution: At this moment, it is convenient to express Eq. 6.16 in matrix vector form. We adopt a choice for the test function as v = η i , and the trial function p = p i η i . Linear shape functions η i are used. The matrix-vector form of Eq. 6.16 becomes:

Q e = 3 k=1 Q k = 3 k=1 ∂ Ωk
Ne e=1 [[K e ][p e ] + [q]] = Ne e=1 [f e ] + Ne e=1 [g e ] (6.17)
where:

[K e ] = 1 µ Ωe h 3 e 12
∇η T ∇η dx (6.18) This assembly operation yields a global matrix vector system of equations in the form:

  K Ω1     p   +   q   =   b Ω1   (6.19)
where [b] represents the right-hand side of Eq. 6.17

Yarns domain

The terms in Eq. 6.6 that apply to the yarns domain are:

Ω 2 ∇ • u v dx = Ω 2 f v dx (6.20)
In this model, Darcy's law is adopted as the constitutive law for the resin velocity inside the yarns domain Ω2 . This choice is justifiable: since the presence of the interface is not explicitly modeled in the Ω2 domain, Brinkman equation simply reduces to Darcy's law.

By the volume averaging theorem on operators [START_REF] Pillai | Governing equations for unsaturated flow through woven fiber mats. Part 1. Isothermal flows[END_REF], we can then replace the pointwise velocity in Eq. 6.20 by the corresponding volume-averaged velocity u , plus a transverse flow term q (last term in Eq. 6.8) which will be treated separately as an interface flux term. Darcy's equation in the yarns reads:

u = - K µ ∇p (6.21)
The Ω2 skeleton mesh is composed by classic triangle elements which are identical to the elements in Ω1 . The weak formulation of Eq. 6.20 on Ω2 is achieved by following the exact same steps than in Section 6.4.1 and replacing the definition of Reynolds velocity from Eq. 6.9 with Darcy velocity from Eq. 6.21. The same interface flux term [q], accounting for the contribution of the last term of Eq. 6.8, is added even here to the weak form of Eq. 6.20. Therefore, the procedure will not be repeated here. Following the same manipulation we obtain matrix-vector weak form: where this time the local stiffness matrix of each element is defined as:

[K e ] = 1 µ Ωe h e k∇η T ∇η dx (6.23)
and k is the transversely isotropic permeability tensor assigned to the yarn. In the current implementation of the model, we assign a homogeneous, isotropic permeability at all points in the yarn, corresponding to the longitudinal permeability K of the yarn:

k = K e i ⊗ e i = K 0 0 K (6.24)
The transverse permeability of the yarn will be taken care of separately later, by the interface elements Ω3 . This assembly operation yields a global matrix-vector system of the same form as Eq. 6.19:

  K Ω2     p   +   q   =   b Ω2   (6.25)
where [b] represents the right-hand side of Eq. 6.20. Once again, as in domain Ω1 , we will consider that no volumetric source terms or Neumann boundary conditions are assigned, hence [b] = [0].

Interface domain

The skeleton subdomains Ω1 and Ω2 are interconnected by a dense array of 1D elements representing the interface in between. This set of elements is labeled as Ω3 . It is important to remark that these elements do not hold any kind of physical meaning (i.e. they cannot be related to some physical space), but merely represent a logical connection between interacting nodes among which an exchange of fluid is allowed. Therefore, special care and a set of reasonable assumptions must be observed when implementing this subspace in the model. Fig. 6.4 shows a typical interface element from Ω3 . If a flow of resin is allowed in i j q int q i q j A e (6.26) where u is the axial velocity of the flux, in direction of the rod, and A e is the area of the cross-section associated to the element, as defined in section 6.3.2. By adopting outwards flux from the element as positive, we can express the local conservation of mass in terms of balance of fluxes at the inlet and outlet nodes of the element: node (i) : q i + q int = 0 node (j) : q j -q int = 0 (6.27) To close the system, a definition of flux velocity u in the element is needed. The element does not represent a physical space where some constitutive law is well defined, therefore some modeling assumptions must be made.

q int = (u • n)A e = uA e

Resin velocity in the interface

The interface element carries the flow of resin between the mesh of the channel domain Ω1 and the mesh of the yarn domain Ω2 . Furthermore, the 1D interface elements are always oriented in the direction normal to the channel boundaries. The lubrication approximation is adopted in the channels, as a consequence the pressure gradient in the direction of the interface element is null in the region Ω 1 . On the other hand, since Darcy's law is adopted inside the yarns, the pressure gradient in the direction of the interface element is not null in the region Ω 2 . In the hybrid skeleton representation, the 1D interface element spans inside both regions, always sharing one node with Ω1 and one node with Ω2 . Therefore, a choice is made to model the flow velocity in the interface using Darcy, but only in the portion of the element inside Ω 2 , and considering u = 0 in the portion of the element inside Ω 1 (no transverse flow in the channel region).

Darcy's velocity (scalar form) in the 1D interface element is expressed as:

u = - K t µ ∂p ∂s (6.28)
where K t is the yarn transverse permeability and s denotes the curvilinear coordinate of the 1D element. To be consistent with the choice made, it is necessary to properly define the pressure derivative ∂p/∂s of Eq. 6.28 so that only the yarn region is taken into account. This can be achieved by defining the pressure derivative as a linear pressure gradient:

∂p ∂s = p i -p j L e δ e (6.29)
where L e is the total length of the element, and δ e is a scalar coefficient corresponding to the fraction of the element length which is immersed in the yarn domain Ω 1 in the original geometry, as defined in section 6.3.2. Adopting this expression, the linear system from Eq. 6.27 becomes:

q i -Kt µ p i -p j
Le δe A e = 0 q j + Kt µ p i -p j Le δe A e = 0 → q i q j = K e p i p j (6.30)

with

[K e ] = K t A e µ L e δ e 1 -1 -1 1 (6.31)
where A e represents the cross-sectional area associated to the element. A schematic representation of the pressure field as modeled on the interface elements is shown in Fig. 6.5. This system expresses the local net flux balance for a 1D interface element ∈ Ω3 between nodes i and j. A global system of equations can be built by assembling all Ω3 elements, to create a system of nDOFs equations:

  q   =   K Ω3     p   (6.32)

Global assembly

Local matrix-vector forms have been written in Eqs. 6.19, 6.24 and 6.32 for all hybrid skeleton subdomains Ω1 , Ω2 and Ω3 . In order to write a global assembly it is sufficient to equate the [q] term in these expressions. A global matrix-vector system is obtained:

  K Ω1     p   +   K Ω2     p   +   K Ω3     p   =   b   (6.33)
Finally by assembling all terms in a global form:

  K Ω     p   =   b   (6.34)
where

[K Ω] = [K Ω1 ] + [K Ω2 ] + [K Ω3 ] is a n × n global stiffness matrix, [p]
is the n × 1 nodal values vector and [b] is the n × 1 right-hand side vector of Eq. 6.5, n being the number of DOFs. The system of equations in Eq. 6.34 can be solved for the nodal coefficients [p] provided that some boundary conditions are prescribed. Without loss of generality, we remind that in all following cases we will adopt [b] = [0], and Dirichlet boundary conditions on pressure will be prescribed on some boundary Γ of the hybrid skeleton mesh.

Model test and validation

In order to validate the dual-scale skeleton model, a reference model is needed. As detailed in Section 6.2, the dual scale flow in a clear-fluid/porous medium domain can be conveniently modeled by Stokes-Brinkman equations. Therefore in the following, we solve the DSS model and compare it to the Stokes-Brinkman solution for a simple geometry. The following procedure was adopted: -the same dual-scale flow is solved using the DSS model;

-the velocity obtained from both solutions is compared

Test case geometry

A simple rectangular channel geometry is created using Gmsh. The channel is composed by two layers: one channel domain Ω 1 on top, and a porous domain Ω 2 on bottom. Both layers have the same constant thickness h = 0.2 m. The geometry is shown in Fig. 6.6a. Inlet and outlet boundaries are indicated. Geometrical and material parameters are provided in Table 6.1. We remark that the porous medium permeability value is quite high, and as a consequence the boundary layer thickness becomes important. These are not realistic values, however the choice was made in order to have a significant effect of the boundary layer on the solution. The weak form of Brinkman's equations, obtained by multiplying Eq. 6.2 by test function v and integrating by parts, can be written as: u plug = 2.00E-05 m/s on Ω 2 for s 1 = 0, s 2 = 0.15 (porous domain away from boundary layer).

One can see the non-negligible difference between u Ω 2 and u plug as a consequence of the high value of permeability K.

Numerical solution: DSS model

A numerical solver for the DSS model described in Section 6.4 is written in-house using Matlab scripting to solve Eq. 6.34 for solution of pressure p. For this test case, to be coherent with the SB model: µ = µ and K = K t = K are selected.

Channel domain

The DSS model yields a velocity value of u Ω1 = 6.75E-04 m/s on Ω1 .

From the SB model we obtain a higher value of u Ω 1 = 7.72E-04 m/s on Ω 1 . Overall, the difference is reasonable, with a percentage relative error (ξ) between the two values of ξ = 12.5%. The DSS model underestimates the fluid velocity in the channel region. We can attempt to provide an explanation: the DSS model is written adopting the no-slip condition, hence a Poiseuille velocity profile, in the channel. In the SB model the velocity at the tow interface is not null, since there is a transition between the Poiseuille flow in the channel region and the boundary layer in the porous region near the interface. This interface-slip velocity (see Beavers-Joseph boundary condition [START_REF] Saffman | On the boundary condition at the surface of a porous medium[END_REF][START_REF] Mikelic | On the interface boundary condition of Beavers, Joseph, and Saffman[END_REF][START_REF] Cao | Coupled Stokes-Darcy model with Beavers-Joseph interface boundary condition[END_REF]), which is absent in the DSS model, leads to higher velocity in the region close to the interface in the Stokes solution.

Porous domain

The DSS model yields a velocity value of u Ω2 = 2.02E-05 m/s on Ω2 , compared to the value of u Ω 2 = 2.90E-05 m/s on Ω 2 from SB model (ξ = 30.3%). However, when the strong velocity change in the interface boundary layer is neglected, the "plug" Darcy velocity inside Ω 2 becomes u plug = 2.00 E-05 m/s. This is a very close value to the one obtained from the DSS model (rel. perc. error ξ = 1%). The meaning is, no boundary layer is accounted for in the DSS model. Instead, the model limits the interaction between regions to pressure. This is clearly a simplification from the full-dimensional dual-scale Stokes-Brinkman model. On the other hand, the size of the boundary layer is small with respect to the overall domain, especially when the porous medium permeability is low [START_REF] Tan | Finite element implementation of stress-jump and stresscontinuity conditions at porous-medium, clear-fluid interface[END_REF]. Away from the interface, when the Darcy velocity is stabilized, the DSS model agrees well with the SB model.

Conclusion

The CPU time needed to compute solution according to the models is different. Values are reported in Table 6.3. One must remark that this comparison is arbitrary, since the number of elements is very different (the analysis was carried out with validation purpose, so a high number of elements were used in the SB solution). Even so, the DSS model shows a large advantage in CPU time.

Based on all the considerations drawn, we conclude that the DSS model is a viable first attempt at modeling the dual-scale flow taking advantage of the skeletonized version of a domain. The main drawback is the difference in the channel region velocity, which is underestimated by ≈ 12% in the test, due to the lack of slip-velocity at the interface. This condition could be addressed in future refinement of the DSS model. However, at the current level of development, and in reason of the low CPU time trade-off, this approximation is deemed acceptable especially if the goal is to predict dual-scale flow in large complex geometries.

It is worth mentioning that the simplification of the Stokes-Brinkman model by the Reynolds-Darcy model is often used in the context of the homogenization of the permeability of porous media with bidisperse porosity [START_REF] Tran | Unified microporomechanical approach for mechanical behavior and permeability of misaligned unidirectional fiber reinforcement[END_REF]. For example, the interactions between the inter-tow and intra-tow regions, obtained here with 1D interface elements, are reflected by the Eshelby tensor in the context of the mean field homogenization. It has been shown that this simplified model nevertheless provides satisfactory results in the context of composites reinforced with fabrics made of fiber tows.

In the remainder of the chapter, efforts are made to implement the DSS model in a more powerful solver than the in-house Matlab solution, to take advantage of its full potential in complex scenarios.

One final remark is due: the single test case that was carried out for validation should not be considered exhaustive for the quantification of the DSS model approximation. Both the geometry and the boundary conditions used resulted in a somewhat "simple" test case. The comparison between DSS and SB models should be carried out over more complex test cases in order to properly assess the potential and limitation of the model we propose. However, in this first implementation, the objective is shifted on building the complete workflow that was stated in the introduction of the chapter. Therefore the simple validation presented will be deemed sufficient for the task at hand.

LIMS implementation

The final objective of this chapter is to use the dual-scale skeleton model to simulate the unsaturated flow of resin during the filling of a textile preform. To do so, we resort to the same numerical tool that was used in chapter 5: LIMS. This solver was developed to model Darcy flow only, and is not capable of solving a dual-scale problem on a mixed elements skeleton mesh. However, this is not entirely true: a model mixing 2D and 1D elements to simulate the dual-scale impregnation of tows was already implemented in LIMS by its developers [START_REF] Simacek | A numerical model to predict fiber tow saturation during liquid composite molding[END_REF]. Because of this work, LIMS already possesses all the core functionalities needed for the implementation of the dual-scale skeleton model. Since the equations differ, some modification is required. By acting on the material and process variables assigned via the *.dmp input file, the same form of equations as in section 6.4, hence the same pressure solution, can be obtained. However, because of the redundancy of elements in the hybrid skeleton mesh (1D interface elements), a special care will be needed to correctly compute the filling. The next section outlines how the input data can be adapted so that the equations solved by LIMS and the flow advancement correctly reflect the DSS model.

LIMS solution scheme

The solution scheme of the LIMS solver is detailed in chapter 5. The two steps solved at each time increment (pressure solution, fill factor update) are isolated and the modifications needed to replicate the DSS model are highlighted. The baseline is the mixed 2D-1D elements solution scheme already implemented in LIMS in [START_REF] Simacek | A numerical model to predict fiber tow saturation during liquid composite molding[END_REF].

Pressure solution

LIMS solves the classical Darcy-Continuity equations at every time step:

u = - K µ ∇p (6.37) ∇ • u = 0 (6.38)
The objective is to obtain the same solution as the one yielded by the skeleton model. To do so, it is convenient to express Eq. 6.37 in matrix-vector weak form:

[K][p] = [f ] (6.39)
When the mixed 2D-1D elements scheme is used, the global stiffness matrix can be split in the 2D ( ) and 1D ( ) contributions:

[K ][p] + [K ][p] = [f ] → Ne e ([K ] e [p] e + [K ] e [p] e = [f ] e ) (6.40)
Each term on the left-hand side of Eq. 6.40 is taken care of separately.

2D elements

The 2D elements contribution of the global stiffness matrix is represented by summation of all local terms, each expressed as:

[K ] e [p] e (6.41)
In the case of triangular linear elements as used by LIMS, the local stiffness matrix is obtained according to classical FE theory as:

[K ] e = 1 µ
Ωe h e k e ∇η T ∇η dx (6.42)

In the dual-scale skeleton model, the 2D elements are found in the channels domain Ω1 and in the yarns domain Ω2 . The local stiffness matrix forms of the corresponding elements are given in Eqs. 6.18 and 6.23. By comparing them with Eq. 6.42, one can see that the form is identical and the same quantity can be expressed by simply acting on the definition of the element permeability tensor k e . The implementation of the local 2D elements contribution from the dual-scale skeleton model into the LIMS solver is then achieved by simply defining:

k e = h 2
k e = K e i ⊗ e i = K 0 0 K on Ω2 (6.44) (6.45)
where K is the longitudinal permeability of the yarns.

1D elements

The 1D elements contribution of the global stiffness matrix is the summation of all local element stiffness matrices:

[K ] e [p] e (6.46)

LIMS uses rod elements to compute the local stiffness matrix of its 1D elements. According to classical FE formulation, this is expressed as:

[K ] e = 1 µ
Ωe A e k e ∇η T ∇η dx (6.47)

Linear shape functions (η) are adopted by LIMS. If we denote by x the arc length coordinate within a generic line element delimited by nodes i and j, and by L e the length of the element, that yields:

η i η j = -1 Le x + 1 2 1 Le x + 1 2 (6.48) ∇η i ∇η j = 1 L e -1 1 (6.49)
Developing Eq. 6.47 we obtain the form:

[K ] e = k e A e µ L e 1 -1 -1 1 (6.50)
This is the form of the local stiffness matrix computed by LIMS for each 1D element. Admitted input parameters element-wise are: k e , A e The fact that linear shape functions are used is very convenient, because it can be adapted to the 1D interface elements of the skeleton model. Comparing Eq. 6.50 with the 1D element of the skeleton model from Eq. 6.31, one can see that the same quantity can be computed by simply setting:

A e = A e δ e (6.51)
and assigning the element-wise value A e to the 1D interface elements in the skeleton mesh input file. We also assign the yarn transverse permeability (K t ) as the element permeability:

k e = K t (6.52)

Filling solution

The filling step is straightforward: once the pressure field is correctly computed at all nodes with the modifications presented above, the velocity at the resin front is computed to prime the filling of the elements. LIMS computes Darcy's phase averaged velocity, which is defined element-wise as:

u e = - k e µ ∇p (6.53)
2D elements Since the definitions of k e detailed in Eq. 6.43 are being assigned, the proper velocity is already computed on the 2D elements ( Ω1 and Ω2 ). After computing the velocity, the volumetric fluxes going into the elements are computed using Eq. 5.27

q =

Γe u e • ndx (6.54)

and the fill factor in the elements is updated. If the fill factor value in an element is denoted as σ e , the fill factor at time step t is calculated as:

σ e(t)
= σ e(t-1) + q A e h e φ e ∆t (6.55)

The filling of the 2D elements is computed correctly by assigning:

φ e = 1 in Ω1 (channel region which is empty space)

-

φ e = (1 -V tow f
) in Ω2 (space available is a function of yarn's porosity)

The proper porosity value is assigned element-wise to 2D elements in Ω1 , Ω2 using the mesh *.dmp input file.

1D elements

The 1D elements need special care. In fact, the volume of the entire domain, i.e. empty channel and porous region, is completely described by the corresponding Ω1 and Ω2 skeleton meshes. The volume represented by the 1D interface elements is an extra volume that wrongly adds to the total. The interface elements serve the important purpose of connecting the DOFs of the Ω1 and Ω2 and ensuring the pressure solution is properly achieved. After that, they become redundant and, even more so, hinder the proper execution of the filling step.

One remedy adopted to address the issue is to define an area correction factor, a large scalar value that shall be named λ. We shall adopt, for instance, λ = 1E + 05. Eq. 6.50 is modified to be multiplied and divided by λ. In the specific, integrating also the modified version of Eq. 6.51:

[K ] e = k e A e µ L e λ λ 1 -1 -1 1 (6.56)
We proceed to define the following quantities:

k e = k e λ (6.57
)

A e = A e λ = r 2
e π δ e λ (6.58)

Both quantities k e and A e can be assigned element-wise according to their definitions from Eqs. 6.57 to the 1D interface elements using the mesh *.dmp input file. By adopting this approach results in: 6.7. Test cases -the pressure solution that was defined is not altered, since the value of Eq. 6.31 remains the same. The contribution of the 1D interface elements to the pressure field has not been modified.

-the velocity computed using Eq. 6.53 in the 1D interface elements is much higher, by a factor λ, than it should be.

-the volume associated with the 1D interface elements is much smaller, by a factor λ. This condition together with the previous one results in the 1D interface elements being filled almost instantaneously.

-the volumes Ω1 and Ω2 , associated to channels and yarns, are filled normally.

-The extra volume associated with the 1D interface elements that becomes filled (quickly) is very small, and the error committed in the total volume to be filled is of order

O = λ -1
The final result which is obtained by adopting this modification is that the correct pressure solution is conserved step-wise, at the cost of a very small approximation in the filling step.

Validation of LIMS Dual-Scale Skeleton solution

The validity of the entire filling simulation cannot be directly compared with any other reference solution. Such validation would require a full-3D Stokes-Brinkman FE/CV filling simulation software. These solutions do exist (for example, GeoDict by Math2Market), but their accessibility and implementation are out of scope for this thesis. It is undeniable however that a rigorous validation should be eventually performed, because at the current time it is impossible to estimate the degree of approximation of the DSS model.

One validation that can be carried out, instead, is the correct implementation of the DSS model in the LIMS solver. As such, we will proceed in identical fashion as what was done for the implementation of the single-scale skeleton model in chapter 5, i.e. ensuring the pressure solution is correctly maintained.

The validation is detailed in Appendix C. The results from the analysis show that the difference between the LIMS and the custom-written reference solution is only 0.005 %. The result substantiates the correct implementation of the DSS model within LIMS.

Test cases

The remainder of the chapter presents a few representative test cases in which the dualscale filling of a multi-layer layup is simulated at the mesoscale using the DSS model. The procedure followed for all cases is the same:

-an input textile geometry is generated using either Gmsh or TexGen; 6.7. Test cases -the hybrid skeleton mesh is generated; -the mesh is imported in LIMS and the DSS model is implemented by following the procedure from Section 6.6;

-the simulation is run, and the dual-scale filling is visualised;

For all presented cases, a resin viscosity of µ = 0.3 Pa•s and a yarn fiber volume fraction of V yarn f = 0.5 are used. Tow permeability values K and K t are provided case by case.

Boundary conditions

The same boundary conditions from Eq. 5.20 of inlet pressure (p = p in on Γin ) and outlet pressure (p = p out on Γout ) are imposed for all cases. For all cases, the values of inlet and outlet pressure are set to p in = 1.0E + 03 Pa and p out = 0 Pa, respectively. The only exception is represented by case DS/C, where a value of p in = 1.0E + 05 Pa is used. The proper pressure difference is chosen to avoid a too fast filling of the mesh, since the domain lengths are in the order of centimeters.

Workstation specifications

All test cases are run on the same workstation detailed in chapter 4, Section 4.5.2, running LIMS solver on Windows operative system.

DS/A: Single yarn in channel

The purpose of the first case is to focus on the dual-scale filling mechanism using an elementary geometry. The DSS model is applied to the simplest scenario: one yarn is immersed in an empty channel running above and below. The hybrid skeleton mesh used for the simulation is the same that was used in Section 6.3 to introduce the hybrid skeleton mesh, re-proposed here in Fig. 6.10a. Dimensions are representative of situations encountered in fabrics, except for the reduced length L: the geometry is created for illustrative purpose. Geometrical and material parameters are reported in Table 6.5. Boundary conditions are applied, the DSS filling is computed by LIMS and results are exported.

Discussion of results

The filling pattern is shown as a time sequence in Fig. 6.11. The interface 1D elements are not displayed due to format constraints, which is fine because their presence would clutter the visual. One can notice the faster flow quickly filling the channels Ω1 first, and the slower flow front following in the yarn Ω2 . We remark that the entire filling is rather quick because the domain is short and p in = 1.0E+03 Pa. The role of the interface elements is evidenced by the partially-saturated region that precedes the saturated flow front in the yarn. This region is the effect of the transverse flow coming into the yarn from the already filled channel regions above and below. One can also notice the concave shape of the flow front in the yarn, which is due to the higher thickness of the yarn in the middle (more volume to fill). Finally we remark that the CPU cost needed to compute the dual-scale filling is in the order of seconds (Table 6.4).

In conclusion, the DSS model provides qualitatively satisfactory results, although it is impossible at the current stage of development to express any quantitative comparison.

DS/B: Array of yarns with variations

The purpose of this test case is to study how the dual-scale filling of an array of yarns is affected by local features: in particular, the yarn permeability and the presence of racetracks. An array of four yarns are lined up in a channel defined by two parallel plates, as shown in Fig. 6.12a. The geometry is modeled using Gmsh. A flow is prescribed along the yarns. Three different scenarios are considered:

-scenario a): the yarns are assigned the same permeability value as in case DS/A. -scenario b): the yarns permeability is decreased by one order of magnitude.

-scenario c): a racetrack gap of size d is introduced in the middle of the yarns.

For each of these scenarios, a different input mesh is created and the corresponding hybrid skeleton mesh is extracted. Geometry and permeability values for the three scenarios are provided in Table 6.6. The hybrid skeleton mesh from scenario a) is shown in Fig. 6.12b as an example.

Discussion of results

The results from the three scenarios are displayed in Fig. 6.14. In the interest of the study, the top of the mesh is hidden to show both channels and yarns underneath at the same time. This visualisation cut is illustrated in Fig. 6.13. All three scenarios show the dualscale flow behaviour: the flow front proceeds faster in-between the yarns, and follows after some lag distance in the yarns. The flow front in the yarns is preceded by a partiallysaturated region, similar to what was observed in case DS/A, caused by the transverse flow coming from the filled channels. The variation of isolated parameters has an observable effect on the flow. In Scenario b), the yarn permeability, both longitudinal and transverse, is one order of magnitude lower than in Scenario a). As a consequence, the resin advances more slowly in the yarns and the partially-saturated region is longer. With the yarns being less permeable, the flow rate of the resin flowing transversely from the channels to the yarns is lower and hence the channels are filled faster. In Scenario c) the larger racetrack between the two central yarns has a trailing effect on the entire dual-scale flow, distorting the flow front pattern. As the central channel is filled first, the transverse impregnation of the yarns originates in prevalence from the center. The partially-saturated region in the yarns is comparable in length to the one of Scenario a), which is justified by the same value of yarn permeability. All simulations require a higher CPU time to carry out than case DS/A in virtue of the much higher number of elements. Values are reported in Table 6.4. The CPU time across all Scenarios is fairly similar, given that the number of elements is comparable, with the notable exception of Scenario c). This might be due to the fact that, since the front is distorted, by the time the lateral unsaturated regions are filled the resin velocity is very low and many more time steps are necessary. This is only an interpretation since the flow front position with CPU time was not tracked, but this investigation goes beyond the scope of the study. Considering the size of the domain (20 × 15 mm) and the hardware used, the values of CPU time remain affordable. The hybrid skeleton mesh is extracted and shown in Fig. 6.16. The two main components of the hybrid skeleton, the channels mesh Ω1 and the yarns mesh Ω2 , have been isolated in the figure. It is possible to see the presence of the racetrack channel in the middle of the Ω1 mesh. In order to contain the memory usage during the filling simulation, a second, coarser mesh is generated from the one in Fig. 6.16 by merging smaller elements together. The h field is conserved by interpolation of the original field over the coarser mesh. This allows one to decrease noticeably the number of elements and the CPU cost. Inlet and outlet regions are defined as in a channel flow experiment ( Γin and Γout in Fig. 6.16). Unlike the previous test cases, the value of inlet pressure has been increased to p in = 1.0E+05 Pa because of the larger domain.

Discussion of results

The computed flow patterns in layers are displayed as a time sequence in Fig. 6.17. The mesh is cut for visualisation at every layer, so that the flow progress can be tracked individually. At each level, it is possible to see the yarns being progressively filled with the flow front preceded by the usual partially-saturated region as a result from the transverse infiltration. As expected the flow front in the channels is more advanced. We remark that the filling is carried out in a very short time since the imposed pressure in this case was increased to p in = 1.0E+05 Pa. The effect of the racetrack on the flow front can be observed in layer 2. Its influence can be observed also in the adjacent layers 1 and 3, less in layer 4 (the most distant). In particular, in Fig. 6.17 we have highlighted a yarn region in layer 3 which is located just above the racetrack. An unfilled area in the yarn is entrapped between the advancing flow front saturating the yarn and the transverse flow infiltrating from the racetrack. A scenario like this represents a potential risk of air entrapment and dry spot 6.8. Conclusions formation.

Regarding the computational cost of the simulation, CPU time with respect to number of elements is reported in Table 6.4. The time needed, slightly higher than one hour, is deemed reasonable considering that a dual-scale filling simulation is being carried out on a 35 × 35 mm domain composed by 4 layers. We remark that the CPU time needed was comparable to test case DS/B, even though the number of elements is lower. This result was checked by running the simulation 3 times. The reason is unknown, although it is most certainly related to the time stepping and solution scheme of LIMS.

Conclusions

In this chapter we have proposed a new Dual-Scale Skeleton (DSS) model which can handle the dual-scale flow onto a novel hybrid skeleton mesh, both in saturated and transient form. The model was developed as an extension to the Single-Scale Skeleton model that was introduced in chapters 4 and 5. The accuracy of the DSS model was assessed by comparison against a full-dimensional Stokes-Brinkman solution. However, in this regard, an additional validation over more complex scenarios is still necessary.

The main objective of the chapter, which brings the entire thesis to a fitting conclusion, was the application of the DSS model to the simulation of the dual-scale filling of composite parts. This was achieved by implementing the DSS into the software LIMS. Using this new tool, a series of representative test cases were presented in which the filling of arbitrary dualscale porous domains was carried out on their hybrid skeleton representation. At the current level of development, the results require further validation, which could be achieved only by either experimental measurement or by solving the same scenarios numerically using a high-computing Stokes-Brinkman transient flow solver. It is worth to note that experimental validation is still very challenging. The dual scale flow in fabrics is very difficult to monitor and is subjected to experimental errors. However, all the results show good qualitative agreement with how the flow is expected to behave in presence of their distinctive features (racetracks, low/high transverse yarn permeability...).

The approximations introduced by the DSS model are well balanced by the computational efficiency: filling simulations for all test cases were carried out in a CPU time of about one to two hours, a reasonable result considering that: i) no solving optimisation was made; ii) a modest workstation was used; iii) domain sizes were larger than typical RVEs. The actual efficiency of the model should be further assessed by obtaining a term of comparison from a full-dimensional Stokes-Brinkman transient filling solution.

In summary, the methodology presented in this chapter can be considered a first successful attempt at applying the novel Dual-Scale Skeleton model to the direct meso-scale simulation of RTM processes.

Chapter 7

Conclusions

Summary and contributions

Driven by the aim to bypass the need for traditional permeability characterization while capturing the essential geometric features of fibrous preform, we have presented a contribution to the modeling of RTM processes. The first part of the thesis has focused on the acquisition of relevant geometrical data from fibrous textiles. We have proposed a novel measurement methodology, based on an affordable pressure-sensitive film, to interpret fine mesoscale features of a textile in its compacted state. In particular, we have shown how the theory of Moiré interference, combined with image spectral analysis, is able to detect the precise orientation and spatial frequency of yarns of individual layers in a multi-layer stack. This represents a novel application of the theory, introducing a new non-destructive tool of geometric investigation to go alongside with or replace the traditional ones (µ-CT [START_REF] Barburski | Micro-CT analysis of internal structure of sheared textile composite reinforcement[END_REF], photography [START_REF] Gommer | Quantification of mesoscale variability and geometrical reconstruction of a textile[END_REF]) . This technique has shown promising potential: the necessary data is the compaction pressure field, and this measurement could be integrated in current manufacturing workflow, since the preform compaction is already an integral step in the production process. Furthermore unlike µ-CT, using a cost-effective pressure film removes any limitations to the size of measurement.

In the second part of the thesis, the geometrical data obtained was used as input to develop an efficient flow model based on skeletonization. Adopting the lubrication approximation, the creeping flow typical of RTM processes was solved on a Skeleton mesh representing the meso-scale architecture of the material at a reduced computational cost. Although the use of lubrication theory in textiles is not novel, for the first time the Reynolds flow model was solved on complex three dimensional skeleton geometries that represent the inter-yarn channels of a preform. This method allows one to compute the upscaled permeability of the domain, but more importantly allows one to directly solve the mold filling flow problem on an entire preform. The stated objective is to achieve the direct simulation of an RTM flow, at the mesoscale, on a domain representing the whole part. While the current level of implementation is still far from the goal, we were able to solve domains on the order of several centimeters of length at a reasonable computational cost and without any numerical optimization, substantiating the potential of the methodology in modeling local flow effects.

Finally, in the last chapter of the thesis, we have extended the model to propose the first development of a Dual-Scale Skeleton (DSS) model. Based on a specifically designed hybrid skeleton mesh, the model is capable of solving the dual-scale transient flow impregnation of a textile preform taking into account the filling of both inter-yarn and intra-yarn pore spaces. The DSS model was implemented in the existing solver LIMS to tackle a series of representative cases, and the results have shown good agreement with known flow behaviours.

Unlike with macroscale models, the Single and Dual-Scale Skeleton models are able to directly account for local defects or variability in the material, relaxing any assumption of periodicity, which makes them valuable for a future application in the simulation of RTM processes.

Future prospects

Along with the answers provided, this research opens many new questions. Potential future developments that spawn from this work are numerous, and can be conceptually divided in two groups: geometrical data acquisition and skeleton flow simulation.

Geometrical data acquisition

The first part of the work was dedicated to using a pressuresensitive film to measure and interpret geometrical data from textiles. At the current level of development, only flat geometries were investigated using this sensor.

-While the doubly-curved surfaces are a different challenge in terms of implementation (a deformable sensor is needed), it would be interesting to extend the analysis to single-curvature shapes, which are commonly encountered in RTM parts.

The Prescale pressure system has served well its purpose to lay the basis for a new methodology: however, this pressure film was not developed for the application we envisioned, and presents some strong limitations in resolution, data readability and fluid resistance.

-In light of an industrial application, it would be wise to resort to a different system for pressure field measurement: possibly a digital solution (section 2.2.1) directly integrated in a "smart mould", so that the pressure data from closing a preform in the mould could be directly read and analysed, and the simulation of flow could be carried out, before the resin is injected into the actual part.

In terms of the methodology itself, the spectral Moiré analysis of pressure paves the way for future studies. In chapter 3 we have detailed the complete theory on how to relate the position of impulses in the FFT to the orientation and spacing of tows. The focus of the analysis was put on the detection of orientations only. 7.2. Future prospects -It would be worthwhile to study also the variation of inter-yarn spacings using the same technique. For instance, we can hypothesize that one layer in a stack with yarns closer together than nominal (for instance, because it has been sheared) would show up on the FFT as an impulse with higher frequency than expected (outside of the orientation circle), indicating an anomaly.

-Another development is possible in the textiles architecture: although we focused on quasi-UD textiles (gratings), we have remarked that Moiré theory can also be applied to periodic grids [19] and non-periodic structures [START_REF] Amidror | The Theory of the Moiré Phenomenon[END_REF]. Based on these premises, one could explore if spectral Moiré analysis based on pressure would be able to characterize woven fabrics (grids) or more complex features, such as curved yarn paths in automated fiber placement technologies.

-Finally, one interesting extension to what discussed in chapter 3 would be to couple spectral Moiré analysis with a mechanical model, able to use the pressure intensity information to extract some information about the yarn cross-sections. This would add to the realism of the geometric reconstructions (see section 3.6) made possible by the analysis of the compaction pressure field.

Skeleton flow simulation

The Skeleton model for flow, both in its single-scale and dualscale (DSS) versions, has been applied to simulate the infusion of representative test cases at mesoscale level. Three-dimensional flows are solved on two-dimensional skeleton meshes, with good potential in terms of computational savings when compared to full 3D solutions. The model trades off some accuracy to the advantage of the size of domain that could be solved. However, at the current level of implementation, we were only able to scale up to domain sizes of a few centimeters (case RT/DT in chapter 5 is 10×10 cm) which, although larger than traditional RVEs, are still far from the target of a whole part. The cause of this limit is not in the model itself, but rather in the preliminary generation of the skeleton mesh. Using the Voronoi diagram method, the memory requirements for the generation of the skeleton mesh increase steeply with the domain size L (O(L 2 )). Furthermore, the large scale separation between the inter-yarn channel dimensions and the domain size implies that a very high density of generating points must be used in order to properly compute the channel gap values. This results in a very fine skeleton mesh, independently of the domain size.

-With the objective to generate large skeleton meshes at the scale of entire parts, the first step of improvement is to abandon the Voronoi diagram method for the definition of the skeleton nodes. Alternative methods such as the "grassfire" algorithm were discussed in section 4.4.1: it would be worth to test if they could yield better performance with large geometries. However it is the author's belief that the best approach would be to take advantage of the parametric description of a textile (yarn paths, cross sections) to directly generate a "parametric skeleton mesh" without using a set of input points.

Finally, some future prospect should be dedicated to the Dual-Scale Skeleton (DSS) model. Applying the model to simulate flow in classical scenarios we have obtained results which are qualitatively satisfactory. However, a rigorous validation by comparison with a full three-dimensional Stokes-Brinkman unsaturated flow is still lacking.

-A validation by comparison with non-skeleton solution should be carried out. Depending on the domain size, such validation could have a high computational cost but would be useful to quantify the approximation introduced by the DSS model.

-The DSS model is already capable of treating yarns as transversely isotropic by assigning different values of longitudinal and transverse permeability to the yarn elements ( Ω2 ) and interface elements ( Ω3 ) respectively. This scenario was not tackled, but it could be interesting to explore the effect of varying transverse permeability on the flow.

-Capillary effects inside fiber tows are not modeled in the DSS model and hence it would be worthwhile to include this behavior by identifying regions where the capillary flow is not negligible.

-So far, fibrous structures were considered as non-deformable porous media. Flowinduced decompaction in liquid resin infusion could be developed, expanding the applicability of the model to a wider range of LCM processes [START_REF] Hemmer | Unloading during the infusion process: direct measurement of the dual-scale fibrous microstructure evolution with X-ray computed tomography[END_REF]. The DSS model is suitable because it includes the h metric of local thickness, which should be modified to evolve in time with the flow progression.

work however we do not rely on geometric Moiré, which can be used as an alternative tool of interpretation. where K is the permeability tensor of the material, v is the test function and µ is the viscosity. The permeability tensor K is constant on the element, and LIMS allows to assign an element-wise different permeability tensor to the mesh. By comparing the weak forms of the Darcy Eq. B.1 solved by LIMS with Reynolds Eq. 4.41 solved by FEniCS, it is clear that the equations can be reduced to the same form by assigning to each element the permeability tensor: The tensor k e is computed for each element in the Skeleton mesh and a *.dmp file is written to import the mesh into LIMS. Viscosity from Table B.1 and and boundary conditions from Eq. 5.20 are applied. The transient filling simulation is carried out and a a snapshot file in *.tp format is saved at every time interval ∆t = 10 s. It should be noted that the only time when the pressure solutions from FEniCS (saturated flow) and LIMS (transient flow) are comparable is when the mesh is completely filled, i.e. at the last computed time step. The pressure field form the last time step is shown in

1 Introduction 1 .

 11 sponding hybrid skeleton mesh . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Triangle element with interface flux . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Interface element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 1D element Ω3 connecting the domains Ω1 , Ω2 . The assumed pressure gradient is represented . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 The geometry used for DSS model test and validation . . . . . . . . . . . . . . 6.7 DSS model validation: Stokes-Brinkman solution (velocity) . . . . . . . . . . xiii 6.8 DSS model validation: extracted SB velocity plot along s . . . . . . . . . . . . 6.9 DSS model validation: hybrid skeleton solution (velocity) . . . . . . . . . . . 6.10 DS/A: Input geometry and hybrid skeleton mesh . . . . . . . . . . . . . . . . 6.11 DS/A: result fill sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.12 DS/B: Input geometry and hybrid skeleton mesh . . . . . . . . . . . . . . . . 6.13 DS/B: close-up view of the solution features and visualisation cut of the interlayers channel domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.14 DS/B: result fill sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.15 DS/C: Input geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.16 DS/C: Hybrid skeleton mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.17 DS/C: result fill sequence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 Validation test geometry: tapered channel . . . . . . . . . . . . . . . . . . . . B.2 Validation test geometry: tapered channel . . . . . . . . . . . . . . . . . . . . B.3 Visual comparison of FEniCS and LIMS solutions. Mesh distorted by p • e 3 . . . B.4 Solution comparison and error . . . . . . . . . . . . . . . . . . . . . . . . . . C.1 Validation test geometry: dual-scale channel . . . . . . . . . . . . . . . . . . . C.2 Comparison of the pressure solutions obtained from the in-house Matlab solver (left) and LIMS implementation (right) of the DSS model . . . . . . . . C.3 Comparison of pressure solutions from DSS model plotted on the medial line of skeleton meshes Ω1 and Ω2 . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv List of Tables 2.1 Measured reference cross-section area A 0 from an 8 samples population. . . . 2.2 Comparison of reconstructed vs measured yarn geometries [mm] . . . . . . . 3.1 Impulses collected from the spectrum F G . . . . . . . . . . . . . . . . . . . . . 3.2 List of test cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Two Ply Stack case: impulses in the spectrum (Fig.3.7a) . . . . . . . . . . . . 3.4 Three Ply Stack case: impulses in the spectrum . . . . . . . . . . . . . . . . . 3.5 Four Ply Stack case: impulses in the spectrum . . . . . . . . . . . . . . . . . . 3.6 4PPm: impulses in the spectrum (Fig.3.13a) . . . . . . . . . . . . . . . . . . . 3.7 4PPa: impulses in the spectrum (Fig.3.13b) . . . . . . . . . . . . . . . . . . . 3.8 9PPm: impulses in the spectrum (Fig.3.15a) . . . . . . . . . . . . . . . . . . . 3.9 9PPa: impulses in the spectrum (Fig.3.15b) . . . . . . . . . . . . . . . . . . . 4.1 Geometrical parameters for Gebart's expressions . . . . . . . . . . . . . . . . . 4.2 List of test cases for 1D Medial Skeleton model. ξ indicates the perc. rel. error, as defined in Eq. 4.42 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Boundary labels for the two-dimensional test cases . . . . . . . . . . . . . . . 4.4 Problem parameters of test case SR2/A . . . . . . . . . . . . . . . . . . . . . . 4.5 Channel flow: Stokes vs Reynolds solution comparison . . . . . . . . . . . . . 4.6 Results from SR2/B test case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Parameters of strongly varying gap problem . . . . . . . . . . . . . . . . . . . 4.8 Results from SR2/C test case . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Parameters for SR2/D test case . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Results from SR2/D test case: permeability and error . . . . . . . . . . . . . . 4.11 Boundary labels for the three-dimensional test cases . . . . . . . . . . . . . . 4.12 Parameters of 3D Strongly varying channel flow . . . . . . . . . . . . . . . . . 4.13 Strongly varying channel 3D: Stokes vs Reynolds solution comparison . . . . . 4.14 Geometrical parameters of doubly slanted plates test case . . . . . . . . . . . 4.15 Results from SR3/B test case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Geometrical parameters of SR3/C test case . . . . . . . . . . . . . . . . . . . . xv 4.17 Results from of SR3/C test case . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.1 Yarn geometrical parameters for all test cases . . . . . . . . . . . . . . . . . . 5.2 RT/A: Textile geometrical parameters . . . . . . . . . . . . . . . . . . . . . . . 5.3 Results from test case RT/A . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 RT/B: Textile geometrical parameters . . . . . . . . . . . . . . . . . . . . . . . 5.5 Results from test case RT/B . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 RT/C: Textile geometrical parameters . . . . . . . . . . . . . . . . . . . . . . . 5.7 Results from test case RT/C . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Characterised permeability tensor for test case RT/C . . . . . . . . . . . . . . 5.9 Fill and CPU time results from the unsaturated flow test cases . . . . . . . . . 5.10 RT/AT: Textile geometrical parameters . . . . . . . . . . . . . . . . . . . . . . 5.11 RT/BT: Textile geometrical parameters . . . . . . . . . . . . . . . . . . . . . . 5.12 RT/CT: Textile geometrical parameters . . . . . . . . . . . . . . . . . . . . . . 5.13 RT/DT: Textile geometrical parameters . . . . . . . . . . . . . . . . . . . . . . 5.14 RT/DT: Evolution of flow front position and principal permeability values . .Chapter

1 1. 2 .

 12 Figure 1.1: a) Schematic of the RTM process; b) dual-scale flow impregnating a quasi-UD textile[START_REF] Facciotto | Modeling of anisotropic dual scale flow in RTM using the finite elements method[END_REF]

Figure 2 . 1 :

 21 Figure 2.1: Examples of Prescale prints of textures

Figure 2 . 2 :

 22 Figure 2.2: Linear fit for calibration of LLW and LLLW Prescale ® films

Figure 2 .

 2 Figure 2.3: UDT400 fabric

Figure 2 . 5 :

 25 Figure 2.5: Single yarn compression device Figure 2.6: Camera recorded yarn compression

  (a) Finite element yarn compression, σ 2 stress field on a quarter of yarn. (b) pnum database for imposed final h = 0.25mm, half yarn normalized width (vertical symmetry). Units in legend are [mm]

Figure 2 . 9 :

 29 Figure 2.9: Numerical simulation of yarn compression using the hyperfoam material model

2. 5 .Figure 2 . 11 :

 5211 Figure 2.11: Pressure field analysis for the textile reconstruction

Figure 2 . 12 :

 212 Figure 2.12: Synoptic comparison of real and reconstructed geometry of a 98 × 98mm UDT400 textile

3. 2 Figure 3 . 1 :

 231 Figure 3.1: C-Weave UT710P textile

Figure 3 . 2 :

 32 Figure 3.2: Examples of superposed pairs of periodic structures that originate Moiré effect

3. 3 .Figure 3 . 3 :

 333 Figure 3.3: The function g 1 and its spectrum F 1

3. 3 .Figure 3 . 4 :

 334 Figure 3.4: The function g 2 and its spectrum F 2

3. 3 .Figure 3 . 5 :

 335 Figure 3.5: The product function G and its spectrum F G

3. 4 Figure 3 . 6 :Figure 3 . 7 :

 43637 Figure 3.6: Two Ply Stack case: pressure field

( a )Figure 3 . 8 :

 a38 Figure 3.8: Two Ply Stack case: isolated impulses and associated gratings

Figure 3 . 9 :

 39 Figure 3.9: Photographic validation of two ply stack with different orientations case

Figure 3 . 10 :

 310 Figure 3.10: Three Ply Stack case

Figure 3 . 11 : 36 3. 5 .

 311365 Figure 3.11: Four Ply Stack case

Figure 3 . 12 : 38 3. 5 .Figure 3 . 13 :

 312385313 Figure 3.12: Four Ply Periodic (4PP): pressure field at 30 kN force

3. 5 .Figure 3 . 14 :

 5314 Figure 3.14: Nine Ply Periodic (9PP): pressure field at 30 kN

Figure 3 . 15 :

 315 Figure 3.15: Nine Ply Periodic (9PP): F (u, v) spectrum data (enhanced contrast)

Figure 3 . 16 :

 316 Figure 3.16: Digital reconstruction of the Four Ply Stack preform

4. 2 .

 2 Predictive methods for permeability characterisation (a) Square packing (b) Hexagonal packing

Figure 4 . 1 :

 41 Figure 4.1: The array packings used by Gebart

50 4. 3 .

 503 Model equations along with the continuity equation, which for a constant density fluid (∂ρ/∂t = 0) is:

Figure 4 . 2 :

 42 Figure 4.2: µ-CT acquisition of a multi-layer quasi-ud layup[START_REF] Hemmer | Unloading during the infusion process: direct measurement of the dual-scale fibrous microstructure evolution with X-ray computed tomography[END_REF]

Figure 4 . 3 :

 43 Figure 4.3: Effective channel gap b d definition as used in [101]

Figure 4 . 4 :

 44 Figure 4.4: The original concept of medial skeleton from Blum's 1967 paper [103]. The skeletons are highlighted. a) Medial Axis Skeleton. b) Medial Surface Skeleton.

  26) (a) Voronoi diagram in two dimensions of a set of 7 generator points.(b) A pair of 3D Voronoi regions[START_REF] Ledoux | Interpolation as a tool for the modelling of three-dimensional geoscientific datasets[END_REF]

Figure 4 . 5 :

 45 Figure 4.5: Examples of Voronoi diagrams. Generator points are marked as p i . Lines (2D) or planes (3D) separate Voronoi regions (V i ), bounded by Voronoi vertices (v i ).

30 )

 30 The dimension (d -1) of the boundary hyperplane calls for clarification: if the original geometry is two-dimensional (d = 2), each Voronoi boundary is a line (Fig.4.5a). Likewise if the Voronoi regions are three-dimensional (d = 3), then each Voronoi boundary is a plane delimiting the region (Fig.4.5b

59 4. 4 .

 594 Medial Skeleton For the case d = 3: -C3/0: associates each Voronoi region with the bounding Voronoi vertices. Elements: three-dimensional polyhedra of varying number of vertices -C2/0: associates each Voronoi plane with its Voronoi vertices. Elements: planar polygonal regions of varying number of vertices -C1/0: associates each Voronoi ridge with its Voronoi vertices. Elements: lines with 2 Voronoi vertices

Figure 4 . 6 :

 46 Figure 4.6: Flowchart of Methodology

Figure 4 . 7 :

 47 Figure 4.7: Stokes flow constant channel domain

Figure 4 . 8 :

 48 Figure 4.8: Solution of Stokes flow in channel

Figure 4 . 9 :

 49 Figure 4.9: Mesh convergence analysis on SR2/A

  (a) Velocity profile (b) Numerical solution error with mesh size

Figure 4 . 10 :

 410 Figure 4.10: Convergence study and validation of the FEniCS Stokes solver.

Figure 4 . 11 :Figure 4 . 12 :

 411412 Figure 4.11: Channel thickness h on the medial axis mesh

Figure 4 . 13 :

 413 Figure 4.13: Domain for smoothly varying gap size test case

Figure 4 . 14 :

 414 Figure 4.14: Different inlet sizes h 0 used for the parametric study

2 .

 2 An example of the result is shown in Fig. 4.15. Results are summarised in Table 4.6.

Figure 4 . 15 :

 415 Figure 4.15: Solution fields for test case SR2/B: Stokes flow

Figure 4 . 16 :

 416 Figure 4.16: Solution fields for test case SR2/B: Skeleton model

Figure 4 . 18 :

 418 Figure 4.18: Geometry and Skeleton mesh for test case SR2/C

Figure 4 . 19 :Figure 4 . 20 :

 419420 Figure 4.19: Solution fields for test case SR2/C: Stokes flow

  Several meshes are generated for different values of V f for both the square and hexagonal packing configurations. A high number of elements are used to remove solution bias due to elements size. On each mesh, Stokes flow problem is solved following the procedure detailed in Section 4.5.2. In Figs. 4.21 and 4.22 we show an example of solutions, for the square and hexagonal packing scenarios respectively. Results are summarised in Table4.10.

Figure 4 . 21 :Figure 4 . 22 :

 421422 Figure 4.21: Solution fields for test case SR2/D: Stokes flow (Square packing, V f = 0.78)

  2. An example of results is shown in Figs. 4.23 and 4.23, for the square and hexagonal packing scenarios respectively. Results are summarised in Table 4.10.

Figure 4 . 23 :

 423 Figure 4.23: Solution fields for test case SR2/D: Skeleton model (Square packing, V f = 0.78)

Figure 4 . 24 :

 424 Figure 4.24: Solution fields for test case SR2/D: Skeleton model (Hexagonal packing, V f = 0.73)

Figure 4 . 25 :Figure 4 . 26 :

 425426 Figure 4.25: Square packing: models comparison, permeability (left) and error (right)

Figure 4 . 28 :

 428 Figure 4.28: Solution fields for test case SR3/A: Stokes flow

Figure 4 . 29 :

 429 Figure 4.29: Solution fields for test case SR3/A: Skeleton model

Figure 4 . 30 :

 430 Figure 4.30: Domains for SR3/B test case

Figure 4 . 31 :

 431 Figure 4.31: Solution fields for test case SR3/B: Stokes flow (h 1 /h 2 = 4)

Figure 4 . 32 :

 432 Figure 4.32: Solution fields for test case SR3/B: Skeleton model (h 1 /h 2 = 4)

4. 7 .

 7 Three-dimensional test cases (a) Stokes flow geometry (slice) (b) Skeleton mesh

Figure 4 . 33 :

 433 Figure 4.33: Streamlines from velocity field, h 1 /h 2 = 2

  Channel between two yarns (b) Problem geometry (c) Medial Surface Skeleton mesh

Figure 4 . 34 :

 434 Figure 4.34: Domains for SR3/C test case

Figure 5 . 2 :

 52 Figure 5.2: Schematic of the channel flow test case scenario

Figure 5 . 3 :

 53 Figure 5.3: Schematic of radial flow measurement apparatus [1]

5. 3 .Figure 5 . 4 :

 354 Figure 5.4: Schematic of yarn geometrical parameters

Figure 5 . 5 :

 55 Figure 5.5: An example of layup generated using TexGen (elliptical yarn cross sections, 4 layers)

Table 5 . 2 :Figure 5 . 6 :

 5256 Figure 5.6: RT/A: Textile geometry

5. 4 .Figure 5 . 7 :

 457 Figure 5.7: RT/A: Skeleton mesh

Figure 5 . 8 :

 58 Figure 5.8: Solution fields for test case RT/A: Longitudinal flow

5. 4 .Figure 5 . 9 :

 459 Figure 5.9: Solution fields for test case RT/A: Transverse flow

Figure 5 .

 5 Figure 5.10: RT/B: Textile geometry

5. 4 .Figure 5 . 12 :Figure 5 . 13 :

 4512513 Figure 5.12: Solution fields for test case RT/B: Longitudinal flow

Figure 5 . 15 :

 515 Figure 5.14: RT/C: Generic layup geometry and skeleton mesh

LIMS uses a

  Darcy solver to compute the pressure solution at every time step. With no source terms or Neumann conditions imposed, the equation solved by LIMS becomes:

  Figure 5.17: RT/AT: Textile and skeleton mesh

  , a series of considerations arise. One can observe that the flow front location is different in each layer, because the flow at each level of the mesh is affected by the connectivity of the N layers Orientations L [mm] W [mm] H [mm] Packing factor p in [Pa] 4 [0.0 • , 18.08 • , 4.48 • , 76.5 • ]

Figure 5 . 19

 519 Figure 5.19

Figure 5 .

 5 Figure 5.20: RT/BT: top and bottom resin filling sequences

Figure 5 .

 5 Figure 5.21: RT/CT: comparison between racetrack and no racetrack filling scenarios

  5.22).

Figure 5 .

 5 Figure 5.22: RT/DT: Textile with inlet gate (left) and skeleton mesh (right)

First, a newFigure 6 . 1 :

 61 Figure 6.1: Example of dual-scale flow in a quasi-UD fabric Interglas 92146[START_REF] Facciotto | Modeling of anisotropic dual scale flow in RTM using the finite elements method[END_REF]

Figure 6 . 2 :

 62 Figure 6.2: Yarn immersed in a channel: three-dimensional original geometry and corresponding hybrid skeleton mesh

where

  Γe indicates the boundary subset where Neumann boundary conditions are prescribed by the boundary flux term g = h 2 e /12(∇p • n). This is the classical finite element formulation for the solution of an elliptical PDE.

Figure 6 . 3 :

 63 Figure 6.3: Triangle element with interface flux

1 =

 1 internal element flux balance A.2 = interface flux term B = external source term C = Neumann flux boundary condition on element boundary ∂Ω e

  . The vector [b] contains the prescribed volumetric source terms and Neumann boundary conditions. Without loss of generality, we will consider that [b] = [0] in the test cases addressed in Section 6.7 (no source or external flux terms imposed).

Figure 6 . 4 :

 64 Figure 6.4: Interface element

Figure 6 . 5 :

 65 Figure 6.5: 1D element Ω3 connecting the domains Ω1 , Ω2 . The assumed pressure gradient is represented

6. 5 .

 5 Model test and validation-a dual-scale flow is solved using the Stokes-Brinkman model;

Figure 6 . 6 : 1 -Brinkman equation on Ω 2 -

 6612 Figure 6.6: The geometry used for DSS model test and validation

Figure 6 . 7 :Figure 6 . 8 :

 6768 Figure 6.7: DSS model validation: Stokes-Brinkman solution (velocity)

  Figure 6.10: DS/A: Input geometry and hybrid skeleton mesh

Figure 6 .

 6 Figure 6.11: DS/A: result fill sequence

Figure 6 . 12 :

 612 Figure 6.12: DS/B: Input geometry and hybrid skeleton mesh

Figure 6 .

 6 Figure 6.14: DS/B: result fill sequence

Figure 6 .

 6 Figure 6.15: DS/C: Input geometry

Figure 6 .

 6 Figure 6.16: DS/C: Hybrid skeleton mesh

Figure B. 1 :

 1 Figure B.1: Validation test geometry: tapered channel

  Figure C.1: Validation test geometry: dual-scale channel

  Figure C.3: Comparison of pressure solutions from DSS model plotted on the medial line of skeleton meshes Ω1 and Ω2

  Figure 2.4: Schematic of Yarn cross section where w is the width and h is the thickness subjected to normal distributed traction vector, p film e 2

				2.4. Methodology
	-p film e 2
	Γ film	O	Ω	h
	e 2
	e 1	w

Table 2 .

 2 1: Measured reference cross-section area A 0 from an 8 samples population.

	Yarn sample	1	2	3	4	5	6	7	8
	A 0 [mm 2 ]	1.10 1.11 1.11 1.09 1.12 1.08 1.04 1.12

Table 2 . 2

 22

: Comparison of reconstructed vs measured yarn geometries [mm] Specimen w meas 0 h meas 0 w calc 0 (err) h calc 0 (err) 1 2.94 0.35 3.192(8.5%) 0.34(2.8%) 2 3.39 0.33 3.392(0.05%) 0.32(3%) 3 2.7 0.41 2.65(1.8%) 0.405(1.2%) 4 2.7 0.41 2.77(2.6%) 0.4(4.8%) 5 3.42 0.34 3.39(1%) 0.32(5.8%) 6 2.75 0.39 2.855(3.8%) 0.38(2.6%)

Table 3 .

 3

		1: Impulses collected from the spectrum F G
	Vector	Lin. comb. Index Order Position (u, v) Position (f, θ) Amplitude	Source

Table 3 .

 3 3: Two Ply Stack case: impulses in the spectrum (Fig.3.7a)

Index Lin. comb. Order Position (u, v)

Table 3 .

 3

	2 -

4: Three Ply Stack case: impulses in the spectrum Index Lin. comb. Order Position (u, v)

Table 3 .

 3

5: Four Ply Stack case: impulses in the spectrum Index Combination Order Position

Table 4 .

 4

	1: Geometrical parameters for Gebart's expressions
	Packing Square Hexagonal π/(2 V f max π/4 √ 3) 53 16/9π A B 57 16/9π	√ √	2 6

Table 4 .

 4 2: List of test cases for 1D Medial Skeleton model. ξ indicates the perc. rel. error, as defined in Eq. 4.42 representing the medial surface skeleton of domain Ω, with associated scalar field h. An example of Medial Surface Skeleton mesh is shown in Fig.4.27b.

	Label	Description	ξ range
	SR2-A	Poiseuille flow in straight channel 0.3%
	SR2-B	Slowly varying gap channel	(0.1%, 15.7%)
	SR2-C	Strongly varying gap channel	0.6%
	SR2-D	Replica of Gebart's study [70]	(1.0%, 10.3%)
	SR2/E-H Gallery of 2D test cases	(3.2%, 6.6%)
	SR3/A	Strongly varying gap channel	4.5%
	SR3/B	Doubly slanted plates	(0.2%, 2.3%)
	SR3/C	Inter yarns channel	8.2%
	SR3/D-E Gallery of 3D test cases	(1.3%, 5.0%)
	elements mesh Ω,

Table 4 .

 4 and outlet boundaries are labeled as Γ in and Γ out respectively. All other solid boundaries are indicated as Γ w . Domains and boundaries belonging to the Skeleton domain are indicated with the |•| superscript. Depending on the geometrical dimension of the domain, these boundaries assume a different topological dimension. A summary of boundary labels and dimensions is listed in Table4.3. For the comparison to be consistent, the same set of boundary conditions must be applied to both the Stokes and Reynolds problems: 3: Boundary labels for the two-dimensional test cases

	61

Table 4 .

 4 4: Problem parameters of test case SR2/A

Table 4 .

 4 5.

	N e	CPU t [s] RMSE
	232	0.06	0.188
	1062	1.85	0.044
	3332	4.11	0.029
	9552	16.30	0.026
	20332 21.52	0.025
	26312 29.62	0.025

(a) Mesh convergence analysis: RMSE

Table 4 .

 4 5: Channel flow: Stokes vs Reynolds solution comparison Reynolds lubrication approximation does not introduce any error: the analytical expression for velocity profile from Eq.4.20 is exact and the pressure gradient in the transverse direction, as expressed in Eq.4.19, is null. This initial test case was carried out with the purpose of validating the Stokes and Reynolds solvers. However, there is little interest in comparing solutions on channels of constant cross section geometry. The method will now be extended to non-constant cross section cases.

	Scenario p in [Pa] Ψ [W]	N e .	CPU t [s] ξ
	Stokes Reynolds 3.54E03 3.538E03 298 3.55E03 3.55E03 26312 29.62 0.035	0.3%

Table 4 .

 4 6: Results from SR2/B test case

				Stokes			Skeleton
	h 0 /λ h 1 /h 0	N e CPU t [s]	Ψ [W]	N e CPU t [s]	Ψ [W]	ξ
		0.03	138202	212	1.06E+07 515	0.003	9.25E+06 12.74
		0.07	38360	55	1.88E+06 517	0.003	1.81E+06	4.03
	0.15	0.1	64354	88 69.07E+04 690	0.003 68.45E+04	0.87
		0.2	75348	119 12.67E+04 691	0.003 12.66E+04	0.10
		0.3	86864	137 48.09E+03 692	0.003 48.12E+03	0.1
		0.016	88740	160	7.50E+06 572	0.003	6.85E+06	8.72
		0.03	92000	139 13.30E+05 577	0.003 13.03E+05	2.00
	0.3	0.23	132946	213 11.11E+03 580	0.003 11.18E+03	0.61
		0.33	155252	229	4.77E+03 281	0.002	4.80E+03	0.62
		0.66	226970	349 10.01E+02 282	0.002 10.06E+02	0.50
		0.008 126184	240 53.29E+05 455	0.005 47.43E+05 10.98
		0.016 129052	110 94.47E+04 455	0.005 91.86E+04	2.76
		0.05	139900	276 61.76E+03 457	0.004 61.97E+03	0.34
	0.6	0.08	150802	257 17.59E+03 459	0.004 17.76E+03	0.92
		0.13	168150	262 56.11E+02 463	0.004 56.97E+02	1.53
		0.25	209708	367 12.44E+02 475	0.005 12.77E+02	2.64
		0.5	303550	296 24.46E+01 519	0.004	25.4E+01	3.88
		0.005 156544	190	4.14E+06 321	0.004	3.49E+06 15.77
		0.01	156544	190	7.37E+05 321	0.004	7.06E+05	4.22
	1.0	0.02 0.04	160864 168962	132 122	1.32E+05 319 2.39E+04 319	0.004 0.004	1.31E+05 2.42E+04	0.46 0.96
		0.08	185304	188	4.45E+03 319	0.004	4.55E+03	2.34
		0.2	238930	272	5.00E+02 321	0.004	5.25E+02	5.04

Table 4 .

 4

		7: Parameters of strongly varying gap problem
	L [m] H [m] µ [Pa s] Q in [m 3 /s] p out [Pa]
	0.4	0.2	1	1	0
					73

Table 4 .

 4 Section 4.2.1 the model proposed by Gebart[START_REF] Gebart | Permeability of unidirectional reinforcements for RTM[END_REF]. It is therefore interesting to see how the Medial Skeleton solution compares to the analytical results of Gebart for transverse permeability from Eq. 4.4. when the same scenario is simulated. We present a parametric study where arrays of cylinders are packed in both square and hexagonal configurations as to replicate exactly the cells used by Gebart. For each packing configuration, a series of domains are generated with different sizes of cylinders radius R, as to vary the V f values of Eq. 4.4. Geometrical parameters are summarised in Table4.9. The flow problem is then solved numerically on the full domain (Stokes) and on the Medial Skeleton (Reynolds) and the equivalent macro permeabilities are computed from the result. Values of permeability from both models are compared with the values obtained from Gebart's analytical expressions (Eqs. 4.43, 4.44). In this study the full two-dimensional result (Stokes) is the intended benchmark solution, because it is rigorous with least assumptions.

			8: Results from SR2/C test case
	Scenario	p in [Pa]	Ψ [W]	N e	CPU t [s]	ξ
	Stokes Reynolds 19.00E+03 19.00E+03 18.87E+03 18.88E+03 16365 227	16 0.004	0.6%

Table 4 .

 4

9: Parameters for SR2/D test case Packing

Table 4 .

 4 10: Results from SR2/D test case: permeability and error

Table 4 .

 4

		13: Strongly varying channel 3D: Stokes vs Reynolds solution comparison
	Scenario	p in [Pa]	Ψ [W]	N e	CPU t [s]	ξ
	Stokes Reynolds 1.788E+04 17.94E+03 71897 1.896E+04 18.78E+03 300758	1810 0.37	4.5%

Table 4 .

 4

14: Geometrical parameters of doubly slanted plates test case

Table 4 .

 4

				15: Results from SR3/B test case
			Stokes			Skeleton
	h 1 /h 2	N e	CPU t [s] Ψ [W]	N e	CPU t [s] Ψ [W]	ξ
	2	53237	175	5943 28898	0.75	5954 0.20
	3	68363	218	2552 28747	0.73	2578 1.02
	4	86260	277	1349 28359	0.73	1380 2.30

Table 4 .

 4

16: Geometrical parameters of SR3/C test case

Table 4 .

 4

			17: Results from of SR3/C test case
	Scenario	p in [Pa]	Ψ [W]	N e	CPU t [s]	ξ
	Stokes Reynolds 27.0E+03 27.0E+03 12678 24.5E+03 25.0E+03 79672	220 0.4	8.2%

Table 5 .

 5

	1: Yarn geometrical parameters for all test cases
	w [mm] e [mm]	n	d [mm] g [mm]
	4.8	0.56	0.6	0.3	0.05
	the geometry but that have not been implemented in the current library are:
	-local distortion of yarn paths
	-local variability in cross sections (non-constant through the yarn length)
					100

Table 5 .

 5

3: Results from test case RT/A

Table 5 .

 5 • , 18.08 • , 4.48 • , 76.5 • ] 7: Results from test case RT/C

			Orientations	L [mm] W [mm] H [mm] Packing factor
	4	[0.0 50	20	2.7	0.73
				Table 5.6: RT/C: Textile geometrical parameters
		θ	Ψ [W]	N elem. CPU t [s] u [mm/s]	K [m 2 ]
		0 •	1.38E-03 596768	3.76	2.33E-05	1.16E-09
		45 • 1.04E-03 497346	2.83	1.64E-05	8.20E-10
		90 • 8.08E-04 550482	3.54	1.44E-05	7.17E-10

Table 5 .

 5 09 -6.99E-11 7.13E-10 1.18E-09 7.03E-10 [0.98, 0.15] [0.15, -0.98] 8.50 • 8: Characterised permeability tensor for test case RT/C

Table 5 .

 5 9: Fill and CPU time results from the unsaturated flow test cases

	-Boundary conditions assignment
	-Solution and visualisation
	These steps are detailed below:

Table 5 .

 5 10: RT/AT: Textile geometrical parameters N layers Orientations L [mm] W [mm] H [mm] Packing factor p in [Pa]

Table 5 .

 5 11: RT/BT: Textile geometrical parameters

Table 5 .

 5

12: RT/CT: Textile geometrical parameters

Table 5 .

 5 13: RT/DT: Textile geometrical parameters

	N layers	Orientations	L [mm] W [mm] H [mm] Packing factor r 0 [mm] p in [Pa]
	4	[0 • , 45 • , 0 • , 45 • ]	100	100	2.7	0.73	0.2	5E+03

Table 5 .

 5

	1	0.008	0.006	8.0910E-10 5.1008E-10
	2	0.010	0.008	7.0262E-10 5.0820E-10
	3	0.012	0.010	7.2201E-10 5.3348E-10
	4	0.014	0.011	7.2685E-10 5.2419E-10
	5	0.015	0.012	7.3003E-10 5.0425E-10
	6	0.016	0.013	7.3722E-10 5.1644E-10
	7	0.017	0.013	7.3732E-10 5.0086E-10
	8	0.018	0.014	7.4340E-10 4.9976E-10
	9	0.019	0.015	7.4170E-10 5.1044E-10
	10	0.020	0.016	7.2829E-10 5.1543E-10
	11	0.021	0.016	7.3185E-10 4.9627E-10
	12	0.022	0.017	7.2304E-10 4.9680E-10
	13	0.023	0.018	7.2488E-10 4.9910E-10
	14	0.023	0.018	7.0768E-10 5.0158E-10
	15	0.024	0.019	7.2082E-10 5.0165E-10
	16	0.025	0.019	7.2162E-10 4.9305E-10
	17	0.025	0.020	7.2388E-10 4.9159E-10
	18	0.026	0.020	7.3372E-10 4.9657E-10
	19	0.027	0.021	7.1950E-10 5.0376E-10
	20	0.027	0.021	7.2514E-10 4.9615E-10

14: RT/DT: Evolution of flow front position and principal permeability values

Table 6 .

 6 2: Results from dual-scale validation test case u Ω 2 = 2.90E-05 m/s on Ω 2 for s 1 = 0, s 2 = 0.2 (porous domain).

	6.5. Model test and validation

Table 6 .

 6 3: Mesh and CPU time values from validation test case

	Mesh	Solver	N e	CPU t [s]
	Original domain FEniCS (SB) 63316	97
	Hybrid skeleton Matlab (DSS) 9354	0.02

Table 6 . 4

 64

	Test Case Description	N e	t fill [s] CPU t [s]
	DS/A	Yarn immersed in channel	19456	1.14	20
	DS/B a)	Yarns array -high yarn perm.	383273	19	3300
	DS/B b)	Yarns array -low yarn perm.	383273	16	3700
	DS/B c)	Yarns array -central racetrack	402404	16	6500
	DS/C	Multi-layer layup -central racetrack 124818	1	3700

: Fill and CPU time results from the DSS filling test cases

Table 6 . 5

 65

	1	0.5	0.4	0.2	0.1	0.125	1.0E-11 1.0E-11

6.7. Test cases L [mm] W [mm] H [mm

] h 1 [mm] h 2 [mm] h 3 [mm] K [m 2 ] K t [

m 2

] : DS/A: geometrical and material parameters

Table 6 .

 6 6: DS/B: geometrical and material parameters

Acknowledgements

Three-dimensional test cases

SR3/A: Steep variation in gap size

We replicate here, in its three-dimensional form, the test case SR2/B from Section 4.6.3. The domain consists in a rectangular cross-section channel defined by parallel plates and a transverse cylindrical obstacle (Fig. 4.27a). The width of the domain is set to unity to mimic quantitatively the results from the two-dimensional case. Geometrical and problem parameters are provided in Table 4.12.

Stokes solution

The domain is modeled using Gmsh. Boundary conditions 4.36, and 4.46 on Γ sym are prescribed. The solution fields p, u are plotted in Fig. 4. [START_REF] Mannsfeld | Highly sensitive flexible pressure sensors with microstructured rubber dielectric layers[END_REF]. Results are summarised in Table 4.17.

Reynolds solution

The Medial Skeleton mesh is generated from the 3D yarn domain. Given the planar symmetry of the geometry, the skeleton appears as a flat surface traversing the centerline of the channel, on which the varying gap size is defined as the field h (Fig. 4.34c). Solution is shown in Fig. 4.36. Results are summarised in Table 4.17.

Validation

There will be no three-dimensional validation of the skeleton flow problems by Stokes model comparison. This is because the large size of the textile RVEs compared to the size of the inter-yarn gaps would require a high number of elements, and consequently a high computational cost, to solve a 3D Stokes flow problem on the full domain. This kind of simulations falls within the competence of high performance computing, and the challenge is beyond the scope of this research. The extensive analysis presented in chapter 4 served specifically to validate the Medial Skeleton model, and specific features were isolated for the purpose. The set of test cases in this chapter presents geometrical features that can be assimilated to already analysed features, rendering another validation unnecessary. Furthermore, one can note that the difficulty to perform full three-dimensional flow simulations on large domains is itself the core motivation behind the Medial Skeleton model.

Governing equations

The governing equations, both in strong and weak form, are identical to the ones used in chapter 4. They are summarised here for convenience. The creeping flow of the fluid according to the Skeleton model is governed by Reynolds gap-averaged velocity on the medial surface skeleton:

Eq. 5.21 must be coupled with continuity equation:

Numerical solution

Once again the FEniCS finite elements library is used to solve the weak form of Reynolds equation and numerically compute a solution p kinematically admissible with boundary conditions 5.20. A direct LU decomposition solver is used in all the cases. From the pressure solution, velocity is computed on the elements using relation 4.22. The inlet and outlet total volumetric fluxes are also computed from the solution, to ensure that the mass is conserved:

where Γ denotes the inlet/outlet edge boundary of the Skeleton mesh (topological dimension: 1).

Numerical solution: Skeleton model

In order to perform three distinct channel flow measurements on the Skeleton mesh, three different meshes are prepared. Slices of the layup of size L=50 mm × W=20 mm are cut out from the Skeleton mesh of Fig. 5.14b, each oriented with side L aligned to directions

Skeleton flow problem is solved following procedure detailed in Section 4.5.2 with boundary conditions from 5.20. The velocity field solution is shown for all three orientations of the mesh in Fig. 5.15, with velocity streamlines highlighted.

Permeability determination

Values of averaged velocity in the direction e 1 are computed from the three solutions using Eq. 5.6. Effective permeability values are computed using Eq. 5.9. Results are reported in Table 5.7. The permeability tensor in the basis [e 1 , e 2 ] can be calculated by solving the system of equations 5.15 for K 11 , K 22 , K 12 . Using the measurement directions θ 1 = 0 • , θ 2 = 45 • , θ 3 = 90 • this yields the system of equations:

(5.25)

Solving the system of Eqs. 5.24-5.26 yields the three components of the permeability tensor of the material in the [e 1 , e 2 ] basis:

FE/CV method

Finite Elements / Control Volume methods (FE/CV) are hybrid methods that incorporate traits of both FEM and Finite Volume (FV) techniques, used to simulate filling of a domain by a fluid. In the FE/CV method the filling of a numerical mesh is simulated by discretising time in a series of time increments. At each time increment, the method consists in looping through two essential steps. In the first step, the pressure field is solved at every node using the finite elements method. The constitutive equation to relate velocity and pressure is usually Darcy's law and the continuity equation is used to write a Laplace equation of the same form as the one in chapter 4. For efficiency, the solution must be sought only in the yet unsaturated region (p = p in in the saturated region). This model of flow front progression is also known as slug flow or abrupt front [START_REF] Michaud | A review of non-saturated resin flow in liquid composite moulding processes[END_REF]. In the second step, the pressure solution is used to compute the velocity field for the current time step using Darcy's law. The fill factor of the elements at the front is updated as shown in Fig. 5.16 by calculating the net fluxes of resin entering the control volumes:

where s denotes the boundary of the control volume. After, the loop is repeated until the entire part is filled. The FE/CV method is used by commercial softwares like LIMS developed by University of Delaware [START_REF] Bruschke | A finite element/control volume approach to mold filling in anisotropic porous media[END_REF][START_REF] Simacek | Desirable features in mold filling simulations for liquid composite molding processes[END_REF].

LIMS: filling simulation software

In the present work, LIMS is adopted as framework to implement the Skeleton model. This solution was chosen because of the possibilities it offers in terms of non-standard implementations. LIMS is distributed with a Graphical User Interface (GUI) that allows to perform all

Skeleton mesh importation

To read a mesh into the solver, LIMS uses a *.dmp input file. dmp is a text format designed to contain all the necessary data for the simulation: nodal data, elements data and resin viscosity value. The *.dmp file does not contain any boundary conditions information.

The Skeleton mesh data is converted in the necessary format and written into a customcreated *.dmp file via a Python script. In addition to the mesh basic information, each element is assigned three key values for the correct implementation of the Skeleton model:

-thickness: the element thickness will affect the volume to be filled. The gap value h e for each element is assigned.

-volume fraction: this value affects the effective volume associated to the element. In the context of the pure-fluid domain this value has no meaning. The very low value V f = 1E -6 is assigned. That means the entire volume must be filled.

-principal permeability values: the values K 1 = K 2 = h 2 e /12 are assigned to each element, according to the relation discussed in Eq. 5.29. Furthermore, in the *.dmp the resin viscosity value µ is defined. For all test cases in the remainder of the chapter, the value of µ = 0.3 Pa•s is used.

Boundary conditions assignment

The same set of boundary conditions used in the first part of the chapter are maintained here. An inlet pressure is imposed at the boundary defined as Γin . An outlet pressure is imposed at the boundary defined as Γout . To define these boundary conditions, LIMS uses a scripting *.lb file, in LBASIC scripting language. The LBASIC scripting language is a customised version of the traditional BASIC language, which has been expanded with a set of functions for use in LIMS. Once all inlet and outlet nodes of the skeleton mesh have been identified, the BC assignment is carried out by writing the *.lb file is written via a Python script.

Solution and visualisation

Once input file and boundary conditions are assigned, the LIMS solver is called to compute the solution. One key factor when dealing with filling simulations is the visualisation of results. To properly monitor the process, the evolution of the flow front must be tracked. To to so, the *.lb script was tweaked to save a snapshot file of the current solution, in *.tp format, at every given interval ∆t of filling time. In Table 5.9 we report some numerical values extracted from the solutions. The filling time t fill is defined as the time when the most advanced region of the flow front reaches the outlet boundary Γout . The CPU time (CPU t) includes the time spent to write multiple *.tp snapshot files of the solution (about 100 tecplot files per simulation). The actual CPU time dedicated to the solution might be sensibly lower, however at the current implementation it was not tracked.

Values used for the test case are provided in Table 6.1. Eq. 6.35 is added with the weak forms of Stokes and Continuity equations defined in Eq. 4.38 to the system of equations to solve for p and u.

Boundary conditions A pressure difference between inlet and outlet is imposed as boundary condition, along with no-slip condition on the solid walls.

p = p in on Γ in p = p out on Γ out u = 0 on Γ w Numerical values are provided in Table 6.1

Numerical solution

The Stokes-Brinkman flow problem is solved once again using FEn-iCS. To reduce the computational cost without affecting the result of interest, only a slice of the domain is modeled and a planar problem is solved. This is possible because the geometry and the problem variables do not vary in the e 2 direction. Adopting this reduction allows to mesh the domain very finely, as to remove any mesh resolution bias from the solution. Mesh parameters are provided in Table 6.2. The velocity field solution is shown in Fig. 6.7 and the velocity in the channel direction u 1 is extracted along a section line s and plotted in Fig. 6.8. One can observe the typical velocity profile from Stokes-Brinkman model: the parabolic profile in the channel region and the plug flow profile in the porous region, connected by the characteristic boundary layer near the interface. To extract some values that can be comparable with the Skeleton model, we compute the average velocity in the e 3 direction for both fluid and porous regions. Since in the porous domain the value is affected by the (fast) boundary layer, we compute also the plug value of Darcy's velocity, indicated in Fig. 6.8 as u plug . The e 3 -averaged velocity on an interval [s 1 , s 2 ] is obtained as:

The following averaged values were obtained:

u Ω 1 = 7.72E-04 m/s on Ω 1 for s 1 = 0.2, s 2 = 0.4 (channel domain) Values are provided in Table 6.1

Numerical solution

The system of Eqs. 6.34 is solved after applying inlet/outlet boundary conditions using the Matlab solver, by direct inversion. Nodal solution p is exported and written back onto the mesh *.xdmf file to be able to visualise the results in Paraview. Phase averaged velocities are computed from the pressure solution on Ω1 and Ω2 using Eqs. 4.22 and 6.21, respectively. The velocity solution on the hybrid Skeleton obtained from the DSS model is shown in Fig. 6.9

The following velocity values are obtained:

u Ω1 = 6.75E-04 m/s on Ω1 (channel domain)

u Ω2 = 2.02E-05 m/s on Ω2 (porous domain).

The values are constant on the respective skeleton domains since the thickness is constant.

Discussion of results

A series of considerations can be drawn from the synoptic view of the Dual-Scale Skeleton and Stokes-Brinkman solutions. The domains will be examined individually.

Appendix A

Geometric Moiré

In the present technique the information contained in the O1 impulses (located on the orientation circle) was used to detect orientation and spatial frequency of plies. One remark can be made on the apparent redundancy of information carried by the spectrum: in fact the same quantities (orientations, frequency) can be obtained by applying the principles of geometric Moiré to the low-frequency O2 impulses. To demonstrate this the Two Ply Stack (2PS) case will be used as an example. The spectrum in Fig. 3.7 shows a pair of O2 lowfrequency impulses: (1,-1) and (-1,1). Under the assumption that the period T of the tow center-lines is known (or measured from O1 impulses) and remains the same in both layers (T 1,0 = T 0,1 = T ), the following relations are valid [19]:

where ∆θ and Σθ denote respectively the difference and sum of the orientation angles θ i for the two layers. The sole measurement of the polar coordinates f, θ of the second order impulse (1,-1) is sufficient to determine the tow center-line period of the base layers by combining Eqs. A.1 and A.2:

Here the first layer was viewed as the reference orientation (θ 1,0 = 0 • ), so that ∆θ = Σθ. These relations can be used to confirm the validity of the values in Table 3.3. For example plugging T = 1/f 1,0 = 1/f 0,1 = 5.12mm and ∆θ = θ 1,0 -θ 0,1 = 12.16 • in Eq. A.1 yields f 1,-1 = 0.0415, proving that the measurement of the impulse positions on the spectrum was consistent. More generally, second order impulses can complement first order impulses in providing the necessary information, especially considering that low frequency second order impulses like (1,-1) are often easy to detect in periodic Moiré patterns [START_REF] Amidror | The role of Fourier theory and of modulation in the prediction of visible moire effects[END_REF]. In the present

Appendix B

Validation of LIMS solution:

Single-Scale Skeleton model

The purpose of this section is to demonstrate that the LIMS solver is able to correctly solve Reynolds equation on the Skeleton mesh, replicating the same results as the FEniCS solver.

To do so we will be solving a simple test case, with a non-linear pressure field solution. The skeleton problem is first solved using the established FEniCS solver, which is now adopted as the reference solution. Consequently, the mesh input data is adapted to the Darcy solver of LIMS in order to solve the same equation. The solution is computed using LIMS and the result is compared with the reference.

Problem domain

In order to compare the solutions, we need a non-trivial test case that present some degree of non-linearity in the solution. To that purpose we generate a 3D tapered channel, with the gap height larger on one end. The geometry is shown in Fig

The geometric parameters are given in Table B.1. The medial skeleton mesh appears as a flat surface on which the channel gap field h e is assigned to each element. In order to have a comparable solution, we choose the set of boundary conditions from Eq. 5.20. These BCs impose a given pressure difference between inlet and outlet and are used for all test cases in chapter 5.

FEniCS solution Reynolds equation 4.41 with assigned boundary conditions from Eq. 5.20 is solved in FEniCS. Solution is plotted in Fig. B.2a. We remark that the gradient of the solution is not constant due to the varying channel gap. The purpose of this validation is to ensure that the implementation of the DSS model into LIMS is correct. We set up a test case and ensure that the exact same problem is solved: the pressure solution obtained from the custom written dual-scale skeleton FE solver (reference solution) is compared with the pressure solution obtained from the implementation of the dual-scale skeleton model in LIMS. The aim is to verify that the pressure solution is identical within reasonable tolerance, to confirm that the tweaks introduced in chapter 6 to the *.dmp input file to implement the DSS model actually do not affect the pressure solution.

The test case geometry is a dual-layer channel, in which the top layer is a clear-fluid region Ω 1 and the bottom layer represents a porous yarn (Ω 2 ). To make the pressure solution non-trivial, the thickness of the Ω 1 region is not constant. The input geometry is shown in Boundary conditions are imposed as follows:

p = p in = 1.0E + 05Pa on Γin ; -p = p out = 0Pa on Γout ;

Notice that to make the solution less trivial the inlet boundary condition is only imposed on the skeleton boundary of the channel domain ( Γ1 ), while the yarn inlet boundary is left unconditioned (natural Neumann condition).

The pressure solution is computed for the two cases using the custom-written Matlab solver and LIMS. In the case of LIMS solution, as was done for the single-scale validation, we extract the pressure field at the last time step, when the front has reached the end of the part.

Solutions from both solvers are displayed side by side in In the first part, the focus is on the acquisition of geometrical data from preforms: we propose a novel methodology based on the analysis of the pressure field experienced by a dry preform under compaction. By using a commercial pressure-sensitive film, the pressure field exerted by a stack of layers against mould walls is captured and analyzed. Taking advantage of the periodic morphology of textiles, geometric patterns revealed by the pressure field are interpreted according to spectral Moiré analysis to recover the orientation and spatial distribution of each individual layer in the stack.

In the second part, the reconstructed digital architecture of the preform is used to carry out numerical flow simulations at the scale of the yarns, to characterize permeability of the stack or directly perform filling simulations. The stack geometry is replaced by a skeletonized representation of the same, on which a two-dimensional flow problem can be solved numerically, greatly reducing the computational cost when compared to a full 3D approach. This "medial skeleton" model is first formulated in its single-scale version (flow in channels) and then extended to dual-scale (flow in channels and yarns). The model potential is illustrated through several test cases. This research establishes a pathway going from the non-destructive acquisition of data to the simulation of the dual-scale flow inside a multi-layer layup of textiles.