%~Mouliné par MaN_auto v.0.34.0 (b0f8651d) 2025-06-26 09:18:38
\documentclass[CRMECA,Unicode,biblatex,published]{cedram}

\addbibresource{CRMECA_Girard_20250239.bib}

\usepackage{subcaption}
\usepackage{calc}
\usepackage{pgfplots}
\usepackage{siunitx}

\pgfplotsset{compat=newest}
\usepgfplotslibrary{fillbetween}
\pgfplotsset{width=7cm,compat=1.3}
\usepgfplotslibrary{polar}
\usetikzlibrary{plotmarks}
\usetikzlibrary{patterns}
\usepgfplotslibrary{groupplots}
\usetikzlibrary{pgfplots.groupplots}

\newcommand*{\dd}{\mathrm{d}}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\bracks}{[}{]}

%%%---------------------------------------------------------------------------

\graphicspath{./figures/}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}

%%%---------------------------------------------------------------------------

\makeatletter
\def\@setafterauthor{%
  \vglue3mm%
\begingroup\hsize=10.7cm\advance\hsize\abstractmarginL\raggedright
\noindent
   \leftskip\abstractmarginL
  \normalfont\Small
  \@afterauthor\par
\endgroup
\vskip2pt plus 3pt minus 1pt
}
\makeatother

%%%---------------------------------------------------------------------------

\title{Towards the determination of optimal 3D initiation crack shapes}
\alttitle{Vers la détermination des formes optimales de fissures à l'amorçage en 3D}

\author[H. Girard]{\firstname{Hugo} \lastname{Girard}\CDRorcid{0000-0002-0478-5971}\IsCorresp}
\address{Univ. Lyon, INSA Lyon, Université Claude Bernard Lyon 1, CNRS, MATEIS, UMR5510, 69621 Villeurbanne, France}
\email{hugo.girard@insa-lyon.fr}

\author[A. Doitrand]{\firstname{Aurélien} \lastname{Doitrand}\CDRorcid{0000-0003-1000-3507}}
\address[1]{Univ. Lyon, INSA Lyon, Université Claude Bernard Lyon 1, CNRS, MATEIS, UMR5510, 69621 Villeurbanne, France}
\email{aurelien.doitrand@insa-lyon.fr}

\author[M. Tosti Umemura]{\firstname{Milena} \lastname{Tosti Umemura}}
\address[1]{Univ. Lyon, INSA Lyon, Université Claude Bernard Lyon 1, CNRS, MATEIS, UMR5510, 69621 Villeurbanne, France}
\email{milena.umemura@insa-lyon.fr}

\author[B. Koohbor]{\firstname{Behrad} \lastname{Koohbor}\CDRorcid{0000-0002-5787-4644}}
\address{Department of Mechanical Engineering, Rowan University, 201 Mullica Hill Rd., Glassboro, NJ 08028, United States of America}
\email{koohbor@rowan.edu}

\author[R. G. Rinaldi]{\firstname{Renaud} \middlename{G.} \lastname{Rinaldi}\CDRorcid{0000-0001-7066-891X}}
\address[1]{Univ. Lyon, INSA Lyon, Université Claude Bernard Lyon 1, CNRS, MATEIS, UMR5510, 69621 Villeurbanne, France}
\email{renaud.rinaldi@insa-lyon.fr}

\author[N. Godin]{\firstname{Nathalie} \lastname{Godin}\CDRorcid{0000-0003-3813-1594}}
\address[1]{Univ. Lyon, INSA Lyon, Université Claude Bernard Lyon 1, CNRS, MATEIS, UMR5510, 69621 Villeurbanne, France}
\email{nathalie.godin@insa-lyon.fr}

\author[J. Bikard]{\firstname{Jérôme} \lastname{Bikard}\CDRorcid{0000-0003-1364-8415}}
\address{Axel’One (Syensqo), 87 Avenue des Frères Perret, CS 70061, 69192 Saint Fons, France}
\email{jerome.bikard@syensqo.com}

\keywords{\kwd{Coupled criterion} \kwd{crack} \kwd{shape} \kwd{3D} \kwd{fiber-matrix} \kwd{debonding}}
\altkeywords{\kwd{Critère couplé} \kwd{fissure} \kwd{forme} \kwd{3D} \kwd{interface fibre-matrice} \kwd{décohésion}}

\begin{abstract}
An approach is proposed to determine the optimal initiation crack shapes when using the Coupled Criterion~(CC) to assess 3D debonding initiation at a fiber-matrix interface. The optimal crack shapes are determined among potential crack shapes described by their length along the fiber-matrix interface, their angle at the free surface and the shape of the debonding front. The optimal crack shapes actually correspond to the crack shapes maximizing the energy criterion, regardless of the interface fracture properties. Stress-isocontours fail to replicate the debonding shapes that were experimentally observed on the side of a single glass fiber-transparent epoxy sample, unlike the energy-based crack shapes. Both methods yield similar interface strength identification, however, the energy-based shapes result in larger identified critical energy release rates.
\end{abstract}

\begin{altabstract}
Une méthode est proposée pour déterminer les formes optimales de fissures à l'amorçage, en utilisant le Critère Couplé~(CC) pour prédire l'amorçage d’une décohésion à l’interface fibre-matrice en 3D. Les formes optimales de fissures sont déterminées parmi un ensemble de formes potentielles, décrites par leur longueur le long de l’interface fibre-matrice, leur angle à la surface libre et la forme du front de décohésion. Les formes de fissures optimales correspondent en réalité à celles qui maximisent le critère en énergie, indépendamment des propriétés à rupture de l’interface. Les iso-contours de contrainte ne parviennent pas à reproduire les formes de décohésion observées expérimentalement sur le flanc d’un échantillon composé d’une seule fibre de verre dans une matrice d’époxy transparente, contrairement aux formes de fissures basées sur l’énergie. Les deux méthodes permettent une identification similaire de la résistance de l’interface\tralicstex{}{ }; cependant, les formes basées sur l’énergie conduisent à des taux de restitution d’énergie critiques plus élevés.
\end{altabstract}

\dateposted{2025-07-25}
\begin{document}

\maketitle

\section{Introduction}

The Coupled Criterion~(CC), introduced by Leguillon~\cite{leguillon_strength_2002}, is an approach for predicting crack initiation and associated crack size across various configurations~\cite{doitrand_review_2024}. It is particularly efficient under linear elastic assumption, as its solution can be obtained from a single calculation, independent of the material fracture properties. This approach allows for the inverse identification of fracture properties to be performed as a post-processing step~\cite{doitrand_determination_2020}, unlike computationally expensive approaches such as the Cohesive Zone Model~(CZM) \cite{barenblatt_formation_1959,dugdale_yielding_1960} or the Phase-Field~(PF) method for fracture~\cite{bourdin_numerical_2000}. Moreover, the CC provides a deeper understanding of the fracture phenomenon since it involves studying configurations that are not favorable in order to explain the transition to initiation configurations, based on the combined analysis of the energy and stress criteria. Indeed, further analysis is needed to clarify the role of each criterion in the initiation process, whether it is predominantly triggered by energy, stress, or both~\cite{leguillon_prediction_2016}.

Yet, implementing the CC requires the crack path to be pre-determined, which is not necessary with PF, for example. Several possibilities exist for describing the potential initiation crack shapes. For instance, following the stress isocontours which is computationally efficient as it only requires an elastic calculation without cracks~\cite{leguillon_attempt_2014,doitrand_3d_2018}. Mittelman and Yosibash~\cite{mittelman_energy_2015} studied mode III crack initiation and found that the crack path maximizing the energy release rate, along the V-notch bisector, did not correspond to the experimental one, whether the crack shape is circular or elliptical. Instead, they noted that the crack path observed experimentally matched the maximum normal stress location~\cite{yosibash_3_d_2016}. Duminy et al.~\cite{duminy_anisotropic_2023} used stress isocontours to describe interface crack path in nacre-like alumina. In different configurations, stress-isocontours based crack shapes may not be optimal for the energy criterion, which can drive the initiation of the crack, see~\cite{mantic_interface_2009,doitrand_damage_2017,martin_competition_2008,girard_numerical_2023}. Another solution is to assume the shape of the crack, choosing for instance a geometric shape, such as in~\cite{doitrand_damage_2017,garcia_numerical_2016}, to obtain a configurable description of the cracks. Papšík et al.~\cite{papsik_prediction_2023} compared CC solutions of indentation induced-cracking obtained using different crack paths, starting from the contact surface between the indenter and the sample, thus identifying the optimal one. Elliptical cracks were assumed in~\cite{papsik_elip_2023} to study crack initiation in layered ceramics. Nevertheless, the above-mentioned approach does not consider all possible crack geometries. As a consequence, it is yet not established which initiation crack shape among all possible crack shapes is the most favorable, i.e., results in the minimum imposed loading at initiation. In our previous work~\cite{girard_influence_2024}, all possible geometrical debonding shapes were compared using a 2D model consisting of two neighboring fibers in model macro fiber composites. Stress-based and energy-based shape solutions, corresponding respectively to the shapes maximizing either stress or energy conditions, were compared to the optimal CC solution. Stress-based shapes provided the optimal CC solutions for relatively small interface brittleness numbers. For larger brittleness numbers, energy-based shapes yielded the optimal CC solution. There exists a transition zone though, for intermediate brittleness numbers, where neither the stress nor the energy-based shapes were able to provide the optimal CC solution. A similar methodology for determining the optimum initiation crack shape in 3D has yet to be developed.

Determining the crack path in 3D is challenging because the third dimension increases drastically the number of possible crack shapes. To overcome this drawback, García et al.~\cite{garcia_numerical_2016} used configurable rectangles whose dimensions can be varied using two variables. Burhan et al.~\cite{burhan_three_dimensional_2024} investigated configurable semi-elliptical shapes to study interlaminar crack initiation in composite materials, optimal crack shapes were selected as those minimizing the CC solution. Interlaminar strength and critical energy release rate were then identified by minimizing the difference between the CC and experimental test results for both failure stresses and crack initiation length. Leguillon~\cite{leguillon_attempt_2014} attempted to determine the crack shapes at an interface corner based on stress isocontours. Fracture properties identified in 3D were similar to those obtained with 2D simulations. However, Doitrand et al.~\cite{doitrand_damage_2017} found that stress isocontours-based shapes did not fulfill the energy criterion in woven composites when the CC solution was found to be driven solely by the energy condition. In such a configuration, idealized crack surfaces spanning the whole ply thickness and having straight crack fronts were assumed. These surfaces may not be the optimal ones compared to energy-based shapes to describe the 3D potential crack shapes when dealing with energy-driven configurations. With regard to the fiber-matrix debonding configuration studied previously, stress isocontours-based shapes did not allow the identification of shear fracture properties because they are mainly subjected to mode I loading, see Girard et al.~\cite{girard_comparison_2025}. In addition, they provided concave debonding shapes and large debonding angles that were not observed experimentally. In this sense, energy-based shapes or intermediate shapes lying in between stress isocontours-based or energy-based ones could better describe the experimentally observed debonding shapes and allow for the shear fracture properties to be determined. However, their determination is computationally expensive, as an infinite number of possible debonding shapes need to be investigated.

A methodology is therefore presented in this article to efficiently determine the optimal 3D debonding initiation shapes provided that debonding initiates from the fiber pole located at the free surface as observed experimentally. The experimental configuration consisting of a single glass fiber embedded in a transparent epoxy matrix is presented in Section~\ref{sec:exp_sample}, together with the obtained experimental results. Then the CC predictions obtained with the stress-based and energy-based shapes are compared for a wide range of interface fracture properties in Section~\ref{sec:det_opt_shape} to assess the relevance of the stress-based or energy-based shapes in adequately describing the optimal CC solution. Optimal debonding shapes are then compared to the experimental observation. Finally, inverse identification of the fracture properties of the fiber-matrix interface is performed in Section~\ref{sec:INV} using both approaches to quantify the differences they yield.

\section{Experimental characterization of the fiber-matrix debonding}\label{sec:exp_sample}

\subsection{Sample preparation and testing}

The configuration consists in a single glass fiber embedded into an epoxy matrix tensile specimen. The geometry and dimensions of the samples are given in Figure~\ref{single_batch2}.

\begin{figure}[h!]
\begin{subfigure}[c]{0.49\textwidth}
\centering
\newcommand*\svgscale{.9}
\begingroup%
  \makeatletter%
  \providecommand\color[2][]{%
    \errmessage{(Inkscape) Color is used for the text in Inkscape, but the package 'color.sty' is not loaded}%
    \renewcommand\color[2][]{}%
  }%
  \providecommand\transparent[1]{%
    \errmessage{(Inkscape) Transparency is used (non-zero) for the text in Inkscape, but the package 'transparent.sty' is not loaded}%
    \renewcommand\transparent[1]{}%
  }%
  \providecommand\rotatebox[2]{#2}%
  \newcommand*\fsize{\dimexpr\f@size pt\relax}%
  \newcommand*\lineheight[1]{\fontsize{\fsize}{#1\fsize}\selectfont}%
  \ifx\svgwidth\undefined%
    \setlength{\unitlength}{200.35049589bp}%
    \ifx\svgscale\undefined%
      \relax%
    \else%
      \setlength{\unitlength}{\unitlength * \real{\svgscale}}%
    \fi%
  \else%
    \setlength{\unitlength}{\svgwidth}%
  \fi%
  \global\let\svgwidth\undefined%
  \global\let\svgscale\undefined%
  \makeatother%
  \begin{picture}(1,0.69421002)%
    \lineheight{1}%
    \setlength\tabcolsep{0pt}%
    \put(0.22129638,0.46810339){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}Front surface\end{tabular}}}}%
    \put(0.71231633,0.48627317){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}Side surface\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=1]{./figures/single_batch2.pdf}}%
    \put(0.65758052,0.3748182){\color[rgb]{0,0,0}\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}Fiber\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=2]{./figures/single_batch2.pdf}}%
    \put(0.64,0.2065629){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$h = 5$ mm\end{tabular}}}}%
    \put(0.33020888,0.1950131){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}5 mm\end{tabular}}}}%
    \put(0.61239345,0.01546815){\color[rgb]{0,0,0}\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$-\frac{U}{2}$\end{tabular}}}}%
    \put(0.61239345,0.65184974){\color[rgb]{0,0,0}\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$\frac{U}{2}$\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=3]{./figures/single_batch2.pdf}}%
    \put(0.3821865,0.39015747){\color[rgb]{0,0,0}\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$\sigma^{\infty}$\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=4]{./figures/single_batch2.pdf}}%
    \put(0.04074895,0.34418079){\rotatebox{90}{\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}40 mm\end{tabular}}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=5]{./figures/single_batch2.pdf}}%
    \put(0.87577601,0.27983894){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$\mathrm{y}$\end{tabular}}}}%
    \put(0.9636233,0.18323685){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$\mathrm{x}$\end{tabular}}}}%
    \put(0.75757853,0.1593795){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$\mathrm{z}$\end{tabular}}}}%
  \end{picture}%
