Jacques-Eric Bergez -Inrae
email: jacques-eric.bergez@inrae.fr

J E Bergez

E Ceschia

K Klumpp

S Reynders

From carbon neutral cropping systems To Climate Neutral Farms

Demonstrate that

A consortium of 33 partners will interactively integrate and improve existing solutions to achieve economically viable business models in farming systems through a multi-actor approach.

The project at a glance GRA workshop -Sevilla

Agriculture is part of the problem but is also part of the solution!

Innovative

Systemic Solution Spaces

17/01/2023

Overall concept

The goal of I3S is to develop business models that ensure the financial sustainability of the solutions, with an upscaling methodology.

MRV and predicting

• When designing an MRV scheme, it is important to consider several factors, including data availability and feasibility.

• Data availability refers to the amount and quality of data that is available for use in the MRV scheme.

• It is important to ensure that there is sufficient data available to accurately measure and report on the initiative, and that the data is of high quality and can be trusted.

• Feasibility refers to the practicality and viability of implementing the MRV scheme.

• Other factors that may need to be considered when designing an MRV scheme include the goals and objectives of the initiative, the target audience for the MRV reports, and any legal or regulatory requirements that must be met.

MRV requirements that should be considered

There are a number point to be considered when designing an MRV:

1. The number of models (tools) being necessary/ used to accurately "measure and report " -> it may be necessary to use multiple models or approaches.

2. Field sampling versus modeling: Depending on the goals of the initiative and the resources available, it may be necessary to use soil sampling techniques or modeling techniques (or a combination of both) to measure and report on progress.

Data quality and accuracy:

it is important to ensure that the data used in the MRV scheme is of high quality and accuracy, as this will impact the reliability of the measurements and reports.

Verification procedures:

It is important to have clear and transparent verification procedures in place to ensure the accuracy and reliability of the MRV reports.

Stakeholder engagement:

Engaging with stakeholders such as local communities and government agencies can help to ensure the success and relevance of the MRV scheme.

Legal and regulatory requirements:

It is important to ensure that the MRV scheme is in compliance with any relevant legal and regulatory requirements.

7. Using a decision tree can help to make the MRV scheme selection process more structured and systematic, and can help to ensure that all relevant factors are considered.

  Demonstrate that innovative systemic solutions have the potential to generate positive impacts by 2030  Achieving climate neutrality of farms and farming systems  Reducing GHG emissions  Increasing carbon sequestration and storage  Consider other climatic effects (albedo change, surface energy partitioning...)

  Ideal farm using map and post-its -mixed groups Activity 2 -Ideal farm using solution cards (uptake line and matrix) -mixed groups workshop 2 -How can stakeholders help farmers to implement impactful solutions Activity 3 -Barriers and pre-requisites of high impact solutions Activity 4 -What do farmers need; what can stakeholders offer Plenary Activity 5 -Plenary discussion on how stakeholders can help farmers implement impactful solutions Activity 6 -How DF can support high impact, low uptake solutions  Create collective intelligence  Empowerment of farmers  Spark discussion between stakeholder, an appreciation for each other's points of view  Gather knowledge on how we can help get solution onto farms  Lessons learnt from the first creative arena can be provided to help other I3S plan their creative for carbon sequestration & GHG reduction in agriculture is under construction at European & international level (certification rules, sources of financement...) • In this context, and among other things, a methodological framework for MRV is under construction A short focus on the MRV .L.ITY: QUantification, Additionality and baselines, Long-term storage and sustainabilITY

-

  Analyse the current/existing frameworks for MRV at different production systems (duration of the projects, definition of the baselines, levers considered...) ==> towards the proposition of a unified framework, -Contribute to the development of digital tools for monitoring by : • Analysing the compliance of the current digital tools (e.g. models) with the certification methodologies (literature study) : analysis of strength and weaknesses of current models, needs for further development, accounting or not for key climate mitigation practices ==> proposition of improvement of the current certification methodologies, • Testing a range of models and implementing them at Demo and Lead commercial farms ==> analysis of minimum dataset needed, data collection, sensitivity and accuracy analysis, ensemble modelling, • Developing dedicated protocols for data collection at farm level (e.g. soil and vegetation sampling) and best use of new data steams (e.g. remote sensing) for model's input & validation, • Developing prototypes of digital Monitoring tools for the MRV taking the advantage of new data streams (e.g. remote sensing) => increase the accuracy, the scalability and reduce de cost of implementation of Monitoring for MRV schemes 17/01/2023 14 Therefore, in strong collaboration with the OrcaSa and Marvic Horizon projects and building on the conceptual framework for Monitoring by Smith et al. (2020)

  ClieNFarms scope is based on a demonstration approach through the creation of I3S

	and
	cover crops
	17/01/2023

Farms that are experimental/research sites on which a range of existing solutions will be tested, to reduce GHG emission and increase Carbon storage and sequestration, account for other climatic effects. Commercial farms who are pioneers in innovation testing and are well connected to the DFs, contextualizing prototypes in an operational environment. Network of the surrounding farms of LCFs to dissimenate project results to empower farmers for adoption. Includes the replication of successful systems through the supply chains that are either local, national or even international.

	17/01/2023  Technico-economic blueprints Business models 

Carbon certification platform Upscaling the I3S Communicatio, dissemination & exploitation of project results 6 Synergies with other EU projects, policies and initiatives

			Policy solutions 
		hub
	Existing solutions	Knowledge	Criteria
		consolidated
		data	Models &
		Protocols 	Tools
		 Protocols validation
		Data curation

I3S

Quantified solutions

(videos, narratives, reports, socio economic, values)

Bring concept to live  Results gathered in the reservoir

Catalogue of solutions

Existing New Applicable Already applied Not applicable Could be applied Not applicable Applicable Already applied Could be applied WP5 WP6 WP4 WP3 WP2

WP1 17/01/2023 A holistic approach to climate-neutral and climate-resilient farming 9 Livestock (

including feeds management, animal management, housings, grass management…) Carbon sequestration (including soil, humic balance, impact of farming practices on soil carbon sequestration, effect of hedgerows, agroforestry, biochar, albedo…), reduction of GHG emissions from soils, albedo effects… Integral Environmental sustainability (including biodiversity, risk of nutrients leaching, ammonia emissions, soil erosion…) Crops and special crops (including fertilization, soil management, crops management, specific mitigation practices, crops rotation…) Low carbon energy production and consumption (including fuel and electricity consumption, biogas plant, wooden chips from hedgerows, photovoltaic…) Other approaches (including circular organization, governance, territorial approaches, value chain, collaboration with other stakeholders…)

  The I3S network is structured around 6 six major themes to investigate and test solutions:

	1	2
	3	4
	5	6

17/01/2023

• This includes considerations such as the resources required to implement the scheme and whether those resources are available, as well as any technical or logistical challenges that may arise.

From farm to larger scale • Two approaches

• Spatial modeling (ex Cland)

• Research approach

• Developing sound business models from farms archetypes In this informal workshop, a small group of invited policy makers and stakeholders will discuss the proposal of the legislative framework for the certification of carbon removals. First, the European Commission will present the proposal. Then ClieNFarms will provide feedback on the proposal based on the experiences in the project. Finally, we will open the discussion for all invited participants to bring in their different perspectives.

International Advisory Board / Int. Initiatives

General outputs

ClieNFarms knowledge Hub

Innovative solutions

Farmers' Network and Multi-actors ecosystem Providing sufficient, safe, nutritious, healthy and affordable food for all.

Improving the overall sustainability of food systems.

Improving the resilience of food systems to shock and stress.

Expected Impacts

Follow us on our digital channels!


Keywords: Volume integral method, air gap, magnetostatics, facet element

The presence of air gaps in magnetic modeling problems is well known. When the air gap is thin, its numerical treatment is complex because the integrals arising on the model lead to inaccuracy. This paper presents an air gap model for thin air gaps using the magnetostatic volume integral method for linear and non linear problems. The thin air gap is considered as a face region with a numerical treatment to improve the accuracy. An application to an academic case and a current transformer is reported to validate the model.

I. INTRODUCTION

V OLUME integral method (VIM) is a powerful approach to solve magnetic problems. Contrary to the finite element method (FEM), only the active regions have to be discretised, thus avoiding to mesh the air. The VIM started in 1970s with the magnetic moment method, based on uniform magnetization on each element [START_REF] Mj Newman | GFUN: An interactive program as an aid to magnet design[END_REF] but it wasn't until the last decade that it became more popular as a result of three main improvements: advanced integral formulations such as [START_REF] Carpentier | Resolution of Nonlinear Magnetostatic Problems With a Volume Integral Method Using the Magnetic Scalar Potential[END_REF] [3] (H-conforming formulations) or [START_REF] Le-Van | A Volume Integral Formulation Based on Facet Elements for Nonlinear Magnetostatic Problems[END_REF] (B-conforming formulation), matrix compression algorithms such as the fast multipole method (FMM) [START_REF] Greengard | A Fast Algorithm for Particle Simulations[END_REF] or hybrid cross approximation (HCA) [START_REF] Börm | Hybrid cross approximation of integral operators[END_REF] and the increase of computers RAM memory. A VIM based on facet elements for linear and non linear problems [START_REF] Le-Van | A Volume Integral Formulation Based on Facet Elements for Nonlinear Magnetostatic Problems[END_REF] is considered in this paper.

A frequent feature of magnetic problems is the presence of air gaps, that can be very thin, arising for example when two parts are welded. Modeling air gaps presents difficulties, there are two options to do it. The first is to not mesh the air gap, then the proximity of the two air gap faces leads to inaccurate integral computations of Green's kernel. The second option is to mesh the air gap as a volume element, this case is more accurate than not meshing the air gap but it also leads to difficult integral calculations on the mesh elements of the air gap because the two air gap faces are very close.

A third alternative is presented in this paper: a VIM model of thin air gaps considered as a face region. It avoids the numerical treatment difficulty that arises when the air gap is not meshed. It also allows to avoid a thin mesh of the air gap volume region that is computationally expensive, meshing only the face of the air gap and keeping the precision of computations. This strategy has already been done on a FEM approach for a B-conforming formulation [START_REF] Nakata | 3D magnetic field analysis using special elements[END_REF] and a Hconforming formulation [START_REF] Guérin | Three dimensional magnetostatic finite elements for gaps and iron shells using magnetic scalar potentials[END_REF], [START_REF] Ospina | Modeling of thin structures in eddy current testing with shell elements[END_REF]. Nevertheless, for the VIM case, special attention has to be given to the integral equations Manuscript received XXX; revised XXX. Corresponding author: M. Hernandez Alayeto (email: mayra.hernandezalayeto@se.com).

that model the air. The face air gap development is applied to an academic case and an industrial example and it is compared with a converged volume air gap model to show its accuracy.

