The Theory of Fracture — Deriving the Theoretical Strength, Griffith’s Equation, Fracture Toughness, and Weibull Statistics

This page is the companion explanation to fracture-simulator-en.html. The simulator showed “a gap of tens to hundreds of times between the theoretical strength and the actual strength,” “strength that falls as the −1/2 power of the crack length,” and “strength that scatters even when specimens are made under the same conditions, and the Weibull modulus m.” This page derives where these come from, without skipping intermediate steps. All you need is Hooke’s law, the elastic strain energy σ²/2E, finding a maximum by differentiation, and the first rule of probability that “the probability of independent events is a product.” The page is written with second- and third-year undergraduates in mind. The physics of the phenomena that the simulator makes visible — the necking of a tensile specimen, the stress field around an elliptical hole, the plastic zone at a crack tip and the atomic picture, the distribution of strengths that builds up as specimens are broken one after another, the stress distribution in a bending test, the temperature and thermal-stress distributions in thermal shock, and the Vickers indentation with the cracks running out from its corners — is also derived in §1.1–§1.5, §3.2–§3.5, §5.5–§6.7, §8.5, §9.5–§9.6, §10.2–§10.5, and §12 (for the complex functions of §3.2 and the numerical method of §10.4, it is enough to understand the idea and how to use the result).

1. Introduction — How Ceramics Break

Drop a teacup on the floor and it breaks. Drop a metal spoon and it may bend, but it will not break. A tensile test makes the difference clear: a ceramic breaks suddenly at some stress while stress and strain are still proportional. The strain just before fracture is only 0.092 % for alumina and 0.071 % for soda-lime glass. A metal, after yielding, stretches a great deal by plastic deformation and fails only after necking (§1.1–§1.5).

The root of this difference is whether dislocations can move. Metallic bonds have no directionality, and atomic planes slip over one another relatively easily. In ionic crystals, ions of the same sign come face to face partway through slip, which costs a lot of energy, and in covalent crystals the bonds are directional, so slip requires the bonds to be switched over. As a result, in polycrystalline ceramics at room temperature dislocations hardly move and there are too few independent slip systems (single crystals such as MgO do slip a little even at room temperature), so the crack tip is never rounded off by plastic deformation. A crack whose tip stays sharp therefore governs the strength.

The elastic energy stored per unit volume up to fracture is the area of the triangle under the stress–strain line, σ_f²/2E. It is 0.16 MJ m⁻³ for alumina, 0.018 MJ m⁻³ for glass, and even for Y-TZP only 2.4 MJ m⁻³. In a metal the work of plastic deformation is added, so it absorbs far more energy than this before it fails. This page uses the values in Table 1 as examples.

1.1 The tensile test — nominal stress and true stress

Tab ① of the simulator carries out a student-laboratory tensile test on the screen. A specimen whose parallel section has length (gauge length) L₀ and cross-sectional area A₀ is pulled with a load F, and the elongation ΔL is measured. What the testing machine records directly is F and ΔL, so the first quantities we obtain are these divided by the original dimensions.

However, as the specimen is pulled it becomes thinner, so the stress actually carried by the current cross section A is F/A. The strain, likewise, is defined by adding up the increments dL/L as the length grows from L to L + dL. These are the true stress σ_t and the true strain ε_t.

In the plastic deformation of a metal the atoms merely rearrange, and the volume hardly changes. Setting AL = A₀L₀, we have A = A₀L₀/L = A₀/(1 + e), so

In the elastic range the volume does change slightly, but the strain is only about 0.1 %, so the difference can be neglected. After necking begins, the cross section is no longer uniform along the specimen, so the true stress and true strain are taken from the cross section at the thinnest place A_min, as σ_t = F/A_min, ε_t = ln(A₀/A_min) (the readouts of the simulator use this definition, too).

1.2 Work hardening — Hollomon’s equation

Even after a metal yields, dislocations become tangled with one another and harder to move, so a larger stress is needed to keep deforming it (work hardening). The relation between true stress and true strain is well described for many metals by the following power function (Hollomon, 1945).

n is called the work-hardening exponent and K the strength coefficient (the K and n in this section are different quantities from the stress intensity factor K_I in §5 and the number of flaws n in §7). In the simulator, after the elastic line σ_t = Eε_t reaches the yield stress σ_y, we use σ_t = max(σ_y, Kε_tⁿ). Applying Eq. (4) as it stands at the yield strain σ_y/E gives 209 MPa for structural steel, which falls below σ_y = 250 MPa, so the stress is held constant at σ_y up to the strain ε₁ = (σ_y/K)^1/n at which Kε_tⁿ catches up with σ_y (0.36 % for steel, 0.47 % for the Al alloy). This is a formal stand-in for the yield-point elongation seen in mild steel (the flat part while a Lüders band travels along the specimen). n and K are determined, using the result of §1.3 below, from the tensile strength in Table 1 and the uniform elongation (18 % for structural steel and 8 % for the aluminum alloy; rough values for teaching).

1.3 When does necking begin? — The Considère condition

Work hardening makes the material stronger, but at the same time the cross section gets thinner. How the load F = σ_tA changes is decided by the competition between these two effects. Differentiating the product gives

Differentiating the constant-volume condition AL = constant gives L dA + A dL = 0, that is, dA/A = −dL/L = −dε_t. Putting this in, we get

The load is at its maximum (dF = 0) when

(Considère, 1885). The meaning is as follows. Suppose that one place somewhere in the specimen happens to become slightly thinner. The stress rises there, but that place is also strained extra, so it becomes stronger through work hardening. As long as the work-hardening rate dσ_t/dε_t is greater than σ_t, “getting stronger” wins, the thin place stops deforming, and the deformation spreads over the whole specimen (uniform elongation). When the work-hardening rate falls below σ_t, “getting thinner” wins and the deformation concentrates at that one place. This is necking. As Eq. (6) shows, this coincides with the point where the load turns from increasing to decreasing, that is, with the highest point of the nominal stress–nominal strain curve.

For Hollomon’s equation (4), dσ_t/dε_t = nKε_tⁿ⁻¹ = (n/ε_t)σ_t, so Eq. (7) becomes

(the second is obtained by putting ε_t = n into σ_n = σ_te^−ε_t from Eq. (3)). This is a simple and important result: the true strain at which necking begins is equal to the work-hardening exponent n. Conversely, if the uniform elongation e_u (the nominal strain at the highest point) is known, n = ln(1 + e_u), and if the tensile strength σ_UTS is known, K = σ_UTS/(nⁿe⁻ⁿ). For structural steel (e_u = 18 %, σ_UTS = 400 MPa), n = 0.166 and K = 636 MPa, and the true stress at the onset of necking is 472 MPa, and the work-hardening rate there is indeed 472 MPa as well. For the aluminum alloy (8 %, 310 MPa), n = 0.077 and K = 408 MPa, and both quantities are 334.8 MPa. The Al alloy, with its small work-hardening exponent, necks earlier. Rewriting in terms of the nominal strain gives dσ_t/de = σ_t/(1 + e), so on the graph of σ_t versus e, the point of tangency of the tangent line drawn to the curve from the point e = −1 is the highest point. This is the Considère construction found in many textbooks. In the simulator, if you check “Also show true stress–true strain and the work-hardening rate (Considère condition),” the same point is shown as the intersection of the dσ_t/dε_t curve (green) and the σ_t curve.

1.4 Stress inside the neck — Bridgman’s correction

Once a neck forms, the outer profile at the thinnest cross section is curved, and that cross section is sandwiched from above and below by the thicker parts. The thin part tries to contract sideways, but the thicker material around it holds it back, so inside the cross section tensile stresses (radial and circumferential) arise in addition to the stress along the axis. As a result, the axial stress σ_zz is no longer uniform across the cross section but higher at the center and lower at the surface. Bridgman (an analysis collected in his 1952 book) derived, for the neck of a round bar, with a the radius of the smallest cross section and R the radius of curvature of the neck profile, the following distribution (r is the distance from the center and σ_eq is the flow stress (equivalent stress) of the material).

Let us compute the average stress σ̄ = F/πa², the load F divided by the minimum cross section πa², from this distribution. Substituting u = a² − r² (du = −2r dr), ∫₀^a ln(1 + (a² − r²)/2aR) 2r dr = ∫₀^a2 ln(1 + u/2aR) du = 2aR[(1 + x) ln(1 + x) − x] (with x = a/2R; we used ∫ ln t dt = t ln t − t), so

The coefficient on the right-hand side, B = (1 + 2R/a) ln(1 + a/2R), is at least 1 and returns to 1 as a/R → 0 (no neck). In other words, the F/A_min measured after necking makes the flow stress of the material look higher by the factor B, and dividing it back, σ_eq = σ̄/B, is Bridgman’s correction. For the structural steel model just before fracture (a/R = 0.646), B = 1.146, and in Eq. (9) the value at the center is 1.12 times σ̄ and the value at the surface is 0.87 times σ̄. The fact that the stress is highest at the center and is tensile in all three directions is the reason that ductile fracture starts at the center (see §1.5 below).

1.5 Ductile fracture and brittle fracture

