Dynamic propagation of mode III cracks in a Lattice Boltzmann method for solids

This work presents concepts and algorithms for the simulation of dynamic fractures with a Lattice Boltzmann method (LBM) for linear elastic solids. This LBM has been presented previously and solves the wave equation, which is interpreted as the governing equation for antiplane shear deformation. Besides the steady growth of a crack at a prescribed crack velocity, a fracture criterion based on stress intensity factors has been implemented. This is the first time that crack propagation with a mechanically relevant criterion is regarded in the context of LBMs. Numerical results are examined to validate the proposed method. The concepts of crack propagation introduced here are not limited to mode III cracks or the simplified deformation assumption of antiplane shear. By introducing a rather simple processing step into the existing LBM at the level of individual lattice sites, the overall performance of the LBM is maintained. Our findings underline the validity of the LBM as a numerical tool to simulate solids in general as well as dynamic fractures in particular.


Introduction
Regarding solid bodies and structures, not only the deformation under external loads is of interest.Different fields of science and engineering, such as geophysics or civil engineering, are also concerned with the mechanism of fracture and the study thereof.Thus, a number of different numerical techniques for the simulation of dynamic crack propagation and other fracture related phenomena have emerged.More prominent among these are the finite element method (FEM) [1,18] the boundary element method [7,31], or more recently peridynamics [26,10] and phase field methods [4,25].
Lattice Boltzmann methods (LBM) [12] are another approach to simulations in engineering and are now widely used in computational fluid dynamics.The principle, however, can be adapted to different problems in various disciplines of science [28,27].The usage of LBM in computational solid mechanics is recently developed [16,21,19,8] and the application to problems involving fractures has been considered in [5,22,23].LBMs are mesoscale methods, employing elements of statistical mechanics, and work with a rather simple transport mechanism on a regular grid.This promises good computational efficiency, especially when the algorithm is parallelized [14,17].This work regards the dynamics of antiplane shear, which reduces the Navier-Cauchy equation to a 2D wave equation.Different LBMs [5,30,9] have been proposed to solve the problem of wave propagation.Chopard and Luthi [6] even proposed an LBM to simulate crack growth in a simplified solid.However, this model is not consistent with linear elastic solid mechanics, nor does it include a fracture criterion based on theories of classical fracture mechanics.
In a previous work [22], we have shown the LBM for antiplane shear, based on the formulation of Yan's [30] LBM for waves, and applied this to a stationary crack with mode III opening.We then improved thereon with the introduction of non-lattice conforming boundary conditions [24], which are highly relevant for the utilization of the LBM in the field of fracture mechanics.2018) [22] .The outline of the body in the deformed configuration is indicated by dashed lines.
Now we directly build on this by introducing the extension to dynamic crack propagation.It is our intention to demonstrate that crack growth can be modeled with the LBM in a way that is in agreement with classical fracture mechanics.However, the comparison to competing approaches such as finite element simulations of phase field models for fracture [13,3] and the X-FEM [18] is beyond the scope of this work.In essence, handling dynamic cracks is reduced to a post-processing of lattice sites, without a need for remeshing due to the fixed grid.The criterion for crack propagation that we employ is based on the concept of stress intensity factors (SIFs) by Irwin [11].The SIF characterizes the load on the crack tip and is evaluated via the elastic fields in the vicinity of the crack tip in our approach.
Our work is structured as follows.First a short background on continuum and fracture mechanics is given in Sec. 2, followed by a review of the LBM for waves in Sec. 3. Next, in Sec. 4, we present the algorithm and discuss details regarding the implementation of crack propagation and fracture criteria.In Sec. 5 numerical results are shown for two problems.In order to demonstrate the agreement with classical fracture mechanics and to increase the reliability of our approach, we deliberately chose relatively simple problems for which an analytical benchmark solution exists.In the first example we do not employ any fracture criterion to decide whether a crack will propagate, but simply force the crack to grow at a certain speed.For this case, we demonstrate that the LBM is able to predict the elastic fields surrounding a moving crack tip, i.e. the SIF, in agreement with analytic solutions.
In the second numerical example, the implementation of the fracture criterion is tested and it is shown that, in our algorithm, a crack propagates if the fracture criterion is fulfilled.
Lastly, Sec. 6 summarizes the results.It also discusses the algorithm and its generalization to different LB schemes, such as plane strain, or different propagation criteria.

