Numerical Modelling of The Quasi-Brittle Behaviour of Refractory Ceramics by Considering Microcracks Effect

In the steelmaking industry, the inner lining of ladles is made of refractory ceramics, which are constantly subjected to thermal shocks during their service. Experimentally, it is observed that pre-existing microcracks could significantly increase the thermal shock resistance of these ceramics. The presence of such microcracks network within the refractory microstructure could lead to a non-linear quasi-brittle mechanical behaviour. To model this quasi-brittle behaviour, a suitable numerical approach is the Discrete Element Method (DEM), which can circumvent the limitations of more conventional continuum approaches in capturing microstructural effects required to simulate multi-fracture propagation. Here, it is aimed to simulate such quasi-brittle behaviour by initial well-distributed damages, with a strength dispersion following a Weibull distribution. In this way, the microcracks effect on the quasi-brittle behaviour of a numerical sample under uniaxial and cyclic tensile tests is investigated. Ultimately, a quantitative DEM model to simulate such a complex behaviour is


Introduction, motivation of the study
The ceramic materials currently used as refractory linings of industrial vessels required to produce steel, glass, and cement are regularly subjected to thermal shocks and cycles. For example, during the steelmaking process, due to the filling and emptying of melted iron in the steel ladle, refractories face severe thermal cycles. Therefore, these refractories should be resistant to such thermal loadings. A high number of research works has pointed out that the microstructure design of refractories has a significant influence on their ability to sustain thermal shocks during their service time [1][2] [3]. However, understanding this phenomenon is not an easy task due to the complex microstructures of these heterogeneous ceramic materials. Industrial and academic studies have then already been managed in order to better understand how thermal shock resistance can be improved [4] [5]. Such an improved resistance was first correlated to an increase in mechanical toughness. Then, it became clear that the thermal expansion mismatch between constituents (large aggregates within a matrix), entailing microcracking around the aggregates, was mainly responsible for an improved thermal shock resistance [6]. Indeed, the presence of these microcracks will induce a non-linear quasi-brittle mechanical behaviour [7].
It is thus clear that numerical modelling can also help to build a better understanding of such complex quasi-brittle behaviour resulting from this extensive microcrack network.
For such simulations, it should be highlighted that due to a very high number of evolving discontinuities, it could be challenging to anticipate the correct value for the elastic parameters of the material by using Finite Element Method (FEM) simulations [8]. To solve this problem and to introduce arbitrary crack propagation, the Extended Finite Element Method (XFEM) was introduced by Belytschko et al. and Moës et al. in 1999 [9] [10] [11]. Many studies about microcracks and crack propagations using XFEM can be cited. Souza et al. have developed a multiscale model for the transition of local cracks to the macrocracks by using XFEM to model cohesive zones [12]. As another example, Goyal et al. used XFEM to model crack propagations in an aluminium alloy in a uniaxial tensile test [13]. However, XFEM is not able to manage multiple crack initiation and propagation at the same time.
For this purpose, the Discrete Element Method (DEM) is an adequate alternative numerical method as it can simultaneously manage initiation and propagation of many cracks. In fact, using DEM to model fracture mechanisms in such a pseudo-continuum media, especially in the ceramics field, is a new and promising approach under active development .
Based on a study by Sharafisafa et al. [14], for simulation of crack propagation in the brittle materials, DEM could simulate realistic secondary crack branching (propagating from the main crack path). In contrast, the secondary cracks were not observed in XFEM models. Furthermore, DEM could simulate the entire process of crack initiation, propagation, and, more importantly, the coalescence of the fracture. In another study by Hedjazi et al. [15], for simulating crack branching in a dense vitreous material, it was shown that DEM results were closer to the analytical models in comparison to FEM, and it was more in accordance with experimental visual observation for crack branching in certain conditions.
Simulating quantitatively a continuum media with DEM is not as usual as with FEM because the continuum mechanics cannot be introduced directly into the DEM models. On the other hand, as mentioned, simulating discontinuities such as fractures or damages are more straightforward to be reproduced in DEM frameworks because of the discrete nature of discontinuous phenomena [16] [17].
The present study investigated the randomisation of local fracture criteria within a virtual sample using a Weibull distribution. The key idea is here to verify if the non-linear mechanical behaviour of a refractory ceramic can be reproduced by initial well-distributed damages with a strength dispersion following a Weibull distribution. This approach could allow studying the effect of pre-existing microcracks on the non-linear quasi-brittle behaviour of a numerical sample (representing a refractory ceramic) under uniaxial and cyclic tensile tests. Ultimately, this proposed approach will lead to a qualitative and quantitative DEM model to simulate the non-linear quasi-brittle behaviour of refractory ceramics.
In this paper, prior to focus on the numerical aspect, the reference material, Alumina Spinel brick classically used as a refractory lining of steel ladle [18], and the technics to monitor its microcracking during the thermal cycles will be presented. Afterwards, the numerical foundation for this study will be introduced. Later, the randomisation process of bond strengths will be described. Then, the influence of this randomisation process on the apparent mechanical behaviour and microcracking will be studied. Finally, a quantitative DEM numerical approach will be proposed to model the non-linear quasi-brittle behaviour.

