<figure><img src="image_1.png" /></figure>

HAL
open science

# Convergence analysis of crack features extraction using conjugate work integral

Héloïse Dandin, Nicolas Chevaugeon, Julien Réthoré

► To cite this version:

Héloïse Dandin, Nicolas Chevaugeon, Julien Réthoré. Convergence analysis of crack features extraction using conjugate work integral. 2025. hal-05018756v2

HAL Id: hal-05018756

https://hal.science/hal-05018756v2

Preprint submitted on 5 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Convergence analysis of crack features extraction using conjugate work integral

Héloïse Dandin¹, Nicolas Chevaugeon², Julien Réthoré²

¹Univ Brest, Bretagne INP, IRDL, UMR CNRS 6027, F-29200 Brest

²Nantes Université, Ecole Centrale Nantes, CNRS, GeM, UMR 6183, F-44000 Nantes

November 5, 2025

**Abstract**

**Background:** The extraction of Williams higher-order terms ($n > 1$) is crucial, since these terms control the crack stability and propagation, among other phenomena. Nevertheless, current numerical methods strongly depend on the projection domain and/or are unstable, and thus not suited to analysing experimental data.

**Purpose:** This paper aims to investigate the method based on the Bueckner-Chen conjugate work integral, which provides a bilinear form in which Williams series terms are orthogonal, to retrieve all Williams higher-order coefficients. We reformulate the original path-independent integral into an Equivalent Domain Integral (EDI), more robust when dealing with experimental data or finite element discretisation, and compare the convergence of both approaches.

**Methods:** Williams fields are imposed on a meshed body with decreasing element sizes. The different error sources, i.e. numerical integration, stress interpolation and Finite Element (FE) discretisation, are decoupled by applying relevant boundary conditions to the linear elastic problem.

**Results:** Predictably, the error due to FE discretisation is higher than the other types of errors. Bueckner-Chen integral converges markedly faster than the *J*-integral on the same domain and provides a better approximation of the Stress Intensity Factor (SIF). As expected, the EDI formulation is more accurate and efficient to retrieve Williams higher-order terms.

**Conclusion:** This work presents a convergence analysis on a robust and efficient method based on the Bueckner-Chen integral to retrieve Williams higher-order coefficients.

**Keywords:** Williams coefficients, higher-order terms, Bueckner-Chen integral, equivalent domain integral, convergence analysis

Licence: CC-BY @TheAuthors

# 1 Introduction

In Linear Elastic Fracture Mechanics (LEFM), the asymptotic displacement field in the vicinity of the crack tip is accurately described by Williams Eigenfunction Expansion Form (EEF) [Williams, 1957]. These infinite series are parametrised by coefficients denoted $A_n, B_n$ referring to modes I, resp. II, that have various physical meanings. The coefficients at order $n = 1$ are called *singular* and are proportional to the Stress Intensity Factors (SIF) $K_I$ and $K_{II}$. They are used to effectively predict crack propagation, e.g. Paris' law in cyclic fatigue, maximum hoop stress criterion for crack orientation in brittle materials, Kanninen's law in dynamic fracture [Kanninen and Popelar, 1985]. Moreover, the energy release rate $G$ can be computed from $K_I$ and $K_{II}$ or from the *J*-integral [Rice, 1968]. Consequently, in pure mode I the SIF $K_I$ can be derived from the *J*-integral: $J = G = K_I^2/E^*$ (with $1/E^* = (1 - \nu^2)/E$ in plane strain).

The influence of some higher-order or *subsingular* terms has also been explored in a number of different situations: brittle failure of polymers [Cotterell, 1966, Eid et al., 2023], fatigue crack propagation [Bellecave et al., 2014, Bahaj Filali et al., 2025], failure of concrete [Nguyen et al., 2020], among others. For instance, the coefficient $A_2$ is related to the *T*-stress, which influences the stability of the crack's direction, while the third coefficient refers to the *B*-stress and can be
responsible for the stability of the crack propagation. The fourth coefficient describes the sign of the maximum shear stress gradient along the crack axis.

Consequently, the extraction of Williams coefficients is a question that still captivates the fracture mechanics community. Most of the time, this extraction is to be performed using displacement field data obtained either from finite element numerical simulation or more recently from displacement field measured by digital image correlation during experiments.

The kinematic approach consists in minimising the distance between a displacement field and Williams EEF. To ensure the convergence of the results, the minimisation is performed over several orders of the series or over the so-called K-dominance zone if only SIF are sought for. One can mention the work by Kim et al. [2010] in computational mechanics or that of McNeill et al. [1987], Hamam et al. [2007] in experimental mechanics. It has also been proposed to directly extract Williams coefficients by considering their amplitude as kinematic degrees of freedom [Roux and Hild, 2006, Xiao and Karihaloo, 2007, Passieux et al., 2010, Roux et al., 2009]. On the other hand, energetic approaches are based on path-independent integrals, such as the M-integral [Attigui and Petit, 1997] or the interaction integral [Chen and Shield, 1977, Bui, 1983] which allow extracting only one Williams coefficient at a time (usually the SIF). All these methods can be written under the form of a scalar product between the actual displacement field and an extractor or auxiliary field that can be associated with a kinematic or energetic norm [Réthoré et al., 2011]. Bueckner-Chen's conjugate work integral [Bueckner, 1973, Chen, 1985] stands out because the terms of the Williams series are orthogonal in the sense of the bilinear form it provides.