Ductile fracture (metals). Metals contain scattered inclusions such as oxides and sulfides, and second-phase particles. As plastic deformation proceeds, particles crack or the interfaces between particles and the matrix peel apart, and small cavities (voids) form. The voids are stretched and grow with plastic deformation, and they grow faster where the stress is tensile in all three directions. The center of the neck is, as Eq. (9) shows, exactly where those conditions are met, so the voids grow and join at the center and become a central crack perpendicular to the axis. Finally the outer ring slips along a surface at about 45° and tears off, producing the fracture surface called a cup-and-cone. Seen in an electron microscope, the fracture surface is covered with pits left by the voids (dimples). A model that explains the magnitude of ductility from the amount of inclusions and the nucleation and growth of voids was presented by Le Roy et al. (1981). The enlarged (schematic) view in the simulator draws this sequence, but the way the voids grow is shown only as a shape; the growth equation is not solved.

Brittle fracture (ceramics). Dislocations hardly move, so the stress rises along a straight line, and at the instant the stress intensity factor at the tip of a small crack that was there from the start reaches the fracture toughness (§5), the crack crosses the cross section while accelerating. The fracture surface is a flat surface perpendicular to the direction of the load, with no trace of plastic deformation. Taking the values in Table 1 and a surface crack (Y = 1.12), the crack depth a = (K_IC/1.12σ_f)²/π (Eq. (43)) that breaks at the strength σ_f is 25 μm for alumina, 16 μm for Y-TZP, and 57 μm for soda-lime glass. The enlarged view in the simulator draws a crack of this depth.

The energy absorbed before fracture (the area under the nominal stress–nominal strain curve) is, in the simulator’s model, 94 MJ m⁻³ for structural steel and 35 MJ m⁻³ for the aluminum alloy, which are about 590 times and about 215 times that of alumina (0.16 MJ m⁻³), respectively. Although the strength itself is higher for alumina (350 MPa) than the yield stress of steel (250 MPa), what differs decisively is the energy that can be absorbed before fracture.

2. Estimating the Theoretical Strength (Orowan)

First we estimate the stress needed to pull apart atomic planes in a perfect crystal with no cracks or flaws — the theoretical strength. Let the spacing of the atomic planes be a₀ and the displacement of a plane from its equilibrium position be x. As x increases, the stress first rises according to Hooke’s law and, after reaching a maximum σ_th, returns to 0 as the bonds break one after another. We approximate this hump by a half period of a sine function (Orowan’s approximation).

λ/2 is the displacement at which the stress returns to 0 (the bond is completely broken). There are two unknowns, σ_th and λ, so we prepare two conditions.

2.1 Condition 1: Hooke’s law at small displacements

For x ≪ λ, sin θ ≈ θ, so Eq. (11) gives σ ≈ σ_th·2πx/λ. On the other hand, the strain is ε = x/a₀, so Hooke’s law reads σ = Eε = Ex/a₀. Equating the coefficients of x, we get

2.2 Condition 2: The work of pulling apart equals the energy of the two new surfaces

The work done by the stress in pulling an atomic plane of unit area completely apart is the area under the curve of Eq. (11). Pulling apart creates two surfaces, so this is set equal to 2γ_s.

(We used cos π = −1 and cos 0 = 1.) Therefore λ = 2πγ_s/σ_th.

2.3 Eliminating λ

Substituting this λ into the right-hand expression of Eq. (12), we get

Putting in the alumina-equivalent values E = 380 GPa, γ_s = 1.0 J m⁻², and a₀ = 0.2 nm gives

The strength of alumina in Table 1 is 350 MPa, so the theoretical strength is about 125 times the actual strength. Here λ = 2πγ_s/σ_th = 0.144 nm, and the displacement λ/2 at which the bond is completely broken is 0.072 nm, which is 0.36 times a₀ — a reasonable size: the bond breaks when the interatomic distance is stretched by a little under 40 %.

3. Stress Concentration (Inglis)

The first clue to explaining this large gap is that stress concentrates at the tip of a hole or crack. Inglis (1913) obtained the elastic solution for an elliptical hole of major axis 2a and minor axis 2b in an infinitely wide plate under a uniform tensile stress σ applied perpendicular to the major axis, and showed that the stress at the end of the major axis is

For a circular hole (a = b) this is 3σ, the well-known “stress concentration factor of 3.” Near the end (a, 0) of the major axis of the ellipse (x/a)² + (y/b)² = 1, the ellipse can be approximated by the parabola y² ≈ 2(b²/a)(a − x); comparing this with y² ≈ 2ρ(a − x) for a circle of radius ρ gives the radius of curvature ρ = b²/a. From this, b = √(aρ), and therefore 2a/b can be rewritten as 2a/√(aρ) = 2√(a/ρ).

What determines the degree of concentration is not the size of the hole itself but “the ratio of the length to the sharpness of the tip.” For example, if the tip of a crack with a = 10 μm is as sharp as a single atom (ρ = a₀ = 0.2 nm), then √(a/ρ) = √(5×10⁴) = 223.6, and the concentration factor is 1 + 2 × 223.6 ≈ 448. Merely applying σ = 100 MPa from outside brings the stress at the tip to 44.8 GPa — which exceeds the theoretical strength of 43.6 GPa in Eq. (16).

3.1 Combining Orowan and Inglis

Assume that the specimen breaks when the stress at the tip reaches the theoretical strength, that is, σ_max = σ_th, and set ρ = a₀. From the approximate form of Eq. (18) and from Eq. (15), we get

a₀ has dropped out. For a = 10 μm, σ_f ≈ 97 MPa, so the order of magnitude of real strengths comes out. Even the form “the strength is inversely proportional to the square root of the crack size” appears.

3.2 The stress field around an elliptical hole — conformal mapping and complex potentials

Tab ② of the simulator draws not only the “value at the single point at the tip” of Eq. (17) but the stress all around the elliptical hole as a color map. Here is an outline of how that is computed. Consider an infinitely wide plate (the in-plane stresses do not depend on ν, in either plane stress or plane strain) with an elliptical hole of semi-axis a in the x direction and semi-axis b in the y direction, with no force on the edge of the hole (a free surface). Far away, uniform stresses σ_xx = S_x, σ_yy = S_y, σ_xy = 0 are applied. In ②, the load “vertical” is (S_x, S_y) = (0, σ), “horizontal” is (σ, 0), and “equibiaxial” is (σ, σ).

Step 1: Express the stress with two complex functions. It is known that a two-dimensional stress state satisfying the equilibrium equations and the compatibility condition can be written, using analytic functions φ(z) and ψ(z) of z = x + iy, as follows (the Kolosov–Muskhelishvili equations; for the derivation, see textbooks on elasticity theory: Muskhelishvili, 1953).

For example, if φ = Γz and ψ = Γ′z (Γ and Γ′ real), then Φ = Γ, Φ′ = 0, and Ψ = Γ′, so we get the uniform stress σ_xx + σ_yy = 4Γ, σ_yy − σ_xx = 2Γ′, σ_xy = 0. To match the stress far away, it suffices to choose Γ and Γ′ as follows:

With the hole present, φ and ψ approach this form far away and are determined so that the force at the edge of the hole is 0.

Step 2: Map the outside of the ellipse onto the outside of a circle. Imposing conditions on the edge of an ellipse is awkward, so we use a function that maps the outside of the unit circle |ζ| ≥ 1 in the plane of another complex number ζ onto the outside of the elliptical hole.

Let us check. At a point ζ = e^iθ on the unit circle, 1/ζ = e^−iθ, so z = R[(1 + m) cos θ + i(1 − m) sin θ]. Since R(1 + m) = (a + b)/2 × 2a/(a + b) = a and likewise R(1 − m) = b, z = a cos θ + ib sin θ — exactly the edge of the ellipse. If b = a, then m = 0 and the hole is circular; as b → 0, m → 1 and it becomes a crack of length 2a. A function that maps shapes while preserving angles in this way is called a conformal mapping.

Step 3: Determine the functions for which the force on the edge is 0. In the ζ plane, express the functions that take the form of Step 1 far away as powers of ζ and 1/ζ, and determine the coefficients so that the force on the edge is 0 on the unit circle; this gives the following form.

Inglis (1913) solved the same problem using elliptic coordinates. Eq. (23) is a rewriting of that solution in the form of complex potentials, and the simulator inverts the mapping to find ζ from z (of ζ = (z ± √(z² − a² + b²))/2R, the root with the larger |ζ|), and substitutes Φ = φ′(ζ)/ω′(ζ) and Ψ = ψ′(ζ)/ω′(ζ) into Eq. (20) to compute the stress at each point. That the force is 0 at the edge of the hole, that the stress returns to uniform far away, and that the equilibrium and compatibility equations are satisfied have been checked numerically (Appendix A.9).

Step 4: The stress at the edge of the hole. At a free edge the stress normal to the edge σ_nn is 0, so the stress along the edge σ_tt equals the sum σ_xx + σ_yy = σ_nn + σ_tt, which does not depend on the orientation of the coordinates, that is, 4 Re Φ. At ζ = e^iθ, Φ = φ′/ω′ = (Γ − Ae^−2iθ)/(1 − me^−2iθ). To make the denominator real, multiply the numerator and denominator by (1 − me^2iθ); the denominator becomes 1 − 2m cos 2θ + m², and the real part of the numerator is Γ + mA − (A + mΓ) cos 2θ. Using A + mΓ = −Γ′ and Γ + mA = Γ(1 − m²) − mΓ′, we obtain