\endgroup%
\caption{}\label{single_batch2}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.49\textwidth}
\centering
\newcommand*\svgscale{1.}
\begingroup%
  \makeatletter%
  \providecommand\color[2][]{%
    \errmessage{(Inkscape) Color is used for the text in Inkscape, but the package 'color.sty' is not loaded}%
    \renewcommand\color[2][]{}%
  }%
  \providecommand\transparent[1]{%
    \errmessage{(Inkscape) Transparency is used (non-zero) for the text in Inkscape, but the package 'transparent.sty' is not loaded}%
    \renewcommand\transparent[1]{}%
  }%
  \providecommand\rotatebox[2]{#2}%
  \newcommand*\fsize{\dimexpr\f@size pt\relax}%
  \newcommand*\lineheight[1]{\fontsize{\fsize}{#1\fsize}\selectfont}%
  \ifx\svgwidth\undefined%
    \setlength{\unitlength}{132.55549215bp}%
    \ifx\svgscale\undefined%
      \relax%
    \else%
      \setlength{\unitlength}{\unitlength * \real{\svgscale}}%
    \fi%
  \else%
    \setlength{\unitlength}{\svgwidth}%
  \fi%
  \global\let\svgwidth\undefined%
  \global\let\svgscale\undefined%
  \makeatother%
  \begin{picture}(1,1.08072764)%
    \lineheight{1}%
    \setlength\tabcolsep{0pt}%
    \put(0,0){\includegraphics[width=\unitlength,page=1]{./figures/real_sample.pdf}}%
    \put(0.77941111,0.60145114){\color[rgb]{0,0,0}\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}Glass fiber\end{tabular}}}}%
    \put(0.17046609,0.80918667){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}Polished\\surfaces\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=2]{./figures/real_sample.pdf}}%
    \put(0.15585408,0.33958708){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}5 mm\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=3]{./figures/real_sample.pdf}}%
    \put(0.69747022,0.43834498){\color[rgb]{0,0,0}\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}Epoxy\end{tabular}}}}%
  \end{picture}%
\endgroup%
\caption{}\label{real_sample}
\end{subfigure}
\caption{(a) Single fiber specimen geometry, dimensions and testing. (b) As prepared sample with the glass fiber embedded into an epoxy matrix.}\label{fig:sample}
\end{figure}

The isotropic elastic properties of the epoxy matrix (SD2505/SR1500, Sicomin, France) are determined in-house, using 5 cured samples with the same epoxy mixture as the single-fiber samples. Tensile tests are carried out on the samples at room temperature under quasi-static loading with a cross-head speed of 0.18 \unit{mm/min}. The average elastic properties obtained from the 5 samples are given in Table~\ref{tab:prop}, together with those of the 1 \unit{mm} glass fiber supplied by the manufacturer (quartz glass fiber, McMaster-Carr, USA).

\begin{table}[h!]
\caption{Linear elastic properties of the epoxy matrix and glass fiber.}\label{tab:prop}
\begin{tabular}{|c|c|c|}
\hline
&&\\*[-2.7ex]
\hline
\textbf{Properties} & \textbf{Epoxy} & \textbf{Glass fiber} \\
\hline
&&\\*[-2.7ex]
\hline
$E$ [\unit{GPa}] & 2.7 $\pm$ 0.2 & 72 \\
\hline
$\nu$ & 0.41 $\pm$ 0.04 & 0.17 \\
\hline
\end{tabular}
\end{table}

Three single fiber samples are fabricated by embedding a 1 \unit{mm} diameter glass fiber in an epoxy matrix, denoted sample 1, 2 and 3. The dimensions of the samples are defined for use in a micro tensile testing machine (see Figure~\ref{single_batch2}). The fiber is held in place at the intersection of the 2 symmetry axes using a 3D printed support. The surface of the fiber is cleaned with isopropanol. The matrix is then poured into a silicone mold, where the mold is drilled at the location of the fiber to allow it to slightly protrude from the final sample. This eliminates the need to polish the end of the fiber beforehand. Curing time is one month at room temperature. After curing, the front surface of the sample is lightly polished with fine sandpaper (10 \unit{\micro\metre} grit). The side surfaces are also polished, first with fine sandpaper (8 \unit{\micro\metre} grit), then using a rotating cotton disc coupled with a polishing solution to obtain a translucent and smooth side surface. A speckle pattern is then applied to the front surface to enable 2D Digital Image Correlation (DIC). A dual-vision system consisting of two perpendicularly positioned cameras (FLIR Grasshopper GS3-U3-41C6M-C) is installed~\cite{uddin_simultaneous_2024}. The front camera allows DIC close to the fiber-matrix interface, while the second side camera is used to observe the fiber side.

\subsection{Experimental characterization of the debonding initiation and propagation}

The geometric characteristics that can be extracted from the two cameras after debonding initiation are depicted in Figure~\ref{fig:param_geom}.

\begin{figure}[h!]
\newcommand*\svgscale{1.1}
\begingroup%
  \makeatletter%
  \providecommand\color[2][]{%
    \errmessage{(Inkscape) Color is used for the text in Inkscape, but the package 'color.sty' is not loaded}%
    \renewcommand\color[2][]{}%
  }%
  \providecommand\transparent[1]{%
    \errmessage{(Inkscape) Transparency is used (non-zero) for the text in Inkscape, but the package 'transparent.sty' is not loaded}%
    \renewcommand\transparent[1]{}%
  }%
  \providecommand\rotatebox[2]{#2}%
  \newcommand*\fsize{\dimexpr\f@size pt\relax}%
  \newcommand*\lineheight[1]{\fontsize{\fsize}{#1\fsize}\selectfont}%
  \ifx\svgwidth\undefined%
    \setlength{\unitlength}{148.47850003bp}%
    \ifx\svgscale\undefined%
      \relax%
    \else%
      \setlength{\unitlength}{\unitlength * \real{\svgscale}}%
    \fi%
  \else%
    \setlength{\unitlength}{\svgwidth}%
  \fi%
  \global\let\svgwidth\undefined%
  \global\let\svgscale\undefined%
  \makeatother%
  \begin{picture}(1,0.604431)%
    \lineheight{1}%
    \setlength\tabcolsep{0pt}%
    \put(0,0){\includegraphics[width=\unitlength,page=1]{./figures/param.pdf}}%
    \put(0.6619271,0.29468318){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$h$\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=2]{./figures/param.pdf}}%
    \put(0.3775722,0.53912818){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$\ell_\mathrm{d}$\end{tabular}}}}%
    \put(0.28662991,0.39972342){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$\delta_\mathrm{nn}$\end{tabular}}}}%
    \put(0.30595191,0.33989623){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$\theta_\mathrm{d}$\end{tabular}}}}%
    \put(0.72405797,0.49035256){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}Debonding\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=3]{./figures/param.pdf}}%
    \put(0.19093693,0.56011065){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}Pole\end{tabular}}}}%
    \put(0.0628363,0.37563355){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}Equator\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=4]{./figures/param.pdf}}%
    \put(0.83237764,0.20833688){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$\mathrm{y}$\end{tabular}}}}%
    \put(0.95091485,0.07798618){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$\mathrm{x}$\end{tabular}}}}%
    \put(0.63443735,0.02614933){\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}$\mathrm{z}$\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=5]{./figures/param.pdf}}%
  \end{picture}%
\endgroup%
\caption{Set of geometrical parameters used to describe the debonding geometry.}\label{fig:param_geom}
\end{figure}

The debonding can be described by its angle $\theta_\mathrm{d}$ at the front surface, corresponding to the angle between the two tips of the debonding, see~\cite{girard_experimental_2023} for more information about its experimental determination using the side camera. The quantity $\delta_\mathrm{nn}$ corresponds to the maximum normal opening of the debonding located at the pole of the fiber with respect to the loading direction~\cite{girard_comparison_2025}. The interface opening reflects the distance between the two lips of the debonding and it can be determined by performing DIC at the front surface. The debonding length through the thickness is described by $\ell_\mathrm{d}$, which corresponds to the maximum length in the z direction at the pole of the fiber~\cite{uddin_simultaneous_2024}. Figure~\ref{fig:exp_SFII-4} illustrates the variation of the different geometric characteristics as a function of the remote loading applied ($\sigma^\infty$) to the sample 3. The results obtained for sample~3 are representative of all three samples tested, using the same approaches as described above.

\begin{figure}[h!]

\pgfplotstableread{OP/BATCH_II/F009_top.dat}{\ninetopp}

\begin{tikzpicture}
\pgfplotsset{set layers}

\begin{axis}[width = 0.5\textwidth,
height = 0.4\textwidth, axis y line*=left, xlabel={Remote loading ($\sigma^\infty$) [\unit{MPa}]}, xmin = 5., xmax = 16., xtick distance = 2., ylabel={Normalized length ($\ell_\mathrm{d}/h$) [-]}, legend style={at={(0.02,0.9)},anchor=west},]
\addplot[green] file[] {LEN/BATCH_II/F009.dat};
\node[rotate=0] at (axis description cs:0.47,0.75) {\rotatebox{0}{Initiation}};
\draw[->] (axis description cs:0.47,0.7) -- (axis description cs:0.47,0.57);
\legend{$\ell_\mathrm{d}/h$}
\end{axis}

\begin{axis}[width = 0.5\textwidth,
height = 0.4\textwidth, axis y line*=right, axis x line=none, xmin = 5., xmax = 16., xtick distance = 2., ytick distance = 3., ylabel={Opening ($\delta_\mathrm{nn}$) {[}\unit{\micro\metre}{]}}, legend style={at={(0.02,0.75)},anchor=west},]
\addplot[dashed,red] table [x = {S}, y = {Op}] {\ninetopp};
\legend{$\delta_\mathrm{nn}$}
\end{axis}

\begin{axis}[width = 0.5\textwidth,
height = 0.4\textwidth, axis y line*=right, axis x line=none, xmin = 5., xmax = 16., xtick distance = 2., ylabel= Half angle ($\theta_\mathrm{d}/2$) {[}deg.{]}, legend style={at={(0.02,0.6)},anchor=west},]
\pgfplotsset{every outer y axis line/.style={xshift=2cm}, every tick/.style={xshift=2cm}, every y tick label/.style={xshift=2cm}}
\addplot[dashdotted,blue] file[] {ANGLE/BATCH_II/F009.dat};
\legend{$\theta_\mathrm{d}$}
\end{axis}
\end{tikzpicture}
\caption{Variation of the geometrical characteristics of the debonding as a function of the remote loading for the sample 3.}\label{fig:exp_SFII-4}
\end{figure}

All the debonding characteristics enable the determination of a similar initiation remote stress around 10 \unit{MPa}. They evidence a step-wise debonding process. Unstable debonding propagation phases are followed by either stable propagation phase or debonding arrest. During the unstable propagation phases, both the debonding length and angle increase, as well as the opening. During the debonding arrest phase, both the length and angle remain constant contrary to the opening which increases with an increase in the tensile loading. During the stable propagation phase, the debonding angle increases while its length remains constant. Table~\ref{tab:exp_results} summarizes the different characteristics magnitude after debonding initiation with the associated remote loading for the 3 samples.

\begin{table}[h]
\caption{Fiber-matrix debonding geometric parameters obtained experimentally from the three samples after initiation.}\label{tab:exp_results}
\begin{tabular}{|c|c|c|c|}
\hline
&&\\*[-2.7ex]
\hline
& \textbf{Sample 1} & \textbf{Sample 2} & \textbf{Sample 3} \\
\hline
&&\\*[-2.7ex]
\hline
$\sigma^\infty$ [\unit{MPa}] & 7.25 $\pm$ 0.04 & 8.03 $\pm$ 0.04 & 10.04 $\pm$ 0.04 \\
\hline
$\theta_\mathrm{d}/2$ [deg.] & 29 & 37 & 40 \\
\hline
$\ell_\mathrm{d}/h$ [-] & 0.03 & 0.02 & 0.05 \\
\hline
$\delta_\mathrm{nn}$ [\unit{\micro\metre}] & 1.5 $\pm$ 1.0 & 2.1 $\pm$ 0.8 & 2.1 $\pm$ 1.0 \\
\hline
\end{tabular}
\end{table}