Description of the reference material and monitoring the microcracking process
The chosen reference material for this study is an Alumina Spinel brick because it promotes a strong non-linear quasi-brittle behaviour that suits this paper's overall goal. This particular material is commonly used in the working lining of steel ladle. The purpose of this short section is to illustrate (by a couple of key experimental results) where the microcracks (which constitute the central point of the present study) come from and how these microcracks lead to non-linear quasi-brittle behaviour.
Kaczmarek et al. has investigated this material [19]. This Alumina Spinel brick is composed of tabular Alumina grains (up to ~3 mm), Alumina Spinel matrix and 19.7 % of porosity. A typical microstructure of this Alumina Spinel is shown in Fig. 1.  [19] The influence of the microcracking process on the elastic properties can be examined by using the ultrasonic method and acoustic emission. The ultrasonic method was used to detect and record the evolution of Young's modulus versus temperature. This low-frequency method uses a long bar configuration to perform measurements at high temperature. This method has high precision in measuring the impact of microcracks on the macroscopic elastic properties.
The acoustic emission method records and analyses acoustic waves coming from the sample. An acoustic emission event results from the propagation of an elastic wave coming from a localised source within the sample. Each event is counted, leading to the results as accumulative numbers of events versus temperature [20].
Linking the ultrasonic measurements and acoustic emission events could help to better understand the relationships between elastic properties and microcracking process. In this way, Young's modulus and acoustic events evolutions of the referenced Alumina Spinel is plotted in Fig. 2 versus temperature. Just focusing on the cooling part of this curve, Young's modulus increases from 1500°C to 700°C and then decreases from 700°C to ambient temperature (quite abnormal behaviour). In fact, this decrease is due to microcrack initiation and propagation resulting from the Coefficient of Thermal Extension (CTE) mismatch between constituents, which leads to simultaneous increase in acoustic events. Note that although theses cracks are not obviously visible in Fig. 1, experimental results from ultrasound and acoustic emission clearly demonstrate the occurrence of those microcracks during cooling stage of the thermal treatment. Fig. 2. Young's modulus and acoustic events evolutions during thermal treatment [19] Now, to focus on the mechanical behaviour of such an Alumina Spinel, a uniaxial tensile test was performed to characterise this material at ambient temperature. Fig. 3 is showing Alumina Spinel brick experimental results during incremental cyclic mechanical loadings. a. a progressive decrease of its rigidity after each cycle; b. a hysteresis loop between each unloading and reloading; c. a residual strain after each unloading. In order to model this experimental behaviour accurately, a calibration process for simulating with DEM such a non-linear quasi-brittle behaviour should be defined. In the following section, such a quantitative calibration process will be proposed.

Numerical approach for simulating non-linear quasi-brittle behaviour with DEM
In this section, the numerical basis required for simulating non-linear behaviour using Weibull distribution within DEM will be introduced. For the present study, the main DEM framework was the Particle Flow Code (PFC). At first, the calibration of the Flat Joint Model input parameters for a perfect linear elastic brittle material will be presented. Then, the apparent behaviour given by this calibration process will be compared to the experimental uniaxial test of Alumina Spinel. Later, a randomisation process of bond strengths will be introduced following Weibull distributions. Then, the influence of this randomisation process on the apparent mechanical behaviour and microcracking will be studied.