This paper is organised as follows: section II explains the volume integral formulation using B-facet elements, section III describes the face air gap model developed for the previous formulation, section IV presents two applications of the face air gap model and finally a conclusion in section V completes the paper.

II. VOLUME INTEGRAL FORMULATION

We consider Maxwell's equations under magnetostatics assumptions in a domain composed of a ferromagnetic region Ω and coils. The volume magnetostatic integral equation in terms of induction B [START_REF] Le-Van | A Volume Integral Formulation Based on Facet Elements for Nonlinear Magnetostatic Problems[END_REF] is:

νB + 1 4π ∇( Ω ((ν 0 -ν)B) • ∇( 1 r )dΩ) = H 0 ( 1
)
where ν is the reluctivity of the ferromagnetic region that depends on B, ν 0 is the vacuum reluctivity, ∇ is the differential operator, r is the distance between the integration point and the observation point and H 0 is the source field created by source currents. Interpolating B with facet elements (Whitney 2-form [START_REF] Bossavit | Whitney forms: A class of finite elements for threedimensional computations in electromagnetism[END_REF]) and applying a Galerkin projection with facet functions as in [START_REF] Le-Van | A Volume Integral Formulation Based on Facet Elements for Nonlinear Magnetostatic Problems[END_REF], the resulting discretized integral system reads

(R + L){Φ} = {∆φ} + {S} ( 2
)
where Φ is the magnetic flux across each facet of the mesh, ∆φ is the difference of magnetic potential of two mesh elements, S is a term coming from the external sources, R is the finite element matrix

R ij = Ω w i • νw j dΩ (3)
where w i is the facet shape function of element i and L is the integral full matrix

L ij = 1 4π Γext i 1 s i Γj 1 s j δν j 1 r dΓ j dΓ exti (4)

An Air Gap Model for Magnetostatic Volume Integral Formulation

Mayra Hernandez Alayeto where s i is the surface of the facet element i, δν is the reluctivity jump between two adjacent elements that have a shared facet, Γ is the ensemble of facets of the problem and Γ ext refers only to the facets that are on the boundary of the ferromagnetic region. Let us notice that for the non linear case the hypothesis of ν being constant at each mesh element is considered, which is already true for tetrahedral mesh elements because in that case B is constant at each element. The non linear resolution of the system is carried out with the Newton-Raphson method.

In practice, eq. ( 2) is a representation of the problem with an equivalent electric circuit approach expressed with a dual mesh as in Fig. 1 [START_REF] Nguyen | An integral formulation for the computation of 3-D eddy current using facet elements[END_REF], where the branches of the dual mesh correspond to the facets of the primal mesh and the nodes of the dual mesh are the elements of the primal mesh. Then R matrix can be seen as the internal reluctance matrix for the magnetic material and L matrix as the one taking into account flux in air region. To solve the problem, either the Kirchhoff's nodal rule or the Kirchhoff's mesh rule can be used. This procedure ensures the free divergence of B. In order to speed the computations, a matrix compression technique needs to be used, such as FMM [START_REF] Greengard | A Fast Algorithm for Particle Simulations[END_REF] on matrix L. To that aim, we apply an iterative solver with GMRES that can be accelerated with the use of a preconditioner of type LU on the finite element matrix R.

III. FACE AIR GAP MODEL

Let us consider the following assumptions:

• Thickness of the air gap is small compared to the dimension of the other elements. • Flux leakage from the external faces of the air gap is null.

• Flux inside the air gap is perpendicular to the air gap. We can assume this because the permeability of the material is much higher than the permeability of the vacuum and the air gap is thin. Due to these assumptions, the air gap can be considered as a face region instead of a volume region (Fig. 2). The equivalent magnetic circuit can be built adding the air gap contribution to matrices R and L (eqs. ( 3) and ( 4)). For the finite element matrix R, an additional reluctance on the mesh branches that cross the air gap is considered, given by µ0e s , where e is the thickness of the air gap, s is the surface of the facet element and µ 0 is the vacuum permeability. This can be seen as the relutance of a parallelepiped of surface s and depth e. For matrix L, since the air gap is characterized by having two regions that are near, the evaluation on two parallel faces Γ 1 and Γ 2 that are in front of each other and separated by a distance equal to the thickness of the air gap (Fig. 3) is given by:

1 4π Γext i 1 s i Γ1 j 1 s j (ν 1 -ν 0 ) 1 r 1 dΓ 1j dΓ exti - 1 4π Γext i 1 s i Γ2 j 1 s j (ν 2 -ν 0 ) 1 r 2 dΓ 2j dΓ exti (5)
Fig. 3: Air gap representation.

where Γ 1j are the facet elements j that are on face Γ 1 and minus sign between the two previous terms comes from the orientation of the faces, going from Γ 1 to Γ 2 . Now, expressing ν 1 and ν 2 as

ν 1 = ν avg + δν 2 , ν 2 = ν avg -δν 2 where ν avg = ν1+ν2 2
and δν = ν1-ν2 2 , the previous expression becomes:

1 4π Γext i 1 s i δν 2 1 s j Γ1 j 1 r 1 dΓ 1j - Γ2 j 1 r 2 dΓ 2j dΓ exti - 1 4π Γext i 1 s i (ν avg -ν 0 ) Γ1 j 1 r 1 dΓ 1j - Γ2 j 1 r 2 dΓ 2j dΓ exti .
(6) The latter equation shows that when the air gap is thin r 1 and r 2 are almost the same value and then the difference between the two large terms involving r 1 and r 2 can lead to inaccuracy. This happens when the air gap is meshed as a volume element and when it is not meshed at all. Nevertheless, when it is meshed as a volume element the inaccuracy is only given on matrix L whereas if there is no air gap mesh, inaccuracy is given on both matrices R and L.

For the face air gap model, the first term can be approximated as:

1 4π Γext i 1 s i Γj 1 s j δν r dΓ j dΓ exti (7)
using a first order Taylor expansion

1 r1 + 1 r2 = 2 r + ε(O 2 ) ≈ 2 r
, where ε(O 2 ) is a term of order two.

The second term is approximated as:

1 4π Γext i 1 s i Γj (ν avg -ν 0 ) e s j ∇( 1 r ) • ndΓ j dΓ exti (8) given that ∇( 1 r ) • n = d dn ( 1 r ) ≈ lim e→0 1 r 1 -1 r 2 e
, therefore

1 r1 -1 r2 ≈ e∇( 1 r ) • n.
Assembling the previous information, matrix L can be written as:

L ij = 1 4π Γext i 1 s i Γj 1 s j δν r dΓ j dΓ exti + 1 4π Γext i 1 s i Γj (ν avg -ν 0 ) e s j ∇( 1 r ) • ndΓ j dΓ exti . ( 9
)
We observe that the first term of L corresponds to (4), which doesn't consider the air gap, being the second term the additional contribution given by the air gap.

When the thickness of the air gap is small the correction of R matrix is enough to model the air gap behaviour. Otherwise, the additional correction to L matrix must be considered.

IV. APPLICATIONS

A. Academic case

The proposed air gap model is applied to a magnetic part (linear isotropic material with permeability µ r = 1, 000) with an air gap and a coil of 1,000 turns and 1A of imposed current as Fig. 4 shows. In order to validate the face region air gap method, the calculation of the magnetic flux through the coil following [START_REF] Huang | General Integral Formulation of Magnetic Flux Computation and its Application to Inductive Power Transfer System[END_REF] is compared with the values obtained when the air gap is considered as a volume region. We create a mesh with enough number of elements to have a converged solution on the air gap volume model so that it will be the reference. This mesh has 42,334 tetrahedral elements in the magnetic region and around 1,900 elements in the volume air gap (varying depending on the air gap size). We will mesh the face air gap region with Figure 6 shows that the relative difference between the air gap as a volume region and the air gap as a face region using the contribution of matrices R and L is very low, inferior to 0.2% and it is stable for the three air gaps tested. The difference between the air gap as a volume region and the air gap as a face region using only the contribution of matrix R remains inferior to 2% for the three air gaps tested, but it increases as the air gap thickness enlarges.

The conclusion for this application is that adding the contribution to both matrices R and L provides more accurate and more stable results.

B. Industrial case

The second validation case is a current transformer composed of two coils and a ferromagnetic region with a non linear material that follows an isotropic analytic saturation of arc tangent type with two coefficients, µ r =100 and saturation magnetization = 1.2T as Fig. 7 shows.

The air gap is produced when welding the magnetic part, which is opened in order to introduce the coil inside the device, thus producing and air gap on a face region of the current transformer. The primary coil has 1 turn and an imposed current of 10A, the secondary coil has 980 turns and no current imposed. In order to validate the face region air gap method, we follow the same method as the academic case, comparing the flux through the secondary coil for a model of a volume air gap with a converged mesh (reference) and the face air gap model. The converged mesh has 39,165 tetrahedral elements in the magnetic region and around 500 elements in the volume air gap. The face air gap region has 160 elements. Results for three different air gap thicknesses of 20µm, 30µm and 50µm are given in Fig. 8 and Fig. 9. Figure 9 shows that the face and volume models provide closer results for thin air gaps, increasing the difference when the air gap is thicker . It is also exhibited that the contribution of both matrices R and L presents lower relative difference to the volume model than considering the contribution of matrix R only.

In terms of computation time of the problem, the face air gap model is less time consuming than the volume air gap model, taking around 180 seconds to solve the face air gap model (R and L contributions) and roughly 330 seconds for the volume one. The time difference might be due to the higher number of mesh elements in the volume air gap region compared to the number of elements of the face air gap region.

This application provides a second validation of the VIM face air gap model.

V. CONCLUSION

A model of thin air gaps using the VIM is proposed. It considers the air gap as a face region with a given thickness, avoiding the integration difficulties that produces not meshing the air gap or meshing it as a volume region. It is made for a B-facet volume integral formulation for linear and non linear magnetostatic problems [START_REF] Le-Van | A Volume Integral Formulation Based on Facet Elements for Nonlinear Magnetostatic Problems[END_REF].

It adds new terms on matrices R and L of the system that defines the problem, showing that for thin air gaps the contribution on matrix R can be enough, but for thicker air gaps both contributions are needed.

An application to an academic case and an industrial case are shown for different air gap thicknesses, keeping the accuracy of the model more thoroughly when the air gap is thinner.

Fig. 1 :

 1 Fig. 1: Primal mesh and dual mesh of a region.

Fig. 2 :

 2 Fig. 2: Equivalent magnetic circuit of a device for an air gap as volume region (left) and as face region (right).

Fig. 4 :

 4 Fig. 4: Academic case geometry.

Fig. 5 :

 5 Fig. 5: Flux through coil.

