A discrete element model of brittle damages generated by thermal expansion mismatch of heterogeneous media

. At the macroscopic scale, such media as rocks or ceramics can be seen as homogeneous continuum. However, at the microscopic scale these materials involve sophisticated micro-structures that mix several phases. Generally, these micro-structures are composed by a large amount of inclusions embedded in a brittle matrix that ensures the cohesion of the structure. These materials generally exhibit complex non linear mechanical behaviors that result from the interactions between the different phases. This paper proposes to study the impact of the diffuse damages that result from the thermal expansion mismatch between the phases in presence. The Discrete Element Method (DEM) that naturally take into account discontinuities is proposed to study these phenomena.


Introduction
Rocks, ceramics or refractories are heterogeneous materials that exhibit a multi-phase composition and involve different sizes of aggregates, bonding phases and various additives. The description and the prediction of thermomechanical behavior of such materials present a real difficulty due to their complex microstructure. Hashin & Shtrikman (H&S) have developed an analytical approach that predicts thermo-mechanical properties of perfectly cohesive multi-phase materials free of damages. However, most of real materials present numerous micro-cracks at room temperature resulting from thermal expansion mismatches between their phases. These damages highly influence the thermo-mechanical properties of such materials.
For instance, Young's modulus is strongly affected by the presence of micro-cracks and the measured values are often in disagreement with H&S's prediction [1]. The study proposed here proposes a numerical method in order to predict the occurrences of these damages and their influences on macroscopic properties such as Young's modulus. The Finite Element Method (FEM) is often used to study mechanical properties of materials. However, at the microscopic scale, this method is not adapted to describe discontinuities (micro-cracks) without assumptions on their localizations, their paths, their growths and their coalescences.
The DEM could be an interesting alternative to study multi-damaged materials because it takes naturally into account these discontinuities. The DEM implements a group of distinct elements (also named discrete element) that are in interaction though contacts or cohesive laws. This model consists of an assembly of discrete elements, deformable or not, linked up by simple mechanical laws to mimic the behavior of the material. The advantages are the description, in a natural way, of the crack initiations, the crack propagations and coalescences. Researchers have used this method to study damages in solids, such as rocks or ceramics. However, the main difficulty is to tackle quantitative results due to the necessity to find the relations between unknown microscopic laws, at the discrete element scale, and the macroscopic properties, at the structure scale, such as Young's modulus, Poisson's ratio or Coefficient of Thermal Expansion (CTE).

Model material
The reference 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.

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 micro-crack network. The main thermo-mechanical properties of both individual materials are given in the table 1. For details about the experimental processing of model materials, readers may refer to [1,2].

Induced residual thermal stresses
The introduction of spherical particles in the matrix leads the occurring of thermal stresses during the cooling stage of the sample processing. Because of the CTE mismatch, the matrix is under tensile mode and the inclusions are subjected to compressive stresses. More precisely, The matrix is subjected to tensile orthoradial stresses σ orth and radial compression σ rad . Brittleness of the glass matrix induces orthoradial cracks that occurs and propagates in the matrix. Figure 1 shows the resulting micro-crack network with several inclusions on a model material.

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 because this analytical approach is well suited to describe inclusional materials. The apparent Young's modulus can be predicted from the For details about this analytical approach, reader may refer to [3]. As shown in Figure 2, the plots of the experimental Young's modulus (noted E) are lower than H&S lower bound (HS − ) prediction. The Young's modulus was measured through an ultrasonic technique [4]. 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 rigidity 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.
The following parts of this paper will introduce the discrete element model used to predict apparent Young's modulus that can not be treated with standard H&S models and FEM models.

The discrete element model
The discrete element approach used here is a mix between the lattice models and the particle models as it was first proposed by POTYONDY in [5]. Unlike continuous approaches, the main difficulty for DEM is to simulate quan-  Figure 3 shows the main simulation steps implemented in this study. The first step consists in calibrating the thermomechanical parameters of the discrete element model (for details, see the next sections). 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 limit σ μ f and CTE α μ . Considering the values reported in Table 1, both borosilicate glass and alumina parameters are calibrated separately. After building the initial cubic discrete domains, the second step consists in inserting the spherical alumina inclusions. This step involves a simple geometrical algorithm. Then, the virtual sample is cooled from 450 • C (the glass transition temperature) to 20 • C (room temperature). This step, that involves thermo-mechanical simulations, leads to cracks initiation and propagation in the virtual samples (see Figure 4). In this study, the temperature is supposed constant within the sample. As a first approach, thermal conduction is neglected. Finally, the damaged virtual samples given by the last step are frozen. In this case, only the elastic behavior is taking into account in order to forbid crack extensions. These frozen samples are submitted to virtual tensile tests to evaluate their apparent Young's modulii E M .

Elastic behavior
The cohesive beam bond model [7] 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 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. 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 of the simulated material.

Brittle behavior
Within the discrete element approach, cracks are generally simulated by breaking the cohesive bonds if a criterion is reached. The main existing approaches are based on the computation of bond strains or stresses [5]. However, these approaches are not able to simulate complex crack paths such as the hertzian cone crack that exhibits on fused silica glass [8]. 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. To compute this hydrostatic stress, an equivalent Cauchy stress tensor for each discrete element is computed as : 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 i j 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 i j is the relative position vector between the center of the two bonded discrete elements i and j. Finally, the developed criterion assumes that fracture occurs if the computed tensile hydrostatic stress is higher than a threshold value σ f :

Thermal expansion behavior
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 macroscopic thermal expansion of the discrete domain. To study the relation between the microscopic and the macroscopic thermal expansion, some simulations were carried out. These simulations conduct to the following simple result : So, the thermal expansion coefficient can be introduced directly without any calibration.

Results and conclusion
Following the approach described in the last sections, discrete domains were subjected to cooling tests and then, to tensile tests in order to quantify their apparent Young's modulii. During the tensile tests, the micro-cracks are "frozen" and are not allowed to extend. Experimentally, the Young's modulii are measured by ultrasonic pulse echography technique that do not damage the materials. Finally, the results of numerical simulations and experimental observations are plotted in Figure 5, where : the absciss Volume fraction is the volumic fraction of inclusion in samples, the Hv-and Hv+ curves denotes bounds of the H&S predictions, the Numerical 3D model without cracks curve corresponds to the undamaged DEM samples, the Numerical 3D model with cracks curve corresponds to the damaged DEM samples and the Experimental curve corresponds to the experimental observations. These results shows that the Numerical 3D model without cracks match perfectly with the H&S predictions. This result confirm that the H&S predictions are not able to take into account the damages inside the material. However, the Numerical 3D model with cracks match with the experimental results measured by ultrasonic pulse echography technique. It allows to validate quantitatively the numerical 3D DEM model. This work related to the simulation of thermo-elastic behaviour using DEM presents a significant improvement of the discrete element methods applied to the simulation of continuum media. This model has been applied to study the influence of the damages generated during thermal cooling of multiphase materials. These damages highly influence the rigidity of materials and is of high level of importance for common engineering applications. The proposed method seems to be adapted to predict this rigidity and allows to consider further studies to improve the understanding of more complex multiphase materials.