# A FDM-based experimental benchmark for evaluating quasistatic crack propagation in anisotropic linear elastic materials

Xinyuan Zhai, Thomas Corre, Véronique Lazarus

► To cite this version:

Xinyuan Zhai, Thomas Corre, Véronique Lazarus. A FDM-based experimental benchmark for evaluating quasistatic crack propagation in anisotropic linear elastic materials. Engineering Fracture Mechanics, In press, 10.1016/j.engfracmech.2025.111175. hal-05042969

HAL Id: hal-05042969

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

Submitted on 22 Apr 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
A FDM-based experimental benchmark for evaluating
quasistatic crack propagation in anisotropic linear elastic
materials

Xinyuan Zhaiª, Thomas Correᵇ, Véronique Lazarusª

ªIMSIA, CNRS, EDF, ENSTA, Institut Polytechnique de Paris, 828 boulevard des Maréchaux, Palaiseau, 91762 cedex, France

ᵇÉcole Centrale de Nantes, 1 Rue de la Noë, Nantes, 44300, France

**Abstract**

Anisotropic materials gain importance in industry pushed forward by new appli-
cations, from wood construction, single crystals for photovoltaic applications or
aerospace turbines, to Additive Manufactured materials. But due to the lack of
certified damage tolerance models to guarantee their resistance to failure, their use
in aerospace or nuclear sensitive components whose failure must be avoided at all
costs is limited. In this context, we present herein a new experimental benchmark,
instrumented by Digital Image Correlation, aiming to test quasistatic crack propaga-
tion models for anisotropic linear elastic materials, for instance phase field variational
approaches. The benchmark is based on Compact Tensile and Compact Tensile Shear
samples printed by Fused Deposit Modeling of polycarbonate threads. The flexibility
of the printing process is used to produce specimens that exhibit directional depen-
dence of the fracture toughness. Fracture experiments with these samples have been
designed to enter in the framework of two-dimensional (2D) Linear Elasticity Frac-
ture Mechanics (LEFM) with strong fracture anisotropy, while avoiding additional
complexities by preserving isotropic elasticity and quasistatic propagation. Close
comparisons with the experiments demonstrate that the generalized maximum en-
ergy release rate criterion accurately predicts crack propagation. While serving as
a benchmark, these experiments also lay the foundation for a better fundamental
understanding of crack propagation in anisotropic materials, beyond FDM printed
parts.

Keywords:

Linear Elastic Fracture Mechanics (LEFM), Strong fracture anisotropy, Additive
Manufacturing, Experimental benchmark, Fused Deposit Modeling (FDM)
Nomenclature

<table>
  <tr>
    <td>a</td>
    <td>Crack length</td>
  </tr>
  <tr>
    <td>a<sub>I</sub><sup>n</sup></td>
    <td>Asymptotic coefficients for Mode I</td>
  </tr>
  <tr>
    <td>a<sub>II</sub><sup>n</sup></td>
    <td>Asymptotic coefficients for Mode II</td>
  </tr>
  <tr>
    <td>&beta;</td>
    <td>Loading angle</td>
  </tr>
  <tr>
    <td>e</td>
    <td>Sample thickness</td>
  </tr>
  <tr>
    <td>E</td>
    <td>Young's modulus</td>
  </tr>
  <tr>
    <td>E'</td>
    <td>Effective modulus (plane stress/strain)</td>
  </tr>
  <tr>
    <td>F</td>
    <td>Applied force</td>
  </tr>
  <tr>
    <td>G</td>
    <td>Energy Release Rate</td>
  </tr>
  <tr>
    <td>G*</td>
    <td>Energy Release Rate at kink initiation</td>
  </tr>
  <tr>
    <td>G<sub>c</sub></td>
    <td>Critical surface energy</td>
  </tr>
  <tr>
    <td>K<sub>I</sub></td>
    <td>Mode I Stress Intensity Factor</td>
  </tr>
  <tr>
    <td>K<sub>II</sub></td>
    <td>Mode II Stress Intensity Factor</td>
  </tr>
  <tr>
    <td>K<sub>IC</sub></td>
    <td>Fracture toughness in Mode I</td>
  </tr>
  <tr>
    <td>L, W</td>
    <td>Sample length and width</td>
  </tr>
  <tr>
    <td>r, &theta;</td>
    <td>Polar coordinates with the origin located at the crack tip</td>
  </tr>
  <tr>
    <td>s, &phi;</td>
    <td>Kink extension length</td>
  </tr>
  <tr>
    <td>&phi;</td>
    <td>Kink angle</td>
  </tr>
  <tr>
    <td>S&#770;</td>
    <td>Sample stiffness</td>
  </tr>
  <tr>
    <td>T</td>
    <td>T-stress</td>
  </tr>
  <tr>
    <td>u</td>
    <td>Displacement</td>
  </tr>
  <tr>
    <td>&Gamma;</td>
    <td>Fracture energy</td>
  </tr>
  <tr>
    <td>&Gamma;<sub>0</sub></td>
    <td>Fracture energy of the weakest plane</td>
  </tr>
  <tr>
    <td>&Gamma;<sub>max</sub></td>
    <td>Maximum fracture energy</td>
  </tr>
  <tr>
    <td>&Delta;&Gamma;</td>
    <td>Fracture anisotropy (&Gamma;<sub>max</sub> &minus; &Gamma;<sub>0</sub>)</td>
  </tr>
  <tr>
    <td>&mu;</td>
    <td>Shear modulus</td>
  </tr>
  <tr>
    <td>&nu;</td>
    <td>Poisson's ratio</td>
  </tr>
  <tr>
    <td>&theta;<sub>0</sub></td>
    <td>Material axis direction of surface energy</td>
  </tr>
  <tr>
    <td>&epsilon;<sub>p</sub></td>
    <td>Potential energy</td>
  </tr>
  <tr>
    <td>&sigma;<sub>c</sub></td>
    <td>Critical stress</td>
  </tr>
  <tr>
    <td>AM</td>
    <td>Additive Manufacturing</td>
  </tr>
  <tr>
    <td>COD</td>
    <td>Crack Opening Displacement</td>
  </tr>
  <tr>
    <td>CT, CTS</td>
    <td>Compact Tension sample, Compact Tension Shear sample</td>
  </tr>
  <tr>
    <td>DIC</td>
    <td>Digital Image Correlation</td>
  </tr>
  <tr>
    <td>FDM</td>
    <td>Fused Deposition Modeling</td>
  </tr>
  <tr>
    <td>FEM</td>
    <td>Finite Element Method</td>
  </tr>
  <tr>
    <td>GMERR</td>
    <td>Generalized Maximum Energy Release Rate</td>
  </tr>
  <tr>
    <td>LEFM</td>
    <td>Linear Elastic Fracture Mechanics</td>
  </tr>
  <tr>
    <td>PC</td>
    <td>Polycarbonate</td>
  </tr>
  <tr>
    <td>P-&sigma;</td>
    <td>Plane stress</td>
  </tr>
  <tr>
    <td>P-&epsilon;</td>
    <td>Plane strain</td>
  </tr>
  <tr>
    <td>SIF</td>
    <td>Stress Intensity Factor</td>
  </tr>
</table>
## 1. Introduction

Determining how and when a pre-existing crack will propagate is of crucial importance for sensitive components whose failure would endanger life, particularly in the aerospace and nuclear industries. Damage tolerance approaches based on Linear Elastic Fracture Mechanics [1] has shown its efficiency in this respect. The idea, based on the fact that almost no material is perfect, is to assess that any possible defect will not cause catastrophic failure taking the worst-case scenario of a crack propagating in a brittle material.

Well-established for 2D linear elastic problems, involving isotropic and homogeneous materials, simple loads and geometries, fatigue or quasistatic slow crack propagation [2], the determination of crack propagation still presents scientific challenges in taking into account heterogeneities [3], 3D aspects [4, 5], dynamic propagation [6, 7, 8] or material anisotropy [9, 10, 11].

Herein, we will focus on the effect of material anisotropy on quasistatic propagation, proposing benchmark experiments designed to avoid the other complexities, as heterogeneities, 3D aspects, dynamics, for instance.

Anisotropic materials gain importance in industry pushed forward by new applications, from wood [12], single crystals for photovoltaic applications [13] or aerospace turbines [14], to Additive Manufactured (AM) materials [15, 16, 17, 18, 19]. However, the lack of validated tolerance damage models for strongly anisotropic materials limits their broader adoption. Our benchmark aims to provide a robust experimental framework for validating fracture models while enabling a deeper understanding of anisotropic fracture problem.

