R Gantois
email: renaud.gantois@mines-albi.fr

A Cantarel
email: barthur.cantarel@iut-tarbes.fr

J.-N Félices

N Pirc

F Schmidt

BEM Simulation Of 3D Updated Resin Front For LCM Processes

Keywords: Numerical simulation, Boundary Element Method, Moving Mesh, Liquid Composite Molding, Darcy's law

The use of composite materials in large structures has increased in recent years. Aircraft industry has recently begun to investigate the field of Liquid Composite Molding (LCM) through research programs because of its ability to produce large parts at low cost. The present paper focuses on modeling a 3D radial impregnation through an anisotropic fibrous perform. As a preliminary work, it is assumed an isothermal flow and no hydro-mechanical coupling. Governing equations are Darcy's law and mass conservation. Simulation is performed combining Boundary Element Method (BEM) with a lagrangian moving mesh method. An analytical solution is developed to assess the numerical model.

Introduction

Composite materials have been used in the aircraft industry over the past three decades for their excellent specific mechanical properties. Large structural parts made of resin/fiber composites require specific forming processes. The traditional prepreg autoclave process produces high quality parts, but involves heavy manufacturing cost due to autoclave investment and use of high value added raw materials. As an alternative, recent technologies in Liquid Composite Molding (LCM) have become established recently under the designation of infusion processes. Liquid Resin Infusion (LRI) is one of them. It consists in forcing the liquid resin to circulate transversally into the dry preform using a distribution medium, consecutively to the vacuum application [START_REF] Correia | Analysis of the vacuum infusion moulding process: I. Analytical formulation[END_REF].

Prediction of resin flow during the infusion filling stage requires a full three-dimensional model in order to take into account the lag due to the distribution medium. Previous works revealed the adequacy of Boundary Element Method (BEM) solving the Darcy's flow [START_REF] Um | A study on the Mold Filling Process in Resin Transfer Molding[END_REF][START_REF] Schmidt | Numerical Simulation of Resin Transfer Molding Using Linear Boundary Element Method[END_REF][START_REF] Soukane | Application of the level set method to the simulation of resin transfer molding[END_REF]. It was performed assuming a two-dimensional flow, neglecting the impregnation through the thickness.

The present paper is focused on a three-dimensional numerical modeling of the Newtonian isothermal flow using constant Boundary Elements. Combining Darcy's law with incompressibility equation leads to Laplace equation which is solved using BEM at each step time on an updated lagrangian mesh. Preliminary results are presented, i.e. BEM simulations performed in the case of central injection. Computed results are compared with an analytical solution.

Governing equations

Modeling Assumptions. Impregnation is modeled as a Darcy's flow. Fibrous perform is regarded as a porous medium. Hydro-mechanical couplings are neglected, so that the medium is assumed to be not deformable and of a constant permeability. The flow is considered as isothermal and resin viscosity is assumed to be constant, neglecting chemical reactions effects. Consequently, governing equations reduce to Darcy's law and incompressibity equation.

Darcy's law. Darcy's law [START_REF] Darcy | Les Fontaines Publiques de la Ville de Dijon. V. Dalmont[END_REF] states that the average velocity v ! is proportional to the pressure gradient according to the following equation:

[ ] p K v ∇ - = ! ! µ (1)
where [ ] K is the permeability tensor, p the resin pressure and µ the liquid resin viscosity.

Combining Eq. 1 and incompressibility equation gives:

[ ] 0 . =         ∇ - ∇ p K ! ! µ (2)
For sake of simplicity, we assume in this paper that the principal orthogonal axis of the reinforcement coincides with the reference coordinates system, such as [ ] K reduces to diagonal matrix. Eq. 2 can then be reduced to Laplace equation in the isotropic equivalent domain by stretching the coordinates. In the Cartesian coordinate system, the equivalent coordinates ) , ( 21 e e e x x x = " may be deduced from the physical coordinates ) , ( 21 x x x = ! using the following transformation:

i i e e x K K x i = (3)
where i K is a principal value of [ ] K and e K the equivalent permeability given as

3 3 2 1 K K K K e = .
In the isotropic equivalent domain, Eq. 1 and Eq. 2 reduce simply to:

p K v e e ∇ - = ! ! µ (4) 0 = ∆p (5)
Analytical solution. Following the approach introduced by Adams et al. [START_REF] Adams | Radial penetration of viscous liquid into a planar anisotropic porous medium[END_REF], an analytical solution is used to calculate the evolution of the flow front versus time. Computations are performed in the isotropic equivalent system, using the equivalent permeability The evolution of the flow front is governed by the following differential equation: ε µ

( )         - - =      
- = +       -         (9)
Eq. 9 gives immediately the flow front which is spherical in shape for an isotropic reinforcement. For the orthotropic case, the reciprocal transformation is applied to give the three semi-axis of the flow front which is ellipsoidal in shape.

=umerical Simulation

Boundary Element Method. The calculation domain is delimited by its boundary

∪ q p Γ Γ = Γ
A know pressure p (resp. normal pressure gradient q ) is prescribed on the boundary surfaces p Γ (resp. q Γ ). Eq. 4 is multiplied by the Green function * p . Integrated twice over the calculation domain leads to the well-known Somigiana's equation [START_REF] Brebbia | Boundary Elements : an introductory course[END_REF][START_REF] Pirc | BEM-based cooling optimization for 3D injection molding[END_REF]:

π θ ∫ ∫ Γ Γ i i i c d q p d pq p c (10)
Where i p is the value of the pressure at a point i M located on the boundary (which is supposed to be smooth), * q the pressure gradient associated with * p and θ is the internal angle of the corner in radians. For a three-dimensional domain, * p and * q are given as:

∫ ∑ ∫ ∑ Γ = Γ = Γ = Γ + j j d p q d q p p c 1 j j 1 j j i i (12)
Introducing the 21 matrix H and G , Eq. 12 can be transformed into:

ij 1 1 G q H p 1 j j ij 1 j j ∑ ∑ = = = (13) where ∫ Γ Γ + = j d q H ij * ij 2 1 δ and ∫ Γ Γ = j d p G ij * .
Finally Eq. 13 is reordered to take the form of a linear system as:

F AX = (14)
Where X is the vector of 1 unknowns, A is a 21 matrix obtained by reordering H and G . F is the known vector computed from the boundary conditions, H and G matrix.

Fluid domain updating.

The fluid domain is continuously expanding as it flows, constraining to update the resin flow front. The filling stage is regarded as a succession of quasi-steady states, so that the geometry is redefined at each time step. Unknown pressures p and normal pressure gradients q are calculated using Eq. 14. Nodes on the moving boundary are relocated using a backward Euler integration scheme:

n n v t t x t t x ! ! ! ! ! ) . ( ) ( ) ( ∆ + = ∆ + (15)
Combining Eq. 15 and Eq. 4 yields to:

) ( ) ( ) ( n q K t t x t t x e ! ! ! µ ∆ + = ∆ + (16)
procedure updates the resin flow front, but induces confined mesh degeneration. To overcome this difficulty, an additional adjustment is performed relocating the nodes in contact with the mold.

Algorithm. The implemented algorithm is summarized below:

• Initial computation domain is meshed.

• Stationary boundary conditions p or q are applied on each BEM element.

• Nodes coordinates are transformed according to Eq. 3.

• Normal pressure gradients q are calculated according to Eq. 14.

• Normal pressure gradients are prescribed on each node averaging values from surrounding elements. • Nodes on the moving boundary are relocated according to Eq. 16. Nodes in contact with the upper mold wall are relocated according to a homogenization procedure.

• New normal pressures p and normal pressure gradients q are calculated next steps. The three last steps are repeated until the allotted mold-filling time is reached. • Nodes coordinates are transformed according to the reciprocal transformation.

Validation of numerical simulation

Boundary conditions. The simulation consists in propagating the resin front into an anisotropic reinforcement characterized by its permeability tensor [K] and its porosity ε . A constant pressure 0 p is prescribed at the spherical inlet of radius 0 r . On the flow front, the pressure f p is constant and set by the vacuum pump. Processing and rheological parameters assigned are referenced in Table 1. They keep close to the LRI process conditions but the distribution medium is not taken into account as a first approach.

Fig. 1 shows the boundary conditions prescribed on the boundary of the initial fluid domain. The non penetration condition is assumed assigning 0 = q at the upper mold wall. The initial mesh generation is performed on the equivalent isotropic domain. In order to assess the mesh effect on the numerical method, some simulations have been carried out meshing the physical domain. Deviations are all the more important that the permeability values are different. This is due to the transformation from Eq. 3 which stretches the mesh, degenerating elements before any calculation.

Flow front at t=0.5 s Flow front at t=10 s

Conclusion

The 3D radial front propagation is accurately simulated. The proposed technique is able to give a high accuracy of the front position, making it suitable for front tracking in permeability measurement [START_REF] Nedanov | A method to determine 3D permeability of fibrous reinforcements[END_REF]. It is a great advantage comparing to Volume Of Fluid techniques. However, the proposed method is not directly applicable to complex shapes, because there is no contact algorithm implemented. In order to perform it, current works [START_REF] Gantois | Numerical Simulation Of Resin Transfer Molding Using BEM and Level Set Method[END_REF] focus on fixed grid approach, more suitable to handle contact (in particular merging and dividing fronts). Future works will include hydro-mechanical couplings occurring during Liquid Composite Molding process, such as fibrous perform deformations.

  e K and the spherical injection gate of radius e r 0 . A constant pressure injection e p 0 is applied, so that the macroscopically resin front shape is spherical of radius e f r . In such configuration the resin pressure varies only in the radial direction according to: the expression of the gradient pressure into Darcy's law and introducing the system porosity ε gives the radial component of the front velocity e

  time dependant flow front position can be deduced by integrating the previous equation during the filling stage. This leads to the following solution:

  the point of application of the Dirac delta function and P the point under consideration. Meshing the boundary into 1 constant boundary elements and applying Eq.

Figure 1

 1 Figure 1: Boundary conditions

  Comparison with analytical solution. For symmetry reasons, the simulation is performed on a quarter domain meshed using 540 elements as shown in Figure2. The step time is set at 10 -2 s. CPU time for the complete simulation is around one hour on Intel® Core™2 Duo CPU (2.26 GHz, 1.93 Go of RAM) using Matlab® environment. The computed flow fronts are in fair agreement with the analytical one as shown in Figure3. At 10 s elapsed time, the deviations between the numerical and analytical flow fronts positions do not exceed 1%. Prolongating the simulation until 60 s gives deviations around 3%, due to the gradual mesh degeneration as the iterations go along. Figure4illustrates the flow front evolution during the filling stage, comparing the front shapes at t=0.5 s and t=10s.

Figure 2 :

 2 Figure 2: Initial mesh domain Figure 3: Comparison with analytical solution

Figure 4 :

 4 Figure 4: Flow front computed at t=0.5 s and t=10 s

Table 1 :

 1 Processing parameters

with * * = Γ = Γ +

Acknowledgments

This work was supported by DAHER-Socata and CRCC-Resource and Competencies Center for Composite Materials. The authors are especially grateful to M. Niquet for its assistance during this works.


Sabine Le Roux

Farid Medjedoub

Gilles Dour

Farhad Rezai
email: rezai@mines-albi.fr

Farhad Rézaï-Aria

Image analysis of microscopic crack patterns applied to thermal fatigue heat-checking of high temperature tool steels

Keywords:

come

Image analysis of microscopic crack patterns applied to thermal fatigue heat-checking of high temperature tool steels

Introduction

One of the major damage mechanisms occurring in hot forming tool steels under thermal fatigue solicitations is the formation of a network of interconnected cracks, often named "heat-checking" [START_REF] Jean | An investigation on heat checking of hot work tool steels[END_REF][START_REF] Medjedoub | Endommagement par faïenc ¸age des moules de fonderie sous pression en fatigue thermique: origines, mécanismes et approches[END_REF][START_REF] Norström | Thermal fatigue behaviour of hot work tool steels[END_REF][START_REF] Persson | Simulation and evaluation of thermal fatigue cracking of hot work tool steels[END_REF]. This kind of crack pattern is generally observed when a surface is subjected to the removal of a diffusing heat (in cooling) or liquid (in drying), as for example in thermally quenched ceramic materials [START_REF] Korneta | Topological and geometrical properties of crack patterns produced by thermal shock in ceramics[END_REF], dried mud [START_REF] Velde | Structure of surface cracks in soils and muds[END_REF][START_REF] Vogel | Studies of crack dynamics in clay soil. I: experimental methods, results, and morphological quantification[END_REF], drying of thin layers and films prepared by a sol-gel process [START_REF] Bockmeyer | Crack formation in TiO2 films prepared by sol-gel processing: quantification and characterization[END_REF][START_REF] Shorlin | Development and geometry of isotropic and directional shrinkage-crack patterns[END_REF]. [START_REF] Persson | Simulation and evaluation of thermal fatigue cracking of hot work tool steels[END_REF] have related the heat-checking developing on hot work tool steels to a strong thermal gradient that provides tensile stresses during the cooling phase. In brittle materials, the tensile stresses are relieved by the development of a crack pattern on the surface [START_REF] Korneta | Topological and geometrical properties of crack patterns produced by thermal shock in ceramics[END_REF].

Different techniques are used to analyze surface crack damage. [START_REF] Norström | Thermal fatigue behaviour of hot work tool steels[END_REF] evaluated and ranked the steel resistance to heat-checking by a comparative method based on in-house standard charts. This technique allows a qualitative classification of the steels, but is limitative. Indeed, such method does not give a very precise quantification of the damage grading, and it is only suitable for a specific type of damage at a given scale of observation. Moreover, this evaluation is rather subjective, since the results can vary from an observer to another. More recently, quantitative investigations have been proposed using image analysis software to evaluate the "crack density" [START_REF] Andersons | Advanced fragmentation stage of oxide coating on polymer substrate under biaxial tension[END_REF][START_REF] Bockmeyer | Crack formation in TiO2 films prepared by sol-gel processing: quantification and characterization[END_REF][START_REF] Jean | An investigation on heat checking of hot work tool steels[END_REF][START_REF] Lemoine | Etude de la fatigue thermique superficielle d'un acier inoxydable austénitique[END_REF][START_REF] Maillot | Thermal fatigue crack networks parameters and stability: an experimental study[END_REF][START_REF] Wu | Computer aided evaluation of thermal fatigue cracks on hotwork tool steel[END_REF]. These methods constitute a more precise tool to qualify the damage rate. Different definitions (and unities) are proposed for the crack density, which are not always clearly expressed. It can be stated that the crack density is a global parameter, which solely represents the severity of the crack network, independently of its morphology. The fractal dimension, measured by the "box-counting" method, is sometimes used to quantify the complex geometry of soil cracking patterns [START_REF] Velde | Structure of surface cracks in soils and muds[END_REF] or heterogeneous materials under compressive loading [START_REF] Yan | A quantitative study on the surface crack pattern of concrete with high content of steel fiber[END_REF]. This approach, which leads also to represent the irregularity of the crack pattern by a global coefficient, generally requires different scales of observation of the same pattern, and is furthermore difficult to interpret. Some geometrical aspects, like the fragment area [START_REF] Andersons | Advanced fragmentation stage of oxide coating on polymer substrate under biaxial tension[END_REF], the crack segments length [START_REF] Lemoine | Etude de la fatigue thermique superficielle d'un acier inoxydable austénitique[END_REF][START_REF] Velde | Structure of surface cracks in soils and muds[END_REF], or the orientation of the crack segments [START_REF] Lemoine | Etude de la fatigue thermique superficielle d'un acier inoxydable austénitique[END_REF], have been in addition examined in different works.

Some studies focus specifically with topological properties [START_REF] Korneta | Topological and geometrical properties of crack patterns produced by thermal shock in ceramics[END_REF][START_REF] Shorlin | Development and geometry of isotropic and directional shrinkage-crack patterns[END_REF], and the hierarchical structure of the crack patterns [START_REF] Bohn | Hierarchical crack pattern as formed by successive domain divisions. I: temporal and geometrical hierarchy[END_REF]. These approaches, applicable to any cellular structure, provide interesting information on the structural organization of a crack network. For example, the distribution of the number of sides or neighbors of the The present paper deals with a semi-automatic analysis of the superficial crack pattern from SEM micrographs, using image processing and mathematical morphology tools. This method is applied to characterize the microscopic heat-checking patterns, appearing generally on tool steels under thermal fatigue tests. The oxidation and cracking mechanisms are described. The heatchecking density is measured and some quantitative features are determined on the crack branches (length, orientation), the branching points (density, class junctions, angles) and the cells (size, shape, number of first neighbors). A complete description of the geometrical and topological features of the heat-checking pattern is reported for a unique test condition.

Thermal fatigue experiments

Test specimens

Thermal fatigue tests are performed on X38CrMoV5 (AISI H11) tool steel (see Table 1), quenched and double tempered to achieve a tempered martensitic structure with 47 HRC hardness [START_REF] Jean | An investigation on heat checking of hot work tool steels[END_REF]Medjedoub, 2004). Specimens are of hollow cylindrical shape with a central part of 40 mm length, a 30 mm external diameter, and an internal axial hole of 10 mm diameter, for circulation of cooling water (Fig. 1a). The specimens have a wall thickness of 10 mm in the central gauge area. The external surface is mechanically polished to reach an arithmetic mean roughness (R a ) of 0.02 m.