DEM to model microstructure properties relationship using Flat Joint Model (FJM)
DEM numerical method was initially designed to reproduce the mechanical behaviour of intact rocks accounting for their complex microstructure in 1971 [21] [22]. By bonding particles in DEM models, it was possible to simulate pseudo-continuum media, such as intact rocks. DEM has been developed further and became more practical for simulating the fracturing mechanism and crack propagation, especially for brittle materials [23] [24] [25]. Besides, different Bonded Particle Models (BPM) [26] were introduced. However, based on the studies of Cho et al. in 2007, two common bonded models, linear contact bond and parallel bond models, exhibit certain limitations. Namely, the tensile strength to compressive strength ratio was not following the laboratory test results [27]. Hence, in 2012, the Flat Joint Model (FJM) was introduced as one of the BPM models [28], which improved the mentioned limitation. Different studies were done using FJM to show its advantages over the other BPM models, for example, matching the compressive to tensile strength ratio to the real granite sample [29] and the perforation failure in the sandstone [30]. Considering these improvements, the FJM was chosen as the contact model for the present study.
In the FJM model, the behaviour of the bonded model is linear elastic until a strength limit [31]. Table 1 shows the list of FJM independent modifiable parameters that were used in this study. These FJM parameters were applied to the contacts who their distance is smaller than an initial gap ( ). Also, for this study, the normal stiffness and shear stiffness of the contacts were set using local Young's modulus and stiffness ratio.
To be more explicit, the local input parameters of the contact model will be marked by 'loc' in superscript. If the parameters refer to the apparent elastic properties of the specimen, it will be shown by 'ap' in superscript. It is important to note that these different parameters should be calibrated in order to fit material apparent elastic properties such as Young's modulus or Poisson's ratio. Further information about theses parameters are detailed in [28].

Flat Joint Model and the calibration process for perfect brittle linear elastic material
One of the main disadvantages of the DEM approach is the lack of well-defined relationships between the contact properties (the local parameters of the bonds between discrete elements) and the apparent brittle elastic properties, which are related to the simulated material. Hence, the FJM requires a calibration process to assign the correct values for the contact properties to reproduce the desired apparent brittle elastic properties (initial linear part of stressstrain curve). Generally, the calibration process can be done either by a rationalised trial and error [32] or by a direct calibration process [33]. However, due to the high number of local input parameters of the FJM and the complex interdependencies of these parameters on the apparent behaviour, no practical, direct calibration process has been proposed for FJM yet. Nevertheless, some calibration algorithms have been recently proposed by [32] and [34]. From these works and for the general goal of this paper, a modified indirect calibration algorithm is proposed. The proposed algorithm (shown in Fig. 4) considers the different calibrations steps for the most influential FJM local parameters.
The very preliminary step for the calibration is to fix the numerical sample geometry, which here is a cube. Then, the modelled sample should be mechanically loaded in order to obtain the apparent mechanical behaviour targeted for the material. As the present paper is dedicated to refractory materials (which are well-known to be more sensitive in tension, mode I), the calibration process proposed in Fig. 4 is exclusively focused here on tensile behaviour. The proposed calibration of FJM contact model parameters follows a systematic step by step approach where each step is detailed below.
1. Initial bonding gap ( ): assigning a small value for the initial bonding gap will result in an asymmetric rigidity between compression and tension, which is unwanted for the purpose of this study. Thus, here, this initial gap is set to the 40% of the average discrete element's diameters ( ̅ ) to satisfy this condition: Local stiffness ratio ( ): in the second step, the local stiffness ratio ( ) is calibrated. As mentioned earlier, this parameter greatly influences the apparent Poisson's ratio (ν ap ). Hence, the local stiffness ratio ( ) should be set to a value that produces the targeted apparent Poisson's ratio (ν ap ). It should be mentioned that this step must be done before calibrating the local Young's modulus ( ) as the apparent Young's modulus also depends on the parameter [34].

Local Young's modulus (
): after calibrating and fixing the local stiffness ratio ( ), the local Young's modulus ( ) should be set to a value that produces the targeted apparent Young's modulus ( ) [34]. 4. Local tensile strength ( ): finally, after calibrating the elastic properties, the targeted apparent tensile strength ( ) is calibrated by changing the local tensile strength ( ). The starting guess for can be the value itself. Generally, the apparent tensile strength of the material is not so far from the local tensile strength. It is essential to highlight that the proposed calibration algorithm is for perfect brittle linear elastic materials. Therefore, in the case of Alumina Spinel (which is a quasi-brittle ceramic), these calibration steps were used for the initial linear part of the tensile behaviour of the ceramic, which is measured at the very beginning of the stress-strain curve (Fig. 3).