Mechanics of Linear Elastic Solids and Fractures
This section summarizes the concepts of solid mechanics and fracture mechanics that the proposed new LBM builds upon.

Antiplane Shear Deformation of Linear Elastic Solids
A linear elastic body with shear modulus µ in a Cartesian x, y, z-coordinate system is considered, for which the displacement field is restricted to u = w(x, y)e z .Here, x/y are the in-plane coordinates and e z is the unit vector in the out-of-plane direction, see Fig. 1.This deformation, in which the out-of-plane displacement is the only nonzero displacement component and is a function of the in-plane coordinates, is commonly referred to as antiplane shear deformation.For a linear elastic body with density ρ, the governing equation of the displacement field is the wave equation

Dynamic Linear Elastic Fracture Mechanics
For sufficiently brittle materials, the debonding of the material during fracturing is determined by the mechanical fields directly in the vicinity of the so-called process zone in which material separation actually takes place.It is commonly assumed that in such a case, the fields surrounding the process zone can still be determined with sufficient accuracy by the theory of linear elasticity.Thus, fracture criteria are constructed by characterizing the elastic fields in the vicinity of crack tips and comparing these to critical values, which need to be determined in experiments for a given crack loading mode and material.
In this work, attention is restricted to one such fracture criterion in order to model crack growth in brittle materials.The elastic fields surrounding the crack tip are obtained by the LBM and are fed into the fracture criterion that eventually determines whether a crack propagates.
The criterion employed in this work is Irwin's stress intensity factor (SIF) [11] criterion, which states that for a given crack deformation, i.e. a particular crack opening mode, the elastic fields in the vicinity of the crack field are a dominated by a universal function in a polar coordinate system attached to the crack tip, see Fig. 2, in which the only unknown is a single scalar parameter.This parameter determines the 'intensity' of the loading in the process zone and is referred to as the SIF K.For pure antiplane shear loading the displacement field in the vicinity of the crack tip is approximated in terms of K by Equation ( 2) can be used to determine K for a given, i.e. simulated, displacement field.For example, in a problem with a characteristic length scale L the displacement jump across a crack at a small distance r L away from the crack is given by which can be solved for K as The current stress intensity can subsequently be compared to a material specific critical SIF K C at which crack growth occurs in experiments in order to obtain a crack growth criterion given by i.e. a crack will grow if the SIF reaches the critical value.3 Lattice Boltzmann Method for the Wave Equation The LB scheme for wave equations, as proposed by Yan [30], was already presented in [22] and summarized in [24] as well.Therefore only a short overview is given here, which aims at readers that are already familiar with the LBM in general.
The body B, is discretized by a lattice representation B with uniform spacing ∆h, see Fig. 3.A D2Q5 lattice scheme is used, meaning that five distribution functions f α per lattice site are regarded, one for each lattice velocity c α , α ∈ {0, 1, 2, 3, 4}, defined as where c = ∆h /∆t is the lattice speed with a time step ∆t and c is closely related to the shear wave speed c s .
The lattice Boltzmann equation with the BGKW-collision operator [2,29] is given by For the wave equation, the distribution functions are connected to the particle velocity and the equilibrium distributions for the wave equation are defined as where and τ = ∆t.
Finally, the displacement w is computed from the particle velocity ẇ by an Euler integration scheme via