Accurate debonding topology is also assessed after initiation for one specimen using Computed Tomography (CT), the equipment is manufactured by the Phoenix X-ray company. A voxel size of 5 \unit{\micro\metre} is achieved. The approach is more challenging than using the dual-vision system since it requires adding contrast penetrant to highlight the debonding location. Contrast penetrant similar to that in~\cite{yu_2d_2015} is employed. As a consequence, only one scan is performed, just after debonding initiation observed using the side camera. The loading is maintained constant while the penetrant is added. The debonding observed probably has a larger surface area than that at initiation, having potentially propagated when the penetrant is added. Nevertheless, the shape observed is close to that of the initiation and can therefore still be used as a target for the experimental-numerical dialogue. Figure~\ref{fig:tomo} presents the experimental debonding topology after initiation. The debonding front, indicated by the white line, exhibits a convex shape.

\begin{figure}[h!]
\begingroup%
  \makeatletter%
  \providecommand\color[2][]{%
    \errmessage{(Inkscape) Color is used for the text in Inkscape, but the package 'color.sty' is not loaded}%
    \renewcommand\color[2][]{}%
  }%
  \providecommand\transparent[1]{%
    \errmessage{(Inkscape) Transparency is used (non-zero) for the text in Inkscape, but the package 'transparent.sty' is not loaded}%
    \renewcommand\transparent[1]{}%
  }%
  \providecommand\rotatebox[2]{#2}%
  \newcommand*\fsize{\dimexpr\f@size pt\relax}%
  \newcommand*\lineheight[1]{\fontsize{\fsize}{#1\fsize}\selectfont}%
  \ifx\svgwidth\undefined%
    \setlength{\unitlength}{129.81447204bp}%
    \ifx\svgscale\undefined%
      \relax%
    \else%
      \setlength{\unitlength}{\unitlength * \real{\svgscale}}%
    \fi%
  \else%
    \setlength{\unitlength}{\svgwidth}%
  \fi%
  \global\let\svgwidth\undefined%
  \global\let\svgscale\undefined%
  \makeatother%
  \begin{picture}(1,0.9996557)%
    \lineheight{1}%
    \setlength\tabcolsep{0pt}%
    \put(0,0){\includegraphics[width=\unitlength,page=1]{./figures/tomo.pdf}}%
    \put(0.22841637,0.64182178){\color[rgb]{1,1,1}\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}Matrix\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=2]{./figures/tomo.pdf}}%
    \put(0.6859407,0.72138292){\color[rgb]{1,1,1}\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}Fiber\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=3]{./figures/tomo.pdf}}%
    \put(0.50356942,0.1920117){\color[rgb]{1,1,1}\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}Debonding\end{tabular}}}}%
    \put(0,0){\includegraphics[width=\unitlength,page=4]{./figures/tomo.pdf}}%
    \put(0.85133833,0.42976083){\color[rgb]{1,1,1}\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}x\end{tabular}}}}%
    \put(0.78713805,0.54183943){\color[rgb]{1,1,1}\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}y\end{tabular}}}}%
    \put(0.71532049,0.30027174){\color[rgb]{1,1,1}\makebox(0,0)[c]{\lineheight{1.25}\smash{\begin{tabular}[t]{l}z\end{tabular}}}}%
  \end{picture}%
\endgroup%
\caption{Debonding topology observed after initiation using X-ray tomography. The debonding zone is highlighted using a contrast penetrant while the debonding front is depicted by a white thick line. The volume of the scan is 3 \unit{mm} in the x-direction and y-direction, 5 \unit{mm} in the z-direction.}\label{fig:tomo}
\end{figure}

\section{Determination of the optimal initiation crack shape}\label{sec:det_opt_shape}

This section provides a methodology to determine the optimal initiation debonding shape to be used when implementing the CC.

\subsection{The fiber-matrix interface model}

A 3D model is employed to account for the singularity at the free surface of the sample, caused by the fiber-matrix elastic mismatch. The model is established to compare stress-based and energy-based crack shapes to describe the potential debonding path. The model consists of either 1/4 or 1/8 simplification of the real geometry, using ten-node quadratic elements under the assumption of small strain and linear elasticity. A 1/4 simplification is used for the calculation of energy-based shapes, as it allows a single initiation site to be studied. A second 1/8 simplified model is used with stress-based shapes as they provide a two-site initiation, i.e., at both free surfaces, so that additional symmetry can be added to the model~\cite{girard_comparison_2025}. Two mesh convergence studies are performed since the CC involves both stress and energy conditions.
Overall, a difference of less than 2\% is achieved, compared to a converged solution, on the stress fields at the fiber-matrix interface and the elastic strain energy released for a fixed debonding surface.
A mesh size of 20 \unit{\micro\metre} is thus adopted at the fiber-matrix interface.

\subsection{Implementation of the CC}\label{section_3_2}

The debonding initiation is predicted by means of the CC similar to Manti{\v c}~\cite{mantic_interface_2009} and Girard et al.~\cite{girard_numerical_2023}. The CC combines two conditions to provide the initiation loading, the debonding surface and shape. The energy condition is based on an energy conservation principle before and after crack initiation. It is a necessary condition but it is not sufficient. It is similar to the Griffith criterion but considers an incremental balance instead of a differential one. 
As a result, the variation in the external work forces has to be equal to the variation of other sources of energy such as:
\[
\Delta W_\mathrm{k} + \Delta W_\mathrm{el} + \Delta W_\mathrm{d} + \overline{G_\mathrm{c}} S  = \Delta W_\mathrm{ext},
\label{eq:Energy_balance}
\]
where $\Delta W_\mathrm{ext}$ denotes the external work forces, $\Delta W_\mathrm{k}$ the kinetic energy, $\Delta W_\mathrm{el}$ the elastic strain energy, $\Delta W_\mathrm{d}$ the energy dissipated into other mechanisms such as plasticity, diffused damage or friction for instance and $\overline{G_\mathrm{c}} S$ the crack surface creation energy, where $S$ is the crack surface and $\overline{G_\mathrm{c}}$ is the average critical ERR whose calculation is detailed in Section~\ref{sec:ener_shape}.
A rigorous description of the energy balance in the thermodynamic context in presented in Dormieux et al.~\cite{dormieux_remarks_2025}.
Different thermal effect occurring at crack initiation could also be accounted in the energy balance, for instance local heating zone resulting from the cracking process, sliding, viscous effects.
These effects could be accounted through $\Delta W_\mathrm{d}$.
However, they are difficult to estimate and would require the use of a dedicated numerical model or a thermal camera during the experiments.
As a result, crack initiation will be assumed to be isothermal and all dissipated energy will not be taken into account.
Furthermore, although the dissipated energies are not taken into account in the energy balance, they are in fact accounted for phenomenologically through the critical ERR during its determination, under the assumption that these effects do not vary with the crack surface.
Dormieux et al.~\cite{dormieux_remarks_2025} also proposed rigorous bounds for the loading level at crack initiation and associated observable crack length.
They discussed conditions under which quasi-static crack initiation can be considered.
The validity of the quasi-static crack initiation assumption depends on the variation of the Incremental Energy Release Rate (IERR) $G_\mathrm{inc}$:
\[
G_\mathrm{inc}(\Gamma,\sigma^\infty) = \frac{\Delta W_\mathrm{ext} - \Delta W_\mathrm{el} - \Delta W_\mathrm{k}}{S} = \overline{G_\mathrm{c}},
\]
where $\Gamma$ denotes the debonding path and $\sigma^\infty$ the remote tensile stress.
Noting that the ERR can be obtained from the IERR using:
\[
G = \underbrace{G_\mathrm{inc}}_{=\overline{G_\mathrm{c}}} + S \frac{\mathrm{d} G_\mathrm{inc}}{\mathrm{d} S},
\]
and that $G_\mathrm{inc} = \overline{G_\mathrm{c}}$ at initiation.
So to ensure a stable crack initiation with $G \leq \overline{G_\mathrm{c}}$ at initiation, $\frac{\mathrm{d} G_\mathrm{inc}}{\mathrm{d} S}$ has to be lower or equal to $0$.
This case is detailed in Dormieux et al.~\cite{dormieux_remarks_2025}.
It can be obtained when $G_\mathrm{inc}/\overline{G_\mathrm{c}}$ exhibits a maximum \cite{girard_numerical_2023} or in the presence of strong singularities \cite{aranda_prediction_2023} for which deceasing variations of $G_\mathrm{inc}$ and $G$ as a function of the crack surface are obtained.
In this case, crack initiation can be studied within the framework of a quasi-static hypothesis with neglected kinetic energy.
On the contrary, if $\frac{\mathrm{d} G_\mathrm{inc}}{\mathrm{d} S} > 0$, $G > G_\mathrm{c}$, unstable crack initiation takes place and the quasi-static assumption is no longer valid.
It would be necessary to take kinetic energy into account in the energy balance, see \cite{doitrand_dynamic_2022,chen_crack_2023}.
Quasi-static crack initiation is assumed in this work since it can not be quantified experimentally, so the IERR reverts to:
\[
G_{\mathrm{inc}}(\Gamma,\sigma^\infty) = \frac{W(0,\sigma^\infty)-W(\Gamma,\sigma^\infty)}{S}
\]
with $W$ denotes the elastic strain energy of the model. As a consequence, only stable crack initiation can be considered.

The energy condition solely is generally adopted to study the propagation of an already existing semi-infinite crack in an infinite medium.
The latter needs to be complemented by a stress condition to predict crack initiation.
The stress condition must be satisfied at any location over the entire debonding path $\Gamma$ before initiation, according to the condition proposed in \cite{leguillon_strength_2002}.
This non-local stress condition thus enables its application even if the stress fields are non-monotonic.
A quadratic stress criterion $\sigma_\mathrm{eq}$ is employed to account for the mode mixity, i.e., the normal stress $\sigma_\mathrm{nn}$, the in-plane and out-of-plane shear stresses, respectively, $\sigma_\mathrm{nt}$ and $\sigma_\mathrm{nz}$:
\[
\sigma_\mathrm{eq} (\vec{x},\sigma^\infty) = \sqrt{\sigma_\mathrm{nn}^{2}(\vec{x},\sigma^\infty) + \tfrac{1}{\mu^2} \parens[\big]{ \sigma_\mathrm{nt}^2(\vec{x},\sigma^\infty) + \sigma_\mathrm{nz}^2(\vec{x},\sigma^\infty)}},
\]
where $\vec{x}$ is the location vector and $\mu$ the shear to tensile strength ratio.
As for the energy condition, this criterion is necessary but not sufficient.
These two conditions are combined as follow:
\[
\begin{cases}
\sigma_\mathrm{eq} \geq \sigma_\mathrm{c} 									&\forall \vec{x} \in \Gamma, \\
G_{\mathrm{inc}}(\Gamma,\sigma^\infty) \geq \overline{G_\mathrm{c}}(\Gamma),
\end{cases}
\]
where $\sigma_\mathrm{c}$ is the interface tensile strength.

One of the unknown to be determined in the CC is the initiation crack path. The potential crack path can be chosen in order to minimize the imposed loading that fulfills (i) the stress condition, (ii) the energy condition or (iii) both stress and energy conditions.

\subsection{Stress-based debonding shapes}

Figure~\ref{fig:stress_iso} illustrates the debonding stress-based shapes whose fronts are defined based on the normal stress isocontours, superimposed on one-quarter of the fiber-matrix interface. Further details about the determination and implementation of the stress-isocontours based debonding shapes can be found in Girard et al.~\cite{girard_comparison_2025}.

\begin{figure}[h!]

\pgfplotstableread{CHAP_6/shape_on_demand/stress_surf.dat}{\table}

\begin{tikzpicture}
\begin{axis}[view/az=50,
xlabel={x [mm]}, zlabel={y [mm]}, ylabel={-z [mm]}, y = {(0.5*0.3\textwidth,0.5*0.075\textwidth)}, x = {(0.5*0.2\textwidth,-0.075*0.5\textwidth)}, z = {(0.0\textwidth,0.5*0.315\textwidth)}, domain=0:2.375, y domain=0:pi, xmin=-0.6, xmax=0.6, ymin=0, ymax=2.375, zmin=0.0, zmax = 0.6, samples=30, ytick distance = 0.5, xtick distance = 0.5, ztick distance = 0.5, colorbar horizontal, colorbar style={at={(0.5,1.1)},anchor=south,xticklabel pos=upper,title={$S$ [\unit{mm^2}]},title style={at={(0.5,1.4)}}}, point meta min = {0.0044}, point meta max = {1.7}, colormap/jet,]
\addplot3 [surf,fill=lightgray,faceted color=lightgray, line width = 0.2pt] ({0.5*cos(deg(y))},{x},{0.5*sin(deg(y))});
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt] table [x = {x}, y = {y}, z={z}] {\table};
\end{axis}
\end{tikzpicture}
\caption{Potential debonding shape fronts obtained using normal stress isocontours superimposed on one-quarter of the fiber-matrix interface for varying $S$.}\label{fig:stress_iso}
\end{figure}