Thermal fatigue rig

The thermal fatigue rig uses the induction heating of specimens by a CELES-25 kW high frequency generator (100-400 kHz). The induction frequency is about 115 kHz, resulting in a very fast heating of the specimen. During the test, continuous cooling is ensured by cold water circulating through the internal hole with a flow rate of 20 l/min. The external surface is air-cooled by natural convection. The temperature of the specimen is measured by a K-type thermocouple spot welded to the specimen surface (Fig. 1b).

Test conditions

Different thermal cycles, with a maximum temperature (T max ) ranking from 600 to 685 • C and a heating time (ht) ranging from 1.0 to 1.6 s, are applied (Fig. 2). For all the experiments, the cooling time (ct) is adjusted to achieve a minimum temperature (T min ) of 100 • C on the specimen surface (Table 2). In order to investigate the surface heat-checking and cracking evolution, the thermal fatigue experiments were regularly interrupted at different intervals: every 500 or 1000 cycles until 3000 cycles, then every 2000 or 2500 cycles up to 15,000 cycles, and every 10,000 cycles right up to the end of the test [START_REF] Medjedoub | Endommagement par faïenc ¸age des moules de fonderie sous pression en fatigue thermique: origines, mécanismes et approches[END_REF].

Crack pattern analysis procedure