Uniaxial tensile test simulation with non-randomised (uniform) strength of bonds
As described in the previous section, the calibration algorithm for a perfect brittle linear elastic material was used in order to achieve apparent Young's modulus ( ), apparent Poisson's ratio (ν ap ) and apparent direct tensile strength ( ) of the ceramic material. It should be emphasised that, at this stage, the has a uniform single value of 1.33 MPa for all FJM bonds. Then, by using these FJM input parameters, a uniaxial tensile test was simulated and compared to the experimental Alumina Spinel results in Fig. 6 As shown in Fig. 6, this preliminary calibration step leads to a perfect brittle elastic behaviour that targets the proper value of the macroscopic failure strength (1.58 MPa). At the beginning of the stress-strain curve, the simulated curve tangent the experimental one, which means that the initial Young's modulus value was successfully calibrated. It should be mentioned that the initial oscillation at the very beginning of the stress-strain curve is due to applying a small displacement rate (to model a direct tensile test) and changing the condition of the sample from the fully static to the quasi-static. In this way the discrete elements at the loading boundaries (see Fig. 5) would have vibrations, leading to the mentioned oscillation. This is an artefact due to the initial sudden acceleration of the system at the very starting point of loading. These oscillations would not affect the overall output and of such mechanical simulations.
As expected, the simulation is not showing any quasi-brittle behaviour since no parameter was introduced into the model for that purpose. Also, at this stage, the simulation simply exhibits an abrupt failure that is not in line with the real material behaviour. Therefore, it is necessary to introduce an additional technique capable of reproducing this quasi-brittle behaviour. For such a purpose, a random distribution of local tensile strength is introduced to model such non-linear behaviour quantitatively. In the following parts, a dedicated quantitative calibration procedure will be proposed and explained step by step.

Introduction of a Weibull distribution within bond strengths
By using Weibull distribution, the probability of failure for a solid material under a tensile load Pf (σ) is given as [35]: where σ, m and σ0 are the tensile stress (MPa), the Weibull modulus, and the scaling parameter, respectively. 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 The Weibull statistical distribution is commonly used as a macroscopic concept to describe the occurrence of random defects within a 3D solid sample, which could cause a variation in the macroscopic tensile strength of different specimens for a given material. As a material parameter, the Weibull modulus (m) could be experimentally obtained from different mechanical tests (such as the uniaxial tensile test or 3-point bending test) conducted on the same sample geometries.
In the present study, this concept is used in the context of refractory ceramics under tension by varying the local bond tensile strength ( ). Hence, it is proposed to reproduce the microscopic defects (coming from Coefficient of Thermal Expansion (CTE) mismatch during cooling after sintering) within the DEM numerical sample by randomising the local bond tensile strength ( ). In such a perspective, the Weibull distribution law is chosen here to set randomly values, following Eq. 2, where σ corresponds thus to local bond tensile strength values ( ). Moreover, the scaling parameter 0 for these local tensile strengths is set (at this stage) to the value previously determined without any Weibull distribution (corresponding in fact to m = ∞ and thus will be noted later as =∞ 0 = =∞ ). In this way, the impact of this local randomisation will be investigated on the macroscopic non-linear behaviour of the Alumina Spinel refractory under tensile loadings. In reality, this non-linear effect comes from defects originated by thermal expansion mismatch between constituents within the microstructure.