For isotropic linear elasticity, the Maximum Energy Release Rate (MERR) together with Griffith's energy balance has shown its efficiency to determine the propagation of a crack [20]. Both can be rewritten in a more general framework using variational approaches [21]. For numerical purposes, regularisation by Phase-Field (PF) models has giving rise to a huge number of formulations [22, 23, 24], whose convergence to the variational approach was assessed in the sense of $\Gamma$-convergence [25]. However, the final selection toward the so-called AT1 model could only be done by comparison with experiments [26]. For anisotropic materials, extensions of the MERR have been used to predict crack paths in various contexts [27, 28, 29, 30, 10]. The formal reformulation in term of variational approaches and phase-field models is quite straight-forward, with some notable differences: the presence of a kink in the crack path poses some theoretical issues [31], no-proof of $\Gamma$-convergence is available yet (to the best of our knowledge) and the selection of the correct phase field formulations among the proposed ones [32, 33, 34, 9, 35, 36] in comparison to experiments is only in its early stages [30, 10, 13, 37].

More precisely, Takei et al. [30] is considering tearing of weakly anisotropic sheets, highlighting the existence of forbidden directions, but those experiments are not pertinent for planar elasticity problem since the loading is out of plane and the
70

71

72

73

74

75

76

77

78

79

80

81

82

elastic energy could be neglected [38]. Mesgarnejad et al. [10] considered the fracture
of fiber composites highlighting the necessity to take T-stress and a finite length
into account. But no measure of the critical load was reported and failure in part
of the experiments occurred dynamically, that is beyond the quasistatic framework.
Nagaraja et al. [13] study crack propagation in single Silica crystal, comparing closely
to several PF models, but those experiments are complex to set up and comparison to
quasistatic models may be questioned since the setup gave rise to abrupt propagation.
Lampron et al. [37] compared experimental results of FDM-printed PLA samples with
a multi-damage phase-field model. In their work, crack propagation exhibited diffuse
damage, with significant plastic deformation, making it unsuitable for validating
models within the LEFM framework.

83

84

In this paper, we present a new experimental benchmark for testing quasistatic
crack propagation models in anisotropic linear elastic material.

It uses Compact Tensile (CT) and Compact Tensile Shear (CTS) samples of
polycarbonate printed by

85

86

87

88

89

90

91

Fused Deposit Modeling (FDM) as proposed by Corre and Lazarus [39]. These experiments are shown to enter in the framework of 2D Linear Elasticity Fracture Mechanics (LEFM) with strong fracture anisotropy, while avoiding additional complexities by preserving isotropic elasticity and quasistatic crack propagation. The aim here is not to focus on practical application in the widely studied context of FDM, but to propose low-cost experiments designed for a better understanding and modeling of crack propagation in anisotropic materials.

92

93

94

95

96

97

98

99

100

An extensive study of the crack propagation is presented for different mixed mode
I+II loading and anisotropy orientation, combining experiments (section 2) and a
numerical analysis based on LEFM (section 3). The results (section 4) are well de-
scribed by the Generalized Maximum Energy Release Rate criterion (GMERR) based
on Griffith energy balance. This demonstrates the relevance of these experiments for
validating models that aim to simulate crack propagation within the LEFM frame-
work, in particular variational phase field approaches that aim to predict the crack
path according to the GMERR criterion without explicitly meshing a sharp crack
during propagation.

101

## 2. Fracture experiments

102

### 2.1. *Sample printing*

103

104

105

106

107

108

109

The samples were produced by Fused Deposition Modeling with Polycarbonate (PC) filament, specifically, *PC – Plus*TM from Poly-maker, with a filament diameter of 1.75 mm. The digital model of the sample was imported into IdeaMaker to be sliced and generate a G-code file that lists the commands of the 3D printer. The machine used is the Raise3D Pro2. To prevent warping of the sample during the process, the build platform was preheated to 100 °C for one hour and a raft was added, ensuring a proper calibration of the first layer.
The printing process begins by heating the nozzle to 260 °C and extruding a continuous filament of PC through a 0.4 mm diameter nozzle. The material is heated to a semi-liquid state and then deposited onto the platform or previous layers by moving the nozzle. Over time, the deposited material is gradually cooling down and solidifying to bond with the previously deposited thread. The printing speed for infilling the sample was set to 50 mm/s, within the range recommended by Poly-maker for PC (30 - 50 mm/s). The thickness of each layer was approximately 0.15 mm and each layer was filled with an infill rate of 100 %. Yet, due to polymer retraction during printing, the samples exhibit a small porosity, between 1.3 % and 2.5 % as estimated from the measured mass and the theoretical volume.

To get a benchmark for studying crack propagation in presence of strong anisotropic fracture energy, within 2D isotropic elasticity, the printing strategy described in Corre and Lazarus [39] was used. In this strategy, each layer was filled with a selected infill angle or *raster angle* (angle of deposited material with respect to the notch), and successive layers were arranged orthogonally to each other. Fig. 1 shows the principle of FDM and the printing strategy for the CT sample.

The fracture behavior of this 3D printed material is investigated with Compact Tensile samples and Compact Tensile Shear samples (Fig. 2). The displacement is imposed through pins inserted in the holes. The symmetry of the CT sample ensures mode I conditions, while the inclination of the holes in a CTS sample induces an additional amount of mode II. The *loading angle* denoted $\beta$, is defined as the inclination of the sample versus the axes of the tensile machine, that corresponds also to the inclination of the holes axis.

<figure><img src="image_2.png" /><figcaption>Figure 1: (Color online) (A): FDM printer. (B): Preview of the printing process for the outer shell (in red) and infill threads (in yellow) of the initial layer (the shell is not printed at the notch extremity). (C): Orthogonal arrangement of the printing threads.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: (Color online) a) Dimensions of the CT sample, b) CT sample with a speckle pattern for Digital Image Correlation (DIC) analysis, c) CTS sample with adjoined pin holes used to create a pre-crack by mode I fatigue cyclic loading and d) CTS sample after fatigue pre-cracking, mode I pin holes are removed. Green and blue zones aligned with the pin holes in (B) and (C) are the domains between which the Crack Opening Displacement (COD) is measured by DIC.</figcaption></figure>

## 2.2. Initial notch and fatigue pre-cracking

To facilitate the crack initiation, the outer shell near the crack tip was not printed. To achieve this, we modified systematically the G-code file using a Python script to order the printing of the red line on Fig. 1B,C. Nevertheless, the microstructure just ahead of the notch is different from the criss-cross pattern, because of the complex movements of the nozzle is this area. Fig. 3(A) illustrates the details of this microstructure, with the outer shell highlighted in red and the raster arrangement of the bulk underlined in yellow.

To ensure that the crack propagate directly into the intended anisotropic bulk material, a fatigue pre-crack was added to the notched specimen by applying mode I
loading using an Instron ElectroPuls 3000 Dynamic tensile machine. The specimens were subjected to a sinusoidal force at 3 Hz, with predefined maximum and minimum force amplitudes. After 20,000 cycles, a pre-crack of approximately 1 to 1.5 mm length was initiated. Besides ensuring that the initial crack is within the bulk anisotropic material (Fig. 3(A)), we observe that, as awaited, the presence of the pre-crack significantly helps to stabilize crack propagation, particularly under Mode I+II loading conditions: a quasi-static crack propagation was observed for most of the experiments, although without pre-cracking, a brutal failure was observed in most cases. Note that structures resembling ‘Mickey ears’, shown in fig. 2(C), were used to initiate a mode I fatigue pre-crack in the CTS sample, then removed before the fracture test. Fig. 3(B) shows the crack tip of a CT sample with a 25° raster angle once propagation occurred, highlighting the kinked extension occurring at the tip of the fatigue pre-crack.

<figure><img src="image_4.png" /><figcaption>Figure 3: (A) Details near the notch at the outer shell of the printed sample. (B) Zoom in on the crack tip as crack kink occurs (with speckles pattern)</figcaption></figure>

## 2.3. *Fracture experiments*

Fracture tests were carried out using an Instron 5596 tensile machine, with a displacement imposed at the upper pin hole at a rate of 0.5 mm.min⁻¹. The experiment is monitored with Nikon D300s camera mounted with an AF-S Micro NIKKOR 60 mm lens, taking 4288x2848 pixel images every 4 seconds (physical size of pixel is approximately 20 µm). Data collection, load and displacement, and image capture were synchronised by an Arduino system.

As benchmark, we choose to study extensively CT tests and CTS tests with a 15° loading angle, varying for each case, the raster angle among 0°, 5°, 15°, 25°, 35°, and 45°.

Each test was repeated three times, resulting in a total of 36 experiments. Each one was instrumented by Digital Image Correlation (DIC). A few other loading angles were also tested, but without systematic DIC post-processing. In the sequel, we
169 denote CT\_R15 for CT samples with a 15° raster angle and CTS\_L15\_R25 for
170 CTS with a 15° loading angle and 25° raster angle, and so on.

## 2.4. Experimental data

171 In this section, we present the raw experimental results. Section 2.4.1 presents the force-displacement curves and section 2.4.2 presents representative pictures of the crack paths observed in CT and CTS samples with various raster angles.

