1 Introduction

Ordinary refractory ceramics possess coarse microstructures comprising aggregates of different sizes, often up to 5 mm, with fine particles and defects. They are typically used in complex thermomechanical environments, such as steel ladles in the iron and steel industry. Tensile fracture is one of the main failure mechanisms caused by severe cold thermal shocks during service [1].

Theoretically, uniaxial tensile testing is the ideal approach to detecting fractures in ordinary refractory ceramics under pure tensile loading conditions. Nevertheless, their inhomogeneous microstructures and larger ratios of defect to aggregate size demand sufficiently large sample dimensions to gain representative results [2,3,4,5]. However, direct tensile tests have strict requirements for the preparation of testing configurations, and special measures should be taken to avoid bending or torsion [6,7,8]. For these reasons, direct tensile tests have been conducted less frequently in the testing of ordinary refractory ceramics than three-point bend tests, Brazilian tests, and wedge splitting tests [9,10,11,12,13]. However, the samples tested with the latter approach often experience different loading conditions, such as a tensile load at the flexural area and a compressive load at the residual area close to the loading point for a three-point bend test [14, 15].

How do ordinary refractory ceramics respond under purely uniaxial tensile loading conditions? What are the fracture mechanisms and patterns in ordinary refractory ceramics with various properties? These questions cannot be easily answered through experimentation. Therefore, the present study used the discrete element method (DEM) to establish a two-dimensional numerical model with a minimum ratio of specimen width to aggregate size of 8 for the numerical uniaxial tensile test.

The stochastic spatial distribution of particles and statistically distributed interface properties at the grain level were considered in the DE modelling, with an emphasis on the inhomogeneity of ordinary refractory ceramics. Furthermore, the crack propagations appearing in the samples during testing were visualised and quantified to provide a failure mechanism analysis at the mesoscale.

2 Methods: discrete element model configuration

The Discrete Element (DE) approach exhibits considerable potential in investigating the behavior of geomaterials at the individual particle level [16,17,18,19,20,21,22]. The approach was first developed by Cundall and Strack [23] and applies rigid circular discs, spherical spheres, and walls as basic components. The movement of each element is determined by Newton's second law, and the contact force is adjusted based on force–displacement relations [24].

For the tensile test simulation, a 2D representation of the refractory samples with a size of 80 mm × 40 mm was established, as shown in Fig. 1. The Dinger-Funk equation (Eq. 1 [25]) with n equal to 0.37 described the grain grading of a refractory ceramic material, where the largest diameter of the particles was 5 mm, and the smallest diameter was 0.088 mm. The cumulative volume percentage curve with respect to particle diameter is shown in Fig. 2 as introduced in the ref [26].

$$ P_{d} = \frac{{d^{0.37} - d_{min}^{0.37} }}{{d_{max}^{0.37} - d_{min}^{0.37} }} $$
(1)
Fig. 1
figure 1

Two-dimensional representation of a sample for tensile test simulation

Full size image
Fig. 2
figure 2

Exemplary cumulative volume percentage curves of a Dinger-Funk equation and from a DE model [26]

Full size image

Pd denotes the cumulative volume distribution of particles with diameters less than or equal to d. The particle size was classified into five ranges: 0.088–1 mm, 1–2 mm, 2–3 mm, 3–4 mm, and 4–5 mm. The particle size was distributed evenly in each group.

The samples were mechanically loaded in the vertical direction to simulate tensile testing. As shown in Fig. 1, two sets of elements were selected from the top and bottom rows of the samples, which are marked in red and green, respectively. A loading strain rate of 0.25 s−1 was chosen as the optimal rate [27] to allow for quasi-static conditions. Consequently, a constant velocity was applied to the two rows of elements, as indicated by the arrows in Fig. 1. The row of elements placed at the top travelled upward at a constant velocity. In contrast, the row of elements located at the bottom moved downward with the same constant velocity. The strain was measured from the particles at the two ends. The stress was measured by creating a single measurement circle with a radius of 20 mm at the centre of the sample. Stress was computed by averaging the contact forces within the circle.