The Bueckner integral [Bueckner, 1973] is based on Betti's reciprocal theorem: let us consider two symmetrical load systems denoted by $Q^{(1)}$ and $Q^{(2)}$. The work accomplished by the tractions of $Q^{(1)}$ through the displacements of $Q^{(2)}$ is equal to the work done by the tractions of $Q^{(2)}$ through the displacements of $Q^{(1)}$. An application was given by Rice [1972]: if the displacements and SIF are known for any arbitrarily chosen load system $Q^{(2)}$, then it is possible to determine the SIF for any other load system $Q^{(1)}$. A decade later, Chen [1985] continued this work and demonstrated, for a carefully chosen auxiliary load case $Q^{(2)}$, the orthogonality property of Williams EEF with respect to the Bueckner integral. Thanks to this property, he was able to retrieve the SIFs and he proved the equivalence of his approach with the J-integral in mode I [Rice, 1968]. More generally, this orthogonality property lets any coefficient of the Williams series be retrieved regardless of the order or mode. To this end, Melching and Breitbarth [2023] developed an automatic procedure based on the Bueckner-Chen integral to extract any subsingular term. They proved the path-independence of the Bueckner-Chen Line Integral (LI) with varying crack distance and under several loading conditions.

For these reasons, we focussed our research on this particular integral and continued Melching and Breitbarth's work to verify the numerical convergence of the Bueckner-Chen integral. Furthermore, for a robust estimation of SIF, path integrals are usually transformed into Equivalent Domain Integrals (EDI) using the concept of virtual crack extension field [Destuynder et al., 1983, Moran and Shih, 1987]. We thus propose a new formulation of the crack features extraction problem based on the Bueckner-Chen EDI. This paper is divided into two main parts. Section 2 details the analytical methods employed, from Williams series to the determination of its coefficients with Bueckner-Chen LI and EDI. Then, Section 3 presents the numerical study that has been led to verify the convergence of the extraction method based on the Bueckner-Chen integral, both with the LI and EDI formulations. The results are also compared with the J-integral. In this study, we endeavour to separate and quantify the numerical errors due to numerical integration, stress interpolation and FE discretisation.

# 2 Analytical methods

## 2.1 Williams series

We work within the framework of Linear Elastic Fracture Mechanics (LEFM) and consider a cracked plate with a local polar coordinate system centred on the crack tip (see Fig. 1). Williams Eigenfunction Expansion Form provides the following mathematical expression for the 2D stress and displacement fields for a semi-infinite straight crack:

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

where $A_n, B_n \in \mathbb{R}$ are Williams coefficients. The stress trigonometric functions for modes I and II are defined as:

$$ f_{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\} \quad (3a) $$

$$ f_{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\} \quad (3b) $$

$$ f_{I,12}(\theta, n) = f_{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\} \quad (3c) $$

$$ f_{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\} \quad (3d) $$

$$ f_{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\} \quad (3e) $$

$$ f_{II,12}(\theta, n) = f_{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\} \quad (3f) $$

The displacement trigonometric functions are defined as:

$$ g_{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) \quad (4a) $$

$$ g_{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) \quad (4b) $$

$$ g_{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) \quad (4c) $$

$$ g_{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) \quad (4d) $$

These expressions depend on several material parameters, namely,

* the Young's modulus $E$,

* the Poisson ratio $\nu > 0$,

* the Kolosov constant, defined by $\kappa = (3 - \nu)/(1 + \nu)$ under plane stress and $\kappa = 3 - 4\nu$ under plane strain conditions,

* the shear modulus $\mu = E/(2(1 + \nu))$.

The terms of order $n < 0$ in the EEF don't have any physical meaning. However, they can still be considered in regions excluding the crack tip for estimating the position of an equivalent elastic crack tip or the size of a process zone when experimental displacement fields [Roux et al., 2009] or cohesive zone model simulations [Réthoré and Estevez, 2013] are analysed. In the following, we will restrict our study to positive orders of the EEF, i.e. singular and subsingular terms.

## 2.2 Bueckner-Chen Line Integral and determination of Williams subsingular terms

Melching and Breitbarth [2023] implemented the Bueckner-Chen conjugate work integral [Bueckner, 1973, Chen, 1985] to extract the coefficients of Williams subsingular terms. In the following, we recall this procedure for the path-independent Line Integral (LI).

Let $\Gamma$ be a closed path around the crack tip, $n_j$ an outward unit normal to $\Gamma$ and $x_j$ coordinates with $x_3$ in the thickness direction as illustrated in Fig. 1. Let also $(u^{(1)}, \sigma^{(1)})$ and $(u^{(2)}, \sigma^{(2)})$ be pairwise-coupled fields corresponding to linear elastic solutions of two different loadings for a planar crack. The Bueckner integral [Bueckner, 1973] is a conjugate work integral whose path-independence property stems from Betti's theorem. It writes:

$$ I_{\Gamma}(u_i^{(1)}, \dots, u_i^{(2)}, \dots) = \int_{\Gamma} \left( \sigma_{ij}^{(1)} u_i^{(2)} - \sigma_{ij}^{(2)} u_i^{(1)} \right) n_j \, d\Gamma \quad (5) $$

where $d\Gamma$ is the arc length along $\Gamma$ and the stress fields are given by Hooke's law, $\sigma = \mathbb{C} : \nabla^s u$. Chen [1985] established that:

