Analysis of 3D crack patterns in a free plate caused by thermal shock using FEM-bifurcation

Damage to components made of brittle material due to thermal shock represents a high safety risk. Predicting the degree of damage is therefore very important to avoid catastrophic failure. An energy-based linear elastic fracture mechanics bifurcation analysis using a three-dimensional finite element model is presented here, which allows the determination of crack length and crack spacing for a defined thermal load in a free plate. It is assumed that a hierarchical crack pattern is formed due to cooling penetration. The constant growth of the ideal regular pattern of hexagons can change into a pattern with a different symmetry by slightly changing the cooling conditions. This bifurcation point is determined by the second derivative of the mechanical potential with respect to the geometry of the crack front. The very high computational effort for the second derivative is reduced by describing the three-dimensional crack front with a limited number of Fourier coefficients. A one-dimensional transient temperature field at a sufficient distance from the plate edge is assumed. For alumina, the crack length and crack spacing curves are computed for different quenching temperatures and heat transfer coefficients. The corresponding final crack lengths are also calculated as a measure of damage. Comparison with a two-dimensional model confirms the expected 1/2 difference in crack spacing. Data from thermal shock experiments are also presented. However, due to the cracks caused by the strong cooling at the edge, these correspond to the results of the two-dimensional model.


Introduction
The study of thermal shocks in brittle materials and the accompanying damage in form of cracks employed by scientists for decades. Fascinating complex crack pattern formations arise during thermal shock. Beside this it is also important to know the residual strength after a thermal shock case to avoid catastrophic failure in technical systems such as gas turbines. Engineering ceramics are widely used as thermal barriers, as they can withstand very high temperatures. This application is always critical due to the relatively low fracture toughness of ceramics. It was analyzed in Kingery (1955) and Hasselman (1963Hasselman ( , 1969 which material properties affect thermal shock damage. This was investigated experimentally by quenching alumina rods in water (Hasselman 1970). An overview of various failure phenomena and related thermal stress resistance parameters is given in Hasselman (1985). The first investigations of 2D hierarchical parallel crack patterns and the associated analysis with a bifurcation instability criterion were carried out in Nemat-Nasser et al. (1978) and Bažant et al. (1979). Based on this criterion, the 2D crack patterns that occur when glass quartz plates are quenched in water (Bahr et al. 1986) were qualitatively investigated in Bahr et al. (1988) using the Boundary Element Method (BEM). It was found that after thermal shock the final crack length has a significant influence on the strength degradation. Nemat-Nasser et al. (1980) then published another bifurcation instability criterion which was also used in later publications (Bahr et al. 1996(Bahr et al. , 2010Hofmann et al. 2011). In Bahr et al. (1992) another bifurcation criterion based on the mechanical potential and an eigenvalue analysis was shown. This criterion is equivalent to that of Nemat-Nasser et al. (1980). Furthermore, in Bahr et al. (1992) the correctness of the assumption of the hierarchical crack pattern, that one crack stops and another one grows, and the influence of a disturbance on it was shown. Many publications investigated crack patterns by experiments and by various simulation methods, e.g. Shao et al. (2011, and Xu et al. (2016a). However, they all focus on 2D crack patterns.
In Bourdin et al. (2014), a scalar damage variable within a non-local damage model was used to simulate the growth of 3D columnar joints in thermal shocked ceramics. It was shown that this can explain the formation of imperfect polygonal patterns and their selective coarsening during propagation, but an analysis of the 3D crack pattern geometry is still missing. 3D crack patterns of basalt columns were investigated with the Finite Element Method (FEM) in Bahr et al. (2009) by a fracture mechanic bifurcation analysis based on the local energy release rate. This research is similar to the investigated 2D cases in Bahr et al. (1988Bahr et al. ( , 2009Bahr et al. ( , 2010. In order to determine the 3D crack contour, an elaborate gradient method was used for the bifurcation analysis. To avoid this cumbersome iteration of the crack front, a new effective calculation method was developed based on the global energy release rate G and on a Fourier series expansion of the 3D crack front in order to simplify the bifurcation analysis. The application of this method has been demonstrated in Anderssohn et al. (2018). To enable a comparison with Bahr et al. (2009), the growth of 3D crack patterns of basalt columns was also investigated with a steady-state temperature field. The method of Anderssohn et al. (2018) can now be used to investigate the growth of other crack patterns, e.g., generated by drying or, as here, by a transient thermal shock load case for a finite length.
In the present work, we use the ideal hexagonal three column pattern from Bahr et al. (2009) and the Fourier series expansion of the 3D crack front contour from Anderssohn et al. (2018). An FE model for an infinite free plate of thickness 2b is built, which allows for calculating the resulting normalized crack spacing L/b and crack length a 0 /b due to thermal shock. Finally, the goal of the present study is to predict the normalized final crack length a end /b which is crucial for the residual strength of ceramic components that were damaged by thermal shocks. This paper is organised as follows: Sect. 2 presents the analytical temperature field and the basic thermomechanical equations. In addition, the underlying model of the periodic contiuable ideal hexagonal three column pattern is described. The presentation of the theory of eigenanalysis of the global mechanical potential together with the development of the crack front geometry as a Fourier series, which was used to find the bifurcation points, concludes the section. In Sect. 3, the numerical parameters of the FE model and the central Finite Difference Method (FDM) used to determine the derivatives of the mechanical potential are explained. A convergence analysis and the effects of different numbers of Fourier series coefficients on the result are shown. Successful calculations for different temperature differences ΔT and heat transfer coeffi-cients h were performed in Sect. 4. The influence of ΔT and h on the crack length a 0 /b and the crack spacing L/b is discussed. For the validation of the 3D model, thermoshock experiments and calculations for an already existing 2D model were performed, evaluated and compared with the data of the 3D model in Sect. 5. Furthermore, a novel evaluation method based on Computed Tomography (CT) was applied to the thermo-shocked specimens in this work. Note: Similar experiments have already been performed, for example, in Shao et al. (2010) and Xu et al. (2016b). Unfortunately, no material data are given in these publications and thus cannot be used for validation. Finally, in Sect. 6, several conclusions are drawn and an outlook on possible future applications of the model presented in this paper is given.

Temperature profile
First, a plate is heated to the temperature T 0 and then quenched in cooling liquid of temperature T 1 . The area under consideration is the center of the plate far away from the edges. In the center the temperature field can be described by a simple 1D analytical solution, while at the edges a more complex temperature field is expected. It is assumed that the cracking process has no influence on the temperature field. This assumption is confirmed by two points: firstly, the cracks are orthogonal to the isotherms that means that they do not interfere with the heat flux. Secondly the crack opening is very small so heat transport through the air filled cracks is very small compared to the heat transport through the solid. Due to the short duration of the thermal shock process, this is a common hypothesis (Bahr et al. 1988;Li et al. 2013). Therefore, the transient temperature profile is used for an infinite plate with a thickness of 2b and with symmetrical boundary conditions (Tautz 1971) where δ = √ 4Dt, μ n = hb/κ cot(μ n ) and ΔT = T 0 − T 1 . Here, μ n , n = 1...∞ are the positive eigenvalues, κ the thermal conductivity, h heat transfer coefficient and D the thermal diffusivity. In Tautz (1971), the origin of the z coordinate is in the centre of the plate, whereas in the present study it is at the surface. This shift of the coordinate system is implemented in Eq. (1) by the −μ n term in the second cosine function (see also Bahr et al. 1987Bahr et al. , 1988. For short times δ 2 /4b 2 1, (it should be noted that in the literature often τ for δ 2 /4b 2 is used) the number of the used eigenvalues μ n must be high enough. Otherwise Eq. (1) is not valid (Bahr et al. 1987;Martin et al. 2019). 50 and 100 eigenvalues were calculated with the mathematic program MATLAB. For δ/b = 0.01, the maximal difference of the summation term in Eq. (1) was approx 3.6 · 10 −4 . In the simulations the smallest time is δ/b = 0.1, therefore 100 eigenvalues are more than enough.
In Fig. 1 the temperature function Eq. (1) is plotted for different hb/κ and δ/b. The greater hb/κ is, the better is the heat transport between the heated plate and the cooling liquid. An increase of hb/κ leads to an increase of the temperature gradient and therefore the thermal stresses also increase. Because of that, the value of the heat transfer coefficient h is very important for the determination of the damage by the thermal shock. It should be noted that the experimental determination of h is very difficult. In Singh et al. (1981), Zhou et al. (2012), and Jiang et al. (2012) different methods were used to determine h in the case of a thermal shock of alumina in water. The results for the heat transfer coefficient h range from 10 3 Wm −2 K −1 to 10 5 Wm −2 K −1 .

Fundamental thermoelastic equations
Due to the thermal shock, the ceramic plate experiences shrinkage leading to mechanical stresses. In Takeuti and Furukawa (1981) it is shown that the ratio of V = v p b/D decides if acceleration terms have to be considered for the penetration of a temperature field.
. The thermal diffusivity can be calculated with D = κ/ρC p . However, for the materials used in this work, neither D nor the specific heat capacity C p is given, see Table 1. From Zhou et al. (2012), D = 5.4 · 10 −6 m 2 s −1 for alumina with similar purity (99.5%) and density (ρ = 3.98 · 10 3 kg m −3 ) can be taken. For the Al 2 O 3 material with 99.7% purity used in the experiments, with a plate thickness of b = 3.5 · 10 −3 m, V = v p b/D = 6.9 · 10 6 is obtained. For such large values of V , according to Takeuti and Furukawa (1981), p. 118, the ratio of dynamic and quasi-static stresses is equal to one. Thus the acceleration terms can be neglected. Volume forces are not present and the local mechanical equilibrium conditions read where σ i j are the mechanical stresses. Eq.
(2) must be fulfilled in every point of the structure. We assume an isotopic linear thermo-elastic material with linarised kinematics for which the strain ε i j can separated into an elastic and a thermal part: Here, u i denotes the displacement, E is the Young's modulus, ν the Poisson's ratio and δ i j the Kronecker symbol in the elastic strain ε el i j . The third term describes the thermal strain ε th i j . Here, α is the thermal expansion coefficient and T (z, δ) − T 0 the temperature difference according to Eq. (1). T 0 is the reference temperature at which no thermal stresses are present.
For a plate having a 3D crack pattern, to which an instationary temperature field like Eq. (1) is applied, unfortunately no analytical solution of Eq. (2) is available. Therefore, the solution must be obtained numerically. In the present research, this is done by the FEM through the program Ansys. Please note that an analytical solution of the stress field which fulfills Eq. (2) for a plate without cracks is given in Timoshenko and Goodier (2017) and Parkus (1968).

Model development and boundary conditions
In experimental investigations on dried starch-water suspensions, it was found that a regular hexagonal cracking pattern causes three columns to merge into a larger one, see Fig. 2a. A similar type of column could be observed in basalt columns, which were formed by thermal contraction of solidified lava. In both cases, the   tensile stress caused by drying or cooling is the reason for the appearance of the cracks (see also Goehring et al. 2006;Goehring 2008). The process of merging three columns into one is periodically repeatable in the depth (z-direction) and also in the width of the material. This was used to form the model described below. Fig. 2b shows the idealised model with three regions. In region (I) there are three hexagonal columns of equal size with column diameter L. These columns are caused by the steady propagation of the crack front and are referred to as fundamental solution. The model then assumes that the crack stops growing at the (bp) point in the middle of the three columns. So there is another solution besides the fundamental solution. The point (bp) is hence a bifurcation point in the mathematical sense.
In region (II), the crack front between s/p = 0 and s/p = 1 in Fig. 3a stops growing. With this termination, the high symmetry of the hexagons is lost and mixed mode crack propagation with curved crack faces occurs. The three columns merge in region (III) to form a larger column with a new column diameter of √ 3L. It should be noted that region (II) after the bifurcation point is called the post-critical solution.
As in Anderssohn et al. (2018) and Bahr et al. (2009), we use symmetry and periodicity to reduce the computational costs. The representative volume element, used for the calculation in the case of a bifurcation, is one half of a hexagonal column (gray area in Fig. 3b). This is sufficient due to the periodic repetition. In Fig. 3b this is shown by the three points 0 in which the crack stops in the other two surrounding three hexagonal column configurations. This leads to a point symmetry with respect to point 2 with the following coupling of the displacements in the ligament surfaces E and F of the representative volume element in Fig 3b: Here, the index n stands for the normal direction of the ligament surfaces and s is the coordinate of the crack front. Furthermore, the stresses must be consistent at the surfaces E and F. The plate is not subjected to external kinematic constraints, so that the two surfaces G and the ligament surface at D are unloaded in the normal direction, i.e. 0 = σ i j n j d A and u i n i = constant .
This condition was implemented by the Ansys command CE with the requirement that all FE nodes have the same displacement in the surface normal direction. The surface of the plate z = 0 and the crack surface are unloaded, i.e. σ i j n j = 0. The lower face, at z = b, is the symmetry area at which the boundary condition u i n i = 0 and σ i j n j = 0 (i = j) applies.
In the case of uniform column growth, a twelfth of the column is sufficient for the calculation since a higher symmetry is present than in the case of column merging. The ligament face B, shown in Fig. 4, is forcefree in the normal direction according to Eq. (5). The faces A, C, and lower face, at z = b, are symmetry areas for which the boundary conditions u j n j = 0 and σ i j n j = 0 (i = j) hold. The top surface and the crack surface are free of forces σ i j n j = 0 as in the half column model.
for the parameterization of the crack front. Here, a 0 is the crack length, j the number of cosine elements, C i the time-invariant coefficients,p = L/ √ 3 the column side width shown in Fig. 3b, and s the crack front curve variable. To ensure that the Fourier series can be used as the crack front geometry for the fundamental and also for the bifurcation solution, the product 2p is included. The function a(s) is the crack length measured from the surface of the plate, as shown in Fig. 4.
For steady growth, the crack front geometries of the three-column configuration must match the opposite column. Due to the symmetry and periodicity of this configuration, i = 4, 8, . . . are the only coefficients in Eq. (6) that are nonzero. Due to the huge increase in computation time by increasing the number of coefficients and the low effect on the result (see Fig. 9), the highest Fourier coefficient C 4 is selected.
We introduce here for the case of quasi-static crack growth, the change in the global mechanical potential (see for example Kuna (2013, pp. 42-44)) where the elastic strain energy reads and the fracture energy W a is the external work. In Eq. (7), the dot above the symbols represents the change between two states. These states can be two different times or, as in this study, the first state is before the crack extension and the second state is after the crack extension ∂ A. As there are no external forces or supports acting on the free plate, Furthermore, the kinetic energy equals zero. The thermal energy has no influence on the crack growth and thus the change in the thermal energy is negligible. Consequently, Eq. (7) has the character of a potential. The integral over the crack front a(s) in W f r is equal to the full crack area. G C is the critical energy release rate of the material. Integration over the full volume of the strain energy density 1 2 σ kl ε el i j leads to the elastic strain energy U el . The elastic strain energy U el is implicitly dependent on a(s). Using the Fourier series (Eq. (6)), the mechanical potential Π(a(s)), which is a functional, is thus converted into the function Π(a 0 , C i ). As is the length in the 3D case. The fracture mechanics length l 0 is a value that contrasts the fracture energy per crack area G = G c , necessary for crack propagation and the stored elastic energy per volume. According to Irwin (1958, p. 560), the mode I energy release rate for the 3D case can be calculated by l 0 and κ/ h are given by the mechanical and thermal material properties. The geometrical length b is given by half of the plate thickness. We consider the problem quasi-statically for a time δ = √ 4Dt. For this time, a 0 and L can be calculated using the two equations for stationary crack propagation G = G c and the bifurcation criterion, which will be discussed below.
In the following we give a short overview of the bifurcation analysis considerations according to Anderssohn et al. (2018).
For stationary crack propagation, it is necessary that is fulfilled.
The system tends to a minimum. This can be calculated by setting the first partial variation of the mechanical potential Eq. (7) to zero, i.e.
Due to the mechanical potential Π(a 0 , C i ), Eq. (14) can be rewritten as which is the fundamental solution of the problem. This solution leads to j + 1 equations. G c and l 0 are both known in Eq. (11) and if we assume that L is also given, we can use Eq. (15) to determine a 0 and C i at a time δ. Note that changes in the pattern perpendicular to the crack growth direction are excluded in the case of the fundamental solution. This is ensured by the fact that mode I is the only stress intensity factor acting on the entire crack front. The phenomenon of the merging of the three hexagonal columns, as described in Sect. 2.3, can be understood as a mathematical bifurcation problem. This means that there exists a second solution in addition to the fundamental one. Our goal is to find this bifurcation solution for a set of critical characteristic lengths. Starting from Eq. (15), we formed the second derivative, which leads to the Hessian matrix where C 0 = a 0 . The eigenvalue problem has a nontrivial solution when the coefficient determinant is zero This leads to j + 1 eigenvalues λ i . All eigenvalues have real values because of the symmetry of the Hessian matrix, which results in j + 1 real eigenvectors v i computed by (H − λ i I)v i = 0. If Eq. (18) results in positive eigenvalues, then the minimum of the mechanical potential Π (Eq. (7)) is found and the system is stable. This is the case for the fundamental solution, see Sect. 2.3. When one or more of the eigenvalues become negative, the system is unstable and bifurcation occurs. This leads to the merging of the three columns into one. The behavior after the bifurcation can be determined from the eigenvector v in conjunction with the fundamental solution a(s) (Nguyen 1987). A set of characteristic lengths is found where the smallest eigenvalue becomes zero, i.e.
Then, the bifurcation point is reached. The crack grows stationary as long as Eq. (13) is true. It stops growing when the condition (see also Bahr et al. 1988) is fulfilled. This condition states that for a fixed crack length a 0 the energy release rate G does not increase even though the temperature field (the cooling) continues to penetrate. If the condition is fulfilled also with the still further penetrating temperature field, it does not come again to the crack growth. Thus, the fixed crack length a 0 in Eq. (20) is the final crack length a end .

Finite element model
In this section we will first describe the general procedure how to determine a 0 /b and L/b from the FE analysis. In the following we will give more details about the used numerical methods and the influence of the number of Fourier series coefficients. As described in Sect. 2.4 we have six characteristic lengths of which four are given. The remaining two can be determined with the condition for stationary crack growth Eq. (13) and with the condition for the occurrence of a bifurcation point Eq. (19).
For a fixed hb/κ and δ/b, for different crack lengths a 0 /b and crack spacings L/b, the energy release rate G cF E M is calculated according to the Eq. (15) 1 and the smallest eigenvalue λ min is determined according to the Eqs. (16)-(18) by FE analysis. Due to the restriction to linear elastic material and small deformations, normalized material values E = 1, α = 1 and a normalized temperature difference ΔT = 1 can be used in the FE calculations. The dependence of ν will be discussed below. Now, by using G cF E M and the normalized E, α, ΔT and a given ν, the fracture mechanical length l 0F E M , according to Eq. (11), can be obtained. In a further step, the searched a 0 /b and L/b are now  19) is fulfilled. Furthermore, l 0Mat is calculated from a given material, e.g. in Table 1. The fracture mechanical length l 0F E M is then linearly interpolated so that it matches l 0Mat . This corresponds to the condition for the steadystate crack growth, Eq. (13). Thus, the crack length a 0 /b and crack spacing L/b were determined for a given δ/b and a given hb/κ. This process was then repeated for further δ/b and hb/κ.
Note that ∂G/∂δ is always determined with the FE analysis. During the evaluation it can be checked whether the condition for crack stop, Eq. (20), is true. Furthermore, Eq. (11) depends on the Poission's ratio ν. In Fig. 5 the normalised fracture mechanical length l 0 /b is plotted against ν. It can be seen that ν has a nonnegligible influence. The dependence of ν is also given in the analytical solution of the temperature-loaded free plate in Timoshenko and Goodier (2017). Therefore, the Poisson's ratio ν of the material must be taken into account in the FE analysis.
Geometric limits result from the variation of the crack space L/b = 0.05 . . . 1. The change of L causes a change of the model size and due to the meshing, the crack length is thus limited in the range of a 0 /b = 0.1 . . . 0.9. For some combinations of a 0 /b and L/b, crack closures occur when δ/b < 0.5. A more detailed investigation for the values a 0 /b < 0.1, L/b > 1 and δ/b < 0.5 is possible by adapting the FE mesh.
The determination of the derivatives in the fundamental solution Eq. (15), the bifurcation solution Eq. (19), and the change of the energy release rate over δ/b Eq. (20) was performed with the FDM analogous to Anderssohn et al. (2018). The central difference quo-tient with step size h a for crack length and h δ for time was used for this purpose.
Isoparametric hexahedral elements with quadratic shape functions were chosen for the FE mesh. To guarantee the independence of the step size h a in the FDM for the first and second derivatives, small step sizes are required compared to the crack spacing. This means that the strain energies must be determined as accurately as possible. Especially at the crack front, increased stresses and strains occur due to the singularity. Therefore, the crack tip was meshed with singular elements (Barsoum 1976) and a much finer discretisation was carried out in the vicinity of the crack tip (see Fig. 1). It was checked by a convergence analysis that the FE mesh is fine enough and the step sizes h a /b and h δ /b are sufficiently small for the FDM. Fig. 6 shows the results for the fundamental solution Eq. (15) 1 using the twelfth hexagonal model and Fig. 7a shows the results for λ min (Eqs. (16)-(18)) using the half hexagonal model. Differences in the computation of l 0 /b and λ min are not recognized for the presented range of a 0 /b with increasingly finer FE meshes and shorter step sizes. For the calculation of λ min , the second derivative of the strain energy is needed according to Eq. (16). Thus, λ min is more sensitive to numerical inaccuracies compared to l 0 /b. Therefore, for λ min the range in which the abscissa is cut (in which then also the bifurcation point occurs according to Eq. (19)), was shown detailed in Fig. 7b. Because of the many possible configurations for L/b, δ/b, and hb/κ, the following settings were chosen for the further analyses for certainty: 65000 nodes for the twelfth hexagonal model, 380000 nodes for the half hexagonal model and h a /b = 5·10 −4 for the step size of the crack length. For the choice of h δ /b, a similar consideration was made for ∂G/∂δ as for λ min shown in Fig. 8. Again, the differences due to various step sizes h δ /b are negligible. Based on Fig. 8b, the step size h δ /b = 5 · 10 −4 was chosen for the further calculations.
The Hessian matrix Eq. (16) was passed to the mathematical program Matlab, where the eigenvalues λ i were calculated.
The 1D temperature field Eq. (1) was loaded into the nodes of the 3D FE mesh using a special routine.
To determine C 4 /b (according to Eq. (15) 2 ), five FEM calculations were performed with given crack length a 0 /b and given crack spacing L/b using the twelfth column model (fundamental solution). From these calculations with different given values for C 4 /b, Fig. 6 Results of the convergence test calculations l 0 /b according to Eq. (15) 1 with Eq. (11) for given L/b = 0.5, δ/b = 1, ν = 0.22 and hb/κ = 1.0833, testing three different mesh qualities (n. mean number of FE nodes) and three different step sizes h a five values for the elastic strain energy Eq. (8) are obtained. According to the principle of minimum potential energy, the correct value for C 4 /b is the minimum of the least squares interpolation function. As an example, the calculated minimum and maximum values of C 4 /b with the corresponding settings are given in Table 2.
The influence of the number of coefficients j in the Fourier series Eq. (6) is shown in Fig. 9. The calculations were performed for the Al 2 O 3 99.8% plate (Table 1) with h = 10 4 W/m 2 K resulting in hb/κ = 1.0833 and ΔT = 800K . The crack spacing L/b versus the crack length a 0 /b is presented in Fig. 9a. In Fig. 9b, the time variation of the energy release rate ∂G/∂δ versus the crack length a 0 /b is shown. As can be seen, the difference in the results is low. This is in good agreement with the results by Anderssohn et al. (2018). For example, in Fig. 9b, where the curves cross the abscissa, the final crack length a end /b is readable. The relative difference between j = 1 (a end,C1 /b = 0.4839) and j = 4 (a end,C4 /b = 0.4936) is about 2%. To save computational time, j = 1 (a 0 and C 1 ) was chosen for further calculations.

Crack pattern in thermo-shocked free alumina plate
It is assumed that an alumina plate is heated to a temperature T 0 . Then the plate is quenched in a cooling Table 2 Calculated minimum and maximum values of C 4 /b with the corresponding settings for hb/κ = 1.0833 and ν = 0.23  (19), with j = 1 (a 0 and C 1 ) for given L/b = 0.5, δ/b = 1, ν = 0.22 and hb/κ = 1.0833, testing three different mesh qualities (n. mean number of FE nodes) and three different step sizes h a , b detailed plot at the intersection with the abscissa  (Table 1) with h = 10 4 Wm −2 K −1 → hb/κ = 1.0833 and ΔT = 800 K) liquid of constant temperature T 1 . At a sufficient distance from the interface, the temperature field (Eq. (1)) holds in the plate. Table 1 shows the material, thermal and geometrical properties of alumina plates.
Extensive simulations were performed for the Al 2 O 3 99.7% plate ( Table 1). The results are illustrated in Figs. 10 and 11.
If a end /b or L end /b or better both are known from experiments and all parameters listed in Table 1 are also known, it is generally possible to determine the heat transfer coefficient h from diagrams like Fig. 10 between any brittle material and the cooling medium.
Here, the value hb/κ = 1 represents a kind of transition point. For hb/κ < 1, small changes of hb/κ have a very strong influence on the temperature gradient and thus on the stress. The final crack length achieved (Fig. 10a) and the final crack spacing (Fig. 10b) are very different. For hb/κ > 1, only large changes in hb/κ have a noticeable effect on a end /b and L end /b. If hb/κ > 2 the values of a end /b and L end /b very quickly approach the values of the asymptote at hb/κ = ∞. This means for the alumina plate 99.7% (Table 1) that the dependence in the final results is small for values of the heat transfer coefficient in the range h = 16,000Wm −2 K −1 to h = ∞ Wm −2 K −1 . This could explain why in the literature, e.g. Hasselman (1970), Jiang et al. (2012), Li et al. (2013), and Xu et al. (2016a), the heat transfer coefficient h between water and ceramics was determined to vary over a large range.
The progression of L/b over a 0 /b for hb/κ = 1.0833 and hb/κ = ∞ as significant examples is shown in Fig. 11. As can be seen, L/b decreases with increasing ΔT . This could also be observed in the experiments, see Figs. 12, 13, and 14, and by Bahr et al. (1986, Fig . 2) and Xu et al. (2016b, Fig . 2) on the outside of the quenched plates. Furthermore, in Fig. 11a it can be seen that the value of the crack spacing L/b hardly changes with the growth of the crack a 0 /b. This indicates that there is no merging of three columns as described in Sect. 2.3 and that the columns therefore grow in a steady course (fundamental solution). The calculated bifurcation points are not valid because the fixed point 0 in the model (Fig. 3) wants to catch up again. A postcritical analysis could show this (see Bahr et al. 2009, p. 4). However, for hb/κ = ∞ in Fig. 11b the crack space L/b is almost tripled for ΔT = 400 K as well as for ΔT = 450 K. For the other ΔT , L/b changes also significantly with increasing crack length. If the perturbation of the crack pattern is small enough, this will show the same trend as in Fig. 3 and the new diameter should be close to √ 3L. With the increase of ΔT the tensile stress increases according to Eq. (3). This leads to more cracks and thus to a decrease in the crack spacing L/b. This was also found out in 2D by thermal shock experiments with 50 × 10 × 1 mm 3 thin ceramic specimens in Jiang et al. (2012 , Table 1).
Consequently, the competition between short cracks increases. Short cracks have a smaller crack spacing while long cracks have a larger one.

Comparison with 2D model and experiments
To validate the developed model, analyses were also performed with a comparable 2D FE model. Furthermore, alumina plates were thermally quenched for different ΔT and evaluated. In Fig. 12 the results from the thermal shock experiment, the 2D and the 3D FEM bifurcation analysis are plotted. First, a brief explanation of the experiments and their evaluation is given. This is followed by a description of the 2D model with a comparison to the data from the 3D model and the experiments. At the end of this section, possible reasons why the 2D model fits the experimental data better than the 3D model are discussed.
The plates for the experiments were made of DOCE-RAM A-132 Al 2 O 3 99.7% with the properties given in Table 1 with dimensions of 40 × 40 mm 2 . The plates were heated to T = 550 • C, 750 • C, 950 • C and quenched in boiled water of T = 100 • C. The quenched temperatures were thus ΔT = 450 • C, 650 • C, 850 • C. Then the plates were treated with a contrast agent (Schilling et al. 2005) to achieve a better separation of ceramic and crack for the evaluation. Next, the plates were scanned with CT in High Aspect Ratio Tomography mode.
Through the CT scan and the contrast agent, it is possible to produce images of the surface of the sample, but also of the deep interior of the sample with the cracks highlighted, see Figs. 13 and 14.
The crack length a 0 /b and the crack spacing L/b of the thermal shock experiment were measured from Fig. 14. The procedure is exemplified by Fig. 14 (e). The sectional view was loaded into the programme Engauge Digitizer version 10.10. Based on the scale line, the thickness of the sample was determined (2b = 6.97 mm). Then horizontal lines were inserted from the top to the middle. These lines start and end at the outer cracks. The intersections of the lines with the cracks were counted. The crack widths are given by L =(width of the line)/(number of intersections -1). The associated crack length a 0 is the z-position of the line. Then the ratios L/b and a 0 /b can be determined with the measured b. All information is collected in Table 3. The 2D FE model is described in Bahr et al. (1988). For the analysis the bifurcation criterion from Nemat-Nasser et al. (1980) and Bahr et al. (1992) was used. In order to compare the results with the 3D model, the same material, thermal and geometrical settings (see Table 1 Al 2 O 3 99.7%) were used in the 2D model as in the 3D model. For the 2D model, the same temperature field as in Eq. (1) was applied.
The results for the final crack length a end /b with the corresponding final penetration depth of the temperature field δ end /b and the final crack space L end /b in the case of ideal cooling (hb/κ = ∞) are presented in Fig. 15 for the 2D and 3D FE model. The values of δ end /b match quite well. The achieved values for a end /b are slightly larger in the 2D FE model than in the 3D FE model. The difference for L end /b is slightly larger. A factor of 2 or less between the results of the 2D FE model and those of the 3D FE model agrees well with the findings in Bahr et al. (2009). This kind of difference in a end /b and L end /b was also found out in the experiments by Shao et al. (2010). In our investigations, similar results were obtained for other values of hb/κ, as shown in Fig. 12.
As shown in Fig. 12, the results of the 2D analysis between hb/κ = 0.5416 and hb/κ = 1.0833 agree well with the measured values of the samples. Thus, a heat transfer coefficient between h = 4.3 · 10 3 Wm −2 K −1 and h = 8.7 · 10 3 Wm −2 K −1 can be determined. This is slightly lower than the measured h in Zhou et al. (2012). This is due to the fact that nucleate boiling occurs at low sample temperatures due to the lower water temperature in Zhou et al. (2012) of 20 • C, thus increasing the heat transfer coefficient. Quenching in boiling water favours the formation of a vapour film, which decreases the heat transfer coefficient (Singh et al. 1981).
It is the question of why the 2D bifurcation analysis seems to predict the crack width L/b better than the 3D one. Here are two possible explanations.
First, due to the fine grain of the material, the cracks can close again after unloading. Thus, it is possible that the contrast agent does not reach every crack and thus the cracks or their lengths are partially poorly detectable. For more details see Zielke et al. (2021).
Second, the 2D and 3D models do not predict cracks going through the entire body, see Fig. 14. In Fig. 13, especially in Fig. 13b, e, h, it can be seen that the cracks start from the edge and go into the whole body. It can be assumed that there is a surplus of energy at the edge  (Table 1) for different values of hb/κ and the quenching temperatures ΔT of 400 K, 450 K, 650 K, 850 K and 1050 K Fig. 11 Normalised crack space L/b versus normalised crack length a 0 /b (to crack stop) for a hb/κ = 1.0833 and b hb/κ = ∞ for the material parameters of alumina plate 99.7% (Table 1) due to the high temperature gradient. This causes the cracks to grow unstably into the body. The fact that the cracks also branch out is an indication of this (Kanninen and Popelar 1985, pp. 205-207). These cracks relieve the body and convert the 3D stress state in a 2D stress state.
Thus, there is no longer a closed plate when the cooling penetrates further from the surface. The penetrating cracks have cut the plate into many small slices. As has been shown by Shao et al. (2010, Figs. 3e and f), when comparing a whole plate and stacked individual plates at high temperature differences ΔT , the crack pattern on the cooled surface looks similar. So, even for the stacked individual plates, the crack pattern is not invariant in the x-or y-direction (see also Bahr et al. 1986). However, between the interior surfaces of the stacked single plates and the cut surfaces at the whole plates, the patterns differ. Thus, more cracks are found in the whole plates, and thus a shorter crack spacing L/b, than in the interior surfaces of the stacked individual plates. This is not addressed in Shao et al. (2010). Only the fact that crack branching occurs in the stacked single plates and not in the whole plate is explained by a different stress field caused by edge effects in the stacked single plates.
One reason or a combination of both reasons may explain why the 2D bifurcation analysis fits well with the measured data from the experiments.
Nevertheless, in the following aspect the data from the 3D model, as well as the 2D model, match the experiments. Both models predict that the crack spacing will change only slightly. This means that the cracks grow stable and do not recede. This agrees well with the measured data shown in Fig. 12, although it should be noted that the increase in experimental data in L/b is due to the cracks penetrating from the edge.

Conclusion and future work
The model developed by Bahr et al. (2009) was applied and adapted to predict the complex growth of the hierarchical 3D crack pattern of a brittle plate under thermal shock conditions without external constraints. This model is based on the assumption that the cracks form ideal hexagonal columns and that in this idealised crack pattern, the merging of three columns into a larger column is the mechanism for increasing the crack spacing. For cooling penetration, it was assumed that a 1D stationary temperature field (Eq. (1)) exists far enough away from the edges. Based on these assumptions, a 3D FEM bifurcation analysis with the Fourier expansion of the crack front (Anderssohn et al. 2018) was performed. Besides the different type of temperature field, the following differences and extensions are present in the model used in this work compared to Bahr et al. (2009) and Anderssohn et al. (2018).
The model no longer has infinite dimensions in the z-direction but the finite thickness b. The temporal penetration of the symmetrical temperature field Eq. (1), the associated reduction of the temperature gradient and thus the reduction of the mechanical stress, leads to crack arrest and thus to a final crack length a end . The final crack length was determined by evaluating the change in the energy release rate as a function of time (Eq. (20)). a end is the measure of the damage to the structure and can be used later for a residual strength analysis. Because the plate can contract without hindrance, freedom from forces in the lateral direction of the model is required.
Through a convergence study, which includes the optimisation of the FE mesh and the influence of the FDM step size, the computational effort could be reduced. Furthermore, it was shown that the results converge rapidly for an increasing number of Fourier coefficients. To reduce the computational effort even further, the different calculations were carried out with only one Fourier coefficient ( j = 1).
Due to the linear nature of the problem and the use of the characteristic lengths, describing the material and the thermal load, it was effectively possible to determine the evolution of the crack spacing L/b over the crack length a 0 /b by a parametric analysis.
The 3D model was verified by a comparison with an existing 2D model by Bahr et al. (1988). As in Bahr et al. (2009), the expected difference between 2D and 3D FE model in L/b was smaller than a factor of 2. Due to the intense cooling at the edge of the specimens and the associated formation of cracks through the whole specimen, a validation of the 3D model by experimental data was unfortunately not possible. In further research, additional experimental investigations are planned in which the specimen will be isolated at the edge. It should be noted that due to the evaluation by CT scan, the dimensions of the specimens are limited and specimens with larger dimensions cannot be simply used.
Further investigations with the presented 3D model are possible. For example, by modifying the FE mesh towards very short cracks, the critical temperature difference ΔT c for thermal shock could be analysed at which no crack growth occurs. It should be noted that a good knowledge of the mechanical and thermal material data is important for an accurate analysis. Conversely, it is possible to determine the critical stress intensity factor K I C and the heat transfer coefficient h using 3D FEM bifurcation analysis. Both values are difficult to determine using other methods. For this, the curves of a 0 /b and L/b must be known from thermal shock experiments.   Fig. 14 a is I, b is II, c is III, d is IV, e is V and f is VI  Fig. 13 b, e, h Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Cre-ative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.