Following the stress isocontours before initiation, debonding initiates at the fiber pole with respect to the loading direction and propagates in a concave shape, forming a tunnel for sufficiently large cracks ($S > 0.2$ \unit{mm^2}). Although only one ratio between fiber diameter and sample thickness is investigated in this work, the debonding shapes derived from the stress isocontours will remain similar for fibers that are sufficiently long with respect to their diameter. In the case of a very thin sample with respect to the fiber diameter, the singularities of the two free surfaces will influence the isocontours. Irrespective of the crack surface, the potential debonding shapes given by the stress criterion significantly differ from the one observed using CT (Figure~\ref{fig:tomo}). For small crack surface, the debonding angle given by the stress isocontours is larger than the one observed experimentally, and for larger crack surfaces the debonding shape is concave, contrary to the convex shapes observed experimentally. It can thus be concluded that for the studied configuration, the stress-based crack shapes significantly differ from the crack shapes observed experimentally. This discrepancy suggests the possibility of alternative optimal initiation debonding shapes, which are studied in next sections.

\subsection{Energy-based debonding shapes}\label{sec:ener_shape}

Determining energy-based debonding shapes is challenging because it requires computing the elastic strain energy difference between the states with and without crack. It thus requires one Finite Element (FE) calculation for a given crack shape and there is theoretically an infinite number of 3D crack shapes to be evaluated. The determination of energy-based crack shapes consists in defining the debonding shape that maximizes the IERR to critical ERR ratio for a given crack surface increment. However, it is not possible to test an infinite number of debonding shape configurations. Therefore, assumptions about the potential debonding shapes are made in the following. On the one hand, since the stress criterion reaches a peak at the pole of the fiber at the free surface, it is likely that the debonding emanates from this location, as presented in Figure~\ref{fig:tomo}. On the other hand, although experimental observations show convex shapes, both concave and convex shapes are investigated for comparison purposes.

The potential debonding shapes are therefore approximated using a power-law relationship, involving geometric parameters depicted in Figure~\ref{fig:param_geom}. The location of the debonding front at the free surface is determined by the debonding angle ($\theta_\mathrm{d}$), while the debonding length ($\ell_\mathrm{d}$) determines the location at the pole of the fiber through the thickness. The topology of the debonding front, connecting these two points, can therefore be approximated using a power-law whose exponent $n$ can vary from 0.6 to 4, see~\eqref{shape_approx}, where $r$ is the fiber radius.
\begin{equation}
z(x) = - \parens*{\ell_\mathrm{d}^n \parens*{1 - \parens*{\frac{x}{\sin \parens[\Big]{\frac{\theta_\mathrm{d}}{2}} r}}^n}} ^{\frac{1}{n}}\label{shape_approx}
\end{equation}

Figure~\ref{inf_param_shape} shows the influence of the three parameters on the topology of the debonding front.

\begin{figure}[h!]
\begin{subfigure}[c]{0.45\textwidth}
\centering

\pgfplotstableread{CHAP_6/shape_on_demand/var_theta.dat}{\table}

\begin{tikzpicture}
\begin{axis}[view/az=50,
xlabel={x [mm]}, zlabel={y [mm]}, ylabel={-z [mm]}, y = {(0.3\textwidth,0.075\textwidth)}, x = {(0.3\textwidth,-0.075\textwidth)}, z = {(0.0\textwidth,0.315\textwidth)}, domain=0:1.2, y domain=0:pi*0.5, xmin=0, xmax=0.6, ymin=0, ymax=1.2, zmin=0.0, zmax = 0.6, samples=30, ytick distance = 0.5, xtick distance = 0.5, ztick distance = 0.5, colorbar horizontal, colorbar style={at={(0.5,1.1)},anchor=south,xticklabel pos=upper,title={$\theta_\mathrm{d}/2$ [deg.]},title style={at={(0.5,1.4)}}}, point meta min = {9}, point meta max = {91}, colormap/cool,]
\addplot3 [surf,fill=lightgray,faceted color=lightgray] ({0.5*cos(deg(y))},{x},{0.5*sin(deg(y))});
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt] table [x = {x}, y = {y}, z={z}] {\table};
\end{axis}
\end{tikzpicture}

\caption{$\ell_\mathrm{d} = 1 \ \unit{mm}, \ n = 2$}\label{inf_angle}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.45\textwidth}
\centering

\pgfplotstableread{CHAP_6/shape_on_demand/var_length.dat}{\table}

\begin{tikzpicture}
\begin{axis}[view/az=50,
xlabel={x [mm]}, zlabel={y [mm]}, ylabel={-z [mm]}, y = {(0.3\textwidth,0.075\textwidth)}, x = {(0.3\textwidth,-0.075\textwidth)}, z = {(0.0\textwidth,0.315\textwidth)}, domain=0:1.2, y domain=0:pi*0.5, xmin=0, xmax=0.6, ymin=0, ymax=1.2, zmin=0.0, zmax = 0.6, samples=30, ytick distance = 0.5, xtick distance = 0.5, ztick distance = 0.5, colorbar horizontal, colorbar style={at={(0.5,1.1)},anchor=south,xticklabel pos=upper,title={$\ell_\mathrm{d}$ [mm]},title style={at={(0.5,1.4)}}}, point meta min = {0.1}, point meta max = {1}, colormap/cool,]
\addplot3 [surf,fill=lightgray,faceted color=lightgray] ({0.5*cos(deg(y))},{x},{0.5*sin(deg(y))});
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt] table [x = {x}, y = {y}, z={z}] {\table};
\end{axis}
\end{tikzpicture}

\caption{$\theta_\mathrm{d} = 80$ deg., $n = 2$}\label{inf_length}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.45\textwidth}
\centering

\pgfplotstableread{CHAP_6/shape_on_demand/var_order.dat}{\table}

\begin{tikzpicture}
\begin{axis}[view/az=50,
xlabel={x [mm]}, zlabel={y [mm]}, ylabel={-z [mm]}, y = {(0.3\textwidth,0.075\textwidth)}, x = {(0.3\textwidth,-0.075\textwidth)}, z = {(0.0\textwidth,0.315\textwidth)}, domain=0:1.2, y domain=0:pi*0.5, xmin=0, xmax=0.6, ymin=0, ymax=1.2, zmin=0.0, zmax = 0.6, samples=30, ytick distance = 0.5, xtick distance = 0.5, ztick distance = 0.5, colorbar horizontal, colorbar style={at={(0.5,1.1)},anchor=south,xticklabel pos=upper,title={$n$ [-]},title style={at={(0.5,1.4)}}}, point meta min = {0.5}, point meta max = {6}, colormap/cool,]
\addplot3 [surf,fill=lightgray,faceted color=lightgray] ({0.5*cos(deg(y))},{x},{0.5*sin(deg(y))});
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt] table [x = {x}, y = {y}, z={z}] {\table};
\end{axis}
\end{tikzpicture}

\caption{$\ell_\mathrm{d} = 1 \ \unit{mm}, \ \theta_\mathrm{d} = 80$ deg.}\label{inf_order}
\end{subfigure}
\hfill
\caption{Possible debonding shapes superimposed on the fiber-matrix interface obtained by varying (a) the debonding angle, (b) the debonding length and (c) the exponent of the power-law description.}\label{inf_param_shape}
\end{figure}

The influence of the debonding angle on the debonding front shape is shown in Figure~\ref{inf_angle}, while the exponent and length are arbitrarily fixed. The debonding half-angle $\theta_\mathrm{d}/2$ can be varied between 10 and 80 degrees to assess a wide range of debonding shapes. The influence of debonding length is shown in Figure~\ref{inf_length} and the power-law exponent in Figure~\ref{inf_order}, the other parameters being fixed. The debonding length varies between 0.05 and 1 \unit{mm} in Figure~\ref{inf_length}. The variation in exponent (Figure~\ref{inf_order}) induces a transition from a concave shape, e.g., for an exponent smaller than 1, similar to that provided by stress-based shapes, to a convex shape, e.g., for an exponent larger than 1, tending towards a squared debonding shape for a large exponent. Using the above geometrical parameters, a wide range of debonding shapes can be evaluated to determine the optimal energy-based debonding shapes.

A procedure is established to determine the energy-based debonding shapes by exploring all possible configurations for each exponent $n$ of the power-law description. The first step consists in identifying the couples of debonding angle and length leading to the same debonding surface for a fixed exponent. Figure~\ref{surf_n_1_5} shows the variation in debonding surface as a function of both debonding angle and length for another fixed $n$~=~1.5 exponent.

\begin{figure}[h!]
\begin{subfigure}[c]{0.49\textwidth}
\centering
\begin{tikzpicture}
\begin{axis}[view/az=300,
xlabel={$\ell_\mathrm{d}$ [mm]}, ylabel={$\theta_\mathrm{d}/2$ [deg.]}, zlabel={Debonding surface ($S$) [\unit{mm^2}]}, width = .8\textwidth, height = .8\textwidth, ytick distance = 10., xtick distance = 0.25, title = {$n=1.5$}]
\addplot3[surf,mesh/rows=40, mesh/cols=40] file[] {CHAP_6/n_1_5/surf_n_1.5.dat};
\end{axis}
\end{tikzpicture}

\caption{}\label{surf_n_1_5}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.45\textwidth}
\centering

\begin{tikzpicture}
\begin{axis}[view={0}{90},
width = 0.8\textwidth, height = 0.8\textwidth, xlabel={$\ell_\mathrm{d}$ [mm]}, ylabel={$\theta_\mathrm{d}/2$ [deg.]}, ytick distance = 10., xtick distance = 0.2, title = {$n=1.5$}, colorbar, colorbar style={title={$S$ [\unit{mm^2}]},title style={at={(0.5,0.97)}}}, point meta min = {0.008}, point meta max = {0.41},]
\addplot3[contour prepared={contour label style={/pgf/number format/fixed,/pgf/number format/precision=2,}},line width = 1 pt] file[] {CHAP_6/n_1_5/iso.dat};
\end{axis}
\end{tikzpicture}

\caption{}\label{iso_n_1_5}
\end{subfigure}
\caption{(a) Variation of the debonding surface as a function of the debonding angle and length. (b) Debonding surfaces isocontours highlighting the ranges of debonding angle and length couples leading to a similar debonding surface.}\label{fig:surf_iso_n_1_5}
\end{figure}

Increasing both the debonding angle and the length leads to increasing debonding surfaces. Consequently, several couples of angle and length, leading to the same debonding surface, can be determined using the debonding surface isocontours as shown in Figure~\ref{iso_n_1_5}. The minimum debonding surface considered is 0.001 \unit{mm^2}, with surfaces smaller than 0.01 \unit{mm^2} are not presented in the sequel for the sake of clarity. Ten ($\ell_\mathrm{d}$, $\theta_\mathrm{d}$) couples are selected for each debonding surface isocontour. The different debonding shapes obtained for two selected surfaces of 0.18~\unit{mm^2} and 0.3 \unit{mm^2} are shown in Figure~\ref{shape_debond_pairs}, noting that only five couples of angles and lengths are displayed for the sake of clarity.

\begin{figure}[h!]
\begin{subfigure}[c]{0.49\textwidth}
\centering

\pgfplotstableread{CHAP_6/shape/n=1.5_S=0.18_1.dat}{\tableun}
\pgfplotstableread{CHAP_6/shape/n=1.5_S=0.18_2.dat}{\tabledeux}
\pgfplotstableread{CHAP_6/shape/n=1.5_S=0.18_3.dat}{\tabletrois}
\pgfplotstableread{CHAP_6/shape/n=1.5_S=0.18_4.dat}{\tablequatre}
\pgfplotstableread{CHAP_6/shape/n=1.5_S=0.18_5.dat}{\tablecinq}

\begin{tikzpicture}
\begin{axis}[view/az=50,
xlabel={x [mm]}, zlabel={y [mm]}, ylabel={-z [mm]}, y = {(0.3*1.1\textwidth,0.075*1.1\textwidth)}, x = {(0.3*1.1\textwidth,-0.075*1.1\textwidth)}, z = {(0.0\textwidth,0.315*1.1\textwidth)}, domain=0:1., y domain=0:pi*0.5, xmin=0, xmax=0.6, ymin=0, ymax=1., zmin=0.0, zmax = 0.6, samples=30, ytick distance = 0.5, xtick distance = 0.5, ztick distance = 0.5, colorbar horizontal, colorbar style={title={$S$ [\unit{mm^2}]},title style={at={(0.5,0.97)}}}, colorbar to name={colorbarshape}, point meta min = {0.008}, point meta max = {0.41}, title = {$n=1.5$},]
\addplot3 [surf,fill=lightgray,faceted color=lightgray] ({0.5*cos(deg(y))},{x},{0.5*sin(deg(y))});
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt] table [x = {x}, y = {y}, z={z}] {\tableun};
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt,dashed] table [x = {x}, y = {y}, z={z}] {\tabledeux};
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt,densely dotted] table [x = {x}, y = {y}, z={z}] {\tabletrois};
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt,dashdotted] table [x = {x}, y = {y}, z={z}] {\tablequatre};
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt,densely dashdotted] table [x = {x}, y = {y}, z={z}] {\tablecinq};
\end{axis}
\end{tikzpicture}

\caption{$S=0.18$ \unit{mm^2}}\label{shape_debond_0.18}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.49\textwidth}
\centering

