# Stress intensity factor determination along a kinked crack path by DIC analyses

Thomas Corre, François Hild, Veronique Lazarus

► To cite this version:

Thomas Corre, François Hild, Veronique Lazarus. Stress intensity factor determination along a kinked crack path by DIC analyses. International Journal of Fracture, 2025, 249, pp.42. 10.1007/s10704-025-00862-2. hal-05085982

HAL Id: hal-05085982

https://hal.science/hal-05085982v1

Submitted on 27 May 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Stress intensity factor determination along a kinked crack path by DIC analyses

T. Corre¹*, F. Hild² and V. Lazarus³

¹Nantes Université, Ecole Centrale Nantes, CNRS, GeM, UMR 6183, 44000 Nantes, France.

²Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS–Laboratoire de Mécanique Paris-Saclay, 91190 Gif-sur-Yvette, France.

³IMSIA, CNRS, EDF, CEA, ENSTA Paris, Institut Polytechnique de Paris, 91120 Palaiseau cedex, France.

*Corresponding author(s). E-mail(s): thomas.corre@ec-nantes.fr;

## Abstract

Sharp kinks may be observed under shear loading or in materials containing weak directions, such as those produced by additive manufacturing. A better understanding of the fracture of these materials, both theoretically and experimentally, is required to deploy them in structural applications. This study focuses on the measurement of stress intensity factors (SIFs) around a sharp kink using digital image correlation (DIC). The performances of two DIC-based techniques, namely, integrated-DIC and post-processing of DIC-measured displacement fields, are assessed on a benchmark test using fused deposit modeling capabilities, and are compared to a reference finite element solution. It is shown that Williams’ expansion remains valid on a large enough region around the crack to extract reliable SIFs even very close to the crack kink. Both techniques are very trustworthy, provided the SIF identification zone is carefully defined to exclude the kink zone of influence.

**Keywords:** Angular crack; Crack tip position; Fused Deposition Modeling (FDM); Interaction integral; Linear Elastic Fracture Mechanics (LEFM); Williams’ series
1 Introduction

Additively Manufactured materials and crystals may have anisotropic fracture prop-
erties leading to unexpected crack propagation [1, 2]. Expecting an increasing use of
such materials (e.g., aeronautics/space industry [3]), it becomes essential to extend lin-
ear elastic fracture mechanics tools from isotropic to anisotropic cases. Several reasons
may be invoked. On the one hand, prerequisite to their use in sensitive components is
to assess in a damage tolerance framework that a potential preexisting crack will not
lead to catastrophic failure. On the other hand, a better knowledge of the effect of
anisotropy on crack propagation will open the way to optimal design of the printing
strategy that triggers such anisotropy [4].

While Linear Elastic Fracture Mechanics (LEFM) has shown its efficiency to deal
with crack propagation in isotropic materials [5], its extension to anisotropic settings
still requires additional developments, notably to be validated by close comparisons
with experiments. Strong anisotropy with preferential weak fracture planes may lead
to crack kinking [6]. It is thus essential to be able to deal with kinked crack paths, in
particular with the change of Stress Intensity Factors (SIFs) or Energy Release Rate
(ERR) across the angular point. From a theoretical point of view, their values have
been shown to be discontinuous from one side to the other [7]. In the case of elastic
isotropy, their changes crossing the kink position is obtained using Amestoy-Leblond's
formula [8].

While Linear Elastic Fracture Mechanics (LEFM) has shown its efficiency to deal
with crack propagation in isotropic materials [5], its extension to anisotropic settings
still requires additional developments, notably to be validated by close comparisons
with experiments. Strong anisotropy with preferential weak fracture planes may lead
to crack kinking [6]. It is thus essential to be able to deal with kinked crack paths, in
particular with the change of Stress Intensity Factors (SIFs) or Energy Release Rate
(ERR) across the angular point. From a theoretical point of view, their values have
been shown to be discontinuous from one side to the other [7]. In the case of elastic
isotropy, their changes crossing the kink position is obtained using Amestoy-Leblond's
formula [8].

From an experimental point of view, the previous fracture mechanics parameters
are assessed with Digital Image Correlation (DIC) measurements [9]. Very early on,
it was shown that DIC could be used to extract SIFs in experiments under static [10]
and dynamic [11] loading conditions. Other parameters such as crack tip opening angles [12] or crack opening displacements [13–15] were estimated. One key quantity for assessing fracture mechanics parameters is the crack tip location, which has a direct consequence on the amplitude of SIFs [16]. This problem was, for instance, addressed by using the first supersingular field [9]. This criterion will be used herein.

SIF measurement techniques using DIC measurements are commonly based on Williams' series [10, 16]. However, the theoretical displacement field they provide relies on the assumption of a straight crack in an infinite medium, which thus is valid only close enough to the crack tip in comparison to the distance to the boundaries or to any wavy part of the crack. Its zone of validity close to a kink needs to be assessed, especially as a kink defines two angular boundaries and is therefore also a stress concentrator [17] that may perturb Williams' expansion around the crack tip. This observation questions until which distance of the tip reliable SIF measurements may be accurately obtained by DIC if any.

To answer this question, fused deposit modeling (FDM) polycarbonate samples [6] were considered as benchmarks, as they had the advantages of (i) having isotropic linear elastic properties (hence consistent with the framework of isotropic LEFM), while (ii) being endowed with anisotropic fracture properties enabling the crack to kink toward one of the weak printing directions. In a first study [6], a classical Pacman-shaped extraction zone with a straight mouth was considered, thereby preventing the authors to obtain reliable SIFs when the kink was in the identification zone. The present study seeks to go much closer to the kink by analyzing the influence of the identification parameters and ensuring theoretical consistency of the obtained results.
It focuses on two DIC-based techniques, namely, integrated DIC and post-processing of DIC-measured displacement fields.

Section 2 of the paper recalls the definitions of SIFs and their expressions across a kink according to Linear Elastic Fracture Mechanics. The experiment and the DIC methods utilized herein are described in Section 3. This section also presents the finite element methods used to provide reference values for SIFs along the whole experiment. This solution near the kink has been carefully validated against closed-form solutions (Appendix A). Section 4 discusses the results about the influence of the main parameters of the SIF analysis, namely, shape and size of the identification zone and modeling choice of the crack path. The measured SIFs with different angular openings of the identification zone are compared for the two techniques with the reference finite element solution to identify good practice in a reasoned way.

