Crack dynamics and crack tip shielding in a material containing pores analysed by a phase ﬁ eld method

Many naturally occurring materials, such as wood and bone, have intricate porous micro-structures and high sti ﬀ ness and toughness to density ratios. Here, the in ﬂ uence of pores in a material on crack dynamics in brittle fracture is investigated. A dynamic phase ﬁ eld ﬁ nite element model is used to study the e ﬀ ects of pores with respect to crack path, crack propagation velocity and energy release rate in a strip specimen geometry with circular pores. Four di ﬀ erent ordered pore distributions are considered, as well as randomly distributed pores. The results show that the crack is attracted by the pores; this attraction is stronger when there is more energy available for crack growth. Crack propagation through pores also enables higher crack propagation velocities than are normally seen in strip specimens without pores (i.e. homogeneous material), without a corresponding increase in energy release rate. It is further noticed that as the porosity of an initially solid material increases, the crack tip is increasingly likely to become shielded or arrested, which may be a key to the high relative strength often exhibited by naturally occurring porous materials. We also ﬁ nd that when a pore is of the same size as the characteristic internal length then the pore does not localise damage. Since the characteristic internal length only regularises the damage ﬁ eld and not the strain end kinetic energy distributions, crack dynamics are still a ﬀ ected by small pores.


Introduction
Dynamic crack propagation in heterogeneous materials isin some sensesprofoundly different from that in the special case of homogeneous materials.A striking example is that crack propagation velocities close to or even above the Rayleigh velocity (c R ) have been observed, both experimentally and numerically, in materials with a weak zone embedded in an otherwise homogeneous material [1][2][3].Previous studies suggest that crack dynamics is primarily governed by the amount and distribution of energy, and the rate at which this energy can be released [3,4].Elastic heterogeneities, primarily variations in local Young's moduli, affect the strain energy distribution and influence the crack path by either attracting or repelling a crack.Dynamic heterogeneities, which include variations in both local density and local stiffness, affect the dynamic response and wave propagation velocities of the material, and can influence the kinetic energy distribution in a specimen, which in turn influences the crack dynamics [4].
In a material containing pores, the energy distribution is affected by the microstructure and these micro-structural effects can inhibit crack growth and limit the crack propagation velocity, or arrest the crack entirely [5].This is often the case with biological materials such as wood or bone, where evolutionary processes have favoured a combination of toughness and lightness, giving rise to materials with intricate micro-structures which are extremely well adapted to the load cases they are subjected to; e.g.wood shows a ratio of density to stiffness, strength and fracture toughness similar to or even higher than most engineering composites [6,7].
In a qualitative experimental evaluation of crack propagation, a difficulty is that the crack propagation velocity for most brittle materials is too fast for the naked eye to capture (and alsodepending on material and geometryfor moderate-priced high-speed cameras to capture with good resolution).Quantitative experimental evaluation of dynamic crack propagation is often complicated, as the energy release rate varies with the velocity of the moving crack.The relation between crack tip velocity and energy release rate is typically given by the dimensionless universal function of crack speed, ĝ v ( ), cf.[8], where G is the instantaneous energy release rate, G c the critical energy release rate of the material, v is the crack tip velocity and c R , as mentioned, is the Rayleigh velocity of the material.This relation has proved useful for prediction of crack propagation velocity in strip specimens of homogeneous material in several studies, cf.[9,3,10], but has shown to be less accurate when the crack interacts with waves, typically reflected at the boundary of the specimen [11,10].A benefit of the strip specimen is that the strip is in quasistatic equilibrium before the onset of crack propagation, thus the only wave phenomena are due to the propagation.
The study of dynamic crack propagation has primarily been focused on comparatively homogeneous materials such as steel, glass and amorphous polymers like polymetylmetakrylat (PMMA) [12][13][14][15].Interesting exceptions exist primarily within the area of geosciences, where dynamic crack propagation has been studied in rock and rock-like materials [16][17][18][19].Quasi-static crack propagation in heterogeneous (though not necessarily porous) materials has been studied extensively in e.g.civil engineering where dynamic crack propagation in concrete has attracted a lot of attention [20][21][22][23].The authors have not been able to find any appropriate references concerning rapidly growing cracks and crack dynamics in porous or cellular materials; a few studies concerning fracture toughness were found, e.g.[24][25][26].
Here, the influence of pores in a material on crack dynamics in a brittle material is studied by means of a dynamic finite element (FE) phase field model.Phase field models for fracture have received increasing attention over the last decade.In this work we are using a phase field method based on the variational approach to fracture (other methods exist, e.g.methods based on the Ginzburg-Landau phase transition are also common [27]).The variational approach to fracture was suggested by Francfort and Marigo [28] (see [29] for a thorough review), and is based on the original works of Griffith [30].The variational approach is thus closely related to Griffith's ideas, that crack propagation follows the principle of minimum potential energy.According to the works of Griffith, crack growth is irreversible, the fracture toughness is the upper limit of the energy release rate and a crack will not grow unless the energy release rate is critical.While the first observation requires some additional treatment (as will be outlined in the next section), the latter two follow naturally from the variational approach.As remarked by Francfort and Marigo [28], the minimisation principle is strongly reminiscent of the Mumford-Shah functional [31] (it has also been proven that a crack tip is indeed a stationary point of the Mumford-Shah functional [32]).This was used by Bourdin et al. [33], who used the Ambrosio-Tortorelli approximation of the Mumford-Shah functional using elliptic functionals for efficient numerical implementation [34].This approximation introduces a new variable, which can be thought of as a regularised crack phase field.The crack phase field depends in turn on a regularisation parameter of dimension length, often referred to as an internal characteristic length.Since the minimisation covers all possible crack states in the body there is no need for additional criteria for crack initiation, path, bifurcation etc.This makes the phase field method an ideal candidate for studies on porous and heterogeneous materials since a crack is able to "jump" across a heterogeneity or pore if energetically favourable.Also, since the introduced internal characteristic length can be related to the ultimate tensile strength and amount of work done on a specimen before fracture (cf.[35]), an energetically consistent phase field model prescribes a value of the internal length.
The phase field method can be straightforwardly extended to dynamics, which has been shown by, primarily, Bourdin, Larsen and Richardson [36,37] who proved the existence and convergence of solutions to dynamic phase field problems for fracture, as well as [4] and Kuhn and Müller [40].While our approach to evaluate energy release rate directly from dissipated energy is similar to that of the former two, the latter uses an equivalent concept of the near-field -integral, , tip which is evaluated in a contour around the crack tip [40].