\pgfplotstableread{CHAP_6/shape/n=1.5_S=0.3_1.dat}{\tableun}
\pgfplotstableread{CHAP_6/shape/n=1.5_S=0.3_2.dat}{\tabledeux}
\pgfplotstableread{CHAP_6/shape/n=1.5_S=0.3_3.dat}{\tabletrois}
\pgfplotstableread{CHAP_6/shape/n=1.5_S=0.3_4.dat}{\tablequatre}
\pgfplotstableread{CHAP_6/shape/n=1.5_S=0.3_5.dat}{\tablecinq}

\begin{tikzpicture}
\begin{axis}[view/az=50,
xlabel={x [mm]}, zlabel={y [mm]}, ylabel={-z [mm]}, y = {(0.3*1.1\textwidth,0.075*1.1\textwidth)}, x = {(0.3*1.1\textwidth,-0.075*1.1\textwidth)}, z = {(0.0\textwidth,0.315*1.1\textwidth)}, domain=0:1., y domain=0:pi*0.5, xmin=0, xmax=0.6, ymin=0, ymax=1., zmin=0.0, zmax = 0.6, samples=30, ytick distance = 0.5, xtick distance = 0.5, ztick distance = 0.5, point meta min = {0.008}, point meta max = {0.41}, title = {$n=1.5$},]
\addplot3 [surf,fill=lightgray,faceted color=lightgray] ({0.5*cos(deg(y))},{x},{0.5*sin(deg(y))});
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt] table [x = {x}, y = {y}, z={z}] {\tableun};
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt,dashed] table [x = {x}, y = {y}, z={z}] {\tabledeux};
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt,densely dotted] table [x = {x}, y = {y}, z={z}] {\tabletrois};
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt,dashdotted] table [x = {x}, y = {y}, z={z}] {\tablequatre};
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt,densely dashdotted] table [x = {x}, y = {y}, z={z}] {\tablecinq};
\end{axis}
\end{tikzpicture}

\caption{$S=0.3$ \unit{mm^2}}\label{shape_debond_0.3}
\end{subfigure}
\caption{Debonding shape with varying couples of debonding angle and length obtained for a similar debonding surface of (a) 0.18 \unit{mm^2} and (b) 0.3 \unit{mm^2}.}\label{shape_debond_pairs}
\end{figure}

This illustrates the variation in debonding angle and length for a similar debonding surface. Noticeably, the debonding shapes corresponding to a surface of 0.18 \unit{mm^2} show large variations in length and angle (see Figure~\ref{shape_debond_0.18}). Since a maximum value is set on the length range to conform to experimental observations, smaller variations are encountered for a surface of 0.3 \unit{mm^2} (see Figure~\ref{shape_debond_0.3}).

The IERR corresponding to each couple can be evaluated by FE using a mesh that includes a debonding profile corresponding to the area delimited by a given debonding front. The debonding shape is projected onto the fiber-matrix interface so that the mesh topology matches the geometrical configuration. As a consequence, one FE mesh is created for each debonding front configuration.

The average critical ERR is calculated using the Hutchinson and Suo~\cite{hutchinson_mixed_1991} relationship, which depends on the mode mixity $\psi$. A condition is applied to the mode mixity to consider only shear when the interface undergoes compression, see Girard et al.~\cite{girard_numerical_2023} for more details. Mode mixity is evaluated locally which gives a local value of $G_\mathrm{c}$ that is averaged on the anticipated crack path, see~\eqref{GC_barre}:
\begin{equation}
\overline{G_\mathrm{c}}(\Gamma) = \frac{1}{S} \int_{}^{\Gamma} G_{\mathrm{IC}}\bracks[\Big]{1 + \tan^2\bracks[\big]{(1 - \lambda) \psi(\vec{x})}} \,\dd x \,\dd y \,\dd z.\label{GC_barre}
\end{equation}

The parameter $\lambda$ is related to a mode II to mode I critical ERR ratio, see~\cite{doitrand_review_2024}.

The IERR to critical ERR ratio variation is shown in Figure~\ref{fig:Ginc_Gc_all_n_1_5} as a function of $\theta_\mathrm{d}$ for each selected ($\ell_\mathrm{d}$, $\theta_\mathrm{d}$) couples and varying debonding surfaces, which values are differentiated using the same color code as Figures~\ref{fig:surf_iso_n_1_5} and~\ref{shape_debond_pairs}.

\begin{figure}[h!]
     \begin{subfigure}[c]{0.49\textwidth}
         \centering
        
\pgfplotstableread{CHAP_6/n_1_5/Ginc_Gc_max.dat}{\table}        

\begin{tikzpicture}
\begin{axis}[view={0}{90},
width = 0.95\textwidth,
height = 0.855\textwidth,
ylabel={$G_\mathrm{inc}/\overline{G_\mathrm{c}}$ [-]},
xlabel={$\theta_\mathrm{d}/2$ [deg.]},
title = {$n=1.5$, $\lambda = 0.2$},
colorbar horizontal,
colorbar style={title={$S$ [mm$^2$]},title style={at={(0.5,0.97)}}},
colorbar to name={colorbarGinc},
point meta min = {0.008},
point meta max = {0.41},
xmin = 25,
xmax = 90,
ymin = 400,
ymax = 850,
xtick distance = 10,
ytick distance = 100,
]
\addplot3[contour prepared={labels=false},line width = 1 pt,mark size = 1pt] file[] {CHAP_6/n_1_5/Ginc_Gc_S.dat};
\addplot[scatter,only marks,point meta=\thisrow{c}] table [x = {x}, y = {y}] {\table};
\end{axis}
\end{tikzpicture}
         \caption{}
         \label{fig:Ginc_Gc_all_n_1_5}
     \end{subfigure}
     \hfill
     \begin{subfigure}[c]{0.49\textwidth}
         \centering
         
\pgfplotstableread{CHAP_6/n_1_5/Ginc_Gc.dat}{\table}

\begin{tikzpicture}
\begin{axis}[
width = 0.95\textwidth,
height = 0.855\textwidth,
ylabel={$\max \left( G_\mathrm{inc}/\overline{G_\mathrm{c}} \right)$ [-]},
xlabel={$S$ [mm$^2$]},
title = {$n=1.5$, $\lambda = 0.2$},
%colorbar,
%colorbar style={title={$S$ [mm$^2$]},title style={at={(0.5,0.97)}}},
point meta min = {0.008},
point meta max = {0.41},
ymin = 400,
ymax = 850,
tick label style={/pgf/number format/fixed}
]
\addplot[scatter,point meta=\thisrow{c}] table [x = {x}, y = {y}] {\table};
\end{axis}
\end{tikzpicture}
                
        \caption{}
     	\label{fig:Ginc_Gc_max_n_1_5}
     \end{subfigure}

\includegraphics[width=.35\textwidth]{./figures/Ginc_Gc_n_1_5}
\caption{(a) IERR to critical ERR variation obtained for different debonding surfaces as a function of the debonding angle. The maximum is indicated by a marker for each surface. (b) Maximum IERR to critical ERR obtained for each debonding surface.}
        \label{fig:Ginc_Gc_n_1_5}
\end{figure}

Overall, for a given debonding surface, there exists a ($\ell_\mathrm{d}$, $\theta_\mathrm{d}$) couple that maximizes the energy criterion, i.e., the IERR to critical ERR ratio, indicated by a round marker. Consequently, these couples lead to the geometrical configuration that is the most favorable from an energy point of view. The maximum of IERR to critical ERR ratio obtained from each surface is shown in Figure~\ref{fig:Ginc_Gc_max_n_1_5}, with a similar marker color, as a function of the corresponding debonding surface. It yields the maximum value of IERR that can be attained as a function of the crack surface for given exponent $n = 1.5$ and $\lambda = 0.2$ parameter.

The debonding angle and length couples that maximize the IERR to critical ERR ratio are shown in Figure~\ref{fig:path_n_1_5}. These couples therefore correspond to the most energetically favorable geometrical configuration. Since only the debonding initiation is studied here, the angle and length of debonding do not necessarily have to increase monotonically. For instance, a decrease in the debonding angle is observed as the debonding surface increases, see Figure~\ref{fig:path_n_1_5} for surfaces close to 0.4 \unit{mm^2} ($\ell_\mathrm{d}$ close to 0.8 \unit{mm}).

The configuration corresponding to the favorable debonding shape is superimposed onto the fiber-matrix interface in Figure~\ref{fig:iso_W_n_1_5}. The debonding remains localized near the free surface fiber pole for small surface areas. For a larger debonding area, there is a more pronounced increase in debonding length with respect to the angle.

\begin{figure}[h!]
\begin{subfigure}[c]{0.49\textwidth}
\centering

\pgfplotstableread{CHAP_6/n_1_5/Path_lamb=0.2_n=1.5_with_surf.dat}{\table}

\begin{tikzpicture}
\begin{axis}[
width = 0.8\textwidth, height = 0.72\textwidth, ylabel={$\ell_\mathrm{d}$ [mm]}, xlabel={$\theta_\mathrm{d}/2$ [deg.]}, xtick distance = 10., ytick distance = 0.2, title = {$n=1.5$, $\lambda = 0.2$}, colorbar horizontal, colorbar style={title={$S$ [\unit{mm^2}]},title style={at={(0.5,0.97)}}}, colorbar to name={colorbariso}, point meta min = {0.008}, point meta max = {0.41},]
\addplot[scatter,only marks,point meta=\thisrow{c}] table [x = {x}, y = {y}] {\table};
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:path_n_1_5}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.49\textwidth}
\centering

\pgfplotstableread{CHAP_6/n_1_5/iso_lamb=0.2_n=1.5.dat}{\table}

\begin{tikzpicture}
\begin{axis}[view/az=50,
xlabel={x [mm]}, zlabel={y [mm]}, ylabel={-z [mm]}, y = {(1.2*0.3\textwidth,1.2*0.075\textwidth)}, x = {(1.2*0.3\textwidth,-0.075*1.2\textwidth)}, z = {(0.0\textwidth,1.2*0.315\textwidth)}, domain=0:1, y domain=0:pi*0.5, xmin=0, xmax=0.6, ymin=0, ymax=1, zmin=0.0, zmax = 0.6, samples=30, ytick distance = 0.5, xtick distance = 0.5, ztick distance = 0.5, point meta min = {0.008}, point meta max = {0.41}, title = {$n=1.5$, $\lambda = 0.2$},]
\addplot3 [surf,fill=lightgray,faceted color=lightgray] ({0.5*cos(deg(y))},{x},{0.5*sin(deg(y))});
\addplot3[contour prepared={labels=false}, point meta=\thisrow{level}, line width = 1pt] table [x = {x}, y = {y}, z={z}] {\table};
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:iso_W_n_1_5}
\end{subfigure}

\includegraphics[width=.3\textwidth]{./figures/path_iso_n_1_5}
\caption{(a) Debonding angle and length couples corresponding to the maximum in IERR to critical ERR ratio in Figure~\ref{fig:Ginc_Gc_all_n_1_5}, with the colors corresponding to the debonding surfaces. (b) Optimal potential debonding shapes determined from an energy point of view with an exponent of 1.5, superimposed onto the fiber-matrix interface (gray surface).}\label{fig:path_iso_n_1_5}
\end{figure}

The debonding shapes shown in Figure~\ref{fig:iso_W_n_1_5} are those that maximize the energy condition. However, they are obtained for a fixed power-law exponent $n$. Thus, the optimum energy conditions obtained for different power-law exponents can be compared to determine the configurations which maximize the energy conditions. Figure~\ref{fig:Ginc_n} shows the variation in the IERR to critical ERR ratio obtained for different exponents. It is found that a $n=1.5$ exponent yields the most favorable debonding shapes from an energy point of view whatever the debonding surface and $\lambda$ values explored.

\begin{figure}[h!]
\begin{tikzpicture}
\begin{axis}[
width = 0.45\textwidth, height = 0.36\textwidth, ylabel={$G_\mathrm{inc}/\overline{G_\mathrm{c}}$ [-]}, xlabel={$S$ [\unit{mm^2}]}, title = {$\lambda = 0.2$}, legend style={at={(1.05,1.)},anchor=north west}, cycle list name=black white, tick label style={/pgf/number format/fixed}]
\addplot+[] file[] {CHAP_6/Ginc/Ginc_lamb=0.2_n=0.6.dat};
\addplot+[] file[] {CHAP_6/Ginc/Ginc_lamb=0.2_n=0.8.dat};
\addplot+[] file[] {CHAP_6/Ginc/Ginc_lamb=0.2_n=1.0.dat};
\addplot+[] file[] {CHAP_6/Ginc/Ginc_lamb=0.2_n=1.5.dat};
\addplot+[] file[] {CHAP_6/Ginc/Ginc_lamb=0.2_n=2.0.dat};
\addplot+[] file[] {CHAP_6/Ginc/Ginc_lamb=0.2_n=4.0.dat};
\legend{$n=0.6$,$n=0.8$,$n=1.0$,$n=1.5$,$n=2.0$,$n=4.0$}
\end{axis}
\end{tikzpicture}
\caption{Variation of the optimal IERR to critical ERR ratio as a function of the debonding surface for different power-law exponents, emphasizing the optimal exponent of 1.5 for the defined $\lambda = 0.2$.}\label{fig:Ginc_n}
\end{figure}

The fiber-matrix interface undergoes both opening and shear modes simultaneously. Consequently, $\lambda$ can influence the loading required to fulfill the energy condition. Figure~\ref{fig:Ginc_lamb} shows the influence of $\lambda$ on the magnitude of the IERR to critical ERR ratio.