## 2 Linear Elastic Fracture Mechanics

The aim of this section is (i) to introduce the definition of the key quantities describing stress fields around a crack tip, namely, the SIFs $K_I$, $K_{II}$ and the T-stress together with the underlying assumptions, and (ii) to describe their discontinuous change across a kink. The analyses are restricted to 2D problems as DIC analyses were conducted to measure 2D displacements of the sample surface. Further, only surface measurements were performed. No information was available on the crack front shape at the vertex, or on displacement fields in the bulk [18, 19]. Consequently, 3D effects, such as out-of-plane displacements correlated to Poisson effects [20] or ill-defined SIFs at that point [21, 22] are not addressed herein. A plane strain assumption is adopted
throughout the paper. While using a plane stress assumption would slightly alter the SIF values, it would not affect the main conclusions since the same assumption is consistently applied for the Williams' fields and the finite element simulations.

## 2.1 Definition of SIFs and T-stress

Let us consider a crack (Figure 1(a)) in an isotropic and homogeneous elastic medium subjected to some loading (i.e., prescribed displacements and/or forces).

<figure><img src="image_2.png" /><figcaption>Fig. 1: Zoom around the tip of a crack embedded in a linear elastic medium and associated notations. (a) Initial crack. (b) Kinked extension of length s in direction $\varphi$.</figcaption></figure>

In the vicinity of its tip ($r \to 0$), the stress field $\sigma$ takes the following form [23] in the Cartesian frame defined in Figure 1(a)

$$ \sigma_{ij}(r, \theta) = \frac{K_I}{\sqrt{2\pi r}} f_{ij}^I(\theta) + \frac{K_{II}}{\sqrt{2\pi r}} f_{ij}^{II}(\theta) + T\delta_{ix}\delta_{jx} + O(\sqrt{r}) \quad (1) $$

where $r, \theta$ denote the polar coordinates defined with respect to the crack tip and the crack surface defining the $x$ axis. The fields $f_{ij}^p(\theta)$ are universal (in the sense that they do not depend on the geometry of the body and the boundary conditions) and depend on $\theta$ only [5, 24].
The prefactors $K_I$ and $K_{II}$ quantifying the intensity of the stress singularity in the vicinity of the crack tip are the so-called SIFs. Note that in the 2D plane case, there is no mode III contribution ($K_{III} = 0$). The first non singular term is referred to as the T-stress [25].

## 2.2 SIF discontinuity across a kink

Let us assume that the loading of the medium is kept constant between the two configurations of Figure 1(a) (just before kink) and Figure 1(b) (same configurations with a kinked portion of length $s$ in direction $\varphi$). Let $K_p^0$ and $T^0$ denote the SIF and T-stress before the kink, and $K_p(s, \varphi)$ the SIF in the second configuration, taking the new tip and crack direction as reference for the polar coordinate system (Figure 1(b)). Then, in the limit $s \to 0$ [7]

$$K_p(s, \varphi) = F_{pq}(\varphi)K_q^0 + \sqrt{s}T^0G_p(\varphi) + O(s), \quad (2)$$

where $F_{pq}(\varphi)$ and $G_p(\varphi)$ are universal fields depending on the sole value of $\varphi$. This result means that knowing the values of $K_q^0$ and $T^0$ before the kink, $K_p(s, \varphi)$ is obtained after the kink with no need for solving the entire elasticity problem. In the case of isotropic elasticity, the fields were provided by Amestoy and Leblond [8] for any value of $\varphi$, either in the form of some integral equations coupled with Anderson's formula (Equations (34)$_1$, (35)$_1$, (36), (39) therein) or of a series (Equation (66) therein).

Let us introduce the limit $K_p^*(\varphi)$ for $s \to 0$

$$K_p^*(\varphi) = \lim_{s \to 0} K_p(s, \varphi) \quad (3)$$
and using Equation (2)

$$K_p^* = F_{pq}(\varphi)K_q \quad (4)$$

Despite the continuity of the displacement and stress fields when $s \to 0$, the SIFs are discontinuous crossing the kink, that is $K_p^*(\varphi) \neq K_p$ and their derivative are $1/\sqrt{s}$-singular for $s \to 0$. Conversely, if the crack path is smooth, then $\varphi = 0$, $F_{pq} = \delta_{pq}$ and $G_p = 0$, such that the change of $K_p(s)$ and its derivative are continuous. The term $K_p(s, \varphi = 0)$ is then of order $s$ for $s \to 0$, and depends on the curvature of the path and on the whole geometry of the elasticity problem [7].

## 3 Analysis Methods

This section details the experimental and finite element methods used to analyze the evaluation of the SIFs close to a kink. The experiment providing a practical case is described, followed by the DIC techniques employed to measure SIFs during this test. Last, the construction of a reference solution using finite elements is discussed.

### 3.1 Experiment

This study is based on a new analysis of one of the crack propagation experiments published by Corre et al. [6]. In that study, kinked crack propagation was observed in a polycarbonate sample produced by fused deposition modeling (FDM). As a full description of the experimental procedure is available in that reference, only its main features are recalled hereafter.

Compact tension (CT) specimens were fabricated to monitor crack propagation (Figure 2(a)). Samples (50 × 48 × 10 mm³) were produced with a Raise 3D N2 plus
124 FDM machine along the width direction, aiming for a 100% filling. The 0.3 mm notch 125 was directly printed as part of the geometry. The outer shell was not printed at the 126 notch root to promote crack propagation within the bulk of the material. No specific 127 pre-cracking procedure was applied. The filaments were deposited in straight lines 128 on each 0.15 mm layer, with each of them oriented orthogonally with respect to the 129 previous one. The raster angle was defined with respect to the CT notch (smallest 130 deposit angle of the two types of layers). The filament width *d* was about 0.2 mm.

<figure><img src="image_3.png" /><figcaption>**Fig. 2:** (a) Crack propagation in the studied sample (last image of the test; raster angle of 25°). The inset shows a zoom on the kink. (b) Loading curve of the test. Each circle indicates a picture acquired for DIC measurements.</figcaption></figure>