For the vertical load (Γ = σ/4, Γ′ = σ/2, A = −σ(2 + m)/4), calculating at the tip θ = 0 gives a numerator of 4(1 − m)(Γ − A) = (1 − m)(3 + m)σ and a denominator of (1 − m)², so σ_tt = σ(3 + m)/(1 − m). Putting in 1 − m = 2b/(a + b) and 3 + m = (4a + 2b)/(a + b) gives σ(4a + 2b)/2b = σ(1 + 2a/b) — Eq. (17) has come out. At the top and bottom edges θ = π/2 (cos 2θ = −1), the numerator is −(1 + m)²σ and the denominator is (1 + m)², so σ_tt = −σ (compression).

3.3 It depends on the direction of the load — comparison with Kirsch’s solution for a circular hole

If the direction of the load is changed in Eq. (24), the place on the edge where the stress concentrates changes accordingly. The results are summarized in the table (values divided by σ).

Even if a long thin hole like a crack is pulled in the direction along the crack (horizontal), the stress hardly concentrates (1 + 2b/a → 1 at the top and bottom edges). Only the stress in the direction that opens the crack faces is dangerous. Under equibiaxial loading, the concentration at the tip, 2a/b, is almost as large as for the vertical load.

For a circular hole, the solution for a uniform tension σ in the x direction, in polar coordinates (r, θ) with the center of the hole as origin (θ is measured from the x axis, the direction of the load), was found in closed form by Kirsch (1898).

At the edge r = a, σ_rr = σ_rθ = 0 (free surface), and σ_θθ = σ(1 − 2 cos 2θ). It is 3σ at θ = 90°, perpendicular to the load, and −σ at θ = 0, in the direction of the load. As r → ∞, σ_rr → (σ/2)(1 + cos 2θ) = σ cos²θ, returning to the uniform tension. This also agrees with the value obtained by setting m = 0 in Eq. (24). Under a uniform compression of the opposite sign, a tension of +σ appears at the edge in the direction of the load. This is why small cracks can grow from the edge of a hole or flaw even when only compression is applied, and it corresponds to the explanation of Misconception 1 in §1.

For a long thin ellipse close to a crack (b ≪ a), the stress ahead of the tip also takes a simple form. In the limit b → 0 (a crack of length 2a), the stress at a point a distance r from the tip on the extension of the crack is σ_yy = σ(a + r)/√((a + r)² − a²), and for r ≪ a, (a + r)² − a² = 2ar + r² ≈ 2ar, so σ_yy ≈ σ√(a/2r) = σ√(πa)/√(2πr) — this is K_I/√(2πr) in Eqs. (37) and (38) of §5 (derived in §5.5). For example, at the point b/a = 0.005, r = 0.01a, the solution for the elliptical hole gives 7.128σ, the exact solution for the crack 7.124σ, and the K-field approximation 7.071σ, which nearly agree. That the region of high σ₁ at the tip looks like a kidney shape bulging above and below the crack line in the color map can be explained by Eq. (47) in §5.5.

3.4 The direction of the principal stress and the direction in which the crack advances

When the stress at a point (σ_xx, σ_yy, σ_xy) is viewed on a plane whose normal is rotated by an angle ϑ from the x axis, the stress normal to that plane is

(this follows from writing the force balance on a small triangular element). Differentiating with respect to ϑ and setting the result to 0 gives tan 2ϑ_p = 2σ_xy/(σ_xx − σ_yy), and putting the resulting cos 2ϑ and sin 2ϑ back gives the maximum and minimum values. These are the principal stresses.

The default display of simulator tab ② is the maximum principal stress σ₁, and the short gray lines show its direction (ϑ_p). A crack in a brittle material advances perpendicular to the direction of the maximum principal stress near its tip. Erdogan and Sih (1963) compared the way an inclined crack turns with experiments under the condition that the crack extends in the direction in which the circumferential stress around the tip is maximum (the maximum circumferential stress criterion), and Gol’dstein and Salganik (1974) presented the principle of local symmetry, that a crack advances in the direction in which the field at the tip has no shear component (mode II). Both are expressions of the same intuition that “a crack advances perpendicular to the tension that opens its tip.” The green arrow in ② is drawn in this direction (the outward normal of the edge) from the point on the edge where σ₁ is maximum. For the vertical load the crack extends in the x direction from the tip, and for the horizontal load it extends in the y direction from the top and bottom edges. With equibiaxial loading on a circular hole the whole edge has the same 2σ, so where the hole breaks is decided by the small flaws present there, and no arrow is drawn.

Right next to a free edge, the principal stress directions are always along the edge and perpendicular to the edge. This is because at the edge both the normal stress to the surface and the shear stress are 0, so the stress along the edge σ_tt and the zero normal stress are the principal stresses. If σ_tt > 0, σ₁ runs along the edge, and a crack advancing perpendicular to it meets the edge at a right angle. If you check “Overlay the direction of σ₁ as short lines” in ②, you can confirm that the lines near the edge of the hole line up along the edge.

3.5 Discussion: why do glass fragments seem to have so many triangles?

Fact 1: A crack advances perpendicular to the local σ₁ (§3.4). The shapes of the fragments therefore trace the pattern of principal-stress directions in the plate at the moment of fracture.

Fact 2: A crack that arrives later meets an earlier crack at almost a right angle. The surface of an earlier crack is a free surface, so near it σ₁ runs along the surface, as in §3.4, and a later crack turns to meet the surface perpendicularly. When cracks form one at a time, in sequence and divide a region, as in the cracking of a ceramic glaze or of drying mud, T-shaped junctions are common, and it has been reported that the number of sides of the resulting cells averages 4 (Bohn, Douady, Couder, 2005; Bohn, Pauchard, Couder, 2005). This means that “if the region is simply divided in sequence, the fragments tend to be quadrilaterals rather than triangles.”

Fact 3: A fast-running crack branches. When a material such as glass breaks while still holding a large amount of stored elastic energy, the crack accelerates and splits into two, and the branches split further. Yoffe (1951) showed that at the tip of a fast-running crack the direction of maximum circumferential stress deviates obliquely from the direction of travel, and this is regarded as one of the classical explanations of branching (the cause of branching is still under study). The fractography literature summarizes that the branching angle varies with the stress state, being about 30–45° for uniaxial tension (with large scatter) and larger for biaxial tension (Quinn, 2020). It is also well known that when a point impact or bending is applied to a plate, cracks that extend radially from the load point (radial cracks) and circumferential cracks that cross them form (Quinn, Swab, Patel, 2021; Bradt, 2011).

What the facts seem to imply. When a third crack crosses the narrow wedge between two cracks that emerge from one point at a sharp angle (the two branches of a branching crack, or adjacent radial cracks), a roughly triangular fragment is formed. That the branching angle is sharp, about 30–45°, is consistent with triangular fragments with a sharp corner at the tip of the wedge being easily formed. On the other hand, if the T-shaped junctions of Fact 2 continue, the fragments approach quadrilaterals. Therefore, if it does feel as though there are “many triangles,” the natural interpretation is that many branches and radial cracks were formed in a fast, high-energy fracture. However, whether triangles really are in the majority, and by how much in which kinds of breakage, cannot be said without counting the fragment shapes. Please also note that the elliptical-hole figure in simulator tab ② is a static elastic solution and does not represent fast-running cracks or branching themselves. A question for students: In ②, change the direction of the load, predict from the pattern of σ₁ lines “where the crack will turn and in which direction,” and compare your prediction with photographs of actually broken glass.

4. Griffith’s Energy Balance [Core]

Griffith’s (1921) idea was not to ask what happens at the tip of the crack but to look only at the energy bookkeeping when the crack extends a little. When a crack extends, (i) new surfaces are created, which costs energy, and (ii) the material around the crack relaxes, so elastic strain energy is released. If (ii) exceeds (i), the crack keeps extending without anything being done from outside.

We set the situation as follows. In an infinitely wide plate of thickness 1 (unit thickness) there is a through-crack of full length 2a, and a uniform tensile stress σ is applied far away, perpendicular to the crack. Plane stress is assumed. We find (i) and (ii) as functions of a for the case in which the crack is introduced, and add them.

4.1 The shape the crack opens into

According to elasticity theory (the limit b → 0 of the Inglis elliptical hole in §3), this crack opens into an elliptical shape. Taking the center of the crack as the origin and the x axis along the crack, the displacement of the upper crack face is

(the lower face moves by −v(x)). This result is borrowed from elasticity theory. The opening is greatest at the center, 4σa/E, and even for a = 10 μm and σ = 156 MPa it is only 16 nm.

4.2 The elastic energy released when a crack is introduced

Consider the following procedure. Starting from a plate without a crack that carries σ, we cut along the surface where the crack is to be introduced, but on the cut surface we leave the stress σ that was originally acting there as a “holding force.” While it is held, nothing changes. If this holding force is slowly reduced from σ to 0, the crack faces open from 0 to v(x). The work that the material does against the holding force at this time is the energy released by introducing the crack.