Non-Lattice Conforming Boundary Conditions
Boundary conditions, which are able to accommodate arbitrary boundary geometries, are important for the accurate representation of a crack within the body B. For the LBM used in this work, two strategies have been proposed in [24].A summary of the macroscopic strategy is given here, since boundary handling is relevant to the algorithm proposed in Sec.4.1.
The LB equation ( 7) is able to determine the evolution of the distribution functions, and thus also the evolution of the macroscopic displacement, in the interior of the domain.However, lattice points at the boundary miss one or more lattice links, i.e. neighbor lattice points.These lattice points are denoted as boundary lattice points x B .For boundary lattice points, not all distribution functions can be obtained from collision and streaming, i.e. by the lattice Boltzmann equation (7).Instead, the missing distribution functions have to be determined by the boundary conditions.
In this work, we employ a macroscopic algorithm that is capable of handling non-lattice conforming boundary conditions.The algorithm approximates the displacement field w in the vicinity of a particular boundary lattice point by means of a quadratic polynomial that is determined by the value of the boundary condition as well as the displacement field at two interior points x I and x II , see Fig. 4. The value of the displacement field at the interior points is determined by bilinear interpolation inside the cells C I and C II .Evaluating the resulting polynomial at x B yields a linear equation in terms of the unknown displacements at boundary lattice points for each x B .These equations are assembled in a linear system of equations of the form where S only contains information depending on discretization and geometry, and R also involves the current value of the boundary conditions.Note that S is time dependent if geometry and lattice change over time as is the case during crack propagation.
The system of equation is subsequently solved for the unknown displacements at the boundary lattice points w B , which is then used to determine the value of the missing distribution functions according to where F x B is the set of distribution functions that can be determined by the LB equation (7) and n miss is the number of missing distribution functions at x B that cannot.
The part of the algorithm that deals with the implementation of the boundary conditions determines the macroscopic field w(x B , t + ∆t) at each boundary lattice point such that it is consistent with the boundary conditions on the macroscopic scale.This is why, we refer to the algorithm as a 'macroscopic' algorithm for the treatment of the boundary conditions.

Implementation
This section describes the algorithm and considerations for the implementation in general terms.Further details can be found in Sec. 5, where numerical models and results are discussed.

Concepts and Algorithm of Crack Propagation
The propagation of dynamic cracks is handled mostly in a rather simple geometric manner.An initial crack is needed.
It is modeled as a line and independent of the lattice.This also defines the initial crack tip, viewed as a point, and the boundary conditions.Additionally the direction ĉ of crack growth is needed.This restricts problems to straight cracks, but in turn the lattice can be easily aligned to single cracks.Our implementation allows two types of simulations.First, it is possible to prescribe crack growth at a constant rate ȧ = ∆a /∆t.The second option is to let crack growth be determined by the K-criterion, based on the SIF.For this case the critical value K C at which a crack will grow must be specified.
A time step 2 of the LBM ends with the time integration (10), thus the displacement field is fully updated at every point of the lattice B. The handling of crack propagation is then appended to the end as an additional step.This step itself is subdivided, as described in Alg. 1.First the criterion is evaluated.This is trivial for steady growth.For the K-criterion the SIF is computed according to Eq. ( 4) and compared to the critical value K C , see also Sec. 2.2.This procedure is repeated for multiple crack tips, if needed.
If the crack grows, the crack tip is moved along the direction ĉ.Subsequently, a new segment is created, see Fig. 5a, between the previous and the new crack tip position, with length where v = ȧ/cs is the relative speed.Since this new segment acts as a boundary within the computational domain, the adjoining lattice points need to be processed accordingly.Generally, these points must be found first.For each point, every associated lattice link is examined, see Fig. 5b.These links can be treated as lines connecting the point and its respective neighbor and are checked for an intersection with the crack segment.Any intersected link is subsequently severed and no information is exchanged along it in further time steps.Both associated points need to be processed for the boundary conditions and their implementation, as for the initial crack.For the macroscopic implementation proposed in [23,20] and used throughout the numerical models in Sec. 5, this entails extending and inverting the boundary coefficient matrix and expanding the vector of interpolation coefficients by the newly found boundary points, see Fig. 5c.This task is computationally expensive and requires more time, the longer the crack grows.
While the growth of the crack and determining severed lattice links is a universally applicable concept, the processing of new boundary points needs to be adapted to different techniques of boundary handling, e.g. when a different LBM for solids is employed.
For this LBM, the lattice wave speed c should surpass the shear wave speed c s , i.e. c = κ c s , where κ 1.Since ȧ < c s , it follows that ∆a < ∆h /κ ∆h.Thus only one pair of new boundary points is expected, at most, for straight cracks.Thus the number of links to be checked can be severely reduced by regarding only the immediate vicinity of the 2 The steps of the LBM, including boundary conditions, are summarized in Alg. 1 of [24].return B 31: end function crack tip.A queue is generated from the neighbors of the last pair of boundary points that have been found.For every point in the queue, the associated links need to be checked.Once a severed link is found, the neighbors of the linked points are added to the queue.The number of checks to be performed can reduced further by marking points that have been completely visited, i.e. all associated links have been checked.This type of processing also works for different kinds of LB methods, e.g. on D2Q9 lattices, and for cases in which the direction of crack propagation is not prescribed.