\begin{figure}[h!]
\begin{subfigure}[c]{0.49\textwidth}
\centering

\begin{tikzpicture}
\begin{axis}[
width = 0.9\textwidth, height = 0.81\textwidth, ylabel={$G_\mathrm{inc}/\overline{G_\mathrm{c}}$ [-]}, xlabel={$S$ [\unit{mm^2}]}, title = {$n = 1.5$}, cycle list name=black white, legend columns=-1, legend entries={$\lambda=0.01$,$\lambda=0.1$,$\lambda=0.2$,$\lambda=0.4$,$\lambda=0.6$,$\lambda=0.8$}, legend to name=Ginc_lamb, ytick distance = 100, tick label style={/pgf/number format/fixed}]
\addplot+[] file[] {CHAP_6/lamb/Ginc_lamb=0.01_n=1.5.dat};
\addplot+[] file[] {CHAP_6/lamb/Ginc_lamb=0.1_n=1.5.dat};
\addplot+[] file[] {CHAP_6/lamb/Ginc_lamb=0.2_n=1.5.dat};
\addplot+[] file[] {CHAP_6/lamb/Ginc_lamb=0.4_n=1.5.dat};
\addplot+[] file[] {CHAP_6/lamb/Ginc_lamb=0.6_n=1.5.dat};
\addplot+[] file[] {CHAP_6/lamb/Ginc_lamb=0.8_n=1.5.dat};
\node[] at (axis description cs:0.52,0.4) {$\lambda_\searrow$};
\draw[->] (axis description cs:0.25,0.85) -- (axis description cs:0.45,0.3);
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:Ginc_lamb}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.49\textwidth}
\centering

\begin{tikzpicture}
\begin{axis}[
width = 0.9\textwidth, height = 0.81\textwidth, ylabel={$\sigma^\infty_\mathrm{req}$ [-]}, xlabel={$S$ [\unit{mm^2}]}, title = {$n = 1.5$}, cycle list name=black white, ytick distance = 0.5, tick label style={/pgf/number format/fixed}]
\addplot+[] file[] {CHAP_6/lamb/Ginc_Gc_lamb=0.01_n=1.5.dat};
\addplot+[] file[] {CHAP_6/lamb/Ginc_Gc_lamb=0.1_n=1.5.dat};
\addplot+[] file[] {CHAP_6/lamb/Ginc_Gc_lamb=0.2_n=1.5.dat};
\addplot+[] file[] {CHAP_6/lamb/Ginc_Gc_lamb=0.4_n=1.5.dat};
\addplot+[] file[] {CHAP_6/lamb/Ginc_Gc_lamb=0.6_n=1.5.dat};
\addplot+[] file[] {CHAP_6/lamb/Ginc_Gc_lamb=0.8_n=1.5.dat};
\node[] at (axis description cs:0.57,0.45) {$\lambda_\searrow$};
\draw[->] (axis description cs:0.2,0.1) -- (axis description cs:0.5,0.55);
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:Ginc_Gc_lamb}
\end{subfigure}
\ref{Ginc_lamb}
\caption{Influence of $\lambda$ on (a) the IERR to critical ERR ratio and (b) the loading required to fulfill the energy criterion ($\sigma^\infty_\mathrm{req}$) for a fixed exponent of $n=1.5$.}\label{fig:inf_lamb_energy_condi}
\end{figure}

The opening critical ERR $G_\mathrm{IC}$ is fixed while the shear critical ERR $G_\mathrm{IIC}$ increases, leading to a decrease in $\lambda$. The critical ERR in mode II appears to influence the profile of the energy criterion over the full range of the debonding surface, suggesting that small variations in mode mixity are encountered on the different debonding shapes. Overall, increasing $G_\mathrm{IIC}$ leads to a decrease in the level of the energy criterion. Consequently, a larger loading is required to fulfill the energy criterion as $G_\mathrm{IIC}$ increases, see Figure~\ref{fig:Ginc_Gc_lamb}. However, the difference in required loading induced by a variation in $\lambda$ is small, less than 0.5 \unit{MPa} for a fixed debonding surface. Consequently, the shear component appears to be negligible on the energy-based shapes determined.

Although shear has little influence on the loading required to fulfill the energy condition, it can influence the ($\ell_\mathrm{d}$, $\theta_\mathrm{d}$) couples that describe the energy-based shapes. Figure~\ref{fig:path_inf_lamb} shows the ($\ell_\mathrm{d}$, $\theta_\mathrm{d}$) couples corresponding to energy-based shapes for two extreme $\lambda$ of 0.01 and 0.8.

\begin{figure}[h!]
\begin{subfigure}[c]{0.49\textwidth}
\centering

\pgfplotstableread{CHAP_6/lamb/Path_lamb=0.01_n=1.5_with_surf.dat}{\tableone}
\pgfplotstableread{CHAP_6/lamb/Path_lamb=0.8_n=1.5_with_surf.dat}{\tabletwo}

\begin{tikzpicture}
\begin{axis}[
width = 0.9\textwidth, height = 0.81\textwidth, ylabel={$\ell_\mathrm{d}$ [mm]}, xlabel={$\theta_\mathrm{d}/2$ [deg.]}, xtick distance = 10., ytick distance = 0.2, title = {$n=1.5$}, colorbar horizontal, colorbar style={at={(0.5,-0.35)},anchor=north,xticklabel pos=lower,title={$S$ [\unit{mm^2}]},title style={at={(0.5,-1.1)},anchor=north}}, point meta min = {0.008}, point meta max = {0.41}, legend pos = north west,]
\addplot[scatter,only marks,mark=triangle*,point meta=\thisrow{c}] table [x = {x}, y = {y}] {\tableone};
\addplot[scatter,only marks,mark=square*,point meta=\thisrow{c}] table [x = {x}, y = {y}] {\tabletwo};
\legend{$\lambda=0.01$,$\lambda=0.8$}
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:path_inf_lamb}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.49\textwidth}
\centering

\pgfplotstableread{CHAP_6/lamb/mesh_3D_mode_mixity.dat}{\table}

\begin{tikzpicture}
\begin{axis}[view/az=50,
xlabel={x [mm]}, zlabel={y [mm]}, ylabel={-z [mm]}, y = {(0.2*1.8\textwidth,0.05*1.8\textwidth)}, x = {(0.2*1.8\textwidth,-0.05*1.8\textwidth)}, z = {(0.0\textwidth,0.21*1.8\textwidth)}, ymax = 1, xtick distance = 0.5, ztick distance = 0.5, colormap/bluered, colorbar horizontal, colorbar style={at={(0.5,1.2)},anchor=south,xticklabel pos=upper,title={$\psi$ [rad.]},title style={at={(0.5,1.4)}}}, point meta min = {0}, point meta max = {1.57}, legend columns = 2, legend style={at={(0.5,-0.7)},anchor=north},]
\addplot3[patch, patch type=rectangle, shader=faceted interp,forget plot] table [x = {x}, y = {y}, z = {z}, point meta=\thisrow{c}] {\table};
\addplot3[line width = 1pt] file[] {CHAP_6/lamb/iso_select/iso_lamb=0.8_n=1.5_S=0.05.dat};
\addplot3[line width = 1pt,dashed,gray] file[] {CHAP_6/lamb/iso_select/iso_lamb=0.01_n=1.5_S=0.05.dat};
\addplot3[line width = 1pt] file[] {CHAP_6/lamb/iso_select/iso_lamb=0.8_n=1.5_S=0.18.dat};
\addplot3[line width = 1pt,dashed,gray] file[] {CHAP_6/lamb/iso_select/iso_lamb=0.01_n=1.5_S=0.18.dat};
\addplot3[line width = 1pt] file[] {CHAP_6/lamb/iso_select/iso_lamb=0.8_n=1.5_S=0.3.dat};
\addplot3[line width = 1pt,dashed,gray] file[] {CHAP_6/lamb/iso_select/iso_lamb=0.01_n=1.5_S=0.3.dat};
\legend{$\lambda = 0.8$~~~~,$\lambda = 0.01$}
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:path_lamb_mixity}
\end{subfigure}
\caption{(a) Debonding geometrical parameters ($\ell_\mathrm{d}$, $\theta_\mathrm{d}$) corresponding to the energy-based debonding shapes for two extremes values of $\lambda$. (b) Energy-based shapes obtained for three surfaces and two $\lambda$ values highlighting the influence of the mode mixity (whose variation is superimposed on the fiber-matrix interface) on the debonding shapes.}\label{fig:inf_mixity_shape}
\end{figure}

For a relatively small debonding area ($\ell_\mathrm{d} < 0.3$ \unit{mm} and $\theta_\mathrm{d}/2 < 45$ deg.), $\lambda$ has no influence on the energy-based debonding shape, which mainly undergoes opening mode. However, for a larger debonding surface, smaller $\lambda$ results in both decreasing $\theta_\mathrm{d}$ and increasing $\ell_\mathrm{d}$. This difference can be explained by the stress shear component that is higher close to the fiber equator, i.e., at large debonding angles. Consequently, the increase in the critical ERR in mode II makes debonding initiation for a large angle less favorable than for a small $G_\mathrm{IIC}$. This is illustrated in Figure~\ref{fig:path_lamb_mixity}, where the optimum debonding shapes obtained for three different debonding surfaces and the two extreme values of $\lambda$ are superimposed with the field of variation in mode mixity at the fiber-matrix interface. Overall, mode mixity is similar to that observed with a front 2D modeling (see Girard et al.~\cite{girard_numerical_2023}) for a sufficiently large distance from the free surface and it is affected by the singularity at the free surface. As with debonding shapes based on stress isocontours, the fiber diameter to sample thickness ratio should not affect energy-based shapes for a fiber that is sufficiently long in relation to its diameter.

\subsection{Energy-based vs stress-based: CC solution}

Once the potential initiation debonding shapes based on stress isocontours or energy are determined, the CC can be implemented. The loading required to fulfill one or other of the criteria ($\sigma^\infty_\sigma$, $\sigma^\infty_\mathrm{G}$) is compared in Figure~\ref{fig:CC_compare_sol} on the basis of the stress isocontours-based and energy-based debonding shapes for two brittleness numbers:
\begin{equation}
\gamma = \frac{1}{\sigma_\mathrm{c}} \sqrt{\frac{G_{\mathrm{IC}}E^*}{r}},\label{eq:gamma}
\end{equation}
where $E^*$ denotes the harmonic mean of the effective elastic moduli of the two constituents. The $\lambda$ and $\mu$ parameters are selected on the basis of those identified in~\cite{girard_numerical_2023}. The minimum loading for which both criteria are fulfilled is the initiation loading and the corresponding debonding surface is the initiation crack surface, which are indicated by the round marker in Figure~\ref{fig:CC_compare_sol}.

\begin{figure}[h!]
\begin{subfigure}[c]{0.49\textwidth}
\centering

\begin{tikzpicture}
\begin{axis}[width = 0.9\textwidth,
height = 0.81\textwidth, ylabel={$\sigma^\infty_\mathrm{req}$ [\unit{MPa}]}, xlabel={$S$ [\unit{mm^2}]}, title = {$\lambda = 0.13$, $\mu = 3$, $\gamma=0.5$}, legend columns=2, legend entries={$\sigma^\infty_\sigma$ energy-based~~~~,$\sigma^\infty_\mathrm{G}$ energy-based,$\sigma^\infty_\sigma$ stress-based~~~~,$\sigma^\infty_\mathrm{G}$ stress-based}, legend to name=CC_sol_leg, xmax = 0.1, ytick distance = 5., xtick distance = 0.025, tick label style={/pgf/number format/fixed}]
\addplot[blue] file[] {CHAP_6/CC/ener_iso_CC_stress_gamma=0.5.dat};
\addplot[blue,dashed] file[] {CHAP_6/CC/ener_iso_CC_energy_gamma=0.5.dat};
\addplot[red] file[] {CHAP_6/CC/stress_iso_CC_stress_gamma=0.5.dat};
\addplot[red,dashed] file[] {CHAP_6/CC/stress_iso_CC_energy_gamma=0.5.dat};
\addplot[only marks,red,mark size=1.5pt] coordinates {(0.0224,6.906)};
\addplot[only marks,blue,mark size=1.5pt] coordinates {(0.0039,6.13)};
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:gamma_0_5}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.49\textwidth}
\centering

\begin{tikzpicture}
\begin{axis}[width = 0.9\textwidth,
height = 0.81\textwidth, ylabel={$\sigma^\infty_\mathrm{req}$ [\unit{MPa}]}, xlabel={$S$ [\unit{mm^2}]}, title = {$\lambda = 0.13$, $\mu = 3$, $\gamma=3$}, xtick distance = 0.5, ytick distance = 5.,]
\addplot[blue] file[] {CHAP_6/CC/ener_iso_CC_stress_gamma=3.0.dat};
\addplot[blue,dashed] file[] {CHAP_6/CC/ener_iso_CC_energy_gamma=3.0.dat};
\addplot[red] file[] {CHAP_6/CC/stress_iso_CC_stress_gamma=3.0.dat};
\addplot[red,dashed] file[] {CHAP_6/CC/stress_iso_CC_energy_gamma=3.0.dat};
\addplot[only marks,red,mark size=1.5pt] coordinates {(1.05,5.45)};
\addplot[only marks,blue,mark size=1.5pt] coordinates {(0.107,4.48)};
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:gamma_3}
\end{subfigure}
\vspace{.5\baselineskip}
\ref{CC_sol_leg}