### 2.4.1. Force-displacement curves

175 Fig. 4 illustrates the force-displacement curves for the CT (left) and CTS (right) samples, for a selection of representative raster angles. They all exhibit a peak-shaped curve, with the crack starting to advance at the maximum. Before the peak, a nearly linear behavior and similar slopes are observed regardless of the printing orientation, in accordance with a 2D linear elastic isotropic constitutive law. After the peak, the force gradually decreases as the crack propagates quasi-statically.

<figure><img src="image_5.png" /><figcaption>Figure 4: (Color online) Force-displacement curves for CT (left) and CTS (right) sample with different raster angles.</figcaption></figure>

### 2.4.2. Crack paths

182 Fig. 5 shows images of the crack path for a selection of raster angles for CT and 183 CTS samples. Regardless of the loading conditions or raster angle, the crack propagated exclusively along the printing direction. This was consistently observed across 184 all tested samples. This suggests that crack propagation in directions other than 185 the printing direction is impossible in our material, in line with strongly anisotropic 186 fracture properties [30]. 187
<figure><img src="image_6.png" /><figcaption>Figure 5: Crack paths in CT and CTS samples with various selected raster angles</figcaption></figure>

## 2.5. Digital Image Correlation

To perform DIC analysis, the specimens were coated with a speckle pattern, as illustrated in Fig. 2 (B) and (D). Tracking the speckle pattern on the specimen allowed for the measurement of displacement and strain on the specimen surface. The displacement field was then used to determine the Crack Opening Displacement (COD) (§2.5.1), as well as the crack tip position and the Stress Intensity Factors (SIFs) through projection onto Williams' series (§2.5.2). The open-source software UFreckles was used for these tasks [40, 41].

### 2.5.1. Crack Opening Displacement

The COD is defined as the difference in y-axis displacement between two points aligned with the center of the CT and CTS sample pin hole (blue and green dots in Fig. 2 (B) and (C)).

To enhance measurement accuracy, displacement values were averaged from multiple points along the x- and y-axis in two distinct regions (about 32px × 32px) near the crack lips, rather than relying on a single point.
### 2.5.2. Crack tip position and Stress Intensity Factors

Fracture mechanics parameters were extracted by DIC, as detailed in previous works [42, 40, 39]. The post-processing involves projecting the extracted displacement field $u$ near the crack tip onto the Williams' series [43]:

$$u(r, \theta) = \sum_{n=-\infty}^{\infty} a_I^n \psi_I^n(r, \theta) + a_{II}^n \psi_{II}^n(r, \theta) \quad (1)$$
where $a_I^n$ and $a_{II}^n$ are the asymptotic coefficients and $\psi_I^n$ the terms of the solution fields series.

The extraction domain is a disk-shaped region with a radius denoted as $r_{ext}$, centered on the crack tip as shown in Fig. 6. A predefined crack path was proposed based on the discontinuity in the measured displacement field to assist in detecting the crack tip position. Initially, the coefficient $a_I^n$ was calculated using an estimate of the crack tip position. To more accurately determine the crack tip location, the assumed crack tip was then shifted by a distance $d$ along the predefined crack path to eliminate the most dominant supersingular term $a_I^{-1}$ [44]:

$$d = \frac{2a_I^{-1}}{a_I^1} \quad (2)$$
The crack tip detection process is considered complete after several iterations, once the shift distance $d$ is reduced to less than one pixel. Subsequently, the displacement field of the extracted domain at the detected crack tip is projected onto the Williams' series to derive the asymptotic coefficients. These coefficients are then used to compute fracture mechanics parameters, such as SIFs:

$$K_I = 2\mu\sqrt{2\pi a_I^1}, \quad K_{II} = 2\mu\sqrt{2\pi a_{II}^1} \quad (3)$$
where $\mu$ is the shear modulus.

## 2.6. Material characterisation

### 2.6.1. Elastic constants

The elastic behavior of the material was determined using classical tensile tests on dog-bone samples until failure. It was found that the material is nearly isotropic in the plane of printing, with $E = 1890 \pm 111$ MPa, $\nu = 0.34 \pm 0.07$ and the failure stress $\sigma_c = 42 \pm 2$ MPa. Section 4.2 will show that the elastic constants are coherent with the stiffness of the fracture samples (Fig. 4). Further details are provided in Appendix A.

### 2.6.2. Fracture energy $\Gamma(\theta)$ shape

Due to the criss-cross printing strategy, which results in an orthotropic arrangement of the threads, we choose the form of the anisotropic fracture energy presented in Zhai et al. [11]:
<figure><img src="image_7.png" /><figcaption>Figure 6: (Color online) Left: measured displacement field obtained by DIC. Right: extracted domain near the crack tip, the red line represents the proposed crack path along which the crack tip position is determined</figcaption></figure>

$$ \Gamma(\theta) = \Gamma_0 \left( 1 + \frac{\Delta\Gamma}{\Gamma_0} \left| \sin 2(\theta - \theta_0) \right| \right) \qquad (4) $$

Here, $\Gamma_0$ represents the critical fracture energy of the weakest planes, and $\Delta\Gamma$ is defined as $\Gamma_{max} - \Gamma_0$, where $\Gamma_{max}$ is the maximum fracture energy in the four threaded directions. The directions corresponding to $\Gamma_0$ and $\Gamma_{max}$ are related by a $\pi/4$ rotation. Consequently, $\frac{\Delta\Gamma}{\Gamma_0}$ represents the degree of anisotropy of the surface energy. $\theta_0$ denotes the material axis of orientation-dependent surface energy. In practice, $\theta_0$ corresponds to the raster angle in our samples, which defines the direction of the closest to the crack, weakest plane. Fig. 7 presents a polar plot of the reciprocal four-fold anisotropic surface energy with $\theta_0 = 45^\circ$. The surface energy becomes non-convex, hence entering in the category of strongly anisotropic materials [9], when $\frac{\Delta\Gamma}{\Gamma_0}$ exceeds 0.4.

### 2.6.3. Determination of $\Gamma_0$

The critical fracture energy of the weakest planes, $\Gamma_0$, can be directly estimated using force-displacement curves. Following Triclot et al. [45],

the potential energy per unit thickness at each time step can be calculated by:

$$ \mathcal{E}_p^i = \frac{1}{e} \left( \frac{1}{2} F_i u_i - \int_0^{u_i} F du \right) \qquad (5) $$

where $u_i$ and force $F_i$ are the measured displacement and force, $e$ is the sample thickness.

Plotting the potential energy against the crack length $a$ determined through DIC, revealed that $\mathcal{E}_p$ decreases nearly linearly.
<figure><img src="image_8.png" /><figcaption>Figure 7: Polar plot of 1/Γ(θ) from Eq. (4) with θ₀ = 45° and different values of ΔΓ/Γ₀</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Plot of potential energy versus crack length for different samples</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Fracture energy of the weakest planes $\Gamma_0$, measured using 18 CT samples and 18 CTS samples with different raster angles and 15° loading angle. The red dashed line represents the average value, and the grey shaded area indicates the standard deviation.</figcaption></figure>

This linearity is consistent with the Griffith criterion $G = \Gamma_0$ so that the slope of the plot provides an estimate of $\Gamma_0$, as shown by the dashed lines in Fig. 8.

Fig. 9 reports the critical fracture energy of several tested CT and CTS samples plotted against the raster angle. The calculated value of $\Gamma_0$ is approximately 4.6 kJ/m², with a standard deviation of 0.66 kJ/m².

This measure agrees with the values reported in the literature for PC, which range from 1 kJ/m² to 5 kJ/m² [46, 47]. Notice that the value of $\Gamma_0$ is independent of the raster angle and loading device CT vs CTS, showing that the fracture energy of the weakest planes is independent of the fracture mode. The measured mean material parameters are summarized in Table 1.

<table><caption>Table 1: Average material properties of the printed samples</caption><thead><tr><th colspan="4">Material properties</th></tr><tr><th><i>E</i> (MPa)</th><th><i>ν</i></th><th>σ<sub>c</sub> (MPa)</th><th>Γ<sub>0</sub> (kJ/m<sup>2</sup>)</th></tr></thead><tbody><tr><td>1890 (std 111)</td><td>0.34 (std 0.07)</td><td>42 (std 2)</td><td>4.6 (std 0.66)</td></tr></tbody></table>
### 3. Crack propagation determination by Linear Elastic Fracture Mechanics

In this section, we use the GMERR criteria and a Finite Element Method (FEM) model within the framework of LEFM to reproduce the experimental results. Our objective is to determine the kink direction and the related critical force at equilibrium.

Our approach is as follows. First, we perform a dimensional analysis and solve by FEM the dimensionless equilibrium elasticity problem to determine the SIFs and Energy Release Rate (ERR) for (i) the initial geometry and (ii) the case with a given crack extension of dimensionless length $\hat{s}$ in the direction $\varphi$ (section 3.1.1). To supplement FEM for a vanishing small crack extension length, Amestoy and Leblond [48] formula will be used (section 3.1.2).

