A Swindlehurst

A Saxena
email: aksaxena@uci.edu

A Mezghani
email: amezghan@uci.edu

I Fijalkow
email: inbar.fijalkow@ensea.fr

MINIMUM PROBABILITY-OF-ERROR PERTURBATION PRECODING FOR THE ONE-BIT MASSIVE MIMO DOWNLINK

Keywords: bit quantization, massive MIMO, linear precoding, perturbation precoding, minimum probability of error precoding

Linear precoders have been shown to perform reasonably well at low SNR when the basestation of a MIMO downlink employs one-bit digital-to-analog converters to quantize the precoder outputs. However, at medium-to-high SNRs, an error floor is encountered due to the coarse quantization. This paper examines methods for slightly perturbing the transmitted signal prior to quantization in an effort to improve downlink performance at higher SNRs. The perturbation is performed with the goal of minimizing the worst-case probability of error among the user terminals, and assumes that the symbols to be transmitted are drawn from a finite alphabet constellation. Two different types of perturbations are studied, and it is found via simulation that the methods can provide dramatic gains in downlink performance.

INTRODUCTION

The use of one-bit quantization has been proposed to reduce the cost, complexity and energy consumption associated with massive MIMO systems. Rather than reducing the number of RF chains and decreasing beamforming flexibility as in hybrid digital-analog approaches [START_REF] Han | Large-scale antenna systems with hybrid analog and digital beamforming for millimeter wave 5G[END_REF]- [START_REF] Alkhateeb | Mimo precoding and combining solutions for millimeter-wave systems[END_REF], one-bit quantization maintains one RF chain per antenna and reduces the cost per RF chain by replacing each high-precision analog-to-digital converter (ADC) or digital-to-analog converter (DAC) with a simple one-bit implementation. The savings are considerable since energy consumption and cost associated with the ADCs and DACs grow exponentially with the number of quantization bits, and other RF components can be simplified as well. This is especially important for the downlink, where a one-bit DAC obviates the need for linear amplifiers operating with a high back-off and hence low energy efficiency.

Most research on one-bit quantization for MIMO systems has focused on the uplink (e.g., see [START_REF] Liang | A mixed-ADC receiver architecture for massive MIMO systems[END_REF]- [START_REF] Mezghani | Capacity lower bound of MIMO channels with output quantization and correlated noise[END_REF]). The downlink case, with one-bit DACs implemented on each RF chain, has received considerably less attention [START_REF] Mezghani | Transmit processing with low resolution D/A-converters[END_REF][START_REF] Usman | MMSE precoder for massive MIMO using 1-bit quantization[END_REF][START_REF] Saxena | On one-bit quantized ZF precoding for the multiuser massive MIMO downlink[END_REF][START_REF] Jedda | Minimum BER precoding in 1-bit massive MIMO systems[END_REF]. The impact of one-bit quantization on the performance of linear precoding was studied in [START_REF] Saxena | On one-bit quantized ZF precoding for the multiuser massive MIMO downlink[END_REF], and it was shown that the signalto-quantization-and-interference ratio at the receivers grows approximately linearly with the ratio of the number of antennas to the number of users. Using one-bit DACs at the transmitter is equivalent to restricting the transmitted signal from each antenna to be QPSK. This observation was exploited in [START_REF] Jedda | Minimum BER precoding in 1-bit massive MIMO systems[END_REF] to directly design the transmitted signal for each possible combination of symbols desired at the user terminals. Rather than using a linear precoder, a nonlinear minimum probability-of-error criterion was used to directly determine the QPSK symbols for each antenna.

A key observation in much of the work for both the uplink and downlink is that at low signal-to-noise ratio (SNR), the loss in performance due to one-bit quantization is relatively small, only about 2dB (a factor of 2/π) in most cases. This is encouraging since the low SNR regime is a likely operating point for massive MIMO systems that take advantage of array or beamforming gain to make up for low signal power. However, it is still of interest to study methods that minimize the impact of one-bit quantization for medium and higher SNRs.

In this paper we focus on making linear precoding methods more effective in conjunction with one-bit quantization by perturbing the linearly precoded signal to favorably impact the probability of decoding error at the user terminals. We assume that the basestation is transmitting symbols drawn from a finite alphabet, and hence that each user needs only to make a correct detection of its symbol rather than receive something "close" to the symbol itself. With channel state information available at the transmitter, we can determine in advance the effect of the one-bit quantization at each user, and thus we choose an appropriate perturbation that increases the minimum distance to a decision boundary over all the users. Instead of searching for an optimal perturbation which would be prohibitively expensive, we show via simulation that making a relatively small and fixed number of perturbations to improve the performance of the worst-case user can lead to a significant improvement in overall system performance.

The next section introduces the data model assumed in the paper and the minimum symbol error probability criterion that we will use to update the perturbation for the transmitted symbol. Section 3 presents the perturbation methods considered in this work and how they are computed, and then simulation results are given in Section 4 to illustrate the performance gains compared with a transmitter with full-precision quantization and one with no perturbation.

QUANTIZED LINEAR PRECODING

Data Model and Assumptions

To illustrate the basic idea of our approach, we assume a flat fading downlink channel involving an M -antenna basestation (BS) and K single antenna user terminals. We will focus on a massive MIMO scenario in which M K, although this is not strictly necessary for the proposed methods. We make no assumptions regarding the K × M channel matrix H other than it is generically full rank. We let the M -element vector x represent the transmitted signal at the BS, so that the downlink data model can be represented as

r = Hx + n , (1)
where the k-th element of r, denoted by rk , represents the signal received by user k, and n is a vector whose K elements are independent identically distributed (i.i.d.) Gaussian noise. We will denote r = Hx as the noise-free version of r, and we will let H km denote the element of H in row k and column m.

In [START_REF] Saxena | On one-bit quantized ZF precoding for the multiuser massive MIMO downlink[END_REF], it is assumed that the vector x is the one-bit quantized output of a linear precoder:

x = Q(Ps) , ( 2
)
where s is the desired K-vector of symbols to be decoded at the users, P is an M × K precoder, and the one-bit quantizer Q(•) operates element-wise on its argument:

[Q(z)] m = 1 √ 2 (sign(Re{z m }) + jsign(Im{z m })) ,
where sign(•) returns the sign of the argument, and Re{•} and Im{•} represent the real and imaginary parts of their arguments, respectively. As explained below, in our approach we will perturb the signal in (2) in order to improve the probability of correctly detecting s at the user terminals. We will focus on perturbed versions of the maximal ratio transmission (MRT) scheme in which P = H H , and the zero-forcing (ZF) precoder where P = H H HH H -1 . In order to focus on the impact of the one-bit quantization, we assume in this paper that the channel matrix H is known at the transmitter. In practice the channel will be subject to estimation errors that will lead to degraded detection performance.

Minimum Probability of Error Precoding