$$ I_{\Gamma} \left( \frac{r^{n/2}}{2\mu} g_{J,i}(n), \dots, \frac{r^{m/2}}{2\mu} g_{K,i}(m), \dots \right) = \begin{cases} -\frac{\pi(\kappa+1)}{\mu}(-1)^{n+1}n & \text{if } m+n=0 \text{ and } J=K \\ 0 & \text{else} \end{cases} \quad (6) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Integration domain and virtual crack extension field</figcaption></figure>

what states for the orthogonality of Williams eigenfunctions with respect to the updated conjugate work integral. In practice, the solution fields for which we seek to extract the unknown coefficients $A_n, B_n$ is denoted by the subscript (1) and the auxiliary fields by the superscript (2). According to Chen, the solution displacement $\mathbf{u}^{(1)}$ is defined by Eq. 2 and the auxiliary displacement field $\mathbf{u}^{(2)}$ writes:

$$u_i^{(2)} = u_i^{(2)}(m) = \frac{r^{m/2}}{2\mu} (\hat{A}_m g_{I,i} + \hat{B}_m g_{II,i}) \quad (7)$$

Inserting into Eq. 5, we denote:

$$I_{\Gamma}(m) = I_{\Gamma}(u_{i}^{(1)}, \dots, u_{i}^{(2)}(m), \dots) \quad (8)$$

Due to the orthogonality property, we get the relations between the sought coefficients $A_n, B_n$ at order $n$ when $m = -n$. For pure mode I ($\hat{B}_m = 0$), resp. II ($\hat{A}_m = 0$), they write:

$$I_{\Gamma}^{I}(-n) = -\frac{\pi(\kappa + 1)}{\mu} (-1)^{n+1} n A_n \hat{A}_m \quad (9a)$$

$$I_{\Gamma}^{II}(-n) = -\frac{\pi(\kappa + 1)}{\mu} (-1)^{n+1} n B_n \hat{B}_m \quad (9b)$$

Finally, the values for $A_n, B_n$ are obtained by setting $\hat{A}_m = \hat{B}_m = 1$ [Melching and Breitbarth, 2023]:

$$A_n = -\frac{\mu}{\kappa + 1} \frac{1}{\pi n (-1)^{n+1}} I_{\Gamma}^{I}(-n) \quad (10a)$$

$$B_n = -\frac{\mu}{\kappa + 1} \frac{1}{\pi n (-1)^{n+1}} I_{\Gamma}^{II}(-n) \quad (10b)$$

## 2.3 Bueckner-Chen Equivalent Domain Integral

In this section, we rewrite the Bueckner-Chen LI into an Equivalent Domain Integral (EDI). The demonstration is similar to those of other path-independent integrals in fracture mechanics [Destuynder et al., 1983, Moran and Shih, 1987].

Let us define an area $A$ delimited by a closed integration path $C = \Gamma + \Gamma^+ - \Gamma_\varepsilon + \Gamma^-$, as illustrated in Fig. 1. The Bueckner-Chen integral along $\Gamma_\varepsilon$ writes:

$$I_{\Gamma_{\varepsilon}} = -I_{C} + I_{\Gamma} + I_{\Gamma^{+} + \Gamma^{-}} \quad (11)$$
where $I_{\Gamma\varepsilon} = I_{\Gamma}$ due to the path-independence property of the integral. We then define a scalar, continuous and differentiable virtual crack extension field such that:

$$ q = \begin{cases} 0 & \text{on } \Gamma \\ 1 & \text{on } \Gamma_{\varepsilon} \end{cases} \qquad (12) $$

Upon insertion into Eq. 5, the contribution of this virtual field transforms Eq. 11 into:

$$ I_{\Gamma\varepsilon}(u_i^{(1)}, \dots, u_i^{(2)}, \dots) = - \oint_C (\sigma_{ij}^{(1)} u_i^{(2)} - \sigma_{ij}^{(2)} u_i^{(1)}) n_j q \, d\Gamma + \int_{\Gamma^+ + \Gamma^-} (\sigma_{ij}^{(1)} u_i^{(2)} - \sigma_{ij}^{(2)} u_i^{(1)}) n_j q \, d\Gamma \quad (13) $$

where the term over $\Gamma$ vanished. We also assume tension-free crack faces, which translates into $\sigma_{ij}^{(1)} n_j = \sigma_{ij}^{(2)} n_j = 0$, thereby simplifying the last term of the equation:

$$ I_{\Gamma^+ + \Gamma^-} = 0 \qquad (14) $$

Applying the divergence theorem transforms the closed contour integral over $C$ into a domain integral over $\mathcal{A}$:

$$ \begin{aligned} I_{\Gamma\varepsilon} &= - \int_{\mathcal{A}} \frac{\partial}{\partial x_j} \left[ (\sigma_{ij}^{(1)} u_i^{(2)} - \sigma_{ij}^{(2)} u_i^{(1)}) q \right] dA \\ &= - \int_{\mathcal{A}} \left[ \frac{\partial}{\partial x_j} (\sigma_{ij}^{(1)} u_i^{(2)} - \sigma_{ij}^{(2)} u_i^{(1)}) q + (\sigma_{ij}^{(1)} u_i^{(2)} - \sigma_{ij}^{(2)} u_i^{(1)}) q_{,j} \right] dA \quad (15) \\ &= - \int_{\mathcal{A}} \left[ (\sigma_{ij,j}^{(1)} u_i^{(2)} - \sigma_{ij,j}^{(2)} u_i^{(1)} + \sigma_{ij}^{(1)} u_{i,j}^{(2)} - \sigma_{ij}^{(2)} u_{i,j}^{(1)}) q + (\sigma_{ij}^{(1)} u_i^{(2)} - \sigma_{ij}^{(2)} u_i^{(1)}) q_{,j} \right] dA \end{aligned} $$