131 A quasi-static fracture experiment was then performed on a universal tensile 132 machine with a prescribed stroke rate of 0.5 mm/min. With this experimental pro- 133 cedure, crack propagation always followed one of the deposit directions [6, 26]. An
example of a sample with a 25° raster angle is shown in Figure 2. After a short straight propagation, the crack kinked to follow the deposit direction. Thus, fracture occurred at the interface of two filaments in every other layer, and perpendicular to the filaments in the others (see broken filaments in Figure 2(a)). This surface acted as a weak plane in the material and prescribed a kink angle equal to that of the raster pattern. The propagation remained stable, as illustrated by the load-displacement curve of Figure 2(b). The material just ahead of the notch did not exhibit the same behavior because of the locally complex printing trajectory. However, the crack path in that area, before the kink, was close to a straight line. This material exhibited a strong anisotropic fracture pattern. However, it was shown that its elastic behavior remained close to isotropy in the printing plane with a Young's modulus $E = 1980$ MPa, and Poisson's ratio $\nu = 0.3$ [6].

The experiment presented in Figure 2 is the case study for this paper. In the following, the kink angle will be equal to 25°. To perform DIC analyses, the test was monitored with a Nikon D300s camera acquiring grayscale images of definition 4288 × 2848 pixels every 4 s around the kinking event, and 6 s otherwise. The camera was mounted with an AF-S Micro NIKKOR 60 mm lens, resulting in a physical pixel size of 19.3 µm.

## 3.2 SIF measurements with digital image correlation

Two different routes were followed to assess SIFs via DIC. First, FE-based DIC was run and subsequently post-processed to extract fracture mechanics parameters. Second, integrated DIC, which merges the two previous steps into a single one, was also carried
out. In both cases, the reference image was always the very first one, namely, in the unloaded configuration.

### 3.2.1 FE-based DIC

Finite-element based DIC was used herein [27–29]. The sum of squared differences $\rho$

$$ \rho = \sum_{\text{ROI}} \rho^2(\underline{x}) \qquad (5) $$

is minimized over the whole region of interest (ROI), where $\rho$ denotes the pixel-wise gray level residuals

$$ \rho(\underline{x}) = f(\underline{x}) - g(\underline{x} + \underline{u}^{Q4}(\underline{x})) \qquad (6) $$

with $f$ the image of the reference configuration, $g$ that of the deformed configuration.

The sought displacement field $\underline{u}_{Q4}$ is written as

$$ \underline{u}_{Q4} = \sum_{i} v_i \underline{N}_i^{Q4}(\underline{x}) \qquad (7) $$

where $v_i$ refers to the nodal displacement associated with shape function $\underline{N}_i^{Q4}$ of a Q4 discretization (i.e., 4-noded quadrilaterals with bilinear shape functions).

An iterative Gauss-Newton scheme was implemented, which leads to linear systems

$$ [\underline{H}^{Q4}]\{\delta\underline{v}\} = \{\underline{h}^{Q4}\} \qquad (8) $$

where $[\underline{H}^{Q4}]$ is the DIC Hessian matrix, $\{\underline{h}^{Q4}\}$ the right hand side column vector, and $\{\delta\underline{v}\}$ the corrections to the current estimate of the nodal displacement gathered in the column vector $\{\underline{v}\}$.

### 3.2.2 SIF extraction domain

The first method used in this study to estimate SIFs consists in post-processing measured displacement fields [10, 16, 30] and will only briefly be recalled herein. In the
vicinity of its tip ($r \to 0$), the displacement field $\underline{u}$ takes the following form in polar coordinates [31]

$$ \underline{u}(r, \theta) = \sum_{n=n_{\min}}^{\infty} a_I^n r^{\frac{n}{2}} \underline{\psi}_I^n(\theta) + a_{II}^n r^{\frac{n}{2}} \underline{\psi}_{II}^n(\theta), \quad (9) $$

where $n_{\min}$ is generally equal to 1. The fields $\underline{\psi}_I^1(\theta)^1$ and $\underline{\psi}_{II}^1(\theta)$ are universal; their expressions can be found in Refs. [9, 30]. The extraction of Williams' amplitudes was carried out on a tailored ROI, referred to as Pacman in the present study (Figure 3(a)), via least squares minimization with respect to the measured displacement fields.

The Pacman domain is a disk of external diameter $2R_{\max}$ centered about the crack tip and deprived from an inner disk of diameter $2R_{\min}$. A section of thickness $R_{\min}$ around the mouth was also removed. This domain could also be opened by an angle $\alpha$, meaning that the corresponding angular sector was removed from the domain. The proposed crack path was defined as either two straight lines that mark a sharp kink or a curved line that smooths the path about the kink. The test was analyzed using direct or reverse extraction strategies along the crack path (Figure 3(a)), namely, starting from either the first or last image of the full acquired sequence. Such choice was dictated by the presence of a kink that makes both options not necessarily equivalent, as will be shown later on.

As the Pacman ROI does not cover the sample surface, the Williams' series (9) are truncated. In the present case, $n_{\min} = -3$ and $n_{\max} = 7$ were selected. The mode I field corresponding to $n = -1$ was used to detect the position of the crack tip at each step of the test [9]. An iterative procedure was followed to cancel out the amplitude of the first mode I supersingular field. The radius $R_{\max}$, the opening angle $\alpha$, the shape of the crack path, and the analysis strategies are the four parameters whose
influence on SIF estimates is assessed in this study. The minimum order $n_{min} = -3$ is standard practice [9]. The maximum order $n_{max}$ depends on the size of the pacman. A convergence analysis was performed on the crack tip position and SIFs to select its value [32].

The software UFreckels [33] was used for the analysis of the test with the post-processing technique. DIC analyses were performed with 16 × 16 px Q4 elements. The radius $R_{max}$ varied between 100 px and 300 px, and $R_{min}$ was set to 32 px. Figure 3(b,c) show examples of the extraction domain with Q4-DIC elements (the extraction domain is defined at the node level). This post-processing method is denoted W-DIC in the following.

### 3.2.3 Integrated DIC