Second, we use these values in the GMERR criterion, that puts in balance the fracture costs with the elastic energy available to advance the crack, in order to get (i) the crack propagation direction $\varphi^*$, (ii) the force $F_c(\varphi, s)$ required to get quasi-static crack propagation ruled by $G = \Gamma_0$ from initiation (section 4.3.1) to the quasistatic propagation phases (section 4.3.2).

#### 3.1. Determination of the sample stiffness, the SIFs and the Energy Release Rate

##### 3.1.1. Two dimensional elasticity problem

<figure><img src="image_11.png" /><figcaption>Figure 10: (Color online) Elasticity problem for CT (left) and CTS (right) sample, where $s$ represents the kink length and $\varphi$ denotes the kink direction. The difference in the $y$-axis displacement between blue (upper side) and green point (lower side) is defined as the Crack Opening Displacement. These two points are aligned with the center of the pin along the $y$-axis.</figcaption></figure>
The experiments are modeled by considering the 2D elastic problem depicted in Fig. 10. A force is applied at the top and the displacement is fixed at the bottom. The experimental sample shape (Fig. 2) is preserved, taking as reference unit crack length $a_0 = 26$ mm corresponding to the mean crack length before kinking in the experiments. The kink direction $\varphi$ and the extension length $s$ are obtained from the crack tip positions measured using DIC.

For a given crack geometry parameterized by $s$, $\varphi$, the problem is solved using the Finite Element Method. From the displacement field solution, we extract (i) the COD, defined as the $y$-direction displacement jump between the green and blue dots in Fig. 10, (ii) the ERR $G$ by applying the $G - \theta$ method [49], and the SIFs $K_I$, $K_{II}$ by using interaction integrals to decompose the different fracture modes [50, 51]. In practice, the geometry is meshed using the open-source mesh generator Gmsh [52], and the FEM software Cast3M [53] is used to solve the elasticity problem and calculate $K_I$, $K_{II}$, and $G$. Details and validation of the FEM method are provided in Appendix B.

Thanks to these calculations, we compute the following dimensionless quantities:

$$ \hat{S}(\varphi, \hat{s}) \equiv \frac{1}{E} \frac{F}{e \text{COD}} \qquad (6) $$

$$ \hat{K}_p(\varphi, \hat{s}) \equiv K_p \frac{e\sqrt{a_0}}{F}, \quad p = I, II \qquad (7) $$

$$ \hat{G}(\varphi, \hat{s}) \equiv G \frac{Ee^2 a_0}{F^2} \qquad (8) $$

for each value of $\varphi$ and $\hat{s} = s/a_0$.

Advantage of applying a force, as sketched in Fig. 10, instead of upper bolt displacement and of considering a stiffness $\hat{S}$ based on COD and not the pin displacement is that, thanks to Saint-Venant's principle [54], $\hat{S}$, $\hat{K}_I$, $\hat{K}_{II}$, and $\hat{G}$ remain unaffected by how the force is applied in the numerical simulation whereas the results is sensitive to the way the displacement would be modeled at the pins [45].

### 3.1.2. Stress Intensity Factors evolution across a kink

To correctly predict the kink initiation, the value of $\hat{G}^*(\varphi) \equiv \lim_{\hat{s} \to 0} \hat{G}(\varphi, \hat{s})$ is required. However, at the kink position, the SIFs are not continuous [55] and an accurate estimation by FEM is impossible. Instead we use the Amestoy and Leblond [48] formula. It provides a

relationship between SIFs before and for a vanishing kink, assuming that the loading remains unchanged. It can thus be written on the dimensionless stress intensity factors $\hat{K}_p$ since they preserve the loading by construction:

$$ \hat{K}_p^*(\varphi, s) = F_{pq}(\varphi) \hat{K}_q^0 + G_p(\varphi) T\sqrt{s} + O(s) \qquad (9) $$

where $\hat{K}_p^*(\varphi) = \lim_{\hat{s} \to 0} \hat{K}_p(\varphi, \hat{s})$, $\hat{K}_q^0$ are the initial dimensionless SIF without a kinked extension, and $T$ denotes the T-stress. The function $F_{pq}(\varphi)$ and $G_p(\varphi)$ are
313 universal, independent of the elasticity problem, and their values can be found in [48].

315 By substituting the SIFs after kinking $\widehat{K}_p^*(\varphi)$ into Irwin's formula [56], we can obtain the energy release rate for all potential propagation directions:

$$ \hat{G}^*(\varphi) = \frac{1}{E'} (\widehat{K}_I^*(\varphi)^2 + \widehat{K}_{II}^*(\varphi)^2) = \frac{1}{E'} F_{pq}(\varphi) F_{pr}(\varphi) \widehat{K}_q^0 \widehat{K}_r^0 \quad (10) $$

317 where Einstein summation notation is used in last expression, and $E' = 1$ holds for plane stress (P-$\sigma$) assumption while $E' = \frac{1}{1-\nu^2}$ for plane strain (P-$\epsilon$) assumption.

## 3.2. The Generalized Maximum Energy Release Rate criterion

### 3.2.1. At initiation

320 Denote $G^*(\varphi) = \lim_{s\to 0} G(\varphi, s)$. The propagation direction and force that initiate propagation can be determined using the GMERR criterion [31, 30]. For an increasing loading yielding a growing $G^*(\varphi)$, it states that:

324 • as long as $G^*(\varphi) < \Gamma(\varphi)$, $\forall \varphi$ no propagation occurs;

325 • as soon as $\exists \varphi^*$, such as $G^*(\varphi^*) = \Gamma(\varphi^*)$ propagation becomes possible.

<figure><img src="image_12.png" /><figcaption>Figure 11: (Color online) Polar plot of $\Gamma_0/\Gamma(\varphi)$ versus $\Gamma_0/G^*(\varphi)$ for (A) CT 25° raster angle and (B) CTS sample 15° loading angle, 45° raster angle. The fracture energy is taken as $\Gamma(\theta) = \Gamma_0 (1 + \frac{\Delta\Gamma}{\Gamma_0} |\sin 2(\theta - \theta_0)|)$ with $\Delta\Gamma/\Gamma_0=1$.</figcaption></figure>
A convenient way to reformulate it, is to use Frank plot [9] (Fig. 11) to put $1/G^*$ into balance with $1/\Gamma$ or, equivalently, under dimensionless form $\Gamma_0/G^*$ with $\Gamma_0/\Gamma$. Increasing the loading from $G^* = 0$, $\Gamma_0/G^*$ decreases from $\infty$ toward 0. As long as the loading is small enough, $\Gamma_0/\Gamma < \Gamma_0/G^*$, that is $G^*(\varphi) < \Gamma(\varphi)$ in all directions $\varphi$: the energy released by crack propagation is too weak to promote any crack propagation. As soon as $\Gamma_0/G^*$ meets $\Gamma_0/\Gamma$ at one point $\varphi^*$, propagation may start in this direction and the corresponding critical load $F_c$ verifies $G^*(\varphi^*) = \Gamma(\varphi^*)$. Using the dimensional quantities of Eq. (8), we obtain

$$ \varphi^* = \operatorname{argmax}_{\varphi} \left[ \frac{\hat{G}^*(\varphi)\Gamma_0}{\Gamma(\varphi)} \right] \qquad (11) $$

$$ F_c = \left( E e^2 a_0 \frac{\Gamma(\varphi^*)}{\hat{G}^*(\varphi^*)} \right)^{1/2} \qquad (12) $$

By employing Eq. (10) and the FEM to calculate the initial dimensionless Stress Intensity Factor $\hat{K}_p^0$, we can determine $\varphi^*$ and $F_c$ for a given fracture anisotropy $\Gamma(\varphi)$.

### 3.2.2. Influence of the anisotropic degree on the kink angle

Consider the surface energy of Eq. (4). Applying the GMERR using Eq. (11) together with Eq. (10), we get the kink angle $\varphi^*$ as function of the anisotropic degree $\Delta\Gamma/\Gamma_0$, the raster angle $\theta_0$ and the mode mixity ratio $K_{II}^0/K_I^0$.

Fig. 12 shows the predicted kink angle as a function of $\Delta\Gamma/\Gamma_0$ for the raster angle $\theta_0 = 45^\circ$ and different values of $K_{II}^0/K_I^0$ (corresponding to the range of loading angles $\beta = 0-45^\circ$). In pure mode I, as shown by the solid line, the crack propagation follows the direction of $0^\circ$ when the anisotropic level is low and switches to the weak plane $\theta_0$ when $\Delta\Gamma/\Gamma_0 > \Delta\Gamma/\Gamma_0)_c$. Such a threshold above which $\varphi^* = \theta_0$, can be defined for all $\theta_0$ and $K_{II}^0/K_I^0$. Fig. 13 gives its evolution with $K_{II}^0/K_I^0$ for different values of $\theta_0$. For $\theta_0 = 45^\circ$, it decreases with $K_{II}^0/K_I^0$ in coherence with Fig. 12.