Influence of Weibull distribution on the mechanical behaviour and microcracking process during loading
The impact of applying randomisation on the local bond tensile strength ( ) is studied on Alumina Spinel ceramic as the reference material. As mentioned, Alumina Spinel exhibits a quasi-brittle behaviour with an apparent tensile strength of 1.57 MPa, which corresponds to the  These samples were submitted to a uniaxial tensile test in order to study the influence of Weibull modulus on the stress-strain curve. The results, which include more cases for Weibull modulus, are shown in Fig. 8 (a). This figure shows that using Weibull distribution to randomise has a strong influence on the mechanical behaviour of the numerical samples. Low values of Weibull modulus (m) promote non-linearity in the behaviour of the sample during tension. During the progressive loading of the samples, the weaker bonds created by the randomisation process breaks earlier than stronger ones, as expected. These weaker bonds can be considered as resulting from micro defects. While decreasing Weibull modulus (m), it increases the number of these micro defects that induces a progressive decrease in the apparent rigidity during loading (non-linear mechanical response). In addition to the stress-strain curves, Fig. 8 (a), the number of detached FJM elements related to microcracks is shown in Fig. 8 (b) for different m values. It highlights that by decreasing m: (a) the microcracking process starts for lower strains; (b) the complete failure is achieved with a higher number of microcracks.  It could be noticed that a significant change in the apparent mechanical response and the progression of microcracking can be observed for the lowest m value (m = 1.75). Even if no explanation can be given at this stage of the study, it probably indicates that it could be interesting to investigate further this range of low values of Weibull modulus ( = [0. 5,2]). This point could be particularly important because the fracture energy (surface below the green curve of m = 1.75 in Fig. 8 (a)) seems to become larger Fig. 9 shows the position of microcracks within the numerical sample during these direct tensile tests for three different Weibull moduli (m = ∞, 5 and 1.75) at the same loading state corresponding to 40% of maximal stress after each peak (post-peak region). Note that the complete detachment of an FJM contact is considered one microcrack for plotting these figures (shown in yellow). On these figures, the red dashed areas show the macroscopic failures 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 resulting from the localisation of diffuse microcracks. The sample with a uniform value of 1.33 MPa for the distribution (where m = ∞, (a)) is showing a single fracture, which leads to an abrupt failure of the sample. While the sample with low Weibull modulus (m = 1.75, (c)) shows diffuse microcracking, resulting in quasi-brittle behaviour. Hence, lowering the Weibull modulus (m) promotes non-linear mechanical behaviour by increasing the occurrence of diffuse microcracks. This fracturing process is quite similar to the phenomena that explain the quasi-brittle behaviour of certain refractory materials [20]. In addition to the influence on non-linear mechanical behaviour, the Weibull modulus also strongly influences the ultimate apparent tensile strength (as shown in Fig. 8 (a)). By decreasing m, the ultimate apparent tensile strength decreases. Considering these points, in order to achieve quantitative simulations of the nonlinear behaviour of such real material, a robust calibration method is proposed in the next section.

Quantifying experimental non-linear behaviour
The first step to calibrate the Alumina Spinel reference material is to quantify its non-linear behaviour. In this regard, the secant modulus ( ) at the peak have been used, which corresponds to the slope between the origin and the stress-strain curve for maximum stress, described as follows: where, and are the stress and strain at the peak of the curve, respectively. In fact, using the secant modulus for non-linear materials, such as soils, especially under cyclic tests, is common in the field of civil engineering for characterising and predicting the behaviour of such materials [39].
By calculating this secant modulus ( ) and having the initial Young's modulus ( ), which has been measured at the early stages of loading, it is possible to quantify Young's modulus decreasing rate. These two parameters were calculated for Alumina Spinel reference material, as shown in Fig. 10. As shown in Fig. 10, the initial Young's modulus is 10.37 GPa. After subjecting to six cycles of tensile loadings, the secant modulus of the peak was measured equal to 6.47 GPa. Instead of using absolute values, normalised values of the decreasing Young's modulus rate 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 are used to be more general and case-independent. Hence, this Young's modulus normalised decreasing rate is defined as: In the case of Alumina Spinel brick, this normalised decreasing rate was 62.4%. From this experimental rate, the way to mimic this decrease quantitatively in DEM will be explained in the next section.

