On dynamic crack propagation in a lattice Boltzmann method for elastodynamics in 2D

In recent years, the development of lattice Boltzmann methods (LBMs) for solids has gained attention. Fracture mechanics as a viable application for these methods has been presented before, albeit for mode III cracks in the context of an LBM for antiplane shear deformation. The performance of the LBM itself is promising, while the usage of a regular lattice simplifies the modelling of fractures significantly. Recent advancements in LBMs for solids, especially the description of Dirichlet‐ and Neumann‐type boundary conditions, now make it possible to extend the LBM simulation of crack propagation to the plane strain case with modes I and II crack opening, including growth with non‐uniform speed in arbitrary directions. For this, the configurational force acting on a crack tip is utilised. The definition of the moments of the LBM, which are based on the balance laws of continuum mechanics, render the evaluation of macroscopic fields in the configuration straightforward. In this work, the general in‐plane case of dynamic crack propagation is shown and necessary considerations for the implementation are discussed. Lastly, numerical examples showcase the capabilities of the proposed method to model dynamic fractures and establish a proof‐of‐concept.


INTRODUCTION
Fractures are regarded in a multitude of different fields of science and engineering.The most general case, propagating cracks under dynamic loading, is difficult to describe analytically and also present a demanding task for simulations.
The change in topology during the runtime of simulations typically requires either a remeshing of the computational domain in a standard finite element method (FEM) or more specialised discretization methods, such as the extended FEM (X-FEM) or scaled boundary FEM (SB-FEM).Free fracture models [3] can solve the problem of crack growth, such as the phase field method, due to the regularisation of the crack interface, or the increasingly popular peridynamics.However, for these methods, the numerical performance can be an issue.
With regards to the computational performance, the lattice Boltzmann method (LBM) is rather promising, as it is easily parallelised and scales well [4].While the method was initially conceived for fluid dynamics, it is also actively developed to be applied in solid mechanics.Recent advancements describe a numerical method [5] for plane strain deformation of elastic solids, which is inherently dynamic.Furthermore, newly developed implementations of boundary conditions [1] now allow the extension of the initial modelling of crack propagation to the more general case of mixed-mode loading.
In addition to its computational performance, the LBM also offers the simple geometrical modelling of its underlying regular lattice.This is also favourable for the simulation of cracks.An LBM for solids has already been applied to crack propagation [6,7], where the proof-of-concept was given for mode III cracks, albeit in rather reduced mechanical and numerical setting.The same algorithm is used now with modifications due to the different boundary conditions.

NUMERICAL METHOD
An in-depth description of the principles underlying the numerical method would greatly surmount the scope of this work.The reader is referred to textbooks [8,9] and previous publications [1] and references therein.This section restricts itself to those definitions, that are necessary for understanding the crack propagation algorithm presented in Section 3.
The LBM as a method resorts to approaches from statistical mechanics -solving transport problems via Boltzmann's equation [10] -by regarding the evolution of probability distributions in phase space.This is modelled as a relaxation towards a local equilibrium state, through collisions of (fictitious) particles, with a subsequent streaming of distributions, facilitating the exchange of information about the configuration throughout the lattice.Thus the LBM works on a mesoscopic scale.The connection to the macroscopic scale, that is, continuum mechanics, is given by a decomposition into moments of a probability distribution [11], where the moments are identified with mechanical quantities.
The discretisation of the continuous domain is undertaken with a regular lattice (lattice spacing Δ).Additionally, a set of discrete velocities is needed, which is defined by the lattice scheme, see Figure 1, in relation to the lattice velocity .In the method presented in his work, a D2Q9 lattice scheme is used with , where Δ is the time step and velocities are numbered with index  ∈ {0, … , 8}.This connects each lattice point to its eight surrounding neighbours.Note that the zeroth component resides at the lattice point.

LBM for elastodynamics
The aim is to describe the dynamic behaviour of linear-elastic solids, which can be modelled through the Navier-Lamé equation for the displacement field   , where  and  are the Lamé parameters, characterising the material.The method is based on Murthy et al. [12], using crystallographic lattices, and Escande et al. [5], who developed this further for the more common D2Q9 lattice.Subsequently, boundary conditions have been extended [1], with some additional modifications to the method.
The starting point is a chain of balance equations from continuum mechanics for the mass density , the linear momentum density , and a constitutive law for a stress measure , Spatial derivatives d  are taken with respect to the reference configuration.Equations ( 2)-( 4) are then modelled as a chain of moments of probability distributions , where describes a source term, where the gradient is obtained from second-order finite differences and is called the Poisson stress tensor.Both  and  contain contributions that vanish for a Poisson material, that is, materials with  ≡ , giving Poisson's ratio as  = 0.25.In that case,   only consists of external body forces  and  reduces to the (negative) Cauchy stress .
The evolution of the probability distributions  is given by the Lattice-Boltzmann equation which includes the source term   with weights   per velocity    .The equilibrium state is determined by distributions which are defined in terms of the moments ,  and .
As a post-processing step, the deformation  is computed from the momentum  =  0 u via time integration

