M Georges

Dargaud

Prof Pascal Forquin

P Forquin
email: pascal.forquin@3sr-grenoble.fr

M Blasone
email: maria.blasone@3sr-grenoble.fr

D Georges
email: david.georges@3sr-grenoble.fr

M Dargaud

Continuous and discrete methods based on X-ray computed-tomography to model the fragmentation process in brittle solids over a wide range of

Keywords: Dynamic fragmentation, A: Voids and inclusions, B: Brittle materials, C: Probability and statistics, C: X-ray Computed Tomography

Introduction

Achieving a good understanding and modelling of the tensile behaviour of brittle materials until failure remains a major challenge for their use in many application areas [START_REF] Weerheijm | Understanding the tensile strength of concrete[END_REF][START_REF] Forquin | A probabilistic damage model of the dynamic fragmentation process in brittle materials[END_REF][START_REF] Forquin | Brittle materials at high-loading rates: an open area of research[END_REF]. The main reason for this is the fact that their dynamic tensile strength and fracturing energy are strongly sensitive to the applied loading-rate. For the past two decades, such behaviour have been extensively studied for several types of brittle materials (ceramics, rocks, concretes and polycrystalline ice) by means of spalling tests conducted using Hopkinson bar apparatus [START_REF] Schuler | Spall experiments for the measurement of the tensile strength and fracture energy at high strain rates[END_REF][START_REF] Weerheijm | Tensile failure of concrete at high loading rates: new test data on strength and fracture energy from instrumented spalling tests[END_REF]Erzar andForquin, 2010, Saletti et al, 2019) and plate-impact experiments [START_REF] Murray | The spall strength of alumina ceramics[END_REF][START_REF] Antoun | Spall fracture[END_REF][START_REF] Grote | Experimental characterization of the dynamic failure behavior of mortar under impact loading[END_REF][START_REF] Paris | The spall strength of silicon carbide and boron carbide ceramics processed by spark plasma sintering[END_REF][START_REF] Dandekar | Spall Strength of Silicon Carbide Under Normal and Simultaneous Compression-Shear Shock Wave Loading[END_REF] or high pulsed power current facilities (Zinszner et al, 2015a) for higher strain-rates. Over a large range of strain-rates, two regimes can be distinguished as shown on Figure 1a and1b. A small strain-rate sensitivity of the tensile strength to strain-rate is noted in the quasi-static regime or at low strain-rates, whereas a pronounced increase of strength with strain-rate is observed above a strain-rate threshold value [START_REF] Erzar | Analysis and modelling of the cohesion strength of concrete at high strain-rates[END_REF][START_REF] Weerheijm | Chapter 7: Response mechanisms of concrete under impulsive tensile loading[END_REF]. The postulate of two regimes interrogates on whether this transition strain-rate is an intrinsic property of the material or depends on size and shape of the considered sample. Moreover, since [START_REF] Reinhardt | Concrete under impact loading, tensile strength and bond[END_REF] many authors have used to express their experimental results as an evolution of the DIF (Dynamic Increase Factor), which corresponds to the dynamic strength to static strength ratio, according to the applied loading-rates. Again, it must be found out if this representation (i.e. DIF versus strain-rate) can be considered as an intrinsic material property or not.

Figure 1. Influence of strain-rate on the tensile strength of brittle solids. (a) Case of dry and wet microconcrete and common concrete [START_REF] Erzar | Analysis and modelling of the cohesion strength of concrete at high strain-rates[END_REF][START_REF] Weerheijm | Chapter 7: Response mechanisms of concrete under impulsive tensile loading[END_REF], (b) Case of two SiC ceramics [START_REF] Zinszner | Strain-rate sensitivity of the tensile strength of two silicon carbides: experimental evidence and micromechanical modelling[END_REF], (c) Post-mortem analysis of dry common concrete after a spalling test at 120 s -1 [START_REF] Erzar | Experiments and mesoscopic modelling of dynamic testing of concrete[END_REF], (d) Microscopic view of the SPS-S damaged specimen recovered after a test performed at 2,000 s -1 [START_REF] Zinszner | Strain-rate sensitivity of the tensile strength of two silicon carbides: experimental evidence and micromechanical modelling[END_REF].

As first answer, it is well established that, under quasi-static regime, the mean failure strength of most brittle materials is sensitive to the size of the tested sample, i.e. the larger the sample size the lower the mean failure strength, and to the stress heterogeneities [START_REF] Weibull | A Statistical Distribution Function of Wide Applicability[END_REF][START_REF] Freudenthal | Statistical Approach to Brittle Fracture[END_REF]Hild et al, 1992). This size effect and the scatter of failure stress in tension can be related to flaws population (Hild and Marquis, 1992;[START_REF] Alava | Statistical Models of Fracture[END_REF]. On the other hand, at high strain-rates (respectively above 100 s -1 and 1,000 s -1 for concretes and ceramics), the size effect is supposed to subside. In such high-rate conditions, as concretes and ceramics have small tensile failure strains, their loading times to failure are limited to a few µs and a few tenths of µs, respectively [START_REF] Hild | On the Probabilistic-deterministic Transition Involved in a Fragmentation Process of Brittle Materials[END_REF][START_REF] Forquin | A probabilistic damage model of the dynamic fragmentation process in brittle materials[END_REF]. Moreover, given the limited crack speed in these materials (about few km/s according to [START_REF] Kanninen | Advanced fracture mechanics[END_REF][START_REF] Freund | Dynamic fracture mechanics[END_REF])), such small loading times imply that damage and failure processes should be driven by the simultaneously propagation of a large number of cracks in a considered centimetre volume. Therefore, in relatively large volume loaded at high strain-rates, what happens in a given location only has an influence on its neighbourhood, which could explain why the size effect is expected to vanish in dynamic conditions.

The way the strain-rate transition between both regimes is expected to evolve, according to tested sample volume and the strain-rate sensitivity of the considered material, is illustrated on Figure 2. Assuming a quasi-static tensile strength σQS and a dynamic strength σDyn, provided by Equations ( 1) or (2), the DIF and transition strain-rate appear to be dependent of the size of the considered volume Z, which means that, in a such framework, these concepts of DIF and transition strain-rate cannot be considered as intrinsic properties of the brittle solid, (

) 2
where / , / are the exponents traducing the strain-rate sensitivity in "quasi-static" and "dynamic" regimes and (! , ! ) are the reference strain-rates defined for each regime. In addition, it can be observed from Equations ( 1) and (2) that the lower the strain-rate sensitivity in the dynamic regime / , the higher the sensitivity of the transition strain-rate to the quasi-static strength and the volume size Z. However, the study of brittle materials tensile strength over a large range of strain-rates and considering different sample sizes remains a difficult task due to the many measurement bias inherent to dynamic testing. For this reason, the size effect at high strain-rates stays a topic hardly explored.

Figure 2. Influence of the volume zone on the transition strain-rate. Cases with (left) and without (right) strain-rate sensitivity in the quasi-static regime.

In parallel to the experimental characterization of the macroscopic tensile strength provided via different spalling techniques, extended efforts have been developed to study the growth of damage and final fracturing patterns in various types of brittle materials at high-strain-rates. Among the numerous existing works, one may mention the use of the edge-on impact testing technique applied to ceramics [START_REF] Strassburger | Visualization of impact damage in ceramics using the edge-on impact technique[END_REF], concretes, mortars and rocks [START_REF] Forquin | Dynamic Fragmentation of an Ultra-High Strength Concrete during Edge-On Impact Tests[END_REF][START_REF] Grange | On the dynamic fragmentation of two limestones using edge-on impact tests[END_REF]Erzar andForquin, 2011, Saadati et al, 2014). In 'open configuration', an ultra-high speed camera is used to visualize the growth of damage 'in real time' [START_REF] Riou | Visualization of the damage evolution in impacted silicon carbide ceramics[END_REF]. In the sarcophagus configuration, a metallic or polymeric casing is used to hold the fragments in place to be able to examine the fragmentation patterns via post-mortem observations [START_REF] Forquin | Effect of aluminum reinforcement on the dynamic fragmentation of SiC ceramics[END_REF]Zinszner et al, 2015b) or determine the fragments size distribution from X-ray tomography analysis [START_REF] Forquin | Application of micro-tomography and image analysis to the quantification of fragmentation in ceramics after impact loading[END_REF]. More recently, transparent sarcophagi were used to perform both analyses (ultra-high speed imaging and post-mortem analysis) in one single test [START_REF] Forquin | Experimental Study of the Dynamic Fragmentation in Transparent Ceramic Subjected to Projectile Impact[END_REF]. Ultra-high speed photography and post-mortem analysis were also applied to spalling tests performed using Hopkinson pressure bar apparatus (cf. example in Figure 1c). This type of analysis is much more difficult for plate-impact experiments, for which the spalled sample is more difficult to recover. However, a post-mortem analysis was carried out in spalling tests performed with the pulse-power technology (cf. example in Figure 1d). From these numerous experimental studies, it can be concluded that the cracking density increases and the mean fragments size decreases as the applied strain-rate increases. Moreover, several studies have demonstrated the major role played by the microstructure of brittle materials, mainly the flaws population, on their fragmentation properties [START_REF] Hild | On the Probabilistic-deterministic Transition Involved in a Fragmentation Process of Brittle Materials[END_REF][START_REF] Zhou | A cohesive-model based fragmentation analysis: effects of strain rate and initial defects distribution[END_REF]. For instance, the fragmentation process induced in a porous SiC ceramic and in the same ceramic after being filled with an aluminium alloy though its open porosity was investigated in [START_REF] Forquin | Effect of aluminum reinforcement on the dynamic fragmentation of SiC ceramics[END_REF]. It was concluded that the difference of cracking pattern was due to the effect of the metal infiltration on the inception stresses of cracks. In (Forquin et al, 2018), the fragmentation properties of four SiC grades were compared by means of edge-on impact, normal impact tests and fragments size