Mimicking the apparent non-linear behaviour
As already presented in Fig. 8 (a), the Weibull modulus (m) strongly influences the non-linear behaviour in tension. In order to go a step further, it could be valuable to find a reliable relationship between m and this non-linear behaviour through the secant modulus, as defined in the previous section. In this regard, eight uniaxial tensile test simulations on a cubic sample with different values of Weibull modulus (m = 1, 2, 3, 5, 7, 10, 20, 30) have been performed. Therefore, all initial elastic parameters (Young's modulus and Poisson's ratio) for this series of simulations were constant. Again, to remain independent from initial elastic values, the results were normalised and plotted in Fig. 11. For a better understanding of this curve (Fig. 11), value close to 100% means that the numerical sample exhibits minimal damage before its complete rupture. On the opposite, value very far from 100% (e.g. 35%) indicates that a large quantity of damage occurs within the numerical sample before its complete rupture (in the pre-peak region).
Thus, very limited damage is observed for the infinite value of m, and the final rupture takes place in a brittle way. This brittle behaviour remains similar by lowering m down to m = 8. For further decrease of m, a larger quantity of damage can be observed before the final rupture, which becomes, in these cases, quasi-brittle. From this curve, which can be used as a calibration chart, it becomes straightforward to set the m value (m ≅ 1.78) in accordance with the targeted value of ( = 62.4%) in the case of Alumina Spinel reference material. In such a way, as shown in Fig. 12, the result of the simulation for this value of m (m = 1.78) perfectly fit the experimental results for both initial Young's modulus and peak secant modulus . Of course, at this stage, the ultimate tensile strength still needs to be adjusted. Therefore, the model should go through another calibration step to reproduce the same ultimate apparent tensile strength of the Alumina Spinel reference material.

Calibrating ultimate apparent tensile strength
In Fig. 12, it was shown that despite the accordance of the non-linear behaviour, the ultimate apparent tensile strength of the simulation was not equal to the experimental result. The reason is that considering Fig. 8 (a), decreasing the Weibull modulus m will decrease, at the same time, the ultimate apparent tensile strength ( ) of the numerical model. To quantify the impact of m on , eight uniaxial tensile test simulations on a cubic sample, with different values of Weibull modulus (m = 1, 2, 3, 5, 7, 10, 20, 30), have been performed and plotted in Fig. 13 (a). Instead of using absolute values of , the normalised ultimate strength values ( ) have been again preferred, which is defined as follows :   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62   Therefore, globally, a decrease in m allows, on the one hand, to generate a non-linear behaviour in the stress-strain curve, but on the other hand, it also leads to a strong decrease in the apparent tensile strength value. Thus, this decrease in apparent tensile strength should be compensated by an increase in the scaling parameter of the Weibull formulation. In this way, the inverse of has been plotted in Fig. 13 (b) in order to obtain a pertinent correction factor that need to be applied to the scaling parameter of the Weibull formulation initially set to =∞ 0 = =∞ . The objective here is to recover the same apparent tensile strength as =∞ by applying the following scaling parameter: In the present case of Alumina Spinel reference material, for which m = 1.78, Fig. 13 (b) allows estimating this correction factor to 174% leading to a scaling parameter ( 0 ) of 2.3 MPa. In order to validate the proposed approach, the numerical results obtained by the calibrated parameters have been compared to the real experimental results for Alumina Spinel brick in Fig. 14. Concerning Fig. 14, after calibrating the non-linear behaviour, the simulation shows a perfect qualitative and quantitative accordance with the pushovers of the experimental results of Alumina Spinel brick. The stress-strain curve of the simulation starts with the same initial Young's modulus. Afterwards, it promotes non-linear behaviour, following the pushovers of the cyclic loadings of the experimental curves. Then, it fails at the same ultimate tensile strength as the real material, and in the end, it potentially predicts a reasonable post-peak failure behaviour of the material.
Overall, it can be said that the proposed procedure for calibrating the DEM model to reproduce the non-linear behaviour of the quasi-brittle materials is validated by comparing it to the real experimental results. This procedure has been summarised in the next section.

Influence of cyclic loadings on non-linearity
In this section, the simulation of cyclic loadings will be reviewed. Firstly, the non-linear stress-strain curve during cyclic loading will be presented. Then the microcracking and elastic properties evolution during cyclic loading will be investigated. Finally, at the end of this section, the cyclic and monotonic uniaxial tensile test simulations will be compared.