Fig. 6 :

 6 Fig. 6: Relative difference for different air gap thickness.

Fig. 7 :

 7 Fig. 7: Geometry of magnetic region (light blue), primary coil (dark blue) and secondary coil (yellow).

Fig. 8 :

 8 Fig. 8: Flux through secondary coil.

Fig. 9 :

 9 Fig. 9: Relative difference for different air gap thickness.

  1, 2 , Gérard Meunier 1 , Loic Rondot 2 , Olivier Chadebec 1 , Matthieu Favre 2 and Jean-Michel Guichon1

Univ. Grenoble Alpes, CNRS, Grenoble INP, G2ELab, Grenoble, France

Schneider Electric, Global Technology Eybens, France


Julien Francken
email: julien.francken@arcelormittal.com

Didier Maquin

José Ragot

Bertrand Bèle
email: bertrand.bele@arcelormittal.com

Robust data reconciliation for the determination of the control system set-points of a basic oxygen furnace

Keywords: Data reconciliation, Parameter estimation, Set-point control, Batch Process, Estimation techniques

In the steel industry, the determination of the control system set-points of batch processes is a common problem. It consists in adjusting the set-points in order to reach the given product specifications thanks to a process model. Small changes in operating conditions may impact final product quality. This is particularly true for the Basic Oxygen Furnace (BOF) where the information collected during a specific batch serves to adjust the set-points of the next batch. For being able to control that type of process, measurements must be made coherent and it may be convenient to use data reconciliation procedure. The proposed paper describes a method allowing simultaneous data reconciliation and model parameter estimation. Parameter estimation results can either be used to update the process model or to detect abnormal parameter variations due, e.g. to fouling, corrosion, degradation of parts of the process.

INTRODUCTION

Context

In the steel industry, the determination of the control system set-points of batch processes is a common problem. It consists in adjusting the set-points in order to reach the given product specifications thanks to a process model. Batch processes are characterized by prescribed processing of raw materials into refined products. The objective is to produce products of the desired quality. However, small changes in operating conditions may impact final product quality. Moreover, batch processes with the same trajectory generally exhibit some degree of batch-to-batch variation. This is particularly true for the Basic Oxygen Furnace (BOF) where the information collected during a specific batch serves to adjust the set-points of the next batch. Usually measurements collected from an industrial process contain random errors, systematic biases or gross errors due to the quality of the instrumentation chain and its environment. Based on first principles (mass and energy balances, chemical equilibrium relationships, etc), modelfiltering approach, commonly called data reconciliation may be used to provide coherent data. However, the main hypothesis underlying these methods is that the model describing the system is perfectly known which is not the case in real world. The proposed paper describes a method allowing simultaneous robust data reconciliation and model parameter estimation.

The Basic Oxygen Furnace

Basic oxygen furnace is a method of primary steelmaking in which carbon-rich molten iron is made into refined steel. The vast majority of steel manufactured in the world is produced using the basic oxygen furnace. Modern furnaces will take a charge of iron of up to 350 tons and convert it into steel in less than 40 minutes.

By blowing oxygen through molten pig iron, the carbon content of the alloy is lowered and changes the material into low-carbon steel. The basic oxygen furnace actions are scheduled as follows: molten iron from a blast furnace is poured into a large refractory-lined container. Then the container is filled with the required ingredients which quantities are computed thanks to a setup charge balance model. As the required thermal energy is produced during the process, it is relevant to maintain a proper charge balance. A water-cooled lance is so lowered down into the vessel to blow pure oxygen onto the steel and iron, igniting the carbon dissolved in the steel and burning it to form carbon oxide, causing the temperature to rise to about 1700 • C. This melts the scrap, lowers the carbon content of the molten iron and helps remove unwanted chemical elements. Other fluid components are added to form slag which absorbs impurities of the steelmaking process and emulsifies to make easier the refining process. At the end of the blowing cycle, the temperature is measured and samples are taken to give a chemical analysis.

For the BOF, the setup problem consists to determine the quantity of iron ore to add and the oxygen volume to blow to reaching the target of Carbon rate in steel and the temperature defined for each heat (batch) by the given product specification. This task is done by a static charge computation based on a model formed by comprehensive heat and mass balances. Notice that for each heat, the system state variables may have a different operating point.

Main objectives and proposed method

Data reconciliation and parameter estimation are fundamental components to real time optimization of industrial processes. Great efforts have been made to develop models for various processes. However, as mentioned by [START_REF] Rolandi | Simultaneous dynamic validation/identification of mechanistic process models and reconciliation of industrial process data[END_REF], physical and chemical process phenomena are complex and difficult to model conceptually and mathematically. Indeed, thermodynamic and transport properties and reaction rates are difficult to characterize experimentally and are subject to parametric uncertainties. Moreover, even when measurements are sufficiently numerous, raw process-instrumentation data are also corrupted by different errors. Indeed, data sets taken at a certain operation point are not model-consistent. Online optimization provides a mean to run a process near its optimum operating condition by providing real-time computed optimal set-points to the control system [START_REF] Faber | An optimization framework for parameter estimation of largscale systems[END_REF]. However, the implemented algorithms need to be fed by coherent data. The more classical way to deal with that problem consists to reconcile process data based on the minimization of measurement errors subject to satisfying the model constraints and next estimate the model parameters from these reconciled values. An alternative to this two-step approach consists to develop simultaneous strategies for data reconciliation and parameter estimation [START_REF] Rod | Iterative estimation of model parameters when measurements of all variables are subject to error[END_REF][START_REF] Tjoa | Simultaneous strategies for data reconciliation and gross error detection of nonlinear systems[END_REF][START_REF] Joe | A robust strategy for joint data reconciliation and parameter estimation[END_REF]. The proposed paper is dedicated to the presentation of such method.

On a mathematical point of view, data reconciliation is generally based on the assumption that the measurement errors have Gaussian probability density function (pdf) with zero mean. As some model parameters are uncertain, they are also modeled as random Gaussian variables. Simultaneous state variable and parameter estimates are obtained using a maximum likelihood estimation approach applied on a sliding time window of observation data. To make less sensitive the parameter estimation to measurement uncertainties, the parameters are considered as constant value on a given time window. Parameter estimation results can either be used to update process model or to detect abnormal variation due, e.g. to fouling, corrosion, degradation of parts of the process.

SIMULTANEOUS DATA RECONCILIATION AND PARAMETER ESTIMATION

The proposed paper describes a method allowing simultaneous data reconciliation and model parameter estimation in order to provide model-consistent data. In real process, model parameters are subjected to smooth evolution/degradation, that's why simultaneous state variable and parameter estimates will be extracted from a sliding time window of observation data. For each observation data window, parameters are assumed to be constant whereas state variables evolve around different operating points. This computation artifice will make less sensitive the parameter estimates to measurement uncertainties. The method of simultaneous data reconciliation and parameter estimation on a given window of length N is firstly presented and the corresponding algorithm is outlined.

Next, the practical implementation using a sliding window is described

The process model

Consider the system described by:

• a set of nonlinear equations in relation to a vector of state variables x * i and a parameter vector with true value a * .

F (x * i , a * ) = 0, i ∈ [1, N ] (1) x * i ∈ IR v , a * ∈ IR p , F : IR v+p → IR n • observation equations x i = x * i + ε xi , i ∈ [1, N ] (2) • a priori knowledge of the parameters a = a * + ε a (3)
The batch process dynamics are managed as a sequence of static constraints F , without any time correlation between the different state variable x i . The x i measurements of the true values x * i are available for the N observations. For the parameter a * , the a priori knowledge is formulated as a kind of "observation equation". The parameter is defined on the basis of a nominal value and a "pseudo-noise". This form allows to express in the same manner, the state variable and parameters probability density functions, given as follows:

p xi = 1 (2π) v/2 |V | 1/2 exp - 1 2 (x * i -x i ) T V -1 (x * i -x i ) (4) p a = 1 (2π) p/2 |W | 1/2 exp - 1 2 (a * -a) T W -1 (a * -a)
(5) with V and W the variance matrices of the state variables and parameters respectively, |.| is used for the determinant of the matrix. The measurement and parameter error distributions are assumed independent random variables.

State and parameter estimation method

The problem issue concerns the estimation of the true value of the state variables and the model parameters based on the knowledge of the measurements (2) on the given horizon [1, N ], the a priori knowledge (3) of the parameter and the process model (1). The estimation can be performed by the maximum likelihood principle. Taking into account that measurement errors are independent, the likelihood function V is the product of the probability density functions:

V = N i=1 p xi p a (6)
Estimations xi , â of the true values maximize the likelihood function V in relation to x * i and a * with respect to the constraint (1) satisfaction.

Let Φ be the Lagrange function relative to the loglikelihood function and the constraint that have to be satisfied:

Φ = ln V + N i=1 λ T i F (x * i , a * ) (7)
where λ i ∈ IR n are the Lagrange parameters.

Generally, this nonlinear optimisation problem must be solved using an iterative algorithm solution. In this case, instead of a global system linearization arised from Lagrange regularity conditions, it is advisable to previously linearize the constraint equations.

For all the observations i ∈ [1, N ], at the iteration j, we assume that a solution pair {x i,j , âj } is available. From this solution pair, estimates can be improved based on a first order Taylor series expansion in the neighborhood of this solution pair. Assuming the convergence of the estimation process, the successive solution pairs constitute a convergent series toward the problem solution.

At the iteration j + 1, for a solution pair {x i,j+1 , âj+1 }, in the neighborhood of {x i,j , âj }, the constraint first order Taylor series is:

F (x i,j+1 , âj+1 ) =F (x i,j , âj ) + G x (x i,j , âj )(x i,j+1 -xi,j ) + G a (x i,j , âj )(â j+1 -âj ) (8
) with the constraint Jacobians as follows:

G x (x i,j , âj ) = ∂F (x * i , a * ) ∂x * T i x * i =xi,j ,a * =âj (9a) G a (x i,j , âj ) = ∂F (x * i , a * ) ∂a * T x * i =xi,j ,a * =âj (9b)
For more readability, these last expressions will be denoted F i , G ix , G ia omitting the dependancy with regard to the quantities at iteration index j. At iteration j + 1, estimation problem reduces to the search of the Φ j+1 Lagrange function extremum in relation to x * i,j+1 , a * j+1 , and λ i :

Φ j+1 = Φ 0 - 1 2 (a * j+1 -a) T W -1 (a * j+1 -a) - 1 2 N i=1 (x * i,j+1 -x i ) T V -1 (x * i,j+1 -x i ) (10) + N i=1 λ T i F i + G ix (x * i,j+1 -xi,j ) + G ia (a * j+1 -âj )
where Φ 0 is the constant component of Lagrange function.