SEM image acquisition

The heat-checking surface morphology is observed using a scanning electron microscope (SEM) with a back-scattering electrons detector. A magnification of 250 times, sufficient to distinguish the micro-cracks, is selected. A cartography of a randomly selected area of 0.825 mm 2 , chosen in the central region of the specimen where the temperature is uniformly distributed ("gauge area"), is systematically reconstructed by assembling a matrix of (3 × 3) image-fields (Fig. 3). To facilitate this reconstruction, an overlapping region representing approximately 15% of the image field is taken between two adjacent images. The resulting cartography is saved in the "Tiff" format, and transferred on the image analysis unit. The size of this digital image is about (1740 × 1210) pixels 2 , leading to a resolution of 0.5 m/pixel. This mapping technique is optimized to provide detailed information on a large representative area, with a sufficiently precise resolution to detect the microcracks.

Image processing

The image analysis software Aphelion [ADCIS] is used to assess and evaluate the heat-checking density and morphology on the cartographies. Two algorithms are employed depending on the image contrast.

When the contrast between the micro-cracks and non-cracked oxide layer is sufficiently high (as in Fig. 4a), the cartography analysis is performed automatically. The treatment consists of several steps including pre-processing techniques, Boolean operations and mathematical morphology transformations [START_REF] Coster | Précis d'analyse d'images[END_REF][START_REF] Russ | The Image Processing Handbook[END_REF][START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]:

1. Median filtering using a (3 × 3) pixels 2 mask, to reduce the noise in the initial gray-scale image (Fig. 4b). 2. Contrast enhancement of "b" by the "top-hat" transformation using a (5 × 5) pixels 2 structuring element, followed by a linear expansion of the gray scale range (Fig. 4c). This allows to extract thin black lines (representing the cracks) placed on an heterogeneous background (oxide surface). 3. Gray double thresholding of image "c" followed by a region growing routine, to extract the darkened micro-cracks (Fig. 4d).

The obtained binary cracks network may contain some discontinuities resulting from an incomplete thresholding of the finest cracks. 4. Complementation of the thresholded image "d" in order to split the resulting heat-checking cells, using a "watershed segmentation" which reconstructs most of the cells boundaries (Fig. 4e). This operation suppresses the small discontinuities in the cracks network. 5. Thickening of the binary image "e" using a "skeleton by zone of influence", followed by a complementation of the resulting image (Fig. 4f).

This algorithm leads to a skeletonized representation (reduced to a single pixel in width) of the fully inter-connected cracks network prior to quantification analysis.

