5.3 Crack tip field

Williams expansion

Linear elastic fracture mechanics is based on stress intensity factors, specifically the $r^{-1/2}$ stress singularity related to the first term of Williams' series. The Williams eigenfunction expansion form can have positive and negative terms, as it is an infinite series. The displacement field eigenfunction expansion for pure mode I is as follows:

$$u_{i}(r,\theta) = \underbrace{\sum_{n=-\infty}^{-1} A_{n} g_{{\rm I},i}(\theta, n) \ r^{\frac{n}{2}}}_{\text{supersingular terms}} + \underbrace{A_0 g_{{\rm I},i}(\theta, 0)}_{\text{translations}} + \underbrace{A_1 g_{{\rm I},i}(\theta, 1) \sqrt{r}}_{\text{singular}} + \underbrace{\sum_{n=2}^{\infty} A_n g_{{\rm I},i}(\theta, n) \ r^{\frac{n}{2}}}_{\text{subsingular terms}}$$

The formulation implemented in CrackPy (from crackpy.fracture_analysis.crack_tip import williams_displ_field, williams_stress_field) is as follows [Melching23, Code, arXiv]:

$$\sigma_{ij}(r,\theta) = \sum_{n=1}^{\infty}{r^{\frac{n}{2}-1}\ \left( A_n f_{\text{I},ij}(\theta,n) + B_n f_{\text{II},ij}(\theta,n) \right)}$$ $$u_{i}(r,\theta) = \sum_{n=0}^{\infty} \frac{r^{\frac{n}{2}}}{2\mu} \ \left( A_n g_{\text{I},i}(\theta,n) + B_n g_{\text{II},i}(\theta,n) \right)$$

The parameters $A_n, B_n \in \mathbb{R}$ are called Williams coefficients and depend on the loading conditions. The trigonometric eigenfunctions $f_{\text{I},ij},f_{\text{II},ij}, g_{\text{I},i},g_{\text{II},i}$ for the stress tensor are defined as follows:

$$f_{\text{I},11}(\theta,n)=\frac{n}{2} \left\{ \left[ 2+(-1)^n+\frac{n}{2} \right] \cos \left(\left(\frac{n}{2}-1\right) \theta \right) - \left( \frac{n}{2}-1 \right) \cos \left(\left(\frac{n}{2}-3\right) \theta \right) \right\}$$ $$f_{\text{II},11}(\theta,n)=\frac{n}{2} \left\{ \left[ -2+(-1)^n-\frac{n}{2} \right] \sin \left(\left(\frac{n}{2}-1\right) \theta \right) + \left( \frac{n}{2}-1 \right) \sin \left(\left(\frac{n}{2}-3\right) \theta \right) \right\}$$ $$f_{\text{I},22}(\theta,n)=\frac{n}{2} \left\{ \left[ 2-(-1)^n-\frac{n}{2} \right] \cos \left(\left(\frac{n}{2}-1\right) \theta \right) + \left( \frac{n}{2}-1 \right) \cos \left(\left(\frac{n}{2}-3\right) \theta \right) \right\}$$ $$f_{\text{II},22}(\theta,n)=\frac{n}{2} \left\{ \left[ -2-(-1)^n+\frac{n}{2} \right] \sin \left(\left(\frac{n}{2}-1\right) \theta \right) - \left( \frac{n}{2}-1 \right) \sin \left(\left(\frac{n}{2}-3\right) \theta \right) \right\}$$ $$f_{\text{I},12}(\theta,n)=f_{\text{I},21}(\theta,n)=\frac{n}{2} \left\{ \left( \frac{n}{2} - 1 \right) \sin \left(\left(\frac{n}{2}-3\right) \theta \right) - \left[ \frac{n}{2}+(-1)^n \right] \sin \left(\left(\frac{n}{2}-1\right) \theta \right) \right\}$$ $$f_{\text{II},12}(\theta,n)=f_{\text{II},21}(\theta,n)=\frac{n}{2} \left\{ \left( \frac{n}{2} - 1 \right) \cos \left(\left(\frac{n}{2}-3\right) \theta \right) - \left[ \frac{n}{2}-(-1)^n \right] \cos \left(\left(\frac{n}{2}-1\right) \theta \right) \right\}$$

and the functions for the displacement part are defined as

$$g_{\text{I},1}(\theta,n)=\left[ \kappa + (-1)^n + \frac{n}{2} \right] \cos \left( \frac{n}{2} \theta \right) - \frac{n}{2} \cos \left( \left( \frac{n}{2} - 2 \right) \theta \right)$$ $$g_{\text{II},1}(\theta,n)=\left[ -\kappa + (-1)^n - \frac{n}{2} \right] \sin \left( \frac{n}{2} \theta \right) + \frac{n}{2} \sin \left( \left( \frac{n}{2} - 2 \right) \theta \right)$$ $$g_{\text{I},2}(\theta,n)=\left[ \kappa - (-1)^n - \frac{n}{2} \right] \sin \left( \frac{n}{2} \theta \right) + \frac{n}{2} \sin \left( \left( \frac{n}{2} - 2 \right) \theta \right)$$ $$g_{\text{II},2}(\theta,n)=\left[ \kappa + (-1)^n - \frac{n}{2} \right] \cos \left( \frac{n}{2} \theta \right) + \frac{n}{2} \cos \left( \left( \frac{n}{2} - 2 \right) \theta \right)$$

Here, $\kappa$ is the Kolosov constant, namely, $\kappa = (3-\nu) / (1+\nu)$ under plane stress conditions and $\kappa=3-4\nu$ for plane strain, where $\nu>0$ is the Poisson ratio and $\mu = E/(2(1+\nu))$ is the shear modulus.