Integrated DIC consists in performing the previous two steps (i.e., DIC analysis and subsequent extraction of fracture mechanics parameters) at the same time, namely, the selected kinematic basis is directly Williams' series. Two routes can be followed. Either a standalone DIC code is implemented, as was the case in earlier applications [16], or FE meshes can still be used but tailored to Williams' series [34] (i.e., pacman-like ROIs, see Figure 3(d,e)). Williams' displacement fields are then gathered in a rectangular matrix [$\Psi_W$] that enables each nodal displacement to be related to the considered Williams' field

$$ \{v\} = [\Psi_W]\{a\} \qquad (10) $$

where {$a$} is the column vector gathering the sought amplitudes. Such strategy enables for the implementation of any reduced basis in a non intrusive way.
<figure><img src="image_4.png" /><figcaption>Fig. 3: (a) Parameters describing the identification ROI of SIFs. Two extraction strategies are implemented to estimate the fracture mechanics parameters. The first one (direct) follows the crack path as experimentally. The second one (reverse) starts from the last acquired image to the first one. Examples of pacmen with Rmax = 250 px. (b,c) Post-processing cases (α = 0° and 25°, before the kink). (d,e) Integrated DIC ROIs (α = 0° and 25°, after the kink).</figcaption></figure>

The integrated DIC Hessian matrix then reads

$$[\mathbf{H}_W] = [\mathbf{\Psi}_W]^\top [\mathbf{H}^{T3}] [\mathbf{\Psi}_W] \qquad (11)$$
and right-hand side vector

$$ \{h_W\} = [\Psi_W]^T \{h^{T3}\} \quad (12) $$

where $[H^{T3}]$ is the DIC Hessian associated with an FE discretization of thetmeman with 3-noded triangular (T3) elements, which allows for more faithful descriptions of pccman-like ROIs (Figure 3(d,e)).

As for the post-processing route, the crack tip position is iteratively estimated by canceling out the amplitude of the first supersingular field. Therefore, there are two iterative processes (for DIC and for crack tip detection) that are nested. This imple- mentation implies that the integrated DIC mesh (pacman) is iteratively moved, and integrated DIC is performed on the new ROI at each change of the crack tip position. Integrated DIC was implemented in the Correli 3.0 framework [35]. C++ kernels deal- ing with all the heavy computations (i.e., node to pixel projector, computation of the T3-DIC Hessian matrix $[H^{T3}]$, right-hand side vector $\{h^{T3}\}$, and gray level residuals $\rho$) are based upon T3 discretizations. These kernels are then called by Matlab via MEX files, and the minimization is still performed within the Matlab environment such that any choice of reduced basis can be made. The same set of parameters was investigated as in the previous case. Integrated DIC will be designated as I-DIC in the sequel.

## 3.3 Finite Element computations as reference

To evaluate the performances of DIC to measure SIFs close to kinks, reference solutions were obtained by running 2D Finite Element (FE) computations for each picture of the test. As accurately modeling the boundary conditions at the pins of CT tests is challenging [36], the FE model was restricted to the ROI of Q4-DIC and the Dirichlet
boundary conditions were the nodal displacements directly measured during the test.

A plane strain model of the ROI was created, to be consistent with Williams' fields, assuming linear elasticity and using the experimentally determined Poisson's ratio (Section 3.1).

Figure 4 shows the ROI on the sample where DIC was run and the corresponding mesh for FE analyses. The position of the nodes on the outer layer of the model matched those of Q4-DIC mesh. Then, for each recorded image, the reference value of $K_I$ and $K_{II}$ was obtained through the following procedure:

• the crack geometry and in particular the crack tip position, the nodal displacements at the boundary of the ROI were extracted from DIC analyses;

• the displacement field was computed via FE analyses of elasticity problems with the measured displacements prescribed on the outer boundary of the model;

• the SIFs were extracted from this solution using the G_theta method [37, 38].

The use of finite elements enables for the implementation of a sufficiently refined mesh to apply the G_theta method to a kinked extension of arbitrarily small length while excluding the kink from the computation zone. As a result, the calculation of SIFs using finite elements is not affected by the presence of the kink, unlike the values extracted from the less spatially resolved DIC fields. The validity of this numerical approach near the kink is assessed by comparing the computed SIFs for various kink angles with the theoretical predictions given in Equation (2). The method used to compute these reference SIF values, along with its validation, is detailed in Appendix A.
<figure><img src="image_5.png" /><figcaption>**Fig. 4:** (a) Region of interest where Q4-DIC was performed with the displacement field expressed in pixels (norm). (b) Mesh for the computation of the reference value. Example with the crack tip located after the kink.</figcaption></figure>

The reference values were obtained once and for all with one set of DIC and Pacman parameters to get the nodal displacements and the crack geometry. While the nodal displacements along the boundary of the mesh (used as Dirichlet boundary conditions) are insensitive to these parameters, the successive positions of the crack tip depend on them. This is especially the case when the Pacman parameters were clearly not adapted (when $R_{max}$ was too small for instance). That is the reason why the reference value was set with one set of parameters that proved to be well suited, namely, I-DIC and direct analysis with a kinked path and $R_{max}$ = 200 px. The locations of crack tips were nearly identical for similar parameters, but the instant when the crack kinked
277 could not be accurately determined due to experimental inaccuracies within one or
278 two images.

279 The post-processing procedure used in W-DIC may also be applied to finite-element
280 displacement fields, using them as synthetic fields. In that case, the position of the
281 crack tip was set and only the amplitudes of Williams' series were identified. The
282 use of FE fields avoids experimental constraints related to DIC, such as measurement
283 uncertainty or minimum pattern size, to concentrate on the performance of the iden-
284 tification of SIFs close to the kink. In practice, the minimum size of the element was
285 set to 2 px around the kink. This size was similar to that of the smallest element used
286 here in the I-DIC method and allowed a relevant comparison to be maintained with
287 practical DIC measurements. This type of SIF identification on synthetic fields will
288 be denoted as W-FE, by analogy with the post-processing route of DIC fields.

## 4 Results

289 This section presents the influence of the DIC technique and extraction domain on
290 SIF estimates. The focus is only put on the salient features of each technique, then
291 followed by a comparison with the reference solution. The analysis is based on four
292 methods:
293

294 • *W-DIC*: SIFs were obtained by post-processing DIC measured displacement fields
295 using Williams' series.