When micro-cracks are not sufficiently opened (specially at the beginning of the heat-checking process), the contrast is often not high enough to allow an automatic detection of the heat-checking pattern (Fig. 5a). In this case, the cartography is carefully drawn on a transparent paper and digitized using a scanner at 256 gray levels and 150 dpi. The resulting gray scale images have a lower resolution (estimated at 1 m/pixel) than original SEM, due to the width of the manual line tracing. Although this acquisition step is time-consuming, it greatly facilitates the image processing. Indeed, the binary thinned crack network can be directly obtained from the scanned cartography by simply applying an automatic thresholding (using the "entropy maximization" algorithm developed by [START_REF] Kapur | A new method for gray-level picture thresholding using the entropy of the histogram[END_REF] followed by a skeletonization of the cracks (Fig. 5). The manual drawing of the cracks induces certainly a subjectivity, as very short edges may be ignored or smoothed. But it remains nevertheless preferable to a direct analysis of the initial cartography which would lead to uncompleted segmentation and produce biased results. The binary skeletonized structure (Figs. 4f and5d) represents in fact a thinned median axis of the cracks, which is a powerful tool containing both geometrical and topological information. The most direct way to describe such a crack network is to calculate the crack density. We define this global parameter, mentioned hereafter as the heat-checking density " hc " (in mm -1 ), as the cumulated crack length "L c " (in mm) divided by the area of the reference surface "A ref " (in mm 2 ):

hc (X) = L c (X) A ref
(1)

Geometrical parameters

The heat-checking network can be regarded as a planar mosaic constituted of quasi-convex polygons (called "cells"). To characterize the morphology of these cells, it is necessary to invert the binary image of the skeletonized crack network by the complementation operation (Fig. 6a). The cells which are partially included in the image field are removed (Fig. 6b), and the remaining cells are labeled (Fig. 6c).

For each cell c i , the following calibrated geometrical parameters are calculated: area A(c i ), Crofton perimeter P(c i ), length (evaluated by the minimum Ferret diameter F min (c i )), and width (evaluated by the maximum Ferret diameter F max (c i )) [START_REF] Coster | Précis d'analyse d'images[END_REF][START_REF] Wang | Image analysis of particles by modified Ferret method-best-fit rectangle[END_REF]. The shape of the cells is characterized by the compactness, defined as:

C(c i ) = 16.A(c i ) P(c i ) 2 (2)
and the elongation, defined as:

E(c i ) = F max (c i ) -F min (c i ) F max (c i ) + F min (c i ) (3)
both of which range between 0 and 1. The compactness value is 1 for a perfect square, and smaller for shapes with irregular boundaries. The elongation equals to 0 in the case of a circle or a square, and approaches 1 for an ellipse that is long and narrow. All these data are calculated on the entire field, and used to obtain frequency distributions describing the cells geometry and statistical parameters (such as the average size and the standard deviation). Because the heat-checking cartographies are taken on a large area, the number of cells included inside the images is very high (between 500 and 4000 cells depending on the test conditions). The classical corrections generally required for planar sampling, like the Miles-Lantuéjoul statistical method used to compensate the edge effects [START_REF] Lantuéjoul | Grain dependence test in a polycrystalline ceramic[END_REF], are therefore not needed.

On the skeletonized crack network, the pixels having at least three neighboring pixels are considered as cracks junctions or crack intersections (called "branching points" or "nodes"). These nodes delimit the cell sides called "crack branches". To analyze the morphological characteristics of crack branches, one need first to remove the nodes from the skeletonized crack network (Fig. 7a-c). This results in individualization of the crack branches. The crack branches intersecting the field border are eliminated and the remaining ones labeled.

The length and the orientation of the crack branches are measured respectively by the maximum Ferret diameter (F max ) and its corresponding angle  [START_REF] Wang | Image analysis of particles by modified Ferret method-best-fit rectangle[END_REF]. The possible range between 0 • and 180 • for the cracks orientation is divided into 36 classes having a width of 5 • . The cumulated crack length is calculated as the centered sum of all branches having an orientation in the range of ±2.5 • around the prescribed direction, and the normalized frequency distribution of cracks orientation P(Â) is presented as polar graphs (called "rose plot").

The crack junction orientations are determined by the following procedure. The nodes of the crack network are extracted one by one on the labeled image of the branching points (Fig. 7d), and dilated twice using a square structuring element (Fig. 7e). The intersection between the dilated nodes and the crack branches of Fig. 7c provides markers in the branches radiating from the considered node. A geodesic dilatation (also called "reconstruction") from these markers reconstitutes the complete radiating branches (Fig. 7f), of whom the orientations are measured using the maximum Ferret angle (defined above) with a precision of 1 • . The angle (ˇ) between each radiating branch is then calculated, and the frequency distribution of the cracks junction orientation is achieved considering 18 classes having a width of 10 • .

The previous procedure allows also to classify the crack junctions with respect to the number of radiating branches, simply by counting the individualized branches in Fig. 7f. For the sake of simplicity, we call "class 3 junction" a node connected to three crack branches, "class 4 junction" a node connected to four crack branches, "class 5 junction" a node connected five crack branches, etc.

Topological parameters

Topological parameters are features which remain invariant when the analyzed object is modified by a geometric transformation (like stretching or distortion). From a topological point of view, the crack network can be characterized by the number of crack branches (N b ), the number of nodes (N n ) whereas the cellular structure (i.e., the inverted crack network), by the number of cells (N c ). The density of crack branches (d b ), nodes (d n ) and cells (d c ) are then expressed in mm -2 , as the ratio of the number of each feature (corrected to take into account the edge effects) to the reference surface A ref .

Another interesting topological feature of the cellular structure is the number of first neighbors n of each cell, and specially the average value ⟨n⟩ and the frequency distribution P(n). Two cells are considered neighbors if they share at least one side. According to the Lantuéjoul's procedure [START_REF] Lantuéjoul | Grain dependence test in a polycrystalline ceramic[END_REF], the determination of n requires four steps. First, the individualized cells are extracted one by one from the cellular mosaic (Fig. 8a). Then, the selected cell is dilated twice using a square structuring element (Fig. 8b). The intersection between the dilated cell and the mosaic structure provides markers both in the analyzed cell and in the first neighboring cells, which are reconstituted by a geodesic reconstruction (Fig. 8c). It is straightforward to remove the central cell and count then the number of remaining cells constituting the first neighbors (Fig. 8d). This algorithm is applied to all the cells of the network, and leads to a label image in which the cells are classified according to their number of first neighbors (Fig. 8e). The n-neighbored cells can be selected by a simple thresholding of the label image, and the average number of neighbors of order 2 (i.e., the average number m(n) of neighbors of the cells that are adjacent to n-neighbored cells) can be determined by applying again the Lantuéjoul's procedure [START_REF] Fortes | Average topological properties of successive neighbours of cells in random networks[END_REF]. These topological numbers are then used in some semi-empirical laws describing the correlation between neighboring cells in planar cellular structures.

According to the classic Aboav-Weaire's law, the average number of neighbors of cells surrounding an n-neighbored cell m(n) is linearly proportional to 1/n [START_REF] Korneta | Topological and geometrical properties of crack patterns produced by thermal shock in ceramics[END_REF][START_REF] Zsoldos | Appearance of collectivity in two-dimensional cellular structures[END_REF], which is usually expressed by the equation:

n • m(n) = (⟨n⟩ -a)n + (⟨n⟩ • a + 2 ) (4)
where "⟨n⟩" is the average number of first neighbors of cells, "a" is a system-constant, and " 2 " is the variance of P(n). This relationship indicates that many-neighbored cells tend to have few-neighbored neighbors, and vice versa [START_REF] Fortes | Average topological properties of successive neighbours of cells in random networks[END_REF][START_REF] Korneta | Topological and geometrical properties of crack patterns produced by thermal shock in ceramics[END_REF]. In Eq. ( 4), "a" is related to the slope of the regression line, ranging for most of cellular structures between 0.5 and 1.5. " 2 " characterizes the rate of disordering in the cellular structure: the more it increases, the more the cells are randomly distributed. Different values are reported for " 2 ": 0.6-1.1 in biological cells, 0.24-2.86 in soap foam, 1.72 in arsenic-selenium glass, 1.07 in cork patterns, 2.75 ± 0.07 in crack patterns produced by thermal shock in ceramics, etc. [START_REF] Zsoldos | Appearance of collectivity in two-dimensional cellular structures[END_REF].

Results and discussion

Thermal fatigue cracking mechanism

Detailed SEM observations have revealed that the heat-checking pattern is composed of thin interconnected cracks of different lengths and orientations, due to the oxide-scale cracking under the thermal stresses (Figs. 3 and9). This is due to cyclic transition thermal gradients, induced by repeated heating and cooling of the surface. The morphology of oxide layer changes with the testing conditions, specially the maximum temperature T max of the thermal cycle (Medjedoub, 2004). For T max up to 600 • C, the oxide layer is made of juxtaposed spangles, and when T max increases beyond 650 • C, the structure of the oxide layer become alveolar and more adherent (Fig. 9). Post-mortem observations of specimen cross-sections revealed the formation of a duplex oxide-scale: a superficial layer rich in iron oxide (constituted mainly of hematite oxide Fe 2 O 3 ), and an internal compact layer rich in Cr 2 O 3 . These two layers present quasi-identical thicknesses (Fig. 10a). As can be seen in Fig. 10b, the micro-cracks occur in the outer-layer of the oxide scale and often stop at the interface between the two layers. As the external oxide layer is more porous, one can suppose that the micro-cracks originate on structural defects like porosities. It leads to a homogeneous crack pattern on the oxide surface (Fig. 3), which conforms to the heat-checking density statistical analysis (see hereafter).

Unfortunately, the crack width, which in fact represents the crack lips opening, is not large enough to be precisely quantified considering the cartography resolution. However, SEM observations seem to show that the crack width tends to grow while cycling, specially at the first stages of the heat-checking cell formation (Fig. 11a). In comparison to our case, the crack opening is much easier to observe in desiccated soils [START_REF] Vogel | Studies of crack dynamics in clay soil. I: experimental methods, results, and morphological quantification[END_REF].

Initiation and evolution of the heat-checking density

A typical evolution of heat-checking pattern and density as a function of the number of cycles is shown in Fig. 11, for the reference thermal fatigue test (100-650 • C). It is observed that there is an incubation period where no crack is detected (it must be remembered that the detection limit does not allow to detect the crack with opening less than 0.5 m). Then, micro-cracks initiate, propagate and coalesce rapidly to form a network of interconnected cracks (Fig. 11a andb). Fig. 11c shows a sigmoidal trend for the evolution of the heat-checking density hc versus number of thermal cycles. As soon as heat-checking formation process is initiated, hc raises abruptly to reach a maximum value, prior to stabilizing to a asymptotic limit (named sat ) beyond 5000 cycles. This limit remains quasi-constant (or saturated) for further thermal fatigue cycles. The number of cycles necessary to achieve sat is named N sat (Fig. 11c). The slight decrease of hc following the peak value is explained by an oxide healing of some very shallow cracks. To assess the homogeneity of the heat-checking in the gauge area of the specimen, five cartographies were taken in different zones at the end of two thermal fatigue tests (T max = 600 and 650 • C). A standard deviation ( ) of about 3.5 was obtained for a hc of 135 mm -1 . The small variation of hc values in the saturated stage (Fig. 11c) is of the same order of magnitude as (2.4-4.4, depending on the testing conditions), proving a well stabilized unmodified heat-checking pattern. Furthermore, we take for our careful observations that indeed the crack network remains quasi the same once sat has been reached. At the most, some cracks tend to open more, but follow the initial path.

According to the definition of [START_REF] Glenny | Thermal and High-strain Fatigue[END_REF], the number of cycles for the heat-checking initiation (N i ) is the arithmetic value between the number of cycle when the heat-checking is first detected and the last previous thermal fatigue test interruption cycle, respectively named N 2 and N 1 (Fig. 11d). Because of the drastic sigmoidal evolution of the heat-checking density, the accuracy of the N i value depends principally on the thermal fatigue interruption intervals and on crack detection limit. As in major experiments, tests were interrupted after every 500 cycles prior to heat-checking cells formation, we achieve a relatively fine estimation of N i .

Morphological and topological features of the heat-checking network

The evolution of the morphological and topological characteristics of the heat-checking cells resulting from the oxide-scale cracking have been investigated during the cycling, for the reference thermal fatigue test (100-650 • C). The size and the morphology of the heat-checking cells become more regular and homogeneous with increasing number of cycles (Fig. 12a). The cell compactness and elongation reach stabilized mean values (respectively 0.6 and 0.4), characterizing polygonal cells with irregular boundaries (Fig. 12 b).

Fig. 13, which presents the rose plot of the individualized crack branches orientation, indicates that the micro-cracks develop in all directions. A preferential cracking along 0 • and 90 • , corresponding to the maximal thermal stresses and strains solicitation directions, is however observed. The heat-checking network remains quasiisotropic along with the number of cycles. It can be noted that this isotropic trend is enhanced with increasing the number of cycles. Indeed, at the beginning of the fragmentation process (i.e., for low number of cycles), crack branches oriented along 0 • and 90 • are less predominant (Fig. 13a). When the heat-checking is stable, we can clearly observe four symmetry axes on the rose plot along 0 • , 45 • , 45 • , 90 • and 135 • directions (Fig. 13c).

All along the fragmentation process, the crack junctions are predominantly trivalent in the heat-checking network (Fig. 14a). The same phenomenon was noted by [START_REF] Korneta | Topological and geometrical properties of crack patterns produced by thermal shock in ceramics[END_REF] on quenched ceramics. [START_REF] Bohn | Hierarchical crack pattern as formed by successive domain divisions. I: temporal and geometrical hierarchy[END_REF] have observed that triplets of cracks are generally formed on a surface defect. Frequency distributions of crack junction angles show that they range mostly between 100 • and 150 • , with a peak at 120 • appearing when the network reaches the stabilized regime (Fig. 14b). This angle corresponds, in fact, to "Y-shaped" junctions, related to an hexagonal type crack network. One can also note that 90 • and 180 • angles are weakly represented, and so-called "T-shaped" junctions are less present. A "T-shaped" junction is generally formed when a pre-existing cell is divided in two. Observing similar frequency distributions (with a strong peak around 120 • ) in desiccated clay soils, [START_REF] Vogel | Studies of crack dynamics in clay soil. I: experimental methods, results, and morphological quantification[END_REF] concluded that the cells remains stable after their early formation. [START_REF] Shorlin | Development and geometry of isotropic and directional shrinkage-crack patterns[END_REF] have reported that in drying layers, 120 • junction angles are formed in the early stage of the crack pattern development, whereas 90 • ones correspond to later crack propagation (due to the local stress field around cracks). The fact that "Y-shaped" junctions are predominantly observed on the oxide layer of the steel investigated here proves again that the heat-checking network is also very rapidly stabilized following its initiation.

From the initiation of the heat-checking network to its stabilization, the cells mosaic are constituted of a majority of cells with 5, 6 or 7 first neighbors (Fig. 15). The average number of first neighbors of cells ⟨n⟩ is very closed to 6, which is the ergodic theoretical value expressed by Euler's theorem [START_REF] Bohn | Hierarchical crack pattern as formed by successive domain divisions. I: temporal and geometrical hierarchy[END_REF][START_REF] Lantuéjoul | Grain dependence test in a polycrystalline ceramic[END_REF]. The frequency distributions P(n) follow a quasi-normal law, except when the heat-checking network is not saturated (i.e., for low number of cycles). We can also note that P(n = 6) > P(n = 5) > P(n = 7) is verified whatever the number of cycles. The same characteristics have been observed in quenched ceramics [START_REF] Korneta | Topological and geometrical properties of crack patterns produced by thermal shock in ceramics[END_REF], whereas crack networks formed in a drying layer are mostly four sided [START_REF] Shorlin | Development and geometry of isotropic and directional shrinkage-crack patterns[END_REF]. 1998). The decreasing evolution of 2 versus the number of cycles reveals also that the heat-checking tends to reach a more ordered state in the stabilized regime (Fig. 16b).

Another topological correlation, known as the Lewis' law, states the fact that the average area of the cells (A m ) increases linearly with their number of first neighbors n [START_REF] Zsoldos | Appearance of collectivity in two-dimensional cellular structures[END_REF]. The heat-checking network follows perfectly this law, specially in the saturated regime (Fig. 17a). The slope of the linear variation decreases to a pseudo-saturation (Fig. 17b), and the correlation coefficient, always higher than 0.96, increases slightly with the number of cycles (Fig. 17a). One can note in Figs. 16a and 17a that when n > 10, the scattering increases, due to a lower number of cells for this population.

This investigation of the topological properties of the heatchecking pattern indicates that the stabilized cellular mosaic is a near hexagonal structure. This kind of organization characterizes many random tessellations observed in natural phenomena, for example crack patterns in ceramics [START_REF] Korneta | Topological and geometrical properties of crack patterns produced by thermal shock in ceramics[END_REF], soap foams, cooling-induced patterns in lava flows, shrinkage crack patterns in mud, etc. [START_REF] Zsoldos | Appearance of collectivity in two-dimensional cellular structures[END_REF].

Conclusions

Heat-checking of hot work tool steel (X38CrMoV5) has been investigated under thermal fatigue experiments. The specimen surface is subjected to cyclic induction heating and natural cooling, under various test conditions. The thermal fatigue cycle is periodically interrupted, and the surface of the specimen is examined by scanning electron microscopy. Detailed observations have revealed that a multi-directional connected crack network is mainly formed on a bi-layer oxide scales. The heat-checking formation and evolution has been investigated by image analysis. Automated or semi-automated procedures have been developed, depending upon the quality of the SEM images, in order to describe the morphological and topological features of the microscopic heat-checking network.

The heat-checking density follows a sigmoidal evolution, rapidly reaching an asymptotic value indicating the crack pattern saturation. The heat-checking cells size and shape homogenizes with increasing the number of cycles. The saturated crack pattern is quasi-isotropic, with two dominant orientations corresponding to the maximal thermo-mechanical loading directions. Most of the crack junctions are trivalent, "Y-shaped", with an angle distribution centered around 120 • . This cellular pattern is closed to a trivalent polygonal structure. According to the Euler's law, the average number of neighbors for the cells tends to be 6. An investigation of the correlations between neighboring cells shows that both Aboav-Weaire's and Lewis' laws are applicable for heat-checking patterns, even in the non-saturated regime. All the characteristics of the microscopic heat-checking pattern (aspect, density, morphological and topological parameters) tend to stabilize after a critical number of thermal cycles.

The method develop here could be adapted to quantify the morphology of any kind of interconnected crack network, at any scale of magnification. For example, the crack patterns appearing in the glaze of ceramics [START_REF] Bohn | Hierarchical crack pattern as formed by successive domain divisions. I: temporal and geometrical hierarchy[END_REF], in thin films prepared by sol-gel process [START_REF] Bockmeyer | Crack formation in TiO2 films prepared by sol-gel processing: quantification and characterization[END_REF], and in oxide coatings on polymer substrate under biaxial tension [START_REF] Andersons | Advanced fragmentation stage of oxide coating on polymer substrate under biaxial tension[END_REF] present a quite similar aspect and development to those observed on our specimens.

Fig. 1 .

 1 Fig. 1. Thermal fatigue specimen (a) and rig (b).

Fig. 2 .

 2 Fig. 2. Temperature-time profiles measured in the central area of the specimen for the (100 • C -Tmax) thermal fatigue tests (Tmax = 600, 650 and 685 • C).

Fig. 3 .

 3 Fig. 3. Typical cartography (1.1 mm × 0.75 mm) obtained by assembling 9 SEM images of the specimen surface observed after a thermal fatigue test (BSE detector, 250×).

Fig. 4 .

 4 Fig. 4. Image processing procedure applied to a well-contrasted cartography (only a detailed part of images is represented for a better visualization): (a) initial SEM gray scale image; (b) median filtering of the image a; (c) contrast enhancement of the image b; (d) double thresholding of the image c; (e) segmentation of the binary image of heat-checking cells; (f) skeletonized cracks network.

Fig. 5 .

 5 Fig. 5. Image processing procedure applied to a depicted cartography (only a detailed part of images is represented for a better visualization): (a) initial SEM gray scale image; (b) depicted cartography; (c) automatic thresholding of image b; (d) skeletonization of image c.

Fig. 6 .

 6 Fig. 6. Individual analysis of the heat-checking cells: (a) binary image of cells (complementary image of Fig. 4f or Fig. 5d); (b) elimination of all the cells interrupted by the image border; (c) individualization of the cells by the labeling procedure.

Fig. 7 .

 7 Fig. 7. Detection of the nodes and crack branches (only a detailed part of images is represented for a better visualization): (a) skeletonized crack network; (b) detection of branching points (in gray) on the skeletonized structure; (c) crack branches resulting from the removing of branching points on the crack network; (d) selection of one node; (e) selected node dilated twice; (f) reconstructed radiating crack branches from the selected node; (g) determination of ˇangles between radiating branches from the selected node.

Fig. 8 .

 8 Fig. 8. Procedure to determine the number of first neighbors of the heat-checking cells: (a) selection of one cell in the cellular mosaic; (b) dilatation of the selected cell; (c) neighboring cells reconstruction; (d) labeling and counting of the first neighbors; (e) image of the cellular mosaic in which the cells are classified with a gray level related to their number of first neighbors.

Fig. 9 .

 9 Fig. 9. SEM micrographs of the superficial oxide layer, showing various morphologies depending on the maximal temperature of the thermal fatigue test: Tmax = 600 • C (a); Tmax = 650 • C (b); Tmax = 685 • C (c).

Fig. 10 .

 10 Fig. 10. SEM cross-section observations of the specimen after the reference test (100-650 • C cycle), showing the oxide layer structure (a), and micro-cracks in the oxide outer-layer (b).

Fig. 11 .Fig. 12 .Fig. 13 .

 111213 Fig. 11. Evolution of the heat-checking during the reference thermal fatigue test (100-650 • C cycle): heat-checking aspect on SEM micrographs (a) and cartographies (b) at 1500, 5000 and 15,000 cycles; (c) evolution of the heat-checking density hc versus number of cycles; (d) method to estimate the number of cycles (N i ) for heat-checking initiation.

Fig. 14 .

 14 Fig. 14. Frequency distributions of the crack junction classes P(j) (a), and angles P(ˇ) (b).

Fig. 15 .

 15 Fig. 15. Frequency distributions P(n) of the number of first neighbors n of the heat-checking cells (a), and corresponding classified cellular mosaic images (b), for increasing number of cycles N (1000, 5000 and 20,000 cycles).

Fig. 16 .Fig. 17 .

 1617 Fig. 16. Application of the Aboav-Weaire's law to the heat-checking pattern at different number of cycles N: (a) linear correlation nm(n) versus n (m(n) is the average number of neighbors of cells surrounding a n-neighbored cell); (b) evolution of 2 versus N.

  proposed by[START_REF] Lauschmann | Fractography of crack networks in thin layers[END_REF] to analyze the fractography of crack networks in thin layers. It consists in detecting and classifying the crack elements in finishing, crossing or passing branches in knots of the network. This method can be used to model the formation of a crack pattern in relation with fracture mechanics.

	Nomenclature
	a	system-constant in Aboav-Weaire's law
	A	area of the heat-checking cells
	A m	average area of the heat-checking cells
	A ref	reference area of the digital image
	c i	individual cell of the cellular mosaic
	ct	cooling time of the thermal cycle
	C	compactness factor
	d b	number of crack branches per unit area
	d c	number of cells per unit area
	d n	number of branching nodes per unit area
	E	elongation ratio
	F min	minimum Ferret diameter
	F max	maximum Ferret diameter
	ht	heating time of the thermal cycle
	j	crack junctions class
	L c	total cracks length
	m(n)	average number of first neighbors of n-neighbored
		cells
	n	number of first neighbors of a heat-checking cell
	⟨n⟩ nsd	average number of first neighbors of the cells variance of the number of first neighbors of the cells
	N	number of cycles
	N c	number of heat-checking cells
	N b	number of crack branches
	N i	number of cycles necessary for heat-checking initi-
		ation
	N n	number of nodes
	N sat	number of cycles necessary to achieve sat
	P	Crofton perimeter of the heat-checking cells
	P(A)	frequency distribution of the area of the heat-
		checking cells
	P(j)	frequency distribution of the crack junction classes
	P(n)	frequency distribution of the number of first neigh-
		bors of the cells
	P(ˇ)	frequency distribution of angles between the radi-
		ating branches from a node
	P(Â)	frequency distribution of the crack branch orienta-
		tions
	R a	arithmetic mean roughness
	t	time
	T	temperature
	T min	minimum temperature of the thermal cycle
	T max	maximum temperature of the thermal cycle
	ˇangle between the radiating branches from a crack
		junction
	2	variance of the frequency distribution P(n)
	hc	heat-checking density
	sat	saturated heat-checking density
	Â	orientation of crack branches
	polygonal fragments or crack junction angles can help to under-
	stand the formation mechanism of the crack network, as long as
	the experimental observations are carried out appropriately. Vogel
	et al. (2005) exploit Minkowski functions to investigate the geom-
	etry of the crack patterns formed during the desiccation of clay
	soils. Based on mathematical morphology tools and especially the

distance function, this generic method allows to describe the morphological and topological features of any kind of crack pattern, even not fully-connected ones. But it requires that the area of the crack can be precisely evaluated, and therefore, images with a sufficiently high magnification or resolution are necessary. A complex method, based on automatic edge detection and recognition algorithms, is

Table 1

 1 Chemical composition of X38CrMoV5 (AISI H11) steel.

	Element	C	Cr	Mo	V	Si	Mn	S	P	Fe
	Weight%	0.396	5.05	1.25	0.47	0.94	0.36	<0.003	0.009	bal.

Table 2

 2 Thermal cycles used in the thermal fatigue experiments (the test with Tmax = 650 • C is considered as the reference test).

	Minimum	Maximum	Heating time,	Cooling
	temperature,	temperature,	ht (s)	time, ct (s)
	T min ( • C)	Tmax ( • C)
		600	1.0	15.0
	100	650	1.2	17.0
		685	1.6	22.0

= 1.50 ± 0.23) correspond to the case of a trivalent polygonal system[START_REF] Zsoldos | Appearance of collectivity in two-dimensional cellular structures[END_REF]. The rather low 2 value indicates that the rate of disordering is lower than in crack patterns produced by thermal shock in ceramics, where 2 = 2.75 ± 0.07(Korneta et al.,

Acknowledgments

The authors would like to thank the CTIF (Centre Technique des Industries de Fonderie), in particular Dr. P. Hairy, for supporting a part of this activity. The French National Program on Forging (ACRII and Simulforge) and CETIM are also acknowledged for the financial support.