On Crack Dynamics in Brittle Heterogeneous Materials

(1)

UNIVERSITATISACTA UPSALIENSIS

UPPSALA

Digital Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 1917

On Crack Dynamics in Brittle Heterogeneous Materials

JENNY CARLSSON

ISSN 1651-6214 ISBN 978-91-513-0906-4

(2)

Dissertation presented at Uppsala University to be publicly examined in Polhemsalen, Ångströmlaboratoriet, Lägerhyddsvägen 1, Uppsala, Friday, 15 May 2020 at 09:15 for the degree of Doctor of Philosophy. The examination will be conducted in English. Faculty examiner: Associate Professor Ralf Denzer (Lund University).

Abstract

Carlsson, J. 2020. On Crack Dynamics in Brittle Heterogeneous Materials. Digital Comprehensive Summaries of Uppsala Dissertations from the Faculty of Science and Technology 1917. 55 pp. Uppsala: Acta Universitatis Upsaliensis. ISBN 978-91-513-0906-4.

Natural variation, sub-structural features, heterogeneity and porosity make fracture modelling of wood and many other heterogeneous and cellular materials challenging. In this thesis, fracture in such complicated materials is simulated using phase field methods for fracture. Phase field methods have shown promise in simulations of complex geometries as well as dynamics and require few additional parameters; only the material toughness and a length parameter, determining the width of a regularised crack, are needed.

First, a dynamic phase field model is developed and validated against experiments performed on homogeneous brittle polymeric materials, wood fibre composites and polymeric materials with different hole patterns. Then, a high-resolution model of wood is developed and related to experiments, this time without considering fracture. Attention is finally focussed on high- resolution numerical analyses of fracture in wood and other cellular microstructures, considering both heterogeneity and relative density.

The phase field model is found to reproduce crack paths, velocities and energy release rates well in homogeneous samples both with and without holes. In more complicated heterogeneous and porous materials, the model is also able to simulate crack paths, but the interpretation of the length scale is complicated by the inherent lengths of the micro-structural geometry. In sum, the thesis points to possibilities with the proposed method, as well as limitations in our current understanding of both quasistatic and dynamic fracture of heterogeneous and cellular materials. The findings of this thesis can contribute to an improved understanding of fracture in such materials.

Jenny Carlsson, Department of Materials Science and Engineering, Applied Mechanics, Box 534, Uppsala University, SE-751 21 Uppsala, Sweden.