296 • *I-DIC*: Williams' fields were the kinematic basis in integrated DIC, which directly
297 measured SIFs.
298 • Reference: Extraction of SIFs using the G_theta method on finite element computations, performed at each step (i.e., image) with boundary conditions and crack tip locations obtained from the DIC analyses.

300 • W-FE: Application of the same W-DIC post-processing method to synthetic finite element displacement fields using DIC-based boundary conditions.

## 4.1 Shape of identification zone: Study on synthetic fields

Equation (1) assumes a straight crack in an infinite medium, which is unsuitable when the kink is too close to the identification zone, or even within it. Using one synthetic field (W-FE), Figure 5 shows the norm of the displacement difference between the FE simulation and its projection onto Williams' basis. The area in which the displacements are very different from Williams' solution is very small and located around the kink (when $-3 \le n \le 7$). In Figure 5(a,b), the largest differences are around the mouth of the crack and before the kink. These sub-figures illustrate the two possibilities used to remove elements around the crack from the identification zone, namely, a straight zone in the extension of the kinked part of the crack, and removing the elements at a fixed distance of the mouth along the crack. There is a small difference in the quality of the fit between these two ways of defining thehti (only the second one will be used in the remainder of this section). In both cases, the discrepancies were located in the angular sector formed by the initial crack and the kinked branch. That is the reason why removing elements from thehti in that area significantly reduced the fit discrepancies, as shown in Figure 5(c,d) for ahti opening angle of 25°.
<figure><img src="image_6.png" /><figcaption>**Fig. 5:** Norm of displacement differences (expressed in px) for each node between the synthetic FE simulation and its projection onto Williams' basis. Two types of Pacman opening with two different angles were studied. Straight Pacman opening (left column) and opening following the crack path (right column). Axes in px; the dashed line marks the crack path. Note that the scales of top and bottom panels are different.</figcaption></figure>

Figure 6 details the influence of the Pacman opening angle on the relevance of the calibrated Williams' fields with respect to all analyzed synthetic displacement fields.
As shown in Figure 6(a), the RMS error of the identified displacement is significantly reduced for angles $\alpha$ of about 20°. This result is presented for the first image after the kink, but similar results were obtained for all images for which the distance between the kink and the crack tip was less than $R_{max}$. For a kink angle of 25°, an angle $\alpha$ of 25° removes the kink from the identification zone (Figure 3(a)). Removing the kink and the wake area from the Pacman was sufficient to improve significantly the identified displacement field.

<figure><img src="image_7.png" /><figcaption>Fig. 6: (a) RMS displacement error for various sizes and Pacman opening angles (first step after the kink in the W-FE method). (b) RMS displacement error along the propagation path with various angles $\alpha$ ($R_{max}$ = 250 px).</figcaption></figure>

As illustrated in Figure 6(b), large RMS errors are observed close to the kink for angles $\alpha$ less than 25°, until the identification zone was far enough from the kink. Except in the kink zone, the angle $\alpha$ did not significantly influence the RMS error.
This observation is linked to the large number of points available for identification purpose, thanks to the finesse of the finite element mesh. A more significant influence may be observed on a smaller number of DIC measurement nodes.

The observations on the RMS displacement error are confirmed by the direct comparison of SIF levels. Figure 7 compares the identified $K_I$ and $K_{II}$ values to the reference FE solution along the whole propagation path. On these synthetic fields, the results are in very good agreement except close to the kink for small Pacman opening angles. This outcome confirms that to achieve good trustworthiness on SIFs close to the kink, opening the Pacman mouth is necessary. In addition, it does not hinder the robustness far from it. For a given crack tip position, the decrease of the difference with increasing angle $\alpha$, which is particularly visible on $K_{II}$ (Figure 7(d)), confirms that only a small part of the identification zone is responsible for most of the identification error. These findings also suggest that a small external radius $R_{max}$ allows for a better identification by getting closer to the kink. This effect is possible with the use of synthetic fields and very fine meshes. Conversely, for DIC measured fields, there may be an optimum between having enough measurement points to achieve a robust identification and getting close enough to the kink.
<figure><img src="image_8.png" /><figcaption>**Fig. 7:** Mode I and II SIFs (W-FE) along the crack path for different opening angles $\alpha$ in specimen with $R_{max} = 100$ px (a and c) and $R_{max} = 250$ px (b and d).</figcaption></figure>

## 4.2 Crack path modeling

The crack path was positioned from the correlation residual map corresponding to the last acquired image for which the crack extension was the longest. This positioning is faithful for the straight parts of the crack, but the exact location of a singular kink point is more ambiguous. Experimentally, the existence of such sharp kink is debatable. Even with the mesostructure of the studied material creating the kink, it occurs over
a distance of a few filaments (filament width = 10 px). Damage in this zone may also lead to a blurred residual map in addition to the spatial resolution of DIC itself.

The choice between a sharp or smooth description of the crack path is illustrated with I-DIC. Figure 8 shows the SIF profiles measured along these two crack paths. As expected, significant differences arose close to the kink. There was a jump of SIFs for the sharp path while their change was continuous for the smooth path. This difference is particularly visible for the mode II SIF $K_{II}$ (Figure 8(b)). The jump corresponds to the theoretical prediction of Equation (2).

In the case of the sharp path, the choice of the direction of analysis did not impact the positions of the crack tip (Figure 8(c)), except for very few time steps around the kinking event. This observation shows that the exact detection of kinking is delicate. In that region, the usual criterion for pinpointing the crack tip (i.e., canceling out the amplitude of the first supersingular field) did not work very well and oscillations occurred. This effect was in part due to the fact that Williams' fields are only approximations in that region. Thus, in the very close vicinity of the kink, the detection criterion was based on the lowest residual $\varrho$, which indicated that the selected kinematics was the most faithful.
<figure><img src="image_9.png" /><figcaption>**Fig. 8:** Dimensionless SIFs along the crack path. (a) Comparison between smooth and sharp crack paths with direct and reverse I-DIC analyses. (b) Zoom on $K_{II}$ history around the kink. (c) Positions of the crack tip identified along the proposed trajectories.</figcaption></figure>

<footnote n="371">371 All these results prove that the exact occurrence of the kink (i.e., determining the first image after kinking) remains delicate for such experimental study. It is also</footnote>