In a linear elastic body the holding force and the opening change in proportion, so the work is “initial force × final displacement × 1/2” (the triangle under the force–displacement line). Adding up the upper and lower faces, we get

The last integral is the area under y = √(a² − x²) (the upper half circle of radius a), that is, the area of a half circle, πa²/2. Therefore

4.3 The energy of the newly created surfaces

The crack length is 2a, the thickness is 1, and there are two faces, upper and lower.

4.4 Finding the maximum of the total energy

The change in the energy of the system relative to the state without a crack is obtained by subtracting the released part of Eq. (30) and adding Eq. (31):

The first term decreases as a², and the second increases as a. For small a the second term wins, and for large a the first wins, so U(a) has the shape of a hill. Differentiating with respect to a, we get

As long as dU/da > 0, the energy increases when the crack extends, so it does not extend. Once dU/da < 0, the energy decreases the more it extends, so it does not stop. The boundary dU/da = 0 is

Reading this as the stress (strength) at which a crack of a given a breaks, or as the half-length of the crack that starts to extend under a given σ, gives two forms.

The second derivative is d²U/da² = −2πσ²/E < 0, so this point is a maximum. The height of the hill, using πσ²a_c/E = 2γ_s from the left-hand expression of Eq. (34), is U(a_c) = −2γ_sa_c + 4γ_sa_c = 2γ_sa_c. With the alumina-equivalent values and a = 10 μm, we get

This is 1/280 of the theoretical strength of 43.6 GPa. For plane strain (E → E/(1 − ν²), ν = 0.22) it becomes 159 MPa. Conversely, to explain the strength of alumina, 350 MPa, with γ_s = 1.0 J m⁻² unchanged, it suffices to have a crack with a = 2.0 μm.

5. Stress Intensity Factor and Energy Release Rate

Griffith looked at the energy of the whole system. Irwin (1957) moved the viewpoint to the crack tip and expressed the “strength” of the stress field near the tip with a single number. They are two faces of the same fracture condition.

5.1 The stress field at the tip and the stress intensity factor K

In a linear elastic body, at distances r from the tip that are sufficiently smaller than a, the tensile stress on the extension of the crack plane is

It diverges like 1/√r as r → 0, but the form of the divergence is the same regardless of the crack length or the shape of the specimen, and all of the differences go into the coefficient K_I. This is called the stress intensity factor (the subscript I denotes mode I, the opening mode, in which the crack faces are pulled apart). For a wide plate,

Y is the geometry factor. A crack with a straight front that enters from the edge of a plate (a free surface) opens more easily because there is no constraint on one side, so Y is 1.12, greater than 1. For a small semicircular crack in the surface, Y runs along the crack front from about 0.66 (at the deepest point) to 0.73 (where it meets the surface), and the representative value often used is Y ≈ 0.71, the coefficient 2/π of a penny-shaped crack multiplied by the surface correction 1.12. In any case it is smaller than 1. On this page and in the simulator, for simplicity, a surface crack is represented by Y = 1.12 (a conservative treatment that overestimates K_I). The unit of K_I is MPa m^1/2.

5.2 The energy release rate G

The elastic energy released when the crack advances by a unit area is defined as the energy release rate G. A crack of total length 2a has two tips, and when a increases by da each end advances by da, so the crack area (for unit thickness) increases by 2da. Using Eq. (30), we get

Squaring Eq. (38) for Y = 1 gives K_I² = πσ²a, so this can be written directly as

The energy of the whole system (G) and the stress field at the tip (K) are now connected by a single equation.

5.3 The two fracture conditions are the same

The surface-energy side is counted by the same rule. When the crack area increases by a unit amount, the surface increases by two sheets, upper and lower, so the energy required is 2γ_s. Writing the condition “fracture occurs when G reaches 2γ_s,” we get

This agrees with Eq. (35). In the language of K, the condition is “fracture occurs when K_I reaches the critical value K_IC.” If we replace G in Eq. (40) by the critical value G_c, the relations between the two are as follows.

G_c is called the fracture energy (resistance to crack growth), and γ_eff = G_c/2 the effective surface energy (some references call the γ_eff per fracture surface on one side the fracture energy γ_f, so watch out for the factor of 2 when comparing values). K_IC is a material value that can be measured directly in a test, and this is the fracture toughness. For an ideally brittle body γ_eff = γ_s, and with alumina-equivalent values we would expect K_IC = √(2Eγ_s) = 0.87 MPa m^1/2. The 3.5 MPa m^1/2 in Table 1 is 4.0 times that, or 16 times in terms of G_c. The origin of this gap is the subject of §6.

5.4 The tolerable crack size

Solving the left-hand equation of Eq. (42) for a gives the largest crack that can be tolerated under a stress σ.

For a design stress σ = 300 MPa and a surface crack (Y = 1.12), this is 34.5 μm for alumina, 1.6 μm for soda-lime glass, and 180 μm for Y-TZP. Because a_c is proportional to the square of K_IC, a 2.3 times larger K_IC (alumina → Y-TZP) allows a crack 5.2 times larger. The 1.6 μm for glass is a size that even a light scratch on the surface exceeds. Compared with the estimate seen at the end of §4, that “for γ_s = 1.0 J m⁻², a crack of 2.0 μm is enough to give 350 MPa” (a through-crack), the crack size at which real alumina (K_IC = 3.5 MPa m^1/2) fractures at 350 MPa is 32 μm for the same through-crack (16 times, the same as the ratio of G_c), and 25 μm for a surface crack (Y = 1.12).

5.5 The whole stress field around the crack tip — Westergaard’s solution and the K field

Eq. (37) held only on the extension of the crack plane. The upper figure of simulator tab ③ colors σ_yy over the whole region around the tip. We derive that field by the method of Westergaard (1939). For a wide plate containing a crack of length 2a (y = 0, |x| < a) under a remote stress σ, we use the complex function

With this function, √(z² − a²) is purely imaginary on the crack faces (y = 0, |x| < a), so that σ_yy = σ_xy = 0, and Z → σ far away (this Z is a biaxial solution that also makes σ_xx equal to σ far away, so for uniaxial tension σ is subtracted from σ_xx; σ_yy and σ_xy do not change). On the extension of the crack (y = 0, x = a + r), σ_yy = σ(a + r)/√((a + r)² − a²), which is the same as the exact solution for a crack in §3.3.

Near the tip we use polar coordinates centered on the tip, z = a + re^iθ (r ≪ a). z² − a² = (z − a)(z + a) ≈ re^iθ × 2a, so we get

(We put K_I = σ√(πa); Z′ is the derivative of Z ∝ (z − a)^−1/2.) Re Z = (K_I/√(2πr)) cos(θ/2), y Im Z′ = r sin θ × (K_I/2√(2π))r^−3/2 sin(3θ/2), and putting these together with sin θ = 2 sin(θ/2) cos(θ/2), we obtain the mode I K field.

For θ = 0 this returns to Eq. (37). The angular functions depend on neither a nor σ, which enter only through K_I — this is what is meant by “K expresses the strength of the field at the tip with a single number.” Substituting into Eq. (27), we get (σ_xx + σ_yy)/2 = (K_I/√(2πr)) cos(θ/2), ((σ_xx − σ_yy)/2)² + σ_xy² = (K_I²/2πr) cos²(θ/2) sin²(θ/2)[sin²(3θ/2) + cos²(3θ/2)], so

Putting s = sin(θ/2), the angular part becomes √(1 − s²)(1 + s); setting its derivative with respect to s equal to 0 gives s = 1/2, that is, θ = ±60°, where the maximum value is 1.30. At the same distance, σ₁ is 30 % higher in the direction 60° off to the side than straight ahead of the crack (1.00), so in the color map of σ₁ in ② the high-stress region has a kidney shape bulging above and below. The σ_yy painted in the upper figure of ③ (the component that opens the crack) also has a shape bulging above and below the line of the crack, since the angular part of Eq. (46), cos(θ/2)[1 + sin(θ/2)sin(3θ/2)], is 1.30 at θ = ±60° (1.3 times its value straight ahead). The differences from σ₁ are that σ_yy is only the component in the direction that opens the crack, and that it decreases faster from the side of the tip toward the region behind it (at θ = 90°, σ_yy is 1.06 and σ₁ is 1.21; at θ = 120°, σ_yy is 0.50 and σ₁ is 0.93).

5.6 The size of the plastic zone — Irwin’s estimate

Eq. (46) goes to infinity at the tip, but in a metal plastic deformation occurs where the stress exceeds the yield stress σ_y. Finding the distance along the extension at which K_I/√(2πr) = σ_y, we get

However, the part r < r_p cannot support a stress larger than σ_y, so the corresponding force is pushed outward, and the plastic zone spreads to about twice that size, 2r_p. Irwin corrected for this by regarding the crack as “behaving as if it were longer by r_p” (the derivation of this correction is summarized in Anderson, 2017). In simulator tab ③ this plastic zone is shown by a dashed circle of diameter 2r_p centered a distance r_p ahead of the tip. In plane strain (the interior of a thick plate), triaxial constraint makes yielding more difficult, and the plastic zone becomes about 1/3 as large.

