Numerical 3D‐bifurcation analysis of star‐shaped crack patterns using the energy method

The present research deals with a three‐dimensional (3D) Finite Element Method (FEM) bifurcation analysis based on the global mechanical potential that can be used to find the parameters at which a crack pattern changes. In our case we want to analyze at which point the star‐shaped shrinkage cracks in an aqueous colloidal suspension filled in a glass cylinder change from four to three or two cracks growing. The driving force for the crack growth is shrinkage caused by diffusion controlled drying. The 3D crack front geometry is described efficiently by using a Fourier series approach. Based on steady‐state crack growth, the Fourier coefficients are determined in a first step using an optimization algorithm. As a result, the time dependent crack growth can be determined. In a second step, the bifurcation point is determined by an eigenvalue analysis of the second order derivatives of the potential energy of the system. If the lowest eigenvalue reaches zero the fundamental solution becomes unstable and a transition will occur. Our analysis shows that the transition from four to two cracks is preferred over the transition from four to three cracks.


INTRODUCTION
Crack patterns due to in-homogeneous shrinkage caused by either drying or thermal shock often occur in brittle materials such as wood, clay, concrete, glass or technical ceramics.The resulting crack patterns depend on the material and process and on the boundary conditions.While this phenomenon is fascinating in nature as the formation of basalt columns shows, it might become critical in engineering applications as for example, the use of ceramics as heat shield in high temperature environments like gas turbines.Therefore, a deep understanding of how such crack patterns evolve is crucial to minimize or to control cracking and thus improve the reliability.
During the penetration of the crack pattern into the material one can distinguish between regions, where the crack pattern grows mainly self-similar and regions where a transition of the crack pattern occurs mostly by leaving parts of the crack front behind.This causes the crack spacing to change.This has been studied for crack patterns which have a two-dimensional (2D) nature for example, in ref. [1] for Glass-quartz plates for thermal shock experiments and for drying caused crack patterns in ref. [14].Research on 3D crack patterns were carried out in refs.[2,3,12,13].The crack spacing depends on parameters like shrinking, velocity or temperature gradient and can be determined by bifurcation analysis [2,4].
F I G U R E 1 (A) Transition from four to three cracks (from ref. [6]), (B) scheme of bifurcation point in the tube and (C) the moisture concentration field (Equation (1)) for different penetration depths   ∕.
Unidirectional drying experiments were performed in refs.[5,6].In these, glass cylinders (diameter  ≈ 1 mm, length  = 12 cm) were filled with an aqueous colloidal suspension.The cylinders were closed on one side to ensure unidirectional drying.The suspension was dried by placing the bottom of the cylinders in desiccant [5] or by placing the cylinder in a drying chamber [6].Drying leads to shrinkage and thus to tensile stresses in the suspension.The tensile stresses become so high that cracks develop.In the experiments, it was observed that four cracks grow.These cracks arrange themselves star-like in the cylinder.Then, as the cracks continue to grow, a transition occurs where one crack [6] (see Figure 1A and B) or even two cracks [5] arrest.Until now, despite its 3D nature, the phenomenon has only been considered using 2D models [7].We now want to present a 3D Finite Element Method (FEM) bifurcation analysis to find the transition point.

Moisture concentration field
In refs.[5][6][7] it is assumed that drying, and thus shrinkage, is driven by capillary pressure and governed by Darcy's law.This is assumed, for example, in ref. [7], from a distance from the bottom of the cylinder of  ≈ 40, see Figure 1B.Here, however, the behavior of the cracks from  = 0 is to be investigated.Therefore, it is assumed that the shrinkage occurs due to an in-homogeneous drying field.To describe this field, the analogy between diffusion controlled drying and heat conduction via thermal shrinkage is used.Assuming that the crack surfaces have only a small influence on the drying, the moisture concentration field depends only on the spatial coordinate  and the time dependent penetration depth   ().Thus, the concentration (,   ) can be described similar to the well-known half-space field with ideal cooling [4] (,   ) =  0 − ( 0 −   )erfc(∕  ), see Figure 1C.Here   denotes the drying diffusion coefficient,  0 the reference concentration at the beginning of drying and   the final concentration.