This Lagrange function presents an extremum for x * i,j+1 = xi,j+1 and a * j+1 = âj+1 if:

∂Φ j+1 ∂x * i,j+1 = 0, ∂Φ j+1 ∂a * j+1 = 0 ∂Φ j+1 ∂λ i = 0 (11)
Equation ( 11) can be transformed into the following system:

-V -1 (x i,j+1 -x i ) + G T ix λ i = 0 (12) -W -1 (â j+1 -a) + N i=1 G T ia λ i = 0 ( 13
)
F i + G ix (x i,j+1 -xi,j ) + G ia (â j+1 -âj ) = 0 (14)
From ( 12), ( 13) and ( 14), one can deduce: xi,j+1 =

x i + V G T ix λ i (15) âj+1 = a + W N i=1 G T ia λ i (16) G ix xi,j+1 = G ix xi,j -F i -G ia (â j+1 -âj ) (17)
The Lagrange parameters can be deduced from ( 15) and ( 17):

λ i = (G ix V G T ix ) -1 (G ix (x i,j -x i ) -F i -G ia (â j+1 -âj )) (18)
Substituting ( 18) in ( 16) gives:

âj+1 = a + W N i=1 S i (G ix (x i,j -x i ) -F i -G ia (â j+1 -âj )) (19
) with S i defined by:

S i = G T ia (G ix V G T ix ) -1 (20)
Let us define R as follows:

R = I + W N i=1 S i G ia (21)
Assuming that the matrix R is regular,

âj+1 = R -1 W N i=1 S i (G ix (x i,j -x i ) -F i + G ia âj ) + R -1 a (22)
Finally the expression of the state estimation is obtained from ( 15) and ( 18): xi,j+1 =

x i + V G T ix (G ix V G T ix ) -1 G ix (x i,j -x i ) -F i -G ia (â j+1 -âj ) ( 23
)
The estimations xi and â reliant to x i and a are obtained using successive iterations provided algorithm convergence (the analysis of this latter can be based on the norms of the Jacobian matrices G ix and G ia ).

Algorithm

For an observation data window of length N , an algorithm based on the method developped in the previous section could be formulated in the following way:

• Initialization j = 0, ∀i ∈ [1; N ], xi,j = x i and âj = a • Repeat
Compute F i , G ix and G ia using ( 8), (9a) and (9b). Compute âj+1 using ( 22) Compute xi,j+1 using (23) j = j + 1 until norms of G ix and G ia are below a given threshold.

• Update the parameter and variable estimates

â = âj+1 xi = xi,j+1 , ∀i ∈ [1; N ]
This algorithm gives coherent parameter and state estimates for a given observation data window of length N .

Let us now examine the practical implementation of this algorithm.

PRACTICAL IMPLEMENTATION

Estimation on a sliding window

In order to be able to monitor the time evolution of the model parameters, the proposed algorithm is implemented using a sliding window. Let us consider a first observation window of length N . Let us denote a (N ) the a priori nominal value of the parameter vector (given by the user) for this observation window. From the knowledge of x 1 , . . . , x N and a (N ) , the proposed algorithm provides the estimates x1 , . . . , xN and â(N) . The observation window is then slid from one observation. Based on the knowledge of the measurements x 2 , . . . , x N +1 and the previously estimated value of parameter vector â(N) which serves as the nominal value for this new observation data window, the algorithm is used to provide the estimates x2 , . . . , xN+1 and â(N+1) . However, in a real-time context, only the last state variable estimate xN+1 is exploited. This process, which is depicted on figure 1 is re-iterated based on the availability of new sets of measurements along the time.

!"#$ % ! & !'% ( !'& x 2 , . . . , x N +1 , â(N) → xN+1 , â(N+1) x 1 , . . . , x N , a (N ) → x1 , . . . , xN , â(N) x 3 , . . . , x N +2 , â(N+1) → xN+2 , â(N+2)

Fig. 1. Principle of estimation on a sliding window

As previously said, the determination of the control system set-points of batch processes in the steel industry is difficult to solve. Indeed, small changes in the operating conditions may impact final product quality and some degree of batch-to-batch variation exists. Therefore, model adaptation algorithm are frequently implemented for that kind of processes. Clearly, the performances of the existing model adaptation algorithm will be enhanced when fed by the coherent estimated values provided by the proposed method.

Variable scaling

In many industrial process, instrumentation chains provide access to a wide variety of different scaled process data measurements. Applications for on-line data reconciliation and optimization methods must be efficient and numerically robust. To avoid numerical difficulties and make this type of application robust, a suitable scaling method is required before any efficient estimation methods.

Indeed, the performance in term of time and convergence of algorithms based on estimation methods depends crucially on how the problem is formulated. As mentioned by [START_REF] Lid | Data reconciliation and optimal operation with applications to refinery processes[END_REF], a proper scaling ensures an important issue to this problem. Note that the scaling is a preliminary step and can be performed off-line. Thus, the computational complexity of the scaling itself is not important. Rather, the objective of the scaling is to minimize the computation time in improving the convergence speed and precision of the algorithms required by estimation or optimization methods.

An unconstrained optimization problem is said to be poorly scaled if a change in x i (resp. a) in one direction produces a much larger change in F = F (x i , a) than in another direction [START_REF] Nocedal | Numerical optimization. Springer series in operations research[END_REF]. The measure of poor scaling is not so clear in constrained optimization, but a poorly scaled model is likely to generate larger rounding errors which may degrade the performance of the algorithm [START_REF] Lid | Data reconciliation and optimal operation with applications to refinery processes[END_REF].

The scaled process model for the proposed method of data reconciliation and parameter estimation is written as: F (x i , ã) = 0 (24) with the scaled variable xi = S -1 vx x i and the scaled parameter ã = S -1 va a. The scaled model can be written as F (x, ã) = S n F (S vx xi , S va ã) where S vx , S va and S n are fixed diagonal scaling matrices.

The equation scaling factors are the reciprocal of an integer power of 10 of the value of a given term or group of terms, normally related to the scale factors of the relevant variable and parameter: S vx (j, j) = 10 -kj , k j = int(log 10 |(x) j |), j = 1, ..., v (25) where x is a mean value of a set of measurement vectors x i , ( ) j , the jth component of a vector, and:

S va (j, j) = 10 -kj , k j = int(log 10 |(a) j |), j = 1, ..., p (26)
For the S n scaling factor, we have chosen to use an integer power of 2, because residual values are already close to each other.

S n (j, j) = 2 -kj , k j = int(log 2 |(f (x, a)) j |), j = 1, ..., n (27)
This scaling method have notable impact on Jacobian matrices: Gix = S n G ix S vx , Gia = S n G ia S va (28) and variance-covariance matrices:

Ṽ = S -1 vx V S -1 vx , W = S -1 va W S -1 va (29)

Gross error rejection

The proposed method is well adapted for filtering random measurement errors. However, measurements are frequently subject to systematic errors also named "gross errors". These errors may be due to instrument malfunction, miscalibration or drift, leakage or even poor sampling. If some measurements have gross errors, precautions have to be taken to avoid biased measurement adjustments or estimates. Generally, gross errors are detected and identified using statistical tests [START_REF] Romagnoli | Rectification of process measurement data in the presense of gross errors[END_REF] or by the use of robust objective functions [START_REF] Tjoa | Simultaneous strategies for data reconciliation and gross error detection of nonlinear systems[END_REF][START_REF] Alhaj-Dibo | Data reconciliation: a robust approach using contaminated distribution[END_REF].

For the concerned application, we have just implemented an empirical approach consisting to modify dynamically the variance of the measurement suspected to contain a gross error. Indeed, the influence of each measurement in the optimization criteria is weighted by its variance.

It is then possible to adjust this variance in order to minimize the influence of abnormal data. Let us consider a given observation data window which comprises the N measurement vectors x j-N +1 , ..., x j . First, a correction ratio vector r j is calculated. Each component (r j ) k is defined by:

(r j ) k = (|x j -xj |) k V (k, k) , k = 1, ..., v (30)
The components (r j ) k are then compared to a given threshold T . When one or several components exceed the threshold, let us denote by m the index corresponding to the greater correction:

(r j ) m = max k (r j ) k , k = 1, ..., v (31)
In that case, the variance of the corresponding measurement (x j ) m is set to a bigger value (let say ten times its original value) in order to simulate the lack of that measurement for the next sliding windows.

APPLICATION TO A SIMPLIFIED BOF PROCESS

The proposed simultaneous state and parameter estimation method has been applied to a simplified nonlinear process model of a basic oxygen furnace (BOF).

The process model

Due to the random nature of chemical reaction in the non homogeneous slag or to the wear-out refractory-lined container or other many causes, the process cannot accurately be described with equilibrium balances. Therefore some of the model parameters need to be estimated simultaneously with state variables.

The model is a set of mass and heat balances. Let the nonlinear equation system be a simplified model of basic oxygen furnace described as next:

(c 1 + x * 3 )x * 2 + (c 2 -x * 6 )x * 5 + a * 1 = 0 c 3 x * 1 + (c 4 + c 5 x * 3 )x * 2 + (c 6 -a * 2 x * 6 )x * 5 + c 7 = 0 (c 8 x * 3 x * 4 + a * 3 x * 3 + c 9 x * 4 + c 10 )x * 2 +(c 11 x * 6 x * 7 + c 12 x * 6 + c 13 x * 7 + c 14 )x * 5 + c 15 = 0
The first equation represents the iron mass balance where x 2 , x 5 are material quantities, x 3 and x 6 are iron mass percentage and a 1 is a parameter corresponding to the iron loss (projection due to the oxygen blowing and to the deformation of the converter) which evolves along the time. In the same manner, the second equation describes an oxygen volume balance and the third a heat balance.

The variable measurement ranges and the measurement accuracies (assumed constant) are given in the table 1. The filtering capacity of the proposed estimator is pointed out. The more the window length is bigger, the more the parameter estimation filters the noise measurement, but, in the same time, the more the estimation is delayed. Figures 4 and5 show the measurements of the two state variables x 1 and x 6 and their estimations. As the measurements are unbiased, the estimations are closed to the measurements.

Figure 6 shows the measurements and the state estimation of the variable x 1 as well as the normalized corrective term in the presence of a gross error. Between instants 30 and 40, and between 110 and 120 the magnitude of the measurement has been multiplied by 1.3 (this situation can correspond to an error on the sensor gain). The proposed algorithm together with the gross error rejection explained in section 3.3 allow a good estimation of the

CONCLUSION