<footnote>372</footnote>
strongly impacted by other experimental parameters, such as the spatial and temporal
resolutions of the images and the damaged area on the surface of the sample.

## 4.3 Validation of DIC-based SIF measurements

Figures 9 and 10 provide the comparison of W-DIC and I-DIC methods with the reference values for mode I and II SIFs. The present benchmark focuses on a few sets of parameters ($R_{max}$ = 250 px and a sharp kink with reverse analyses). It is representative of the main observed features. First, far from the kink, the results are very consistent with the reference values for both DIC techniques.

The robustness of the identification close to the kink improved for W-DIC when opening the Pacman mouth. Wrong positions for the crack tip were found when $\alpha = 0$ and 5° as this parameter was clearly not suited to that case. As a consequence, the values of SIFs for small Pacman opening angles were also far from the reference close to the kink. However, as the Pacman mouth was opened, both kink position and SIF levels became virtually identical to I-DIC and the reference levels. This result is remarkable since the Q4-DIC mesh was rather coarse, namely, the element size was 16 px, which was much larger than for the computed solution (i.e., 4 px on average and around 2 px close to the crack tip).

Simply opening the Pacman mouth was sufficient to provide very accurate SIFs. Interestingly, I-DIC was a lot less sensitive to the angle $\alpha$. This difference may be due to the fact that the Pacman mesh was a lot finer (i.e., the smallest element size was less than 2 px) in comparison to Q4-DIC. Even though Williams' series are not fully
394 suited when used in the immediate vicinity of the kink, they enabled very consistent

395 SIF estimates to be obtained when compared to the FE reference.

<figure><img src="image_10.png" /><figcaption>**Fig. 9:** Dimensionless mode II SIF measured along the propagation path with I-DIC and W-DIC with four opening angles $\alpha$. Comparison with the FE reference.</figcaption></figure>

396 From a mechanical point of view, it is worth noting that:

397 • before propagation, as awaited from the top/down symmetry of the problem, $K_{II} \approx 0$ and $K_I$ increased with the applied displacement;

398 • after propagation, $K_{II} \neq 0$, meaning that the propagation was not governed by the classical Principle of Local Symmetry [39]. This phenomenon is due to the anisotropy

400
of the fracture energy. It induced kinking in the crack propagation path toward the weakest plane corresponding to one of the deposition directions [26].

* Crack kinking is accompanied by a discontinuity in SIF values, in agreement with LEFM (Section 2.2).

* The SIF jump can be measured accurately provided that the Pacman allows the angular sector of the kink to be excluded from the identification zone (i.e. mouth of the Pacman is open enough not to touch the kink).

<figure><img src="image_11.png" /><figcaption>Fig. 10: Dimensionless mode I SIF measured along the propagation path with I-DIC and W-DIC with four angles $\alpha$. Comparison with the FE reference.</figcaption></figure>
In this study, the abilities of DIC-based methods to estimate stress intensity factors was assessed for a crack path with a sharp kink. The benchmark was based on the analysis of a Compact Tension test on a 3D-printed polycarbonate sample. Despite being loaded in mode I, the manufacturing process induced anisotropy of the fracture properties, thereby causing the specimen to fracture with a rather sharp kink. Two DIC techniques were used to assess SIFs along the propagation path, namely, post-processing DIC measurements and integrated DIC. The results were compared to reference values obtained with an FE model driven by measured Dirichlet boundary conditions.

The good agreement of both mode I and II SIF values close to the kink proved the relevance of these techniques based on the use of Williams' series. Even if this expansion is only valid for straight cracks, the deviation in terms of displacements due to the presence of a kink remained located within a very small area around the main branch. The simple removal of the area in the extraction domain, approximately the angular sector of the kink angle, sufficed to provide very trustworthy results as observed on both computed and experimentally determined data. This property was especially important for the post-processing technique, in which the positioning of the crack tip strongly depended on the extraction zone close to the kink. It was not the case for integrated DIC, which was more robust for any Pacman opening angle.

The discontinuity of SIFs predicted by LEFM was retrieved in FE computations and DIC analyses, provided the considered crack path had a sharp kink. Defining the position and shape of a crack path may be challenging, especially when focusing on
small scales. In that case, the choice of a smooth trajectory, and thus of a continuous change of SIFs may reflect the limits of experimental resolutions regarding the propagation direction changes.

When proper care was taken in the choice of the extraction zone, both DIC-based techniques provided equivalent and very reliable results. To further probe the robustness of the proposed framework, other king angles may be analyzed. The limits of accuracy are likely to mainly stem from other experimental parameters such as spatial resolution of displacement fields around the kink and temporal resolution to properly pinpoint the kink occurrence.

## Acknowledgments

The authors thank Drs. Julie Triclot and Xinyuan Zhai, for discussions about FE modeling. The authors acknowledge the support of Agence Nationale de la Recherche (ANR-23-CE51-0054, project 3FAM). The work was also partially supported by Agence de l'Innovation de Défense – AID (AID-2018 60 0071 00 470 75 01) and via Centre Interdisciplinaire d'Etudes pour la Défense et la Sécurité – CIEDS - (projects 2022 - FracAddi).

## A Finite Element computation of SIFs

### A.1 Method

To compute SIFs, the G_theta procedure implemented in the FE code Cast3M [40] is used to obtain SIFs from the solution of the elasticity problem<fnref n="1" />. The procedure

<footnote n="1">¹Similar procedures are implemented in most industrial FE codes dealing with LEFM</footnote>
451 relies on the concepts of virtual crack extension [37, 38] and interaction integrals [41].
452 Although recognized for providing reliable SIF values when the crack paths are
453 smooth [42], its soundness in the vicinity of a kink remains to be assessed.

454 In practice, two contours Γ₁ and Γ₂ embedding the crack tip are chosen by counting
455 the number of element layers from the tip and a virtual vector field Θ is defined
456 (Figure 11(a)). The SIFs are then obtained by computing surface integrals on a crown
457 between Γ₁ and Γ₂. The integrals involve the displacement and stress elasticity solution
458 and Θ. They are, in principle, independent of Θ and of (Γ₁, Γ₂), provided the part of
459 cracked surface is straight in this zone. However Γ₁ must be far enough away from the
460 crack tip to reduce numerical inaccuracies due to the singularity of the stress field [42].
461 Hence, for kinked crack paths, attention must be paid to the choice of Γ₁ to exclude
462 the kink, while being far away for numerical accuracy. In the vicinity of a kink, this
463 condition becomes delicate due to lack of space, and the results have to be assessed.
464 It is performed in Section A.2 by comparison to exact values given in Equation (2)
465 near the kink.