\vskip .5\baselineskip

\begin{subfigure}[c]{0.49\textwidth}
\centering

\begin{tikzpicture}
\begin{axis}[view/az=50,
xlabel={x [mm]}, zlabel={y [mm]}, ylabel={-z [mm]}, y = {(0.3*1.43\textwidth,0.075*1.43\textwidth)}, x = {(0.3*1.43\textwidth,-0.075*1.43\textwidth)}, z = {(0.0\textwidth,0.315*1.43\textwidth)}, domain=0:0.5, y domain=0:pi*0.5, xmin=0, xmax=0.6, ymin=0, ymax=0.5, zmin=0.0, zmax = 0.6, samples=6, ytick distance = 0.5, xtick distance = 0.5, ztick distance = 0.5, legend columns=2, legend entries={Stress-based~~~~,Energy-based}, legend to name=CC_surf_leg,]
\addplot3 [surf,fill=lightgray,faceted color=lightgray,forget plot] ({0.5*cos(deg(y))},{x},{0.5*sin(deg(y))});
\addplot3[fill = red,opacity=0.5, mark=none, draw=none,area legend] file[] {CHAP_6/Sc/stress_surf_gamma_0.5.dat};
\addplot3[fill = blue,opacity=0.5, mark=none, draw=none,area legend] file[] {CHAP_6/Sc/ener_surf_gamma_0.5.dat};
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:surf_0_5}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.49\textwidth}
\hspace{-.5cm}
\begin{tikzpicture}
\begin{axis}[view/az=50,
xlabel={x [mm]}, zlabel={y [mm]}, ylabel={-z [mm]}, y = {(0.3\textwidth,0.075\textwidth)}, x = {(0.3\textwidth,-0.075\textwidth)}, z = {(0.0\textwidth,0.315\textwidth)}, domain=0:2.375, y domain=0:pi*0.5, xmin=0, xmax=0.6, ymin=0, ymax=2.375, zmin=0.0, zmax = 0.6, samples=30, ytick distance = 0.5, xtick distance = 0.5, ztick distance = 0.5,]
\addplot3 [surf,fill=lightgray,faceted color=lightgray] ({0.5*cos(deg(y))},{x},{0.5*sin(deg(y))});
\addplot3[fill = red,opacity=0.5, mark=none, draw=none] file[] {CHAP_6/Sc/stress_surf_gamma_3.0.dat};
\addplot3[fill = blue,opacity=0.5, mark=none, draw=none] file[] {CHAP_6/Sc/ener_surf_gamma_3.0.dat};
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:surf_3}
\end{subfigure}
\hfill

\vspace{3mm}~\ref{CC_surf_leg}
\caption{Required loading to fulfill both the stress (solid lines) and energy (dashed lines) criteria obtained using the stress-based and energy-based shapes. Two interface brittleness numbers are investigated, (a) a relatively small one \smash{$\gamma = 0.5$} and (b) a relatively large one \smash{$\gamma = 3$}. Comparison of the debonding initiation surfaces ($S_\mathrm{c}$) obtained in (a) and (b), respectively (c) and (d).}\label{fig:CC_compare_sol}
\end{figure}

For both stress-based and energy-based debonding shapes, the stress condition increases monotonically since the equivalent stress decreases with the debonding surfaces. The energy condition obtained with the energy-based debonding shapes monotonically decreases, unlike the one obtained using the stress-based shapes which provides local minima. Overall, the stress-based shapes result in the most favorable solution in terms of stress conditions, whatever the debonding surface considered, consistently with their definition. On the contrary, the energy condition provided by energy-based shapes is smaller than that of stress-based shapes, whatever the debonding surface. For a relatively small brittleness number (Figure~\ref{fig:gamma_0_5}), the energy-based shapes provide the most favorable configuration compared to that provided by the stress-based shapes. This differs from the results obtained in 2D~\cite{girard_influence_2024}, where the stress-based shapes are the optimal debonding shapes for a small brittleness numbers, i.e., smaller than 1. In 3D, the energy condition, which is significantly larger for the stress-based shapes, appears to control debonding initiation even for a small brittleness number.
The stress-based and energy-based shapes provide debonding initiation surfaces with relative difference of 140\%.
For a relatively large brittleness number (Figure~\ref{fig:gamma_3}), the energy-based debonding shapes also lead to the most favorable configuration for predicting debonding initiation. The stress-based shapes require a larger loading to fulfill the energy condition. A similar result was obtained in 2D for large brittleness numbers~\cite{girard_influence_2024}. The debonding initiation configuration (surface $S_\mathrm{c}$) obtained for the stress-based and energy-based debonding shapes and the two brittleness numbers are compared in Figure~\ref{fig:CC_compare_sol}. For a small brittleness number, the stress-based surface remains close to the free surface with a debonding angle spanning over the whole fiber circumference. The energy-based surface has a smaller angle and length, thus a smaller surface than stress-based surface. The difference between the two debonding initiation surfaces becomes more pronounced as the brittleness number increases, as shown in Figure~\ref{fig:surf_3}. The stress-based surface provides a tunnel between the two free surface debondings. In contrast, the energy-based surface shape is similar to that obtained with a small brittleness number, with increased surface area.

Similarly, the debonding initiation surfaces obtained using stress-based and energy-based shapes are compared for a range of interface brittleness numbers, see Figure~\ref{fig:Sc_gamma}. The brittleness number range is selected to ensure that the CC solution is provided by the intersection of both stress and energy conditions.

\begin{figure}[h!]
\begin{subfigure}[c]{0.49\textwidth}
\centering

\begin{tikzpicture}
\begin{semilogxaxis}[
xmode=log, ymode=normal, width = 0.9\textwidth, height = 0.81\textwidth, ylabel={$S_\mathrm{c}$ [\unit{mm^2}]}, xlabel={$\gamma$ [-]}, title = {$\lambda = 0.13$, $\mu = 3$}, log ticks with fixed point, legend entries={Energy-based~~~~,Stress-based}, legend to name=CC_sol_gamma, legend columns = -1,]
\addplot[blue] file[] {CHAP_6/Sc/ener_iso_Sc.dat};
\addplot[red] file[] {CHAP_6/Sc/stress_iso_Sc.dat};
\node[rotate=0] at (axis description cs:0.78,0.35) {\rotatebox{0}{Tunneling}};
\draw[->] (axis description cs:0.6,0.35) -- (axis description cs:0.45,0.35);
\end{semilogxaxis}
\end{tikzpicture}

\caption{}\label{fig:Sc_gamma}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.49\textwidth}
\centering

\begin{tikzpicture}
\begin{semilogxaxis}[
xmode=log, ymode=normal, width = 0.8\textwidth, height = 0.81\textwidth, ylabel={$\sigma^\infty_\mathrm{req}$ [\unit{MPa}]}, xlabel={$\gamma$ [-]}, title = {$\lambda = 0.13$, $\mu = 3$}, log ticks with fixed point,]
\addplot[blue] file[] {CHAP_6/CC/ener_iso_CC_gamma.dat};
\addplot[red] file[] {CHAP_6/CC/stress_iso_CC_gamma.dat};
\draw[<->] (axis description cs:0.43,0.23) -- (axis description cs:0.43,0.51);
\draw[-] (axis description cs:0.44,0.37) -- (axis description cs:0.6,0.6);
\node[rotate=0] at (axis description cs:0.6,0.65) {\rotatebox{0}{ $\Delta \sigma^\infty \simeq$ 30\%}};
\end{semilogxaxis}
\end{tikzpicture}

\caption{}\label{fig:gamma_inf}
\end{subfigure}
\ref{CC_sol_gamma}
\caption{(a) Debonding initiation surfaces and (b) corresponding loading predicted by the stress-based and energy-based shapes.}\label{fig:gamma_inf_CC_sol}
\end{figure}

For interface brittleness numbers smaller than 2, the initiation surfaces predicted by the two approaches are in the same order of magnitude. However, for large brittleness number, the initiation surfaces predicted by the stress isocontours approach is larger than those provided by the energy-based debonding shapes. This phenomenon is caused by the tunneling effect of the stress-based shapes, which induces a sharp increase in the debonding surface for solutions triggered by the energy condition. The remote stress required to fulfill the CC is assessed in Figure~\ref{fig:gamma_inf}, for the same range of brittleness number as Figure~\ref{fig:Sc_gamma}, so that the relevancy of the stress isocontours approach in predicting the optimal CC solution can be evaluated. The stress-based shapes yield larger initiation remote stress than the energy-based ones whatever the brittleness number over the range 0.3--8.
In fact, it is observed that the energy-based debonding shapes are close to the ones that allow fulfilling the CC for a minimum imposed loading, as they yield differences smaller than 0.1\% in the initiation remote loading compared to these optimum solutions.

For relatively small interface brittleness numbers ($\gamma < 1$), limited difference of less than 15\% is made by the stress-based shapes.
For intermediate brittleness numbers ($1 < \gamma < 6$), the difference reaches a peak to 30\% caused by energy driven solutions.
The difference decreases to 15\% for larger brittleness numbers ($\gamma > 6$).
The differences in initiation loading are larger than in 2D, where it remains smaller than 5\% whether stress-based or energy-based crack shapes are used~\cite{girard_influence_2024}.

\subsection{Comparison with experimental observations}\label{sec:compare_exp}

Initiation debonding shapes are now compared to the ones observed experimentally just after initiation. The debonding side views are shown in Figure~\ref{fig:SFII_2_without} and Figure~\ref{fig:SFII_4_without} just after initiation for sample 1 and 3, respectively. Note that these two samples represent two independent cases subjected to the same experimental conditions. The debonding zone in each case is highlighted with white pixels, similar to Uddin et al.~\cite{uddin_simultaneous_2024}.

\begin{figure}[h!]
\begin{subfigure}[c]{0.4\textwidth}
\centering

\begin{tikzpicture}
\begin{axis}[enlargelimits=false,
width = \textwidth, ticks=none, xlabel = {z-direction}, ylabel = {y-direction}, title = {Sample 1}]
\addplot graphics [xmin=0,xmax=202,ymin=0,ymax=150] {CHAP_6/exp/SFII-2.png};
\node[rotate=0,white] at (axis description cs:0.5,0.5) {\rotatebox{0}{Fiber}};
\node[rotate=0,white] at (axis description cs:0.5,0.83) {\rotatebox{0}{Matrix}};
\node[rotate=0,white] at (axis description cs:0.35,0.07) {\rotatebox{0}{0.5 mm}};
\node[rotate=0,white] at (axis description cs:0.75,0.15) {\rotatebox{0}{Free surface}};
\draw[->,white] (axis description cs:0.75,0.2) -- (axis description cs:0.72,0.4);
\addplot[white,line width = 3pt] coordinates {(5,8) (32,8) };
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:SFII_2_without}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.59\textwidth}
\centering

\begin{tikzpicture}
\begin{axis}[enlargelimits=false,
width = 2.*0.24\textwidth, height = 2.*0.36\textwidth, ticks=none, xlabel = {z-direction}, ylabel = {y-direction}, title = {Stress}, legend style={at={(0.5,-0.14)},anchor=north},]
\addplot[forget plot] graphics [xmin=0,xmax=24,ymin=0,ymax=36] {CHAP_6/exp/SFII-2_zoom.png};
\addplot[red, line width = 1pt] file[] {CHAP_6/exp/SFII-2_stress_gamma_1.dat};
\addplot[red,dotted, line width = 1pt] file[] {CHAP_6/exp/SFII-2_stress_gamma_0.56.dat};
\addplot[red,dashed, line width = 1pt] file[] {CHAP_6/exp/SFII-2_stress_gamma_0.3.dat};
\addplot[fill=red, opacity=0.3,area legend] file[] {CHAP_6/exp/SFII-2_stress_gamma_1_area.dat};
\addplot[white,line width = 3pt] coordinates {(2,34) (8,34)};
\node[rotate=0,white] at (axis description cs:0.68,0.925) {\rotatebox{0}{10 \unit{\micro\metre}}};
\legend{$\gamma = 1$,$\gamma = 0.56$,$\gamma = 0.3$}
\end{axis}
\end{tikzpicture}
\begin{tikzpicture}
\begin{axis}[enlargelimits=false,
width = 2.*0.24\textwidth, height = 2.*0.36\textwidth, ticks=none, xlabel = {z-direction}, ylabel = {y-direction}, title = {Energy}, legend style={at={(0.5,-0.14)},anchor=north},]
\addplot[forget plot] graphics [xmin=0,xmax=24,ymin=0,ymax=36] {CHAP_6/exp/SFII-2_zoom.png};
\addplot[blue, line width = 1pt] file[] {CHAP_6/exp/SFII-2_ener_gamma_1.dat};
\addplot[blue,dotted, line width = 1pt] file[] {CHAP_6/exp/SFII-2_ener_gamma_0.6.dat};
\addplot[blue,dashed,line width = 1pt] file[] {CHAP_6/exp/SFII-2_ener_gamma_0.4.dat};
\addplot[fill=blue, opacity=0.3,area legend] file[] {CHAP_6/exp/SFII-2_ener_gamma_1_area.dat};
\addplot[white,line width = 3pt] coordinates {(2,34) (8,34)};
\node[rotate=0,white] at (axis description cs:0.68,0.925) {\rotatebox{0}{10 \unit{\micro\metre}}};
\legend{$\gamma = 1$,$\gamma = 0.6$,$\gamma = 0.4$}
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:SFII_2_with}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.4\textwidth}
\centering