This plot shows that if $\Delta\Gamma/\Gamma_0 \gtrsim 0.4$, $\varphi^* = \theta_0$ whatever the raster angle $\theta_0$ and the $K$-ratio. In other words, if the anisotropic degree $\Delta\Gamma/\Gamma_0$ exceeds 0.4, the predicted kink angle will always align with the closest weakest plane, that is be guided by the weakest fracture cost, regardless of the loading condition or printing direction. And for lower values, the outcome of the GMERR depends on both the elastic and fracture energies, that is not trivially on $\Delta\Gamma/\Gamma_0$, $\theta_0$ and $K_{II}^0/K_I^0$.
<figure><img src="image_13.png" /><figcaption>Figure 12: (Color online) Kink angle $\phi^*$ as a function of $\Delta\Gamma/\Gamma_0$ predicted by GMERR criterion for $\theta_0=45^\circ$, and different mode mixity ratios $K_{II}^0/K_I^0$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: (Color online) The minimum critical value of $\Delta\Gamma/\Gamma_0$ required for crack propagation in the direction of the nearest weakest plane with respect to K-ratio $|K_{II}^0/K_I^0|$ for different printing raster angles $\theta_0$.</figcaption></figure>
### 3.2.3. Quasistatic propagation phase

Quasistatic propagation along the direction $\varphi$ implies that the energy balance:

$$G(\varphi, s) = \Gamma(\varphi) \qquad (13)$$

is satisfied during propagation, that is when $s$ increases.

This equation gives the evolution of $F(\varphi, s)$ required to stay quasistatic along this propagation direction. Using Eq. (8), we obtain:

$$F(\varphi, s) = \left( E e^2 a_0 \frac{\Gamma(\varphi)}{\hat{G}(\varphi, s/a_0)} \right)^{1/2} \qquad (14)$$

Notice that in the limit $s \to 0$, the critical initiation force given by Eq. (12) is retrieved.

From the stiffness, we obtain the corresponding equilibrium $COD(\varphi, s)$ value using Eq. (6):

$$COD(\varphi, s) = \frac{F(\varphi, s)}{eE \hat{S}(\varphi, s/a_0)} \qquad (15)$$

Some years ago, Chambolle et al. [57] stated in their abstract that “2D crack kinking in an isotropic setting is incompatible with continuity in time of the propagation”. Does it imply that, given the presence of a kink, there must be a jump from one quasi-static position to another? The answer is not necessarily affirmative. This statement requires context and is only valid if a kink occurs while $G = G_c$ before the kink. In this case, immediately after the kink $G^* \neq G$ (Eq. 10), that is $G^* \neq G_c$, thus necessitating an abrupt readjustment of the load, that is a snapback to retrieve $G = G_c$. In our case, we start with $G < G_c$ and increase the load until $G^*(\varphi) = G_c$. Under these conditions, there is no reason to expect a sudden jump in propagation once the crack kinks in a new direction. By a matter of fact, we will show numerically that our setup induces no snap-back even in presence of a kink.

## 4. Comparison between Linear Elastic Fracture Mechanics theory and experiments

### 4.1. Two-dimensional plane stress or plane strain hypotheses?

To assess whether the plane stress or plane strain assumption is more suitable for our 2D numerical model in representing the experimental samples, we conducted a comparative analysis of the SIFs between experimental and theoretical data, following the methodology outlined by Triclot et al. [45]. The calculation of the SIFs for standard CT tests was determined using the ASTM-1820 [58] norm.

$$K_I = \frac{F}{e\sqrt{W}} \frac{\left(2 + \frac{a}{W}\right) \left(0.886 + 4.64\frac{a}{W} - 13.32\left(\frac{a}{W}\right)^2 + 14.72\left(\frac{a}{W}\right)^3 - 5.6\left(\frac{a}{W}\right)^4\right)}{\left(1 - \frac{a}{W}\right)^{3/2}} \qquad (16)$$
In accordance with the provided equation, where *F* represents the applied load, *e* denotes the thickness, and *W* represents the width of the CT sample measured from the load line to the right edge, the fracture toughness *K*<sub>IC</sub> is determined by substituting the critical load *F*<sub>c</sub>. This result is independent of the Plane Strain versus Plane Stress hypotheses [59].

For the experimental investigations, three CT samples with a raster angle of 0° were used to ensure mode I loading during crack propagation. Fig. 14 compares the SIF $K_I$ obtained through DIC under different 2D assumptions with the average fracture toughness $K_{IC}$ (represented by the grey zone, which accounts for the standard deviation). Hollow and solid points represent results based on the P-$\sigma$ assumption and the P-$\epsilon$ assumption, respectively. The average fracture toughness $K_{IC}$ was calculated using Eq. (16) for these three CT samples.

Observations reveal that during the crack propagation phase, the measured $K_I$ obtained through DIC under P-$\sigma$ assumption closely matches the computed fracture toughness $K_{Ic}$ by substituting the critical load in the norm Eq. (16).

Thus, the P-$\sigma$ assumption was adopted for further analysis.

<figure><img src="image_15.png" /><figcaption>Figure 14: $K_I$ value under different 2D assumptions: Plane stress (P-$\sigma$) and Plane strain (P-$\epsilon$)</figcaption></figure>

## 4.2. *Stress Intensity Factors evolution*

Figures 15-17 compare the SIFs evolution obtained by FEM calculations following the method explained in Section 3.1.1, with those obtained from the experiments using the method of Section 2.5.2. While the FEM SIFs are independent of the Young modulus since force load is applied (see Eq. 7), DIC SIFs are obtained from the measured displacement field by Eq. 3 using the Young modulus measured on the linear part of the corresponding F-COD curves (see Appendix C for the method).
405
406
407
408
409
410

In the linear part, the agreement is excellent in coherence with the results reported in [45]. Concerning the crack propagation phase, the difference is more significant and comes mostly from the inaccurate numerical modeling of the precise unilateral contact (possibly with friction) boundary conditions at the pin hole. Indeed, as discussed in [45], taking the boundary conditions issued from DIC instead of a point force load as done here, eliminates this biais.

<figure><img src="image_16.png" /><figcaption>Figure 15: Numerical and experimental evolution of SIF values in relation to COD for CT sample with 0° raster angle.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 16: Numerical and experimental evolution of SIF values in relation to COD for CT sample with 25° raster angle.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 17: Numerical and experimental evolution of SIF values in relation to COD for CTS sample with 45° raster angle and 15° loading angle.</figcaption></figure>

### 4.3. Application of the GMERR

In the experiments, the propagation direction is always in the weak direction:

411 $\varphi^* = \theta_0$ whatever $\theta_0$ and $\beta$, that is $K_{II}^0/K_I^0$. Figure 13 implies that for the fracture
costs of Eq. (4), we have to choose $\Delta\Gamma/\Gamma_0)_c$ above 0.4 to get the correct propagation direction in all situations. From initiation and during propagation, the critical force $F$ corresponding to $G = G_c$ is then independent of $\Delta\Gamma/\Gamma_0$ and given by Eq. 14 with $\varphi = \theta_0$:

$$F(\theta_0, s) = \left( E e^2 a_0 \frac{\Gamma_0}{\hat{G}(\theta_0, s/a_0)} \right)^{1/2} \qquad (17)$$

### 4.3.1. Force at initiation

The force at initiation can be obtained theoretically using Eq. (17) in the limit $s \to 0$, computing $\hat{G}^*(\theta_0)$ by Eq. (10) using initial SIF $\hat{K}_I^0, \hat{K}_{II}^0$ values obtained by Finite Element (section 3.1.1). Figures 18 and 19 present the critical load against the raster angle, together with the experimental points for CT and CTS samples respectively. The grey shaded area represents theoretical predictions, accounting for variability in the measurements of $E$ and $\Gamma_0$ (table 1).

This comparison shows that the GMERR criterion within the LEFM framework, provides predictions that align fairly with the experimental data in terms of both the trend and magnitude of the initiation force. However, as the force-displacement curves (Fig. 4) departs from linear elasticity at initiation around the peak force, agreement relying on the sole peak force value can only be indicative and more relevant comparison is awaited by studying the whole quasistatic crack propagation phase.
<figure><img src="image_19.png" /><figcaption>Figure 18: Comparison between experimental and theoretical results for CT samples: each black dot corresponds to one experiment. The grey shaded area represents theoretical prediction based on the energy release rate criterion.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 19: Comparison between experimental and theoretical results for CTS samples: each black dot corresponds to one experiment. The grey shaded area represents the theoretical prediction based on the energy release rate criterion.</figcaption></figure>
4.3.2. Quasistatic crack propagation