466 It is worth mentioning that the G_theta method was also used for post-processing
467 DIC measured displacements fields to extract fracture mechanics parameters [38].
468 Noise-robust evaluations were achieved by further tailoring the Θ field. In the present
469 study, more standard extractors were used (i.e., least squares fits based on Williams'
470 series), even though sub-optimal in terms of noise sensitivity [43].
<figure><img src="image_12.png" /><figcaption>**Fig. 11:** (a) Typical $\Gamma_1$, $\Gamma_2$ contours and $\Theta$ field used in the G_theta procedure. The crack surfaces within the blue zone must be parallel to the propagation direction $e_x$. The pre-kink extension was excluded. (b) Dimensionless $K_I(s, \varphi)$ and $K_{II}(s, \varphi)$ obtained by the G_theta procedure (normalized by mode *I* SIF before the kink, $K_I^0$) for different kink angles. The solid lines correspond to Equation (2) that provides the exact solution near the kink. Radius for the G_theta method: $R_{\Gamma_2} = 24$ µm.</figcaption></figure>

## A.2 Validation close to the kink

The relevance of the numerical method near the kink was tested by modeling CT samples with different kink angles and extension lengths (Figure 11). A 2D model with plane strain condition was developed in CAST3M [40]. The model had the same in-plane geometry as the real sample. The focus was not on the simulation of an actual CT test but on the validation of SIFs obtained by the G_theta procedure after the kink. It was simply chosen to model the boundary conditions by applying a unit
displacement at the centers of the bottom and top pins, respectively. The key point
is to keep the same load while changing the crack length, so that one can compare
the change of SIFs after the kink to Equation (2). The model was meshed with linear
triangular elements except for a 40 µm circular zone around the tip, where a structured
mesh with very fine quadrangular elements was used (the size along the radius was
4 µm). The integration surface Γ₁, Γ₂ was defined as a crown following one layer of
this circular mesh. Convergence of SIF values was reached locating the crown in the
6-th layer. Figure 11(b) shows Kᵢ and Kᵢᵢ values, divided by the mode I SIF before
the kink, for kink angles ranging from 5° to 75° along kinked extensions.

Theoretical values were computed with Equation (2) using the SIFs and *T*-stress
values before the kink obtained from a simulation with no kink. The *T*-stress was
estimated from the stress field, with the difference $\sigma_{xx}(r, \theta=0) - \sigma_{yy}(r, \theta=0)$ averaged
over a length where this value is approximately constant (before the dominance of
higher order terms). This classical method, presented for instance in [44], provided
sufficiently accurate results for the present case. The very good agreement between the
numerical predictions and the theoretical levels of SIFs across the kink validate the
use of the *G_theta* procedure close to the kink, for crack extensions as short as 40 µm
 corresp
 (corresponding to 2 pixels in the experiment). The validated method was applied to
compute the change of SIFs in the experiment by using faithful crack positions and
boundary conditions derived from DIC analyses (Section 3.3).
[1] Deegan, R., Chheda, S., Patel, L., Marder, M., Swinney, H., Kim, J., Lozanne, A.: Wavy and rough cracks in silicon. Physical review E **67**(6, Part 2) (2003)

[2] Takei, A., Roman, B., Bico, J., Hamm, E., Melo, F.: Forbidden directions for the fracture of thin anisotropic sheets: An analogy with the wulff plot. Physical Review Letters **110**, 144301 (2013)

[3] Ngo, T.D., Kashani, A., Imbalzano, G., Nguyen, K.T.Q., Hui, D.: Additive manufacturing (3d printing): A review of materials, methods, applications and challenges. Composites Part B: Engineering **143**, 172–196 (2018)

[4] Gardan, J., Makke, A., Recho, N.: Improving the fracture toughness of 3D printed thermoplastic polymers by fused deposition modeling. International Journal of Fracture **210**, 1–15 (2018)

[5] Bui, H.D.: Fracture mechanics: inverse problems and solutions. In: Solid Mechanics and Its Applications vol. 139. Springer, Dordrecht (2006)

[6] Corre, T., Lazarus, V.: Kinked crack paths in polycarbonate samples printed by fused deposition modelling using criss-cross patterns. International Journal of Fracture **230**, 19–31 (2021)

[7] Leblond, J.-B.: Crack Paths in Plane Situations - I. General Form of the Expansion of the Stress Intensity Factors. International Journal of Solids and Structures **25**, 1311–1325 (1989)
[8] Amestoy, M., Leblond, J.-B.: Crack Paths in Plane Situations - II. Detailed Form of the Expansion of the Stress Intensity Factors. *International Journal of Solids and Structures* **29**, 465–501 (1992)

[9] Roux, S., Réthoré, J., Hild, F.: Digital Image Correlation and Fracture: An Advanced Technique for Estimating Stress Intensity Factors of 2D and 3D Cracks. *Journal of Physics D: Applied Physics: Appl. Phys.* **42**, 214004 (2009)

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

[11] Peters, W., Ranson, W., Kalthoff, J., Winkler, S.: A study of dynamic near-crack-tip fracture parameters by digital image analysis. *Journal de Physique Colloques* **46**(C5), 631–638 (1985)

[12] Dawicke, D.S., Sutton, M.S.: Ctoa and crack-tunneling measurements in thin sheet 2024-t3 aluminum alloy. *Experimental Mechanics* **34**, 357–368 (1994)

[13] Sutton, M.A., Zhao, W., McNeill, S.R., Helm, J.D., Piascik, R.S., Riddel, W.T.: Local crack closure measurements: Development of a measurement system using computer vision and a far-field microscope. In: McClung, R.C., Newman Jr., J.C. (eds.) *Advances in Fatigue Crack Closure Measurement and Analysis: Second Volume, STP 1343*, pp. 145–156. ASTM, West Conshohocken, PA (USA) (1999)