\begin{tikzpicture}
\begin{axis}[enlargelimits=false,
width = \textwidth, ticks=none, xlabel = {z-direction}, ylabel = {y-direction}, title = {Sample 3},]
\addplot graphics [xmin=0,xmax=259,ymin=0,ymax=170] {CHAP_6/exp/SFII-4.png};
\node[rotate=0,white] at (axis description cs:0.5,0.5) {\rotatebox{0}{Fiber}};
\node[rotate=0,white] at (axis description cs:0.5,0.87) {\rotatebox{0}{Matrix}};
\node[rotate=0,white] at (axis description cs:0.75,0.1) {\rotatebox{0}{Free surface}};
\draw[->,white] (axis description cs:0.85,0.15) -- (axis description cs:0.8,0.3);
\node[rotate=0,white] at (axis description cs:0.35,0.07) {\rotatebox{0}{0.5 mm}};
\addplot[white,line width = 3pt] coordinates {(5,8) (48,8) };
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:SFII_4_without}
\end{subfigure}
\hfill
\begin{subfigure}[c]{0.59\textwidth}
\centering

\begin{tikzpicture}
\begin{axis}[enlargelimits=false,
width = 1.*0.47\textwidth, height = 1.*0.69\textwidth, ticks=none, xlabel = {z-direction}, ylabel = {y-direction}, title = {Stress}, legend style={at={(0.5,-0.14)},anchor=north},]
\addplot[forget plot] graphics [xmin=0,xmax=47,ymin=0,ymax=69] {CHAP_6/exp/SFII-4_zoom.png};
\addplot[red, line width = 1pt] file[] {CHAP_6/exp/SFII-4_stress_gamma_1.36.dat};
\addplot[red,dotted, line width = 1pt] file[] {CHAP_6/exp/SFII-4_stress_gamma_0.65.dat};
\addplot[red,dashed, line width = 1pt] file[] {CHAP_6/exp/SFII-4_stress_gamma_0.3.dat};
\addplot[fill=red, opacity=0.3,area legend] file[] {CHAP_6/exp/SFII-4_stress_gamma_1.36_area.dat};
\addplot[white,line width = 3pt] coordinates {(3,65) (12,65)};
\node[rotate=0,white] at (axis description cs:0.5,0.925) {\rotatebox{0}{10 \unit{\micro\metre}}};
\legend{$\gamma = 1.36$,$\gamma = 0.65$,$\gamma = 0.3$}
\end{axis}
\end{tikzpicture}
\begin{tikzpicture}
\begin{axis}[enlargelimits=false,
width = 1.*0.47\textwidth, height = 1.*0.69\textwidth, ticks=none, xlabel = {z-direction}, ylabel = {y-direction}, title = {Energy}, legend style={at={(0.5,-0.14)},anchor=north},]
\addplot[forget plot] graphics [xmin=0,xmax=47,ymin=0,ymax=69] {CHAP_6/exp/SFII-4_zoom.png};
\addplot[blue, line width = 1pt] file[] {CHAP_6/exp/SFII-4_ener_gamma_1.6.dat};
\addplot[blue,dotted, line width = 1pt] file[] {CHAP_6/exp/SFII-4_ener_gamma_1.16.dat};
\addplot[blue,dashed, line width = 1pt] file[] {CHAP_6/exp/SFII-4_ener_gamma_0.8.dat};
\addplot[fill=blue, opacity=0.3,area legend] file[] {CHAP_6/exp/SFII-4_ener_gamma_1.6_area.dat};
\addplot[white,line width = 3pt] coordinates {(3,65) (12,65)};
\node[rotate=0,white] at (axis description cs:0.5,0.925) {\rotatebox{0}{10 \unit{\micro\metre}}};
\legend{$\gamma = 1.6$,$\gamma = 1.2$,$\gamma = 0.8$}
\end{axis}
\end{tikzpicture}

\caption{}\label{fig:SFII_4_with}
\end{subfigure}
\hfill

\caption{(a) Side view of the debonding just after debonding initiation. The debonded area is highlighted using a threshold and the white pixels. (b) Debonding shapes obtained with the CC and several brittleness numbers for both stress-based and energy-based shapes, corresponding to the red and blue shaded areas, superimposed on the side view.}\label{fig:SFII_crack_compare}
\end{figure}

Three brittleness numbers are selected and their corresponding debonding shapes are superimposed on the zoomed experimental side view, see Figure~\ref{fig:SFII_2_with} and Figure~\ref{fig:SFII_4_with}. On the one hand, there are energy-based shapes that correspond to the experimental debonding angle and length after initiation for a given range of brittleness numbers. For sample 1, experimental observations yield an upper limit of 0.6 as an acceptable $\gamma$, and 1.2 for sample 3, otherwise excessive debonding angles are encountered. These restrictions are applied to the inverse identification in the following. Energy-based shapes are able to reproduce the experimental debonding configuration, in terms of angle, length and shape. On the other hand, as pointed out previously, stress-based shapes significantly overestimate the debonding angle, regardless of $\gamma$, while accurately predicting the debonding length, for a range of $\gamma$. The upper limit of the acceptable $\gamma$ range based on the experimental debonding length can be determined as 0.56 for sample 1 and 0.65 for sample 3. Nevertheless, the shape provided by the stress isocontours does not capture the experimentally observed debonding zone, overestimating the debonding area whatever the brittleness number. This corroborates the fact that shapes based on stress isocontours are not optimal for the configuration studied, while not providing the optimum CC solution in terms of debonding shape. Similar observations are obtained for sample 2.

\section{Inverse identification of the interface fracture properties}\label{sec:INV}

The couples of tensile strength and critical ERR in mode I identified leading to equivalent remote stress ($\sigma^\infty$) to the experimental one can be identified. The identification is made on the basis of the energy-based debonding shapes with $\lambda = 0.13$ and $\mu = 1$. Since the mode II has no influence at initiation, the $\lambda$ and $\mu$ values are kept constants for the inverse identification as their influence is negligible. As a results, no shear fracture properties are identified. In addition, stress-based debonding shapes are also used for the inverse identification. Similarly, no shear fracture properties can be identified, see~\cite{girard_comparison_2025}. Once the range of fracture properties are identified, additional restrictions can be applied to the range of properties identified on the basis of the geometrical parameters obtained experimentally, whose values are given in Table~\ref{tab:exp_results}, and shown in Figure~\ref{fig:SFII_crack_compare}. The initiation surface provided by the CC can be either smaller than or equal to the arrest surface observed experimentally, as unstable debonding propagation can occur, see Section~\ref{section_3_2}. The initiation surface reverts to the arrest surface when initiation is driven by the maximum of $G_\mathrm{inc} / \overline{G_\mathrm{c}}$. However, the energy-based debonding shapes yield $G_\mathrm{inc} / \overline{G_\mathrm{c}}$ that increases monotonically for the brittleness numbers identified so that no stable debonding can occur after initiation, see Dormieux et al.~\cite{dormieux_remarks_2025} or Girard et al.~\cite{girard_numerical_2023} for the determination of the crack initiation surfaces bounds. Similarly, when the CC solution is given by the intersection of both stress and energy conditions, further debonding propagation takes place after initiation. Therefore, without the use of a high-speed camera, the experimentally observed debonding surface after initiation actually corresponds to the arrest surface that has potentially already propagated. The geometrical parameters shown in Table~\ref{tab:exp_results} therefore serve as an upper limit to those obtained using the CC and therefore allow the identified couples of fracture properties to be restricted, as explained in Section~\ref{sec:compare_exp}. It is worth recalling that, except for the debonding opening, the determination of the debonding angle and length is limited by the pixel size (see~\cite{uddin_simultaneous_2024}). The results derived from these two quantities may therefore be underestimated as sub-pixel phenomena may have occurred. Consequently, the maximum uncertainties in the calculation of the debonding opening is considered for the restriction of the identified fracture properties. The restricted identified properties are summarized in Table~\ref{tab:restri_prop_inv}.

\begin{table}[h]
\caption{Interface fracture properties identified for sample 1, 2 and 3 and a $\mu$ ratio of 1 and $\lambda = 0.13$.}\label{tab:restri_prop_inv}
\begin{tabular}{c|cc}
\hline
&&\\*[-2.7ex]
\hline
& \multicolumn{2}{c}{\textbf{Sample 1}} \\
\hline
&&\\*[-2.7ex]
\hline
$\sigma^\infty$ [\unit{\unit{MPa}}] & \multicolumn{2}{c}{7.25 $\pm$ 0.04} \\ \hline
& Stress-based & Energy-based \\ \hline $\gamma$ & 0.3 - 0.56 & 0.3 - 0.60 \\
$\sigma_\mathrm{c}$ [\unit{\unit{MPa}}] & 11.9 - 16.8 & 13.0 - 20.7 \\
$G_\mathrm{IC}$ [\unit{J/m^2}] & 2.7 - 4.2 & 3.9 - 5.8 \\ 
\hline
&&\\*[-2.7ex]
\hline
& \multicolumn{2}{c}{\textbf{Sample 2}} \\ 
\hline
&&\\*[-2.7ex]
\hline
$\sigma^\infty$ [\unit{\unit{MPa}}] & \multicolumn{2}{c}{8.03 $\pm$ 0.04} \\ \hline & Stress-based & Energy-based \\ \hline $\gamma$ & 0.3 - 0.5 & 0.3 - 0.74 \\
$\sigma_\mathrm{c}$ [\unit{MPa}] & 13.9 - 18.6 & 12.4 - 22.9 \\
$G_\mathrm{IC}$ [\unit{J/m^2}] & 3.3 - 4.9 & 4.8 - 8.2 \\ 
\hline
&&\\*[-2.7ex]
\hline
& \multicolumn{2}{c}{\textbf{Sample 3}} \\ 
\hline
&&\\*[-2.7ex]
\hline
$\sigma^\infty$ [\unit{MPa}] & \multicolumn{2}{c}{10.04 $\pm$ 0.04} \\ \hline & Stress-based & Energy-based \\ \hline $\gamma$ & 0.3 - 0.65 & 0.3 - 1.16 \\
$\sigma_\mathrm{c}$ [\unit{MPa}] & 14.8 - 23.3 & 11.1 - 29.8 \\
$G_\mathrm{IC}$ [\unit{J/m^2}] & 5.2 - 8.9 & 7.6 - 15.8
\end{tabular}
\end{table}

The debonding angle is not taken into account for restrictions applied to properties identified with the stress-based shapes. Overall, similar strength ranges are identified between the stress-based and energy-based shapes, regardless of the sample, in contrast to the critical ERR in mode I. Smaller critical ERR are identified using stress-based shapes due to less favorable energy conditions. The fracture properties identified are in line with those identified in~\cite{girard_numerical_2023}, although the epoxy system and glass fiber used are different.

\section{Conclusion}

The stress-based or energy-based debonding shapes are used to assess the fiber-matrix debonding initiation shape using the CC on an epoxy single glass fiber specimen under transverse tension. Energy-based shapes require additional effort to be determined because they require all possible shape configurations to be compared. On the basis of experimental observation, restrictions are made about the initiation location and the shape description. A power-law with an exponent of 1.5 is found to be optimal regardless of the critical ERR in mode II. Overall, the critical ERR in mode II has a negligible influence on the magnitude of the IERR to critical ERR ratio. The predicted debonding shape results in a smaller angle with larger length with increasing critical ERR in mode II, due to increased mode mixity in the vicinity of the fiber equator.

Energy-based shapes provide the optimal solution regardless of the brittleness number. The stress-based shapes overestimate the optimal CC solution by up to 15\% for either small ($\gamma < 1$) or large ($\gamma > 6$) brittleness numbers. Larger differences, up to 30\%, are encountered for intermediate brittleness numbers ($1 < \gamma < 6$). The energy condition causes this difference, the latter being less favorable using stress-based shapes, especially for intermediate brittleness numbers, where the energy condition triggers the CC solution. Energy-based debonding shapes provide a debonding topology consistent with the experimental observations, predicting the length, angle and the opening of the debonding compatible with those measured experimentally. Stress-based shapes are unable to provide an accurate shape with respect to experimental observation. It should be noted that the results obtained depend on the configuration studied (e.g., geometry, materials, boundary conditions). Therefore, the approach presented must be reproduced in each configuration in order to determine whether stress-based or energy-based shapes are optimal. In addition, the determination of energy-based shapes must be simplified using mathematical expressions, as their calculation is computationally expensive. Additionally, the debonding surface predicted by the CC corresponds to the initiation surface that is likely to propagate after initiation to the arrest surface. Further study of debonding propagation would be required, which could be approximated using an algorithm that iteratively varies the debonding front to obtain a constant stress intensify factor at the front for different increasing debonding surfaces.

Inverse identification is carried out on the basis of the stress-based and energy-based debonding shapes. Only the opening fracture properties are identified, since shear fracture properties slightly influence the CC solution for small surface such as the ones observed experimentally. The two approaches give similar fracture properties in terms of tensile strength, while the difference becomes more pronounced in terms of critical ERR. Admissible tensile strengths ranging from 11 to 30 \unit{MPa} and a mode I critical ERR of 4 to 16 \unit{J/m^2} are identified with the energy-based debonding shapes for the three samples. Further refinement of the interface fracture properties range could be obtained by determining the critical ERR based on the debonding propagation stage.

\section*{Dedication}

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\printbibliography

\end{document}