Fundamental elastic equations
Because of the slow drying and small dimensions, acceleration terms and volume forces can be neglected, so the local mechanical equilibrium can be written as where   are the mechanical stresses and ( ),  = ( )∕  stands for the partial derivative in   -direction.Despite the partially visco-elastic behavior [5], linear-elastic material behavior with small deformations is assumed here for simplicity.
The strains are thus given by where   , ,  and   label displacements, Young's modulus, Poisson's ratio and Kronecker delta.Furthermore, the first and second terms in Equation (3) correspond to the elastic strain and the third term to the strain due to shrinkage.The last term is analogous to thermal shrinkage.

3D fracture mechanics with series expansion for the crack front
The global mechanical potential without external work is given by [3] The first term represents the elastic strain energy   .The second term is the fracture energy   consisting of the critical energy release rate   and the fracture area   = ∫ () .Where () is the crack length, which depends on the path coordinate  along the crack front as shown in Figure 1B.Now, according to ref. [2], the geometry of the crack front can be described by a truncated Fourier series Here (, ) is measured from the bottom of the cylinder, see Figure 2A, and  0 and   are the time-dependent Fourier coefficients.Steady-state crack growth is assumed.This means that the crack growth condition  =   = .is valid at the whole crack front.Furthermore, only mode I is present at the crack tip.Now, in a first step, for a given  0 and a given time , or a given   , the   must be determined to describe the crack front geometry together with the energy release rate .
This can be done by means of the principle of the minimum of the potential energy Π∕ = 0.This is the first variation of Equation ( 4) and is called the fundamental solution.By using Equation ( 5) the mechanical potential is no longer a functional Π(()) but a function Π( 0 ,   ).Thus, the partial derivative ) can be used.This leads to  + 1 equations.From Π∕ 0 = 0,  =   is obtained for the given  0 and   .By means of the remaining  equations the   are determined.This can be done using an optimization algorithm.

Bifurcation analysis
After determining the Fourier coefficients   of the fundamental solution, the stability of the solution is analyzed.If all eigenvalues of the matrix of the second order derivatives of the potential energy are positive; the system is stable and the crack pattern will grow in a self-similar manner.If the lowest eigenvalue becomes zero a bifurcation of the solution can occur and the crack pattern will change.For this the Hessian matrix  with has to be determined.The indices  and  represent the number of the crack.The cracks are numbered in Roman numerals.Here because of the investigation that four cracks grow and because of the use of symmetry, only two cracks need to be varied.Thus  = .For more information, see Section 2.5 and Figure 4C.By solving the eigenvalue problem the eigenvalues  are obtained.Here,  is the unit matrix,  the eigenvector and  the zero vector.The bifurcation point is characterized by the fact that the lowest eigenvalue becomes zero, that is, At this point (with given  0 and   ) beside the fundamental solution at least another solution with the same mechanical potential exists.In other words, the curvature of the mechanical potential Equation ( 4) is zero.On the other hand, if   > 0, there is a stable state and if   > 0, there is an unstable state.Where stable condition means that only the fundamental solution Equation ( 6) is valid, so here all four cracks grow.Unstable state means that there is a solution with a lower mechanical potential than the one determined by the fundamental solution.