The principle of equilibrium states that $\sigma_{ij,j}^{(1)} = \sigma_{ij,j}^{(2)} = 0$. Furthermore, the major and minor symmetry properties of $\mathcal{C}$ in isotropic linear elasticity lead to:

$$ \sigma_{ij}^{(1)} u_{i,j}^{(2)} = C_{ijkl} u_{k,l}^{(1)} u_{i,j}^{(2)} = \sigma_{ij}^{(2)} u_{i,j}^{(1)} \qquad (16) $$

With these two simplifications, the first member of the integral is dropped. The contour integral over $\Gamma_{\varepsilon}$ is equal to the domain integral over $\mathcal{A}$, which provides the final expression for the Bueckner-Chen EDI:

$$ I_{\mathcal{A}}(u_i^{(1)}, \dots, u_i^{(2)}, \dots) = \int_{\mathcal{A}} (\sigma_{ij}^{(2)} u_i^{(1)} - \sigma_{ij}^{(1)} u_i^{(2)}) q_{,j} \, dA \qquad (17) $$

The orthogonality property in Eq. 9 is still valid: the Williams coefficients from the equilibrated elastic field $(u^{(1)}, \sigma^{(1)})$ can thus be retrieved using the procedure described in Section 2.2 (Eq. 10). In the context of a crete finite element solution over a mesh, the resolution of the extraction problem with the domain integral is simplified compared to the procedure based on contour integration. Indeed, stress values, typically computed at Gauss points of the mesh, do not need to be interpolated at contour points (typically edges of the mesh, where the finite element stress field is discontinuous).

## 3 Convergence study

In this section, we verify the convergence of the Bueckner-Chen LI and EDI with various numerical approximations to decouple the effects of numerical integration, stress interpolation and FE discretisation. We also compare the accuracy and efficiency of the proposed integral to those of the $J$-integral.
## 3.1 Numerical implementation

### 3.1.1 Mechanical problem

We consider a rectangle plate with a single edge crack under plane stress hypothesis. The rectangle is defined by its height $H = 200$ mm and its width $W = 200$ mm; the crack is centred on the edge and its length is $a = 50$ mm. The material behaviour is isotropic linear elastic and is parameterised by the Young's modulus $E = 1$ MPa and the Poisson ratio $\nu = 0.3$.

A coarse mesh is generated thanks to Gmsh's Frontal Delaunay algorithm [Geuzaine and Remacle, 2009] using triangle elements with linear shape functions for the displacement field; the crack is described by twin nodes along the crack path. Features that facilitate the numerical evaluation of the Bueckner-Chen integrals, i.e. a path $\Gamma$ and a domain $A$, are added to the design according to the following specifications illustrated in Fig. 2:

* Discretised circle arcs $\Gamma^h$ and $\Gamma_{int}^h$, as illustrated in Fig. 1, will be used to evaluate the LI and EDI. They are defined by radii $R$ and $R_{int}$ respectively.

* The domain $A^h$ used for EDI is composed of all elements between $\Gamma^h$ and $\Gamma_{int}^h$.

* Element length is set to $h_{out} = 16$ mm on the outer boundaries (crack faces excluded).

* Element length is set to $h = 1.6$ mm on $\Gamma^h$, $\Gamma_{int}^h$ and on the crack faces between the crack tip and $\Gamma$.

Two domains are considered to evaluate the EDI: one is such that $R_{int} = 5$ mm and excludes the crack tip, while the second is defined by $R_{int} = 0$ mm. The contour $\Gamma^h$ used for LI is the one obtained when $R_{int} = 0$ mm.

<figure><img src="image_4.png" /><figcaption>Figure 2: Schematic representation of the mesh near the crack tip and features used for numerical evaluation of the integrals</figcaption></figure>

More refined meshes are built by splitting the coarse mesh four times. The final element lengths and number of elements (lines or triangles) used for integration are given in Tab. 1.

<table><caption>Table 1: Mesh specifications</caption><thead><tr><th></th><th>$h/R$</th><th>0.16</th><th>0.08</th><th>0.04</th><th>0.02</th><th>0.01</th></tr></thead><tbody><tr><td rowspan="3">number of<br>elements:</td><td>$\Gamma^h$ ($R = 10$ mm)</td><td>40</td><td>80</td><td>160</td><td>320</td><td>640</td></tr><tr><td>$A^h$ ($R = 10$ mm, $R_{int} = 0$ mm)</td><td>344</td><td>1 376</td><td>5 504</td><td>22 016</td><td>88 064</td></tr><tr><td>$A^h$ ($R = 10$ mm, $R_{int} = 5$ mm)</td><td>278</td><td>1 112</td><td>4 448</td><td>17 792</td><td>71 168</td></tr></tbody></table>

The convergence study aims to investigate the convergence rate and accuracy of the numerical integration implementation. To this end, we seek to recover the Williams coefficients used to define solution displacements $u^W$ and stresses $\sigma^W$ written under the form of Williams EEF in pure mode I with:

$$A_n^W = R^{1/2-n/2} \quad \forall n = 1 \dots 4 \qquad (18)$$
Thereby, all Williams orders contribute similarly to stresses along $\Gamma^h$, while K-dominance is ensured within the integration domain.