Calculating r_p for K_I = K_IC with the values of Table 1 gives 6.4 mm for structural steel and 2.0 mm for the aluminum alloy. In ceramics dislocations hardly move at room temperature, so a yield stress is hard to define, but if we use the rough relation between the hardness H and the yield stress, σ_y ≈ H/3 (Tabor, 1951; §12.4; since a hydrostatic pressure acts under an indentation, this is a nominal value estimated from the plastic deformation that occurs there), we get 69 nm for alumina (H ≈ 16 GPa, σ_y ≈ 5.3 GPa), 27 nm for soda-lime glass, and even for Y-TZP only 0.6 μm — at most about the size of a grain, orders of magnitude smaller than the crack size (tens of μm), that is, effectively 0. Note that H ≈ 3σ_y holds for metals with a large E/σ_y; in ceramics H/σ_y is smaller, about 1.5–2 (Marsh, 1964). H/3 underestimates σ_y, so the r_p here is a rough upper bound (the actual value is about 1/2 to 1/4 of it). The tip of a ceramic has almost no room to be rounded by plasticity.

Substituting K_I = Yσ√(πa) into Eq. (48) gives r_p/a = Y²σ²/(2σ_y²), which becomes as large as 0.63 at σ = σ_y. Fracture can be discussed in terms of K only when the plastic zone is sufficiently small compared with the dimensions of the crack and the specimen (small-scale yielding); for metals at stresses close to the yield stress it remains a rough guide (Misconception 5 in §6). With a plastic zone the tip opens and is blunted, and in the strip-yield model the crack-tip opening displacement is δ ≈ K_I²/(Eσ_y) (Anderson, 2017). For structural steel at K_I = K_IC, δ ≈ 48 μm, and as seen in Problem 2 of §15, blunting of the tip greatly lowers the stress concentration.

6. Fracture Toughness Values and Toughening Mechanisms

We convert the K_IC of Table 1 into the fracture energy G_c and the effective surface energy γ_eff with Eqs. (40) and (42).

The parentheses around the a_c of the metals are there because 300 MPa exceeds the yield stress in Table 1 (250 MPa for steel, 270 MPa for the Al alloy), so the metal would yield before the crack extends. Repeating the same calculation at 200 MPa, below the yield stress, gives 1.6 cm for steel and 5.7 mm for the Al alloy, so at stresses below the yield stress, the crack that a metal can tolerate is millimeters or larger (although, as in §5.6, the plastic zone is large, so this is only a rough guide).

The γ_eff = 4.0 J m⁻² of glass is only a few times the surface energy (around 1 J m⁻²). Griffith was able to confirm his theory with glass because glass is close to an ideally brittle body. For polycrystalline ceramics, however, the values are 16.1 for alumina, 58.1 for silicon nitride, and 152 J m⁻² for Y-TZP, one to two orders of magnitude larger than the surface energy. The G_c of metals is larger still, about 370 times that of alumina for steel and about 400 times for the Al alloy. This means that a large amount of work other than creating new surface is spent in extending the crack.

6.1 Where toughness comes from — shielding at the tip

When actions that tend to close the crack operate around or behind the crack tip, the stress intensity factor that the tip actually feels is smaller than the externally applied K_I. This is called shielding; writing the shielding part as K_sh, we have

The material at the tip itself fractures when K_tip reaches its intrinsic value K₀ (the toughness needed just to break the bonds), so the fracture toughness seen from the outside is K₀ + K_sh. Because it is hard to strengthen the bonds at the tip, the toughening of ceramics aims to increase K_sh. There are four representative methods (transformation toughening, crack deflection, bridging, and microcracking).

6.2 Transformation toughening — zirconia

When cooled from the high-temperature tetragonal form (t phase), ZrO₂ transforms to the monoclinic form (m phase), and its volume increases by about 4 % (the literature gives a range of 3–5 %). A sintered body of pure ZrO₂ cracks itself during cooling because of this expansion. If a small amount of Y₂O₃ or the like is dissolved in it and the grains are kept small, the t phase can be retained even at room temperature in a metastable state. This is Y-TZP (yttria-stabilized tetragonal zirconia polycrystal).

When a crack approaches such a material, the large tensile stress at the tip triggers the t → m transformation in a band-shaped region around the crack. The transformed grains try to expand but are held in by the surrounding untransformed material, so a compressive stress that tends to close the crack faces arises, and this is the K_sh of Eq. (49). Garvie et al. (1975) showed that a transformation of the same type as the martensitic transformation that occurs in the quenching of steel can be used for toughening, and titled their paper “Ceramic steel?” The mechanism of transformation toughening and the design of zirconia-based materials that exploit it are summarized in the review by Hannink et al. (2000). The K_IC = 8.0 MPa m^1/2 of Y-TZP in Table 1 is 2.3 times that of alumina (the K_IC of zirconia depends strongly on the measurement method, and with the standardized notched-beam method values of about 4–6 MPa m^1/2 may be obtained).

6.3 Crack deflection and bridging

Besides transformation toughening, there are several other ways to create shielding through the microstructure (the review by Evans, 1990, is detailed).

• Deflection. When a crack bends in a zigzag along weak grain boundaries or a second phase, the opening component at the tip decreases and the area of the fracture surface also increases.
• Bridging. When elongated grains, coarse grains, or fibers remain connected behind the crack without breaking, they exert a force that pulls the crack faces together. The high K_IC of silicon nitride is due to the large bridging by the columnar grains that grow during sintering.
• Microcracking. Many microcracks form around the tip, dissipating energy and relaxing the stress at the tip.

6.4 The R-curve — a crack that becomes harder to extend the more it extends

The bridging and transformation zones accumulate behind the crack as it extends, so the shielding grows larger the longer the crack becomes. The curve along which the resistance to crack growth R (with the same dimension as G_c) increases with the crack extension Δa is called the R-curve. Since G = πσ²a/E of Eq. (39) is a straight line proportional to a, the condition for fracture to become unstable is

If R is constant, dR/da = 0, so the crack becomes unstable the moment G reaches R (the sudden fracture of §4). In a material in which R rises, even when G = R the crack extends stably by a small amount and stops, and it becomes unstable only when the stress is raised further and the straight line of G touches the R-curve. As a result, there is the advantage that the strength becomes insensitive to the size of the initial crack and the scatter tends to be small. On the other hand, for small cracks the shielding has not yet developed, so they do not enjoy as much benefit as the steady-state K_IC suggests.

6.5 Emitting dislocations or breaking bonds — the Rice–Thomson picture

The plastic zone of §5.6 was a continuum view. Looking at the same problem on the atomic scale, it becomes a competition over which happens first at a crack tip that is sharp down to a single atom. (i) The bond at the tip breaks and the crack advances by one atom (cleavage). (ii) A dislocation is emitted from the tip and escapes along a slip plane. If (ii) occurs first, each time a dislocation is emitted the tip is displaced by one step and becomes rounded (blunting), the stress concentration drops, and the material behaves in a ductile manner. If (i) comes first, the tip advances while remaining sharp and the material fractures in a brittle manner. This view was put forward by Rice and Thomson (1974), and Rice (1992) quantified it using the Peierls model.

The condition for bond breaking is the same as in §5.3, G = 2γ_s. The condition for dislocation emission is determined by the height of the barrier that has to be crossed when the material above and below a slip plane is displaced by one atom — the unstable stacking energy γ_us — and, for a slip plane inclined at an angle θ to the crack plane (mode I), it is given by

(At θ = 60°, 1 + cos θ = 3/2 and sin² θ = 3/4, so the denominator is 9/8.) If G_e < 2γ_s, dislocation emission comes first (ductile); if G_e > 2γ_s, bond breaking comes first (brittle). To convert to a stress intensity factor, use K = √(GE/(1 − ν²)) for plane strain. The lower figure of simulator tab ③ (the schematic of atoms) and its “Rice–Thomson scale” compare these two values with the G at the tip.

The γ_s of the metals is of the same order as the value that Tyson and Miller (1977) estimated from the surface tension of liquids. For steel, the difference between 2γ_s = 4.8 and G_e = 3.6 is small — close enough that the order could be reversed depending on the conditions. In fact, the γ_us = 0.50 J m⁻² in Table 4 is on the low side of values calculated with interatomic potentials (EAM); if we use a first-principles (DFT) value of about 1 J m⁻² for {110}⟨111⟩ in bcc iron, we get G_e ≈ 7 J m⁻² > 4.8, and this simple criterion classifies bcc iron itself on the brittle side. That the verdict can switch within the uncertainty of the estimate of γ_us corresponds to the fact that the same steel can be either ductile or brittle depending on temperature, constraint, and impurities. Some studies also link the embrittlement of bcc iron at low temperatures (the ductile–brittle transition) to this closeness. In fcc metals such as aluminum the difference is large, and they remain ductile even at low temperatures.

6.6 Why it is hard for dislocations to move in ceramics

We look at the reasons stated in §1 a little more quantitatively.