We assume that the symbols that the BS desires to send to the users are drawn from a finite alphabet with defined symbol decision boundaries. Each user compares its received signal rk with these decision boundaries and takes the closest element of the alphabet ŝk as the symbol that was transmitted to it. As such, it is often not necessary that r k be "close" to s k in any sense, just that r k lie somewhere within the decision region associated with the desired symbol s k . In fact, for some constellations (e.g., PSK), a lower probability of error is obtained when r k moves farther away from s k such that its distance to the decision boundaries is larger. For example, consider the case of QPSK symbols depicted in Figure 1. In this figure, the desired symbol for user k is 1 + j in the upper-right quadrant, but any received signal rk that lies in the upper-right quadrant will be decoded as the correct s k , and the farther the received signal from the decision boundary, the lower the error probability. Received signals outside the correct quadrant have a positive distance to the boundary indicated by the lines with arrows, those inside the quadrant have a negative distance. While this distance is specific to a single user, its value will not only depend on s k , but also on the channel and and on the transmitted signal x, which in turn will depend on the other symbols s. Thus, we will let the function d k (x|s, H), more concisely d k , represent the distance of rk to the nearest decision boundary for user k. As discussed in [START_REF] Jedda | Minimum BER precoding in 1-bit massive MIMO systems[END_REF], there are several ways to aggregate these distances to provide a single figure of merit for the downlink. These include minimizing the sum k d k , prod-uct k d k , or (as we focus on in this paper) minimizing the maximum distance over all users:

x = arg min x max k d k (x|s, H) . (3)
Solving (3) over all 4 M possible x is prohibitively expensive, so in this paper we focus on choosing x to be a perturbed version of a quantized linear precoder. We will consider the two general classes of perturbations outlined below. Direct Perturbation (DP) -Here we choose x = Q(Ps)+ , where is an M -element perturbation vector that must be constrained such that the elements of x are restricted to be ±1 ± j. The idea is that by changing the symbols transmitted by a small number of the antennas, received signals near decision boundaries can be made more robust without affecting decisions made by other users.

Symbol Perturbation (SP) -In this case, the perturbation is made to the symbols: x = Q(P(s + )), to reflect the idea that it is not necessary to force recovery of the exact symbols s, just received signals within the correct decision regions. A small perturbation to one of the user symbols can help move the received signal for another user deeper into the decision region.

PERTURBATION PRECODING

In this section we present the details of the two approaches mentioned above. Even though these approaches work with a constrained set of possible values for x, for large M it is still unlikely that an optimization such as (3) can be implemented over all possible perturbations . Furthermore, the quantization non-linearity makes implementation of gradient-based methods problematic. Instead, we will show that simply considering a small constrained set of can yield a significant benefit rather than optimizing over all possible values of .

Direct Perturbation Precoding

The DP approach takes an initial quantized transmit vector x and strategically perturbs its elements to improve the performance of the worst-case user. The algorithm operates at each step by finding the user with the maximum distance d k , and then finding the element of x that makes the largest contribution to d k . Since

r k = M m=1 H km x m ,
the "worst" element of x from the point of view of user k is the one for which H km x m is the farthest from the decision boundary, or in other words the one for which d k (x m |s k , H km ) is largest. This element of x is thus adjusted to one of the other three possible QPSK symbols if that leads to a reduction in the distance for that user; otherwise it remains unchanged. This process is repeated a fixed number of times, with only one perturbation allowed per element of x. Mathematically, the implementation of the algorithm can be described as follows:

1. Initialize p = 0, number of perturbations N p , the index set M = ∅ and find an initial x (0) using (2) for some linear precoder P. p) and end; otherwise return to step 2.

Find k