Finally, we define the virtual crack extension field $q$ as a linear function between $\Gamma^h$ and $\Gamma_{int}^h$ and its gradient will be computed analytically assuming $\Gamma^h$ and $\Gamma_{int}^h$ are circles. This projection circumvents the need for an interpolation technique and the ratio between the element size and the area lost by discretising the circle is sufficiently small to justify the validity of this hypothesis.

### 3.1.2 Convergence indicators

In this section, we distinguish between two sets of integration points:

* triangle integration points (one point per element) refer to the FE mesh and the Gaussian quadrature method used for solving Boundary Value Problems (BVP),

* edge integration points are used to evaluate the Bueckner-Chen LI, the number of points per edge depends on the chosen order for numerical integration.

We define the relative error between $A_n^W$ (used to define solution fields) and $A_n$ computed with the Bueckner-Chen integral:

$$ \text{err} = \left| \frac{A_n^W - A_n}{A_n^W} \right| \qquad (19) $$

This error is theoretically independent of the integration contour [Chen, 1985, Melching and Breitbarth, 2023]. Nevertheless, it depends on the mesh, the numerical integration approximation and the definition of the input fields $u^{(1)}$ and $\sigma^{(1)}$. Therefore, we characterise the influence of these parameters by evaluating:

(A) the numerical integration error that is obtained when $u^{(1)}$ and $\sigma^{(1)}$ are known exactly, i.e. at every edge integration point we have $u^{(1)} = u^W$ and $\sigma^{(1)} = \sigma^W$,

(B) the error that is due to interpolating solution displacements and stresses from nodes, resp. triangle integration points, to edge integration points,<fnref n="1" />

(C) the FE discretisation error is calculated by applying displacements on the specimen boundaries, then solving the BVP with the FEM as described above.

Each of these errors (see Fig. 3) adds up to the previous ones.

<figure><img src="image_5.png" /><figcaption>Figure 3: Location of known fields for the convergence study: (left) integration error, (centre) interpolation error, (right) FE discretisation error.</figcaption></figure>

Stresses are evaluated from the triangle to the edge integration points thanks to a piecewise linear interpolation method provided by the Scipy library [Virtanen et al., 2020]. The 2D interpolator is based on a Delaunay triangulation constructed with Qhull [Barber et al., 1996].

### 3.1.3 Comparison with J-integral

We wish to compare the Bueckner-Chen integral accuracy and convergence rate with the J-integral [Rice, 1968], which is commonly used for plane stress and plane strain problems.

For a straight crack, the J-integral writes:

$$ J = \oint (w \delta_{1j} - \sigma_{ij} u_{i,1}) n_j d\Gamma \qquad (20) $$

where $w$ is the strain energy density and $\delta$ the Kronecker symbol. This definition extends to an EDI:

$$ J = - \int_{A} (w \delta_{1j} - \sigma_{ij} u_{i,1}) q_j dA \qquad (21) $$

<footnote n="1">¹The linear interpolation of displacements from nodes provides the exact values at edge integration points.</footnote>
with $q$ the virtual crack extension field defined in Eq. 12. Both integrals are computed according to the numerical strategies given in Section 3.1.2. The displacement gradients are evaluated at the same locations as stresses and interpolated to the edge integration points if necessary.

For plane stress and mode I, Williams coefficient $A_1$ can be directly inferred from $J$:

$$A_1 = \sqrt{\frac{JE}{2\pi}} \qquad (22)$$

## 3.2 Line integral: numerical results and discussion

The following section presents the results obtained with the Bueckner-Chen and $J$-integrals when computing the three types of errors introduced in Section 3.1.2.

### 3.2.1 Exact integration

We first evaluate the numerical integration error (A) for the Bueckner-Chen integral with a Gaussian quadrature technique of the Scipy library, which is based itself on Fortran library QUADPACK. This routine adapts the number of integration points required, therefore performing the numerical equivalent of exact integration. Since the numerical fields $u^{(1)}$ and $\sigma^{(1)}$ are known exactly, considering the finite sample is equivalent to solving the extraction problem for an infinite medium. Tab. 2 summarises the relative errors on Williams coefficients obtained for orders $n = 1...4$ and these results demonstrate that exact integration is actually achieved.

<table><caption>Table 2: Relative error on Williams coefficients due to numerical integration with locally enriched Gaussian quadrature</caption><thead><tr><th>n</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>$\left|\frac{A_n^W - A_n}{A_n^W}\right|$</td><td>0</td><td>$1.76 \times 10^{-16}$</td><td>$9.71 \times 10^{-16}$</td><td>$6.58 \times 10^{-16}$</td></tr></tbody></table>

### 3.2.2 Evaluation of convergence indicators

In this section, we evaluate the convergence of the Bueckner-Chen LI with the three indicators defined in Section 3.1.2.