While $A_n$ and $B_n$ have the unit $\text{MPa} \cdot \text{mm}^{1-n/2}$ (from crackpy.fracture_analysis.crack_tip import unit_of_williams_coefficients).

The stress intensity factors and the T-stress are then:

$K_{\rm I} = \sqrt{2\pi} \cdot A_1$

$K_{\rm II} = -\sqrt{2\pi} \cdot B_1$

$T = 4 \cdot A_2$

You can use the following code to plot the Williams field:

'''
This script plots the Williams field for different number of terms.

Needed:
    - Williams expansion coefficients A and B

Output:
    - Plots of the Williams field

'''

import os

from crackpy.fracture_analysis.crack_tip import williams_displ_field, williams_stress_field
from crackpy.structure_elements.material import Material
from crackpy.fracture_analysis.optimization import Optimization
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable

# Set matplotlib settings
plt.rcParams.update({
    "font.size": 20,
    "text.usetex": True,
    "font.family": "Computer Modern",
    "figure.figsize": [10, 10],
    "figure.dpi": 300
})

# Parameters
K_I = 23.7082070388 * np.sqrt(1000)  # MPa.mm^0.5
T = -7.2178311164  # MPa

a_1 = K_I / np.sqrt(2 * np.pi)
a_2 = T / 4
a_3 = -2.8672068762
a_4 = 0.0290671612

# Define the Williams expansion coefficients
coefficients_a = [a_1, a_2, a_3, a_4]
coefficients_b = [0, 0, 0, 0]
coefficients_n = [1, 2, 3, 4]

# Output path
OUTPUT_PATH = os.path.join('Williams_field')

# check if output path exists
if not os.path.exists(OUTPUT_PATH):
    os.makedirs(OUTPUT_PATH)

# Define the material
material = Material(E=72000, nu_xy=0.33, sig_yield=350.0)

# Loop over the number of terms and add one term at a time to the Williams expansion
for index in range(1, len(coefficients_a) + 1):
    print(f'Plotting Williams field for {index} {"term" if index == 1 else "terms"}.')
    for i in range(index):
        print(f'a_{coefficients_n[i]}= {coefficients_a[i]:.2f} N * mm^({-1 - coefficients_n[i] / 2})'
              f', b_{coefficients_n[i]}= {coefficients_b[i]:.2f} N * mm^({-1 - coefficients_n[i] / 2})')

    min_radius = 0.01
    max_radius = 100.0
    tick_size = 0.01

    r_grid, phi_grid = np.mgrid[min_radius:max_radius:tick_size, -np.pi:np.pi:tick_size]
    disp_x, disp_y = williams_displ_field(coefficients_a[:index],
                                          coefficients_b[:index],
                                          coefficients_n[:index],
                                          phi_grid, r_grid, material)
    sigma_xx, sigma_yy, sigma_xy = williams_stress_field(coefficients_a[:index],
                                                         coefficients_b[:index],
                                                         coefficients_n[:index],
                                                         phi_grid, r_grid)
    sigma_vm = np.sqrt(sigma_xx ** 2 + sigma_yy ** 2 - sigma_xx * sigma_yy + 3 * sigma_xy ** 2)

    x_grid, y_grid = Optimization.make_cartesian(r_grid, phi_grid)

    # Matplotlib plot
    number_colors = 120
    number_labes = 5
    legend_limit_max = 100
    Legend_limit_min = 0.0
    cm = 'coolwarm'

    # Define contour and label vectors
    contour_vector = np.linspace(Legend_limit_min, legend_limit_max, number_colors, endpoint=True)
    label_vector = np.linspace(Legend_limit_min, legend_limit_max, number_labes, endpoint=True)
    label_vector = np.round(label_vector, 2)

    # Plot the displacement field
    plt.clf()
    fig = plt.figure(1)
    ax = fig.add_subplot(111)

    # Plot the crack tip field
    plot = ax.contourf(x_grid, y_grid, sigma_vm, contour_vector, extend='max', cmap=cm)
    # Highlight the crack path
    ax.plot([0, np.min(x_grid)], [0, 0], 'k', linewidth=2)

    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.2)
    plt.colorbar(plot, ticks=label_vector,
                 cax=cax,
                 label='Von Mises eqv. stress [$\\mathrm{MPa}$]')
    ax.set_xlabel('$x$ [$\\mathrm{mm}$]')
    ax.set_ylabel('$y$ [$\\mathrm{mm}$]')
    ax.axis('image')
    ax.set_xlim(-50, 50)
    ax.set_ylim(-50, 50)
    title_str_a = ''
    title_str_b = ''
    for i in range(index):
        title_str_a += f'$a_{{{coefficients_n[i]}}} = {coefficients_a[i]:.2f} \\, \\mathrm{{N \\cdot mm^{{{-1 - coefficients_n[i] / 2}}}}}$'
        title_str_b += f'$b_{{{coefficients_n[i]}}} = {coefficients_b[i]:.2f} \\, \\mathrm{{N \\cdot mm^{{{-1 - coefficients_n[i] / 2}}}}}$'
        if i < index - 1:
            title_str_a += ', '
            title_str_b += ', '

    fig.suptitle(f'Williams field with {index} {"term" if index == 1 else "terms"}', y=0.95)
    ax.set_title(title_str_a + '\n' + title_str_b, fontsize=14)

    output_file = os.path.join(OUTPUT_PATH, f'Williams_{index}.png')
    plt.savefig(output_file, bbox_inches='tight')
    plt.clf()

The following plot shows the influence of the Williams series coefficients on the crack tip field:

Christopher-James-Patterson (CJP) model

...