Beyond the initiation point, it is possible to investigate the whole quasistatic
propagation phase. To make the results independent of the way the boundary con-
ditions are modeled, we follow Triclot *et al.* [45] recommendation and plot at each
instant the force-COD curves (instead of force-displacement) for selected represen-
tative examples: the 3 CT samples with 0° (Fig. 20) and 25° raster angles (Fig. 21),
as well as the 3 CTS samples with 15° loading angle and 45° raster angle (Fig. 22).
Experimentally the force is measured by the force cell and the COD by DIC on the
corresponding numerical image. Theoretically, the force $F(\theta_0, s)$ required to ensure
quasistatic propagation along the propagation direction is determined by Eq. (17)
using $\hat{G}(\theta_0, s/a_0)$ obtained by FEM simulations (section 3.1.1) with fixed sample ge-
ometry (Fig. 2) and successive crack tip positions obtained from DIC measurements
on the image taken at the same moment (section 2.5.2). In parallel, the COD is
obtained by Eq. (15) with $\hat{S}$ computed likewise by FEM. To take account of the dis-
persion of material constants from one sample to the other, Young modulus is taken
by extrapolation of the linear part as explained in Appendix C and $\Gamma_0$ is determined
using the method described in section 2.6.3 sample by sample.

The agreement between the experiments and LEFM is overall good, even excellent for the samples of Fig. 21. The larger discrepancies observed for the other samples may be attributed to printing issues or 3D effects in crack propagation observed during the experiments, that are not captured by the 2D numerical model. For instance, we observed that when a crack propagates along a weak plane after a kink, it may jump to another weak plane. While the fracture surface is generally planar, the layer-by-layer printing process can occasionally result in non-uniform crack front evolution. In such cases, the crack propagates along different weak planes across various layers of the material.

As for the transition between the linear part and the propagation phase, we
emphasize that it can not be obtained accurately by the standard tools of LEFM as
it is affected by a stronger deviation from linear elastic behavior, but also sensitive
to the state of the initial crack tip. This is also why calibrating LEFM based-model
is more relevant on the propagation phase than on the peak force.
<figure><img src="image_21.png" /><figcaption>Figure 20: (Color online) Force-displacement curves for CT samples with 0° raster angle, dashed lines represent numerical results within the LEFM framework</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 21: (Color online) Force-displacement curves for CT samples with 25° raster angle, dashed lines represent numerical results within the LEFM framework</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 22: (Color online) Force-displacement curves for CTS samples, dashed lines represent numerical results within the LEFM framework</figcaption></figure>

## 5. Conclusion and perspectives

In this article, we highlight that the FDM method offers a straightforward means of manipulating the microstructure to assess fracture mechanics models. We investigated crack propagation in FDM printed specimens, under different loading modes, including mode I (CT sample) and mixed mode (CTS sample). A criss-cross filling strategy, in which the printing direction is alternate by 90° between each layer, is chosen. Our fracture experiments revealed that cracks consistently propagated along the nearest weakest plane dictated by the filling strategy. By employing the GMERR criteria, DIC and FEM, we demonstrated that the experiments enter in the framework of LEFM with 2D isotropic elasticity but strong fracture anisotropy.

In further work, we plan to use a more configurable FDM 3D printer to change the degree of anisotropy in printed samples by adjusting various printing parameters. Additionally, we currently investigate the potential of phase-field models with anisotropic surface energy to simulate crack propagation numerically (see Zhai [60] for Li and Maurini [9]).

In the context of damage tolerance approaches, this study on quasistatic crack propagation must be supplemented by a study of fatigue propagation. In this respect, we have made the following observation: the crack propagation path is different under cyclic loading than in quasistatic propagation under monotonic loading [11]. This preliminary observation needs to be supplemented by a systematic experimental campaign (in progress).
Besides being of practical interest to get a better understanding of crack propagation in FDM printed materials, we expect that the experimental benchmark introduced in this study will greatly benefit research on crack propagation in anisotropic materials. In addition to the propagation of cracks in FDM-printed PC samples, this study is intended to be transposed to other anisotropic materials: e.g. wood [12], single crystals [13, 14], Additive Manufactured (AM) materials [15, 16, 17, 61].

### CRediT authorship contribution statement

**X. Zhai:** Writing – original draft, Visualization, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **T. Corre:** Writing – review & editing, Supervision, Validation, Methodology, Investigation, Conceptualization.

**V. Lazarus:** Writing – original draft, Supervision, Resources, Project administration, Methodology, Funding acquisition, Conceptualization.

### Declaration of competing interest

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.

### Data availability

Data will be made available on request.

### Acknowledgement

The research at ENSTA Paris was partially supported by Agence de l'Innovation de Défense – AID - via Centre Interdisciplinaire d'Etudes pour la Défense et la Sécurité – CIEDS (AID-2018 60 0071 00 470 75 01 and FracAddi project 2022), by Agence Nationale de la Recherche (ANR-23-CE51-0054 3FAM) and by IP Paris Doctoral Allocation. The authors thank Stella Brach and Andrès Léon Baldelli for scientific discussions, and David Roucou, Nicolas Thurieau and Lahcène Cherfa, for their help in the experimental setup.