Fig. 4 provides the relative error on Williams coefficients obtained for orders $n = 1...4$ with the Bueckner-Chen LI and one integration point per edge (rectangle method, midpoint rule). The upper left graph shows the comparison between Buckner-Chen and $J$-integrals for $n = 1$. Several conclusions stem from these results. The error on Williams coefficients that is due to numerical integration (orange curves with plus markers) decreases like $h^2$ ($h$ being the element length) for all four orders. The results obtained with the $J$-integral for $n = 1$ are very similar. Since the numerical fields $u^{(1)}$ and $\sigma^{(1)}$ are assigned their exact analytical value, considering the finite geometry is equivalent to solving the extraction problem for an infinite medium. Hence, increasing the ratio $h/R$, which comes down to decreasing the number of edge integration points, is essentially the same as reducing the distance to the crack tip. The orange curves in Fig. 4 thus demonstrate the path-independence property of the Bueckner-Chen integral. The numerical scheme (B) adds the contributions of stress interpolation to the relative error. It provides the same decreasing rate of the relative error with regard to element length and a similar accuracy for orders $n = 1...4$. As emphasized by Melching and Breitbarth [2023], although the contribution of the higher-order terms ($n > 2$) in the Williams function vanish when $r \to 0$, the displacement and stress auxiliary functions in the Bueckner-Chen integral behave like $r^{-n/2}$, resp. $r^{-n/2-1}$, and thus diverge in the same area. This leads to numerical discrepancies when calculating higher-order terms too close to the singularity and might explains the unexpectedly low stress interpolation error for $n = 4$. The errors due to numerical integration, stress interpolation and FE discretisation add up to give the performance of the global scheme. The blue curves provide the global error, which decreases like $h$ and is, as expected, higher than the others. Calculating Williams coefficients on a very coarse mesh in an infinite domain comes down to reducing the distance to the crack tip, thereby decreasing the accuracy of the results as observed for the coarsest mesh (largest $h/R$). Interestingly, the $J$-integral provides FE discretisation errors higher than Bueckner-Chen by several orders of magnitude.
<figure><img src="image_6.png" /><figcaption>Figure 4: Relative error on Williams coefficients for an integration path with one integration point per edge (solid lines represent the error obtained with the Bueckner-Chen integral, dotted lines show the error on $A_1$ computed with the $J$-integral, colours refer to the three error types)</figcaption></figure>

Comparatively, Fig. 5 provides the relative error on Williams coefficients with:

* one integration point per edge for the error due to numerical integration, i.e. the same orange curves as in Fig. 4,

* exact integration over the contour for the errors due to stress interpolation and FE discretisation. In other words and contrary to Fig. 4, the numerical integration error does not contribute to the stress interpolation and FE discretisation errors.

The FE discretisation errors (blue curves) is higher than the integration errors (orange curves) for $n = 1...4$, therefore justifying the choice of an approximate evaluation of the integral with one integration point per edge, which provides reasonable accuracy at a reduced computational cost.

In our numerical method, stresses are interpolated from triangle to edge interpolation points with a piecewise linear operator based on a Delaunay triangulation (see Section 3.1.2). Therefore, when each edge counts several integration points, i.e. in the case of exact integration, stresses aren't constant along the edge. This approximation adds up to the uncertainty on the FE fields and might explain the discrepancies in the FE discretisation errors that we observe in Fig. 5 for $n = 1...3$.
<figure><img src="image_7.png" /><figcaption>Figure 5: Relative error on Williams coefficients computed with line integration. The integration error is computed with one integration point per edge, the stress interpolation and FE discretisation errors are computed with exact integration (solid lines represent the error obtained with the Bueckner-Chen integral, dotted lines show the error on $A_1$ computed with the $J$-integral, colours refer to the different error types)</figcaption></figure>
## 3.3 Equivalent Domain Integral: results and discussion

This section focusses on the new EDI formulation for the Bueckner-Chen integral introduced in Section 2.3. We study the convergence of this integral and compare its performance with the J-integral and LI formulations.

EDIs are evaluated on a discretised annulus-shaped domain $A^h$ defined by its outer and inner radii denoted $R$, resp. $R_{\text{int}}$. Two sets of parameters define $A^h$:

(1) $R = 10$ mm, $R_{\text{int}} = 0$ mm, includes the crack tip and therefore considers a worst-case scenario,

(2) $R = 10$ mm, $R_{\text{int}} = 5$ mm, excludes the singularity.

In both cases, $A^h$ is defined within the path $\Gamma^h$ used for LI, since the latter was defined by a radius $R = 10$ mm. Numerical integration is performed at the *triangle integration points*, thus bypassing the need for stress interpolation. For this reason, we subsequently investigate the numerical integration and FE discretisation errors only. The displacements are interpolated onto the integration points using the FE shape functions.

### 3.3.1 Evaluation of convergence indicators

Relative errors on Williams coefficients due to numerical integration are plotted in Fig. 6 for an equivalent domain that includes the crack tip ($R_{\text{int}} = 0$ mm) and compared with the method based on line integration. For $n = 1$, the error with the EDI decreases like $h$, at a slower rate than with the LI. The results provided by the Bueckner-Chen and J-integrals are similar in terms of convergence rate. The lower performance of the EDI compared to LI thus comes from the definition of the integration domain, which includes the singularity. For higher orders, the diverging supersingular fields at the crack tip dominate in the Bueckner-Chen integral, slowing down and eventually preventing convergence. Conversely, the error due to numerical integration decreases

<figure><img src="image_8.png" /><figcaption>Figure 6: Integration errors on Williams coefficients for LI and EDI with $R_{\text{int}} = 0$ mm (solid lines represent the error obtained with the Bueckner-Chen integral, dotted lines show the error on $A_1$ computed with the J-integral, the red curve with square markers is for the EDI and the orange one with cross markers for the LI)</figcaption></figure>

like $h^2$ when the integration domain excludes the singularity ($R_{\text{int}} = 5$ mm), the EDI and LI numerical schemes behaving similarly (Fig. 7).