Evaluation of Stress Intensity Factors
The evaluation of the SIF is carried out according to Eq. 4. The crack opening displacement δ is computed as the difference between w at two lattice points adjacent to, but not on the crack itself, see Fig. 7.For stationary cracks, such as in [22,23], this is straight forward, with a clearly defined distance r from the crack tip.The SIF can be determined in a post-processing step.However, with a growing crack r varies between time steps and for the K-criterion especially.In addition, the SIF has to be evaluated in every time step.The lattice points for evaluation are chosen, such that their distance r to the crack tip lies within an interval [r min , r min + ∆h].As reported in [22,23], the results for the SIF are closer to analytical values when the evaluation occurs at a distance from the crack tip, which is indicated here by the parameter r min .The SIF K is rather sensitive to r, thus for the K-criterion this parameter should be adaptable during   runtime.Preliminary numerical results showed a correlation with v, which can be modeled by This function is shown in Fig. 6a, together with the values of r min used for steady growth in Sec.5.1.It is chosen purely from empirical considerations.

Regularization of the Crack Velocity
With the discretization of time, crack propagation is discretized as well, allowing the crack to grow by a finite length ∆a within ∆t.When using the K-criterion to determine the crack propagation, K is evaluated at the end of each time step and might surpass the critical value K C .This kind of overshooting the critical SIF should be avoided since in real systems, the crack growth would continue to such a state that K K c .
In order to reduce the unphysical overshoot of the SIF, the crack velocity is allowed to increase, up to the maximum crack velocity of v max .Here, the crack velocity is assumed to be a continuous function of K, see Fig. 6b.

Numerical Results
This section shows numerical results to validate the algorithm and its implementation presented in the previous section.
The first examples verifies the evaluation of K in a dynamical model.Since this has not been done before with a propagating crack, a problem with an analytical solution is used for steady crack growth.The second example delivers proof of concept for the K-criterion.The numerical results are assessed with regard to plausibility.

Stress Intensity Factor
We consider a semi-infinite crack that grows at a constant rate ȧ < c s , i.e. v < 1, in a semi-infinite elastic body, as shown in Fig. 8 a).We consider the surface y = h to be subject to a Dirichlet boundary condition and the (moving) crack faces to be traction-free, i.e.
In [15] it is shown that the steady-state stress intensity factor for the problem described above is where The original problem depicted in Fig. 8 a) implies that w = 0 at the lower crack face.In this work, a slightly different but related problem is simulated, see Fig. 8 b).We apply at the top and bottom edges of a strip with width 2L in which a traction-free crack propagates at a steady velocity in x-direction.Although only half the displacement, i.e. 1 2 w 0 , is applied at the top and bottom edge compared to the Table 1: Numerical and statistical data regarding steady crack growth for different values of v; K theo is the expected value from (18).The mean with standard deviation σ and median with differences to the 25th and 75th percentile (cf.Fig. 9) are gathered from experimental values for K, evaluated at rmin /∆h.v rmin /∆h K theo mean ±σ median −25% +75% 0.2 1.50 0.1627 0.1664 0.0083 0.1640 0.0036 0.0064 0.4 2.25 0.1609 0.1642 0.0057 0.1624 0.0035 0.0068 0.6 4.00 0.1569 0.1576 0.0029 0.1584 0.0033 0.0018 0.8 8.00 0.1477 0.1475 0.0007 0.1477 0.0008 0.0005