Boundary conditions
The general problem of boundary handling in LBMs must be considered in relation to missing distribution functions.
There cannot be a streaming of information across the boundary of the domain.Thus, for certain lattice points   at where  is an exchange term, compare Faust et al. [1] for details, which can now be used to emulate Dirichlet-or Neumanntype boundary conditions, given by the prescribed values of the moments  * and  * .This is the connection of macroscopic quantities at the boundary and the mesoscopic distribution functions at the lattice points.For the Neumann boundary conditions, the prescribed traction vector is translated to the Poisson stress tensor  * .For Dirichlet boundary conditions, the momentum  * needs to be prescribed, which follows from the time derivative of a given displacement  * .

CRACK PROPAGATION ALGORITHM
The approach to handle the modelling of crack propagation was published in Müller et al. [6], see also Müller et al. [7].This algorithm is used here without significant changes.In fact, the implementation is more straightforward due to the linkwise boundary conditions, compare Section 2.2, in comparison to the global interpolation matrix needed for the antiplane shear LBM [13].

Outline of the algorithm
The problem of crack growth is handled as the last part within one time step, after the configuration, that is, the deformation field, of the solid is updated.The general approach makes use of the lattice structure, but is otherwise independent of the discretisation.The crack itself is modelled as a polygonal chain by tracking the position of the crack tip.The new segment for one time step is obtained by shifting the crack tip along a direction with a fixed speed, thus defining the increment.The fracture problem is coupled back to the elastodynamics via boundary handling, since a crack acts as a boundary within the domain.The algorithm essentially performs three steps: 1.The fracture criterion is evaluated.This step needs the information about the elastic fields, which couples the fracture problem to the elastodynamical problem.Assuming a growing crack, the increment is computed from the speed of growth and the direction.These can be given by user-input or derived from the criterion.2. The lattice links intersecting the new segment need to be severed, thus inhibiting the exchange of information between lattice points across the crack.In Müller et al. [6], an algorithm is described, that checks the links between lattice points in the vicinity of the crack tip.Only a very small number of points need to be checked.3.Those points with intersected links need to be treated as boundary lattice points.In the present case, only very simple geometric information in relation to the boundary is needed.
In the next time step, the algorithm for the LBM will now apply boundary conditions accordingly, thus coupling the crack propagation algorithm back to the LBM.

Fracture criterion
For the evaluation of a fracture criterion, an extension of the well-known -integral [14] is used.Here,   is a vectorquantity representing the configurational force acting on a crack tip [2].The following formalism is used, starting with the total (Hamiltonian) energy density  and the Eshelby stress tensor Σ  , defined in terms ,   and   , derived from the LBM, ℊ  = Σ , (material force density), ( 17) In comparison to the classical (static) -integral,   has an additional contribution from the time derivative of the canonical momentum, thus taking dynamical effects into account.
As a fracture criterion, the maximum energy release rate (MERR) [15] is used, where the energy release rate  is given by the absolute value | |.Propagation happens, when the criterion | | ≡  >   , for some critical value   , is fulfilled.Moreover, the preferred direction of propagation is given by − .
The growth rate ȧ can be evaluated from the energy release rate as well.Starting with the general relation [16] between the dynamical and stationary value, where   is the Rayleigh velocity.Assuming, that  stat ≡   , the growth rate can be approximated via the linearized overshoot of  over   , by Additionally, a regularised growth rate can be used, which is ensured to be bounded by the Rayleigh velocity, with an empirical scaling factor .

Remarks on the implementation
Equations ( 15)-( 19) are evaluated on each point of a domain  around the crack tip.The integral in Equation ( 19) is assumed to be path-independent.However, in contrast to the more common path-integral over  for the -integral, a volume integral over  is employed, which is easily implemented on a regular lattice. is defined as a rectangle, moving with the crack tip, wherein each point is assigned a volume of (Δ) 2 .This also agrees with the staircase approximation of the boundary conditions.All gradients and divergences are computed from finite differences, which, again, is readily realised on a regular lattice.Lastly, no rules for self-contact of the crack faces are implemented.Thus, the closing of cracks cannot be handled yet.

EXAMPLES
The numerical examples here are meant to exemplify the capabilities of the crack propagation algorithm.However, they do not represent a validation of the method.Two rather simple test cases are considered: quasi-stationary mode II opening and a mode I crack under dynamic loading.In both cases (non-dimensional), material parameters are chosen to achieve a Poisson ratio of  = 0.223, thus slightly deviating from a Poisson material.
For both cases, a square of size 2 by 2 is chosen for the domain, where  is an arbitrary length scale.The domain  around the crack tip is chosen to be 0.125 by 0.125.The initial crack runs from the left edge in -direction, see Figure 3.The crack faces are traction free via the assigned boundary conditions.
Since the aim here is only a proof-of-concept, any parameters defining the numerical model, such as the critical value, are chosen arbitrarily.