Modelling non-linear stress-strain curve during cyclic loading
For verification of the proposed DEM model accuracy, cyclic tensile loading tests are simulated with the same conditions, as explained in section 4. To fit the experimental curve (Fig. 3), five uniaxial tensile cycles were simulated, which target the experimental cyclic peaks at 0.8 MPa, 1.1 MPa, 1.3 MPa, 1.4 MPa and 1.6 MPa (where material failed at 1.57 MPa). In this regard, a series of stress-controlled cyclic tests were done by targeting these five experimental stress peaks. In this simulation, the vertical stresses were gradually applied until it reaches each targeted peak thanks to displacement-driven conditions. After reaching the targeted peak, the displacement directions were inversed to decrease the stress down to the relaxing state equivalent to zero stress. These steps were repeated for each loading cycles up to the failure point. The result of this simulation is quantitatively compared to the experimental test in Fig. 16 .   1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62   By comparing the simulation results and the experimental curve in Fig. 16, it can be observed that the simulation is able to model the progressive decrease of rigidity after each cycle. The overall mechanical behaviour of simulated material and peak of each cycle are quantitatively close to the real material behaviour. This approach could also precisely model the failure peak (+ 0.9% error) of the material after undergoing five cycles. The post-failure behaviour of the material is also modelled reasonably.
On the other hand, the residual strain after each unloading cannot be modelled by this approach. In reality, for the Alumina Spinel ceramic, these residual strains are due to internal stress relaxations, which happen during microcrack extensions. These internal stresses are coming from the cooling phase (after sintering the material), and they are due to the CTE mismatch between the different constituents of the material. At this stage, the proposed DEM model is not able to simulate these internal stresses coming from CTE mismatch and, consequently, the related apparent residual strains. This limitation is because, up to now, the FJM cannot simultaneously model mechanical interactions and temperature variation. However, considering the ongoing development of FJM, this limitation could be solved soon.
The other reason for not mimicking such residual strains could be the perfect fracture closure in the numerical model. This means that after creating a fracture surface during the loading stage, the fracture will close perfectly during unloading without exhibiting any surface mismatch. To simulate a certain level of surface mismatch, an additional artificial gap could potentially be introduced. However, this potential solution was not investigated as a part of this paper.

Microcracking and elastic properties evolution during cyclic loading
To investigate further this cyclic simulation, the stress, the microcracks number and the elastic properties evolutions are plotted versus the simulation time in Fig. 17. To analyse a given cycle, all four plots and the corresponding pictures showing the positioning of cracks should be considered.
Each given cycle starts with a progressive increase of the uniaxial tensile stress up to the targeted peak stress. This progressive increase can be described in two stages: (I) from starting point up to the previous cycle stress peak value corresponding to points a, b, c, and d. (II) from previous cycle stress peak value to the targeted peak corresponding to points 2, 3, 4 and 5. It is worth mentioning that during stage (I), the number of cracks, apparent Young's modulus and Poisson's ratio are constant. While during stage (II), new microcracks are formed, which lead to a non-reversible linear decrease of the apparent Young's modulus. Besides, the slopes of microcrack number (versus time) during the (II) stages are quite similar. It means that the microcracking forming rate is relatively constant. This point is due to the constant displacement rate, which is applied for loading the sample.
After reaching the targeted stress peak, the loading direction is inverted, and the stress decreases to the relaxation state (zero stress). Despite a numerical artefact in elastic properties measurements (which is due to very low values of strain at this point of the curve), the number of microcracks, apparent Young's modulus and Poisson's ratio are quite constant during this stage (I).
For the last cycle, the targeted stress peak was 1.6 MPa, but the sample failed at 1.58 MPa (point 5 to 6) with a sudden drop in the stress. During the failure, the number of microcracks suddenly increases with a slope higher than the observed ones during the previous stages (II). The rising of microcracks number in a short time results in a sharper decrease in apparent Young's modulus.
Finally, in the post-failure part (point 6 to the end of the simulation), the stress tends to zero. During this stage, the microcrack forming rate decreases, and the number of microcracks tends to be constant because complete failure happens. Thus, the apparent Young's modulus and the apparent Poisson's ratio turn to zero in this part.
To get a qualitative visualisation of the microcracking process, pictures of the sample with highlighted cracks are shown for each cycle peak on the top of the figure. These pictures highlight the microcrack progression for each cycle due to the imposed increasing stress after each cycle. In the last cycle, these microcracks merge into a macroscopic fracture (shown in the red dashed area), leading to the complete failure of the sample. This numerical result is in perfect accordance with the localisation phenomenon, which typically is observed experimentally. As expected, this increase of microcracks leads to a progressive decrease in apparent Young's modulus. Also, a noticeable decrease in Poisson's ratio is observed (from 0.16 to 0.13). In fact, this observation is also in accordance with experimental measurements, which in many cases indicate a low value of Poisson's ratio for highly damaged materials. Overall, these different points highlight the pertinence of the proposed model compared to physical phenomena observed experimentally in quasi-brittle ceramics.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 The number of microcracks in simulation could qualitatively be associated with acoustic emission events. To verify this point, the presented simulation is qualitatively compared to experimental observations from [40] of a cyclic tensile test on a quasi-brittle material (marble rock) in Fig. 18. In experimental results, Fig. 18 (a), the drawn KE points refer to Kaiser Effect acoustic emission points. The Kaiser effect is defined as the absence of detectable acoustic emissions until the previous stress peak in a cyclic mechanical test was reached [41]. In fact, the Kaiser acoustic effect in tensile cyclic loading indicates damage accumulation in material under tensile stress [40] and can be used as another verification parameter for the proposed numerical approach.
It can be observed that the evolution of the number of microcracks in this simulation qualitatively resembles the experimental acoustic emission events of a similar experimental test. Moreover, the experimental Kaiser Effect (KE) points (in Fig. 18 (a)) are qualitatively reproduced in the simulation during the forming microcrack stages (as introduced as stage (II)), which is corresponding to the points a, b, c and d in Fig. 18 (b).
This point demonstrates the high potential of DEM numerical models in simulating the acoustic emission events and their evolution. This point could enormously facilitate the interpretation of experimental acoustic emissions results and even, in the future, help to identify the type of defect associated with a given type of event recorded during mechanical tests.

