Thermal shock resistance of Magnesia spinel refractories – Investigation with the concept of configurational forces

In magnesia-spinel refractories, spinel inclusions are embedded in a magnesia matrix in order to increase the thermal shock resistance by micro-crack toughening. The paper investigates damage evolution in the material during the cooling at the end of the production process by a sophisticated numerical analysis. After a stress-strain analysis on a large representative volume element, the development of the damage zone between spinel grains is simulated with the concrete damaged plasticity model of ABAQUS. Application of the concept of configurational forces allows the physically correct determination of the driving force for cracks in the elastic quasi-plastic matrix material. The results demonstrate that inherent matrix cracks are unable to grow during the cooling, but they strongly influence the initiation and development of the damage zone. The pronounced damage initiation, in combination with the low crack driving force, is the main reason for the good thermal shock resistance of magnesia spinel refractories.


Introduction
According to the theories of Kingery [1] and Hasselman [2], two routes appear possible in order to improve thermal shock behaviour. On the one hand crack initiation could be prevented at all; on the other hand the crack propagation could be hindered. For many refractories, the latter route is preferred because stresses are too high to prevent crack initiation completely. Several figures of merit have been proposed for the characterization of the thermal shock resistance. The Hasselman parameter '''' R based on the ratio of elastic strain to work of fracture reads as follows [2], (1) 2 Here, is the specific fracture surface energy, t  the tensile strength and the Young's modulus.
High specific fracture energy f 2 G   and low tensile strength t  is increasing '''' R . The so-called characteristic length [3], is closely related to '''' R in the form, High values for ch l indicate low brittleness because it means a relatively high specific work of fracture compared to the stored elastic strain energy density. This results in more stable crack propagation and higher resistance against crack propagation, i.e. higher thermal shock resistance.
The main idea behind the addition of spinel (MgAl2O3) to pure magnesia (MgO) refractories is to improve the thermal shock behaviour by increasing the crack growth resistance. Due to the thermal mismatch between spinel inclusions and matrix, damage processes are initiated during the cooling at the end of the production process. The damage initiation is intended because it increases the thermal shock resistance due to the reduction of the tensile strength and the increase in the specific fracture energy caused by the formation of a fracture process zone [4]. This has been observed also in the works of Aksel et al. [5][6][7][8][9], who investigated the influence of different spinel contents with grain sizes <30 µm in very fine grained magnesia. In a similar study, Grasset-Bourdel et al. [10] investigated the influence of coarse spinel grains with sizes between 1 and 3 mm in MgO with a grain size varying between 1 µm to 3 mm. The addition of 5 wt% spinel to pure magnesia nearly doubled the characteristic length and the specific fracture energy. At the same time, strength and Young's modulus were significantly reduced. Higher spinel contents of up to 25 wt% yielded similar results for the specific fracture energy, while the strength was further reduced [10]. This paper investigates the damage evolution in the magnesia matrix between the spinel inclusions during the cooling process by a sophisticated numerical analysis. A Finite Element modelling is performed in two steps. First a stress-strain analysis is conducted on a large representative volume element of the refractory. Subsequently, damage evolution is modelled in detail on a fine-meshed sub-model taken from a critical region of the material. Since the presence of defects in the magnesia matrix cannot be excluded, cases without and with inherent cracks between the spinel grains are considered. Additionally, a new, modern tool of mechanics is applied for this investigation, the concept of configurational forces. A short outline of this concept is given in the following section.

The concept of configurational forces
Configurational forces are thermodynamic forces acting on defects in materials, such as cracks or interfaces [11][12][13]. A configurational force on a considered defect arises, if the total energy of the system varies for different positions of the defect in the material. The defect experiences a driving force to move in such a direction so that the total energy decreases. An outline of the concept of configurational forces is found in the books of Maugin [12], Gurtin [13], or Kienzler and Herrmann [14]. The great versatility is shown, for example, in the references [15][16][17][18][19].
A configurational force appears on such positions in the material where the divergence of the configurational stress tensor is different from zero. The configurational force vector is given by the relation, where the symbol ∇ denotes the Lagrangian gradient operator, the strain energy density, I the identity tensor, the transposed of the deformation gradient F, and the first Piola-Kirchhoff stress. The configurational stress is a second rank tensor, which was introduced by Eshelby as "energy momentum tensor" [11]. The crack driving force is a well-known example for a configurational force. The configurational force vector at the crack tip is given by [20][21] The vector is the unit normal vector to a circle % with radius r centred at the crack tip. It is possible to calculate the scalar near tip J-integral from the configurational force at the crack tip, with e as the unit vector in the direction of crack growth. It is seen from Eq. (6) that the vector points into the negative crack growth direction. It is also well-known that a material inhomogeneity influences both the crack driving force and the crack path, see e.g. [21]. The reason is that configurational forces are induced at the boundary of the inhomogeneity, and they influence magnitude and direction of the configurational force vector at the crack tip, . The configurational force along a sharp interface * is determined by the jumps of the strain energy density and the deformation gradient at the interface [23], In Eq. (7), the symbol ⟦ ⟧ denotes the jump and 〈 〉 the mean value of a quantity, calculated from the limiting values on either side of the interface; n is the unit normal vector to the interface Σ. An example for the application of this concept in order to predict the fracture toughness of a ceramic multilayer has been presented in [23].
The main advantage of the use of the configurational force concept is that it allows the evaluation of the crack driving force independent of the constitutive relations of the material [24,25]. If the material is (non-linear) elastic, the total strain energy density is recoverable and the "non-linear elastic configurational force" nlel f is given by Eq. (4). Since no bulk configurational forces exist in a homogeneous elastic material,  nlel f deduced from configurational forces nlel J is identical to the conventional J-integral derived by Rice [26]. Path independence also applies for elastic-plastic materials that are described by deformation theory of plasticity [24][25][26]. In elastic-plastic materials with incremental theory of plasticity, the total strain energy density consists of an elastic and a plastic part,      e p . Only the elastic part of the strain energy density  e is recoverable. Therefore, the "elastic-plastic configurational force" ep f [24,25] is given by Eq. (4), but with  replaced by  e . The scalar, near-tip incremental plasticity J-integral is given by the relation tip J e f , analogously to Eq. (6). Since bulk configurational forces emerge in the plastically deformed regions of the material, the incremental plasticity J-integral is path dependent. The incremental plasticity J-integral for an arbitrary contour  is evaluated from the equation [24,25], where D denotes the region between the crack tip and the contour  . It has been found in [25] that the crack driving force of a stationary crack under monotonic loading is characterised by the parameter ep PZ J , i.e. the incremental plasticity J-integral for an integration path  PZ , which encloses completely the crack tip plastic zone. Note that the conventional Jintegral approach, or the non-linear elastic J-integral nlel J , do, in general, not provide a measure of the crack driving force for elastic-plastic materials [25,26]. However, the incremental plasticity Jintegral ep PZ J and the conventional J-integral are equal under certain conditions [25].
A big advantage of the configurational forces approach is that it enables the evaluation of the crack driving force even for growing cracks in elastic-plastic materials under monotonic or cyclic loads [25,27]. From the foregoing it becomes clear that the elastic-plastic configurational forces and the incremental plasticity J-integral are applicable also to ceramic materials with quasi-plastic behaviour, as described by the Concrete Damaged Plasticity (CPD) model, see below.

The finite element model
In order to simulate the damage evolution in a magnesia spinel refractory during the cooling process, a model similar to that in [30] was used (Fig. 1). The model consists of circular spinel inclusions in a homogeneous magnesia matrix. The 2-dimensional representative volume element (RVE) was very large in order to correctly reflect the size distribution of the spinel grains. The size distribution of the inclusions was chosen according the function of Dinger-Funk [29] where the percentage of grains smaller than a given grain size d is calculated from the relation, The parameter m denotes the distribution modulus; dmin is the minimum and dmax the maximum grain size. A value of m close to 1 gives a uniform grain size distribution; for m < 1 the smaller grain sizes are dominant. In our model we chose a spinel content of 10 vol.%, and the parameters were  al. [28].
Only tensile cracking in the matrix was permitted; cracking of particles and interfaces was not considered. It was decided to neglect the thermomechanical strains above 1200 °C, since the material relaxes to a nearly stress-free state at these temperatures within short time [31]. The material was only thermally loaded. Because of the low cooling rate, a homogeneous temperature distribution was assumed in the RVE. The calculation started with a homogeneous cooling from 1200 °C, assuming a stress free initial state. Hereby small strain theory and generalized plane strain conditions were assumed. After the finite element stress and strain analysis with ABAQUS, the non-linear elastic and elasticplastic bulk configurational forces, nlel f and ep f , and the configurational forces at the interfaces  f were evaluated from the equations given in the preceding section by a self-written postprocessing routine based on [32]. Subsequently, the non-linear elastic and the incremental plasticity J-integrals nlel J  and ep J  were calculated from the configurational forces, Eq. (8). Note that the evaluation of the conventional J-integral in ABAQUS is not possible for generalized plane strain conditions.

Material behaviour
The spinel grains were assumed as linear elastic. The behaviour of the magnesia matrix was modelled with the Concrete Damaged Plasticity (CDP) model. The CDP model in ABAQUS provides the capability of modelling the inelastic behaviour during the cooling of the refractory. It is based on the works of Lubliner [33], Fenves [34] and Hillerborg's fictitious crack model [35].
A possible normal stress versus crack-opening displacement, σ vs. x, response of the material is shown in Figure 2a. The stress increases linearly until the tensile strength σt is reached. A damage process is initiated at this point, which leads to a decrease of the stress with further opening. For an opening equal or larger than xult, the stress is reduced to zero. The σ vs. x relationship is mainly characterised by σt and the specific fracture energy Gf, which is given by For the CDP model, ABAQUS applies a characteristic element length lel to transform the σ vs. x curve of Fig. 2a into a characteristic stress strain curve of the material, Fig. 2b. Hereby the strain is given by the relation The material behaves linear elastic until the tensile strength σt is reached, followed by a strain softening behaviour at larger strains. For a strain equal or larger than ult, the stress is reduced to zero. The used material properties are listed in Table 1; it was assumed that the values are independent from temperature. Data for the Young's modulus E, the coefficient of thermal expansion , the specific fracture energy G f were taken from former works of the authors and from literature

Results and discussion
Due to the lower thermal expansion coefficient of the spinel grains, tensile circumferential stresses develop in the matrix material during the cooling of the RVE, whereas the grains are under compression. Figure 3a shows the distribution of the maximum principal stresses after a cooling of ΔT = 15K. The maximum circumferential stresses in the matrix appear near the interfaces to the grains. They can reach the magnitude of the matrix tensile strength σt. At some locations, where a narrow matrix corridor lies between two spinel grains, high normal stresses appear even distant from the interfaces. Such regions are most critical for damage initiation. One of the most critical regions in the RVE is marked in Figure 3a and selected for further investigation, see below.
Another critical matrix region is seen near the lower left corner of the RVE. More details about the effect of cooling on the mechanical behaviour of the RVE can be found in Fasching et al. [28]. A sub-model is used to investigate in more detail the damage initiation process in the critical matrix region marked in Fig. 3a. The average mesh size in the sub-model is approximately 12 µm, Fig.  3b. The boundary conditions of the sub-model are found by interpolation of the solution from the global model [30]. A basic question that should be answered in the current investigation is, how the presence of an initial defect influences the process of damage initiation and crack growth in the magnesia matrix, in comparison to the undamaged matrix. Therefore, three versions of the sub-model were generated: (i) defect-free matrix, (ii) matrix with short inherent crack, crack length a0 = 34 µm; (iii) matrix with long inherent crack, a0 = 240 µm. The inherent cracks are located in the direction of a line connecting the centres of the two spinel grains, and they start at the surface of the left grain in Fig. 3b, compare also Fig. 5. Figure 4 shows The configurational force vectors emanating from the crack tips are pointing against the crack growth direction, Eq. (6). It is seen in Fig. 4b that the tip f -vector originating from the right tip and pointing to the lower left is larger than the one originating from the left tip. The reason is that the left spinel grain provides a crack-tip shielding effect, since its Young's modulus is higher than that of the magnesia matrix, see e.g. [40,41]. The tip f -vectors are accompanied by a number of small f-vectors, which appear only due to numerical effects. They decrease to zero at a mesh-dependent distance of less than 50 µm. They should be taken into account for the evaluation of the scalar crack driving force Jtip with Eq. (6), since the bulk configurational force f = 0 in a homogeneous elastic material, e.g. [20,25]. For plane strain or plane stress conditions, the Jtip-value of the right tip could be also evaluated from the conventional J-integral concept by employing the virtual crack extension (VCE) method in ABAQUS [30]. Hereby, allowing for numerical inaccuracies, an integration path of 2 or 3 elements around the crack tip should be chosen. Since the evaluation of the conventional J-integral in ABAQUS is not possible for generalized plane strain conditions, the computation of the J-integral from the configurational force concept is important in our case even for elastic material behaviour.
In elastic-plastic materials, or in cases where the CDP-model is used in order to reflect the damage evolution process in the matrix, the crack driving force must be evaluated from the incremental plasticity J-integral computed for an integration path ΓPZ around the (quasi-) plastic zone, ep PZ J (see Eq. 8). The appearance of quasi-plastic strains in the matrix material simulated with the CDPmodel is causing bulk configurational forces ep f in the damage zone. The damage zone is expanding with progressing cooling. Figure 5 shows the maximum principle quasi-plastic strains after ΔT = 10.5K for the three submodels without inherent crack, with inherent short crack and with inherent long crack. Figure 6 illustrates the same after a cooling of ΔT = 15K. For both values of ΔT, the J-integral values for the right crack tips, J PZ ep and J PZ nlel , are listed in Table 2. It is seen from Fig. 5 that after a cooling of 10.5K, a damage process zone is initiated only in the cases with inherent cracks, but not in the defect-free material. It is seen that the J PZ nlel -values are only little larger than J PZ ep ; the reason is that the damage zone ends far from the interface of the right spinel grain. After a cooling of ΔT = 15K, all sub-models exhibit a pronounced damage zone, crossing the whole matrix corridor between the two spinel grains, Fig. 6. Since the damage zone reaches the interface of the right grain, the J PZ nlel -values are significantly larger than the J PZ ep -values and not reliable anymore. The comparison of the three sub-models shows that the pattern of the damage zone is strongly influenced even by the presence of a very small crack. a) b) c) Figure 5: Maximum principle quasi-plastic strains after a cooling of 10.5K for the models a) without inherent crack, b) with inherent short crack, c) with inherent long crack. a) b) c) Figure 6: Maximum principle quasi-plastic strains after a cooling of 15K for the models a) without inherent crack, b) with inherent short crack, c) with inherent long crack.  It is seen from Table 2 that the values of the crack driving force, ep PZ J , are very small, much smaller than the magnitude of the specific fracture energy, f 100 G  J/m 2 , which would be required for crack propagation. Therefore, the inherent cracks are not able to propagate during the initial stages of cooling. Instead, the matrix material is increasingly damaged by the growth of the damage zone, i.e. when the maximum normal stress exceeds the tensile strength σt of the matrix. J for the matrix material with CDP. The contour iDZ represents the "inner damage zone"; it extends 11 elements directly in front of the inherent crack and one row of elements above and below these 11 elements. It can be deduced from the ep J -curves in Fig. 7 that the increasing damage in the matrix material leads to a decrease of the crack driving force so that the condition for crack growth for an inherent crack is never reached. Instead, the quasi-plastic damage zone between the grains consumes energy, without showing pronounced damage localization. This behaviour might be the main reason for the good thermal shock resistance of magnesia spinel materials. Also the experimental results of Harmuth et. al. [4] and Aksel et. al. [9] support this finding, see Introduction. It is clear that the behaviour might change, if the matrix material had a higher tensile strength σt and a lower specific fracture energy f G .
Here it is important to understand the physical meaning of the incremental-plasticity J-integral calculated around a specific contour Γ: The term ep J  gives the driving force for the combined translational movement of the crack tip plus that part of the quasi-plastic zone, which is enclosed by the contour Γ. The term ep tip J would give the driving force for the crack tip alone, which is, in general, not relevant, since the crack tip cannot move without the simultaneous movement of the process-or quasi-plastic zone around the crack tip [25]. The term ep PZ J includes the whole quasiplastic damage zone. This term is reasonable for characterizing the crack driving force for the cases depicted in Fig. 5b and c. However, this is not so for the cases depicted in Fig. 6. The reason is that in Fig. 6 the damage zone already reaches the interface of the right grain. Therefore, a further translational movement of the whole damage zone to the right is not possible anymore. However, a movement of the crack tip within the "inner" or "near-tip" process zone, as reflected by   Fig. 8 it is noticed that the damage zones form "bridges" between the two grains in this temperature range, and a significant increase in quasi-plastic strain near the interfaces is observed. The distribution of the maximum principal quasi-plastic strain after a cooling of ΔT = 200K is shown for the defect-free matrix in Fig. 9a and for the inherent long crack in Fig. 9b. Figure 10 presents the distribution of configurational forces ep f for the two cases. The figures indicate some fundamental differences. Because of the stress concentration at the tip of the inherent crack, damage initiation is happening earlier, i.e. at a lower temperature difference. A region with very high quasi-plastic strain extends in crack growth direction and ends only a few elements from the interface to the right grain. Main part of the damage zone in Figs. 9b and 10b lies in a triangle formed by the crack tip and the two tangents to the right grain. In contrast, damage appears far more distributed in the defect-free matrix, Fig. 9a. In Figs. 9b and 10b very high quasi-plastic strains and large ep f -vectors are seen in the matrix regions close to the interfaces to both grains, whereas this is less significant in Figs. 9a and 10a. This means that matrix fracture near the interfaces of the grains seems to be strongly promoted by the presence of the crack.

Conclusions
The concrete damaged plasticity model and the concept of configurational forces were applied for the investigation of crack initiation of magnesia spinel products during the cooling from the burning temperature. Evaluations show the strong influence of an inherent matrix crack on the damage evolution: Damage initiation happens at lower T for the cases with an inherent crack, the magnitude and shape of the quasi-plastic damage zone are different. Moreover, a strong influence on the possible formation of near-interface damage was demonstrated. For further investigations, it is recommended to consider also heterogeneities in the matrix.
With the concept of configurational forces, the J-integral can be evaluated either according to the theory of deformation plasticity (resulting in the conventional J-integral J nlel ) or according to the incremental theory of plasticity (leading to the incremental plasticity J-integral J ep ). Only the latter term enables the physically correct evaluation of the crack driving force, not only for the concrete damage plasticity model but for arbitrary non-elastic constitutive relations. The results demonstrated that, for the material data given, inherent matrix cracks are unable to grow during the cooling, since the crack driving force is smaller than the specific work for fracture. Instead a quasi-plastic damage zone appears between the grains, which consumes energy and decreases the driving forces of possible inherent cracks, without showing pronounced damage localization. This behaviour is the main reason for the good thermal shock resistance of magnesia spinel materials.