Angle of deflection under mode II crack opening
Pure mode II crack opening is a typical test case in fracture mechanics [17].The initial crack has a length of 1, compare Figure 3A.The lattice is rather coarse with a spacing of 2 −8 , giving 512 × 512 = 262 144 lattice points in total and the domain  is 32 by 32 lattice points.The specimen is held fixed at the bottom edge.The left edge, above the crack, is pushed inwards via a Dirichlet boundary condition.Additionally, a slight pull upwards is applied, in order to counter out any penetration of the crack faces.The crack faces, thus, slide relative to each other, resulting in mode II crack opening.This push is applied quasi-statically, thus suppressing dynamical effects, that is, elastic waves.Eventually, the crack starts to propagate.It grows at an angle, which is measured against the direction of the initial crack, that is, the -axis in this example.This angle is given in literature [17] as approximately −70 • , depending on the fracture criterion.
For this example, the regularised growth rate, Equation ( 22), is used with  = 1, since dynamical effects let the strain around the crack tip rise quickly, once the first lattice connection is severed.However, this cannot be handled correctly by the algorithm at this point.Therefore, the growth rate is limited and only the initial crack propagation is considered here

Dynamic loading of mode I crack
For this example, the specimen is held fixed at the right edge.The crack has an initial length of 0.75.The lattice spacing is 2 −7 , giving 256 × 256 = 65 536 lattice points in total and a domain  of 16 by 16 lattice points.The loading is applied as shearing forces on the left edge, pulling the crack apart, compare Figure 3B.The traction is modulated in time with a sine function, where   lasts for one period, whereas   lasts for two periods, thus introducing an asymmetry in the loading.
Here, the unregularised growth rate according to Equation ( 21) is used.During a first phase, under symmetric loading, a lenticular opening travels towards the crack tip and pulls the crack apart.This eventually results in crack growth in the -direction, which halts at some point, compare Figure 5.During a second phase, the loading is asymmetric, since loading is only applied on the lower left edge.Crack growth starts again, but now following a curved path.

CONCLUSIONS
The numerical examples show plausible results.The angle of deflection of the -initially -pure mode II loading scenario and the subsequent crack path fits the expectations.However, these results depend on the chosen parameters for the loading and fracture resistance.Eventually, it can be concluded, that dynamic crack propagation is indeed possible with a method based on an LBM.While this was already undertaken in a previous work [6], the approach presented there has now been extended to inplane crack opening modes I and II.This extension also includes arbitrary directions of growth and thus allows for curved cracks.This is based on a rather general framework utilising configurational forces.
However, the work on crack propagation is far from finished.Due to the dynamics of the deformation, the growth becomes unstable and an effective control of the growth rate is not possible.This can quickly results in unphysical behaviour.The results of these example were specifically chosen to avoid this.However, additional effects, such as crack branching or the contact of the crack faces, should be included to inhibit some of the unphysical behaviour.
Furthermore, thorough validations of the principles and results are needed, including comparisons with more established methods for fracture mechanics.Studies on the influence of different parameters, notably the size of the domain  in conjunction with the lattice spacing, need to be undertaken as well, since the curvature of the crack could impact the applicability of the configurational force.Also, the lack of boundary rules, that are able to handle the contact of crack faces limits the viability of the method and can lead to unphysical results, especially in the case of dynamic simulations.Furthermore, this method could be extended to incorporate such phenomena as crack nucleation, branching or the interaction of multiple cracks.

A C K N O W L E D G M E N T S
The authors gratefully acknowledge the funding by the German Research Foundation (DFG) within the project 423809639.
Open access funding enabled and organized by Projekt DEAL.

C O D E AVA I L A B I L I T Y
The source code [18], input files for the examples and additional Python modules for the post-processing of simulation data [19] are available under an open license.

F I G U R E 2
The principle of the bounce-back boundary condition, regarding the blue boundary lattice point close to the boundary.The solid is to the left of the dotted line, which marks the (geometrical) boundary.The boundary of the computational domain is approximated by straight lines, halfway between lattice points.Green points are neighbours of the blue point within the domain, while grey points lie outside the domain.The outgoing distribution functions (blue arrows) are send back along the opposite direction (grey arrows).theboundary of the lattice without a full set of neighbours, the values of some incoming distributions need to be determined from the prescribed boundary conditions.Here, (anti-) bounce back rules are used, with modifications for moving boundaries, compare Figure2for details.Outgoing distributions   are send back along opposite direction ᾱ:  ᾱ(  ,  + Δ) = (−)  (  , ) + (,  * ,  * ),

F I G U R E 3
The geometries for (A) the angle of deflection and (B) the dynamic loading.(C) Shows the prescribed boundary values over time.

F I G U R E 4
Result of the example for the angle of deflection: The path of the mode II crack with a line indicating the expected angle of −70 • .The inset shows a magnification of the path.F I G U R E 5 Results of the example under dynamic loading: (A) the path of the mode I crack, (B) the energy release rate , the growth rate ȧ and the total length  of the crack for the mode I example.for the evaluation of the deflection angle.After the onset of crack propagation -and under increased loading -additional effects play an increasing role.The crack path can deviate considerably from the static expectation.Figures4shows the crack path resulting from this simulation.The crack propagates at an angle of −48 • at first, then kinks towards a steeper angle of approximately −72 • and eventually starts to curve rather straight towards the bottom edge.