• The Peierls stress is large. When a dislocation moves by one atomic spacing, the atoms in the core cross a barrier while rearranging their bonds. The shear stress needed to cross that barrier is the Peierls stress, and in the Peierls–Nabarro model, with the spacing of the slip planes d, the magnitude of the Burgers vector b, and the shear modulus μ, it is estimated as τ_P ≈ [2μ/(1 − ν)] exp[−2πw/b], with the core width w = d/(1 − ν) (Hull and Bacon, 2011). The key point is that it depends exponentially on the core width. In metallic bonding, which has no directionality, the core is wide and τ_P is small, whereas in strongly directional covalent bonding the core is narrow and τ_P is larger by orders of magnitude (the numerical values from this formula are only rough estimates).
• In ionic crystals, the directions in which slip can occur are limited. Only planes and directions in which ions of the same sign do not face each other during slip can be used. The main slip system of rock-salt-type crystals is {110}⟨110⟩.
• There are not enough independent slip systems. For the grains of a polycrystal to deform into arbitrary shapes without gaps between them, five independent slip systems are needed (the von Mises condition). The rock-salt-type {110}⟨110⟩ has only two independent ones (Groves and Kelly, 1963). Even if single-crystal MgO slips a little at room temperature, in a polycrystal the deformation does not match at the grain boundaries and cracks form.
• The Burgers vector is long. In crystals with large unit cells such as Al₂O₃, b is long and the energy of a dislocation (∝ μb²) is large.
• Glass has no crystalline dislocations. The arrangement of atoms has no periodicity, so the very notion of a slip plane does not apply.

All of these relax as the temperature rises, and at high temperatures ceramics also show plastic deformation and creep. The fact that dislocations “cannot move” at room temperature is the background of the large G_e in Eq. (51) and of the effectively zero r_p of §5.6.

6.7 The shape of the transformation zone and the shielding

The t → m transformation of §6.2 is one in which the volume increases, so it occurs more easily where the tensile part of the hydrostatic stress (the mean stress σ_m) is large. In plane strain σ_zz = ν(σ_xx + σ_yy), so from Eq. (46) σ_m = (σ_xx + σ_yy + σ_zz)/3 = (1 + ν)(σ_xx + σ_yy)/3. Substituting σ_xx + σ_yy = 2(K_I/√(2πr)) cos(θ/2) and setting it equal to the critical value σ_m^c at which the transformation occurs, the boundary of the transformation zone is

This gives a shape bulging out in front of the tip (McMeeking and Evans, 1982). What is interesting is that a transformation zone ahead of the tip alone does not give shielding; McMeeking and Evans showed that shielding that closes the tip arises only when the crack has extended and the transformation zone remains behind the crack as a band. With the half-width of the band h, the volume fraction transformed V_f, and the volumetric strain of the transformation e^T (about 0.04), the steady-state shielding takes the form ΔK ∝ EV_fe^T√h/(1 − ν), with a coefficient found to be about 0.2 (McMeeking and Evans, 1982; Budiansky et al., 1983). The more the band extends, the more the shielding develops, so Y-TZP shows the R-curve of §6.4. In simulator tab ③, in the atomic picture for Y-TZP, the atoms in a region of the shape of Eq. (52) are colored orange and are left as a band after the crack has advanced (the size is schematic).

7. Deriving the Weibull Distribution from the Weakest Link [Core]

In §4–§6 we saw that the strength is set by the size of the most dangerous crack. But the sizes of the cracks inside specimens are not uniform, even among specimens made by firing the same powder under the same conditions. That is why strength scatters. Moreover, a larger specimen is more likely to contain a larger crack, so the strength decreases on average. Weibull statistics treats these two effects quantitatively.

7.1 The weakest-link idea

A chain breaks as a whole when its weakest link breaks. If each of n links independently withstands a stress σ with probability p_s(σ), the probability that the whole chain withstands it is the product of the probabilities of the independent events.

If we divide a specimen of volume V into V/V₁ small pieces of volume V₁ and make the same statement, P_s(V) = [P_s(V₁)]^V/V₁ holds. Taking the logarithm, ln P_s(V) = (V/V₁) ln P_s(V₁); that is, ln P_s is proportional to the volume. So if we define the “risk” per unit volume, φ(σ) = −ln P_s(V₁)/V₁ (≥ 0), we obtain

P_s is the probability that the specimen survives (does not break) when a stress σ is applied. Everything up to here follows only from the weakest-link assumption that “the small pieces are independent, and if any one of them breaks, the whole breaks.” The form of φ(σ) is not determined by this assumption. Weibull (1939) proposed the following form by taking φ to be a power function, and in a 1951 paper showed that it fits the data of many materials.

m is the Weibull modulus (a measure of how small the scatter is), and σ₀ is the characteristic strength of a specimen of volume V₀. For V = V₀, P_f = 1 − 1/e = 63.2 % at σ = σ₀. The power function was chosen as a simple form that fits experiments well, and Weibull himself did not derive it from theory. Next we show that this form can be derived from the distribution of flaw sizes.

7.2 Deriving it from the distribution of flaws — the model of the simulator

We assume the following four things. Simulator tab ④ creates and breaks specimens exactly according to these assumptions.

1. Number of flaws. The number n of flaws contained in a specimen of volume V follows a Poisson distribution with mean λV (the flaws are scattered independently of one another).
2. Size of flaws. The size a of each flaw is independent of the others and, above a_min, follows a Pareto law (a power law). The larger k is, the less likely large flaws are to occur.
3. Strength of each flaw. A flaw of size a_i becomes the origin of fracture at the stress σ_i = K_IC/(Y√(πa_i)) at which K_I = K_IC in Eq. (42).
4. Weakest link. The strength of the specimen is the smallest of the σ_i of all the flaws.

Step 1: flaws that are dangerous at a stress σ. By assumption 3, σ_i ≤ σ holds when the flaw is at least as large as the following (the same form as Eq. (43)).

Step 2: the probability that a single flaw is dangerous. Put a = a^* into the right-hand equation of Eq. (56). Since a_min/a^* = πa_minY²σ²/K_IC², we get

σ^* is the stress at which the smallest flaw a_min breaks. Eq. (58) holds for σ ≤ σ^*, and for σ > σ^* all flaws are dangerous (q = 1).

Step 3: the probability that no flaw is dangerous. When there are n flaws, the probability that none of them is dangerous is (1 − q)ⁿ (a product, since the sizes are independent). Averaging this over the Poisson distribution of Eq. (56) and using the series of the exponential function Σxⁿ/n! = e^x, we get

Step 4: substituting q. Inserting Eq. (58) gives

Step 5: comparison with the Weibull form. Putting V = V₀ (reading σ₀ as the characteristic strength of a specimen of this volume), Eq. (55) can be written P_s = exp[−(σ/σ₀)^m]. For λV(σ/σ^*)^2k = (σ/σ₀)^m to hold for all σ, the exponents of σ must match, m = 2k, and the coefficients must match, λV(σ^*)^−m = σ₀^−m.

In the range σ ≤ σ^*, the Weibull distribution has come out exactly from the distribution of flaws. Moreover, since σ₀ ∝ V^−1/m, the V/V₀ form of Eq. (55) is also reproduced as it is. With the default values of the simulator (K_IC = 3.5 MPa m^1/2, Y = 1.12, a_min = 5 μm, k = 5, λV = 300), we get

The probability that a specimen whose strength exceeds σ^* appears is e^−λV = e⁻³⁰⁰, which is effectively 0.

7.3 The mean and the coefficient of variation

The probability density of the strength is p(σ) = −dP_s/dσ. Integrating the mean by parts, the boundary term [−σP_s]₀^∞ is 0, so we can write ⟨σ⟩ = ∫₀^∞ P_s dσ (the e^−λV of the part σ > σ^* is regarded as 0). With the substitution t = (σ/σ₀)^m, that is, σ = σ₀t^1/m and dσ = (σ₀/m)t^1/m−1dt, we get

(Here ∫₀^∞ t^z−1e^−tdt is called the gamma function Γ(z). Integration by parts shows that Γ(z + 1) = zΓ(z), and for integers Γ(n + 1) = n!. We used this relation in the last step.) In the same way, integration by parts gives ⟨σ²⟩ = ∫₀^∞ 2σP_s dσ, and with the substitution 2σ dσ = σ₀²(2/m)t^2/m−1dt, we get

The coefficient of variation (standard deviation ÷ mean) is therefore

The approximation is for large m. Using ln Γ(1 + x) = −γ_Ex + (π²/12)x² + ⋯ (γ_E is Euler’s constant), we have ln[Γ(1 + 2/m)/Γ(1 + 1/m)²] ≈ (π²/12)(4 − 2)/m² = π²/(6m²), in which the first-order terms cancel each other, and CoV² = Γ(1 + 2/m)/Γ(1 + 1/m)² − 1 ≈ π²/(6m²).

For the default values (m = 10, σ₀ = 446 MPa), Γ(1.1) = 0.95135 and Γ(1.2) = 0.91817, so the mean is ⟨σ⟩ = 446 × 0.951 ≈ 424 MPa and CoV = √(0.91817 − 0.95135²)/0.95135 ≈ 0.120 (0.128 from the approximate formula). For m = 20 the CoV falls to 0.062.