In this paper, a general methodology for simultaneous robust data reconciliation and parameter estimation for nonlinear model has been proposed. Data reconciliation based on balances was performed to obtain model-consistent measurement data and simultaneously parameter estimates. The use of a sliding window on which parameters are considered as constant values allows the desensitization of their estimations with regard measurement uncertainties. It is interesting to note that the knowledge of parameter distribution errors is not a limiting factor, indeed it would be easy to modify the algorithm in order x 1 state estimation in the presence of a gross error to run without any parameter knowledge. The robustness property with regard gross measurement errors has been introduced. Another approach using a so-called contaminated distribution in the maximum likelihood estimation [START_REF] Arora | Redescending estimators for data reconciliation and parameter estimation[END_REF][START_REF] Alhaj-Dibo | Data reconciliation: a robust approach using contaminated distribution[END_REF] is being developed. Finally the proposed method must be evaluated on real process model with real data in terms of process model adjustment along the successive heat to observe the impact of the set-points adjustment on the successive batch processes.

  Fig. 2. a 1 parameter estimation

  Fig. 3. a 2 parameter estimation

Fig

  Fig. 4. x 1 state estimation and measurement

Fig. 5 .

 5 Fig. 5. x 2 state estimation and measurement state variable x 2 . As this rejection is correct, the estimates of the different parameters a i are little affected by the presence of this bias.

Fig

  Fig. 6. x 1 state estimation in the presence of a gross error

Table 1 .

 1 Measurement ranges and accuraciesFor that example, it is desired to monitor the values of the three parameters a * 1 , a * 2 and a * 3 , the other parameters c i being assumed constant (see table3). A priori knowledge on a 1 , a 2 and a 3 parameters are gathered in table 2.To validate the proposed method, a data base composed of measurement data and true values of evolving parameters has been created for 200 measurement sets. The simultaneous state and parameter estimation algorithm described above was performed based on the measurement data, the knowledge of the model, and the parameter nominal values given in table 2. Because data come from simulated data base, we can easily observe the performance of the parameter estimation by matching parameter estimates with the true values. For the two parameters a 1 and a 2 , results are given in figure2and 3 when the observation data window length is fixed to N = 10, parameter estimates are near to the true value.

	Variable	x 1	x 2	x 3	x 4
	Min value	63593	3.76	4.7	4000
	Max value	96750	4.79	5.74	5232
	Standard deviation	3339	0.19	0.24	202
	Variable	x 5	x 6	x 7
	Min value	3.71	7.27	1.04
	Max value	4.68	12.96	2, 02
	Standard deviation	0.19	0.48	0.052


domain: spi.meca.msmeca chim phys spi

text:
Guillaume Fourrier

Alain Rassineux

François-Henri Leroy

Martin Hirsekorn

Christian Fagiano

Emmanuel Baranger

Automated conformal mesh generation chain for woven composites based on CT-scan images with low contrasts

Keywords: Woven Composites, Computed tomography segmentation, Mesoscopic scale, Manufacturing defects

fiber disorientation and yarn splitting defects

Introduction

Thermoplastic composites are increasingly used in the aeronautic industry for numerous applications. Their thermoforming capability, and the absence of curing step [START_REF] Rozant | Drapability of Dry Textile Fabrics for Stampable Thermoplastic Preforms[END_REF] offers the possibility of automated processes, which is a major goal in order to reduce costs and improve production rates. However, the process automation may create defects like local tow disorientations in the pieces. These defects can impair the material strength of structures.

Modeling tools are therefore needed to gain a better understanding of the impact of these defects on the mechanical behavior.

Composite materials are modelled at three distinct scales : the microscopic scale, at which fibers and matrix are treated separately, the mesoscopic scale, at which the reinforcement architecture is defined, and the macroscopic scale, at which the structure is modeled. Most of the manufacturing defects impact the shape of the reinforcement and thus influences strains and damage evolution at the mesocopic scale [START_REF] Green | Mechanical Modelling of 3D Woven Composites Considering Realistic Unit Cell Geometry[END_REF][START_REF] Varandas | On the Importance of Nesting Considerations for Accurate Computational Damage Modelling in 2D Woven Composite Materials[END_REF]. One of the main challenges in mesoscopic modeling of textile composites is to get a high quality finite element mesh (FEM) describing with accuracy the tow shapes and positions. The different methods of generating mesoscale models of textile composites are distinguished by their fidelity with the real material [START_REF] Gereke | A Review of Numerical Models for 3D Woven Composite Reinforcements[END_REF].

Reinforcement geometries based on analytical descriptions of the tow paths are well suited for parametric analyses of woven composite materials [START_REF] Long | 8 -Modelling the Geometry of Textile Reinforcements for Composites: Tex-Gen[END_REF][START_REF] Lin | Automated Geometric Modelling of Textile Structures[END_REF][START_REF] Wendling | Meshing Preprocessor for the Mesoscopic 3D Finite Element Simulation of 2D and Interlock Fabric Deformation[END_REF]. In these models, idealized tow paths and shapes are used. [START_REF] Wendling | Meshing Preprocessor for the Mesoscopic 3D Finite Element Simulation of 2D and Interlock Fabric Deformation[END_REF]. To improve generated tow shapes, analytical models can use experimental observations to identify the input parameters [START_REF] Hivet | Consistent 3D Geometrical Model of Fabric Elementary Cell. Application to a Meshing Preprocessor for 3D Finite Element Analysis[END_REF][START_REF] Yan | Geometric Modeling of 3D Woven Preforms in Composite T-joints[END_REF][START_REF] Rinaldi | Generating Virtual Textile Composite Specimens Using Statistical Data from Micro-Computed Tomography: 3D Tow Representations[END_REF], or add mechanical regularization as in WiseTex based approaches [START_REF] Lomov | Textile Geometry Preprocessor for Meso-Mechanical Models of Woven Composites[END_REF][START_REF] Verpoest | Virtual Textile Composites Software WiseTex: Integration with Micro-Mechanical, Permeability and Structural Analysis[END_REF][START_REF] Wintiba | An Automated Procedure for the Generation and Conformal Discretization of 3D Woven Composites RVEs[END_REF]. While these methods are fast in generating virtual textile composites, the accuracy of the material representation is limited (inacurrate volume fractions or overestimated material strength) [START_REF] Varandas | On the Importance of Nesting Considerations for Accurate Computational Damage Modelling in 2D Woven Composite Materials[END_REF][START_REF] Sevenois | Avoiding Interpenetrations and the Importance of Nesting in Analytic Geometry Construction for Representative Unit Cells of Woven Composite Laminates[END_REF]. Moreover, interpenetrations between tows are possible and have to be corrected before meshing (like, e.g., in [START_REF] Wendling | Meshing Preprocessor for the Mesoscopic 3D Finite Element Simulation of 2D and Interlock Fabric Deformation[END_REF][START_REF] Rinaldi | Generating Virtual Textile Composite Specimens Using Statistical Data from Micro-Computed Tomography: 3D Tow Representations[END_REF]).

More accurate descriptions of woven reinforcements can be obtained from simulations of the manufacturing process [START_REF] Isart | Internal Geometric Modelling of 3D Woven Composites: A Comparison between Different Approaches[END_REF][START_REF] Wang | Digital-Element Simulation of Textile Processes[END_REF][START_REF] Stig | Spatial Modelling of 3D-woven Textiles[END_REF][START_REF] Grail | Consistent Finite Element Mesh Generation for Meso-Scale Modeling of Textile Composites with Preformed and Compacted Reinforcements[END_REF][START_REF] Mahadik | Finite Element Modelling of Tow Geometry in 3D Woven Fabrics[END_REF][START_REF] Durville | Simulation of the Mechanical Behaviour of Woven Fabrics at the Scale of Fibers[END_REF]. However, the result is still an approximation of the real geometry. Furthermore, these approaches are based on an idealized description of the reinforcement that does not take into account variability of the reinforcement architecture or manufacturing defects.

Another method for producing meso-scale virtual textile composites, is to use nondestructive imaging of the material, especially micro-computed tomography (µCT).These models provide a better representation of the variability of the textile structure and thus allow more accurate mechanical simulations [START_REF] Naouar | 3D Composite Reinforcement Meso F.E. Analyses Based on X-ray Computed Tomography[END_REF]. In most of the published works, representative elementary volumes (REV) of woven composites are reconstructed from a CT-scan in two steps. In the first step, warp and weft tows are separated from each other. The second step consist of the individual labeling of the tows and the identification of the local fiber orientations.

Warp and weft tows can be separated manually [START_REF] Bale | Characterizing Three-Dimensional Textile Ceramic Composites Using Synchrotron X-Ray Micro-Computed-Tomography[END_REF], which is a tedious task if the weaving pattern is complex. Therefore, various automated CT-scan processing methods have been proposed. A common approach is to use the structure tensor [START_REF] Naouar | Meso-Scale FE Analyses of Textile Composite Reinforcement Deformation Based on X-ray Computed Tomography[END_REF][START_REF] Mazars | Damage Investigation and Modeling of 3D Woven Ceramic Matrix Composites from X-ray Tomography in-Situ Tensile Tests[END_REF], which works well if the image contrast is high and the fibers can be seen individually. The required resolution can be reached with synchrotron X-rays or if the scanned volume is small enough, which for standard laboratory tomographs implies that it does not contain a full REV. At this scale, individual fibers can be tracked to reconstruct the fiber bundles composing the tows [START_REF] Auenhammer | Automated X-ray computer tomography segmentation method for finite element analysis of non-crimp fabric reinforced composites[END_REF].

For lower resolutions and for images that present noise and low contrast, the structure tensor is inadequate for tow segmentation. If there is no fibrous pattern in the tows, the Gray Level Co-Occurence Matrix can be used [START_REF] Naouar | 3D Composite Reinforcement Meso F.E. Analyses Based on X-ray Computed Tomography[END_REF]. Recent works have shown the benefits of machine learning algorithms with this process [START_REF] Sinchuk | Variational and Deep Learning Segmentation of Very-Low-Contrast X-ray Computed Tomography Images of Carbon/Epoxy Woven Composites[END_REF][START_REF] Ali | Deep learning based semantic segmentation of micro-CT images for creating digital material twins of fibrous reinforcements[END_REF][START_REF] Blusseau | Instance segmentation of 3D woven fabric from tomography images by Deep Learning and morphological pseudo-labeling[END_REF]. They can be applied to CT-scans with lower resolutions, because they consider larger features of the reinforcement, such as the yarns shapes and boundaries which are difficult to take into account with previously cited methods.

On low contrast CT-images, Sinchuk et al. [START_REF] Sinchuk | Variational and Deep Learning Segmentation of Very-Low-Contrast X-ray Computed Tomography Images of Carbon/Epoxy Woven Composites[END_REF] obtained a more accurate segmentation with a machine learning algorithm than with previously cited methods. However, in order to implement this machine learning segmentation, a manual pre-segmentation is necessary to obtain the required training data. This step is the most time consuming for the user and the number of labelled images directly impacts the accuracy of the segmentation [START_REF] Sinchuk | Variational and Deep Learning Segmentation of Very-Low-Contrast X-ray Computed Tomography Images of Carbon/Epoxy Woven Composites[END_REF].