Fig. 8 and Fig. 9 provide the errors computed with the global scheme for the mesh with $R_{\text{int}} = 0$ mm, resp. $R_{\text{int}} = 5$ mm. Although it is difficult to characterise the convergence in the first case,
<figure><img src="image_9.png" /><figcaption>Figure 7: Integration errors on Williams coefficients for line and domain integrals with $R_{int} = 5$ mm (solid lines represent the error obtained with the Bueckner-Chen integral, dotted lines show the error on $A_1$ computed with the $J$-integral, the red curve with square markers is for the EDI and the orange one with cross markers for the LI)</figcaption></figure>

the latter proves to be more accurate and converges faster than LI. In both cases, the Buckner-Chen integral is more accurate than $J$. Regarding the first integration domain ($R_{int} = 0$ mm), the same numerical discrepancies as in Fig. 6 appear at higher orders. When the singularity is excluded from the integrated fields however (Fig. 9), the error due to FE discretisation decreases like $h$ for $n = 1...4$. Additionally, Bueckner-Chen is always more accurate than the $J$-integral by several orders of magnitude, as demonstrated by the error levels on $A_1$.

### 3.3.2 Comparison of Stress Intensity Factors

This last section summarises the previous results and compares the four methods: the $J$- and Bueckner-Chen LIs and EDIs. To this end, we compute the values obtained for $K_I$, which are directly inferred from $A_1$ using the following relation:

$$K_I = A_1 \sqrt{2\pi} \qquad (23)$$

The pseudo-boundary conditions described for the FE discretisation error are used (see Section 3.1.2, (C)). All numerical schemes behave similarly in terms of convergence rate, since the relative error on $K_I$ decreases like $h$ with all four methods. However, the Bueckner-Chen integral proves to be more accurate than the $J$-integral regardless of the integration method (LI or EDI), even when the integration domain includes the singularity. Moreover, the error on $K_I$ is even lower when the integration domain (line or surface) excludes the crack tip. Therefore and from a strictly numerical point of view, the Bueckner-Chen integral is an accurate, robust and fast converging method to extract the SIF and Williams higher-order terms. The line and domain integrals provide similar results but the latter does not require stress interpolation.

## 4 Conclusion

This work aimed to prove the numerical convergence of the Bueckner-Chen LI and EDI to retrieve Williams coefficients. For this reason, we uncoupled the errors due to numerical integration, stress
<figure><img src="image_10.png" /><figcaption>Figure 8: FE discretisation errors on Williams coefficients for line and domain integrals with $R_{\text{int}} = 0$ mm (solid lines represent the error obtained with the Bueckner-Chen integral, dotted lines show the error on $A_1$ computed with the *J*-integral, the dark blue curve with big markers is for the EDI and the light blue one with small markers for the LI)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: FE discretisation errors on Williams coefficients for line and domain integrals with $R_{\text{int}} = 5$ mm (solid lines represent the error obtained with the Bueckner-Chen integral, dotted lines show the error on $A_1$ computed with the *J*-integral, the dark blue curve with big markers is for the EDI and the light blue one with small markers for the LI)</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10: Comparison of relative errors on $K_I$ computed with line and domain integrals (solid lines represent the error obtained with the Bueckner-Chen integral, dotted lines show the error on $K_I$ computed with the $J$-integral, colours refer to different integration domains)</figcaption></figure>

interpolation and FE discretisation and showed that the latter is always greater by several orders of magnitude. Regarding the EDI, the convergence order of the numerical scheme naturally depends on the integration domain. When it includes the crack tip, the convergence order of the extraction method based on the Bueckner-Chen integral depends on the Williams order. Nevertheless, with an adequate domain (i.e. one that excludes the crack tip) the convergence order is about 2. The comparison with the values obtained for the SIF proves the efficiency and the robustness of the new method compared to the $J$-integral. The results obtained with LI and EDI are similar, although the latter bypasses the need for stress interpolation onto the integration contour.

The next step consists in applying the procedure based on the Bueckner-Chen EDI to experimental data. Adding noise will inevitably degrade convergence, especially as the Bueckner-Chen integral depends on experimental stresses, usually obtained by computing the displacement gradients in the framework of LEFM. In this context, quantifying the orthogonality loss would be an interesting lead to explore: given the approximation due to the numerical integration, does the method conserve the orthogonality property of Williams eigenfunctions? Can we adjust the virtual crack extension field $q$ to correct this orthogonality loss and offset noise? These questions should be addressed in a future work to investigate further the quality of the method based on the Bueckner-Chen integral for the extraction of Williams coefficients.

## Acknowledgments

This research was funded, in whole, by the French national research agency (ANR) through grant ANR-20-CE08-0017.

## Declarations

A CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission, in accordance with the grant's open access conditions.

The code used to perform the simulations presented in this document is available online at: https://github.com/hdandin/adventure.git.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
References

M. Attigui and C. Petit. Mixed-mode separation in dynamic fracture mechanics : New path independent integrals. *International Journal of Fracture*, 84(1):19–36, 1997.

G. Bahaj Filali, M. Coret, A. Leygue, and J. Réthoré. Interpretable crack features for the representation of kinematic fields in the case of fatigue overloads. *International Journal of Fracture*, 249(1):15, Jan. 2025. doi: 10.1007/s10704-024-00830-2.

C. B. Barber, D. P. Dobkin, and H. Huhdanpaa. The quickhull algorithm for convex hulls. *ACM Trans. Math. Softw.*, 22(4):469–483, Dec. 1996. doi: 10.1145/235815.235821.