Conclusions
In this paper, a qualitative and quantitative DEM approach was proposed to model the non-linear quasi-brittle behaviour of a continuum media by introducing a randomisation process that follows a Weibull distribution for the local bond strength. This proposed approach is able to take into account simultaneous microcrack initiations and propagations. Modelling these multi fracturing processes is of a high degree of interest for simulating refractory ceramics because such materials involve a high amount of microcracks which are induced by thermal expansion mismatches between constituents. This fracturing process and 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 microcracking simulations are among the key interests of using DEM in such a case compared to FEM.
In this study, the experimental reference material was Alumina Spinel brick which promotes a typical quasi-brittle behaviour highlighted by a strong non-linear response. To achieve this typical non-linear response, the proposed DEM approach involves a customised calibration process in two steps: 1. a preliminary calibration step capable of reproducing a linear elastic brittle behaviour without considering any randomisation of bond strength and, 2. a complementary calibration step capable to quantitatively reproduce a non-linear quasi-brittle behaviour by introducing bond strength randomisation, which follows an adjusted Weibull distribution. For the linear part, the calibration of the FJM input parameter was done for a perfect linear elastic brittle material. This calibration process led to a perfect brittle linear behaviour, reproducing the initial elastic behaviour of the real material. For the non-linear part, a randomisation process of bond strengths was introduced following a Weibull distribution. Depending on the assigned Weibull modulus (m), this randomisation process influenced the apparent mechanical behaviour resulting from the microcracking process within the simulated material. By decreasing Weibull modulus (m), the microcracking process would increase, and consequently, it induced higher mechanical non-linearities and a lower apparent tensile strength. To achieve quantitative simulations, the secant modulus at the peak was introduced as the parameter to quantify the experimental non-linear behaviour. Thanks to this parameter, finding a relationship between quasi-brittle behaviour and the Weibull modulus (m) was possible. In addition, a relationship between apparent tensile strength and the Weibull modulus (m) was defined. By having these two relationships, it was possible to propose a quantitative calibration process to model a targeted non-linear behaviour. A meta-algorithm was proposed to make this non-linear calibration method more general and applicable for other DEM contact models.
In the last part, a cyclic tensile test was investigated on the reference Alumina Spinel material. In this way, a numerical sample was calibrated by the proposed approach. Results of simulation applied on this numerical sample was compared to experimental observations. The related simulation was designed to reproduce the cyclic loading of the experimental one. The obtained result matched quantitatively to the experimental Alumina Spinel curve. Afterwards, to investigate further the obtained numerical results, the microcrack and elastic properties were monitored. A qualitative comparison with acoustic emissions measurements showed that the simulation could reproduce the experimental Kaiser effect.
Overall, the proposed procedure for calibrating the DEM model to reproduce the non-linear behaviour of the quasi-brittle materials is validated by comparing it to real experimental data. The proposed numerical algorithm could be potentially applied to other DEM bond models and other bond strength distribution laws. Also, this approach can potentially be used to tune the brittleness of material by promoting crack branching processes and diffused damage creation.