Phase field model
In this work, dynamic crack propagation is studied using an explicit time integration finite element scheme in which the displacement and crack phase field are solved in a staggered manner.Considering a two-dimensional problem domain Ω, with exterior ∂Ω (the derivation in 3D is analogous), part of which (∂Ω T ) is subjected to natural boundary conditions (prescribed stress T ) and part of which (∂Ω u ) is subjected to Dirichlet boundary conditions (prescribed displacement U ), and a discrete crack whose surfaces are denoted Γ (Fig. 1a), the total free energy of the system can be written as Here ψ k is the kinetic energy density, ψ e the elastic energy density, and x d denotes integration over spatial coordinates.The line integral over Γ is evaluated over the line segments ds, the integral thus represents the energy consumed by surface creation [29].Through the phase field implementation, the discrete crack Γ is represented as a diffuse regularised crack field d defined over the whole domain Ω, here chosen such that = d 1 represents broken material and = d 0 represents intact material, (Fig. 1b).We have made use of a crack density functional which is linear in d and quadratic in ∇d (cf.[35]), where ∇ is the differential operator and l the regularisation parameter (i.e.internal characteristic length), determining the width of the regularised crack.The crack density functional, while typically chosen as either first-order in d (as in Eq. ( 3)) or second-order in d (which was used in e.g.[9]), can be chosen with great freedom from the family of elliptic functionals approximating the Mumford-Shah functional (cf.e.g.[41]).The first suggestion of the functional used in this work is (to our understanding) made by Francfort and Pham [42,43].Each functional gives a somewhat different stress-strain response, see e.g. the many works of Pham, Marigo and coworkers, e.g.[35,44].The one in Eq. ( 3) has the benefit of having a distinct elastic limit.With Eq. ( 3), the last term of (2) can be approximated by The kinetic energy is where u, and later occurring u ¨, denotes first and second derivatives of displacement with respect to time t and ρ is the density.The kinetic energy is unaffected by the crack phase field d, but d locally degrades the elastic stiffness of the material.In order to account for crack surface contacts, the elastic energy density is split into a positive and a negative part, such that only the tensile-originated strain energy is degraded by the crack phase field d while stiffness is kept in compression to account for crack closure, (cf.[45]) ) .Assuming, as is the case for an isotropic material, that the stiffness tensor can be additively divided into two parts 1) in plane stress and μ and λ are the Lamé parameters and ν is Poisson's ratio, the positive and ne- gative strain energy are taken as and where ε is the Voigt (vector) representation of the linearised strain tensor, ̃= T .The parameter α is determined by the trace of the strain tensor, ε tr and takes the value = α 1 if ̃⩾ ε tr 0 and = α 0 otherwise.Following Miehe et al. [45], a spectral strain decomposition is used, ε is divided into ̃+ ε and ̃− ε by where e i are the eigenvalues of ε n , i are the corresponding eigenvectors, ε / are the Voigt notation vectors of the tensors ̃+ − ε / .It should be noted that, for an isotropic material, the decomposition outlined above is the standard split introduced by Miehe [45], given in (2D) matrix notation.
The Cauchy stress tensor is given as the derivative of the strain energy with respect to strain, ) .
e e e 2 (6) Taking a second derivative gives the positive and negative parts of the consistent tangent stiffness tensor, .  6) and ( 7) are readily calculated using the symbolic toolbox in Matlab [46].
With these relations, the total energy of the system becomes To account for irreversibility of the crack evolution a history field is used in place of + ψ e , such that u t ( , ) is the maximum positive strain energy experienced (at time ⩽ ⩽ τ t 0 ) [45], The history field ensures that it is the largest strain energy experienced in the material during the simulation history that determines the present stiffness.Stresses are however reversible, and are always evaluated according to Eq. ( 6).
The functional ), with integration over Ω, now depends on two independent parameters and their respective deri- vatives: u and u, and d and ∇d.Using the principle of least action, we first obtain the Euler-Lagrange equations for u with respect to time.Then, using the principle of minimisation of potential energy, we obtain the Euler-Lagrange equations for d, for each displacement.This variation of Eq. ( 8) gives the governing equations in the strong form, Here we have neglected any body forces.The equations of motion also require boundary and initial conditions (n is the normal to the boundary),