The labelling of the tows presents the challenge of identifying the contact zones between parallel tows. Topological information on the weaving pattern or the tow shapes may be needed [START_REF] Mazars | Damage Investigation and Modeling of 3D Woven Ceramic Matrix Composites from X-ray Tomography in-Situ Tensile Tests[END_REF][START_REF] Chapoullié | Orientation-Guided Two-Scale Approach for the Segmentation and Quantitative Description of Woven Bundles of Fibers from Three-Dimensional Tomographic Images[END_REF]. The difficulty of the problem increases with the number of contact zones, or the compactness of the reinforcement. A solution is to fit a parametrized description of the geometry to the tomography images using correlation [START_REF] Bénézech | Variational Segmentation of Textile Composite Preforms from X-ray Computed Tomography[END_REF][START_REF] Mendoza | Differentiating 3D Textile Composites: A Novel Field of Application for Digital Volume Correlation[END_REF]. Global descriptors of the textile architecture (like the relative positions of the tows) make these approaches less sensitive to the noise and therefore suitable at lower resolutions. However, convergence is only achieved if the initial model is already close to the real geometry, but even if a manual process may be necessary to determine the initial model, it can be less accurate and thus less labor-intensive compared to a complete manual segmentation. The computational cost is directly proportional to the size of the volume and to the number of parameters chosen to describe the geometry [START_REF] Bénézech | Variational Segmentation of Textile Composite Preforms from X-ray Computed Tomography[END_REF].

Creating high-quality conformal FE meshes of the generated geometries is a complex task [START_REF] Grail | Consistent Finite Element Mesh Generation for Meso-Scale Modeling of Textile Composites with Preformed and Compacted Reinforcements[END_REF]. The main difficulties lie in the multiple contact zones between tows and the sharp edges of the matrix pockets. These edges will prevent most meshing algorithms from converging. This problem can be circumvented by the inclusion of thin matrix layers between the tows, which makes the use of automated meshing tools possible [START_REF] Stig | Spatial Modelling of 3D-woven Textiles[END_REF][START_REF] Lomov | Meso-FE Modelling of Textile Composites: Road Map, Data Flow and Algorithms[END_REF]. Nevertheless, small elements are required to mesh the thin matrix layers, which increases the mesh size and may induce low quality elements. Another solution is to keep a non-conformal mesh between the tows and apply tie constraints to transfer the mechanical loads between them [START_REF] Jacques | Application of Periodic Boundary Conditions on Multiple Part Finite Element Meshes for the Meso-Scale Homogenization of Textile Fabric Composites[END_REF]. Voxel meshes can easily be obtained even for very complex geometries. They are sufficient for the determination of the homogenized elastic behavior, but locally artificial stress concentrations are generated that may interfere with nonlinear phenomena such as damage initiation [START_REF] Doitrand | Comparison between Voxel and Consistent Meso-Scale Models of Woven Composites[END_REF][START_REF] De Carvalho | Reducing the Domain in the Mechanical Analysis of Periodic Structures, with Application to Woven Composites[END_REF].

In a recent work, octree-refined voxel meshes are smoothed to generate conformal meshes [START_REF] Matveev | Efficient Meshing Technique for Textile Composites Unit Cells of Arbitrary Complexity[END_REF]. This technique provides an accurate description of the tow surfaces and reduces significantly the mesh size. Yet, the element size of the voxel mesh is kept on the tow surfaces, which leads to a large number of elements to accurately represent the geometry. Moreover, octree refinement is efficient if the meshed surface is small compared to the enclosed volume, which is not the case for thin tows.

In this article, a complete and mostly automated processing chain is presented to generate representative and realistic mesoscopic FE models of woven composite unit cells with or without defects. The approach is based on a segmentation procedure of µCT-scans, with a novel tow labeling and meshing method. Section 2 introduces the studied material and presents the difficulties of processing the µCT-scans, which present ring artifacts and low contrast images. A strategy is proposed in section 3 to overcome these issues and to segment the matrix and the tow orientations. The tow labeling of the segmented REV is addressed in section 4. The resulting voxelized REV is smoothed, meshed and then coarsened using a technique proposed by Rassineux [START_REF] Rassineux | Robust Conformal Adaptive Meshing of Complex Textile Composites Unit Cells[END_REF] (section 5), which allows for element quality control throughout the whole process. In section 6, the presented approach is applied to the material with defects.

Material

The material used in this work is a woven composite composed of a 5 harness satin weaving (5HS) of carbon fibers embedded in a thermoplastic PEEK matrix supplied by CETIM. The layup sequence is [0/90] 4s . The average ply thickness is 0.3 mm for a total thickness of 2.4 mm, and the unit cell (UC) size, which contains a full REV, is approximately 6.6 × 6.6 mm 2 . Each tow contains 3k fibers. The material was scanned with a NSI X50 tomograph, with a source current of 100µA and a source voltage of 120kV. The reconstruction is based on 2200 projections (average of 20 frames taken at 7fps). Three contrast adjusted slices of the µ-CT are shown in Fig. 1b. The total scanned volume of 25 × 20 × 2.4 mm 3 limits the voxel size to 7µm. The reinforcement geometry is identified on a sub-volume of 1.2 times the UC size, corresponding to 1286x1286x307 voxels. The cell size was chosen in order to treat the tomography with a regular laptop (Macbook Pro with 8Gb RAM). In Fig. 1b we can see that despite the low resolution, the identification of the tow orientation can be identified based on the fibrous texture visible in the in-plane slices. It is easier to identify the tow boundaries on in-plane images rather than on through-thickness slices. Therefore, the in-plane slices of the REV were processed during the segmentation.

The signal to noise ratio (SNR) provides a good descriptor of the image quality. For both matrix and tows this ratio is below 4. In photography, a minimum SNR of 10 is considered an acceptable image quality [START_REF]Photography -Digital Still Cameras -Determination of Exposure Index, ISO Speed Ratings, Standard Output Sensitivity, and Recommended Exposure Index[END_REF]. The gray level distibution of the whole CT-scan (Fig. 1c) shows a significant overlap between the ranges of the matrix and the tows, which makes segmentation by thresholds difficult. Another issue of the tomography is the presence of ring artifacts. These are CT hardware related and inherent to the detector calibration. They emerge as circular lines in the images that interfere with the fiber texture. In our case the ring artifacts are parallel to the weft tows. They are visible in the red rectangle in the through-thickness cross-section normal to the weft direction (Fig. 1b).

The ring artifacts in our images prevent the processing of the tomography with a conventional algorithm based on the structure tensor, which is computed at each pixel of the image I(x, y) by:

T = K σ *   (∇ x I) 2 ∇ x I∇ y I ∇ x I∇ y I (∇ y I) 2   (1)
where ∇ x I, ∇ y I are the image gradients along the x and y axes, respectively. The application of a Gaussian convolution K σ of standard deviation σ makes the structure tensor less sensitive to noise. The eigenvector of the structure tensor associated to its maximum eigenvalue corresponds to the main texture orientation. The coherence, computed with the eigenvalues (λ 1 > λ 2 ), gives a confidence indicator of the direction given by the structure tensor.

Coh = λ 1 -λ 2 λ 1 + λ 2 (2)
In case of a unidirectional texture (λ 1 > 0 and λ 2 = 0) the coherence reaches its maximum of Coh = 1, for an isotropic texture (λ

1 = λ 2 ) it is Coh = 0.
Fig. 2b shows the orientation obtained with the structure tensor calculated on an inplane slice of the CT-scan. The texture quality in the image leads to inaccurate orientations especially in the warp tows. Even after binarizing the orientation map by a threshold at 45 • , the result is not satisfying (Fig. 2d). During the scan, the specimen is rotated around the weft direction (90 • ). Ring artifacts appear in the warp tows as light lines oriented parallel to the weft tows. As a result, there are two perpendicular line textures in the warp tows.

For a more accurate segmentation, directional filtering is needed as shown in section 3.

(a) (b) (c) (d)

Segmentation of the tow orientations

The first step of the proposed segmentation procedure is to separate the matrix from the tows. In-plane slices of the CT-scan are processed one by one. The CT-scan image 2a

presents an illumination gradient that prevents the use of a single gray level threshold to segment the matrix. Chen et al. [START_REF] Chen | Illumination Compensation and Normalization for Robust Face Recognition Using Discrete Cosine Transform in Logarithm Domain[END_REF] considered the illumination gradient as a low-frequency variation of the gray level of the image. They applied a filter based on the discrete cosine transform (DCT) to correct the illumination. The result obtained with this method is shown in Fig. 3a. The corrected image is then binarized at an optimum threshold determined with the Otsu method (Fig. 3b). Due to the overlap between the gray level distributions of the tows and the matrix (Fig. 1c), not only the pixels associated with the matrix are selected.

The segmented image is therefore cleaned by a sequence of operations. First, the holes smaller than 300 pixels are removed in the segmented image and its complement. Then a disk kernel of 2 pixels is used to do a morphological dilation of the image, followed by a gaussian blur (with a sigma of 2 pixel) and a binarization. Finally, the result presented in To separate the warps (0 • ) and the wefts (90 • ), a directional filtering is implemented to deal with the ring artifacts presented in section 2. The aim is to segment one of the two orientations accurately. The second one can then be determined by complementarity with the matrix. The real part of a Fast Fourier transform (FFT) of the image in the frequency domain is shown in Fig. 4a. The two main tow orientations appear as two clearer lines. The 90 • orientation is then filtered from the image using the mask shown in Fig. 4b (white=0, black=1), as it corresponds to the direction of the ring artifacts in the inplane slices. An inverse FFT transforms the image back to the spatial domain (Fig. 4c).

As a result, the ring artifacts are removed from the warp tows, while the texture of the weft tows becomes isotropic noise. As the coherence for an isotropic texture is close to 0 whereas for a unidirectional texture it is close to 1, this choice of filtering maximizes the difference in coherence between the two orientations. The pixels belonging to the warp tows are then those with coherence above a threshold determined with the Otsu method. The pixels belonging to the weft tows are the remaining pixels that are neither warp nor matrix.

The result for one in-plane slice is presented in Fig. 4d.

Filtering the orientation with the warp tows instead would not remove the artifacts, which in this case would be perpendicular to the filtered direction. The result would be a noisy texture with artifacts for the warp tows and an untouched weft texture, which would yield a significantly smaller difference in coherence. After segmenting the in-plane slices one by one, they are assembled to reconstruct the 3D