Validation
The analytical solution for the SIF by Mandal [15] is based on the assumption of a quasi-stationary problem.This implies certain constraints on the geometry of the domain and the prescribed boundary conditions.As mentioned before, we approximate the infinite strip Fig. 8 b) by the domain Fig. 8 c).The left and right edges, as well as the crack faces, are free boundaries.On the upper and lower edge the displacement is prescribed by Since sudden changes in the load can excite spurious waves, the displacement is gradually increased to its final value of 1 2 w 0 at the time t 0 , where t 0 is chosen in relation to the rate of crack growth.The crack continues to grow in a quasi-stationary period until a time of t f > t 0 + 15 cs /L.
The initial crack is rather short compared to r min , with a length of 0.5L.But it grows during the start up period and can be considered to be semi-infinite between t 0 and t f .The position b of the right edge is chosen, such that it does not influence the fields around the crack tip at t f , i.e. b(v) > vc s t f + 2L.Thus both the start up time t 0 and the total length of the strip have to increase with the relative speed v in order to yield comparable results to the infinite strip.This numerical experiment has been conducted for different crack speeds, i.e. v = 0.2, 0.4, 0.6 and 0.8, with a lattice spacing of ∆h = L /16.As described in Sec.4.2, r min is adjusted for each value of v.But the actual value of r differs between time steps within a certain interval, cf.Fig. 7, because the crack tip position changes relative to the lattice points that are used to evaluate K. Due to the sensitivity of K regarding r, this causes K to fluctuate around the expected steady value.Thus a statistical evaluation is undertaken with a total of 300 data points, that are sampled for t ∈ [t f − 15 cs /L, t f ].From this data, the arithmetic mean value with the standard deviation and the median value with the 25th-and 75th-percentile are gathered.This evaluation is compiled in Tab. 1 and the median is also depicted in Fig. 9.The results are close to the analytical curve, each within the margin of error indicated by the percentiles.This margin gets smaller for higher v.The values of K lie above the analytical values, with the exception of v = 0.8, where it is easier to adjust r min .

Validation of the K-criterion
The domain for this example is a rectangle of size 3L by 8L with an initial crack of length 1L along the x-axis, within x = ±0.5L.The crack is located at 1L from the upper and 2L from the lower edge, as depicted in Fig. 10.The lattice spacing is given by ∆h = 2 −6 L, resulting in a total of 98 304 lattice points.
Both crack tips can propagate horizontally along the x-axis.For this example, K C = 0.0055 µ

√
L is chosen, with r 0 = 0.03L ≈ 2 ∆h for r min in Eq. ( 14) and v max = 0.85 as the maximum velocity.For all boundaries the macroscopic implementation of non-lattice conforming boundary conditions, as described in [24], is used.The lower edge is subjected to homogeneous Dirichlet boundary conditions w * = 0.Both lateral edges and the crack faces are tractionfree with t * z = 0, while the upper edge has a time-dependent Dirichlet boundary condition prescribed by which is a half-period of a sine-function.This excites an elastic wave, which then propagates through the domain and is reflected at the outer edges and the crack.
Upon reaching the crack, the incident waves cause the SIF to rise, as can be seen in Fig. 11.When K C is surpassed at a crack tip, the increment ∆a = v ∆t is computed by means of the function v(K; K C ) as defined in Eq. ( 15).Fig. 11 shows the SIF before the crack grows, thus K can be still be higher than K C , in spite of the modification (15).In the next iteration, K should be close to K C .However, since the configuration changes dynamically, K can potentially surpass the critical value again.
The initial wave leads to an increase in K, such that both crack tips initiate crack growth.The subsequent reflected wave also leads to a propagation of the crack tips.This time, K exceeds K C slightly more, due to the changed geometry of the domain, and this results in a higher velocity.During a period of crack growth, K stays at a value close to K C .
After the wave has passed, crack growth halts since K decreases again.