[14] Sutton, M.A., McNeill, S.R., Helm, J.D., Chao, Y.J.: Advances in two-dimensional and three-dimensional computer vision. In: Rastogi, P.K. (ed.) *Photomechanics*,
vol. Topics in Appl. Phys., 77, pp. 323–372. Springer, Berlin (Germany) (2000)

[15] Forquin, P., Rota, L., Charles, Y., Hild, F.: A method to determine the toughness scatter of brittle materials. International Journal of Fracture **125**(1), 171–187 (2004)

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

[17] Dunn, M.L., Suwito, W., Cunningham, S.: Fracture initiation at sharp notches: Correlation using critical stress intensities. International Journal of Solids and Structures **34**(29), 3873–3883 (1997)

[18] Limodin, N., Réthoré, J., Buffière, J.Y., Gravouil, A., Hild, F., Roux, S.: Crack closure and stress intensity factor measurements in nodular graphite cast iron using 3D correlation of laboratory X ray microtomography images. Acta Materialia **57**(14), 4090–4101 (2009)

[19] Rannou, J., Limodin, N., Réthoré, J., Gravouil, A., Ludwig, W., Baïetto, M.C., Buffière, J.Y., Combescure, A., Hild, F., Roux, S.: Three dimensional experimental and numerical multiscale analysis of a fatigue crack. Computer Methods in Applied Mechanics and Engineering **199**, 1307–1325 (2010)

[20] Humbert, L., Valle, V., Cottron, M.: Experimental determination and empirical representation of out-of-plane displacements in a cracked elastic plate loaded in
558

mode I. International Journal of Solids and Structures **37**(39), 5493–5504 (2000)

559

[21] Bažant, Z.P., Estenssoro, L.F.: Surface singularity and crack propagation. International Journal of Solids and Structures **15**(5), 405–426 (1979)

560

561

[22] Zakavi, B., Kotousov, A., Branco, R.: Overview of three-dimensional linear-elastic fracture mechanics. International Journal of Fracture **234**(1), 5–20 (2022)

562

563

[23] Williams, M.L.: On the stress distribution at the base of a stationary crack. ASME Journal of Applied Mechanics **24**, 109–114 (1957)

564

565

[24] Tada, H., Paris, P.C., Irwin, G.R.: The Stress Analysis of Cracks Handbook. Professional Engineering Publishing, New York, NY (USA) (2000)

566

567

[25] Cottrell, B., Rice, J.R.: Slightly curved or kinked cracks. International Journal of Fracture **16**(2), 155–169 (1980)

568

569

[26] Zhai, X., Corre, T., Lazarus, V.: A FDM-based experimental benchmark for evaluating quasistatic crack propagation in anisotropic linear elastic materials. Engineering Fracture Mechanics **324**, 111175 (2025)

570

571

572

[27] Broggiato, G.B.: Adaptive image correlation technique for full-field strain measurement. In: Pappalettere, C. (ed.) 12th Int. Conf. Experimental Mechanics, pp. 420–421. McGraw Hill, Lilan (Italy) (2004)

573

574

575

[28] Sun, Y., Pang, J., Wong, C., Su, F.: Finite-element formulation for a digital image correlation method. Applied Optics **44**(34), 7357–7363 (2005)

576
[29] Besnard, G., Hild, F., Roux, S.: "Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Chatelier bands. Experimental Mechanics **46**, 789-803 (2006)

[30] Réthoré, J.: Automatic crack tip detection and stress intensity factors estimation of curved cracks from digital images. International Journal for Numerical Methods in Engineering **103**(7), 516-534 (2015)

[31] Williams, M.L.: Stress singularities resulting from various boundary conditions in angular corners of plates in extension. Journal of Applied Mechanics **19**, 526 (1952)

[32] Gamot, J., Lasserre, T., Richard, L., Neggers, J., Swiergiel, N., Hild, F.: Calibrating thermoelastic stress analysis with integrated digital image correlation: Application to fatigue cracks. Journal of Strain Analysis for Engineering Design **54**(5-6), 320-330 (2019)

[33] Réthoré, J.: UFreckles. https://doi.org/10.5281/zenodo.1433776

[34] Vargas, R., Neggers, J., Canto, R.B., Rodrigues, J.A., Hild, F.: Analysis of wedge splitting test on refractory castable via integrated DIC. Journal of the European Ceramic Society **36**, 4309-4317 (2016)

[35] Leclerc, H., Neggers, J., Mathieu, F., Hild, F., Roux, S.: Correli 3.0. IDDN.FR.001.520008.000.S.P.2015.000.31500, Agence pour la Protection des Programmes, Paris (France) (2015)
[36] Triclot, J., Corre, T., Gravouil, A., Lazarus, V.: Key role of boundary conditions for the 2d modeling of crack propagation in linear elastic compact tension tests. *Engineering Fracture Mechanics* **277**, 109012 (2023)

[37] Destuynder, P., Djaoua, M.: Sur une interprétation mathématique de l’intégrale de Rice en théorie de la rupture fragile. *Mathematical Methods in the Applied Sciences* **3**(1), 70–87 (1981)

[38] Réthoré, J., Gravouil, A., Morestin, F., Combescure, A.: Estimation of mixed-mode stress intensity factors using digital image correlation and an interaction integral. *International Journal of Fracture* **132**, 65–79 (2005)

[39] Goldstein, R.V., Salganik, R.L.: Brittle fracture of solids with arbitrary cracks. *International Journal of Fracture* **10**, 507–523 (1974)

[40] CEA (French Atomic Commission): Cast3M: Finite Element program. www-cast3m.cea.fr (2017)

[41] Gosz, M., Dolbow, J., Moran, B.: Domain integral formulation for stress intensity factor computation along curved three-dimensional interface cracks. *International Journal of Solids and Structures* **35**(15), 1763–1783 (1998)

[42] Suo, X.-Z., Combescure, A.: On the application of $g(\theta)$ method and its comparison with de lorenzi’s approach. *Nuclear Engineering and Design* **135**(2), 207–224 (1992)
616

617

618

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

619

620

621

[44] Ayatollahi, M., Pavier, M., Smith, D.: Determination of t-stress from finite element analysis for mode i and mixed mode i/ii loading. *International journal of fracture* **91**, 283–298 (1998)