8. The Weibull Plot and Estimating m

8.1 Making it a straight line

If we fix the specimen size at V = V₀, Eq. (55) becomes 1 − P_f = exp[−(σ/σ₀)^m]. Taking the logarithm of the reciprocal gives ln[1/(1 − P_f)] = (σ/σ₀)^m, and taking the logarithm once more gives

If we take x = ln σ on the horizontal axis and y = ln ln[1/(1 − P_f)] on the vertical axis, we get a straight line of slope m whose x at y = 0 is ln σ₀. This is the Weibull plot.

8.2 From the data of N test specimens

Arrange the strengths of N specimens in ascending order, σ₁ ≤ σ₂ ≤ ⋯ ≤ σ_N, and estimate the fracture probability of the i-th specimen by

If we used i/N, the strongest specimen would have F = 1 and y would become infinite, so we shift it by a half. If we fit a straight line to the points (x_i, y_i) = (ln σ_i, ln ln[1/(1 − F_i)]) by least squares, then the slope, together with the fact that the line passes through the mean point (x̄, ȳ) (ȳ = m x̄ − m ln σ₀), gives

Fig. 3 shows the result of breaking N = 100 specimens with the model of §7.2, using the same default values as the simulator and the pseudorandom numbers of seed 20261119. Eq. (68) gives m_est = 9.34 and σ_0,est = 447 MPa, values close to the theoretical ones (m = 10, σ₀ = 446 MPa). Of the 100 specimens, the weakest was 289 MPa and the strongest was 510 MPa.

8.3 The number of tests and the scatter of the estimate

m_est scatters from one set of test data to another. Repeating the test 2000 times with the same model, changing only the seed, gave the following.

We can see that σ₀ is reasonably well determined even from a small number of specimens, whereas m is hard to determine. Even if a test of 10 specimens suggests “a material with m = 15,” the true value may be 10. For N = 10, although the median is 9.9, the distribution has a long tail on the large side, so the mean is 10.5; on average, m is overestimated. Standards for the strength data of ceramics (JIS R 1625, ASTM C1239, and others) specify a method that estimates by maximum likelihood instead of least squares and corrects the bias due to the number of specimens.

8.4 Design strength — choosing the fracture probability and finding the allowable stress

Solving Eq. (55) (V = V₀) for σ, the stress at which a fracture probability P_f is allowed is

(For P_f ≪ 1, −ln(1 − P_f) ≈ P_f.) For P_f = 10⁻⁶ (1 in 1,000,000), σ_d/σ₀ is 0.251 for m = 10 and 0.501 for m = 20. For the default values (σ₀ = 446 MPa, m = 10), the stress is 281 MPa at P_f = 1 % and 112 MPa at 10⁻⁶ — it must be lowered to nearly a quarter of the mean strength of 424 MPa. Using the estimates of Fig. 3 (m_est = 9.34, σ_0,est = 447 MPa) gives 102 MPa, so underestimating m by just 7 % lowers the design stress by 9 %.

8.5 Breaking specimens one at a time brings the histogram closer to the theoretical distribution — the law of large numbers

In simulator tab ④, the specimens lined up on the shelf are put into the testing machine one at a time, or ten at a time, and the strengths are piled up in a histogram. The red curve is the theoretical expected number for the number n of specimens broken so far; with the bin width Δσ and the probability density of the strength p(σ),

(p is the derivative of P_s = exp[−(σ/σ₀)^m], and P_s is the probability of surviving a stress σ.) With the first few specimens the bars do not match the curve at all, but the more specimens are broken, the closer they come to the curve. Let us estimate how they approach it.

Let q (≈ Δσ p) be the probability of falling into a given bin. The quantity X_i, which is 1 if the i-th specimen falls into that bin and 0 if it does not, has mean q, and since X_i² = X_i, its variance is q − q². The specimens are independent, so the mean and variance of the number N = ΣX_i of the n specimens that fall into that bin are the respective sums, so

The relative scatter of the height of a bar decreases as 1/√n. For a bin with q = 0.1, it is 1 ± 0.95 specimens for n = 10 (scatter 95 %), 10 ± 3 for 100 specimens (30 %), and 100 ± 9.5 for 1000 specimens (9.5 %). The fact that the fraction N/n approaches the probability q — this is the law of large numbers, and the reason why the shape of the histogram approaches the red curve. Checking with N = 20000 specimens, the Kolmogorov–Smirnov statistic (the maximum difference between the empirical and theoretical distributions) was 0.0073, below the critical value 0.0096 at the 5 % significance level.

Scatter of the points of the Weibull plot. For any specimen, its strength σ converted to the theoretical fracture probability F(σ) = 1 − P_s(σ) is a uniform random number between 0 and 1 (because the probability that F(σ) ≤ u is exactly u). The F of the i-th of n specimens arranged in ascending order is the i-th value of n uniform random numbers, and its mean and variance are as follows (a standard result for order statistics).

The F_i = (i − 0.5)/n of Eq. (67) uses a value close to this mean as a representative. The vertical axis of the Weibull plot, y = ln ln[1/(1 − F)], is stretched where F is small, as y ≈ ln F, that is, dy ≈ dF/F. For n = 100 and i = 1 (the weakest specimen), the standard deviation of F is 0.0098, the same as its mean, so the scatter of y is 1 or more (1.29 by numerical calculation). It is 0.14 for i = 50 (the middle) and 0.23 for i = 100 (the strongest specimen). On the plot, the vertical position F_i is fixed and the horizontal position ln σ scatters, so the horizontal scatter is the vertical scatter divided by m — for m = 10, about ±0.13 in ln σ for the weakest point and about ±0.014 for the middle point. The few points at the lower left fluctuate the most and swing the slope of the least-squares line, so with a small number of specimens the estimate of m scatters widely (Table 5). If all 100 specimens are broken with the default values of the simulator, the estimate of m is 9.34, that of σ₀ is 447 MPa, the mean is 424 MPa, the weakest is 289 MPa, and the strongest is 510 MPa (the same pseudorandom numbers as in Fig. 3).

9. Volume Effects and Effective Volume

9.1 The volume effect under uniform tension

Apply uniform tension to specimens of the same material with volumes V₁ and V₂, and compare the stresses σ₁ and σ₂ at which they have the same fracture probability. If P_f in Eq. (55) is equal, the contents of the exponent are equal, so

For a component with 1000 times the volume of the test specimen, the strength is 1000^−1/10 = 0.501 times for m = 10 and still only 0.708 times for m = 20. For a volume 10 times larger, it is 0.794 times for m = 10. The median (the stress at P_f = 1/2) is σ₀(ln 2)^1/m, so the ratio of the medians also follows Eq. (73). In simulator tab ⑤ you can confirm this by actually breaking specimens of V₁ and 10V₁ with the model of §7.2 (the meaning of the symbols and a schematic figure are in §9.5).

9.2 When the stress is not uniform — the effective volume

In a bending test the stress differs from place to place. Going back to the fact that ln P_s in §7.1 was the sum of the contributions of the small pieces, we add up the risk for each small piece of volume dV using the stress σ(r) at its location (the compressed part σ ≤ 0 does not break, so it is not included).

Factoring out the maximum stress σ_max in the specimen, we get

V_eff is “the volume that gives the same risk if the maximum stress were applied uniformly.” Because of the m-th power, effectively only the places near the maximum stress are counted.

9.3 Finding the effective volume in 3-point bending by integration

A load is applied at the center of a rectangular bar with support span L, width b, and height h. The bending moment is 0 at the supports, maximum at the center, and varies linearly in between. The bending stress in the cross section is proportional to the distance y from the neutral axis. Measuring the position x from a support and taking the tensile side as y > 0, in the left half (0 ≤ x ≤ L/2) we have

The right half is symmetric, so we take twice the left half, and the compressed side y < 0 is not included. Substituting into Eq. (75), the integrals over x and y separate, and we get

With the substitution u = 2x/L, dx = (L/2)du, and ∫₀^L/2(2x/L)^mdx = (L/2)∫₀¹u^mdu = (L/2)/(m + 1). The integral over y is likewise (h/2)/(m + 1). Therefore

In 4-point bending (with spacing L_i between the inner loading points), the moment is constant in the inner section, so the integral in the length direction is L_i, and the two outer sections together contribute (L − L_i)/(m + 1), so that V_eff = bh[L_i + (L − L_i)/(m + 1)]/[2(m + 1)] = bh(mL_i + L)/[2(m + 1)²]. Putting L_i = L/2 and writing the result with the volume between the supports, V = bLh, we get

9.4 The ratio of bending strength to tensile strength

If the same bar is broken in uniform tension (V_eff = V) and in 3-point bending, the ratio of the σ_max that give the same fracture probability, obtained by replacing the volumes in Eq. (73) with effective volumes, is

For m = 10 it is 242^1/10 = 1.73, and for 4-point bending 40.3^1/10 = 1.45. That is, even for the same material, a strength 1.7 times the tensile strength is measured in 3-point bending and 1.45 times in 4-point bending (3-point bending gives 1.2 times the 4-point value). The larger the scatter, the larger the difference: 2.35 times in 3-point bending and 1.83 times in 4-point bending for m = 5, and 1.40 times and 1.25 times for m = 20, where the scatter is small.