Remarks on Efficiency
To get an impression of the impact the additional crack propagation step has on the efficiency of the LBM, the added computational cost has been examined 3 .For this, the steady growth example with v = 0.4 has been repeated.2 000 time steps were computed with for a total of 10 496 lattice points.The CPU time spent solely on the crack propagation, but also on the computation in total, was measured.By design of the example, crack growth occurs in every time step, but only accounts for about 2.8% of the CPU time.This work introduced a method to simulate dynamic crack propagation using an LBM for solids.In contrast to the more established alternatives, such as finite element methods, it is based on the rather efficient concept of processing very few lattice points in each time step.This keeps the regular lattice unchanged, without the need to re-mesh the domain.While no intensive study on the efficiency has been undertaken so far, it can be projected that the overall efficiency of the LBM is maintained.Only a small number of lattice points is processed for the identification of new boundary points.This processing step is only needed if the crack propagates.Furthermore, boundary conditions are not initialized in each time step that crack growth occurs.Thus little computational effort is added to the LBM due to this crack propagation step.This is corroborated by the comparison of CPU times for the examples in Sec.5.3.Additionally, the matrix inversion necessary for the boundary conditions is a computationally expensive operation, which accounts for some of the needed CPU time.With different boundary conditions, the additional cost could be reduced.
Two cases have been implemented and validated in numerical experiments.In the first one, steady growth of a mode III crack has been compared to analytical results.This shows a good accuracy for the evaluation of the SIF in a dynamical model, as described in Sec.4.2.In the other case, crack propagation with a criterion based on the SIF has been simulated.

Figure 1 :
Figure 1: Antiplane shear deformation of a linear elastic solid with outer normal vector n subjected to Neumann boundary conditions t * and Dirichlet boundary conditions w * , cf.Schlüter et al. (2018)[22] .The outline of the body in the deformed configuration is indicated by dashed lines.

Figure 2 :
Figure 2: Local crack tip coordinate system.The coordinate ξ represents orientation of the crack tip.

Figure 3 :
Figure 3: (a) Lattice representation B of the elastic solid and (b) and the associated lattice velocity vectors (lattice links) for a single lattice point, cf.Schlüter et al. (2018) [22].

Figure 4 :
Figure 4: A section of the considered body where non-lattice conforming boundary conditions are applied at boundary lattice point x B .
(a) The crack is extended by a segment (stippled line) of length ∆a = ȧ ∆t, with P1 as the new crack tip.(b) For finding the intersected link (blue line), all links connected to points in gray boxes are successively checked.(c) The new boundary points are processed.The interpolation involves the points on the blue line, including the closest point on the crack.

Figure 5 :
Figure 5: Handling of a dynamic crack (dark red) with processing of boundary points (filled circles).
(a) The minimum distance of evaluation rmin (14) for r0 = 0.07L (blue line).The green dots represent values chosen for steady growth (see Fig. 9).(b) Regularization of the crack velocity v (15) with respect to the relative overshoot of the SIF K /K C .

Figure 6 :
Figure 6: Functions introduced to define continuous parameters for the K-criterion.

Figure 7 :
Figure 7: Points (blue circles) along the crack (red line) at a distance r ∈ [r min , r min + ∆h] (gray area) are used for the for evaluation of the SIF.

Figure 8 :
Figure 8: (a) Geometry of Mandal's original problem [15] and (b) domain considered for the numerical examination.(c) Numerical realization of this geometry.

Figure 9 :
Figure 9: SIFs for different crack velocities v in comparison to the analytical values of Eq. (18).The values represent the median with errorbars indicating the range from the 25th to the 75th percentile.

Figure 10 :
Figure 10: Domain for the validation of the K-criterion

Figure 11 :
Figure 11: Stress intensity factors, velocity of crack propagation and total additional length for the example of crack growth with the K-criterion, with a maximum allowed velocity v max = 0.85.Due to symmetry, only data for the crack tip propagating in positive x-direction is shown Algorithm 1 Handling of dynamic cracks in the LBMRequire: criterion, direction d, ȧ for steady growth, K C for K-criterion for all n p ∈ N (p) do B ← B ∪ {p, n p } 26: Q ← Q ∪ N (n p )\V