Spatial discretisation
Obtaining the weak forms of Eq. ( 9) is thoroughly described in many publications e.g.[20].Here we only give a brief review since the weak form is a necessary ingredient in the formulation of the FE problem.By multiplying Eq. (9b) with a test function δd, integrating over the domain Ω, making use of the product rule for derivatives, the divergence theorem and the boundary condition of Eq. (10c), we arrive at the weak form of Eq. (9b), Discretising, such that / are some standard finite element basis functions and their first spatial derivative, and observing that the test function is arbitrary, we get the discretised equation for the crack phase field as, where Since a staggered algorithm is used, there is no coupling between the displacement and crack phase fields.
Likewise, for Eq.(9a) we multiply by a test function u δ , integrate over Ω, make use of the product rule for derivatives, the divergence theorem and of the boundary conditions in Eq. ( 10), and we arrive at the weak form of Eq. (9a), Discretising, such that where u i are the nodal values of the displacement field u, and -again -observing that the test function is arbitrary, we get the discretised equation for the crack phase field as, (12) where Here, we have used a lumped mass matrix, obtained by using a quadrature only involving the nodal points (only for M ).

Time integration
In this work, Eq. ( 12) is solved numerically using a fully explicit Newmark time stepping scheme [47].For each time step, the displacement is solved explicitly, then the fracture phase field, Eq. ( 11) is solved, based on the displacements.
In the Newmark scheme, the displacement of the new time step + k 1 is predicted by where t Δ is the length of the time step, and u u , k k and u ¨k are the displacement, velocity and acceleration of the previous time step k.The predicted displacement, + uk 1 , is used to solve the dynamic system of equations for the displacement problem (Eq.( 12)) by where M K , and + f k 1 are the mass-and stiffness matrices and load vector (of the new time step + k 1) of Eqs. ( 13)- (15).The acceleration (Eq.( 17)) is then used to calculate the velocity + uk 1 and updated displacement + u k 1 as Throughout this study, a fully explicit Newmark algorithm is used, with constants = β 1/2 1 and = β 0 2 .