(p) = arg max k d k (x (p) |s, H) 3. Find m (p) = arg max m / ∈M d k (p) (x m |s k , H km ) 4. Add m (p) to M 5. Set element m (p) of x (p) equal to the value ±1 ± j that minimizes max k d k (x (p) |s, H) 6. p → p + 1 7. if p = N p , set x = x (
Since it is relatively simple to update d k (x (p) |s, H) when only one element of x (p ) has changed, for large M most of the computational burden comes in step 3, where M separate values have to be checked for the worst-case contribution to r k (p) . The complexity of the algorithm can be controlled by choosing the maximum number of potential perturbations N p that are allowed.

Symbol Perturbation Precoding

Instead of directly perturbing the elements of x, the SP approach perturbs the values of the K symbols s that are used in the linear precoder. Also, rather than perturbing the symbols one at a time as we did with x, we choose a set of random perturbations { p , p = 1, • • • , N p }, and then select the p that results in the smallest value for d k . There are a number of possible perturbations that can be used; we have found that the following works well:

k = c(u k + jv k ) , (4)
where c is a constant gain and u k and v k are uniform random variables on [0, 1]. The details of the algorithm are given in the following steps:

1. Initialize p = 1, the "size" c of the perturbation, the number of perturbations N p , and find an initial x (0) using (2) for some linear precoder P. Calculate δ(0) = max k d k (x (0) |s, H). 2. Generate (p) according to (4). While this method eliminates the search over M elements for each p, it requires the computation of P(s + (p) ) and Hx (p) for each perturbation. Again, the computational load can be controlled by the choice of N p .

NUMERICAL RESULTS

For the numerical examples presented here, we assume a basestation with M = 128 antennas serving K = 8 singleantenna users. The channel H is composed of circular complex independent identically distributed (c.c.i.i.d) Gaussian random variables with zero mean and variance ρ. The additive noise at the receivers is also c.c.i.i.d Gaussian with zero mean and unit variance, and thus ρ represents the per-antenna SNR. The elements of the signal vector s are independent QPSK symbols. Our performance metric is the symbol error rate (SER) averaged over the K = 8 users.

Figure 2 shows the average SER achieved by several algorithms based on the MRT precoder. The curve labeled MRT refers to the MRT precoder with a full-precision DAC at the transmitter, Q-MRT refers to the quantized MRT precoder with no perturbation, while Q-MRT SP and Q-MRT DP refer to the direct-and symbol-perturbed versions of Q-MRT. Figure 3 shows the same set of results, but for algorithms based on the ZF precoder. The DP approach was implemented with 25 perturbations to the elements of the original x achieved by either Q-MRT or Q-ZF, and the SP approach used 25 different random perturbation vectors with c = 0.25 and at each trial chose the one that resulted in the best performance for the worst-case user.

We see that the perturbation precoders provide a gain in SER when the SNR is 0dB or higher, and the gain becomes more and more significant as the SNR increases. The SP version of Q-MRT eventually achieves the same performance as the unquantized MRT, and the DP version dramatically outperforms the unquantized MRT at high SNR. Both the SP and DP approaches perform similarly for ZF, and come within 3-4dB of the unquantized ZF precoder. Figure 4 shows how performance improves with the number of perturbations considered for the Q-MRT approach with DP, indicating that even better performance may be obtained if N p > 25 is possible.

CONCLUSIONS

This paper has presented methods for improving the performance of a massive MIMO downlink when the output of standard linear precoders is quantized by one-bit DACs at the basestation transmitter. The methods exploit the fact that the impact of the one-bit DACs can be determined prior to transmission, and they operate by making perturbations to the quantized downlink signal prior to transmission in order to improve the worst-case performance over the user terminals. Taking advantage of the assumption that the transmitted symbols are drawn from a finite alphabet, a performance metric for the pertubations is introduced that attempts to directly minimize the probability of decoding error, rather than for example minimizing the distance to the desired symbols. Simulations demonstrate that these perturbation approaches can achieve a dramatic improvement in performance at medium to high SNRs; in fact, the DP approach allows the quantized MRT precoder to significantly outperform standard MRT employing full precision DACs.

Fig. 1 .

 1 Fig. 1. Distance metric for minimizing probability of decoding error for symbol s k ; the distance is shown for several possible received signals rk .

3 .

 3 Calculate x (p) = Q P(s + (p) ) 4. Calculate δ(p) = max k d k (x (p) |s, H) 5. p → p + 1 6. if p < N p ,return to step 2. 7. Find p * = arg min p δ(p) and set x = x (p * ) .

Fig. 2 .

 2 Fig. 2. Average SER performance for MRT algorithms.

Fig. 3 .

 3 Fig. 3. Average SER performance for ZF algorithms.

Fig. 4 .

 4 Fig. 4. Performance improvement with N p for Q-MRT DP.

* This work was supported by the National Science Foundation under ECCS-1547155, and by the Technische Universität München Institute for Advanced Study, funded by the German Excellence Initiative and the European Union Seventh Framework Programme under Grant No. 291763, and by the European Union under the Marie Curie COFUND Program.


Damien André
email: damien.andre@unilim.fr

Bertrand Levraut

Nicolas Tessier-Doyen

Marc Huger

A discrete element thermo-mechanical modelling of diffuse damage induced by thermal expansion mismatch of two-phase materials

Keywords: discrete element method, DEM, thermal expansion, damages, silica, alumina, inclusion, brittle, crack, damage, Young's modulus

At the macroscopic scale, brittle media such as rocks, concretes or ceramics can be seen as homogeneous continua. However, at the microscopic scale, these materials involve sophisticated microstructures that mix several phases. Generally, these microstructures are composed of a large amount of inclusions embedded in a brittle matrix that ensures the cohesion of the material. These materials generally exhibit complex mechanical behaviors resulting from the interactions between the different phases of the microstructure.

As a result, the macroscopic behavior of these media may be predicted considering an accurate knowledge of their microstructures. This paper proposes a model to study the impact of diffuse damage resulting from thermal expansion mismatch between the mixed phases. This type of damage (which is not catastrophic for the integrity of two-phase materials) may appear when heterogeneous materials are subjected to thermal cycles.

This phenomenon involves a high amount of discontinuities and can not be tackled easily with the Finite Element Method (FEM). The Discrete Element Method (DEM) naturally accounts for discontinuities and is therefore a good alternative to the continuum approaches. However, the difficulty with DEM is to perform quantitative simulations because the mechanical quantities are not described in terms of the classical continuum theory. This study describes the approach used here to tackle this fundamental difficulty.

Introduction

Rocks, concretes or ceramics are heterogeneous materials exhibiting multiphase compositions involving different sizes of aggregates, various bonding phases and additives. Description and prediction of thermo-mechanical behavior of such materials present a real difficulty due to their complex microstructure. Hashin & Shtrikman (H&S) have developed an analytic method to predict thermo-mechanical properties of perfectly cohesive (free of damage) multi-phase materials [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF]. However, most of these materials present numerous micro-cracks at room temperature. These micro-cracks result from thermal expansion mismatches between their different phases. These local defects highly influence thermo-mechanical properties and may have a great influence on materials lifetime in service conditions.

For instance, Young's modulus is strongly affected by the presence of microcracks and the measured values are often in disagreement with H&S's prediction [START_REF] Tessier-Doyen | Untypical young's modulus evolution of model refractories at high temperature[END_REF][START_REF] Joliff | Experimental and numerical study of the room temperature elastic modulus of model materials with partly bonded matrix/particles interfaces[END_REF]. The study proposed here focuses on a numerical method able to predict the occurrence of these defects and their overall influence on macroscopic properties such as Young's modulus and Coefficient of Thermal Expansion (CTE). Any materials that exhibit coarse aggregates such as concretes, rocks or ceramics may be addressed by the presented method.

Various numerical approaches can be distinguished. In the mechanical engineering field, the most used is the Finite Element Method (FEM). FEM, which is the most widespread technique is well suitable for materials characterized by reliable stress-strain laws. However, at the microscopic scale, this method is not adapted to describe discontinuities (4) without assumptions on their localization, their paths, their growths and their coalescences. Alternative methods to FEM exist such as X-FEM (Extended Finite Element Method) that keeps the rate of convergence towards singularities (micro-cracks) thanks to form functions. However, this method is adapted to describe the crack propagation in quasi-static mode but is not suitable for the management of opening and closure of numerous microcracks simultaneously.

The DEM is an interesting alternative way to study multi-damaged materials because it takes naturally into account discontinuities. The DEM implements a group of distinct elements (also named discrete element) which are in interaction through contacts or cohesive laws. This model consists in an assembly of discrete elements, deformable or not, linked by simple mechanical laws. The advantage is the natural description of crack initiation, crack propagation and coalescence. Researchers have used this method to study damages in solids, such as concretes [START_REF] Hentz | Discrete element modelling of concrete submitted to dynamic loading at high strain rates[END_REF], rocks [START_REF] Bobet | Numerical models in discontinuous media: Review of advances for rock mechanics applications[END_REF] or ceramics [START_REF] Tan | Discrete element method (dem) modeling of fracture and damage in the machining process of polycrystalline sic[END_REF]. However, the main difficulty is to get quantitative predictions. This is due to the necessity to find relations between microscopic laws, at the discrete element scale, and the macroscopic properties, at the structure scale. In this study the considered macroscopic properties are Young's modulus, Poisson's ratio, coefficient of thermal expansion and failure strength.

In this study, only the thermo-mechanical effects are taken into account. To avoid complex physico-chemical interactions, a model material is preferred in place of an industrial one. The phases of this model material have been chosen to ensure chemical neutrality and to avoid the presence of interphases. This model material is composed of two phases : spherical alumina inclusions embedded in a borosilicate glass matrix. The selected model material exhibits a positive CTE mismatch between the glass matrix and the alumina inclusions.

A short first section is dedicated to model material processing. The second section focuses on the H&S analytic model whose limits are highlighted. The third and fourth section introduces the DEM model. In this part, the method to model thermal expansion with DEM, that presents an innovative feature, is focused. As a preliminary validation, the next section proposes the implementation of the DEM model without cracks. The numerical results are compared to the H&S predictions. In the last part, the micro-cracks are implemented and the given numerical results are compared to the experimental data in terms of macroscopic Young's modulus and CTE. The capability of the presented method to simulate, understand and predict qualitatively and quantitatively the impact of the microcrack density on thermo-mechanical properties is finally discussed.

Experimental study

The material used for this work is a model material. A model material mimics, through a simplified framework, a given discoupled phenomenon observed with real and complex materials. In order to study the impact of thermal expansion mismatches, a two-phase model material, composed of mono-diameter spherical alumina inclusion embedded in a borosilicate glass matrix, is preferred. The thermo-mechanical parameters values for alumina and glass were chosen to produce a micro-crack network during the cooling stage of the sample preparation.

Properties

Processing of two-phase model materials

The model materials used in this study are composed of spherical monomodal alumina inclusions (average diameter equals to 500 µm) which are randomly placed in a borosilicate glass matrix. Alumina (corrundum phase) is a single oxide-based ceramic exhibiting fine grains. The main requirements of the selected glass matrix are homogeneity, isotropy, a rather chemical inertia and the capability to monitor thermal expansion coefficient. In this way, a borosilicate glass has been prepared from the melting of a vitrifiable mixture initially constituted by different raw materials containing silica, boron oxide and other secondary oxides. After grinding and sieving at 40µm, the glass powder was mixed during 3 hours with a slight proportion of organic additives used as binder and lubricant. A perfectly controlled volume fraction of alumina inclusions is incorporated in the mixture and is homogenized during 1 hour to ensure the dispersion of spherical alumina inclusions. Green specimens (80 × 40 × 10 mm 3 ) are shaped by uni-axial pressing (80 MPa) before debinding and sintering under uni-axial pressure (15MPa at 600 • C) to remove residual porosity. Three different volume fractions of inclusion were prepared (15%, 30% and 45%). Figure 1 shows the microstructure of a final two-phase model material highlighting the crack network.

Induced residual thermal stresses

The thermo-mechanical properties of both individual materials are given in Table 1. The selected model material exhibits a thermal expansion mismatch such as :

∆α = α matrix -α inclusion , with ∆α = 4 × 10 -6 K -1 so ∆α > 0 (1)
The introduction of spherical particles in a matrix, which are chemically and physically different, leads the occurring of thermal stresses during the cooling

Inclusions Matrix

Cracks Figure 2 describes this phenomenon for a single inclusion. The matrix is subjected to σ orth and σ rad . The brittleness of the glass matrix induces orthoradial cracks that occur and propagate in the matrix. In the case of a unique inclusion, the radial and orthoradial stresses σ rad and σ orth can be deduced from the hydrostatic pressure p at the interface (8):

p = (α m -α p ) ∆T 1 + ν m 2E m + 1 -2ν p E p -1 (2)
with :

σ rad = -p and σ orth = p 2 (3)
with :

• 'm' and 'p' indices are respectively the matrix and the inclusion,

• α i is the CTE and When multiple inclusions are presents, a complex crack network is created. Figure 1 shows an example of this crack network.

The Hashin & Shtrikman model

Among the various theoretical approaches developed to predict elastic and thermal properties of two-phase materials, the H&S bounds are used here. This analytic approach seems well suited to describe inclusional materials [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF]. The apparent Young's modulus can be predicted from the knowledge of (i) both volume fractions of the matrix (v m ) and inclusional particles (v p ), (ii) bulk and shear modulus of each phase noted K m , G m , K p , and G p respectively. The bounds of the apparent bulk modulus K - v and K + v are computed as :

K - v = K m + v p 1 Kp-Km + 3vm 3Km+4Gm (4)
K + v = K p + v m 1 Km-Kp + 3vp 3Kp+4Gp (5)
The bounds of the apparent shear modulus G - v and G + v are computed as :

G - v = G m + v p 1 Gp-Gm + 6(Km+2Gm)vm 5Gm(3Km+4Gm) (6)
G + v = G p + v m 1 Gm-Gp + 6(Kp+2Gp)vp 5Gp(3Kp+4Gp) (7)
The bounds of apparent Young's modulus E - v , E + v are computed as :

E - v = 9K - v G - v 3K - v + G - v (8)
E + v = 9K + v G + v 3K + v + G + v (9)
As shown in Figure 3, the plots of the experimental Young's modulus and CTE are lower than H&S lower bound (HS -) prediction. The Young's modulus was measured through an ultrasonic technique and the CTE with a standard dilatometer equipment.

The presented results underline the impact of the micro-cracked microstructure on the macroscopic thermo-mechanical behavior. The micro-cracks lead to a significative loss of E and α values compared to the H&S prediction. When the volume fraction of inclusions (15%, 30% and 45%) increases, the results for E values tend to enhance and for α values, they tend to decrease.

The following parts of this paper will introduce the discrete element model used to predict apparent Young's modulus and apparent CTE that can not be treated with H&S models.

The discrete element model

Thanks to the spectacular improvement of computer performances in the last two decades, DEM approach has met a growing interest in the scientific community. Many methods based on the original DEM proposed by Cundall and Strack

HS-

Δα

(a) Young's modulus E v [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] have been proposed in the literature. Among these methods, the lattice model allows to study a continuum as a set of simple physical interactions between discrete elements. For example, Schlangen in [START_REF] Schlangen | Experimental and numerical analysis of micromechanisms of fracture of cement-based composites[END_REF] proposed to study the fracture of concretes with discrete elements bonded by mechanical beams. With such a model, a major problem is that the discrete elements are material points without volume nor shape. An interesting approach is to mix lattice and particle models as it was first proposed by Potyondy in [START_REF] Potyondy | A bonded-particle model for rock[END_REF]. It has been widely used to study phenomena that cannot be treated with the Finite Element Methods (FEM) such as tribological problems with wear phenomena. Unlike continuous approaches, the main difficulty for DEM is to simulate quantitatively the continuum. The purpose of the free DEM software GranOO1 is to face this difficulty. At this time, GranOO embeds models that allow quantitative simulations of continuum materials with DEM [START_REF] Andre | 3D Discrete Element Workbench for Highly Dynamic Thermo-mechanical Analysis[END_REF]. This section introduces the main aspects of the numerical models used in GranOO. In addition, the developments specifically addressed for the presented study are highlighted : the thermal expansion and an improvement of the failure criterion.

Overview of the numerical procedure

Figure 4 shows the main simulation steps implemented in this study. The first step consists in calibrating the thermo-mechanical parameters of the discrete element model. These microscopic parameters, related to the scale of the discrete elements, are denoted by the µ indice. These parameters are the microscopic Young's modulus E µ , radius ratio R µ , tensile strength σ µ f and CTE α µ . Considering the values reported in Table 1, both borosilicate glass and alumina parameters are calibrated separately. Readers may refer to the sections 4.4, 4.5 and 4.6 to get more details about the signification of theses microscopic parameters and the related calibration method.

After building the initial cubic discrete domains, the second step consists in inserting the spherical alumina inclusions. This step involves a simple geometrical algorithm that is not described here.

The aim of this work is to study the influence of the cracks that happen during the cooling stage of the sample preparation. This stage concerns the range of temperature [450 • C -20 • C] from the glass transition temperature to room temperature.

Finally, the damaged virtual samples given by the last step are frozen. In this case, only the elastic behavior is taking into account. To forbid the crack extension, the tensile strength is set to an infinite value. These frozen samples are submitted to virtual tensile tests (see sections 5.1 and 5.2) to evaluate their apparent Young's modulii E M .

Synopsis of the main assumptions

In the following discrete element model the main assumptions are :

• the temperature is considered as constant in the whole domain, the heat transfer is neglected because the studied experimental samples are submitted to a very low thermal loading variation.

• the contacts between the discrete elements are not taking into account, the proposed model is a lattice discrete element model. However, the proposed model must be compatible with discrete element model that manages contacts.

• due to the positive CTE mismatch, only the crack opening is considered. In a first approach, the crack closure is supposed insignificant and is not taken into account. • a crack occurs locally only in a tensile mode. The failure in compression is not allowed locally at the discrete element scale.

• the slight variation of the Young's modulus in respect with the temperature is not considered.

The explicit temporal scheme

The temporal differential equation is solving by using an explicit integration scheme that is well suited to massive DEM simulations [START_REF] Rougier | Numerical comparison of some explicit time integration schemes used in dem, fem/dem and molecular dynamics[END_REF] and high velocity phenomena. The velocity Verlet scheme is chosen for its simplicity and efficiency. The discrete element orientations are described by quaternions. Quaternions provide an efficient way to compute discrete element attitudes (14, §2.5). Quaternions must be normalized for every time step to prevent numerical drift. A pure numerical damping factor may be introduced in the numerical scheme as described by Tchamwa and Wielgosz [START_REF] Mahéo | Damping efficiency of the tchamwawielgosz explicit dissipative scheme under instantaneous loading conditions[END_REF] to increase the convergence rate toward static solutions.

To avoid divergence of the temporal scheme, the chosen time step ∆t must be lower than the smaller oscillation period t os of the whole system. For an elastic interaction the period is computed as:

∆t = 1 k min(t os ) with k ≥ 2π and t os = M ij K ij (10)
where M ij and K ij are respectively the equivalent mass and stiffness of an interaction between a discrete element i and j. Cambou in ( 16) proposes a k value equal to 2π. In this study, the stiffness is related to the cohesive beam interactions (see section 4.4) and is equal to :

K ij = E µ S µ l µ (11)
where E µ , S µ and l µ are respectively the cohesive beam Young's modulus, section and length. The equivalent mass of an interaction between a discrete i and j is :

M ij = M i M j M i + M j (12)
where M i and M j are the mass of the discrete element i and j. The time step used in this work is about the nanosecond.

Initial domain building

In this study, only cubic discrete domains are used. The discrete domains are initially packed. The radius of the discrete elements follows an uniform distribution. The range of the distribution is 25%. This value prevents from creating ordered discrete domains that lead to geometrical anisotropy [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF]. The algorithm used to pack the initial domain comes from the GranOO workbench [START_REF] Jebahi | Discrete Element Method to Model 3D Continuous Materials[END_REF]. The packing process involves three main steps:

1. the given boundary shape (a cube) is filled with discrete elements randomly placed in the domain. The filling process goes on as long as free places are found. 2. even if no free place exists, a few number of discrete element are inserted.

It induces a high level of inter-penetration. A new equilibrium state is achieved by performing a DEM computation. This second step is repeated while the wanted number of discrete element is not reached. 3. finally, the bond network is constructed thanks to a Delaunay triangulation.

The nodes of the Delaunay triangulation are the center of the discrete elements. The edges given by the Delaunay triangulation are used to build the bond network.

Elastic behavior

The cohesive beam bond model ( 17) is used here to simulate perfectly elastic media characterized by a Young's modulus and a Poisson's ratio. The bond network given by the initial domain building (see previous section) is replaced by a beam network.

In such model, the discrete elements are bonded by Euler-Bernoulli beams. These beams are able to work in tensile, bending and torsion modes. The interactions resulting from the cohesive beams at the discrete element scale are denoted microscopic. They are symbolized by the µ suffix. In opposition, the emergent behaviors of the whole assembly network are denoted macroscopic (see Figure 5). The macroscopic scale is symbolized by the M suffix. The macroscopic behavior corresponds to the behavior of the simulated material. The cohesive beams are simply defined by two microscopic parameters : a Young's modulus E µ and a radius ratio R µ . The E µ and the R µ values are quantified thanks to a calibration process. The aim of this process is to reach the required values of macroscopic Young's modulus E M and Poisson's ratio ν M . The calibration process can be achieved in two main steps :

1. the calibration of the microscopic beam radii in order to merge the imposed value of Poisson's ratio ν M and 2. the calibration of the microscopic beam Young's modulus in order to reach the required value of Young's modulus E M .

More details about the calibration process are given in the section 5.3.

Brittle behavior

Within the discrete element approach, cracks are generally simulated by breaking the cohesive bonds when an imposed criterion is reached [START_REF] Tavarez | Discrete element method for modelling solid and particulate materials[END_REF][START_REF] Shiu | Discrete element modelling of missile impacts on a reinforced concrete target[END_REF]. The main existing approaches are based on the computation of bond strains [START_REF] Carmona | Fragmentation processes in impact of spheres[END_REF] or stresses [START_REF] Potyondy | A bonded-particle model for rock[END_REF]. However, these approaches are not able to simulate complex crack paths such as the hertzian cone crack that exhibits on fused silica glass [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF]. In a previous study, a new criterion, based on the computation of an equivalent Cauchy stress tensor and a maximal hydrostatic pressure value was developed. It has been shown that this criterion gives more accurate results than the standard DEM failure criteria. The equivalent Cauchy stress tensor is computed as :

σi = 1 2Ω i 1 2 j r ij ⊗ f ij + f ij ⊗ r ij (13)
where :

• ⊗ is the tensor product between two vectors,

• σi is the equivalent Cauchy stress tensor of the discrete element i,

• Ω i is the volume of the discrete element i,

• f ij is the force imposed on the discrete element i by a cohesive beam that bonds the discrete element i to another discrete element j and,

• r ij is the relative position vector between the center of the two bonded discrete elements i and j.

The developed criterion assumes that fracture occurs if the computed tensile hydrostatic stress is higher than a threshold value σ f :

1 3 trace( σi ) ≥ σ f (14)
In the initial study [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF], all bonds that belong to the discrete element that reach the σ f criterion are destroyed. As a consequence, a debonded discrete element occurs and a debris is created. However, the loss of discrete elements when fracture occurs can be viewed as a loss of mass and volume. In the following study, only the half of bonds that belong to the discrete element is destroyed. Hence, a second computation is necessary to select the bonds to break out.

Because brittle failure occurs in tensile mode (mode I), the breaking bond set must follow a surface normal to the higher tensile stress direction. The higher tensile stress direction d i is assumed to be the direction of the vector formed by the principal stress components d i (σ I i , σ II i , σ III i ) related to the discrete element i. The surface (S i ) created by the fracture is assumed to be normal to the main principal stresses direction d i . To mimic this surface creation, the bonds between the discrete elements i and j are broken assuming the following condition :

r ij • d i > 0 ( 15
)
where • is the scalar product and r ij is the relative position vector between the bonded discrete i and j. This method allows to create a crack path along the surface S i . Figure 6 illustrates this process in 2D. The main advantage of this method is to let the possible management of crack closures by taking into account the contacts between unbonded discrete elements. This last aspect (crack closure) is not treated in this study.

Thermal expansion behavior

All materials subjected to a thermal loading undergo thermal expansion. For a free of crack material, if the temperature is homogeneous in the whole volume, thermal strains do not lead to stress and the mater is still relaxed. As a consequence, imposing forces or torques to discrete elements is not the most valuable way to treat thermal expansion phenomenon using DEM.

The proposed solution consists in incrementing the free length of the cohesive beams l 0 by a ∆l 0 value computed from the temperature variation ∆T (see Figure 7). It leads to the expansion of the whole domain without adding any internal forces. So, free length of cohesive beams become a function of the temperature :

l 0 (T ) = l initial 0 (1 + ∆T × α) (16)
Discrete element i

d i
Fracture surface (S i )

Broken beams Discrete element j To avoid confusion, it is important to note that the current lengths l are not directly affected by these computations. Changing the free length l 0 , for the given cohesive beam, leads to change the reaction forces acting on the bonded discrete elements. If the related discrete elements are unconstrained, the current length l of the beam converges toward the most updated value of its free length l 0 after a given number of iterations.

Numerical study of thermal expansion

As the other mechanical parameters described in the previous sections, the thermal linear expansion α is considered at two levels : microscopically (α µ ) and macroscopically (α M ). The α µ is the input parameter related to the thermal expansion coefficient of the cohesive beams and the α M is the resulting thermal expansion of the bounding shape of the discrete domain.

To study the relation between the microscopic and the macroscopic thermal expansion, some simulations were carried out. The discrete domain used for these simulations is shown in Figure 8. The discrete domain contains around 10,000 discrete elements and the average number of bonds per discrete element is 13.6 (also named coordination number).

l 0 k (T ) = l initial 0 k (1 + ∆T × α µ ) (17)
In this whole study, the temperature is considered as homogeneous in the material. So, in a first approach, heat conduction is not taken into account.

To measure the macroscopic thermal expansion coefficient α M , the length variations of the bounding box are computed and updated during the simulation. This bounding box has six faces : S x+ , S x-, S y+ , S y-, S z+ and S z-.

In this whole study, the first index 'x', 'y' or 'z' means that the considered surface is normal to the related axis : x, y or z. The second index, minus '-' or plus '+', denotes the location of the face according to the axis direction. For example, S x+ denotes the face normal to the x axis and located positively along this axis.

The discrete elements that belong to these faces were identified by the preliminary packing algorithm used to build the initial discrete domain (for details about the packing algorithm, see ( 18), ( 23) and section 4.3).

The lengths L x , L y , L z along the x, y and z axes of the bounding box are computed by averaging the position vectors of the discrete elements that belong to these faces :

L x = p i • x i∈S x+ -p i • x i∈S x- (18)
L y = p i • y i∈S y+ -p i • y i∈S y- (19)
L z = p i • z i∈S z+ -p i • z i∈S z- (20)
where p i is the position vector of the discrete element i.

According to the equation 16, three macroscopic linear thermal expansion coefficients α Mx , α My and α Mz along the x, y and z axes are computed as :

α Mx = ∆L x ∆T × L initial x (21) α My = ∆L y ∆T × L initial y (22) α Mz = ∆L z ∆T × L initial z ( 23
)
where ∆L x , ∆L y , ∆L z are the length variation due to thermal expansion of the bounding box.

Figure 9 shows the results of such numerical simulations. The applied thermal loading corresponds to an increase of temperature of 1/20 Kelvin per iteration. The simulation is achieved after 10,000 iteration that corresponds to a total temperature variation of 500 K. The microscopic thermal expansion coefficient used for these simulations is α µ = 7.6 10 -6 K -1 . The evolutions of the macroscopic linear thermal expansion coefficients α Mx , α My and α Mz versus the iterations number are plotted in Figure 9. The given results show that :

1. the evolutions of the three macroscopic linear thermal expansion coefficients α Mx , α My and α Mz are quite similar and 2. the macroscopic thermal expansion coefficients converge toward the value of 7.6 10 -6 K -1 corresponding to the microscopic one.

These two results conduct to : The first result ensures that thermal expansion is isotropic and the second result shows that thermal expansion coefficient can be introduced directly without any calibration.

α µ = α Mx = α My = α Mz = α (24) E (GPa) ν σ f (MPa) α (10

From discrete to continuous quantities

Continuous quantities such as Young's modulus, Poisson's ratio, tensile failure strength or coefficient of thermal expansion are computed from discrete quantities. This section focuses on the main methods for computing these quantities.

Because continuous quantities are always considered at the scale of the whole discrete assembly (macroscopic scale), these quantities are not directly computed. Continuous quantities are virtually measured during numerical testings. For instance, mechanical quantities such as Young's modulus, Poisson's ratio and tensile failure strength are retrieved through virtual tensile test. The CTE is computed through a virtual thermal expansion test.

A calibration procedure must be achieved to set the microscopic parameter values that match the values of the two simulated constituents : alumina and borosilicate glass. The thermo-mechanical parameter values of these two materials are shown in Table 2.

According to the previous section, the linear thermal expansion coefficient α has not to be calibrated. However, Young's modulus, Poisson's ratio and failure strength limits have to be calibrated in regard of their microscopic quantities thanks to virtual tensile tests. These virtual tests involve cubic discrete domains similar to the presented one on Figure 8.

Virtual tensile test and computational methods of macrosopic quantities 5.1.1. Imposed displacements

To simulate a tensile test, opposite displacements are imposed to the discrete elements S x+ and S x-. These displacements are applied along the X axis. A small constant velocity of 20 nanometers per iteration is used to apply these displacements.

Numerical measurement of the main normal force F

The cohesive normal forces F x-and F x+ resulting from the imposed displacements are considered as the sum of the local forces f i applied by the beams to the S x+ and S x-discrete element sets :

F x-= i∈S x- f i • x (25)
F x+ = i∈S x+ f i • x ( 26
)
If no failure occurs during the simulation, the cohesive normal forces F x-and F x+ are opposite :

F x-≈ -F x+ (27)
The cohesive normal force F is computed by averaging the two opposite forces :

F = 1 2 (F x--F x+ ) (28)

Numerical measurement of the main normal stress σ Mxx

The main normal stress σ Mxx is simply computed as :

σ Mxx = F • x L y × L z ( 29
)
where L y and L z are the lengths of the sample along the y and z axes.

Numerical measurement of the main strains

The macroscopic strain ε Mxx , ε Myy and ε Mzz along the x, y and z axes can be expressed as :

ε Mxx = ∆L x L x ε Myy = ∆L y L y ε Mzz = ∆L z L z (30)

Numerical measurement of the macroscopic Young's modulus

The macroscopic Young's modulus E M is computed as :

E M = σ Mxx ε Mxx (31)

Numerical measurement of the macroscopic Poisson's ratio

The macroscopic Poisson's ratio ν M is computed by averaging the macroscopic Poisson's ratio values along the y and z axes as :

ν My = - ε Myy ε Mxx ν Mz = - ε Mzz ε Mxx (32) ν M = 1 2 ν My + ν Mz (33)
where ν My ≈ ν Mz .

Numerical measurement of the macroscopic failure stress

The macroscopic failure stress σ M f is the value of the macroscopic tensile stress σ Mxx when the failure occurs. The failure is considered when a sudden decrease of the measured macroscopic tensile stress σ M appears (see Figure 10-(c)).

Typical results from virtual tensile tests

Typical behaviors resulting from virtual tensile tests are reported on Figure 10. The (a) and (b) graphs show that the macroscopic Poisson's ratio ν M and the macroscopic Young's modulus E M quickly converge towards constant values : 0.23 and 72 GPa respectively . These converged values are considered as the material properties of the domain in regard of their microscopic quantities. The (c) graph shows the evolution of the macroscopic tensile stress σ Mxx versus the macroscopic strain ε Mxx . A sudden decreasing of the σ Mxx indicates the failure of the virtual specimen. The corresponding value of σ Mxx is considered as the material strength limit σ M f .

The results plotted in Figure 10 correspond to the macroscopic properties of the glass matrix. In other words, the microscopic parameters were calibrated to match the glass values.

Calibration procedure

Deducing microscopic parameters values that match wanted values at the macroscopic scale presents some difficulties. This problem has been discussed in detail by Ostoja-Starzewski [START_REF] Ostoja-Starzewski | Lattice models in micromechanics[END_REF]. The author proposes micro-macro analytic laws for typical ordered lattice configurations. In the last paragraph ( § 6.3) dedicated to the periodic random lattice network, Ostoja-Starzewski proposes numerical tests to calibrate the model. By the way, analytic approaches are limited to ordered configurations as mentioned by Potyondy and Cundall in (11, §3.1).

Another difficulty is the relationship between the sample topology (shape, refinement, etc.) and emergent macroscopic behaviors. It induces that the microscopic parameters are only valid for a single and unique discrete element sample. This problem was underlined by Hentz et al. in [START_REF] Hentz | Identification and validation of a discrete element model for concrete[END_REF] and Liao et al. in [START_REF] Liao | Stress-strain relationship for granular materials based on the hypothesis of best fit[END_REF].

A solution was proposed by Andre and al. in [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF] and [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF]. These studies show that the initial packing processes of the discrete samples is of high level of importance. The packing processes must ensure a stable level of geometrical isotropy, volumic fraction and coordination number between the different created samples. These criteria allows to minimize the level of dependency between the sample topology and the emergent behaviors of these samples. The Cooker program, proposed by the GranOO discrete element workbench, is able to control these packing parameters.

To minimize the number of simulations, the calibration involves a well defined procedural approach described in [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF] and [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF]. The resulting set of microscopic parameters is independent of the size, the shape or the refinement of the discrete domain. It means that this set of microscopic values can be applied to any kind of discrete domains with arbitrary shapes. Figure 11 shows the convergence rate of the macroscopic Young's modulus E M , the macroscopic Poisson's ratio ν M and the macroscopic tensile strength σ f M versus the number of discrete elements. The computations were conducted on discrete samples with an increasing number of discrete element from 2,000 to 20,0000 (see Figure 12). The same set of microscopic parameter values was used for all these different samples. Because the initial packing algorithm involves random procedures, four samples were built for each number of discrete element. It explains the scatters at each given point. These curves underline that the macroscopic quantities converge toward single values. It confirms that a given set of microscopic parameter values gives a single set of macroscopic parameters values for any discrete samples. Table 3: The main microscopic thermo-elastic parameter values of glass and alumina after calibration process where E µ is microscopic Young's modulus, R µ is microscopic radius ratio, σ µ f is microscopic failure strength limit and α µ is microscopic linear thermal expansion coefficient

Two-phase materials

The discrete domains used to obtain the virtual two-phase materials were built in a box shape with approximatively 10,000 discrete elements. Four discrete domains were built to ensure the statistical representativeness of the calculus. After building the discrete domains, a second phase (inclusions) is generated thanks to a geometrical algorithm that randomly distributes spherical inclusions inside the cubic domains. If the inclusions intersect the boundaries of the domains, they are cut by the boundary plane of the cubic domain (see Figure 13a).

For each specimen, several volume proportions of inclusions have been chosen : 0%, 5%, 10%, 15%, 20%, 25%, 30%, 35%, 40% and 45%. The average inclusion radius is 250 µm ±10%. The number of discrete elements that belongs to an inclusion is about 200. The virtual cubic two-phase domains with approximatively 2.2 millimeters length are composed of three types of bonds (see Figure 13b): matrix, inclusion and interface bonds. Because no interphase exists between alumina inclusion and glass matrix (see section 1), the bonds that link a matrix discrete element and an inclusion discrete element are considered as matrix bonds.

The lengths of the virtual domains were chosen to reproduce a statically rep-resentative sample. The imposed number of discrete elements (around 10,000) corresponds to a good compromise between the precision of the results and the computational time cost.

Matrix discret elements

Inclusion discret elements

Unbreakable numerical model

As a preliminary approach, the failure stress of each phase is set to an infinite value. So, no micro-crack appears during thermo-elastic numerical test. The domains are submitted to a numerical cooling from 450 to 20 • C to get the CTE of virtual undamaged two-phase domains. Then, tensile tests are performed to measure macroscopic Young's modulus. The results Figure 14 show that the thermo-mechanical behavior of the simulated cohesive two-phase material are consistent with the predictions of H&S model. It allows to validate the numerical two-phase model. However, the experimental results reported in Figure 3 show that this approach is limited to perfect materials.

Breakable numerical model

In this section, the matrix bonds are considered as breakable. Following the approach described in Figure 4, the discrete domains are subjected to cooling tests (where micro-cracks appear), and then, to a tensile test. Here, the virtual tensile tests are non-destructive. The micro-cracks are "frozen" and they are not allowed to propagate themselves during the tensile test. It allows to measure the apparent Young's modulus of the damaged virtual samples. Experimentally, the Young's modulus is measured by ultrasonic pulse echography technique does not provide any additional damage in normal conditions of using [START_REF] Huger | High-temperature measurement of ultrasonic wave velocity in refractory materials[END_REF]. The CTE is measured with a standard dilatometer equipment.

Discussion

Qualitatively, the evolution of numerical Young's modulus with the increase of volume fraction of inclusions follows experimental observations and remains below H&S predictions. Quantitatively, the numerical results overestimate Young's modulus when volume fraction is lower than 20%. However, Young's modulus estimation is in good accordance with experimental data for volume fraction higher than 20%. In addition, the numerical prediction of CTE matches quantitatively the experimental observations. For numerical samples, the CTE are measured indirectly from length variations. Figure 16 shows the evolution of numerical and experimental thermal strains versus temperature during a cooling stage from 450 • C to 50 • C. These data concern an experimental sample and a numerical sample with 45 vol% of inclusions. Figure 16 shows that the numerical results are in good accordance with the experimental observation. The evolution of numerical strains follow two distinct zones named undamaged and damaged. From this observation, two CTE can be computed from linear regressions applied inside the damaged and undamaged data range. It is observed that the undamaged CTE is inside the H&S bounds. The damaged CTE is lower than the lower H&S bound and is very close to the experimental one. Another value of interest is the number of broken beams during a simulation. Figure 16 shows that the strain slope changing corresponds to the beginning of cracking.

This two distinct regimes are not clearly observed on experimental data. This difference is due to the initial level of damage between the experimental sample and the numerical one. At the initial temperature, i.e, 450 • C, the numerical sample is free of damage. Therefore, the experimental sample is already slightly damaged. So the changing strain slope is smoothed and not clearly visible. That's why the experimental CTE values reported in figure 15a are deduced from linear regression on the complete range from (450 • C to 20 • C) while the numerical CTE are deduced from the range 300 • C to 20 • C.

In addition, only the radial micro-cracks are taken into account in the numerical model. However, real microstructure of such two-phase materials include other discontinuous phenomena, such as debounding around inclusions or defects that may influence the thermo-elastic properties.

As explained in section 2, the system ∆α > 0 promotes the radial microcracks. The agreement between experimental and numerical results while the presence of alumina inclusions increases seems to validate that radial micro-cracks are predominant. For a low volume fraction, the number of radial micro-cracks is low. As a consequence, the other types of damages such as debounding or defects may have a non negligible effects on the apparent Young's modulus.

Radial micro-cracking phenomenon seems to monitor elastic behavior of twophase model materials for high level of inclusions. So, according to the results in Figure 15, the numerical model may be relevant for quantitative studies. The developed model seems to catch and predict the main phenomena that drives the evolution of Young's modulus and CTE of two-phase damaged materials.

Perspectives and conclusions

This work related to the simulation of thermo-elastic behaviour using DEM presents a significant improvement for discrete element methods applied to the simulation of continuum media. The proposed DEM model has been first validated in regards with well established analytic solutions. In the case of undamaged materials, the numerical results show a good level of accordance with H&S prediction. In addition, the developed method allows to predict (i) the damage level generated by thermal cycle on materials that exhibit a thermal expansion mismatch and (ii) the apparent thermo-mechanical parameters such as Young's modulus and CTE. This study confirms that thermal damages provoke a great influence on apparent thermo-mechanical parameters and the developed model may be applied to a large class of heterogeneous brittle materials : geo-materials, civil engineering materials and ceramics.

The prediction of apparent thermo-mechanical parameters such as Young's modulus, Poisson's ratio or CTE is of high level of importance for engineering applications. The proposed method seems to be adapted to predict this parameters and allows to consider further studies to improve the understanding of industrial multiphase materials. Further investigations may address micro-cracks closure, inclusion debounding or self-healing (that appears when the temperature is greater than the glass transition temperature). The main perspective of this study is the

Figure 1 :

 1 Figure 1: Microstructure of a generated cracking network

Figure 2 :

 2 Figure 2: Sketch of the residual stresses and induced damages during cooling stage

  Thermal expansion coefficient α v

Figure 3 :

 3 Figure 3: Comparison between experimental thermo-mechanical properties and analytic prediction at room temperature

1 - 2 -Figure 4 :

 124 Figure 4: General view of the whole simulation steps

Figure 5 :

 5 Figure 5: The cohesive beam bond and the two scales involved by the DEM

Figure 6 :Figure 7 := l initial 0 +

 670 Figure 6: Illustration of the method for selecting broken bonds

Figure 8 :

 8 Figure 8: The cubic discrete domain used to study thermal expansion

- 6 Figure 9 :

 69 Figure 9: Evolution of the macroscopic linear thermal expansion coefficients α Mx , α My and α Mz versus the number of iterations

Figure 10 :

 10 Figure 10: Evolution of (a) the macroscopic Poisson's ratio ν M , (b) the macroscopic Young's modulus E M and (c) the macroscopic tensile stresses σ M versus the macroscopic strain ε Mxx during a virtual tensile test

Figure 11 :

 11 Figure 11: Evolution of (a) the macroscopic (a) Young's modulus E M , (b) the macroscopic Poisson's ratio ν M and (c) the macroscopic tensile strength σ fM versus the number of the discrete element

(a) 2 ,

 2 Figure 12: Snapshot of discrete samples with increasing refinementE µ (GPa) R µ σ µ f (MPa) α µ (10-6 K -1 ) glass 5,200 0.1 60 11.6 alumina 24,560 0.1 360 7.6

Figure 13 :

 13 Figure 13: Virtual sample of the two-phase model material

Figure 14 :Figure 15 :

 1415 Figure 14: Comparison between numerical and analytical H&S models of free of damage twophase materials

Figure 16 :

 16 Figure16: Evolution of (i) numerical and experimental thermal strains and (ii) number of broken beams versus temperature for a sample with 45% of inclusions

Table 1 :

 1 Matrix and inclusions thermo-mechanical properties

	Matrix	Inclusion

  • E i and ν i are respectively the Young's modulus and the Poisson's ratio.

	α matrix > α inclusion and T •
	σ rad
	σ orth
	α inclusion
	σ orth
	Inclusion
	σ rad
	α matrix
	Matrix

Table 2 :

 2 The main thermo-mechanical parameter values of the borosilicate glass and alumina where E is the Young's modulus, ν is the Poisson's ratio, σ f is the failure strength limit and α is the linear thermal expansion coefficient

	-6 K -1 )

see www.granoo.org