REV. The through-thickness section of the reconstructed REV shown in Fig. 5b evidences some remaining segmentation errors. In some zones, the texture quality is so degraded that the directional filtering gives an incorrect segmentation. This is the case for example at the bottom left corner of Fig. 5b, which is the closest to the center of rotation of the specimen during the CT-Scan (red dot). These zones have to be corrected by contouring manually the faulty zone in several cross-sections. They cover less than 0.5% of the total volume. The visual verification of the segmentation and the manual correction takes less than 20 minutes.

Then a post-processing step is applied to correct the smaller topological inconsistencies, which are mainly small included volumes of an incorrect phase visible as white areas in the gray tows in Fig. 5b. The cleaning is a 2D iterative process in the 3 directions. In each iteration, for each class (weft, warp, and matrix) the inclusions are detected and filled automatically with the corresponding material. The result of the corrections is shown in Fig. 5c. Treating the constituent materials separately in this process avoids the removal of small matrix pockets. This step corrects 0.02% of the total volume after 3 iterations and removes all the visible inclusions.

Since the in-plane slices are segmenting independently, small discontinuities between neighboring slices can appear. A 3D Gaussian filter with a standard deviation of 1.5 pixels is applied to smooth the outlines of the tows and to remove the discontinuities. The result after all cleaning steps is shown in Fig. 5c. The segmentation quality can be assessed visually in Fig. 5d. The tows superpose well with the tomography image and their shapes are well respected. The full segmented REV is presented in Fig. 6.

Visually comparing the shapes of the tow cross-sections between the segmentation and the CT-scan gives a first qualitative validation. For a more quantitative comparison, a reference volume or "ground truth" (GT) is determined by manual segmentation of the tomography. Usually, tomographies of woven composites are segmented manually by interpolating the tow surfaces between tow cross-section contours traced manually on several throughthickness slices. However, due to the low image quality, the tow cross-section boundaries are not clear enough to be precisely followed manually, as shown in section 2. The manual segmentation is therefore carried out on the in-plane slices of the REV. Each slice is segmented making the process long but more accurate. A Gaussian filter is then applied on the resulting volume to smooth the contours of the tows.

The quality of the segmentation is quantified by the confusion matrix presented in Table 1. It presents the overlapping of the predicted classes with the manually segmented reference.

The presented segmentation approach tends to overestimate the volume of the matrix in the REV. This is a consequence of the overlapping of the gray-level ranges shown in Section 2. The studied material contains the same amount of warp and weft tows. However, the presented segmentation method yields a lower volume fraction of weft tows compared to the warp tows. This difference is caused by the directional filtering process, which mainly impacts the weft direction. It induces a slight loss of sharpness of the tow edges, which adds to the edge smoothing induced by the regularization term used to compute the coherence.

This slightly reduces the tow cross-sections leading to an overestimation of the matrix volume fraction. The warp orientation, which is extracted by the directional filtering, is slightly oversegmented. Therefore, the weft direction, which is deduced by subtracting the matrix from the complement of the warp tows, is underestimated. Despite this difference in volume fraction, the shape of the tows is respected as shown in Fig. 5, which was the main goal of our approach.

The volume fractions only give a global indication of the segmentation quality. In order to analyse the geometric accuracy locally, the edges of 9 tows are segmented manually on 20 slices in each of the two orientations. The cross-sections are used to extract geometric features of the tows (area, width, thickness). This operation is realised on the automatically segmented volume and on the manually segmented reference. The statistic analysis of these features is presented in Fig. 7a. As expected, the area of the weft tow cross-sections is smaller in the automatically segmented volume (median about 17% lower). Accordingly, the area of warp tows is slighty larger (median increased by 3,5%). The trends are the same for the width and the thickness of the tows with a decrease of 12.2% and 6% respectively.

The differences in volume fractions thus correspond to a global shrinking of weft tows, and a sligth inflation of warp tows compared to the reference model.

Separation of parallel tows

In the case of small contact zones, the tows can be separated by identifying isolated tow sections and by joining these sections across the contact zones [START_REF] Chapoullié | Orientation-Guided Two-Scale Approach for the Segmentation and Quantitative Description of Woven Bundles of Fibers from Three-Dimensional Tomographic Images[END_REF]. Another possibility is to determine the planes between the plies by analyzing the matrix distribution in the REV [START_REF] Sinchuk | Geometrical and Deep Learn-ing Approaches for Instance Segmentation of CFRP Fiber Bundles in Textile Composites[END_REF].

This plane is then used to identify the contact zones between tows belonging to different plies. Separating tows of the same ply is simple if the contact area is small, as for example in a plain weave architecture [START_REF] Sinchuk | Geometrical and Deep Learn-ing Approaches for Instance Segmentation of CFRP Fiber Bundles in Textile Composites[END_REF]. In the segmented REV obtained in section 3, there are more than 1000 distinct large contact zones. Hence, there are only very few isolated sections of tows, which are insufficient for labeling. Moreover, in the material with defects, the plies cannot be separated by planes due to the out-of-plane waviness introduced by the defects. In Fig. 8 the strategy used to separate parallel tows is presented. The first main step is to retrieve the center line of the tows by a combination of tracking and fitting. Then, these lines are inflated within the pre-segmented geometry in order to determine the contact zones between tows and to label each voxel of the REV. To separate the tows, the first step in our process is to retrieve the center path of each tow. This can be done manually by pointing the centers of several cross-sections. However, this operation is time-consuming on a full REV, which contains 140 tows in our case. In order to trace these paths automatically, we use the distance transform on curved space (DTOCS) algorithm [START_REF] Ikonen | Shortest Routes on Varying Height Surfaces Using Gray-Level Distance Transforms[END_REF], which was originally developed to find the shortest route between two points on a 2D gray-level height map. We implemented the DTOCS algorithm in 3D.

In order to reduce the cost of computing the distance transform map, each tow is tracked individually in a box of 3 tows width by 3 tows height centered at the midpoint between the two seeds points used for the tracking. As the algorithm is designed for following the path of the smallest gradient, the matrix and the perpendicular tows are filled with a white noise that increases the gradient. In this way, we ensure that the tracked path stays in the desired tow. The transition between the noise and the inner part of the tow is smoothed by a Gaussian convolution, which increases the robustness of the tracking and further pushes the path towards the tow center. An example of this process is presented in the Fig. 9a,b. The seed points are the manually determined tow centers on each side of the REV.

Corresponding starting points on opposite sides of the REV are identified based on the assumption that the topology (the relative position of the tows compared to each other) is the same on both sides. Each tracking of tows is performed independently, and only the start point of the considerd tow are kept (blue dots on Fig. 9b).

The DTOCS algorithm yields the shortest path through the tows. A distance map is computed based on the gray-levels and the distance between adjacent voxels (Fig. 9c). This map represents, in each voxel, the length of the shortest path between the two end points through the considered voxel. The shortest path shown in Fig. 9d is obtained by selecting in the distance map all the voxels with the minimum value. The algorithm accurately tracks 95% of the tows. However, in the weft direction between the two symmetric center plies, some tracking errors occur, because the contact between tangent tows can be 4 times longer than the tow width in this zone. It is therefore possible that the tracking algorithm finds a shorter path that passes from one tow into another one. These specific paths (6 of the 140 tows) have to be manually corrected by pointing the center of the tows on several crosssections within the incorrectly tracked zone and interpolating the path over the REV. This operation takes approximately 10 min. The final result is shown in Fig. 10. Since the obtained shortest paths do not correspond with the center lines of the tows, they are adjusted by fitting circles to rescaled cross sections of the tows. To do so, the in-plane dimensions of the REV are reduced by a factor of 8.7, which is equal to the average ratio between tow width and thickness in the studied material. This rescaling operation leads to tow cross sections that are approximately circular (Fig. 11a). Fig. 11b shows the tracked tow paths in the rescaled REV surrounded by circles of a diameter equal to the average diameter of the tow cross sections in Fig. 11a. The positions and radii of the circles are then adjusted by minimization of a cost function

f (x) = e corr (x) + w inter e inter (x) (3)
where e corr is similarity error calculated using the Jacquard index between the mask of the idealized circular cross-sections and the segmented slice. e inter is a penalizing term to control interpenetration, which corresponds to the percentage of pixels within the circles that belong to more than one circle. w inter is the weight of the interpenetration correction in the optimization process, which is set to 2 for our case. The result for one slice of the rescaled REV is shown in Fig. 11c. A similar approach was used by Benezech in order to fit an analytical model to a CT-scan [START_REF] Bénézech | Variational Segmentation of Textile Composite Preforms from X-ray Computed Tomography[END_REF]. This procedure is repeated on every fifth slice of the tomography in both warp and weft direction, as a compromise between computational cost and precision. The center line of each tow is then obtained by spline interpolation. The coordinates are then rescaled to the initial REV size.

The center lines are used as starting points of an inflation procedure to fill the REV with labeled tows. To accelerate the procedure, the resolution of the voxel mesh is reduced by a factor of 4 in the in-plane directions and by a factor of 2 through the thickness. This reduces the 470M voxels of the initial REV to 14M. To each voxel V in the coarse grid, we assign the constituent material with the highest number of voxels in V in the fine grid. Through the following steps the voxels are considered cubic despite the anisotropic size reduction of the volume. The initial ratio of the volume is recover at the end of the mehing procedure.

We suppose that the shapes of the tow cross sections are approximately elliptical (which is confirmed by the approximately circular cross sections in the rescaled REV shown in Fig. 11a), and that the deviations are mainly caused by contact interactions with adjacent tows.

The inflation process is therefore based on elliptical cross sections with half axes a and b that are normal to the main direction of the tow (warp or weft). The values of a (1.5mm) and b (0.16mm) are obtained from measurements of the tow width and thickness, respectively, on several slices of the CT-scan. The ellipses are discretized in steps of 0.1 times the voxel size. Likewise, the splines interpolating the tow center lines are also discretized in steps of 0.1 times the voxel size. This gives a regular grid of an idealized tow surface with a constant elliptical cross section. The tows are inflated by iteratively increasing an inflation ratio e starting from 0 in increments of 0.01 (illustrated by Fig. 12a). At each iteration step, for each tow surface grid point, it is checked in which voxel it is located (Fig. 12b). If the voxel

has not yet been assigned to a material, it is attributed to the current tow. If the voxel is already assigned to a material, the inflation process is stopped for the respective point on the surface grid (see Algorithm 1). The inflation ratio is increased simultaneously for all tows up to a maximum value of E = 110%. This accounts for the fact that if a tow is constrained by the contact with another tow in some direction (the inflation is stopped in this direction), the fibers are redistributed in other directions where there is no contact with other tows (the inflation continues for e > 1 in these directions). Note that no complex intersection procedure is applied in this algorithm. A tow surface grid size of less than 0.2 times the voxel size is sufficient to ensure the robustness of the filling process. The procedure is applied to voxel arrays pre-filled with the respective complement, i.e., the warp tows are inflated with the voxels that during the segmentation process were identified as belonging to the matrix or to the weft tows treated as already assigned. The inflation thus stops as soon as a weft tow or a matrix pocket is encountered. Likewise, the weft tows are inflated with the warp and matrix voxels already treated as assigned. This limits the growth of the tows to their actual dimension observed in the tomography images, and the simultaneous growth of all tows leads to a correct identification of the contact zones between parallel tows.