FE model
Finite element analyses are performed on strip specimen geometries with dimensions ( × × = × × L H T 50 25 2 mm) (Fig. 2) [48,15].A state of plane stress is assumed.Material properties are given in Table 1.The length parameter is taken to be of the same order of magnitude as the energetically consistent value, , where σ c is the ultimate tensile strength and/or elastic limit of the material cf.[35].The geometries are discretised using irregular meshes consisting of 40,000-120,000 four-node quadrilateral finite elements (Fig. 2), to ensure an element length of maximum 0.075 mm in the proximity of pores; the regularisation parameter l is three times this value.A crack of length x0 is inserted by specifying a strain history field (cf. [38]).By varying the length of the initial crack, the energy release rate is varied, since a specimen with a short initial crack will require a larger boundary displacement before onset of crack propagation.In the analyses, the strip specimens are first loaded quasi-statically ( = − u Δ 0.5•10 6 m) until fracture nucleates, after which the crack propagation is simulated using an explicit Newmark algorithm.A time step of = − t Δ 2•10 8 s is used; this is about a third of the critical time step (the so-called CFL condition [49]), where h is the smallest element length and c p is the pressure wave (p-wave) speed of the material.Eq. ( 18) ensures that the time step is smaller than the time required for a p-wave to traverse an element.During the dynamic part of the simulation, a vertical velocity of ± 2.5 mm/s is applied to the upper and lower boundary (Fig. 2).Four different geometrical configurations, all with circular pores with diameter = R 2 1.5mm, are studied (Fig. 3).The pores are distributed evenly over the length of the specimen, i.e. = D L/9.The white circles in Fig. 3 represent pores, the grey areas are virgin, unaltered material.Geometry 1 simulates the effect of a crack propagating through a row of holes, where the row of holes is located in a straight line ahead of the crack tip.Geometry 2 is used to investigate the attraction of the holes on the crack tip.This geometry also consists of a row of holes, but now offset by a distance D/2 from the crack tip.Geometry 3 and 4 simulate the effect of a crack propagating between two rows of holes, where the crack is initially between the holes.
In addition, a number of simulations have been run, in which pores have been randomly distributed over the entire sample.In these models the pore radius is 1 mm and the relative density (proportion of intact material to entire volume) is = * ρ 0.86.A small study on the effect of the internal characteristic length parameter, l, in relation to pore radius is also performed.For this purpose, a square model of dimensions × × = × × L H T 25 25 2 mm) is used.The simulation parameters u Δ and t Δ are the same as before.A hole of radius either = R 0.25 mm or = R 0.5 mm is introduced.The parameter l is chosen such that = l 0.25 mm, same as the smaller hole radius.Three different hole placements are studied, a horizontally centered hole located at one quarter of the specimen length, a horizontally centered hole located at half the specimen length and a hole located at one quarter of the specimen length, offset + R l 4 in the vertical direction.
Fig. 2. Geometry and boundary conditions in the numerical simulations.The upper and lower boundary are subjected to a prescribed displacement/ velocity while all other boundaries are free.Insert shows a part of an unstructured mesh.