The interface properties between the elements were defined using the linear parallel bond model in the PFC code [28]. The linear parallel bond model, seen in Fig. 3, describes the interaction between two components by considering two separate interfaces: a bonded interface and an unbonded interface. The unbonded interface is incapable of withstanding tension and relative rotation. To prevent the occurrence of sliding, a restriction on the shear force is imposed, known as the Coulomb limit. The bonded interface has the capacity to endure displacement and relative rotation until it hits its maximum strength threshold. When bond is broken, only the unbonded interface plays role. gs, μ, kn, and ks denote the surface gap, friction coefficient, normal stiffness, and shear stiffness of the unbonded contact, respectively. The parameters that represent the tensile strength, cohesion, friction angle, normal stiffness, and shear stiffness of the bonded interface are symbolised as \(\overline{\sigma }_{c}\),\(\overline{c}\),\(\overline{\phi }\), \(\overline{k}_{n}\)and \(\overline{k}_{s}\), correspondingly.

Fig.3
figure 3

Illustration of linear parallel bond contact in PFC [29] a parallel bond contact between two circular elements, b bonded interface law, c unbonded interface law

Full size image

The determination of the effective moduli and stiffness ratio, as shown in, Eqs 2-6 provides an alternate approach to characterising the features of the interface.

$$ \left\{ {\begin{array}{*{20}ll} { k_{n} = AE^{*} /L , unbonded } \\ {\overline{k}_{n} = \overline{E}^{*} /L , bonded } \\ \end{array} } \right. $$
(2)
$$ \left\{ {\begin{array}{*{20}ll} { k_{s} = k_{n} /k^{*} , unbonded} \\ { \overline{k}_{s} = \overline{k}_{n} /\overline{k}^{*} , bonded} \\ \end{array} } \right. $$
(3)
$$ A = 2rt, \left( {t = 1} \right) $$
(4)
$$ r = \left\{ {\begin{array}{*{20}l} {\min \left( {R^{\left( 1 \right)} ,R^{\left( 2 \right)} } \right), ball - ball } \hfill \\ {R^{\left( 1 \right)} , ball - facet} \hfill \\ \end{array} } \right. $$
(5)
$$ L = \left\{ {\begin{array}{*{20}ll} { R^{\left( 1 \right)} + R^{\left( 2 \right)} , ball - ball} \\ {R^{\left( 1 \right)} , ball - facet} \\ \end{array} } \right. $$
(6)

The effective moduli of the interfaces that are not bonded and that are bonded are represented as E* and \(\overline{E}^{*}\), respectively. The ratios of normal to shear stiffness for the unbonded and bonded interfaces are denoted as k* and \(\overline{k}^{*}\)*, respectively. The symbol A is used to represent the contact area, while r is used to designate the radius of the elements. L is used to signify the distance between the contacted elements, while R(1) and R(2) are used to represent the radii of the contacted elements.

Calibration procedures were performed as described previously [26]. A series of experimental tests of a brittle refractory ceramic material (MgO-based) were conducted, namely cold crushing and wedge splitting tests, and then compared the actual results with simulation curves generated from a 2D DE model of identical testing with homogeneous interface properties defined by a linear parallel bond model in the PFC code. The particle contact parameters were inversely identified by the use of the NL2SOL [30]. The inversely identified parameters are listed in Table 1 as reference.

Table 1 Reference parameter set of the linear parallel bond model [26]
Full size table

The present study focuses on investigating the effects of the particles’ spatial and interface characteristic differences on the mechanical behaviour of tensile tests. For the particle spatial distribution, three models were established with the same minimum element size, grain grading, and homogenous interface features. A random seed was used to assign different particle spatial distributions, and its default value in PFC 5.0 was 10,000. With fixed seeds, the element locations and radii were identical for each assembled model. However, the seed used influences the initial spatial configuration of the elements within a limited space. As shown in Fig. 4, the structures of the samples with random seeds 10,057, 16,857, and 20,057 differed.

Fig. 4
figure 4

Particle spatial distributions resulting from different random seeds a seed number 10057, b seed number 16857, c seed number 20057

Full size image

A case with a constant interface served as a reference for investigating the effects of different contact conditions. Uniform and Weibull distributions with shape parameters 3, 4, and 5 were considered for the interface properties. The median value of each parameter in the cases with uniform and Weibull distributions was the same as the corresponding parameter in the cases with homogeneous interface properties, as shown in Table 1. The parameters were modified identically for the uniform and Weibull distributions; that is, the contact with the highest strength had the highest friction angle. Figure 5 illustrates the assignment of the interface properties using tensile strength as an example. The constant-distributed and uniformly distributed parallel-bond tensile strengths are shown in Fig. 5a and b, respectively. The parallel bond tensile strength assigned to shape parameter 3 of the Weibull distribution is shown in Fig. 5c, and those with shape parameters 4 and 5 are presented in Figs. 5d, e, respectively. Figure 5f depicts the constant case and distribution curves of the tensile strength for the uniformly distributed case and the three Weibull-distributed cases.

Fig. 5
figure 5

Assigned tensile strength (Pa) in the cases with a constant [27], b uniform distribution [27], and c Weibull distribution with a shape parameter of 3, d Weibull distribution with a shape parameter of 4 [27], e Weibull distribution with a shape parameter of 5, and f tensile strength cumulative curves for all the cases

Full size image

3 Results and discussion

3.1 Influence of particle spatial distribution

The number of particles, contacts, and mechanical properties in models assembled with different seeds are listed in Table 2, and the value in brackets shows the relative difference of the respective quantity from that of the reference case with a seed of 10,057.

The stress-strain curves for three cases are shown in Fig. 6. The tensile strength, fracture energy, secant Young's modulus, and brittleness number were calculated using Eqs. 7–10.

$$ \sigma_{t} = \frac{{F_{max} }}{{A_{0} }} $$
(7)

where σt denotes the tensile strength; Fmax is the maximum force; and A0 is the cross-sectional surface area.

$$ G_{f} = \frac{1}{{A_{0} }}\mathop \int \limits_{0}^{{\delta_{ult} }} Fd\delta $$
(8)

where Gf denotes the specific fracture energy; δult is the ultimate displacement; and F and δ are the force and displacement, respectively. For the curves in Fig. 6, the maximum displacement used for the calculation was the point at which the residual force reached a maximum of 15%.

$$ E_{static} = \frac{{\sigma_{fc} }}{{\varepsilon_{fc} }} $$
(9)

where Estatic is the static Young's modulus; σfc is the stress immediately before the first crack occurs; and εfc is the strain immediately before the first crack.

$$ B = \frac{{\sigma_{t}^{2} \cdot L}}{{G_{f} \cdot E}} $$
(10)

where B and L denote the brittleness number and length of the specimen, respectively [31].

Table 2 Number of particles, contacts, and mechanical properties in cases with different seeds
Full size table
Fig. 6
figure 6

Stress–strain curves of cases with different particle spatial distributions

Full size image

The crack density and crack density increment are calculated using Eqs. 11 and 12.

$$ c\rho_{i} = \frac{{n_{i} }}{N} $$
(11)

where cρi and ni are the crack density and number of cracks at the ith loading point, respectively, and N is the number of contacts.

$$ c\rho_{in} = c\rho_{i} - c\rho_{i - 1} $$
(12)

where cρin is the increment in crack density.

The maximum relative difference was observed in the brittleness number of the model with a seed of 20057, which was 4.7%. In general, the seed number has a slight influence on the particle number, contact number, and mechanical properties.

Figure 7 shows the evolution of crack density during the loading process. The total crack density varies with different particle spatial distributions, particularly after reaching 80% of the peak stress in the ascending part, as shown in Fig. 7a. The ratio of tension to shear crack density is shown in Fig. 7b. In all three cases, shear-induced cracking dominated at 80% of the peak stress before the peak load, and tension-induced cracks dominated after the peak stress. The ratio of tension to shear crack density increased from 80% of the peak stress in the ascent to the peak stress for cases with the random seeds 16,857 and 20,057. For stresses higher than the peak value, tensile cracking was dominant. The ratio decreased with decreasing residual stress for the case with a seed of 16,857, whereas the ratio in the case with a seed of 20,057 decreased from the peak stress to the residual stress, reaching 80% of the peak stress, and then increased. In the model with a seed of 10,057, the ratio increased from 80% of the peak stress in the ascending portion to the residual stress, reaching 80% of the peak stress. Tensile cracking dominated after the residual stress reached 80% of the peak stress.

Fig. 7
figure 7

Crack propagation in cases with various particle spatial distributions a Total crack density b Ratio of crack density of tension to shear

Full size image

Figure 8 shows the fracture patterns of the samples with various particle spatial distributions. The main fracture paths were close to the top surface of the samples with seeds 10,057 and 20,057. In contrast, for the sample with a seed number of 16,857, the main fracture path was farther from the top surface. Additionally, two short fracture paths were observed in the samples with seeds 16,857 and 20,057.

Fig. 8
figure 8

Crack patterns of cases with different particle spatial distributions

Full size image

3.2 Influence of statistically distributed interface property

3.2.1 Mechanical properties

Figure 9 shows the stress–strain curves for cases with different particle spatial distributions realised with different seed numbers (10,057, 16,857, and 20,057) and contact characteristics (constant, uniform, and Weibull distributions). Three shape parameters (3, 4, and 5) were considered for the Weibull distribution.

Fig. 9
figure 9

Stress–strain curves of cases with different interface properties

Full size image

In all the cases shown in Fig. 9, at the beginning of the loading procedure, the ascending curves for the cases with different interface properties overlapped. They diverged with increasing strain. The case with a uniform distribution of interface properties displayed the lowest tensile strength and displacement. The cases with constant interface properties showed the highest tensile strength and displacement when the seeds were 10,057 and 16,857, respectively. When the seed was 20,057, the case with shape parameter three exhibited the highest displacement. Generally, the tensile strength is enhanced with an increase in the shape parameter.

Table 3 lists the mean tensile strength, static Young's modulus, specific fracture energy, and brittleness number of the simulated cases. For cases with uniformly distributed interface features, the stress was too small when the first crack occurred. Therefore, the first crack at the beginning of the simulation resulted in an unreasonable value of the static Young's modulus. Due to the fact that the average stress value right before the occurrence of the first crack was around 10% of peak stress in situations with Weibull and constantly distributed interface properties, for cases with uniformly distributed interface properties, Young's moduli were calculated using the stress and strain values from the origin point to 10% of the peak stress in the ascending portion of the stress–strain curves. The case with constant interface properties showed the highest tensile strength, a static Young's modulus, specific fracture energy, and brittleness number. The case with uniformly distributed interface properties possessed the lowest tensile strength, static Young's modulus, specific fracture energy, and brittleness number. The mean tensile strength in the case with uniform distribution was 42.9% of that of the case with constant interface properties; the specific fracture energy of the case with uniform distribution was 31.8%. The mean tensile strength of the cases with a Weibull distribution was 74.5–92.1% of that of the case with constant interface properties, and the mean specific fracture energy was 70.3%–88.2% of that of the case with constant interface properties.

Table 3 Mean tensile strength, static Young's modulus, fracture energy, and brittleness number of simulated cases
Full size table

The Weibull shape parameter has an insignificant effect on the static Young's modulus. Taking the case of Weibull shape parameter 3 as a reference, the other two cases showed relative differences in the tensile strength and specific fracture energy of up to 25.4%. As the Weibull shape parameter increases, the brittleness also increases.

Table 4 displays the mean number of contacts with a given tensile strength interval in cases with a Weibull distribution for interface properties. The mean number of contacts with a strength of no more than 18.12 MPa was 236 for the case with Weibull shape parameter 3; 166 for the case with Weibull shape parameter 4; and 70 for the case with Weibull shape parameter 5. The number of contacts in the cases with Weibull shape parameters 3, 4, and 5 were 4569, 2333, and 1669, respectively, in which the tensile strength was greater than 54.36 MPa. Despite more contacts with high tensile strength in the case with Weibull shape parameter 3, this case had the lowest tensile strength. As expected, the results indicate that weak contacts significantly decrease tensile strength.

Table 4 Mean number of contacts with specified tensile strength intervals
Full size table

3.2.2 Crack density

Figure 10 shows the evolution of the crack density during loading. In Fig. 10a, the mean total crack density of the case with uniformly distributed interface properties is much larger than that of the other instances. When the relative stress reached 80% of the peak stress in the rising portion of the stress–strain curve, the mean total crack density in the cases with a Weibull distribution of the interface properties was larger than that in the case with constant interface properties. The case with Weibull shape parameter 3 showed the largest total density among the three cases with a Weibull distribution. The mean shear crack density exhibits a tendency similar to that of the mean total crack density. The mean tensile crack densities differed in the simulation cases. The difference in the mean tensile crack densities of the cases could not be easily identified.

Fig. 10
figure 10

Crack density with respect to stress relative to the peak stress a mean total crack density, b mean tensile crack density, c mean shear crack density, d mean crack density of tension to shear, e mean increment of tensile crack density, f mean increment of shear crack density

Full size image

In cases with uniformly distributed and Weibull shape parameter 3 distributed interfaces properties, Figs. 10b, c reveal that shear failure dominated the cracking process. Cases with constant and Weibull shape parameters 4 and 5 distributed contact features exhibited prevalent tensile cracks.

Figure 10d shows the mean ratio of the crack density of tension to shear starting from the peak stress because, in some cases, shear cracking did not occur before the peak stress. For the situation with constant and Weibull shape parameter 5 distributed interface properties, tensile cracks dominated the failure mode, although shear cracks increased after the peak stress in the descending portion of the stress–strain curve as the load decreased. Similar to the case with a constant interface property, tensile cracks in the case with Weibull shape parameter 4 dominated the distributed interface property, whereas in the case with Weibull shape parameter 4 distributed interface property, tensile cracks increased more than shear cracks after the peak stress when the load decreased.

The increase in tensile crack density during loading is shown in Fig. 10e. Tensile cracks began to grow at 80% of the peak stress in the ascending portion of the stress–strain curve in the case with constant interface properties and at 60% of the peak stress for the residual cases. Beyond 90% of the peak stress in the descending portion of the stress–strain curve, the increment decreased.

Figure 10f illustrates the increase in shear crack density during loading. For the case with uniformly distributed interface properties, the shear crack density started to increase prior to 40% of the peak stress during the ascending portion of the stress–strain curve. The increment decreased at 60% peak stress but continued to increase until 90% of peak stress during ascending portion of the stress–strain curve. After the residual stress reached 90% of peak stress, the mean increment in the shear density decreased. For the cases with Weibull-distributed interface properties, the shear crack density increased considerably after 80% peak stress in the ascending portion of the stress–strain curve and dropped evidently after 90% peak stress in the descending portion. For the case with constant interface properties, the shear crack density increased significantly at 90% peak stress after the peak load and then decreased.

3.2.3 Crack patterns

Figure 11 illustrates the ultimate fracture patterns of samples with different interface properties at a residual stress of 15% of the peak stress, with the particle spatial distribution controlled by the seed 16,857. The simulation revealed distinct fracture patterns in the five cases. The case with constant interface properties exhibited two noticeable localised fracture bands: the longer band was located at approximately one-third the length of the specimen below its top surface and exceeded the centre of the sample, whereas the shorter band was near the middle of the sample and had a width that was only half of the specimen width. In the case of the uniformly distributed interface, cracks were dispersed throughout the sample.

Fig. 11
figure 11

Crack patterns of cases with various interface properties a constant and b uniform distributions, c Weibull distribution with a shape parameter of 3, d Weibull distribution with a shape parameter of 4, and e Weibull distribution with a shape parameter of 5

Full size image

When the Weibull shape parameter was 3, some cracks occurred randomly in the sample with three distinct fracture bands, and the longest fracture band was close to the top end. The model with Weibull shape parameter 4 had four localised fracture bands, three of which appeared close to the top end and one in the middle of the sample. The longest crack path was approximately one-third that of the specimen’s length beneath its top end. The example with a Weibull shape parameter of 5 had two distinct localised fracture bands; their locations were similar to those of the fracture bands of the case with constant interface properties, one of which was located at approximately 1/3 of the distance beneath the specimen's top end, and the other was close to the middle of the specimen. The lengths of the two fracture bands were greater.

3.3 Comparison between cold crushing test and tensile test simulation results

The simulation results of a cold crushing test and tensile test are compared in this section. The sample was generated with a grain grade obtained using the Dinger-Funk equation with an n of 0.37 and a minimum element size of 0.088 mm. The particle spatial distribution was determined for seed number 10057, and interface properties were assigned using a Weibull distribution considering a shape parameter of 4.

The mechanical properties of ordinary refractory ceramics with the same microstructure exhibited evident asymmetry under compressive and tensile loads, as listed in Table 5. The cold crushing strength was 47.73 MPa, which is nearly four times the tensile strength of 12.56 MPa. The static Young's modulus was measured to be 139.02 GPa from the compression test, but it was 78.62 GPa with tensile loading. The strain values used to calculate Young's modulus for the tensile test were measured from the particles at the two ends, whereas those for the cold crushing test were measured from gauge particles. In the crushing test, the total crack density was approximately ten times that in the tensile test. In addition, for the given interface properties and the given particle spatial distribution, tensile-induced cracks were predominant during the crushing test, whereas shear-induced cracks were predominant during the tensile test. The specific fracture energy calculated from the load–displacement curve of the crushing test was 11.32 times that of the tensile test.

Table 5 Comparison between mechanical properties and crack propagations measured during crushing and tensile tests
Full size table

Figure 12 shows the crack patterns under crushing and tensile loading conditions. The main fracture bands in the sample after the crushing test were inclined at 45º to the horizontal direction, whereas the tensile test cracks were nearly perpendicular to the loading direction. Moreover, the cracks were dispersed in the specimen after the crushing test, whereas only two major crack bands appeared after the tensile testing. The crushing test showed a substantially higher crack density than the tensile test.

Fig. 12
figure 12

Crack patterns of cases with different loading conditions

Full size image

Ordinary refractory ceramics have been characterised in the lab to exhibit anisotropic fracture and creep behaviour under different loading conditions [32,33,34,35,36,37,38,39]. The dynamic Young's modulus is often applied in the characterisation and modelling of refractory behaviours. In contrast, few publications have focused on the determination of the static Young's moduli of ordinary refractory ceramics under tension and compression using mechanical devices and demonstrating their asymmetry under tension and compression. For instance, Simonin et al. [32] reported that the relative difference in the modulus of a high alumina castable between tension and compression determined through a bend test lies between 8% and 30% in the temperature range of 200–1600 °C. Nazaret et al. [33] showed that the static Young's modulus of a geopolymer-based cordierite castable under tension is 1.9 times its static Young's modulus under compression. Schmitt et al. [8] found that the static Young's modulus under tension of a magnesia material with a resin binder that shows a typical nonlinear behaviour is 2.5 times the static Young's modulus under compression, whereas the static Young's modulus under tension of a magnesia material with a pitch binder that shows a typical brittle behaviour is 68.4% of the static Young's modulus under compression. The present numerical simulations also revealed that the ratio of Young's moduli of this MgO-based material under tension to that under compression was 56.6%.

4 Conclusion

A direct tensile testing approach was simulated using a two-dimensional discrete element model to gain an advanced understanding of fracture in ordinary refractory ceramics under uniaxial tensile loading conditions. The stochastic spatial distribution of the particles was controlled using different seed numbers. Statistically distributed interface properties between particles were also considered, and the case with constant interface properties was used as a reference.

The stochastic spatial configuration of the particles had an insignificant influence on the tensile strength, specific fracture energy, static Young's modulus, and brittleness number for the given model setup. The spatial configuration of the particles evidently affects crack propagation, which yields different fracture patterns. A few cracks were dispersed in the samples, and some cracks were localised. No single fracture path was observed in the direct tensile tests. The main non-straight crack path formed was close to one end of the sample and passed through the sample perpendicular to the loading direction.

The bond conditions of the particles in ordinary refractory ceramics cause significant differences in the load–displacement curves, mechanical properties, crack propagations, and patterns under uniaxial tensile loading conditions. The case with uniformly distributed interface properties shows evenly dispersed cracks, whereas the cases with Weibull-distributed interface properties display typical localised fracture paths observed in laboratory mechanical tests. This also indicates that the tensile strength is controlled mainly by weak contacts between the particles, and the brittleness of ordinary refractory ceramics can be reduced when the contact between the components is imperfect.

A comparative study of the mechanical properties obtained from a cold crushing test and tensile test simulation shows asymmetric phenomena in the mechanical responses under uniaxial tension and compression. The crushing strength of a MgO-based brittle material is 3.8 times its tensile strength; the static Young's modulus under tension of a MgO-based brittle material is 56.6% that under compression; and the specific fracture energy under compression is 11.3 times that under tension. With the given Weibull-distributed interface properties defined by shape parameter 4 and the given particle spatial distribution, tensile cracking dominated in the cold crushing test, whereas shear cracking dominated in the direct tensile test.