Numerical realization
The four characteristic lengths are ,   ,  0 and the Griffith length  0 .Here  0 =   ∕(Δ) 2 .For further information see refs.[2,3,6]. is known.Since the energy release rate is an outcome of the used Finite Element (FE) model and due to the linear nature of the problem, we change the variables as explained in refs.[2,3].Thus, for different values of   ∕ and  0 ∕, the corresponding   is calculated using Equation ( 6),  = 0, and also   using Equations ( 7)-( 8) by parametric analysis.Then, for a given material   , the crack growth curve  0 ∕(  ∕) can be determined by linear interpolation.The crack length at which bifurcation occurs is then determined by the condition Equation (9).To solve the Equations ( 6)-( 9) the strain energy and its partial derivatives are required.For this Equations ( 2) and (3) together with Equation (1) as a loading term together with necessary boundary conditions will be solved by using the commercial FE software ANSYS (Ansys Parametric Design Language (APDL)) to get the strain energy   in Equation ( 4) as a result.The partial derivatives were determined for a constant   by changing the   using the Finite Difference Method (FDM), see also ref. [2].The eigenvalues from Equation (8) were determined by using MATLAB.
The FE mesh for the calculation of the fundamental solution is shown with the straight crack front in Figure 2A.Isoparametric hexahedral elements with quadratic shape functions were chosen for the FE mesh.At the crack front, singular elements [8] were used to represent the increased stresses and strains.Due to the symmetry, only an eighth of the cylinder is modeled for the growth of all four cracks.Measured from the ground,  0 is then shown up to the crack front in Figure 2A.The model is assumed to extend infinitely in -direction in the FE model, this was realized as a long cylinder with length 24.It is assumed that the suspension mass adheres perfectly to the glass cylinder.This is ensured by the boundary condition that all displacements at  =  are prevented.This can be seen well in Figure 2C because at the boundary  =  the displacements are zero.The Broyden-Fletcher-Goldfarb-Shanno (BFGS) method was implemented via APDL as the optimization algorithm for the iteration of the fundamental solution.Starting from the straight crack front, the curved crack front was iterated with the BFGS method by a mesh distortion algorithm, see Figures 2B and C.
The material parameters were chosen to be  = 1 MPa and  = 0.3 in accordance to ref. [7].The cylinder radius was selected as  = 1 mm.
For the determination of the bifurcation solution, where one or two cracks arrest,  0 =  0 was varied.The determination of the Hessian matrix is explained for the case that one crack arrests.In Figure 4C, the cracks that continue to grow are highlighted in green and the one that arrests is highlighted in red.Thus, according to the FDM, the three green highlighted cracks are varied in the FE model so that the second derivative  2 Π∕ 2 0 is obtained. 0 stands for the crack length of the green cracks.Then the one red crack is varied according to the FDM and  2 Π∕ 2 0 is calculated.Here,  0 is the crack length of the red highlighted crack.For the mixed derivative, all four crack lengths must then be varied according to the FDM to obtain  2 Π∕ 0  0 .The symmetric Hessian matrix is then given by The procedure is carried out analogously for the case that two cracks arrest.By using the symmetry, only a quarter model had to be created for the case that two cracks arrest and a half model of the cylinder for the case that one crack arrests.
In the quarter model and half model, two and three cracks, respectively, meet in the center.For the determination of the Hessian matrix by FDM, the cracks have a different position due to the change of  0 .In order to be able to merge the FE meshes, three crack tips had to be modeled, which have a distance corresponding to the selected step size in the FDM.

Fundamental solution
As the number of Fourier coefficients  increases, the   should change only slightly.As their influence becomes smaller and smaller, the higher   should approach zero in value.This is not the situation, as can be seen in Table 1.The larger the  is, the larger the   become.Also, the   seem to erase each other by changing the sign.So, it can be said that the Fourier series does not converge.For this reason, the local energy release rates  along the crack front were additionally calculated.This was done by evaluating the stress intensity factor   for plane strain from the crack tip displacement.The energy release rate can then be calculated as  =  2  (1 −  2 )∕.Further information can be found in ref. [9]. along the crack front is shown in Figure 3A for  0 ∕ = 1 and   ∕ = 1.Additionally, the corresponding crack front contours () are shown in Figure 3B.In Figure 3A it can be seen that at least four ( = 4) or better five ( = 5) Fourier coefficients are necessary to fulfill the criterion for steady-state crack growth ( =   in each point at the crack tip).The convergence of the Fourier series with higher  is not yet clarified conclusively.This is because the different singularities appearing at the edges cause convergence problems in the search for the Fourier coefficients.The curved crack, as can be seen in Figures 2B and 3B, approaches for ∕ = 1 the glass cylinder.There the crack is at an interface of two materials with different elastic properties.This changes the nature of the stress singularity that is referred as the Zak-Williams singularity [10].By suppressing the displacements at this edge ( = ) in the model (as explained in Section 2.5), this would correspond to a Young's modulus of  = ∞ for the glass cylinder, which means there is no energy to release.Thus  = 0 should hold for the energy release rate at the glass cylinder (see also ref. [11] p. 133).This is also fulfilled in the FE model as can be seen in Figure 3A on the right hand side at ∕ = 1.
On the left side (∕ = 0) in Figure 3A the effects of a corner singularity can be seen.Here in the center of the cylinder the cracks influence each other.Furthermore, the last singular elements have an angle of 45 degrees which has a negative effect on the local evaluation of .
But still one can see in Figure 3 that in general the difference in contour with increasing  decreases and the differences in energy release rate overall become smaller.