Comparisons with experiments
For comparison, dynamic crack propagation experiments were conducted on pre-notched strip specimens ( × × = × × L H T 50 25 2 mm) of intact PMMA and PMMA with drilled holes, using a universal tensile testing machine (Shimadzu AGS-X) and a high-speed camera (MotionPro Y8) to track crack propagation, cf.[9].The specimens were loaded quasi-statically using a load cell of 10 kN, and a crosshead speed of 5 mm/min until a sudden fracture occurred, and the crack propagated the length of the specimen within less than 100 μs.
The reference crack pathsboth experimental end numericalwhich were obtained for a material without any holes, are shown   in Fig. 4. Experimental and numerical crack paths for geometries 1 -3 and numerical crack paths for geometry 4, for two different initial crack lengths are shown in Fig. 5.
The overall correlation between experimental and numerical crack paths is good.The general trend is that the longer cracksi.e.lower energy release ratesare less affected by the holes than the shorter cracks, i.e. higher energy release rate (Fig. 5).In general, cracks are drawn towards the holes.For example, both in simulations and in experimental results for a short crack ( ̂= x 0.75 0 mm) in geometry 1, there are several successful branching events, which have been eliminated or attracted back to the midline by the holes.The exception to this rule is the short crack in geometry 2, for which it seems like the crack is repelled by the holes.
A study of the (numerical) strain energy distribution during crack evolution (Fig. 6) for the short crack in geometry 2 shows that the main crack branches before the second hole, after which the upper branch is attracted towards the second hole, follows a curved path and becomes shielded by the second hole, while the lower branch follows a straighter (and faster) path and "overtakes" the upper branch.Also, for geometry 1 (short initial crack) we see that while there is actually more than one crack tip propagating simultaneously, the crack tip strain energy singularity does not fully separate; rather it appears as if there was only one crack tip, as opposed to the behaviour of the reference geometry (cf.Fig. 6).

Crack tip velocity
Fig. 7 shows the normalised instantaneous crack tip velocity v c / R versus normalised crack tip position x L / for 52 μs propagation, plotted over the crack path at = t 52 μs, as obtained in numerical simulations for the different geometries considered.The crack tip position x has been determined by considering nodes with damage value over 0.95 as broken.Since the time step, by the CFL condition (cf.Eq. ( 18)), is smaller than the time required for the crack tip to propagate the length of one element (h), the crack tip position evolves in a stepwise manner.To get an accurate estimate of the crack tip velocity, crack tip velocity is obtained by numerically differentiating crack tip position only between time steps in which the crack tip has advanced a distance of about half an element length.Moreover, since crack propagation velocity through a pore is a fictional concept, velocity is calculated only in the second to last time step before the crack tip stops advancing (i.e.enters the pore) and in the second time step after the crack tip advances again (on the other side of the pore); these points are then joined by a dashed line.The crack has been considered to be stationary if it has not moved in more than 60 time steps.The last and first time steps before and after the pore have been excluded in order to make certain that velocity is only calculated between points in which the crack is propagating.In addition to the grey solid line indicating un-smoothed velocity, a smoothed velocity has been added using a solid black line.
For the reference geometry, the crack tip velocity is more or less constant, as expected for the strip specimen [15].For the simple case of a row of holes in front of the crack tip (geometry 1) we see that the crack accelerates sharply before it reaches each hole.When the crack is re-initiated at the other side of the hole, the velocity is similar to the velocity prior to the acceleration.Also, the average crack propagation velocity is higher compared to the reference (longer distance has been traversed in the same time, 52 μs).
Both for geometry 2 (primarily the long initial crack, i.e. ̂= x 2 0 mm) and geometry 3 (primarily the long initial crack) the crack propagation velocity is higher between the holes, i.e. when the crack tip is not shielded by the holes.Geometry 4, with a short initial ), is the exception to this rule as crack propagation velocity is highest at the end of each hole rather than between holes.This velocity increase coincides with the death of one of the crack tipswhat is cause and what is effect is however not obvious.Interestingly, geometries 2, 3 and 4 all show crack propagation velocities slightly lower than the reference geometry for the longer crack tip, i.e. lower energy release rate, and velocities which are slightly higher than the reference geometry for the shorter crack tip, i.e. higher energy release rate.A possible explanation for this behaviour is that for the lower energy release ratewhere we have seen that the crack is less affected by the attraction exerted by the holesthe crack tip shielding of the holes reduces the velocity of the crack tip (Fig. 7, left).For the higher energy release rate this effect (though still present) is dominated by the attraction of the crack to the holes.