For the purpose of Open Access, a CC-BY public copyright licence (https://creativecommons.org/licenses/by/4.0/) 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.
## 514 Appendix A. Elastic properties

The elastic properties of printed samples were determined using 15 tensile specimens with five raster angles (0°, 10°, 15°, 25° and 45°). The tensile test was monitored using the same equipment as previously described, with a displacement rate of 0.5 mm/min. As noted in Corre and Lazarus [39], tests conducted at speeds up to 3 mm/min showed no significant effect on the material properties.

Strain components ($\epsilon_{xx}$, $\epsilon_{yy}$ and $\epsilon_{xy}$) were extracted from a rectangle region centered on the tensile specimens, analogous to attaching a strain gauge to the sample surface, as shown in fig. A.23 (left). The Young's modulus was determined by performing a linear fit of the data up to 0.8% longitudinal deformation $\epsilon_{yy}$. For Poisson's ratio, the transverse strain $\epsilon_{xx}$ was plotted against longitudinal strain $\epsilon_{yy}$. Fig. A.23 (right) shows the results of tensile tests with different raster angle obtained by linear fit of experimental data.

As reported in Corre and Lazarus [39], the printed samples exhibited nearly isotropic behavior in 2D elasticity, with Young's modulus value largely unaffected by the raster angle. However, Poisson's ratio showed greater variability in measurements, possibly due to surface effects from inadequate printing of the first layer.

The measured Young's modulus *E* and Poisson ratio *ν* for the tensile sample were 1890 MPa (standard deviation 111 MPa) and 0.34 (standard deviation 0.07), respectively. The critical stress $\sigma_c$, was determined by measuring the normal stress at the point of rupture and was found to be 42 MPa with a standard deviation of 2 MPa. These results aligns with the values reported in Corre and Lazarus [39] and the material properties provided by the supplier.
<figure><img src="image_24.png" /></figure>

Figure A.23: (Color online) Left shows the measured displacement field along the y-axis (during the linear phase) of the tensile sample with a 45° raster angle. The rectangular area, positioned at the center of the tensile specimen, represents the designated extraction zone for obtaining strain results. Right shows results for tensile test: Young's modulus *E* and critical stress *σ*<sub>c</sub> with respect to the raster angle

## Appendix B. Finite element methods and validation

The aim of this appendix is to give more details and to validate the FE method used to get the elasticity solution (section 3.1.1) and notably $\hat{G}(s, \varphi)$ used in the critical force determination (Eq. 17).

### Appendix B.1. FEM method

To model a kinked crack path, we meshed a full-sized CT sample using Gmsh, a finite-element mesh generator as shown in Figure B.24 for a 25° kink angle.

Figure (A) shows the entire CT sample mesh, which consists of second-order Triangle elements (blue) and second-order Quadrilateral elements near the crack tip (red). Figure (B) zooms in on the mesh near the crack tip. The circular region around the crack tip consists of uniform quadrilateral elements arranged in multiple layers. The radius $R_1$ of this circular region was chosen based on the kink length *s* ensuring that the angular point was outside this region.

The CT specimen was loaded through the bolt in the pin hole in the fracture experiment. A common approach for modeling this contact problem involves filling
<figure><img src="image_25.png" /><figcaption>Figure B.24: (Color online) Mesh generated by Gmsh, CT sample with a kink angle of 25°. (A) shows the mesh of the entire CT sample with 50000 elements, where the blue and red mesh represent the Tri6 elements and Qua8 elements respectively. (B) is a zoomed-in view of the kink region, which consists of 6000 elements in the near crack tip circle. The near crack tip circle has 30 layers, and the kink length is dimensionless with a value of $\frac{s}{L} = 0.06$</figcaption></figure>

the pin holes with Q4 elements [62, 45]. The bolts were assigned the same material properties as the CT sample. We applied a displacement $\delta$ in the y-direction to the pivot of the upper bolt, while the pivot of the lower bolt was fixed.

The finite element analyses were conducted using the open-source Cast3M software. The calculations solved a linear elastic boundary value problem with a fixed crack geometry. The outputs consisted of the displacement field $u$, strain field $\epsilon$, and stress field $\sigma$.

Post-processing using Cast3M Gtheta procedure based on $G - \theta$ method [49] and interaction integrals [50, 51] gives us the values of the ERR and the SIFs. We conducted validation tests to assess numerical convergence concerning the number of elements and layers where the $G - \theta$ method is applied. The convergence analysis for CT samples, both with and without a crack kink, indicated that convergence was achieved with 20000 elements. We conducted simulations with approximately 30 layers near the crack tip and varied the number of layers $n$ in which the $G - \theta$ method was applied. Convergence was notably achieved when $n > 5$. We did not employ specific element like quarter-point elements around the crack tip, due to the fact that the $G - \theta$ method is not sensitive to the element type at the crack tip.

## Appendix B.2. Validation

To assess our numerical procedure, we solve the elasticity problems for different kink length $s$ from 0.01 to 0.5 mm and kink angle $\phi$ from 15° to 65°. We apply a constant load and compare the SIF variation in respect to $s$ with Amestoy-Leblond
formula (Eq. 9). In this formula, we use $\hat{K}_p^0$ and $T$ obtained from a simulation without a kink extension ($s = 0$). Fig. B.25, where dots represent dimensionless SIFs from numerical calculations and lines represent theoretical values, show an excellent agreement for $s \to 0$. As the experimentally observed range of 0.2 to 1.5 mm is larger than 0.01 mm, it validates the values of $\hat{G}(s, \varphi)$ used to determine the critical forces in Eq. 17.

<figure><img src="image_26.png" /><figcaption>Figure B.25: (Color online) Numerical SIFs $K_p(\varphi, s)\sqrt{L}/E\delta$ against kink length $s$ (represented by dots) are compared with theoretical values (represented by lines), obtained using eq. 9 for various kink angles $\varphi$</figcaption></figure>

## Appendix C. Young modulus of the CT and CTS samples

To verify that the Young's modulus of the CT and CTS samples is consistent with the values obtained from tensile tests (as reported in Table 1), we calculated the Young's modulus of the printed CT/CTS samples as follows:

$$E = \frac{1}{e \hat{S}^0} \left( \frac{F}{\text{COD}} \right)_{\text{exp}} \qquad (C.1)$$

$(F/\text{COD})_{\text{exp}}$ was determined by fitting the linear part of the experimental $F$-COD curves reported in Figs. 20-22. The initial stiffness of the sample without a kink $\hat{S}^0$ was computed using the FEM model.
Fig. C.26 presents the Young's modulus values of several printed CT/CTS samples, along with the value obtained from the tensile test. The red dashed line represents the average Young's modulus obtained from the tensile test results, while the grey shaded area indicates the standard deviation. The results show that the Young's modulus of the printed samples, calculated using Eq. C.1, closely matches the value measured in the tensile tests.

<figure><img src="image_27.png" /><figcaption>Figure C.26: Validation of Young's modulus for different printed CT/CTS samples.</figcaption></figure>

## References

[1] H. Tada, P. C. Paris, G. R. Irwin, The Stress Analysis of Cracks Handbook, Del Research Corporation, Hellertown, USA, 1973. doi:10.1115/1.801535.

[2] V. Lazarus, MEALOR II Damage Mechanics and Local Approach to Fracture, Zenodo, 2023. doi:10.5281/zenodo.10125170.

[3] M. Lebihain, M. Vasoya, V. Lazarus, Size effects in the toughening of brittle materials by heterogeneities: a non-linear analysis of front deformations, International Journal of Solids and Structures 280 (2023). doi:10.1016/j.ijsolstr.2023.112392.

[4] M. Wang, M. Adda-Bedia, J. M. Kolinski, J. Fineberg, How hidden 3d structure within crack fronts reveals energy balance, Journal of the Mechanics and Physics of Solids 161 (2022) 104795. doi:10.1016/j.jmps.2022.104795.
604
605
606
607

[5] G. Molnár, A. Doitrand, V. Lazarus, Phase-field simulation and coupled cri-
terion link echelon cracks to internal length in antiplane shear, Journal of the
Mechanics and Physics of Solids 188 (2024) 105675. doi:10.1016/j.jmps.2024.
105675.

608
609
610
611

[6] A. Vasudevan, Y. Lubomirsky, C.-H. Chen, E. Bouchbinder, A. Karma, Oscilla-
tory and tip-splitting instabilities in 2d dynamic fracture: The roles of intrinsic
material length and time scales, Journal of the Mechanics and Physics of Solids
151 (2021) 104372. doi:10.1016/j.jmps.2021.104372.

612
613
614

[7] S. Patnaik, F. Semperlotti, Variable-order fracture mechanics and its application to dynamic fracture, npj Computational Materials 7 (2021) 27. doi:10.1038/s41524-021-00492-x.

615
616
617

[8] Y. Lubomirsky, E. Bouchbinder, Quenched disorder and instability control dy-
namic fracture in three dimensions, Nature Communications 15 (2024) 7494.
doi:10.1038/s41467-024-51573-6.

618
619
620

[9] B. Li, C. Maurini, Crack kinking in a variational phase-field model of brittle fracture with strongly anisotropic surface energy, Journal of the Mechanics and Physics of Solids 125 (2019) 502-522. doi:10.1016/j.jmps.2019.01.010.

621
622
623

[10] A. Mesgarnejad, C. Pan, R. M. Erb, S. J. Shefelbine, A. Karma, Crack path selection in orientationally ordered composites, Physical Review E 102 (2020) 013004. doi:10.1103/PhysRevE.102.013004. arXiv:1811.05891.

624
625
626

[11] X. Zhai, T. Corre, A. Mesgarnejad, A. Karma, V. Lazarus, Path differences between quasistatic and fatigue cracks in anisotropic media, Physical review E 110 (2024). doi:10.1103/PhysRevE.110.L063001.

627
628
629

[12] K. Ostapska, K. A. Malo, Crack path tracking using dic and xfem modelling of mixed-mode fracture in wood, Theoretical and Applied Fracture Mechanics 112 (2021) 102896. doi:10.1016/j.tafmec.2021.102896.

630
631
632

[13] S. Nagaraja, P. Carrara, L. De Lorenzis, Experimental characterization and phase-field modeling of anisotropic brittle fracture in silicon, Engineering Fracture Mechanics 293 (2023) 109684. doi:10.1016/j.engfracmech.2023.109684.

633
634

[14] J. Cormier, G. Cailletaud, V. Maurel, G. Eggeler (Eds.), Nickel Base Single Crystals Across Length Scales, Elsevier, 2021. doi:0.1016/C2018-0-01723-5.

635
636
637
638

[15] T. D. Ngo, A. Kashani, G. Imbalzano, K. T. Nguyen, D. Hui, Additive manufacturing (3d printing): A review of materials, methods, applications and challenges, Composites Part B: Engineering 143 (2018) 172-196. doi:10.1016/j.compositesb.2018.02.012.
639

640

641

642

[16] J. Allum, A. Moetazedian, A. Gleadall, V. V. Silberschmidt, Interlayer bonding has bulk-material strength in extrusion additive manufacturing: New understanding of anisotropy, Additive Manufacturing 34 (2020) 101297. doi:10.1016/J ADDMA.2020.101297.

643

644

645

646

[17] O. Lampron, A. Lingua, D. Therriault, M. Lévesque, Characterization of the non-isotropic tensile and fracture behavior of unidirectional polylactic acid parts manufactured by material extrusion, Additive Manufacturing 61 (2023) 103369. doi:10.1016/j.addma.2022.103369.

647

648

649

650

[18] B. Bahrami, H. Talebi, M. M. Momeni, M. R. Ayatollahi, Experimental and theoretical investigation of the influence of post-curing on mixed mode fracture properties of 3d-printed polymer samples, Scientific Reports 14 (2024) 13216. doi:10.1038/s41598-024-64136-y.

651

652

653

654

[19] M. R. Khosravani, H. Sadeghian, M. R. Ayatollahi, T. Reinicke, Optimization of fracture toughness in 3D-printed parts: Experiments and numerical simulations, Composite Structures 329 (2024) 117766. doi:10.1016/j.compstruct.2023.117766.

655

656

657

[20] Y. Wang, W. Wang, B. Zhang, C.-Q. Li, A review on mixed mode fracture of metals, Engineering Fracture Mechanics 235 (2020) 107126. doi:10.1016/j.engfracmech.2020.107126.

658

659

660

[21] G. A. Francfort, J. J. Marigo, Revisiting brittle fracture as an energy minimization problem, Journal of the Mechanics and Physics of Solids 46 (1998) 1319-1342. doi:10.1016/S0022-5096(98)00034-9.

661

662

[22] B. Bourdin, G. Francfort, J.-J. Marigo, The variational approach to fracture, Journal of elasticity 91 (2008) 5-148. doi:10.1007/s10659-007-9107-3.

663

664

665

666

[23] C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations, International Journal for Numerical Methods in Engineering 83 (2010) 1273-1311. doi:10.1002/NME.2861.

667

668

669

[24] J.-J. Marigo, C. Maurini, K. Pham, An overview of the modelling of fracture by gradient damage models, Meccanica 51 (2016) 3107-3128. doi:10.1007/s11012-016-0538-4.

670

671

672

[25] B. Bourdin, G. A. Francfort, J. J. Marigo, Numerical experiments in revisited brittle fracture, Journal of the Mechanics and Physics of Solids 48 (2000) 797-826. doi:10.1016/S0022-5096(99)00028-9.
673

674

675

[26] E. Tanné, T. Li, B. Bourdin, J. J. Marigo, C. Maurini, Crack nucleation in variational phase-field models of brittle fracture, Journal of the Mechanics and Physics of Solids 110 (2018) 80–99. doi:10.1016/j.jmps.2017.09.006.

676

677

[27] M. Marder, Cracks cleave crystals, Europhysics Letters 66 (2004) 364. doi:10.1209/ep1/i2003-10254-4.

678

679

680

[28] V. Hakim, A. Karma, Crack path prediction in anisotropic brittle materials, Physical Review Letters 95 (2005) 235501. doi:10.1103/PhysRevLett.95.235501. arXiv:0412775.

681

682

683

[29] V. Hakim, A. Karma, Laws of crack motion and phase-field models of fracture, Journal of the Mechanics and Physics of Solids 57 (2009) 342 – 368. doi:10.1016/j.jmps.2008.10.012.

684

685

686

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

687

688

689

[31] A. Chambolle, G. A. Francfort, J. J. Marigo, When and how do cracks propagate?, Journal of the Mechanics and Physics of Solids 57 (2009) 1614–1622. doi:10.1016/j.jmps.2009.05.009.

690

691

692

693

[32] B. Li, C. Peco, D. Millán, I. Arias, M. Arroyo, Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy, International Journal for Numerical Methods in Engineering 102 (2015) 711–727. doi:10.1002/nme.4726.

694

695

696

[33] S. Teichtmeister, D. Kienle, F. Aldakheel, M. A. Keip, Phase field modeling of fracture in anisotropic brittle solids, International Journal of Non-Linear Mechanics 97 (2017) 1–21. doi:10.1016/j.ijnonlinmec.2017.06.018.

697

698

699

[34] T. T. Nguyen, J. Réthoré, M. C. Baietto, Phase field modelling of anisotropic crack propagation, European Journal of Mechanics, A/Solids 65 (2017) 279–288. doi:10.1016/j.euromechsol.2017.05.002.

700

701

702

[35] J.-M. Scherer, S. Brach, J. Bleyer, An assessment of anisotropic phase-field models of brittle fracture, Computer Methods in Applied Mechanics and Engineering 395 (2022) 115036. doi:10.1016/J.CMA.2022.115036.

703

704

705

706

[36] S. Nagaraja, U. Römer, H. G. Matthies, L. De Lorenzis, Deterministic and stochastic phase-field modeling of anisotropic brittle fracture, Computer Methods in Applied Mechanics and Engineering 408 (2023) 115960. doi:10.1016/J.CMA.2023.115960.
707

708

709

710

[37] O. Lampron, D. Therriault, M. Lévesque, Phase-field modeling of fracture in fused filament fabricated thermoplastic parts and experimental validation, Engineering Fracture Mechanics 292 (2023) 109667. doi:10.1016/j.engfracmech.2023.109667.

711

712

713

[38] A. Ibarra, B. Roman, F. Melo, The tearing path in a thin anisotropic sheet from two pulling points: Wulff's view, Soft Matter 12 (2016) 5979-5985. doi:10.1039/C6SM00734A.

714

715

716

[39] T. Corre, V. Lazarus, Kinked crack paths in polycarbonate samples printed by fused deposition modelling using criss-cross patterns, International Journal of Fracture (2021). doi:10.1007/S10704-021-00518-X.

717

718

719

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

720

[41] J. Réthoré, Ufreckles, UFreckles v2.0 (2018). doi:10.5281/zenodo.1433776.

721

722

723

724

[42] S. Roux, J. Réthoré, 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 (2009). doi:10.1088/0022-3727/42/21/214004.

725

726

727

[43] M. L. Williams, Stress Singularities Resulting From Various Boundary Conditions in Angular Corners of Plates in Extension, Journal of Applied Mechanics 19 (1952) 526-528. doi:10.1115/1.4010553.

728

729

730

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

731

732

733

734

[45] J. Triclot, T. Corre, A. Gravouil, V. Lazarus, Key role of boundary conditions for the 2D modeling of crack propagation in linear elastic Compact Tension tests, Engineering Fracture Mechanics 277 (2023) 109012. doi:10.1016/J.ENGFRACMECH.2022.109012.

735

736

[46] M. F. Ashby, Overview No. 80: On the engineering properties of materials, Acta Metallurgica 37 (1989) 1273-1293. doi:10.1016/0001-6160(89)90158-2.

737

738

739

[47] A. Kim, C. P. Bosnyak, A. Chudnovsky, Effects of weathering, scale, and rate of loading on polycarbonate fracture toughness, Journal of Applied Polymer Science 51 (1994) 1841-1848. doi:10.1002/app.1994.070511101.
740

741

742

[48] M. Amestoy, J. B. Leblond, Crack paths in plane situations-II. Detailed form of the expansion of the stress intensity factors, International Journal of Solids and Structures 29 (1992) 465-501. doi:10.1016/0020-7683(92)90210-K.

743

744

745

746

[49] P. Destuynder, P. E. M. Djaoua, L. Chesnay, J. C. Nedelec, 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 (1981) 70-87. doi:10.1002/MMA.1670030106.

747

748

749

[50] B. Moran, C. F. Shih, Crack tip and associated domain integrals from momentum and energy balance, Engineering Fracture Mechanics 27 (1987) 615-642. doi:10.1016/0013-7944(87)90155-X.

750

751

752

[51] X.-Z. Suo, A. Combescure, On the application of *G-theta* method and its comparison with De Lorenzi's approach, Nuclear Engineering and Design 135 (1992) 207-224. doi:10.1016/0029-5493(92)90223-I.

753

754

[52] Gmsh, Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, https://gmsh.info/, 2024.

755

[53] Cast3M, Accueil | Cast3M, http://www-cast3m.cea.fr/, 2024.

756

757

[54] R. V. Mises, On Saint Venant's principle, Bulletin of the American Mathematical Society 51 (1945) 555-562. doi:10.1090/S0002-9904-1945-08394-3.

758

759

760

[55] J. B. Leblond, Crack paths in plane situations—I. General form of the expansion of the stress intensity factors, International Journal of Solids and Structures 25 (1989) 1311-1325. doi:10.1016/0020-7683(89)90094-2.

761

762

763

[56] G. R. Irwin, Analysis of stresses and strains near the end of a crack traversing a plate, Journal of Applied Mechanics 24 (1957) 361-364. doi:10.1115/1.4011547.

764

765

766

[57] A. Chambolle, G. Francfort, J.-J. Marigo, When and how do cracks propagate?, Journal of the Mechanics and Physics of Solids 57 (2009) 1614-1622. doi:10.1016/j.jmps.2009.05.009.

767

768

[58] ASTM-1820, Standard test method for measurement of fracture toughness, Standard, American Society for Testing and Materials. doi:10.1520/E1820-18.

769

770

771

[59] N. I. Muskhelishvili, Some Basic Problems of the Mathematical Theory of Elasticity, Springer Science & Business Media, 1977. doi:10.1007/978-94-017-3034-1.
772
773
774

[60] X. Zhai, Crack propagation in elastic media with anisotropic fracture toughness : experiments and numerical modeling, Ph.D. thesis, Institut Polytechnique de Paris, 2023. URL: https://hal.science/tel-04528865v1.

775
776
777
778

[61] D. Roucou, T. Corre, G. Rolland, V. Lazarus, Effect of the deposition direction on fracture propagation in a duplex stainless steel manufactured by directed energy deposition, Materials Science and Engineering: A 878 (2023) 145176. doi:https://doi.org/10.1016/j.msea.2023.145176.

779
780
781

[62] A. Ortega, P. Maimí, E. V. González, L. Ripoll, Compact tension specimen for orthotropic materials, Composites Part A: Applied Science and Manufacturing 63 (2014) 85-93. doi:10.1016/J.COMPOSITESA.2014.04.012.