Bifurcation analysis
Now the stability of the solution is analyzed by Equations ( 7)- (9).For the fundamental solution five Fourier coefficients ( = 5) were always used.
In Figure 4A and B,   is plotted versus   ∕ for the two assumed cases of the arrest of one crack (blue line) or of two cracks (red line).The bifurcation points Equation ( 9) for the case where  0 ∕ = 1 can be read from Figure 4A with   ∕ ≈ 0.7 for the arrest of two cracks and with   ∕ ≈ 1 for the arrest of one crack.The case, where two cracks arrest occurs before the case where one crack arrest.For longer crack lengths than  0 ∕ = 1 currently no point could be determined with   = 0.The required   ∕ led to such small gradients Π∕  of the energy that the BFGS method could not find a minimum anymore and therefore no Fourier coefficients could be calculated with the fundamental solution.Due to the numerical methods (FEM, FDM) with which Π∕  is determined, the numerical errors have a larger effect.Using more elements in the FE model and a different step size in the FDM did not solve the problem.The use of a gradient free optimization algorithm, such as the Nelder-Mead method, should remedy this.
For a   = 4.77 ⋅ 10 −2 N/mm the crack growth curve could be generated for the case where two cracks stop.To do this,   was determined for intermediate values of  0 ∕ and   ∕, and then linear interpolation was performed.The crack growth curve is presented in Figure 4C.It shows again that four cracks grow up to the point   ∕ ≈ 0.7 and  0 ∕ = 1.Then two cracks arrest.The continuing blue line represents the unstable region and does usually not occur in this way.The case that a transition, where one crack arrests would occur on this blue line.But the arrest of two cracks happens before.This is in contradiction to the experiments in ref. [6].The fact that one crack arrests in the experiment could be explained by the fact that there is an imperfection influence, such as the four cracks not meeting in the center of the cylinder.Furthermore, an explanation is that in the experiment not only diffusion controlled drying occurs, but also capillary effects could play a role [6,7].In ref. [5], the transition from four cracks to two was observed.Two cracks arrest and the other two continued to grow slowly.Then, at some point, the two arrested cracks caught up very quickly and overtook the slow-growing cracks.Then the two slow growing cracks arrest and the cracks exchanged their growth behavior.This cyclic crack growth behavior was explained by visco-elastic effects.

CONCLUSION
The arrest of one or two cracks starting from growing four cracks can be found for a given parameter set of  and   using the given method.The arrest of two cracks occurs before the arrest of one crack.This is in contrary to the experiment from ref. [6], but could be explained by imperfection influence which is difficult to represent in the FE model.Alternate crack growth was observed in the experiments in ref. [5].
With the global approach of the mechanical potential Π and the approximation of the crack front geometry by a Fourier series expansion, the condition for steady-state crack growth  =   can be well satisfied on average.However, the Zak-Williams and the corner singularity at the edges hinder the convergence.But with a local approach, in which an attempt is made to satisfy the condition  =   at each element node in the FE model (see ref. [13]), the two edges would also cause considerable problems.For small   the BFGS method has not found a minimum and therefore no Fourier coefficients.The described problems can be solved by using another optimization algorithm (gradient free) like the Nelder-Mead method.

A C K N O W L E D G M E N T S
The authors would like to acknowledge the financial support of the Deutsche Forschungsgemeinschaft (DFG) under Grant No. HO 4791/1-2.The authors are grateful to the Center for Information Services and High Performance Computing [Zentrum für Informationsdienste und Hochleistungsrechnen (ZIH)] at TU Dresden for providing its facilities for high throughput calculations.
Open access funding enabled and organized by Projekt DEAL.

F
I G U R E 2 (A) FE model with straight crack front for fundamental solution, (B) result of mesh distortion for modeling of the curved crack front and (C) illustration displacements as vector sum for  0 ∕ = 1,   ∕ = 1.5,  = 5.FE, Finite Element.

F
I G U R E 4 (A) lowest eigenvalue   versus   ∕ for  = 5,  0 ∕ = 1, (B)  0 = 8 and (C) crack growth curve for   = 4.77 ⋅ 10 −2 N/mm (the blue line results from the case that one crack arrests and the red that two cracks arrest).
TA B L E 1