urn:nbn:se:uu:diva-406919 (http://urn.kb.se/resolve?urn=urn:nbn:se:uu:diva-406919)

(3)

List of papers

This thesis is based on the following papers, which are referred to in the text by their Roman numerals.

I Carlsson, J. and Isaksson, P. (2018) Dynamic crack propagation in wood ﬁbre composites analysed by high speed photography and a dynamic phase ﬁeld model. International Journal of Solids and Structures 144–145:78–85.

II Carlsson, J. and Isaksson, P. (2019) Crack dynamics and crack tip shielding in a material containing pores analysed by a phase ﬁeld method. Engineering Fracture Mechanics 206:526–540.

III Carlsson, J., Heldin, M., Isaksson, P. and Wiklund, U. (2019)

Investigating tool engagement in groundwood pulping: ﬁnite element modelling and in-situ observations at the microscale. Holzforschung, Epub ahead of print.

IV Carlsson, J. and Isaksson, P. (2020) Simulating fracture in a wood microstructure using a high-resolution dynamic phase ﬁeld model.

Accepted for publication in Engineering Fracture Mechanics.

V Carlsson, J. and Isaksson, P. (2020) A statistical geometry approach to length scales in phase ﬁeld modelling of fracture and strength of porous microstructures. Submitted.

Papers I and II are published with open access. Paper III is republished with permission from the publisher.

(4)

(5)

1. Introduction

In solid mechanics, materials are typically considered to be continuous and homogeneous, which means that we assume that if we take an inﬁnitely small piece of the material, that piece will have the same properties as the original specimen, and also that each small piece will be identical to any other (Cauchy, 1823). What is special about heterogeneous materials is that if we were to take a small piece of the material it will neither have the same properties as the large-scale material, nor will any piece necessarily be the same as any other.

There might be pores in the material, meaning that some such pieces will be empty while others contain material which is both stiffer and heavier than the average of the large-scale material. Another source of heterogeneity is when a large-scale material consists of different materials on the micro-scale, i.e. they are composites. A small piece of the large-scale material will then consist of either one of these materials, and any two pieces will not be alike, nor will they be like the large-scale material.

In many applications it is sufﬁcient to use the average properties at the larger scale (Ostoja-Starzewski, 2007). Fracture, however, is a local process, which occurs in a very small region near the crack tip. For this reason, local properties can play an important role in the fracture behaviour of heterogeneous materials (Hossain et al., 2014; Kuhn and Müller, 2016). Examples of heterogeneous materials are wood, bone, glass- and carbon ﬁbre composites, concrete, solid foams etc. The emphasis in this thesis is on wood, but many results apply to heterogeneous materials in a more general sense.

Wood is a very old – perhaps the oldest – construction material. Wooden water wells, dating more than seven thousand years back, have been found (Tegel et al., 2012; Rybníˇcek et al., 2020). But naturally, the use of wood goes even further back; for example, the numerous founds of axes dating from the stone age imply that wood was being worked regularly in these early societies.

Wood is an extremely versatile material, and some other uses include art, tools, clothes, food products and not least paper. In the 19th century, the invention of the first large-scale process to produce paper from wood, the groundwood process, was the starting point of a rapid decrease in paper prices (Müller, 2014). The groundwood process can thus, at least partly, be credited with the vast increase in information spread during the 20th century: schoolbooks, fiction as well as non-fiction books, letters, newspapers etc. And today, when information to an increasing extent is digital, paper usage in Europe remains near the all-time high because of increased demands for paper for packaging (CEPI, 2019).

(8)

In spite of its widespread use, our understanding of wood is still limited. A part of this can be attributed to the extreme diversity of plants; plants have been on Earth for over a billion years and outnumber animals in weight by a factor of one thousand (Jahren, 2017). Trees are youngsters in the plant family, only 370 million years, but they still make up the major part of the total biomass on land. This abundance makes for great variation. And so it follows that wood is a natural material with large variability: between different species of wood, between different individuals of the same species and even between different parts of the same tree.

But our limited understanding of wood is not only due to this so-called natural variability. Wood, like bone but unlike most man-made materials, is a multiscale heterogeneous material with distinct features from the scale of nanometers to that of a whole tree (Bodig and Jayne, 1982). The cell walls consist of crystalline cellulose microfibrils embedded in hemicellulose and lignin. The proportions of these constituents, the orientations of the cellulose fibrils, the thickness and shape of the cell walls and distribution of cells in different directions all influence the macroscopic properties of wood, making it difficult to model and therefore also difficult to predict its behaviour.

The specific topic of this thesis is modelling of fracture and crack dynamics. In short, fracture is the separation of an object or material into two or more pieces. A more technical definition is that fracture is the process of energy dissipation through the creation of new surfaces; these new surfaces are what we refer to as a crack (Griffith, 1921). This latter definition is very close to the theoretical framework of this thesis, phase field models of fracture based on the variational principle. In this framework, the energy that is dissipated through fracture is considered as a term of the total energy of the system. To- gether with well-established physical principles, such as the minimisation of potential energy or the principle of least action, and mathematical and numerical optimisation methods, such as the finite element method (FEM), fracture is predicted solely on the basis of stationary points of the total energy.

Models of fracture are useful for predicting the behaviour of materials without having to perform a lot of destructive testing, or at scales not typically available for mechanical testing. This knowledge is requested by e.g. the construction industry; wood construction is a hot topic from a sustainability point of view. Another useful application is that of mechanical pulping, where wood is repeatedly fractured until what remains is only a mixture of individual ﬁbres (the pulp). The work also provides insights into fracture in general porous and heterogeneous materials, which can be used to improve the integrity of e.g.

structural panels (sandwich panels), or in the design of medical implants, e.g.

bone implants or bone cements, which require a good understanding of the fracture of the porous cancellous bone structure for optimal strength and dura- bility.

(9)

1.1 Aims of the thesis

The aim of this thesis is to increase the understanding of dynamic fracture in wood and other heterogeneous materials on the cell scale, and more speciﬁ- cally to identify what features need to be considered when modelling this type of fracture. This is done by

• proposing a model for simulating dynamic fracture and validating this model in comparison with experiments;

• proposing a cell-scale model of wood and relating this model to experi- ments;

• identifying key features in the modelling by investigating the inﬂuence of structural features such as internal voids and interfaces on the model predictions; and

• investigating the possibility of generalising the model by going up and down in scales.

(10)

2. Heterogeneous and cellular materials

2.1 Wood

Wood is a cellular, heterogeneous, highly anisotropic material. The woods primarily of interest in this thesis are softwoods (or needlewoods). These consist almost solely of tracheids (longitudinally oriented cells), which are the primary structural building blocks of the wood, while at the same transporting fluids and nutrients, Fig. 2.1. This makes the structure somewhat simpler than that of hardwoods (or broadleaf woods), which have more rays (radially oriented cells) and sap channels (longitudinal cells for fluid transport). Although rays exist also in softwoods they are thin (often one cell wide and a few cells in height) compared to those in hardwoods. Wood without any complicated features, e.g. knots and rays, is referred to as clear-wood. The specific woods considered in the thesis are Scandinavian softwoods: Norway spruce (Picea Abies) and Scotch pine (Pinus Sylvestris). These woods have an average rela- tive density (defined as overall density divided by the density of the cell wall material) of around 0.3–0.4, but since they grow at high latitudes, where the climate is characterised by strong seasonal variation, they also exhibit pro- nounced annual ring structures, meaning that variations in relative density of between 0.1 and 0.6 are common (Thuvander et al., 2000). Both the average and the local relative density affect the overall properties of the wood, such as stiffness, density and strength.

A typical softwood consists of 90% tracheid cells (also referred to as fibres or simply cells). In softwoods, tracheid cell radii (r_cell) are typically around 10 – 40 μm (Bodig and Jayne, 1982). Typical fibre lengths are in the order of a few millimetres. Each tracheid cell consists of several layers, namely the primary wall (P), the three layers of the secondary wall (S) and a void called the cell lumen, Fig. 2.2. The secondary wall, and especially the S2 layer, is the thickest. Due to the high proportion of cellulose, it contributes most to the stiffness of the cell wall (Harrington et al., 1998), Table 2.1–2.2. Between the cells is the middle lamella, the “glue” that holds the cells together, which is weaker due to a high proportion lignin. In addition, it is often difficult to distinguish between the middle lamella and the primary wall, which is why this compound is often referred to as the compound middle lamella (CML).

With 90% of the cells in the longitudinal direction, wood is naturally anisotropic. On a local scale (the scale of a few cells), it can be considered transversely isotropic with the tangential and radial directions approximately equal, Table 2.3. On a global scale, the stiffening effect of rays typically makes the radial direction somewhat stiffer than the tangential.

(11)

Figure 2.1. Rendering of CT micrograph showing the structure of wood and the global, trunk-centered coordinate system (uppercase).

Figure 2.2. The different parts of a wood cell. The lowercase coordinate system refers to the local, cell-centered coordinate system.

Table 2.1. Approximate proportions of the constituents of the different cell wall layers (Harrington et al., 1998).

Cellulose Hemicellulose Lignin

S1 wall 37% 35% 28%

S2 wall 45% 34% 21%

S3 wall 34% 36% 30%

CML 11% 14% 75%

(12)

Table 2.2. Young’s modulus E and ultimate strainεcof the cell wall constituent ma- terials (Gibson et al., 2010).

E εc

Cellulose 165 GPa 20–40%

Hemicellulose 9 GPa

Lignin 3 GPa 3–4%

Table 2.3. Young’s moduli E, Poisson’s ratiosν and shear moduli μ for wood on the cell scale (Harrington et al., 1998). Moduli in GPa. The quantities are given in the cell-centered coordinate system (cf. Fig. 2.2).

El Er Et νlr νlt νrt μlr μlt μrt

S1 wall 53 8 9 0.3 0.3 0.4 3 3 3

S2 wall 64 9 10 0.3 0.3 0.4 3 3 3

S3 wall 50 8 8 0.3 0.3 0.4 3 3 3

CML 18 5 5 0.3 0.3 0.4 2 2 2

2.1.1 Fracture of wood

When loaded in compression, wood does typically not fracture but collapse (Gibson and Ashby, 1997). This is contrary to most brittle materials, e.g.

concrete, which fracture transversely when loaded in compression. In tension however, wood fractures. Due to the highly orthotropic nature of wood, the different fracture directions show different behaviours. Consequently, the different crack propagation directions are often referred to as different fracture systems, Fig. 2.3. The ﬁrst letter denotes the normal direction of the plane of fracture and the second letter the direction of propagation. Sometimes prop- agation in the R direction is divided further into R⁺ and R⁻, depending on whether the crack is propagating away from or towards the pith (the centre of the stem). In this thesis, the fracture systems under consideration are RT and T R.

Figure 2.3. Fracture planes of wood. The ﬁrst letter refers to the direction normal to the crack faces, the second letter to the direction of propagation.

(13)

Figure 2.4. Crack propagation in the tangential (RT) direction (left) and in the radial (TR) direction (right). Adapted from Dill-Langer et al. (2002).

As reported by Ashby et al. (1985), when a crack is propagating in high- density wood, such as latewood, crack growth typically occurs by separating cells in the middle lamella, so-called intercellular fracture. A transition occurs at a relative density of around 0.2. Below this value, cracks tend to propagate primarily through the lumen, i.e. by breaking cell walls, so-called trans-wall fracture; this transition can be related to the relative strengths of the cell wall and the middle lamella (Gibson and Ashby, 1997). The transition between cell separation and cell wall breaking also depends on the architecture of the microstructure. Crack propagation in the pure radial direction is typically straight and occurs by separation of cells in the middle lamella, whereas crack propagation in the tangential direction occurs with a combination of cell separation in the middle lamella and cell wall fracture (Dill-Langer et al., 2002), Fig. 2.4.

Because of the variation in relative density, cracks propagating in the radial direction tend to propagate in a stepwise fashion, with the crack arresting as it approaches the denser late-wood (Stanzl-Tschegg et al., 2011; Thuvander and Berglund, 2000). With further buildup of available strain energy, the crack will break through the late-wood, and due to the excess energy this growth may be unstable. Crack propagation in the RT -direction, on the other hand, is typically stable on a global scale, as the material is more homogeneous in this direction.

2.2 Other cellular materials

The cellular structure of wood, especially the alignment of elongated, stiff cells in a weaker matrix, is similar to that of cortical bone (Gustafsson, 2019).

But also without this very speciﬁc structure, many cellular materials, e.g. cancellous bone or solid foams, have features in common with wood, which are not shared with general homogeneous materials.

2.2.1 Solid foams

Solid foams are typically divided into open-cell and closed-cell foams. Open- cell foams resemble a three-dimensional lattice, or truss structure. Close-cell

(14)

foams typically resemble a union of bubbles, e.g. bread. For both kinds of foams, stiffness and strength depend strongly on the relative density. Like wood, foams do not fracture under compression, but collapse. Only in tension is fracture commonly observed. Due to the large possibility for variation in foams, the results of this thesis do not apply to all foams, but can be most likely be applied to transversely isotropic, open-cell or thin-walled closed- cell foams (i.e. most foams foamed from liquids) made of relatively brittle materials.

2.2.2 Bone

Most bones are made up of an outer shell of relatively dense cortical bone and a core of porous cancellous, or trabecular, bone. The human cortical bone has an architecture similar to wood – longitudinally oriented osteons (bone cells) with a central canal, surrounded by a weaker matrix. However, the relative density of cortical bone is typically much higher than that of wood (>0.7).

Cancellous bone, on the other hand, has a relative density similar to that of wood (0.05–0.7). Cancellous or trabecular bone is made up of a network of connected rods or plates (trabeculae). Lower-density trabecular bone contains primarily rods, whereas higher-density trabecular bone contains more plates (Gibson and Ashby, 1997). In general, the behaviour of trabecular bone is typical of cellular materials; stiffness and strength depend strongly on the relative density and it tends to collapse under compression but can fracture in tension.

2.3 Effective properties of cellular materials

The effective, or macroscale, properties of porous materials typically depend on the relative density of the material. If the solid material making up the walls of the porous material is considered transversely isotropic and linearly elastic, and if the architecture of the material is isotropic, then on the macroscale (con- tinuum scale), we have the macroscale stressσ = Eε, and on the microscale we have the microscale stressσ_μ = E_με_μ, where E is Young’s modulus and ε the strain and the index μ is used to indicate microscale properties. Ac- cording to Gibson and Ashby (1997), wood with its elongated circular, hexag- onal or rectangular cells, can be considered as a slightly perturbed honeycomb material. For a honeycomb material, the effective stiffness and stress on the macroscale depend on the relative density ρ^∗ (deﬁned as the average density of the material divided by micro-scale density, i.e. ρ^∗ = ρ/ρμ) as E_{3rd order} ∝ (ρ^∗)³E_μ. Numerical experiments performed in Paper V suggest that a quadratic scaling is more accurate for the materials considered in this thesis, i.e.

E2nd order∝ (ρ^∗)²E_μ.

(15)

Table 2.4. Scaling relationships for stiffness E, toughness G_c, strainε and stress σ depending the relative density. Numerical experiments performed in Paper V show that the second order model gives the most accurate predictions for the materials studied in this thesis.

E G_c ε σ

Second-order (ρ^∗)²E_μ ρ^∗G_c_μ ε_μ/ρ^∗ ρ^∗σ_μ

The strain energy postulate, stating that the strain energy of the macro- and microscale models coincide, gives expressions for macroscale stress and strain.

In addition, the energy dissipated by fracture must be equal on the macro- and microscales. Since the dissipated energy is proportional to the surface area of the crack on the macro scale, and the surface area is scaled by ρ^∗, the frac- ture toughness must scale as G_c≈ ρ^∗G_cμ. The scaling relationships for the second-order stiffness-scaling model are summarised in Table 2.4.

2.4 Critical crack length of a cellular microstructure

Cellular materials are typically less sensitive to small cracks than continuous materials are. A random porous microstructure can be described by a spatial Poisson point process (Okabe et al., 2000), Fig. 2.5. The average number of cell mid-points in an area A is a spatial Poisson point process with intensity parameter (i.e. number of occurrences per unit area)λA. The probability that the number of points n occupying an area A is exactly equal to k is

P{n = k} = λ_A^k k!e^−λ^A^A.

A lot of useful results have been derived for the Poisson Voronoi tessella- tion, such as the expected total edge length per unit area, L_A= 2√

λA(Okabe et al., 2000).

A typical question that occurs when dealing with fracture of random cellular materials is: How large must a crack or other defect be in order to affect the geometry at all? A straight crack can be described by a sectional diagram of the Poisson Voronoi tessellation, Fig. 2.5. A line sectional diagram of a 2D Poisson Voronoi tessellation is a Poisson line process with intensityλL. Here, λLis the average number of cells encountered in a unit line distance along any direction (Okabe et al., 2000),

λL=2LA

π = 4√ λA

π .

The probability of encountering k cell walls along a distance L is thus P(n = k) =λ_L^k

k!e^−λ^L^L.

(16)

Figure 2.5. A Poisson Voronoi tesselation and its sectional diagram. Crosses indicate the seeds of the spatial process, the average number of seeds per unit area isλA. The line process is given as a sectional diagram of the original tessellation. For the line process, the intensities of walls (rings) and centres are the same,λL.

Suppose we introduce a crack of length a somewhere in the material. The probability that this crack does not intersect any cell wall is then

P{n = 0} = e^−aλ^L.

By setting P= 1/2 we can ﬁnd a critical crack length ac, for which it is equally likely that an introduced crack will interfere with the structure as it is that it will not,

a_c= − 1 λL

ln

1 2

. (2.1)

This is true only when the width of the cell walls is negligible. If this is not the case, then in order for a crack not to intersect any cell wall, it has to be at a distance wπ/4 away from any cell wall, where w/2 is the width of half of the cell wall, and the factor 2/π is the mean of sine over [0,π]. This is equivalent to stating that

P{x = 0} = e^{−(a+wπ/4)λ}^L. Setting P= 1/2, again, we ﬁnd a critical crack length,

ac= − 1 λL

ln

1 2

−wπ 4 .

The critical crack length can be rewritten as a function of the relative density ρ^∗using the expected total edge length per unit area L_A, as,

ac= − 1 λL

ln

1 2

− ρ^∗π 8√

λA

. (2.2)

This critical length can be related to the notch sensitivity of a porous material, as demonstrated in Paper V.

(17)

3. Models of fracture

3.1 Grifﬁth’s theory

According to the theory developed by Griffith (1921), fracture is a process of energy dissipation through creation of new surfaces. Griffith found that during stable crack propagation, there is no net change in total energy. LetΠedenote potential energy in the material, stored in the form of strain energy, andΠsthe surface energy of the crack. Moreover, letF be work done by external forces and let a denote the crack length. Griffith concluded that as the crack advances a distance da, the change in total energy is zero, i.e.

d(Πe+ Πs− F )

da = 0. (3.1)

Since the derivative in (3.1) is equal to zero, the energy in the numerator has a possible extremum with respect to crack length. Grifﬁth thus concluded that “[t]he ‘Theorem of minimum potential energy’ may be extended so as to be capable of predicting the breaking loads of elastic solids, if account is taken of the increase in surface energy which occurs during the formation of cracks”

(Grifﬁth, 1921).

Grifﬁth also showed that the product of the breaking strength σc and the square root of the crack length is constant,

σc

√a≈

2γE

π , (3.2)

The quantity 2γ is twice the surface energy of the material, i.e. the energy required to create two surfaces. It can further be noted that zero crack length predicts inﬁnite breaking strength.

3.2 Later development of linear elastic fracture mechanics

In the 1950’s, Irwin (1957) developed the concept of linear elastic fracture mechanics. Irwin uses the concept of stress intensity factors, which divides fracture into three modes, Fig. 3.1. For each mode, the stress at the crack tip is given by

σi j=

K_I_/II/III

√2πr

fi j(θ)

(18)

Figure 3.1. The three modes of fracture, from left to right: Mode I (opening), Mode II (in-plane shear) and Mode III (out-of-plane shear).

where KI, KII and KIII are the stress intensity factors of the three modes, the function f_{i j}(θ) is a trigonometric function of the angle relative to the crack plane (Williams, 1957), and r is the distance from the crack tip. Irwin and Orowan (1955) also divided the energy dissipated by fracture into two parts, G_c = 2γ + Gp where G_p is energy dissipated by plasticity. As long as the plastic zone, or process zone, remains small, it is possible to replace 2γ in (3.2) by G_c.

The Grifﬁth-Irwin relationship is G=K_I²

E ,

where the stress intensity factor K_Iand energy release rate G are replaced by the critical stress intensity factor KIc and the toughness Gc at the initiation of crack propagation.

In later work on plastic fracture Rice and Rosengren (1968), Cherepanov (1967) and Hutchinson (1968) independently developed the concept of the J- integral which is the equivalent of G.

3.3 Crack kinematics

An issue regarding the aforementioned theories is the inability to predict crack kinematics. In a body Ω ⊆ Rⁿ, a crackΓ ⊂ Ω is a subset of the body, typi- cally of dimension of n− 1. This means that a crack can be parametrised by n− 1 parameters, but require n functions of these parameters. Grifﬁth’s equa- tions are scalar-valued, thus the system is under-determined for n> 1. Irwin’s equations rely on the assumption that a crack is in one of the three modes.

And while there exist expressions for mixed-mode cracks, these do not predict crack propagation direction. So where will the crack propagate?

One hypothesis is that the crack will propagate in the direction of the maximum hoop stress. This hypothesis is usually credited to Erdogan and Sih (1963) but had already been used by Yoffe (1951) to predict crack branching in dynamic fracture. It follows from the assumption that cracks will orient themselves in the mode I direction.

Another common hypothesis is that the crack will propagate in the direc- tion of the maximum energy release rate, G. For the simplest case of pure

(19)

mode II loading, the two assumptions give slightly different predictions on the direction of crack propagation.

3.4 Dynamic fracture

In dynamic fracture, crack growth is not stable, i.e. the energy release rate is not necessarily equal to the change in elastic energy, cf. (3.1). The energy surplus becomes kinetic energy, causing stress waves to propagate through the specimen. Under some conditions, for example when slowly loading speci- mens without or with only a very short initial crack, it is possible to load more elastic energy into the structure than can be released by stable fracture; the resulting fracture must therefor be unstable, or dynamic. Moreover, the energy cannot be dissipated instantaneously, as crack surfaces cannot extend at unlimited velocity. A speed limit for mode I crack propagation is the Rayleigh velocity cR; theoretically, a mode I crack can propagate with velocity v ei- ther below c_Ror above the pressure wave velocity c_p(Broberg, 1989), i.e. the range cR≤ v ≤ cp is “forbidden”. In mode I crack propagation, the velocity typically relates to the energy release rate G as (Freund, 1990),

G

Gc ≈ 1 − v

cR. (3.3)

Equation (3.3) is true only for relatively low velocities, approximately v≤ 0.6cR. At higher velocities, the crack becomes unstable. Instead of accelerat- ing, it is energetically favourable for the crack to branch into two or more crack tips, each propagating with velocities around 0.6cR. This relation between energy release rate and branching was first suggested by Eshelby (1999) and has also been verified numerically using phase field models (Bleyer et al., 2017).

Broberg (1989) also showed theoretically that a mode II (shear) crack can propagate with a velocity under c_Ror over the shear wave velocity c_s. This was veriﬁed experimentally by Rosakis et al. (1999); their experiments have been reproduced using the phase ﬁeld method by Schlüter et al. (2016). Instanta- neous mode I crack propagation velocities as high as 0.9cRhave been observed experimentally (Sharon and Fineberg, 1999), supersonic mode I cracks have been observed only in atomistic simulations (Buehler et al., 2003).

3.5 The variational approach to fracture

The phase ﬁeld models used in this thesis are based on the variational approach to fracture of Francfort and Marigo (1998). The variational approach to fracture is closely related to the works of Grifﬁth, but resolves most of the drawbacks related to crack initiation and direction of propagation. In the variational approach to fracture, a family of possible crack paths is considered.

(20)

For each member of the family, there is an associated energyΠ, depending on the displacement uuu and the crackΓ,

Π(Γ,uuu) = Πe(Γ,uuu) + Πs(Γ), (3.4) where the elastic energy

Πe(Γ,uuu) =

Ωψe(Γ,∇uuu)dxxx, and the surface energy

Πs(Γ) = Gc

Γds.

Throughout this thesis, Π (regardless of index) denotes total energies, inte- grated over the entire bodyΩ whereas ψ denotes energy densities. Also, ∇ denotes the spatial differential operator and dxxx indicates integration over spa- tial coordinates whereas ds indicates integration along line segments given by Γ.

According to Francfort and Marigo (1998), the evolution of the crack and energy must follow some conditions. One such condition is that broken material does not heal, by requiring that the crack surface is non-decreasing. An- other condition is that the “real” crack must give a minimum with respect to energy compared to all other admissible cracks (i.e. cracks for which the previous cracks are enclosed in the current crack).

The crack state is found as the infimum of (3.4), but this infimum is gener- ally not trivial to find, in fact, it is most often impossible. Bourdin et al. (2000) noted the similarity between the variational fracture functional and the Mum- ford and Shah (1989) functional used for image segmentation. Thus, the first numerical experiments in the variational approach to fracture used this similarity to apply the already known Ambrosio and Tortorelli (1990) functional approximating the Mumford-Shah functional to fracture.

3.6 Phase ﬁeld models

In gradient-regularised phase field methods, a sharp crack is represented by a diffuse damage field. The distribution of this damage field is regularised by the gradient of the damage variable, thus dividing the surface energy into two terms, one local and one non-local (gradient-dependent) (Bourdin et al., 2000; Ambrosio and Tortorelli, 1990; Braides, 1998). The damage variable represents the degree of damage in a region of the material, and typically determines the degree of stiffness loss.

Phase ﬁeld methods for fracture have quickly become popular, which can be explained by their extreme versatility. Phase ﬁeld models have been suc- cessfully applied to dynamics (Bourdin et al., 2011; Larsen et al., 2010; Bor- den et al., 2012; Hofacker and Miehe, 2012; Schlüter et al., 2016), anisotropic

(21)

Figure 3.2. Regularised representation of a crack. A sharp crackΓ (left) is represented by a crack density functionγ(d,∇d) as a diffuse damage ﬁeld d (right).

surface energy (Teichtmeister et al., 2017; Li et al., 2015), heterogeneous materials (Hossain et al., 2014; Kuhn and Müller, 2016), plasticity (Hesch et al., 2017; Alessi et al., 2015) and fatigue (Alessi et al., 2018).

Consider a problem domainΩ ⊂ R², with exterior∂Ω, part of which ∂ΩT is subject to natural boundary conditions (prescribed stress TTT ) and part of which

∂Ωuis subject to Dirichlet boundary conditions (prescribed displacement UUU ), with a discrete crackΓ, Fig. 3.2. The total free energy of the system can be written as (Larsen et al., 2010; Larsen, 2010; Bourdin et al., 2011)

Π(uuu,d) =

Ωψk(˙uuu)dxxx −

Ωψe(∇uuu,d)dxxx − Gc

Γds. (3.5) Hereψkis the kinetic energy density. Through the phase ﬁeld implementation, the discrete crackΓ is represented by a crack density function γ(d,∇d), where d is a diffuse regularised crack ﬁeld. Throughout this thesis, it is assumed that d= 1 represents broken material and d = 0 represents intact material. In the works of this thesis, two of the most common crack density functions, the so-called Ambrosio-Tortorelli 1 (AT1) and Ambrosio-Tortorelli 2 (AT2), are used. The AT1 model is linear in d and quadratic in ∇d (Pham and Marigo, 2009; Pham et al., 2011)

Γds≈

ΩγAT1(d,∇d)dxxx =

Ω

3

8l(d + l²∇d · ∇d)dxxx. (3.6) The constant l is a regularisation parameter, a kind of characteristic length, determining the width of the regularised crack. The AT2 model is quadratic in both d and∇d (Bourdin et al., 2000; Miehe et al., 2010),

Γds≈

ΩγAT2(d,∇d)dxxx =

Ω

1 2l

d²+ l²∇d · ∇d

dxxx. (3.7) A number of different crack density functions exists, each one gives a somewhat different fracture behaviour, (cf. e.g. Braides, 1998; Pham et al., 2011;

Marigo et al., 2016).

The evolution of damage must fulﬁll some conditions, closely related to the conditions stipulated by Francfort and Marigo (1998), presented in the

(22)

previous section. Speciﬁcally, damage must be irreversible, the “real” state of the damage must be a minimum of the action integral (or the total energy functional) and the energy balance must always hold.

The kinetic energy is Πk(uuu) =

Ωψk(˙uuu)dxxx = 1 2

Ωρ ˙uuu · ˙uuudxxx,

where ˙uuu, and later occurring ¨uuu, denotes ﬁrst and second derivatives of displace- ment with respect to time t, and ρ is the density. In the case of quastistatic simulations, the velocity and acceleration terms are considered to be so small that they can be neglected, and this term is then zero.

The kinetic energy is unaffected by the crack phase ﬁeld d, but d locally degrades the elastic stiffness of the material. A degradation function of the type (1 − d)² is used in both the AT1 and AT2 models. Other degradation functions give a somewhat different response (cf. Marigo et al., 2016; Kuhn et al., 2015; Sargado et al., 2018).

In order to prevent crack growth under compression, and to account for crack surface over-closure, the strain energy is typically split into a positive part, which is degraded as damage increases, and a negative part, which is unaffected, ψe(∇uuu,d) = (1 − d)²ψe⁺(∇uuu) + ψe⁻(∇uuu) (cf. Amor et al., 2009;

Miehe et al., 2010). There are many different ways to assign the positive and negative strain energy, two are given in Section 3.6.1.

For the AT1 model, the action integral over time t∈ Θ = [0,tend] is J(uuu,d) =

Θ

Ω

1

2ρ ˙uuu · ˙uuu − (1 − d)²ψ_e⁺− ψ_e⁻−3G_c

8l (d + l²(∇d · ∇d))

dxxxdt,

(3.8) and for AT2 model

J(uuu,d) =

Θ

Ω

1

2ρ ˙uuu · ˙uuu − (1 − d)²ψ_e⁺− ψ_e⁻−Gc

2l

d²+ l²∇d · ∇d

dxxxdt.

(3.9) To account for irreversibility of the crack evolution, a history ﬁeld is used in place ofψe⁺(Miehe et al., 2010), so thatH (xxx,t) is the maximum positive strain energy experienced at each point xxx in the time 0≤ τ ≤ t,

H (xxx,t) = max

τ∈[0,t]ψe⁺(uuu(τ)).

The history ﬁeld ensures that it is the largest strain energy experienced in the material during the simulation history that determines the present stiffness.

This kind of history ﬁeld is known to induce a small error in the simulation, especially for small l (Bourdin et al., 2008; Linse et al., 2017). Stresses are however reversible, and are always evaluated asσσσ = ∂ψe/∂εεε, where the lin- earized strainεεε = 1/2(∇uuu + ∇^Tuuu). Taking a second derivative gives the con- sistent tangent stiffness tensorC = ∂²ψe/∂εεε².

(23)

Using the principle of least action, the Euler-Lagrange equations of (3.8–

3.9) are obtained for uuu and d as

⎧⎨

⎩

∇ · σσσ = ρ ¨uuu (3.10a)

2(1 − d)H −3G_c

8l +3G_cl

4 ∇²d= 0 (3.10b)

⎧⎨

⎩

∇ · σσσ = ρ ¨uuu (3.11a)

2(1 − d)H −G_c

l (d − l²∇²d) = 0. (3.11b) For quasistatic cases, we have ˙uuu≈ 0 and the governing equations are still obtained as the Euler-Lagrange equations of Eqs. (3.8–3.9). AssumingΠk= 0 and minimising the total energy rather than the action functional produces the same result. The equations of motion also require boundary and initial conditions (nnn is the normal to the boundary),

⎧⎪

⎪⎪

⎨

⎪⎪

⎪⎩ u

uu= UUU on∂Ωu (3.12a)

σσσnnn = TTT on ∂ΩT (3.12b)

∇d · nnn = 0 on ∂Ω (3.12c)

uuu(uuu,t) = u0(xxx) in Ω (3.12d) uu˙u(uuu,t) = v0(xxx) in Ω. (3.12e)

3.6.1 Strain energy splits

A common strategy in phase ﬁeld modelling is to split the strain energy density such that only tensile strain energy contributes to the damage. This way, materials do not fail under pure compression, and stiffness is kept in the case of crack closure (Bourdin et al., 2008; Amor et al., 2009; Miehe et al., 2010).

Many different split variants have been suggested, (Ambati et al., 2015; Strobl and Seeling, 2016). In the various works of this thesis, two of the most common splits are being used, the spectral split of Miehe et al. (2010) and the hydrostatic-deviatoric split of Amor et al. (2009).

The spectral split (Miehe et al., 2010) uses the fact that for an isotropic material, the principal strain system gives independent modes of deformation which can be split into positive and negative parts,

ψ_e⁺=α

2λtr(εεε)²+ μtr(εεε²₊), ψ⁻= 1− α

2 λtr(εεε)²+ μtr(εεε²₋).

(24)

Figure 3.3. 1D traction model.

where α = 1 if tr(εεε) > 0 and α = 0 otherwise and μ and λ are the Lamé parameters, whereλ= λ in plane strain and λ= λ(1−2ν)/(1−ν) in plane stress. The positive and negative strain contributions are given by a spectral decompostion of the strain tensor,

εεε_±=∑²

1

e_i± |ei|

2 nnn_iiinnn^T_iii, (3.13) where eiare the eigenvalues and nnn_iiithe eigenvectors of the strain tensorεεε.

Another possibility is to use the hydrostatic and deviatoric deformation modes to assign the strain energy to the positive and negative parts (Amor et al., 2009),

ψe⁺= α

2κtr(εεε)²+ μεεεdev:εεεdev. ψe⁻=1− α

2 κtr(εεε)²,

whereκ= λ+ μ is the bulk modulus (in 2D) and εεεdev= εεε −1/2tr(εεε)III is the strain deviator (and III the identity matrix).

3.6.2 Analytical solutions

For some simple quasistatic problems, such as traction of a one-dimensional bar, as well as for the transverse direction across a crack, it is possible to derive analytical solutions.

1D traction of a bar

Consider the barΩ ⊂ R of cross sectional area S, subject to end displacements U(0) = 0 and U(L) = tL (Carlsson, 2019), Fig. 4.17. The stiffness is the scalar stiffness E(1 − d)² and stress σ and strain ε are also scalar-valued. It is possible to drop the nabla notation and instead write d, as d only has one directional derivative. Note also thatε = u= t.

Letting γAT1(d,d) = 3/(8l)(d + l²d²) (cf. Eq. 3.6), i.e. the AT1 model, gives

Π(u,d) = S ^L

0

1

2E(1 − d)²t²+3G_c 8l

d+ l²d² dx.

On Crack Dynamics in Brittle Heterogeneous Materials

On Crack Dynamics in Brittle Heterogeneous Materials

List of papers

Contents

1. Introduction

2. Heterogeneous and cellular materials

3. Models of fracture