Log σ Log ! ! "#$ %&"&' ! "#$ %&"&' 4 4 1 6 ! ! 1 7 / 1 ! Log σ Log ! ! "#$ %&"&' ! "#$ %&"&' 4 4 6 ! ! 7 ! ! 6 ! !
distribution analysis. The difference of cracking patterns was explained by the differences of Weibull parameters related to the flaws population of each microstructure.

For the last three decades, several analytical and numerical models have been proposed to describe and simulate the fragmentation process in brittle solids. Among them, energetic approaches provide analytical predictions of the fragments size distribution generated from dynamic loadings [START_REF] Grady | Geometric statistics and dynamic fragmentation[END_REF], from which the spall strength strain-rate sensitivity of brittle solids can be deduced [START_REF] Grady | The spall strength of condensed matter[END_REF]. This approach was coupled to a statistical distribution of fracture in [START_REF] Grady | A statistics and energy based theory of dynamic fragmentation[END_REF] for predicting the dynamic fragmentation of U6N rings. Based on an irreversible cohesive law, a cohesive model was proposed by [START_REF] Zhou | A cohesive-model based fragmentation analysis: effects of strain rate and initial defects distribution[END_REF] to predict the one-dimensional fragmentation process taking into account for elastic wave propagation and equally-spaced or randomly-distributed point defect distributions. This approach allowed predicting the fragment size distribution in a circular ring that is dynamically expanded [START_REF] Zhou | Analysis of the brittle fragmentation of an expanding ring[END_REF]. The dynamic fracture of brittle materials was numerically investigated considering element interfaces, therefore including parameters related to interfacial strength, fracture energy and crack opening. These approaches allowed the prediction of fragmentation processes at a meso-scale [START_REF] Camacho | Computational modeling of impact damage in brittle materials[END_REF]Ortiz 1996, Espinosa et al, 1998;[START_REF] Maiti | Mesoscale analysis of dynamic fragmentation of ceramics under tension[END_REF] or taking into account for a heterogeneous distribution of defects [START_REF] Levy | Dynamic fragmentation of ceramics, signature of defects and scaling of fragment sizes[END_REF]. However, such approaches present the main limitation of having prohibitive calculation time, as they simulate a fragmentation process that implies up to few tens of thousands of individual cracks generated in the loaded sample. In addition, parameters related to cracks inception (linked to the flaws population) and crack propagation (crack speed) can be difficult to implement.

In parallel to discrete approaches, continuum damage models have been continuously developed to overcome these difficulties and better capture the physics of multiple-fragmentation induced in brittle materials due to the nucleation, growth and coalescence of multiple cracks in mode I [START_REF] Ravi-Chandar | Dynamic fracture of nominally brittle materials[END_REF] or in the case of sliding cracks uniformly distributed in a brittle material subjected to biaxial compressive loading [START_REF] Ravichandran | A Micromechanical model for high strain rate bahavior of ceramics[END_REF]. Based on a continuum approach, [START_REF] Lu | Modelling modelling of dynamic behaviour of concrete materials under blast loading[END_REF] presented an isotropic damage model built on the concept of micro-crack nucleation, growth and coalescence. Even if this model is meant to describe micromechanical behaviours, it requires parameters identified from dynamic testing to properly describe the strainrate sensitivity of concrete. However, the role of microstructural parameters remains insufficiently addressed in the literature. Whereas mesoscopic models constitute a possible way to account for microstructural properties, they present a number of drawbacks such as calculation costs, difficulties in parameters identification and validation procedure. An alternative is the development and implementation of so-called micromechanics-based models that rely on a description of the physical phenomenon involved in the macroscopic mechanical response of the concerned solid, as it is the case in the Denoual, Forquin and Hild model. This anisotropic damage model is based on the concepts of obscuration probability [START_REF] Denoual | A damage model for the dynamic fragmentation of brittle solids[END_REF] and local weakest-link hypothesis [START_REF] Forquin | A probabilistic damage model of the dynamic fragmentation process in brittle materials[END_REF]. Microstructural parameters are introduced through an explicit law of critical defects and through parameters related to crack propagation. This model was extensively used in the last two decades to predict the fragmentation properties of ceramics [START_REF] Denoual | A damage model for the dynamic fragmentation of brittle solids[END_REF]2002;[START_REF] Forquin | Effect of aluminum reinforcement on the dynamic fragmentation of SiC ceramics[END_REF]2018), concretes [START_REF] Forquin | Dynamic Fragmentation of an Ultra-High Strength Concrete during Edge-On Impact Tests[END_REF][START_REF] Erzar | Experiments and mesoscopic modelling of dynamic testing of concrete[END_REF] and rocks [START_REF] Grange | On the dynamic fragmentation of two limestones using edge-on impact tests[END_REF][START_REF] Saadati | A numerical study of the influence of pre-existing cracks on granite rock fragmentation at percussive drilling[END_REF]. In the present work, the "standard DFH model" is presented and applied to predict the fragmentation properties of three different brittle materials. However, these model predictions rely on the Weibull distribution of flaws deduced from a series of quasi-static bending tests. To overcome the limitations that this implies, the density of critical defects of each three material is identified from X-ray tomographic analysis. These data are used as an input of a modelling based on a continuous and a discrete method. The new possibilities offered by both methods, to suitably predict the dynamic strength and final crack density of brittle materials according to the applied strain-rate, are then discussed.

Predictions and limitations of the DFH (Denoual-Forquin-Hild) model

2.1 The obscuration probability concept [START_REF] Denoual | A damage model for the dynamic fragmentation of brittle solids[END_REF][START_REF] Forquin | A probabilistic damage model of the dynamic fragmentation process in brittle materials[END_REF] The fragmentation process induced in brittle materials subjected to low or high strain-rates tensile loading corresponds to the initiation and propagation of a unique, a small or a large number of cracks from volume or surface defects that originally exist in a structure Ω of size Z, with Z being a length, a surface area or a volume. The material flaws are assumed to be randomly distributed in space and activated at random levels of stress. The activation stress of each defect can be expressed as function of its size, as considered in (Jayatilaka and Trustrum, 1977). The number of flaws in Z for which their activation stress is exceeded, named critical defects, can be represented through a density function corresponding to the number of critical defects for a given applied stress level λt (σ). A schematic of a dynamic fragmentation process assuming that the stress level varies linearly or monotonously with time is presented in Figure 3. The defects are distributed along the horizontal axis (spatial scale) and the vertical axis represents the time-scale or stress-level scale. As the stress level raises (σ ° > 0), the number of activated defects increases. However, the first cracks that propagate from the first critical defects lead to the relaxation of tensile stresses in their neighbourhood (σ ° < 0). Such local stress release prevents the activation of any other critical defect in an "obscured domain" of size Zo(T -t) centred on the crack, where T -t corresponds to the time interval between the crack inception t and current T times. The growth of obscuration domains (zones in which σ ° < 0) from activated defects is represented by triangles (obscuration hyper-cones) in Figure 3, considering a constant speed for crack propagation. The fragmentation process ends when the whole domain is obscured.

The fragmentation process can be described by a probabilistic approach, by considering the condition for a point M at a time T not to be obscured. This condition means that none triggered defect is present in the inverted hyper-cones named "Horizon of (M, T )", the top of which is the point M at time T. According to the concept of "local weakest link hypothesis" introduced in [START_REF] Forquin | A probabilistic damage model of the dynamic fragmentation process in brittle materials[END_REF], the probability of non-obscuration Pno is equal to the product of the elementary probabilities of noinception of new crack P i ∉ in each elementary space-time zones of Ω, belonging to the horizon of this point M at time T:

[ ] i T M of horizon no P T M P ∉ Ω ∈ ∏ = ) , ( ) , ( (3)
where P i ∉ is expressed as a function of the density of critical defects λt (x, t). If λt is considered as a continuous function between t and t+dt, P i ∉ can be written as:

      ∂ ∂ - = ∉ dZdt t t x t x P t i ) , ( exp ) , ( λ . ( 4
)
Otherwise, if the density λt (x, t) is discontinuous between ti and ti+dt ( λt(x, ti+dt) -λt(x, ti) = [λt]x,ti ), the elementary probability of no-inception of new crack in the space-time area dZ x dt located at (x, ti) is expressed as:

[ ] ( )

dZ t x P i t x t i i ) , ( exp ) , ( λ - = ∉ .
(5)

Finally, by combining Equations (3), ( 4) and ( 5), the non-obscuration probability at (M, T) takes the following general expression: Given that new cracks are only initiated in the not-obscured part of the domain, the increment of new cracks can be assumed as proportional to the increment of critical defects weighted by the fraction of non-obscured zones. By considering that a fraction of non-obscured zones is equivalent to the probability of non-obscuration, the increment and derivate of the crack density can be expressed as: ( , ) ( , ) ( , )

[ ] [ ] [ ]         - ∂ ∂ - = ∫∫ ∑ ∫ Ω ∈ ∈ Ω ∈ ∈ ) ,
cracks t no M t P M t M t t t λ λ ∂ ∂ = ∂ ∂ , [ ] [ ] i i t M t i no t M cracks t M P , , ) , ( λ λ = (7)
If the initial density of cracks at time (T = 0) is assumed to be equal to zero, the crack density can be obtained by integrating its derivate or summing its increments:

[ ] ∑ ∫ +       ∂ ∂ = i t M cracks T cracks cracks i t d t M t T M , 0 ) , ( ) , ( λ λ λ (8)
In addition, in the situation where the horizon of point (M, T ) does not strongly interact with the boundaries of the domain [START_REF] Forquin | A probabilistic damage model of the dynamic fragmentation process in brittle materials[END_REF], a multiple fragmentation process takes place. In such case, the probability of obscuration can be used as a damage variable and the macroscopic stress Σ can be calculated as a function of the microscopic stress σ, such as:

Σ(Τ ) = (1 -Po(T ))σ (Τ ) (9)
Therefore, the ultimate (maximum) macroscopic strength can be deduced:

) ( ) ( max u no no t u t t Σ = Σ = Σ with 0 ) ( = Σ u no t dt d (10)
Finally, the knowledge of the density of critical defects appears to be essential to predict the macroscopic strength and cracking density in brittle solids, when subject to high strain-rate tensile loadings.

Influence of the domain volume-size on the type of fragmentation process.

The way the size of the considered domain Ω influences the nature of the fragmentation (single or multiple) is illustrated in Figure 4. If the domain Z is large compared to the mean size of the obscuration zones, the horizon of a points M in Ω at a time T (M, T ) is far from interacting with the domain boundaries. Therefore, other defaults are likely to be triggered, thus leading to a multiple fragmentation process (Figure 4a). Contrary, if the domain size Z is small compared to the size of obscuration cone, the horizon of (M, T ) quickly interacts with the domain boundaries, thus preventing the triggering of any other defects. In such case, a single fragmentation process occurs (Figure 4b).

A closed-form solution of the obscuration probability was proposed in [START_REF] Forquin | A probabilistic damage model of the dynamic fragmentation process in brittle materials[END_REF] by splitting the horizon of (M, T ) above and below the time tz, which corresponds to the intersection of the "horizon inverted cone" with the domain boundaries:

        ×         = ∉ ≤ ∉ > ∏ ∏ i t t Horizon i t t Horizon no P P T M P Z Z ) ( ) ( ) , ( (11)
In the specific case of a uniform stress field, the non-obscuration probability is described as:

λ λ = =     = - - × -             ∫ ∫ 14442444 3 144444244444 3 (12)
At low strain-rates, the time tz approaches the current time T and a single fragmentation process occurs. In such conditions, the obscuration probability (Po = 1-Pno) converges towards the failure probability as proposed in the Weibull model (1939;[START_REF] Weibull | A Statistical Distribution Function of Wide Applicability[END_REF]:

( ) ( ) 1 ( ) 1 exp ( ) failure no eff t P T P T Z T λ = - = - - (13)
At high strain-rates, the time tz is equal to zero or small compared to T and a multiple fragmentation process develops. Therefore, the obscuration probability (Po = 1-Pno) tends to the damage evolution law previously proposed in [START_REF] Denoual | A damage model for the dynamic fragmentation of brittle solids[END_REF].

0 ( ) 1 ( ) 1 exp ( ) T t no o d t D P T Z T t dt dt λ   = - = - - -     ∫ (14)
The analytical expressions of the horizon size Zo(T -t) and the density of critical defects λt (σ) are given in the next subsection. In order to be able to express the obscuration probability of Equation 14, two parameters need to be considered: the density of critical defects and the size of the horizon between time T and t (Fig. 3).

In the present subsection, the density of critical flaws is considered to be a power law of the positive (tension) principal stress (Weibull distribution of critical defects):

8 9 = 8 〈 ; 〉 = , ( 15
)
where the constants m and 8 -1/? are respectively the Weibull modulus and scale parameter. The Weibull modulus gives a direct indication on the scatter in failure stresses. The higher the Weibull modulus, the lower the scatter of failure stresses with respect to the average value. The hypothesis are made for each zone of stress release (obscuration zone) to be centred on the point of crack initiation and for all the obscuration zones to be growing in a self-similar way, with a diameter proportional to the crack size. Under these two assumptions, the size of the obscuration zone at a time T, corresponding to a single crack created at time t, can be expressed as:

[ ] ( )
( ) n o Z T t S k C T t - = -     , ( 16
)
where S is a shape parameter of the obscuration volume equal to 3.74 [START_REF] Denoual | Dynamic fragmentation of brittle solids: a multi-scale model[END_REF],

C
is the one-dimensional wave speed (@ = AB/C) and k is a dimensionless parameter related to the crack propagation velocity. Based upon the concept of conservation of energy, it was demonstrated that when the crack length becomes significantly larger than the initial crack size, the cracking velocity tends to a limit close to 0.38C [START_REF] Broek | Elementary engineering fracture mechanics[END_REF][START_REF] Kanninen | Advanced fracture mechanics[END_REF]. It is the reason why, the parameter k is kept constant and equal to 0.38 in the present work. Moreover, by assuming a constant stress rate , the obscuration probability, from Equation ( 14), can be expressed as a function of the space dimension n, the Weibull modulus m, and the dimensionless time T/tc:

! ! ( ) 1 exp ( )! m n o c m n T D T P m n t +       ≡ = - -     +     , ( 17
)
where the characteristic time tc is given by [START_REF] Denoual | A damage model for the dynamic fragmentation of brittle solids[END_REF]:

1 1 0 m n m m n m n m n m n o c L o a d in g r a te C r a c k in c e p tio n C r a c k p r o p a g a tio n S k C t σ λ σ - - + + + -       =             & 1 42 4 3 1 4 2 4 3 1 4 2 4 3 . ( 18
)
For any value of m higher than 3 and considering (n = 3), it can be observed that the obscuration probability is less than 5% at (T = tc) and greater than 95% at (T = 2tc). Therefore, the characteristic time represents the time from which most of the obscuration phenomenon occurs. In addition, a characteristic stress is defined in the following way:

σc = σ °.tc, 1 1 0 m n n m n m n m n m n c o L o a d in g r a te C r a c k in c e p tio n C r a c k p r o p a g a tio n S k C σ σ λ σ - + + + -       =             & 1 424 3 1 4 2 4 3 1 4 2 4 3 , ( 19
)
and a characteristic density is defined as following:

λc = λt (σc), 1 1 0 m n m n m n m n m n m n m n c o L o a d in g r a te C r a c k in c e p tio n C r a c k p r o p a g a tio n S k C λ σ λ σ - - + + + -       =             & 1 424 3 1 4 42 4 4 3 1 4 2 4 3 . ( 20
)
By combining the closed-form solution of the non-obscuration probability (Eq. 17) and the Equations ( 7) and ( 8), the final crack density is expressed as a function of m, n and the characteristic density:

( )! ( ) 1 ! ! m m n cracks c c m n m m n m n λ σ σ λ + +     >> = Γ +     +     (21)
Under the assumption that the non-obscuration probability can be assimilated to the ratio between the obscured and total volumes, a damage variable D and a macroscopic stress Σ may be defined [START_REF] Denoual | A damage model for the dynamic fragmentation of brittle solids[END_REF]:

9 = D ' 9 and Σ 9 = E1 -9 F 9 (22), (23)
In the same way, by merging the closed-form solution of the non-obscuration probability (Equation ( 17)) and Equations ( 9) and ( 10), the ultimate macroscopic strength Σ G (maximum value of the macroscopic stress) and the time H G at which this maximal stress is reached can be expressed as function of the characteristic stress:

Σ G = I J =K L ! =! ! NO and H G = H I =K L ! =! ! NO (24), (25)
Equations ( 21), ( 24) and ( 25) are used in the next section to compare the fragmentation properties of three brittle materials.

Identification of the Weibull parameters and predictions of the DFH model for three brittle materials

The first material considered in the present study is an Ultra-High Performance Concrete (UHPC) called Ductal ® , provided by Lafarge-Holcim Company, characterised by a high compressive strength (170 MPa). This type of concrete benefits from a formulation with a reduced size of grains (≤ 0.6 mm), forming an optimal granular skeleton which provides a good homogeneity of the grains and matrix elastic properties. Such structural features allow decreasing internal stresses [START_REF] Richard | Composition of reactive power concretes[END_REF][START_REF] Cheyrezy | Microstructural analysis of RPC[END_REF]. In addition, the small water to cement ratio (≤ 0.2) enables to minimize the amount of porosity in the final concrete material (water porosity 5-6 %). In the considered commercial composition, steel fibres reinforce the concrete (length: 13 mm, diameter: 0.2 mm). This type of UHPC was studied in [START_REF] Blasone | Ultra-High Performance Fibre-Reinforced Concrete under impact of an AP projectile: parameter identification and numerical modelling using the DFHcoh-KST coupled model[END_REF] to simulate the penetration of a small calibre projectile against a 60 mm thick target. The UHPC mechanical behaviour was investigated using Quasi-Oedometric Compression tests, bending tests and spalling tests at the Hopkinson bar.

The second microstructure studied is an artificial porous polycrystalline granular ice grown in Laboratoire IGE (Institut des Géosciences et de l'Environnement). The specimens are grown from isotropic seeds made of crushed ice, having a maximum particle diameter of 2 mm, and surrounded by water at 0°C. The air trapped between the snow grains during the process leads to the formation of a highly porous microstructure, with about 7 to 10 vol.% of porosity. The final microstructure is characterised by equiaxed grains, an isotropic crystallographic texture and a mean grain size of about 1-2 mm [START_REF] Georges | Investigation of the porosity influence on dynamic tensile strength of polycrystalline ice, 24ème Congrès Français de Mécanique[END_REF].

The third microstructure, provided by Saint-Gobain, is a dense silicon carbide ceramic. This material corresponds to a grade of Hexoloy SA ® with size-controlled porosity. These relatively large and spherical pores are discrete, non-interconnected and dispersed in a controlled manner throughout the body of the material. Their proportion represents less than 2 vol. %, with a maximal diameter below 100 μm. This material has a density of about 3.11 g/cm 3 and is referred to as porous SiC in the present work.

24 bending tests were performed for each material. The tests on UHPC were carried out with 12×16×100 (Height × width × span) mm 3 Ductal ® beams without fibre reinforcement. A servohydraulic Instron testing machine with 100 kN load capacity was used for this material. The failure stress was directly derived from the maximum force. The strain rate of the bending tests was 5×10 -6 s -1 [START_REF] Blasone | Ultra-High Performance Fibre-Reinforced Concrete under impact of an AP projectile: parameter identification and numerical modelling using the DFHcoh-KST coupled model[END_REF].

The bending tests on porous ice were performed with the same testing facility than for UHPC samples (i.e. hydraulic press Instron with a 100 kN loading cell) coupled with a climatic enclosure supplied with liquid nitrogen to regulate the test temperature (set to -30°C) [START_REF] Georges | Etude expérimentale et modélisation de la réponse mécanique de la glace polycrystalline sous chargement de traction dynamique: influence de la porosité[END_REF]. The samples cross-section was 40×20 (height x width) mm ² with a support span of 120 mm. The crosshead speed was set to 4 µm/s.

Results of 3-points quasi-static bending tests performed on the SiC ceramic were provided by Saint-Gobain Research Provence (France). The dimensions of the tested specimens were 3×4×45 mm 3 (height × width × length), with a support span of 40 mm. The loading speed was set to 1 µm/s. As brittle materials are susceptible to crack nucleation at external and internal flaws, each surface of the specimens were polished and the corners chamfered.

The obtained distribution of failure stresses F σ allows building the so-called [START_REF] Weibull | A Statistical Theory of the Strength of Materials[END_REF] diagram in which [ ]

) 1 ln( ln F P - - versus ) ln( F σ
is interpolated by a linear function, the slope of which is the Weibull modulus m. Next, the effective volume is calculated according to the equation provided by [START_REF] Davies | The Statistical Approach to Engineering Design in Ceramics[END_REF], in order to take into account of the stress heterogeneity in the loaded volume:

P JQQ = R 〈 〉 NST = U VW , ( 26
)
where the symbol <.> corresponds to the Macaulay brackets, i.e. the positive value. In the case of three-point bending tests, the effective volume of the loaded structure is given by the analytical solution:

P JQQ = XYZ 4 =K [ , ( 27
)
where h is the height, b is the width of the tested beams and L is the bottom span length. Finally, the Weibull scale parameter 8 L /= is calculated from the following equation:

\ P JQQ = EP JQQ 8 F L N Γ =K = , ( 28
)
where \ is the mean bending failure stress and Γ the Eulerian function of the second kind. The parameters of the DFH model for the three brittle materials are summarised in Table 1.

Table 1. Parameters used in the DFH model for the three brittle materials.

A transition volume Ztransition between single and multiple fragmentation processes can also be defined as a function of the applied strain-rate, considering the equation

( \ "#$ %&"&' = Σ G ! ): . ! ! 1 ( 1)! m m m n transition c Z e m n m Z m n m +   +   = Γ     + -     (29)
The final crack density, the ultimate strength Σ G , the time corresponding to the ultimate strength H G and the single-multiple transition volume, provided by Equations ( 21), ( 24), ( 25) and ( 29), are plotted as functions of strain-rate (defined as ! = /B) for the three materials in Figure 5. To do so, the input parameters of the DFH model provided in Table 1 were considered. Strong differences are observed between the three materials, in terms of predicted ultimate strength (Figure 5a). However, it is interesting to remark that the final crack density, the time corresponding to the ultimate strength and the single-multiple transition volume are found to be relatively similar when compared at same strain-rate.

These predictions rely on the Weibull parameters identified from bending tests, which effective volumes are between the millimetre and centimetre scales. As explained in the next sections, these identification scales may lead to bad predictions when the number of cracks per effective volume significantly exceeds one. 24), ( 25), ( 21) and ( 29)). In this figure, the influence of strain-rate on the inception stresses of pores in UHPC concrete is not taken into account (see next sections).

Material

Identification of the density of critical defects in three brittle materials, based on X-ray tomography CT-scan

In the present work, the microstructure of each of the three studied brittle materials (concrete, porous ice and ceramic) is examined by X-ray micro-tomography, using two different X-ray devices. This method consists of sending an X-ray beam through an object and collecting the released attenuate beam via a panel detector. By placing the studied sample on a rotating platform, the 2D radiography of each slide of the object can be stacked to reconstruct a 3D image of the specimen. The X-ray attenuation level being linked to the atomic number of the crossed material, the grey levels of the resulting 3D image are directly linked to the density gradient in the material.

From this, the 3D image is binarized to highlight the internal structure of the material, such as porosity. A post-processing algorithm is used to label the identified pores (using the multidimensional image processing python package scipy.ndimage) [START_REF] Weaver | Centrosymmetric (cross-symmetric) matrices, their basic properties, eigenvalues, and eigenvectors[END_REF] and evaluate their volumes (using the label toolkit of the python package SPAM (Software for the Practical Analysis of Materials)). The pore equivalent diameter, determined from these analyses, corresponds to the diameter of a sphere having the same volume. Noise limits the ability to reliably identify pores

(a) (b) (c) (d)
smaller than a threshold size. For the SiC ceramic and the UHPC, the lower limit of neighbouring voxels required to adequately represent a single pore was fixed at 5 pixels in diameter and smaller pores were eliminated from the analysis. This criterion is in line with resolution limits presented in the literature [START_REF] Gualda | Quantitative 3D petrography using x-ray tomography: Application to Bishop Tuff pumice clasts[END_REF][START_REF] Tammas-Williams | XCT analysis of the influence of melt strategies on defect population in Ti-6Al-4V components manufactured by Selective Electron Beam Melting[END_REF]. With porous ice the high scanning resolution (7 µm) allowed a characterisation of pores as small as 3 pixels in diameter. Consequently, this threshold size was chosen for this material. Pores are classified in descending order of size and the cumulative density of pores is calculated as the cumulative number of pores divided by the considered volume and can be plotted as a function of the pore equivalent diameter. This methodology was applied to analyse the porosity present in each of the three studied materials.

Two UHPC cylindrical samples (diameter 29 mm, length 38 mm) were scanned using the X-ray micro-CT scanner EASYTOM XL Nanofoyer in the Laboratoire SIMaP. The X-ray source generates a polychromatic cone beam, which is detected by a flat panel detector of 2084 × 2084pixels, each pixel of the panel detector having a size of 27 µm. The X-ray source operated at 150 kV and 119 μA. The UHPC sample is scanned in a helical mode to maximize the resolution, 4,320 projections are acquired during three full platform rotations. This configuration allowed reaching a voxel size of 20 μm. Figure 6a shows a slice of the 3D reconstructed image of one specimen after thresholding. The isolated groups of voxels were identified as pores. The pore size distribution for a volume of 35×40×20 mm 3 , shown on Figure 6b, was obtained by summing up paralepidid crop of the two scanned specimens.

The second scan (Figure 6c) was performed on a cylindrical ice sample (diameter 45 mm, length 120 mm) in a cold room set at -20°C. To do so, the tomography equipment TomoCold DeskTom130 RX Solutions at the Laboratoire CNRM-CEN (Centres d'Etudes de la Neige, Grenoble), specifically adapted to cold temperatures, was used. A current of 238 µA and a voltage of 60 kV powered the X-ray tube. The detector was composed of 1920×1536 pixels with a physical pixel size of 127 µm. A ring filter with a 20-voxel kernel was applied in order to remove ring features, which are artefacts from the acquisition. After reconstruction, the scans are composed of approximately 4500 × 1660 × 1660 voxels. After the threshold step, a binary 3D image is obtained (slice on Figure 6c) and can be segmented using the python algorithm program in order to identify each individual pore with a diameter higher than the chosen cut-off value. The corresponding pore distribution is plotted on Figure 6d. processed 3D image contains 900 × 900 × 900 voxels, and results from the reconstruction of 1856 projections acquired during one full rotation of the specimen around the vertical axis. Such highresolution scan allows reaching a voxel size of 1.33 µm. Threshold and segmentation procedures were applied, similarly to both previous materials, in order to identify each individual pore as shown on Figure 6e. The distribution of pore size obtained is presented on Figure 6f.

Identification of the density of critical defects versus applied stress in three brittle materials

Methodology

In order to use the previous X-ray tomographic analysis in a fragmentation model, the cumulative pores distributions λt(Deq) presented in Figures 6b, 6d and 6f, need to be converted into a function connecting the density of critical defects to the applied microscopic stress. As previously introduced by several authors (Jayatilaka and Trustrum, 1977;[START_REF] Hild | De la rupture des matériaux à comportement fragile[END_REF][START_REF] Forquin | A method to determine the macroscopic toughness scatter of brittle materials[END_REF], the Linear Elastic Fracture Mechanics (LEFM) theory can be used to link the size of a pore to its expected crack-inception stress. According to this theory, the stress intensity factor KI, related to a crack of size (Deq = 2.a) is given by the following equation:

^_ = ` √bc , ( 30
)
where Y is a dimensionless correction factor linked to the crack geometry, often referred to as the geometric shape factor. In the case of a straight crack of length 2a, oriented perpendicularly to the loading direction in an infinite plate subject to a uniform stress field, this correction factor is (Y = 1). In the case of a penny-shaped crack of radius a in an infinite domain under uniaxial tension σ, the geometric shape factor value is (Y = 2/π). From this, the macroscopic stress to trigger an unstable crack is given by:

"#&ddJ# = e fg hAi jk /4 , ( 31
)
with ^_l the critical stress intensity factor of the tested material.

In the present work, it is proposed to determine the parameter Y thanks to the Weibull parameters identified from bending tests. Indeed, according to Equation ( 15), the density of critical flaws is expressed as a power law of the positive (tension) principal stress, with the Weibull modulus being the power value. The Weibull modulus is identified from the distribution of bending failure stresses. The Weibull scale parameter 8 -1/? and mean bending failure stress are determined from the effective volume and Equation (28), respectively. The Weibull linear regression (in a log-log plot) is reported on the plots of Figure 7a, 7b and 7c on the stress range corresponding to the bending tests. Moreover, by considering the equivalent diameter of a flaw to be triggered at a stress level ( ), the previous cumulated flaw density can be expressed as a function of ( ):

Jm = 4 i e fg h. 4 c/V 8 " E Jm F = 8 " E Jm , `, ^_l F , ( 32
), ( 33
)
where ^_l is the critical stress intensity factor provided in Table ( 1) and supposed to be constant whatever the flaw size and ` is the geometric shape factor to be identified from bending tests. From this last equation, it is possible to plot cumulated flaws density 8 " considering several values of the parameter `, as shown on Figure 7a, 7b and 7c. Next, the appropriate shape factor Y is determined by matching the 8 " curve from Equation (33) to the Weibull line (Equation ( 15)), for a stress equal to the mean tensile strength σw measured in bending tests. The calibration of Y also compensates the uncertainty regarding the critical stress intensity factor reported in Table 2. For each material, the identified value of Y is reported in Table 2. It is noted that for the three brittle materials this shape factor value is found to be close to one. In the next sections, this quantity is supposed to be constant for a given material and the Equations ( 32) and ( 33) are used for both the continuous and discrete modelling of the fragmentation. In the case of UHPC concrete a strain-rate sensitivity in the "quasistatic" regime needs to be considered [START_REF] Blasone | Ultra-High Performance Fibre-Reinforced Concrete under impact of an AP projectile: parameter identification and numerical modelling using the DFHcoh-KST coupled model[END_REF], so the 8 " curve varies with the applied strain rate as the crack-inception stress of the pores increases with the applied strain rate:

"#&ddJ# ! = "#&ddJ# , ( 34
)
where / corresponds to an exponent traducing the strain-rate sensitivity in the "quasi-static" regime and ! is a reference strain-rate (Table 1). This equation was already considered in [START_REF] Erzar | An experimental method to determine the tensile strength of concrete at high rates of strain[END_REF][START_REF] Erzar | Analysis and modelling of the cohesion strength of concrete at high strain-rates[END_REF] to account for the strain-rate sensitivity of wet concrete tensile strength due to the presence of free-water [START_REF] Rossi | A physical phenomenon which can explain the mechanical behaviour of concrete under high strainrates[END_REF]. According to the identification proposed in [START_REF] Blasone | Ultra-High Performance Fibre-Reinforced Concrete under impact of an AP projectile: parameter identification and numerical modelling using the DFHcoh-KST coupled model[END_REF], the following parameters can be considered (/ = 0.03 and ! = 5. 10 Lr ). Conversely, no strain rate sensitivity in the "quasi-static" regime is considered for the porous ice and the SiC ceramic. same critical flaws density for a stress equal to the mean tensile strength measured from the bending tests. For concrete, two levels of strain-rate are considered (5×10 -6 s -1 , 100 s -1 ).

Identification of the cumulative flaw density function for each brittle material

The DFH model presents the strong hypothesis that the critical flaws density is considered to be a power law of the positive (tension) principal stress (Equation ( 15)). A piece-wise function matching the cumulative flaw density expressed as a function of obtained from the tomography analysis, was identified with the aim of suitably including the full flaw population in the model. The calibrated interpolation function is given in Table (2) for each material.

Table 2. Parameters considered in Equations ( 32) and ( 33) and interpolation function of the logarithmic cumulative flaw density versus the logarithmic stress

In the case of UHPC (Fig. 7a), the critical stress intensity factor (0.627 sDc√?) was deduced from a Single-Edge Notched Bending (SENB) test performed in [START_REF] Forquin | Endommagement et fissuration de matériaux fragiles sous impact balistique, rôle de la microstructure[END_REF] with a notched sample of UHPC (Ductal ® ) of dimensions 5×3×20 mm 3 (h×w×L, with L the span length). Given the initial notch length a of 1.3 mm (α = a/h) and F the peak force reached in the test (24.2 N), the critical stress intensity factor K1c is deduced according to the following formula (ASTM D5045-14, 2014):

Material Next, a dimensionless parameter Y very close to 1 (Y = 0.83) provided the best fit of the Weibull line for a stress equal to the mean tensile strength of the bending tests (Figure 7a). For stresses close to the mean failure stress of bending tests (σw = 19.42 MPa), the slope of the curve, based on data from tomography analysis, matches quite well the slope of the Weibull straight line, which corresponds to the Weibull modulus (m = 8.9, Table 1). For concrete, as the inception stresses of pores are supposed to be strain-rate dependent according to Eq. ( 34), the critical flaws density versus strain-rate function also depends on the strain-rate. The critical flaws density versus applied stress is plotted in Figure 7a considering two strain-rates (the red curved dash-dotted line is corresponding to 100 s -1 ). The Weibull line of the porous ice, identified from bending tests, is compared to the data extracted from the tomography analysis in Figure 7b. For this material as well, a shape parameter very close to 1 (Y = 0.8) provides the best fit in a stress range close to the mean failure stress value (σw = 1.76 MPa, Table 1), and the slope of the curve is found to fit pretty well with the Weibull modulus (m = 5.3, Table 1).

^_l sDc√?) Y Interpolation function of Λ = tuv 6λ " w N x y 7 + , - z ≤ z → Λ = | z z & < z ≤ z &K → Λ = | & z | &K z | & z m ~+| &K z m ~ /m z > z → Λ = | z •'Hℎ i ⊂ {0, 1, … , n-1} | & z = c & z -z & + `& , c/V z =
The fracture toughness of the SiC ceramic with fuse flaws was measured on a SENB specimen, via a 4-point bending test (crack length 1.6 mm). The dimensions of tested specimen were 4×3×45 mm 3 (height × width × total length). The distances between the internal and external cylindrical supports were set to 20 and 40 mm, respectively. The obtained value, in this configuration, is 2.89 sDc√?. Finally, a shape parameter close to 1 (Y = 0.8) provides a good match between the tomographic data and the Weibull linear regression (Figure 7c) for a range of stresses close to σw (300 MPa, Table 1).

In addition, the slope of the curve coincides very well with the Weibull modulus (m = 16, Table 1). The high Weibull modulus value obtained for this material is due to the good size homogeneity of the large pores (~ 45 µm), which correspond to the flaw population triggered during quasi-static bending tests, therefore resulting in a low scattering in the flexural strength measurement.

The final interpolation function of the cumulative flaw density (Equation ( 33)) expressed as function of the parameters of Table (2) (^ I , `, z ' , `', c ' , ‹ ' ) is plotted according to the applied stress on the Figures 7a, 7b and 7c for each material (curved dash-dotted lines).

Predictions of the fragmentation properties of three brittle materials based on tomographic analysis:

The continuous method

Methodology

The obscuration probability of the DFH model, expressed in Equation ( 14), can be rewritten by using the size of the obscuration zone proposed in Equation ( 16) (with n = 3), but without making any assumption on the density of critical defects λt (σ(t)), except it is a continuous and derivable function.

ln 1 -D ' 9 = -R

•λ ž "

•"

Ÿ @ 9 -H • VH ;

(39)

Considering that ' 0 = 0 (with ' 9 -H = ŸE @ 9 -H F • ), the derivative of the previous equation with respect to time leads to (cf. appendix):

L… ¡ •… ¡ ; •" = R •λ ž " •" 3Ÿ @ • 9 -H 4 VH ; (40)
From the second derivative, the following differential equation is obtained [START_REF] Denoual | A damage model for the dynamic fragmentation of brittle solids[END_REF]:

• [ •" [ L… ¡ •… ¡ ; •" = 6Ÿ @ • λ " 9 (41)
By applying a triple integration of the equation ( 41), the evolution of the obscuration probability can be determined over time for a given loading history σ(t). When the obscuration probability is assimilated to a damage variable (Equation ( 22)), its determination provides the evolution of a macroscopic stress, according to Equation ( 23). From this, the ultimate macroscopic strength of the material (i.e. maximum macroscopic stress) can be deduced. In addition, the final crack density is obtained by integrating Equation ( 7). This methodology was applied to each brittle material, assuming a linear increase of the microscopic stress as a function of time ( 9 = . 9), and considering a large range of strain-rates, supposed proportional to the stress-rate (! = /B). The obtained results are detailed in the next subsection.

Application of the continuous method to the three brittle materials

The damage evolution law of the UHPC material was calculated by introducing the density of critical defects of Equation ( 35) into Equation ( 41). The considered values for the parameters (/ , ! , S, k, C) are listed in Table 1. The obtained ultimate macroscopic strength and final cracking density are plotted as functions of the strain-rate in Figures 8a and8b (solid blue lines). These predictions are compared to the one obtained using a Weibull density of critical defects (power-law function) (Equations ( 24) and ( 21), Figures 5a and5c). Whereas a good match is obtained at a low strain-rate, a higher strength level is predicted with the tomography population of defects at a high strain-rate compared to the DFH closed-form solution (Weibull flaws population). This difference can be explained by the fact that, at the high strain-rates, the defects activated are the smallest ones.

As densities of critical defects diverge at the highest stress levels (Figure 7a), the predictions made from the tomography data and the Weibull power-law function diverge. This precisely demonstrates the interest to consider rather the defects population identified from tomography than the Weibull one, provided by bending tests.

In terms of ultimate strength, the model predictions are compared to experimental results from [START_REF] Blasone | Ultra-High Performance Fibre-Reinforced Concrete under impact of an AP projectile: parameter identification and numerical modelling using the DFHcoh-KST coupled model[END_REF] obtained through spalling experiments using a Hopkinson bar apparatus. As shown on the plot of Figure 8a, the tensile strength given by the DFH model, based on the Weibull power-law function or on the defects distribution identified from tomography, are in good agreement with the experimental data at strain rates in the range of 50-200 s -1 . 1) and a flaw distribution identified from the tomographic analysis (Table 2), respectively.

The same methodology was applied to the porous ice material. The predictions for the ultimate macroscopic strength and final cracking density evolution with respect to the applied strain-rate are plotted on Figures 9a and9b. The ultimate strength predicted by the tomography analysis and Weibull-based DFH solution intersect at two different strain-rates. Below 30 s -1 the predicted strength with the continuous method is higher than with the closed-form solution, which corresponds to activated defects in the part I of Figure 7b. In a range of strain-rates from 30 s -1 to 500 s -1 it is observed that the predicted strength provided by the continuous solution is slightly lower than the one provided by the closed-form solution (Weibull distribution of flaws). In this case, the activated defects mainly correspond to those highlighted in the part II of Figure 7b. Both predictions are in a relatively good agreement in this strain-rate range in term of crack density (Figure 9b). At higher strain-rates, namely above approximately 400-500 s -1 , the flaws of the zone III start to be subsequently solicited resulting in a divergence between the final crack density predicted by the continuous method and Weibull-based solution as the flaw density predicted by a Weibull distribution is no longer relevant (Figure 7b). Similarly to UHPC, the ultimate tensile strength predictions of the DFH model are compared to experimental results (red dots in figure 9a) of spalling tests conducted with a Hopkinson bar apparatus and reported in [START_REF] Georges | Influence of porosity on ice dynamic tensile behavior as assessed by spalling tests[END_REF]. The experiments were performed with porous ice samples presenting a microstructure identical to the one considered in the present work. The increase of experimental tensile strength with strain rate, over the considered range of strain-rates (24-120 s -1 ), is in good agreement with the predictions provided by both approaches (Weibull-based solution and Tomography approach). In terms of tensile strength values, one can notice a slight over-estimation compared to the DFH model but predictions remain satisfying. Experimental data at higher strain rates would be of interest to validate the model predictions of ice tensile strength over a wider range of strain-rates. 1) and a flaw distribution identified from tomographic analysis (Table 2) respectively.

(b) (a)
The ultimate macroscopic strength of the SiC ceramic, with fuse-flaws, is plotted as a function of the strain-rate on Figure 10a. Two different behaviours can be distinguished in the present case, each one corresponding to different domains of the flaw density versus critical stress curve plotted on Figure 7c. For strain-rates below 20,000 s -1 , the continuous method and the DFH Weibull solution match well (Part I and II, Fig. 10a), which is due to a good similarity between the critical flaw density versus stress based on tomographic data and the Weibull slope, in the domain I of Figure 7c. For higher strain-rates, above 20,000 s -1 , the continuous method predicts higher ultimate strengths compared to the Weibull-based solution (Part III, Fig. 10a). This is explained by that fact that no pores of size between 20 -50 µm are available in the material to be activated (plateau identified as part II in Figure 7c) and the loading rate is not high enough to trigger smaller pores. So, the cracking density presents a plateau for a relatively large range of intermediate strain rates, between 10,000 and 40,000 s -1 (Figure 10b). Therefore, in this range of strain rates, the volume is fully damaged by the growth of obscuration zones spreading from the activation of the large sizecontrolled pores. At higher strain rates (> 40,000 s -1 ), the growth velocity of the obscuration zones from the large pores is not sufficient anymore to avoid the trigger of smaller flaws, so the cracking density starts to increase again. Very small pores activated in this high strain-rate domain mainly correspond to sintering defects or impurities, identified as part III in Figure 7c. Finally, the growth of obscuration zones from large flaws prevented the initiation of numerous micro-cracks in the material for loading rates between 10,000 and 40,000 s -1 .

One preliminary experimental result on the porous SiC ceramic is compared to the model predictions in Figure 10a. This result was obtained from a shockless plate-impact spalling technique developed to reach a controlled strain rate in the failure zone. This new experimental technique, based on wavy-machined flyer-plates, was developed in view of evaluating the strain-rate sensitivity of ceramic tensile strength (Dargaud and Forquin, 2020, submitted for publication).

According to this preliminary result, the spall strength of the porous SiC ceramic reaches 410 MPa for a strain rate of 18,150 ± 1,000 s -1 , which is in very good agreement with ultimate strength predictions given by the DFH model. Providing more experimental data for validation is an ongoing work. Compared to the spall strength values given in [START_REF] Zinszner | Strain-rate sensitivity of the tensile strength of two silicon carbides: experimental evidence and micromechanical modelling[END_REF] for a dense SiC ceramic Hexoloy SA (608 MPa at 17,000 ± 1,000 s -1 ), the presence of relatively large size-controlled pores led to a decrease of the material strength, at least for strain-rates lower than 25,000 s -1 .

Figure 10. Predictions of the porous SiC ceramic macroscopic strength (a) and final crack density (b) evolution according to the applied strain-rate. Comparison between the DFH close-form solution and the continuous and discrete methods, based on a Weibull flaw distribution (Table 1) and a flaw distribution identified from tomographic analysis (Table 2), respectively.

Predictions of the fragmentation properties of three brittle materials based on tomographic analysis: The discrete method

7.1 Methodology Discrete numerical simulations have been developed using the programming language Python. The routine is composed of the following steps. Each pore is identified from the tomographic scan using the Python package SPAM.A cut-off value was considered to distinguish an effective pore from potential noise of the tomographic scan. This size limit was set to 3 voxels in diameter for the porous ice and 5 voxels for the SiC ceramic and UHPC concrete materials. The coordinates of each pore are assimilated to the centre of mass position of ideal ellipsoids fitting the pores [START_REF] Ikeda | Three-dimensional study on the interconnection and shape of crystals in a graphic granite by X-ray CT and image analysis[END_REF]. The considered finite volume was obtained by summing up paralepidid crop of the two scanned specimens for the UHPC. The cubic scanned volume is cropped from a large cylindrical sample for the ice, and from beam of millimetric size for the ceramic. A numerical cube of equal size meshed with elementary cube elements of homogenous size is then used as a receptacle for the pores. Detailed parameters for the numerical simulation and the mesh settings are given in Table 3.

In order to simulate the material behaviour resulting from a tensile loading, a constant stress-rate is homogenously applied to the cube. The stress keeps increasing linearly in the whole domain, as long as none of the pore critical stresses is reached. Once the applied stress reaches the critical stress of one pore, a spherical obscuration volume centred on this pore starts to grow with a radius proportional to the constant speed kC. All the neighbouring pores reached by this growing obscuration sphere are obscured, meaning they cannot be activated anymore even if the loading stress happens to exceed their threshold activation stress. Thus, it is possible for a given strain rate to determine how many pores are triggers and how many were simply obscured.

In parallel, once the spherical obscuration zone of an activated pore reach a new mesh element while growing, the corresponding element is considered to be fully damaged. The elements of the cube not yet affected by any spherical obscuration volume keep a stress level increasing linearly.

The damage is defined as the number of damaged elements divided by the total number of elements. The macroscopic stress is calculated from the damage, as presented in Equation ( 9). From this, the macroscopic ultimate strength, which is defined as the maximum value of the macroscopic stress, is deduced. The discrete numerical simulation ends when the whole volume is obscured. With this algorithm, it is also possible to compute the crack density defined as the number of activated pores (for which the stress has reached its critical value without being obscured) during the process divided by the total volume. Finally, this algorithm is applied for a large number of strain-rates in order to evaluate the sensitivity of the macroscopic ultimate strength and crack density to the loading rate. For each brittle material considered, the input parameters and obtained results are detailed in the next subsection.

To better analyse the discrete numerical simulation, bar chart histograms, presented on Figures 11 (a-c-e), show, for three strain-rates and for each material, which are the stresses of the triggered pores at times Tu and Tmax, which corresponds to the time for which the ultimate stress is reached and the final time (volume fully damage) of the calculation, respectively. Histograms corresponding to low strain rates overlap the one of higher loading rates. Moreover, Figures 11 (b-d-f) present the evolution of the crack density as function of the applied microscopic stress that evolves linearly with time ( 9 = . 9), considering the same strain-rates. This figure illustrates which is the population of activated defects involved in the macroscopic strength and final cracking density determination, for various loading-rates.

7.2 Application of the discrete method to the three brittle materials First, the discrete method was applied to the UHPC material, considering the parameters listed in Table 3. A 3D volume 35×40×20 mm 3 was considered. This volume was defined as the largest volume provided by the tomographic scans. 15073pores, of size ranging between 100 µm and 2.5 mm were placed in the volume according to their exact position obtained from the X-ray CT scans (which represents about 2% of the total volume). This porosity represents the entrapped air during the casting process. The mesh element size of 0.5 mm was chosen as a compromised to capture the smallest element size providing a reasonable calculation time. This mesh size is much smaller than the average distance between the considered pores. For strain rates lower than 10 s -1 , the discrete approach predicts an ultimate strength that depends only on the weakest defect in the volume 35×40×20 mm 3 , therefore resulting in a unique crack density and a higher ultimate strength than in the continuous approach (Figure 8). This phenomenon is not captured in the continuous approach as this approach does not involve the size of the considered volume. A good match is observed between the predictions provided by the continuous and discrete methods between 10 s -1 and 1,000 s -1 , in terms of dynamic strength and crack density versus strain-rate evolution. Even if the continuous approach does not consider the exact position of the pores, it leads to results equivalent to the one obtained with the discrete approach. Above approximately 1,000 s -1 , it is noted that the ultimate strength and cracking density predicted by continuous and discrete methods start to diverge. This is due to the limited size of the smallest pores (0.1 mm) captured for the discrete method, whereas for the continuous method the pore distribution is extrapolated above this threshold. Therefore, the continuous approach predicts an increasing number of cracks, whereas the discrete approach reaches a plateau. It results in an estimation of the ultimate strength subsequently higher than the ultimate strength predicted by the continuous approach.

The same methodology was applied to the porous ice material with the parameters listed in Table 3. Similarly, to the UHPC material, a 3D cubic volume (31.32×31.32×31.32 mm³) was considered. 25,926 pores were spatially distributed in the volume at their exact positions, extracted from microCT data. The corresponding air-content is approximately 6.5 %. A mesh element size of about 0.3 mm was found to be satisfying regarding the pore sizes and the pore spatial distribution. The time step was kept constant to 50 ns for every strain-rate applied as, with this value, the total number of time increment was relevant to accurately describe the damage growth in the volume.

The low resolution of the microCT scan did not allow reaching higher strain-rates, the ratio activated pores/total number of pores being close to 1 above 1,200 s -1 . The results from the discrete method are compared to the ones obtained with the continuous method and the analytical solution of the DFH model (Weibull distribution of defects), in terms of ultimate strength and crack density in Figure 9(a) and 9(b). The ultimate macroscopic strength obtained from the discrete method is in excellent agreement with the continuous solution, whatever the strain-rate in-between 2 s -1 and 1,200 s -1 . Below 2 s -1 , a single fragmentation process occurs in the discrete method, whereas the continuous method leads to the crack density of a multiple-fragmentation process that would occur in larger volume. Above 1,200 s -1 , the discrete method is limited by the resolution of the CT scan. A reasonable good match is observed at intermediate strain-rate (around 50 s -1 ) when the fragmentation properties provided by both continuous and discrete methods is compared to the DFH closed-form solution. Indeed, the flaw density based on a Weibull distribution is a relatively good approximation of the actual flaw density extracted from microCT measurements, as long as the pores of the Zone III (Figure 7b) are not solicited. However, beyond this point, the crack density computed with the analytical distribution is not relevant anymore (especially above 400 s -1 ). It results in an estimation of a final crack density subsequently lower than the final crack density predicted by the DFH closed-form solution. However, it can be noted that, even at strain-rate of 200 s -1 , the pores from the Zone III are not significantly activated at the time Tu and have no influence on the values of ultimate strength determined from the continuous and discrete methods.

A series of discrete numerical simulations was conducted on the SiC ceramic with fuse-flaws considering the calibration and input parameters listed in Table 3. About 33,752 pores, with a minimal diameter of 5 pixels were identified via micro-tomography. The scan volume is cropped into a cube of size 1.2 mm 3 . A mesh of 100 elements in each direction is found to provide a sufficient level of discretisation to evaluate the ultimate strength and crack density. A step time of 10 ns was used for all the strain-rates in order to have a fine resolution of the damage growth. The ultimate strength and cracking density sensitivity to the strain-rate obtained with the discrete method are compared to the results of the closed-form solution based on Weibull parameters and the continuous method. It can be noticed that the discrete and continuous methods lead to similar curves on Figure 10, except for strain-rates below 1,000 s -1 and above 300,000 s -1 , for the reasons detailed hereafter. This first result confirms that the continuous function, identified in Table 2, properly describes the flaw population detected by tomography. Even if the continuous method does not take into account the exact position of the pores, it leads to results quasi-equivalent to the one obtained with the discrete method. For strain rates lower than 1,000 s -1 , the discrete method predicts a single fragmentation process in the cube of size 1.2 mm 3 resulting in a unique crack (λcrack × Vcube = 1) and a higher ultimate strength than in the continuous method (Part I in Figure 10a), which corresponds to the critical stress of the weakest defect in the volume (193.3 MPa). This single fragmentation process is not captured in the continuous method as this method does not involve the size of the considered volume. For strain rates between 2,000 and 10,000 s -1 , corresponding to the Part II in Figure 10a, the two methods provide the same estimation of the ultimate strength and cracking density than the closed-form solution based on a Weibull distribution of flaws. Indeed, in this range of strain-rates, the activated defects are the relatively large pores of controlled-size during the manufacturing process, also identified through the bending tests. For higher strain rates, the closed-form solution, based on a Weibull distribution of flaws, underestimates the ultimate strength and over-estimates the cracking density. This is because the Weibull distribution of flaws diverges from the intermediate and small-size flaws populations mainly activated at such strain-rates. Therefore, the closed-form solution of DFH model based on the Weibull parameters is no longer valid at such relatively high strain-rates. The flaw population "plateau" (Part II on Figure 7c) is correctly identified by both methods relying on the tomography analysis, thus both methods predict an increase of the material strength and a reduced number of crack inceptions for loading rates in the range 10,000-200,000 s -1 (Part III in Figure 10a). At strainrates above 200,000 s -1 , the discrete method is limited by the scan resolution, as pores with diameter lower than 6.6 µm (5 pixels) could not be accurately identified for a tomography scan resolution of 1.33 µm/voxel. Therefore, for this range of loading rates, whereas the continuous method predicts an increasing number of cracks thanks to an extrapolation of the pore distribution, the discrete method reaches a plateau (Part IV in Figure 10a).

To better understand which flaws population is triggered according to the applied strain-rate, the number of activated defects at Tu (time at which the ultimate strength is reached) and Tmax (time at which the whole volume is obscured) and the evolution of the crack density are plot as functions of the microscopic stress of the activated pores in Figure 11, for three different loading-rates.

For instance, in the case of UHPC, activated defects at a strain-rate of 100 s -1 have a microscopic stresses below 32 MPa at Tu and 42 MPa at Tmax. The activated defects above 1,000 s -1 correspond to microscopic stresses up to 100 MPa.

In the same way, in the porous ice, activated defects at strain-rates below 50 s -1 correspond to microscopic stresses less than 2.4 MPa at Tu and 4 MPa at Tmax. The activated defects at 500 s -1 are characterised by microscopic stresses up to 7 MPa at Tu and 10 MPa at Tmax.

In the case of porous SiC, the following strain-rates are considered: 2,000, 20,000 and 60,000 s -1 . It is observed that the crack density plateau around 100 cracks/mm 3 results from the plateau in the input flaw distribution. For low strain-rate, only large pores are triggered, therefore leading to a low ultimate strength. For intermediate strain rates, from about 10,000 to 50,000 s -1 , a clear increase of the stress is obtained while maintaining a relatively low number of cracks, as the material presents a limited quantify of intermediate pores (plateau) available to be triggered. The limited number of micro-cracks in this range of strain-rates plays a role on the kinetic of the damage variable, therefore resulting in a great increase of the materials strength. For very high strain rates, above 50,000 s -1 , the growth of obscuration volumes surrounding triggered cracks becomes too slow compared to the loading rate, so a very high number of small cracks are simultaneously nucleated out of these damaged volumes, resulting in a sharp increase of the final crack density. 73.1 (! =100 s -1 ) 6.99 655.9 Max. critical stress value (MPa) 100.9 (! =100 s -1 ) 10.79 717.9 Strain-rates (s -1 ) ; Time increments (ns) Finally, the representation given in the graphs of Figure 11 provides a better understanding of the population of flaws (defined as their activated stresses) governing the ultimate strength and final cracking density in each considered material according to the applied strain-rate. 7.3 Discussion: main advantages and drawbacks of each method Finally, the main advantages and drawbacks of each method (DFH closed-form solution based on Weibull parameters, the continuous and the discrete methods) are listed in the Table 4. Even if the continuous and the discrete methods need a tedious identification of defects from tomographic analysis, it provides the huge benefit to consider a non-Weibull distribution of flaws. So far, this last hypothesis has limited the predictions of the closed-form solutions of DFH model to a restricted range of strain-rates, corresponding to the flaws population activated in bending tests. -Defects have to be detected through X-ray CT scan -Resolution limit in the X-ray tomography (voxel size, ability to differentiate features from their absorption difference, etc.) -Need to identify a constant shape factor based on LEFM Continuous method -Differential equation of the obscuration probability is easy to program in a spreadsheet -Ease to implement in a FE code (non-homogeneous stress field) -Possibility to extrapolate the results -Assumption that the defects are homogeneously distributed in the volume Discrete method -Individual defects are considered at their exact position -Distribution of defects is not necessarily spatially homogeneous -Ease to identify which defects population is involved in the fragmentation process -No implementation in a FE code -Restricted to the scanned defects without possible extrapolation

The continuous and discrete methods were adapted to the study of the three present materials as they all presented a good phase contrast between the flaws and the matrix/grains. The continuous presents the advantage of being based on the obscuration probability differential equation, which can be easily implemented in a spreadsheet or in a Finite-Element code. This offers the opportunity to predict the strength and cracking density even in the case of non-homogeneous stress fields, as it is the case in edge-on impact tests of spalling tests for instance.

Although the discrete method cannot be easily implemented in a FE code as the continuous method, it presents the substantial advantage to consider the exact location and activation stress of each individual defect identified through the CT scan analysis. Predictions of the fragmentation process in brittle solids including non-homogeneous distribution of defects (in terms of size and spatial distribution) could be easily investigated through this method.

Conclusion

In the present paper, a new methodology is proposed to model the multiple fragmentation process induced in brittle solids subject to high-strain-rate loadings. This method is built on the identification of flaws population, provided by micro-tomography analysis. Samples of three different types of brittle materials, an ultra-high performance concrete, a porous polycrystalline ice and a SiC ceramic including size-controlled and relatively large pores have been analysed through an X-ray tomographic analysis. The pore size distribution of each microstructure was extracted from a post-processing of the 3D reconstructed images. By using the linear fracture mechanics theory, the density of critical defects is expressed as a function of the applied stress. Finally, the density of critical defects, obtained from tomography, is used through a continuous and a discrete method. In the continuous method, the differential equation of the obscuration probability of the DFH model is implemented. The calculation of a damage variable provides the evolution of the final cracking density and ultimate macroscopic strength according to the applied strain-rate, without making any assumption on the link between the density of critical flaws and the positive principal stress (power law function according to the Weibull model). In the discrete method, the growth of each spherical obscuration volume centred on its critical defect is simulated considering the exact position and activation stress (size) of each defect identified through the tomographic analysis. For each considered strain-rate, the macroscopic strength is calculated by considering the average value of stresses of all the elements not belonging to any obscuration volumes. The cracking density is obtained by considering the number of triggered defects divided by the total volume.

Both methods (continuous and discrete) were applied to the three brittle materials. In the UHPC concrete, new predictions are obtained at high strain-rates (above 100 s -1 ) that differ from the results provided by the DFH closed-form solution, relying on the Weibull parameters identified from the three-point bending test. For the polycrystalline ice, whereas a reasonably good match with the DFH closed-form solution is observed at intermediate strain-rate (for which the flaw density extracted from microCT measurements do not stray from the Weibull distribution), the predictions in terms of final crack density strongly diverge at the highest strain-rates (especially above 300 s -1 ).

For the SiC ceramic, the high-resolution micro computed-tomographic analysis enabled to capture the presence of two main flaws populations of distinct sizes. Again, it is observed that, above a certain strain-rate, the predictions of the DFH closed-form solution, based on a Weibull distribution of defects which correspond to the largest pores triggered in bending tests, diverge from the predicted results provided by both the continuous and discrete methods.

In addition, the discrete and continuous approaches applied to this porous SiC illustrate how the availability or unavailability of defects may change the strain-rate sensitivity of ultimate strength and crack density. Indeed, these approaches provided information on the valuable effect of the sizecontrolled fuse flaw in the porous SiC ceramic on the cracking density sensitivity to strain rate, for strain rate between 10,000 and 40,000 s -1 . Depending on the loading rate, a competition between the spreading of the obscuration zones and the initiation of new flaws drives the material strength. Such kinetic-dependant crack nucleation on pre-existing flaws explains the high rate-hardening strength, which is characteristic of brittle materials but magnified for this porous SiC ceramic with wellcontrolled pore sizes.

It illustrates how the use of CT scan to identify the flaws population leads to a better understanding of the effect of flaws population on the resulting strength of the material and its fragmentation, as function of the applied loading rate. Finally, the tomography technique is demonstrated to be a powerful tool to accurately take into account the real flaw distribution of the microstructure in order to model and predict the dynamic fragmentation process induced in brittle materials over a large range of strain-rates, as long as the flaws population can be characterised by tomography analysis. It is also observed that each (continuous and discrete) method presents complementary advantages. Indeed, whereas the continuous method is easy to implement and use in a spreadsheet or in a FE code in which transient dynamic loading can be considered (such as impact loading), the discrete method presents the advantage to consider the exact location and activation stress of each individual defects identified through the CT scan analysis, allowing predictions of fragmentation process with non-homogeneous distribution of defects. Finally, both continuous and discrete methods provide a promising way to better understand and model the effects of microstructural features, such as flaws distribution, on the macroscopic behaviour of strain-rate sensitive brittle materials. In particular, these methods are particularly promising with a view to optimizing the microstructure of brittle materials and their manufacturing process (sintering of ceramics, casting of concrete, surface treatment of glass, 3D printing process…) with respect to their use under high strain-rate dynamic loading. For each of the three-material studied in the present work, experimental validations of the modelling predictions constitute a natural prospect.

Figure 3 .

 3 Figure 3. Obscuration phenomenon and horizon of a point (M, T ).

Figure 4 .

 4 Figure 4. Illustration of the influence of the relative sizes of the domain Z and the obscuration zones on the type of induced fragmentation. A large domain Z leads to a multiple fragmentation process (left) because the horizon of a point (M, T ) do not interact with the domain boundaries, whereas a small domain Z (right) leads to single fragmentation process as the horizon of (M, T ) and the domain boundaries interact rapidly.

Figure 5 .

 5 Figure 5. DFH model predictions for the evolution of (a) ultimate strength, (b) time corresponding to the ultimate strength, (c) final crack density and (d) single-multiple transition volume of the three brittle materials as function of the strain-rate (defined as ! = /B). (Equations (24), (25), (21) and (29)). In this figure, the influence of strain-rate on the inception stresses of pores in UHPC concrete is not taken into account (see next sections).

Figure 6 .

 6 Figure 6. X-ray micro-CT scan of a sample of each brittle material. (a-b) UHPC, (c-d) Porous ice sample, (e-f) Dense SiC ceramic with size-controlled porosity. (a, c, e) Slice of the binary 3D reconstructed image. (b, d, f) Cumulative flaw density plotted as a function of the flaw equivalent diameter λt(Deq) (in log.-log. scale).

Figure 7 .

 7 Figure 7. Plot superposition of the critical flaws density, provided by the bending tests (Weibull law, black solid straight line), by tomographic analysis for three different values of the parameter Y, and interpolated function (black curved dash-dotted line) used in the continuous method (Equations (34), (35) and (36)), as a function of the applied stress, for (a) concrete, (b) ice and (c) ceramic. The final shape factor value Y is defined for each material when both first types of curves provide the

Figure 8 .

 8 Figure 8. Predictions of the UHPC concrete macroscopic strength (a) and final crack density (b) evolution according to the applied strain-rate. Comparison between the DFH close-form solution and the continuous and discrete methods, based on a Weibull flaw distribution (Table1) and a flaw distribution identified from the tomographic analysis (Table2), respectively.

  Pictures of the discrete approach at 20, 200 & 2000 s -1

Figure 9 .

 9 Figure 9. Predictions of the porous ice macroscopic strength (a) and final crack density (b) evolution according to the applied strain-rate. Comparison between the DFH close-form solution and the continuous and discrete methods, based on a Weibull flaw distribution (Table1) and a flaw distribution identified from tomographic analysis (Table2) respectively.

Figure 11 .

 11 Figure 11. Results of discrete numerical calculations. (a-b) UHPC, (c-d) Porous ice, (e-f) SiC ceramic with size-controlled porosity. (a, c, e) Overlapped histograms of the number of activated flaws at Tu (time at which the ultimate strength is reached) and Tmax (time at which the whole volume is obscured) for three levels of strain rates. (b, d, f). Evolution of the crack density as function of the microscopic stress for three levels of strain rates. The crack density level corresponding to the time at which the ultimate strength is reached is circled for each loading rate.

Table 3 .

 3 Parameters used for the discrete numerical simulations carried out the three materials

	Parameters	UHPC	Porous ice	SiC ceramic
	Cube size (mm)	35×40×20	31.32×31.32×31.32	1.2×1.2×1.2
	Element size (mm)	0.5	0.3132	0.012
	Number of elements	224000	1000000	10 6
	Lowest critical stress value (MPa) 20.2 (! =100 s -1 ) Average critical stress (MPa)	1.16	193.3

Table 4 .

 4 Advantages and drawbacks of the continuous and the discrete methods compared to the closed-form solution based on Weibull parameters (DFH model).

	Methods	Advantages	Drawbacks
	DFH closed-form	-No need of tomographic analysis and	-Predictions provided in a small range
	solution based on	LEFM	of strain-rate
	Weibull parameters	-No assumption on defects shape	-Need to perform bending tests
		-Closed-form solution	(difficult to perform at small scales)
	Common advantages	-Identification of defects through X-
	and drawbacks of	ray CT scan provides non-Weibull
	continuous and	distribution of flaws
	discrete methods compared to DFH closed-form solution	-Predictions valid over a much wider range of strain-rates

Fundings and acknowledgements:

This work was performed in the framework of the Brittle's CODEX chair, supported by the UGA (Univ. Grenoble Alpes) Foundation and sponsored by the Saint-Gobain and Lafarge-Holcim companies. These sponsors are gratefully acknowledged by the authors.

Appendix

If ones consider the following function:

its derivative reads

(2a)

(2b)

where

is the partial derivative of the function f regards to T. Considering (dT → 0) Equation ( 2) becomes:

Now, if ones consider the following function 9 :

Given that ' 0 = 0 , the first term of Eq. ( 3) vanishes and the derivative of Eq. ( 4) reads:

When assuming the following expression of ' 9 -H = Ÿ @ 9 -H , Eq. ( 5) becomes:

Finally, considering (n = 3) the Equation ( 39) is obtained from the differentiation of Eq. ( 38).