Energy release rate
Fig. 8 shows the normalised instantaneous energy release rate G G / c versus normalised crack tip position x L / for 52 μs propa- gation, plotted over the crack path at = t 52 μs, as obtained in numerical simulations for the different geometries considered.In- stantaneous energy release rate has been evaluated as where E p and E k are the total potential end kinetic energies (since x is crack tip position, Tdx is the change in crack surface area).As was the case with velocity only time steps in which the crack has advanced about half an element length have been considered.If more than 60 time steps have been excluded, the crack is considered not to be propagating, and the jump to the next crack position is indicated in Fig. 8 by a dashed line.In addition to the grey solid line indicating un-smoothed energy release rate, a smoothed energy release rate has been added using a solid black line.
The energy release rates of the geometries with holes are in general lower compared to the reference; by introducing pores we have also reduced the strength of the specimens.In general, the energy release rates of geometry 2, 3 and 4 are slightly lower than the reference.For geometry 1, the energy release rate for the crack propagation between the holes is similar to, or slightly higher than, the reference.
Interestingly, we did not see any reduction in simulated crack tip velocity for the short initial cracks in geometries 2, 3 and 4 compared to the reference, still the simulated energy release rate is lower, cf.Fig. 7. Another interesting feature is that the plots of velocity and energy release rate for geometry 4, short initial crack, are out of phase.When the velocity reaches local maxima at the end of each hole, the energy release rate has local minima, and vice versa at the beginning of each hole.

Relation between energy release rate and crack tip velocity
Based on the observation above we conclude that the approximation in Eq. ( 1) is not fully satisfactory for these geometries.In  1).While neither result coincides exactly with the prediction, the reference follows the principle of Eq. ( 1), that an increase in energy release rate corresponds to an increase in crack tip velocity.On the other hand, those geometries in which the crack has advanced through holes, i.e. geometry 1 (both long and short initial crack) and geometry 3 (short initial crack), deviate significantly from the prediction.The deviation is related to the observation that the energy release rate G appears to drop just before the crack reaches a hole, while the crack tip velocity v increases.1) but for most of the crack propagation they have followed the behaviour of the reference.The out of phase behaviour of geometry 4, short initial crack, is evident from the spiral-like behaviour of the plot for this geometry (dark purple colour).
These deviations from the prediction made by Eq. ( 1) do not imply that Eq. ( 1) is incorrect.Rather it points out that a material containing pores of a certain size does not behave as a continuum at a dominating scale, thus violating the assumption that the crack does not interact with external boundaries, as pointed out in [10].