9.5 The meaning of V₁, V₂, σ₁, and σ₂ — the volume effect seen through the largest crack

Let us confirm the meaning of the symbols once more. V₁ is the volume of the small specimen (for example, the specimen whose strength was measured), V₂ is the volume of a larger specimen or component (V₂ = nV₁), and σ₁ and σ₂ are the stresses at which each breaks with the same fracture probability; in simulator tab ⑤ they are matched at a fracture probability of 50 % (the median). When the stress is not uniform, as in a bending test, read V as the effective volume V_eff (§9.2).

Schematic 1 of ⑤ places a small specimen and a large specimen that contain cracks at the same density λ side by side, and marks the largest crack in each with a red circle. With the model of §7.2, we find the probability that a specimen of volume V contains no crack larger than a. The probability that a single crack exceeds a is (a_min/a)^k, so by the same calculation as for Eq. (59) (averaging over the Poisson distribution), we get

The a at which this becomes 1/2, that is, the median of the largest crack, is obtained by solving λV(a_min/a)^k = ln 2, which gives

This holds for λV ≥ ln 2 (with m = 2k). The strength of a specimen, σ = K_IC/(Y√(πa_max)), is a monotonic function that decreases as a_max increases, so the median of the strength equals the strength calculated from the median of the crack, and their ratio is (a_med,2/a_med,1)^−1/2 = (V₁/V₂)^1/m — the same as Eq. (73). For k = 5 (m = 10), a 10 times larger volume gives a 1.58 times larger largest crack and 0.794 times the strength; 100 times gives 2.51 times and 0.631 times; and 1000 times gives 3.98 times and 0.501 times. Checking with 4000 Monte Carlo runs, the ratios were 1.567 for the largest crack and 0.799 for the strength at 10 times, and 2.511 and 0.631 at 100 times. When the volume is increased, even with the same “density of cracks,” the “largest crack” becomes larger — this is what the volume effect really is.

9.6 The stress distribution in a bending test and the volume under high stress

We confirm from beam theory the bending stress σ ∝ M(x)y that was used in §9.3. The cross section of a bent beam rotates while remaining plane, and the strain at a distance y from the neutral axis is ε = y/ρ, where ρ is the radius of curvature. The bending moment obtained by integrating σ = Ey/ρ over the cross section is M = ∫σy dA = (E/ρ)∫y² dA = EI/ρ (I = bh³/12 for a rectangle of width b and height h), so eliminating E/ρ gives

In 3-point bending with load P and support span L, the reaction at each support is P/2, M = Px/2 (x ≤ L/2), M_max = PL/4 at the center, and σ_max = 3PL/(2bh²). For 4-point bending (P/2 applied at each of two points) in which the spacing of the inner loading points is 1/2 of the span between the outer supports (the ASTM C1161 configuration; simulator tab ⑤ and Eq. (79) use this configuration), the moment is constant in the inner section, M = PL/8, and σ_max = 3PL/(4bh²). In the configuration with an inner spacing of 1/3 (JIS R 1601), M = PL/6 and σ_max = PL/(bh²); putting L_i = L/3 into V_eff = bh(mL_i + L)/[2(m + 1)²] of §9.3 gives V_eff = V(m/3 + 1)/[2(m + 1)²], which for m = 10 is V/55.8, and the ratio of the bending strength to the tensile strength is 1.50 (1.45 for the 1/2 configuration).

Schematic 2 of ⑤ shades this σ = My/I and marks the part where σ ≥ 0.9σ_max in yellow. We now calculate the fraction of the volume that this part occupies. Let ξ ∈ [−1/2, 1/2] be the position measured in units of the span between the supports, and η = y/(h/2) ∈ [−1, 1] the position measured in units of half the height. In 3-point bending, σ/σ_max = (1 − 2|ξ|)η. In the layer of height η ≥ q on the tension side (half of the volume), the fraction of the length where σ ≥ qσ_max follows from 1 − 2|ξ| ≥ q/η as 1 − q/η, so we get

(For 4-point bending (inner spacing 1/2), σ/σ_max = min(1, 2 − 4|ξ|)η, and the fraction of the length becomes 1 − q/2η.) For q = 0.9, compared with 100 % in a tensile test, this is only 0.26 % in 3-point bending and 2.6 % in 4-point bending. Nearly the maximum tension acts only in a very thin layer on the lower surface, so if that layer contains no large crack, a high strength is measured. Counting this with a weight of the m-th power gives the effective volume; for m = 10, V_eff/V = 1/242 (3-point) and 1/40.3 (4-point). When the origin of fracture is a flaw on the surface, the same conversion is made with an effective area instead of a volume (Quinn and Morrell, 1991; Quinn, 2003). Beam theory does not include the contact stresses under the loading points and supports, and it also ignores the stress concentration at the shoulders of the grip section of a tensile specimen.

10. Thermal Shock

A hot ceramic put into cold water can crack. When it is cooled quickly, only the surface cools first and tries to contract, but it is held back by the interior, which is still hot, and cannot contract. As a result, a tensile stress arises at the surface. Ceramics are weak in tension, and the surface also has cracks such as grinding scratches, so the ceramic breaks under this stress.

10.1 Deriving the thermal stress at the surface

Let the surface be the x–y plane and the normal to the surface the z direction. Suppose that the temperature of a thin layer at the surface suddenly drops by ΔT relative to the interior. If it were free, this layer would shrink by αΔT in every direction, but we assume that (i) because it is connected to the thick interior, its in-plane dimensions cannot change (ε_x = ε_y = 0), and (ii) nothing restrains it in the direction normal to the surface (σ_z = 0). Adding the thermal strain −αΔT to Hooke’s law gives the following for the two in-plane directions (a tensile stress σ_y in the y direction contracts the x direction by νσ_y/E through the Poisson effect, so that amount is subtracted).

(Because σ_z = 0, there is no −νσ_z/E term.) The in-plane behavior is isotropic, so σ_x = σ_y ≡ σ. Setting ε_x = 0 gives

A bar constrained in only one direction would have σ = EαΔT, but here the layer is constrained from two in-plane directions, so the stress is multiplied by 1/(1 − ν) (1.28 times for alumina). The temperature difference at which this stress reaches the strength σ_f is the limit at which cracking begins.

We call this R (a quantity different from the R of the R-curve in §6.4) the thermal shock resistance parameter (it was Kingery (1955) who organized the resistance to thermal stress in terms of such a parameter). Calculated with the values in Table 1, it gives the following.

The reason silica glass has an R that is more than 20 times that of soda-lime glass (about 16 times that of alumina) and is especially resistant to thermal shock is that its coefficient of linear thermal expansion is only about 1/16 of that of soda-lime glass. This is why red-hot silica glass does not crack when it is put into water. Silicon nitride has a high strength and a small thermal expansion, so among ceramics it is the representative material that resists thermal shock. Alumina has a high strength, but its Eα is large, so its R stays at 90 K.

10.2 Temperature inside the sample — transient heat conduction and the Fourier number

Eq. (87) was the value at the instant when “only the surface has cooled by ΔT and the interior has not cooled yet.” Simulator tab ⑦ calculates how the temperature spreads through a square sample over time and how the thermal stress changes. Let us start with the temperature. Consider the heat entering and leaving a small volume dV. The heat flux (the heat that flows per unit area per unit time) is given by Fourier’s law q = −k∇T (k is the thermal conductivity), and the net heat that flows in, −∇·q dV, raises the temperature (density ρ, specific heat c), so ρc ∂T/∂t = k∇²T, that is

where κ is the thermal diffusivity. (This is a two-dimensional equation: we consider the cross section of a long square bar and assume that the temperature does not change along its length.) Dividing the coordinates by the side length L (x̂ = x/L, ŷ = y/L) and measuring time by

gives an equation that depends on neither the material nor the size. The Fourier number Fo is time measured in units of “the time L²/κ that heat needs to spread over a distance L.” The problem of ⑦ is a square cross section at a uniform initial temperature T₀; at time 0 the top face is held at T_top and the bottom face at T_bottom, and the side faces are adiabatic (no heat enters or leaves). This corresponds to the most severe quench or sudden heating, in which the surface temperature changes in an instant (the case of an infinite Biot number in the box of §10.1).

Eq. (89) is linear in the temperature, so if we find one solution φ(ŷ, Fo) for “the top face at 1, the initial state and the bottom face at 0,” we can write, by superposition, T = T₀ + (T_top − T₀)φ(ŷ, Fo) + (T_bottom − T₀)φ(1 − ŷ, Fo). Because the side faces are adiabatic and the initial temperature is uniform, the temperature is a function of ŷ only (one-dimensional). φ settles into the straight line 1 − ŷ as time passes, so putting φ = 1 − ŷ − u, u is the solution of heat conduction that is 0 at both ends and starts from 1 − ŷ. By separation of variables we put u = Σ b_n sin(nπŷ) exp(−n²π²Fo), and expanding the initial shape in a sine series, integration by parts gives b_n = 2∫₀¹(1 − ŷ) sin(nπŷ) dŷ = 2/(nπ), so we get