5.3 Crack tip field

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": 16,
    "text.usetex": True,
    "font.family": "Computer Modern"
})

# 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} 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 settings
    plt.rcParams.update({'font.size': 25})
    plt.rcParams['figure.figsize'] = [10, 10]
    plt.rcParams['figure.dpi'] = 300

    # Plot the displacement field
    plt.clf()
    fig = plt.figure(1)
    ax = fig.add_subplot(111)
    plot = ax.contourf(x_grid, y_grid, sigma_vm, contour_vector, extend='max', cmap=cm)
    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', dpi=300)
    plt.clf()

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