A general porous material
Random porous micro-structures were studied numerically by randomly distributing holes in the specimens.Again, only time steps in which the crack has advanced about half an element length have been considered.Since the holes are larger than in the previous examples (1 mm), the crack is considered not to be propagating if more than 80 time steps have been excluded, if so the jump to the next crack position is indicated in Fig. 10 by a dashed line.
These specimens showed somewhat earlier onset of crack propagation compared to the reference material which is expected since we have introduced flaws in the material.In general, the crack propagated by linking holes together, and branching was unusual.A general trend was that these simulations were likely to exhibit crack arrest, as opposed to the general behaviour of a continuum strip specimen, cf.[50].This feature makes the strip specimen a poor choice for simulation if dynamic crack propagation in porous  With the pore size considered here the crack is always influenced by wave interference phenomena caused by the interaction of emitted waves from the crack propagation with the holes, at which point the prediction of Eq. ( 1) appears to fall apart.
In a general random porous geometry, there is a large chance that the crack tip becomes shielded, which can slow down or arrest the crack propagation entirely.This effect may be one key to the high relative strength often exhibited by naturally occurring porous materials.
Within this workas with any modelling worka number of modelling assumptions have been made.One of the most debated is the choice of the strain energy split (Eqs.( 4) and ( 5)).In the split we have used, damage evolution is driven by tensile strain-related energy.Implicitly this corresponds to a tensile mode in crack propagation.Without entering into the debate over which isin general the most correct split, we observe that, for a material containing pores, if we let the porosity increase, the tensile mode will come to dominate over the shear mode (see also the phenomenological damage models in [6]).Another parameter under debate is the physical meaning of the parameter l.It can be considered as a material parameterbut if so it is debated whether it should relate to some micro-structural dimension or to the energy consumed in fracture (ideally they coincide cf.[35]).We have found that when a pore is of the same size as l, then the pore does not localise damage.Since l only regularises the damage field, strain end kinetic energy distribution is unaffected by the choice of l (apart from in the damaged band), and crack dynamics are still affected by small pores.This suggests that l could actually be measured in a laboratory by making small holes in the material and increasing the radius until damage localises at the holes.This should occur when the hole is slightly larger than some characteristic length in the microstructure, which supports the first hypothesis.However, the relation between characteristic length l of the damage regularisation and pore radius R, especially in the context of proper porous materials, should be further investigated.
Combined, these results serve to increase the understanding of crack dynamics and dynamic fracture in naturally occurring porous materials, which in turn can help improve e.g. the integrity of structural panels or concrete.Another important issue is that of how to design medical implants, e.g.bone implants or bone cements, for optimal strength and durability.

2 2
in Eqs. ( where d i are the nodal values of the crack phase field d and N and = ∂ ∂ B N x

Fig. 4 . 2 0
Fig. 4. Crack paths of the reference geometry (without any holes), top: experimental results, bottom: numerical results.The left column corresponds to long initial cracks, ̂= x 2 0

Fig. 6 .
Fig. 6.Strain energy distribution for the reference geometry (top), geometry 1 (middle) and geometry 2 (bottom) at (a) 10, (b) 30 and (c) 50 μs.The colour bar refers to all figures.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Normalised crack tip velocity v c / R versus normalised crack tip position x L / for 52 μs propagation, plotted over the crack path at = t 52 μs, for cracks with initial length, left: ̂= x 2 0

Fig. 8 .
Fig. 8. Normalised energy release rate G G / c versus normalised crack tip position x L / for 52 μs propagation, plotted over the crack path at = t 52 μs, for cracks with initial length, left: ̂= x 2 0

Fig. 9 .
Fig. 9. Inverse normalised energy release rate G G / c versus normalised crack tip velocity v c / R for different geometries and initial crack lengths.The dashed line represents the relation predicted by Eq. (1).The right figure is a closeup of the lower left corner of the left figure.

Fig. 10 .
Fig. 10.Normalised crack tip velocity v c / R (left) and normalised energy release rate G G / c (right) versus normalised crack tip position x L / for 52 μs propagation, plotted over the crack path at = t 52 μs, for cracks with initial length ̂= x 0.75 0 mm for three different random hole configurations.Grey lines are un-smoothed results, black line have 5% smoothing.
[3]38,39]er bound of the crack tip velocity to the elastic wave speeds.Contributions with respect to crack propagation velocities have been made in e.g.[2,38,39].Especially interesting in relation to the current work are the works of Bleyer et al.[3], Ylmaz et al.

Table 1
Material properties used for simulation.