Meshing procedure

The generated voxel mesh is then smoothed using the method developed by Rassineux [START_REF] Rassineux | Robust Conformal Adaptive Meshing of Complex Textile Composites Unit Cells[END_REF]. It is based on a modified Catmull-Clark subdivision technique and converts the voxel mesh into an unstructured tetrahedral mesh, which allows for mesh quality during the voxel deformation. During the smoothing process, a minimum angle between two adjacent element faces of 20 • is imposed, which avoids the presence of flat elements in the final mesh. The smoothing is carried out using the reduced voxel mesh in order to limit the computational costs. When returning to the original scale of the CT-scan, the resulting minimum angle between element faces is 8.2 • . In the studied case, the smoothing combined with the rescaling of the REV increases the volume fraction of the tows by about 4%. This difference increases with the minimum angle. The value was thus chosen as a compromise between mesh quality and respecting the tow shape.

After smoothing, the volume fraction of the matrix in the REV is 16.3%, which is much closer to the real value (see Table 1). The volume fraction of the tows increases accordingly After smoothing, the mesh contains 18.4M nodes and 105.5M elements. Only 0.3M

elements have a quality below 0.1. The quality is measured for a tetrahedron as the ratio of the inscribed sphere of the element to the longest edge times a coefficient chosen such that an equilateral element has a quality of 1. The mesh is then coarsened using the algorithm developed by Rassineux [START_REF] Rassineux | Robust Conformal Adaptive Meshing of Complex Textile Composites Unit Cells[END_REF]. The final mesh is composed of 1.1M nodes and 6.5M elements, i.e., a reduction by a factor of 16.5. The number of elements with a quality between 0.1 and 0.3 is 4313 and there is no element below this quality. The final coarsened mesh is presented in Fig. 14. The low quality elements are distributed over different locations in the mesh and not necessarily found in areas where the geometry is complex, such as in thin matrix layers as could be expected. In order to ensure the convergence of the 3D mesh generator, new points must be added to solve configurations known as Schönhardt polyhedron, especially when the local volume to mesh is non-convex and lower quality elements may be created [START_REF] Schönhardt | Über die Zerlegung von Dreieckspolyedern in Tetraeder[END_REF].

Application to the material with defects

The presented approach is designed to work with local tow disorientations (blue box in Fig. 15a) as induced by the specific manufacturing process of the material. This process also causes tow splitting as highlighted in red in Figs. 15a and15b.

The segmentation procedure described in section 3 does not need to be modified. The path tracking algorithm (section 4) has to be adapted to deal with the split tows. These tows have to be identified by visually inspecting the tomography. In these tows, two distinct If one tow half undulates much more than the other, the tracking algorithm will follow the straighter half for both paths. In the studied volume, two split tows (4 paths in total) require the same manual intervention as presented in section 4 to correct the paths and to obtain a correct representation of the topology with the defects.

The split tows are composed of two strands of different sizes, which is accounted for in the inflation algorithm 1 by multiplication of the inflation ratio r i with the respective size ratio. The strands grow thus at the same speed as the rest of the tows. The coefficient is determined by measuring the strand width on several sections througout the split zone.

The resulting mesh is presented in Fig. 15c and some details of through-thickness slices are superposed to the corresponding CT-images in Fig. 16.

The presented approach deals with the tow orientation defects in the studied material without any problem. In fact, the defect introduces an out-of-plane displacement of less than the thickness of two plies, which does not cause an important modification of the tow curvature. However, tow splitting implies a topological modification of the reinforcement architecture and the tow tracking has to be adapted and some additional manual corrections are needed in order to deal with it. On the other hand, tow splitting is a very specific defect of the manufacturing process of the studied material, which are rarely found in usual textile composites.

Conclusion

In the proposed approach, the process of generating mesoscopic REV models of composites with woven reinforcements has largely been automated. The approach is not optimized for computational cost. The segmentation process takes about 30 minutes and the tracking about 3 hours. Both steps can be implemented in parallel, which would significantly accelerate the processing time. However, the meshing and the coarsening algorithms are implemented in parallel and take about 2 hours for the studied REV. The manual input on the reinforcement topology is limited to a minimum. Only the centers of the tows have to be indicated on the 4 sides of the volume. The main challenge was to process low quality CT-scans with reconstruction artifacts. The ring artifacts are handled by directional filtering. However, the tow texture close to the rotation axis during the CT-scan is too degraded for the proposed method. Manual correction of the tow segmentation is required in these degraded areas of the REV. The work requires less than 20 minutes compared to more than ten hours for a full manual segmentation of the REV.

The tow path tracking also needs some manual correction in some critical areas where parallel tows share a large contact zone. This is the case for 6 out of the 140 tows of the REV. In the case of geometries with fewer contacts, the approach should require even less corrections. The time-consuming manual definition of the reinforcement geometry usually required for geometrical models is avoided. The tow labeling and the meshing are fully automated and do not require any action of the user.

Due to the directional filtering during the segmentation process, the volume fraction of the tows that are not affected by the ring artifacts is slightly increased with respect to the precise manual segmentation. Furthermore, smoothing operations to deal with the low contrast and noise correction reduce the tow cross sections and thus lead to a significant overestimation of the matrix volume fraction. However, the shapes of the tows are correctly determined and superpose well with the CT-images. The smoothing processes used to generate a conformal mesh from the segmented voxel representation compensates to some extent the overestimation of the matrix volume fraction. In mechanical simulations, the quantity of fibers in a tow is known and the local fiber volume fraction has to be adapted to the local tow cross section. In this way, the effect of small variations of the tow volume fractions on the mechanical behavior can be compensated as long as the tow shapes, undulations, and their topology are correctly represented. The fully automated meshing algorithm transforms the voxel mesh into a conformal tetrahedral mesh while avoiding flat elements. High quality elements for accurate mechanical simulations are thus ensured.

The presented approach is able to deal with tow orientation defects as long as the topology of the reinforcement is not significantly impacted, which was not the case with the studied

Figure 1 :

 1 Figure 1: (a) CT-scan of the material, (b) details of in-plane and through-thickness slices, (c) gray level distribution of the whole CT-scan

Figure 2 :

 2 Figure 2: (a) In-plane slice of the CT-scan, (b) orientation map obtained with the structure tensor, (c) coherence, (d) binarized orientation map.

Fig. 3c isFigure 3 :

 3 Fig. 3c is obtained after repeating the hole removal with the same parameters as in the first step.

Figure 4 :

 4 Figure 4: (a) FFT of an in-plane slice of the tomography in the frequency domain, (b) directional filtering mask, (c) slice of the tomography after filtering, (d) result of the segmentation with directional filtering.

Figure 5 :

 5 Figure 5: (a) Through-thickness section of the tomography, (b) result of the segmentation with directional filtering, (c) result after manual correction and cleaning, (d) final segmentation superposed to the tomography image.

Figure 6 :

 6 Figure 6: Final segmented 3D REV.

Figure 7 :

 7 Figure 7: Geometrical variability measured in the CT-scan (a) compared to the variability measured on the segmented volume (b).

Figure 8 :

 8 Figure 8: Diagram of the separation process of parallel tows.

Figure 9 :

 9 Figure 9: Example of implementation of the tow tracking steps with the DTOCS algorithm on a simplified 2D model.

Figure 10 :

 10 Figure 10: Tracked weft and warp tow paths.

Figure 11 :

 11 Figure 11: (a) Rescaled segmented slice of the REV, (b) tracked tow paths with initial average circles, (c) best fitting circles with their centers in red next to the tracked tow paths in blue.

Figure 12 :Algorithm 1 :

 121 Figure 12: Schematization of the inflation process.

to 38 .

 38 0% for the warp and 45.7% for the weft tows. The meshing step doesn't impact significantly the volume fraction inside the REV. The variation of the volume fraction mesured before and after the meshing is less than 1%. The tow shapes are well preserved, because the node displacement during the smoothing process is less than the size of a voxel. Some details of the smoothed REV are shown in Fig.13.

Figure 13 :

 13 Figure 13: (a) Smoothed ply, (b) through-thickness section of the final mesh after smoothing superposed to the CT-scan.

Figure 14 :

 14 Figure 14: Final coarsened mesh of the REV

Figure 15 :

 15 Figure 15: (a),(b) CT-scan of the material with defect highlighted, tow disorientations in blue, tow splitting in red (c) model of the material with defect generated with our approach

Figure 16 :

 16 Figure 16: Detail of the mesh overlapped with the CT-scan

Table 1 :

 1 Confusion matrix for the segmentation approach (SA). The reference used for the comparison is the manual segmentation (ground-truth : GT).

	(%)	Matrix in SA Weft in SA Warp in SA Total Vf in GT
	Matrix in GT	11.2	0.4	0.7	12.3
	Weft in GT	4.1	36.4	3.3	43.8
	Warp in GT	2.5	0.3	41.1	43.9
	Total Vf in SA	17.8	37.1	45.1

material. Likewise, the tow undulations in the studied weaving are not an issue. However, the approach was not tested on composites with 3D woven reinforcements. If there are tows that are oriented locally almost in the through-thickness direction, the inflation process will have to be adapted. Also, the current approach already needs some manual corrections if the tow paths are incorrectly tracked trough parallel tows with large contact areas. In 3D weavings, where due to the high tow waviness some tow parts may be significantly longer than the parts of other tows in the same zone, this may become an issue as the tracking algorithm always chooses the shortest path. Therefore, in the case of 3D weavings with large contact zones between parallel tows, the number of manual corrections of the tow tracking may increase considerably.

The segmentation of split tows requires some adaptations of the algorithm, such as double path tracking for the split tows. The split tows have therefore to be identified by visual inspection of the tomography and taken into account when placing the starting points for the tow tracking on the REV sides. The tow splitting defect is very specific to the manufacturing defect of the studied material.

This application shows that the approach is designed for continuous reinforcements. Tow breaks or ply drops cannot be treated without adaptations. Furthermore, as the approach is based on the analysis of the fibrous texture of the tows, it has to be visible in the CT-images (which is the case in our images taken with a standard laboratory tomograph). However, scans taken with even lower resolutions without any fibrous texture in the tows cannot be treated with the presented approach without significant adjustments.

Declaration of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.