Contents lists available atScienceDirect
Engineering Fracture Mechanics
journal homepage:www.elsevier.com/locate/engfracmech
Crack dynamics and crack tip shielding in a material containing pores analysed by a phase field method
Jenny Carlsson
⁎, Per Isaksson
TheÅngström Laboratory, Uppsala University, Box 534, SE-751 21 Uppsala, Sweden
A R T I C L E I N F O
Keywords:
Dynamic fracture Phasefield method Porous material Crack tip shielding
A B S T R A C T
Many naturally occurring materials, such as wood and bone, have intricate porous micro- structures and high stiffness and toughness to density ratios. Here, the influence of pores in a material on crack dynamics in brittle fracture is investigated. A dynamic phasefield finite ele- ment model is used to study the effects of pores with respect to crack path, crack propagation velocity and energy release rate in a strip specimen geometry with circular pores. Four different 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 propa- gation 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 alsofind 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 damagefield and not the strain end kinetic energy distributions, crack dy- namics are still affected by small pores.
1. Introduction
Dynamic crack propagation in heterogeneous materials is– in some senses – profoundly 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 (cR) have been observed, both experimentally and numerically, in materials with a weak zone embedded in an otherwise homogeneous ma- terial[1–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
https://doi.org/10.1016/j.engfracmech.2018.11.013
Received 29 June 2018; Received in revised form 24 September 2018; Accepted 4 November 2018
⁎Corresponding author.
E-mail address:jenny.carlsson@angstrom.uu.se(J. Carlsson).
Available online 27 November 2018
0013-7944/ © 2018 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/BY/4.0/).
T
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 also– depending on material and geometry – for 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,g v( )̂, cf.[8],
̂ ̂
= ≈ −
G
G g v v
( ) 1 c ,
c
R (1)
where G is the instantaneous energy release rate,Gcthe critical energy release rate of the material,v̂is the crack tip velocity andcR, 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 quasi- static 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–15]. Interesting exceptions exist primarily within the area of geo- sciences, where dynamic crack propagation has been studied in rock and rock-like materials[16–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–23]. The authors have not been able tofind 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–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) phasefield 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 thefirst 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 in- troduces a new variable, which can be thought of as a regularised crack phasefield. 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 het- erogeneity 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 phasefield model prescribes a value of the internal length.
The phasefield 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 phasefield problems for fracture, as well as
Nomenclature
̃
ε,ε+ −̃/ linearised strain tensor/its pos. and neg. parts ε ε, + −/ linearised strain vector/its pos. and neg. parts σ cauchy stress vector
Kd, fd stiffness matrix/force vector for damage problem M, K, f mass matrix/stiffness matrix/force vector for dis-
placement problem
N, B FE basis functions/theirfirst spatial derivative u, u̇, ¨u displacement/first/second derivative w.r.t. time
u t
Δ , Δ displacement increment/time step
̂ ̂
x, v crack tip position/ velocity
0, λ, μ undegraded stiffness matrix/related to normal deformation/related to shear deformation
, + −/ consistent tangent stiffness tensor/its pos. and neg.
parts
historyfield
Ω, Γ problem domain/crack surfaces
∂Ω, ∂Ω ,T ∂Ωu exterior of problem domain/subjected to natural/Dirichlet B.C.
T U, prescribed traction/displacement
ψk, ψe, ψe+ −/ , kinetic/elastic energy density/its pos. and neg. parts
cR, cp Rayleigh velocity/p-wave velocity
d l, Crackfield variable/regularisation parameter E, ν, ρ Young’s modulus/Poisson’s ratio/density λ, μ lamé parameters
Ek, Ep total kinetic energy/total elastic energy
G, Gc instantaneous energy release rate/critical energy release rate
̂
g v( ) universal function of crack speed L, H, T length/height/thickness D, ρ∗ hole spacing/relative density
e.g. the upper 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.
[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].
2. Simulation model
2.1. Phasefield model
In this work, dynamic crack propagation is studied using an explicit time integrationfinite element scheme in which the dis- placement and crack phasefield 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 stressT) 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
∫ ∫ ∫
= x− x−
J ψkd ψed G sd .
Ω Ω Γ c (2)
Hereψk is the kinetic energy density,ψe the elastic energy density, anddx 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 phasefield implementation, the discrete crack Γ is represented as a diffuse regularised crack field d defined over the whole domain Ω, here chosen such thatd=1represents broken material andd=0represents intact material, (Fig. 1b). We have made use of a crack density functional which is linear in d and quadratic in∇d(cf.[35]),
∇ = + ∇ ∇
γ d d
l d l d d
( , ) 3
8 ( 2 · ),
(3) 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 eitherfirst-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]). Thefirst 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
∫
ΓG scd ≈∫
38Gl (d+l (∇ ∇d· d))d .xΩ
c 2
The kinetic energy is
∫
Ωψkdx=12∫
Ωρu u ẋ· ̇d ,whereu̇, and later occurring u¨, denotesfirst 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 phasefield d while stiffness is kept in compression to account for crack closure, (cf.[45])
= − ++ −
ψe (1 d ψ)2 e ψe.
Fig. 1. Regularised representation of a crack; (a) sharp crack and (b) diffuse crack.
Assuming, as is the case for an isotropic material, that the stiffness tensor can be additively divided into two parts0=λ+μ, in the case of 2D, in matrix representation (Voigt notation)
=⎡
⎣
⎢
⎢
′ ′
′ ′ ⎤
⎦
⎥
⎥ λ λ λ λ
0 0 0 0 0
λ
and
=⎡
⎣
⎢
⎢
⎤
⎦
⎥
⎥ μ
μ μ
2 0 0
0 2 0
0 0
μ
where ′ = ⎧
⎨⎩ − −
λ λ
λ ν ν
in plane strain
(1 2 )/(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
= +
+ ε ε ε+ ε+
ψ α
2
1
2 T , and
e T
λ μ
(4)
= −
− ε ε+ ε− ε−
ψ 1 α
2
1
2 T ,
e T
λ μ (5)
whereεis the Voigt (vector) representation of the linearised strain tensor,ε̃ =1/2(∇ + ∇u ( u) )T . The parameterαis determined by the trace of the strain tensor,trε̃and takes the valueα=1iftrε̃ ⩾0 andα=0otherwise. Following Miehe et al.[45], a spectral strain decomposition is used,ε̃is divided intoε+̃ andε−̃ by
̃
∑
= = 〈 〉 ⊗+ −
=
ε e + −n n,
i
i i i
/ 1 2
/
whereeiare the eigenvalues ofε n,̃ iare the corresponding eigenvectors, 〈 〉 =x+ x+| |x
2 and 〈 〉 =x− x−| |x
2 . The vectorsε+ −/ 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,
=∂
∂ = − ∂
∂ +∂
∂
+ −
σ ε ε ε
ψ d ψ ψ
(1 ) .
e 2 e e
(6) Taking a second derivative gives the positive and negative parts of the consistent tangent stiffness tensor,
= ∂
∂ = − ∂
∂ +∂
∂ = − +
+ −
+ −
ε ε ε
ψ d ψ ψ
d
(1 ) (1 ) .
e e e
2 2
2 2
2 2
2
2
(7) The derivatives ∂ψe+/ ,∂ε ∂ψe−/ ,∂ε ∂2ψe+/∂ε2and∂2ψe−/∂ε2in Eqs.(6) and (7)are readily calculated using the symbolic toolbox in Matlab [46].
With these relations, the total energy of the system becomes
∫
∇ = ⎛
⎝ − − − − + ∇ ∇ ⎞
⎠
+ −
u u u u x
J d d ρ d ψ ψ G
l d l d d
( , ̇, , ) 1
2 ̇· ̇ (1 ) 3
8 ( ( · )) d .
e e
Ω
2 c 2
(8) To account for irreversibility of the crack evolution a historyfield is used in place ofψe+, such that u t( , )is the maximum positive strain energy experienced (at time0⩽τ⩽t)[45],
=
∈
u t ψ+ u t ( , ) max ( , ).
τ t e
[0, ]
The historyfield 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 functionalJ( , ̇, ,u u d ∇d), with integration over Ω, now depends on two independent parameters and their respective deri- vatives:uandu̇, and d and∇d. Using the principle of least action, wefirst obtain the Euler-Lagrange equations foruwith respect to time. Then, using the principle of minimisation of potential energy, we obtain the Euler-Lagrange equations for d, for each dis- placement. This variation of Eq.(8)gives the governing equations in the strong form,
⎧
⎨⎩
∇ =
− − + ∇ =
σ ρu
d d
· ¨ (a)
2(1 ) G 0. (b)
l G l 3
8 3
4 2
c c
(9) Here we have neglected any body forces. The equations of motion also require boundary and initial conditions (nis the normal to the boundary),
⎧
⎨
⎪⎪
⎩
⎪⎪
= ∂
= ∂
∇ = ∂
=
= u U σn T
n
u u u x
u u v x
d t t
on Ω (a)
on Ω (b)
· 0 on Ω (c)
( , ) ( ) in Ω (d)
̇ ( , ) ( ) in Ω. (e)
u T
0
0 (10)
2.2. 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),
∫
⎛∫
⎝ + ∇ ∇ ⎞
⎠ = ⎛
⎝ − ⎞
⎠
x x
dδd G l
d δd G
l δd
2 3
4 ( ) d 2 3
8 d .
Ω
c
Ω
c
Discretising, such thatd( )x =N x d( ) i, ∇d( )x =B x d( ) i, δd( )x =N x( )δdiand ∇δd( )x =B x( )δdiwhere diare the nodal values of the crack phasefield d and N and = ∂ ∂B N/ xare some standardfinite element basis functions and their first spatial derivative, and observing that the test function is arbitrary, we get the discretised equation for the crack phasefield as,
=
K dd i f ,d (11)
where
∫
= ⎛
⎝ + ⎞
⎠
K N N G lB B x
2 3
4 d ,
d T T
Ω
c
and
∫
= ⎛
⎝ − ⎞
⎠
f G N x
2 3l
8 d .
d T
Ω
c
Since a staggered algorithm is used, there is no coupling between the displacement and crack phasefields.
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),
∫
Ω[(∇δu) ((1T −d)2 ++ −)∇ +u δ ρu u¨]dx=∫
∂ΩTδuTd( Ω ).∂ TDiscretising, such thatu x( )=N x u( ) i,∇u x( )=B x u( ) i,δu x( )=N x( )δuiand ∇δu x( )=B x( )δuiwhereuiare the nodal values of the displacementfieldu, and– again – observing that the test function is arbitrary, we get the discretised equation for the crack phasefield as,
+ =
Mu¨i Kui f, (12)
where
∫
=
M NTρN xd ,
Ω (13)
∫
= − ++ −
K BT((1 d) ) d ,B x
Ω
2 (14)
and
∫
= ∂
f ∂ N T d( Ω ).T T
ΩT (15)
Here, we have used a lumped mass matrix, obtained by using a quadrature only involving the nodal points (only forM).
2.3. 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 phasefield, Eq.(11)is solved, based on the displacements.
In the Newmark scheme, the displacement of the new time stepk+1 is predicted by
= + + −
ŭ+ u Δ ̇tu 1 β t u 2(1 )(Δ ) ¨ ,
k 1 k k 2 2 k
(16) whereΔtis the length of the time step, anduk, ̇ukandu¨kare the displacement, velocity and acceleration of the previous time step k.
The predicted displacement,ŭk+1, is used to solve the dynamic system of equations for the displacement problem (Eq.(12)) by
= −⎡
⎣ + ⎤
⎦ +
+
−
+ +
u¨ M 1β t K f Ku
2 (Δ ) ( ̆ ),
k 1 2 2 1 k k
1 1
(17) whereM, K and fk+1are 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 velocityu̇k+1and updated displacementuk+1as
= + − +
+ +
u̇k 1 u̇k (1 β1)Δ ¨tuk β1Δ ¨tuk 1
= +
+ + +
u ŭ 1β t u 2 (Δ ) ¨ .
k 1 k 1 2 2 k
1
Throughout this study, a fully explicit Newmark algorithm is used, with constantsβ1=1/2andβ2=0.
2.4. FE model
Finite element analyses are performed on strip specimen geometries with dimensions ( ×L H×T=50×25×2mm) (Fig. 2) [48,15]. A state of plane stress is assumed. Material properties are given inTable 1. The length parameter is taken to be of the same order of magnitude as the energetically consistent value, =l 3G Ec /(8σc2), whereσcis 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 lengthx0̂is inserted by specifying a strain historyfield (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 arefirst loaded quasi-statically (Δu=0.5·10−6m) until fracture nucleates, after which the crack propagation is simulated using an explicit Newmark algorithm. A time step ofΔt=2·10−8s is used;
this is about a third of the critical time step (the so-called CFL condition[49]),
⩽ t h
Δ c ,
p (18)
where h is the smallest element length andcpis 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 diameter2R=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 inFig. 3represent 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 distanceD/2from 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 dimensionsL×H×T=25×25×2mm) is used. The simulation parametersΔuandΔtare the same as before. A hole of radius eitherR=0.25mm orR=0.5mm is introduced. The parameter l is chosen such that =l 0.25mm, 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 4 in the vertical direction.l
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.
3. Results and discussion
3.1. Comparisons with experiments
For comparison, dynamic crack propagation experiments were conducted on pre-notched strip specimens ( ×L H×T=50×25×2mm) 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 paths– both experimental end numerical – which were obtained for a material without any holes, are shown Table 1
Material properties used for simulation.
E ν ρ Gc cR σc
3.24 GPa 0.35 1190 kg/m3 200 J/m2 962 m/s 50 MPa
Fig. 3. Geometries, white circles represent holes: (a) geometry 1, (b) geometry 2, (c) geometry 3 and (d) geometry 4. The dimensionsD=L/9and
=
R L
2 3 /100.
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,x0̂ =2mm, the right column to short initial cracks,x0̂ =0.75mm. (The length of the specimens is 50 mm.)
inFig. 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 inFig. 5.
The overall correlation between experimental and numerical crack paths is good. The general trend is that the longer cracks– i.e.
lower energy release rates– are 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 (x0̂ =0.75mm) 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).
3.2. Numerical results
3.2.1. Crack tip velocity
Fig. 7shows the normalised instantaneous crack tip velocityv c/̂ Rversus normalised crack tip positionx 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 positionx̂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 afictional 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 andfirst 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.x0̂ =2mm) 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
Fig. 5. Experimental and numerical crack paths for the different geometries, left columns: experimental results, right columns: numerical results.
Each row (1–4) corresponds to one geometry (geometry 1–4). Column 1 and 3 corresponds to long initial cracks ( ̂ =x0 2mm), columns 2 and 4 corresponds to short initial cracks (x0̂ =0.75mm). (The length of the specimens is 50 mm.)
crack (x0̂ =0.75mm), 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 tips– what 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 rate– where we have seen that the crack is less affected by the attraction exerted by the holes – the 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.
3.2.2. Energy release rate
Fig. 8shows the normalised instantaneous energy release rate G G/ cversus normalised crack tip positionx 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
̂ ̂
= − −
G dE
Tdx dE Tdx,
p k
whereEpandEkare the total potential end kinetic energies (sincex̂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 inFig. 8by 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.
3.2.3. 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 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 allfigures. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 9the inverse normalised energy release rate G Gc/ has been plotted versus the normalised velocityv c/̂ R. The dashed line shows the prediction by Eq.(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 velocityv̂increases.
Fig. 7. Normalised crack tip velocityv c/̂ Rversus normalised crack tip positionx L/̂ for 52 μs propagation, plotted over the crack path at =t 52 μs, for cracks with initial length, left:x0̂ =2mm and right:x0̂ =0.75mm. Each row (1–5) corresponds to one geometry: 1: reference, 2: geometry 1, 3:
geometry 2, etc. Grey lines are un-smoothed results, black lines have 5% smoothing. Red area indicates crack path. (For interpretation of the references to colour in thisfigure legend, the reader is referred to the web version of this article.)
Fig. 9(right) shows the same plot but with different axes. In fact, all geometries deviate from the prediction of Eq.(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].
Fig. 8. Normalised energy release rate G G/ c versus normalised crack tip positionx L/̂ for 52 μs propagation, plotted over the crack path at =t 52 μs, for cracks with initial length, left:x0̂ =2mm and right:x0̂ =0.75mm. Each row (1–5) corresponds to one geometry: 1: reference, 2: geometry 1, 3:
geometry 2, etc. Grey lines are un-smoothed results, black line have 5% smoothing. Red area indicates crack path. (For interpretation of the references to colour in thisfigure legend, the reader is referred to the web version of this article.)
3.3. 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 inFig. 10by a dashed line.
These specimens showed somewhat earlier onset of crack propagation compared to the reference material which is expected since we have introducedflaws 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
0 1 2 3 4
0 1 2 3 4
0 0.5 1 1.5
0 0.5 1 1.5
Fig. 9. Inverse normalised energy release rate G Gc/ versus normalised crack tip velocityv c/̂ Rfor different geometries and initial crack lengths. The dashed line represents the relation predicted by Eq.(1). The rightfigure is a closeup of the lower left corner of the left figure.
Fig. 10. Normalised crack tip velocityv c/̂ R(left) and normalised energy release rate G G/ c(right) versus normalised crack tip positionx L/̂ for 52 μs propagation, plotted over the crack path at =t 52 μs, for cracks with initial lengthx0̂ =0.75mm for three different random hole configurations. Grey lines are un-smoothed results, black line have 5% smoothing.
materials, but also provides a clue to the high strength often exhibited by porous materials.
3.4. Effect of internal characteristic length l
The effect of the internal characteristic length l with respect to pore radius R is shown inFig. 11. We see that when the hole is located on the reference crack path (i.e. crack path in a geometry without any hole), the crack path is unaffected by the hole, whether its radius isR=lorR=2 . When the hole is placed ahead of the propagating crack, but not in the reference crack path (due to crackl branching), a hole of radiusR=ldoes not affect the crack path; such a small defect is too small to localise damage, due to the regularisation. A hole of radiusR=2 on the other hand, is large enough to localise damage when the crack tips become shielded byl the hole, and the continued crack propagation occurs from this hole. This does however not mean the a crack is unaffected by small disturbances. When the hole is offset vertically, the small hole ( =R l) is still able to attract the crack and cause a crack pattern which is different from the reference. Also crack tip velocity and energy release rates are affected by the small hole (Fig. 11).
4. Conclusions
We have addressed the effect of holes and porosity on the dynamic crack propagation behaviour in a material containing pores. In general, a propagating crack is attracted to holes, and this attraction is stronger when the energy release rate is high, i.e. when there is more energy available for crack growth. When the energy release rate is lower, however, the effect of pores at remote distance from the crack is small.
In addition to the attraction exerted by the holes, the crack is affected by the shielding effect of the holes. When the energy release rate is low (in these cases approximatelyGc<G⩽2Gc) the shielding effect seems to dominate over the attraction exerted by distant holes. When the energy release rate is high (in these cases approximatelyG>2Gc), the attraction of the holes instead dominates over the shielding effect.
Crack propagation velocity in a material containing pores can be both lower and higher compared to a solid reference material, while the energy release rate is typically lower; the latter is unsurprising since the strength is lowered by the introduction of pores.
Fig. 11. Crack paths (upper), normalised crack tip velocityv c/̂ R(middle) and energy release rate G G/ c(bottom) for three different cases: no hole, hole of radiusR=land hole of radiusR=2l(overlaid). Three different geometries are studies, (a) a hole on the reference (no hole) crack path, (b) a symmetrically placed hole not in the reference crack path and (c) a non-symmetrically placed hole not in the reference crack path.
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 work– as with any modelling work – a 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 is– in 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 parameter– but 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 damagefield, 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 mi- crostructure, which supports thefirst 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.
Acknowledgement
The authors wish to thank Adele Wallin for help with the experimental work. This work was supported by the Swedish Energy Agency [Grant No. 37206-2] and the Swedish Research Council [Grant No. 2016-04608]. Per Isaksson greatly acknowledgesfinancial support from the European Union [H2020 MSCR-RISE Proj. No. 734485].
References
[1] Rosakis AJ, Samundrala O, Coker D. Cracks faster than the shear wave speed. Science 1999;284:1337–40.
[2] Schlüter A, Kuhn C, Müller R, Gross D. An investigation of intersonic fracture using a phasefield model. Arch Appl Mech 2016;86(1–2):321–33.
[3] Bleyer J, Roux-Langlois C, Molinari J-F. Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity- toughening mechanisms. Int J Fract 2017;204(1):79–100.
[4] Yilmaz O, Bleyer J, Molinari J-F. Influence of heterogeneities on crack propagation. Int J Fract 2018;209(1–2):77–90.
[5] Persson J, Isaksson P. Modeling rapidly growing cracks in planar materials with a view to micro structural effects. Int J Fract 2015;192:191–201.
[6] Gibson LJ. Biomechanics of cellular solids. J Biomech 2005;38(3):377–99.
[7] Fleck NA, Deshpande VS, Ashby MF. Micro-architectured materials: past, present and future. Proc Roy Soc A: Math Phys Eng Sci 2010;466(2121):2495–516.
[8] Freund LB. Dynamic fracture mechanics. Cambridge: Cambridge Univ Press; 1990.
[9] Carlsson J, Isaksson P. Dynamic crack propagation in woodfibre composites analysed by high speed photography and a dynamic phase field model. Int J Solids Struct 2018;144–145:78–85.
[10]Fineberg J, Bouchbinder E. Recent developments in dynamic fracture: some perspectives. Int J Fract 2015;196(1–2):33–57.
[11]Goldman T, Livne A, Fineberg J. Acquisition of inertia by a moving crack. Phys Rev Lett 2010;104(11).
[12]Sharon E, Fineberg J. The dynamics of fast fracture. Adv Eng Mater 1999;1(2):119–22.
[13] Hauch JA, Marder MP. Energy balance in dynamic fracture, investigated by a potential drop technique. Int J Fract 1998;90(2):133–51.
[14] Ravi-Chandar K, Knauss WG. An experimental investigation into dynamic fracture: II. Microstructural aspects. Int J Fract 1984;26:65–80.
[15] Nilsson F. Crack propagation experiments on strip specimens. Eng Fract Mech 1974;6(2):397–403.
[16] Sagy A, Cohen G, Reches Z, Fineberg J. Dynamic fracture of granular material under quasi-static loading. J Geophys Res: Solid Earth 2006;111(B4).
[17]Durr N, Sauer M, Hiermaier S. Mesoscale investigation of dynamic fracture in quartzite and sandstone and homogenization to macroscale. Int J Solids Struct 2018;144–145:160–79.
[18]Tarasov BG, Sadovskii VM, Sadovskaya OV, Cassidy MJ, Randolph MF. Modelling the static stress-strain state around the fan-structure in the shear rupture head.
Appl Math Model 2018;57:268–79.
[19]Tarasov BG, Guzev MA, Sadovskii VM, Cassidy MJ, Tarasov BG, Cassidy MJ, Guzev MA, Sadovskii VM, et al. Modelling the mechanical structure of extreme shear ruptures with friction approaching zero generated in brittle materials. Int J Fract 2017;207:87–97.
[20]Nguyen TT, Yvonnet J, Zhu QZ, Bornert M, Chateau C. A phasefield method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure. Eng Fract Mech 2015;139:18–39.
[21]Nguyen TT, Yvonnet J, Zhu QZ, Bornert M, Chateau C. A phase-field method for computational modeling of interfacial damage interacting with crack propa- gation in realistic microstructures obtained by microtomography. Comput Meth Appl Mech Eng 2016;312:567–95.
[22]Yilmaz O, Molinari J-F. A mesoscale fracture model for concrete. Cement Concr Res 2017;97:84–94.
[23]Du X, Jin L, Ma G. Numerical simulation of dynamic tensile-failure of concrete at meso-scale. Int J Impact Eng 2014;66:5–17.
[24]Lipperman F, Ryvkin M, Fuchs MB, Lipperman F, Ryvkin M, Fuchs MB. Fracture toughness of two-dimensional cellular material with periodic microstructure. Int J Fract 2007;146:279–90.
[25]Ryvkin M, Aboudi J. Crack resistance in two-dimensional periodic materials of medium and low porosity. Eng Fract Mech 2011;78(10):2153–60.
[26]Leguillon D, Piat R. Fracture of porous materials– influence of the pore size. Eng Fract Mech 2008;75(7):1840–53.
[27]Ambati M, Gerasimov T, De Lorenzis L. A review on phase-field models of brittle fracture and a new fast hybrid formulation. Comput Mech 2015;55(2):383–405.
[28]Francfort G, Marigo JJ. Revisiting brittle fracture as an energy minimization problem. J Mech Phys Solids 1998;46(8):1319–42.
[29]Bourdin B, Francfort G, Marigo JJ. The variational approach to fracture. Springer Verlag; 2008.
[30]Griffith AA. The phenomena of rupture and flow in solids. Philos Trans Roy Soc A: Math Phys Eng Sci 1921;221(582–593):163–98.
[31]Mumford D, Shah J. Optimal approximations by piecewise smooth functions and associated variational problems. Commun Pure Appl Math 1989;42(5):577–685.
[32] Bonnet A, David G. Cracktip is a global mumford-shah minimizer. Société Mathématique de France 274.
[33]Bourdin B, Francfort G, Marigo JJ. Numerical experiments in revisited brittle fracture. J Mech Phys Solids 2000;48(4):1–23.
[34]Ambrosio L, Tortorelli VM. Approximation of functional depending on jumps by elliptic functional via t-convergence. Commun Pure Appl Math 1990;43(8):999–1036.
[35]Pham K, Amor H, Marigo JJ, Maurini C. Gradient damage models and their use to approximate brittle fracture. Int J Damage Mech 2011;20(4):618–52.
[36]Larsen CJ, Ortner C, Süli E. Existence of solutions to a regularized model of dynamic fracture. Math Models Meth Appl Sci 2010;20(7):1021–48.
[37]Bourdin B, Larsen CJ, Richardson CL. A time-discrete model for dynamic fracture based on crack regularization. Int J Fract 2011;168(2):133–43.
[38]Borden MJ, Verhoosel CV, Scott MA, Hughes TJR, Landis CM. A phasefield description of dynamic brittle fracture. Comput Meth Appl Mech Eng 2012;217:77–95.
[39]Hofacker M, Miehe C. Continuum phasefield modeling of dynamic fracture: variational principles and staggered FE implementation. Int J Fract 2012;178(1–2):113–29.
[40]Kuhn C, Müller R. A discussion of fracture mechanisms in heterogeneous materials by means of configurational forces in a phase field fracture model. Comput Meth Appl Mech Eng 2016;312:95–116.
[41]Braides A. Approximation of free-discontinuity problems. Lecture notes in mathematics 1694. Berlin, Heidelberg: Springer; 1998.
[42]Amor H, Marigo JJ, Maurini C. Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments. J Mech Phys Solids 2009;57(8):1209–29.
[43]Pham K, Marigo J-J. Construction and analysis of localized responses for gradient damage models in a 1D setting. Vietnam J Mech 2009;31:233–46.
[44]Marigo J-J, Maurini C, Pham K. An overview of the modelling of fracture by gradient damage models. Meccanica 2016;51(12):3107–28.
[45]Miehe C, Hofacker M, Welschinger F. A phasefield model for rate-independent crack propagation: robust algorithmic implementation based on operator splits.
Comput Meth Appl Mech Eng 2010;199(45–48):2765–78.
[46] Matlab 2017a, the MathWorks, Natick, MA, USA.
[47]Newmark NM. A method of computation for structural dynamics. Proc Am Soc Civil Eng 1959;85(3):67–94.
[48]Nilsson F. Dynamic stress-intensity factors forfinite strip problems. Int J Fract Mech 1972;8(4):403–11.
[49]Courant R, Lewy H, Friedrichs K. Über die partiellen Differenzengleichungen der mathematischen Physik. Math Ann 1928;100:32–74.
[50]Marder M. Molecular dynamics of cracks. Comput Sci Eng 1999;1(5):48–55.