J. Bellecave, S. Pommier, Y. Nadot, J. Meriaux, and J. Araújo. T-stress based short crack growth model for fretting fatigue. *Tribology International*, 76:23–34, 2014.

H. F. Bueckner. Field singularities and related integral representations. In G. C. Sih, editor, *Methods of Analysis and Solutions of Crack Problems*, pages 239–314. Springer Netherlands, Dordrecht, 1973. doi: 10.1007/978-94-017-2260-5_5.

H. Bui. Associated path independent J-integral for separating mixed modes. *Journal of Mechanics and Physics of Solids*, 31:439–448, 1983.

H. Chen and R. Shield. Conservation laws in elasticity of the j-integral type. *Journal of Applied Mathematics and Physics*, 28, 1977.

Y. Chen. New path independent integrals in linear elastic fracture mechanics. *Engineering Fracture Mechanics*, 22(4):673–686, Jan. 1985. doi: 10.1016/0013-7944(85)90131-6.

B. Cotterell. Notes on the paths and stability of cracks. *International Journal of Fracture Mechanics*, 2(3):526–533, Sept. 1966. doi: 10.1007/BF00193691.

P. Destuynder, M. Djaoua, and S. Lescure. Quelques remarques sur la mécanique de la rupture élastique. *Journal de Mécanique Théorique et Appliquée*, 2(1), 1983.

E. Eid, R. Seghir, and J. Réthoré. Crack branching at low tip speeds: spilling the t. *Journal of Theoretical Computational and Applied Mechanics*, 2023.

C. Geuzaine and J.-F. Remacle. Gmsh: A three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. *International Journal for Numerical Methods in Engineering*, 79(11):1309–1331, May 2009. doi: 10.1002/nme.2579.

R. Hamam, F. Hild, and S. Roux. Stress Intensity Factor Gauging by Digital Image Correlation: Application in Cyclic Fatigue. *Strain*, 43(3):181–192, Aug. 2007. doi: 10.1111/j.1475-1305.2007.00345.x.

M. Kanninen and C. Popelar. *Advanced fracture mechanics*. Oxford University Press, New York, 1985.

D.-J. Kim, J. Pereira, and C. Duarte. Analysis of three-dimensional fracture mechanics problems: A two-scale approach using coarse-generalized fem meshes. *International Journal for numerical methods in engineering*, 81(3):335–365, 2010.

S. McNeill, W. Peters, and M. Sutton. Estimation of stress intensity factor by digital image correlation. *Engineering Fracture Mechanics*, 28(1):101–112, 1987.

D. Melching and E. Breitbarth. Advanced crack tip field characterization using conjugate work integrals. *International Journal of Fatigue*, 169:107501, Apr. 2023. doi: 10.1016/j.ijfatigue.2023.107501.

B. Moran and C. Shih. Crack tip and associated domain integral from momentum and energy balance. *Engineering Fracture Mechanics*, 27(6):615–642, 1987.

H. T. Nguyen, M. Pathirage, G. Cusatis, and Z. P. Bažant. Gap test of crack-parallel stress effect on quasibrittle fracture and its consequences. *Journal of Applied Mechanics*, 87(7):071012, 2020.
J. Passieux, A. Gravouil, J. Réthoré, and M. Baietto. Direct estimation of generalized stress intensity factors using a three-scale concurrent multigrid x-fem. *International Journal for Numerical Methods in Engineering*, 85(13):1648-1666, 2010. doi: 10.1002/nme.3037.

J. Réthoré, S. Roux, and F. Hild. Optimal and noise-robust extraction of fracture mechanics parameters from kinematic measurements. *Engineering Fracture Mechanics*, 78(9):1827-1845, 2011.

J. R. Rice. A path independent integral and the approximate analysis of strain concentration by notches and cracks. *Journal of Applied Mechanics*, 35(2):379-386, June 1968. doi: 10.1115/1.3601206.

J. R. Rice. Some remarks on elastic crack-tip stress fields. *International Journal of Solids and Structures*, 8(6):751-758, June 1972. doi: 10.1016/0020-7683(72)90040-6.

S. Roux and F. Hild. Stress intensity factor measurement from digital image correlation: post-processing and integrated approaches. *International Journal of Fracture*, 140(1-4):141-157, 2006.

S. Roux, J. Réthoré, and F. Hild. Digital image correlation and fracture: an advanced technique for estimating stress intensity factors of 2d and 3d cracks. *Journal of Physics D: Applied Physics*, 42(214004), 2009. doi: 10.1088/0022-3727/42/21/214004.

J. Réthoré and R. Estevez. Identification of a cohesive zone model from digital images at the micron-scale. *Journal of the Mechanics and Physics of Solids*, 61(6):1407-1420, 2013. ISSN 0022-5096. doi: https://doi.org/10.1016/j.jmps.2013.01.011.

P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, et al. SciPy 1.0: Fundamental algorithms for scientific computing in Python. *Nature Methods*, 17(3):261-272, Mar. 2020. doi: 10.1038/s41592-019-0686-2.

M. L. Williams. On the Stress Distribution at the Base of a Stationary Crack. *Journal of Applied Mechanics*, 24(1):109-114, 1957. doi: 10.1115/1.4011454.

Q. Xiao and B. Karihaloo. Implementation of hybrid crack element on a general finite element mesh and in combination with xfem. *Computer methods in applied mechanics and engineering*, 196(13):1864-1873, 2007.