\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRMECA,Unicode,biblatex,published]{cedram}
\addbibresource{crmeca20250348.bib}
\newenvironment{noXML}{}{}%
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\sfrac#1#2{{#1}/{#2}}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\def\xsection{}
\def\xtbody{\noalign{\relax}}
\usepackage{etoolbox}
\usepackage{upgreek}
\usepackage[T1]{fontenc}
\usepackage{multirow} 
\providecommand{\algorithmname}{Algorithm}
\usepackage[ruled,vlined]{algorithm2e}
\newenvironment{inftab}{}{}
\def\mathbi#1{\text{\textbf{\textit{#1}}}}
\def\inlinefig#1{\includegraphics{#1}}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\def\jobid{crmeca20250348}
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\def\textemail#1{\href{mailto:#1}{#1}}
\csdef{Seqnsplit}{\\}
\def\botline{\\\hline}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
\def\0{\phantom{0}}
\def\mn{\phantom{$-$}}
\def\refinput#1{}
\def\back#1{}
\def\beqimg#1#2{}
\def\eeqimg#1#2{}
\DOI{10.5802/crmeca.338}
\datereceived{2025-04-25}
\daterevised{2025-09-27}
\dateaccepted{2025-10-28}
\ItHasTeXPublished
}
{\documentclass[crmeca]{article}
\usepackage{upgreek}
\usepackage[T1]{fontenc} 
\usepackage{moreverb}
\usepackage{stmaryrd}
\providecommand{\algorithmname}{Algorithm}
\usepackage[ruled,vlined]{algorithm2e}
\def\xtbody{\noalign{\relax}}
\let\beqimg\relax
\let\eeqimg\relax
\def\CDRdoi{10.5802/crmeca.338}
\def\centering{}
\def\newline{\unskip\break}
\def\hyperlink#1#2{#2}
\def\hypertarget#1#2{#2}
\def\sfrac#1#2{{#1}/{#2}}
}
\makeatother

\dateposted{2025-11-24}
\begin{document}

\begin{noXML}

\CDRsetmeta{articletype}{research-article}

\title{Crack propagation in concrete mortar under quasi-static tensile
loading: experimental insights and phase-field modeling}

\alttitle{Propagation des fissures dans le mortier de b\'{e}ton sous
charge de traction quasi-statique : aper\c{c}u exp\'{e}rimental et
mod\'{e}lisation par champ de phase}

\author{\firstname{Mehdi} \lastname{Boudouh}\CDRorcid{0000-0002-5573-2554}\IsCorresp}
\address{Laboratory of Development in Mechanics and Materials (LDMM),
University of Djelfa, PB 3117, Djelfa, Algeria}
\email[M. Boudouh]{mehdi.boudouh@univ-djelfa.dz}

\author{\firstname{Brahim Elkhalil} \lastname{Hachi}\CDRorcid{0000-0002-6672-746X}}
\addressSameAs{1}{Laboratory of Development in Mechanics and Materials
(LDMM), University of Djelfa, PB 3117, Djelfa, Algeria}
%\email{ek.hachi@univ-djelfa.dz}

\author{\firstname{Boulanouar} \lastname{Azzouzi}\CDRorcid{0000-0002-0732-8071}}
\addressSameAs{1}{Laboratory of Development in Mechanics and Materials
(LDMM), University of Djelfa, PB 3117, Djelfa, Algeria}
%\email{boulanouar.azzouzi@univ-djelfa.dz}

\author{\firstname{Bassem} \lastname{Zouari}\CDRorcid{0000-0001-7373-1099}}
\address{LA2MP Laboratory, National Engineering School of Sfax,
University of Sfax, Sfax 3029, Tunisia}
%\email{bassem.zouari@enis.tn}

\author{\firstname{Mohamed} \lastname{Haboussi}\CDRorcid{0000-0001-8901-6559}}
\address{Laboratoire LSPM Universit\'{e} Sorbonne Paris Nord UPR 3407
CNRS, 93430, Villetaneuse, France}
%\email{haboussi@lspm.cnrs.fr}

\shortrunauthors

\keywords{\kwd{Direct tensile tests}\kwd{Material
characterization}\kwd{Phase field model}\kwd{Crack
propagation}\kwd{Alternated minimization}\kwd{Finite element Methode}}

\altkeywords{\kwd{Essais de traction directs}\kwd{Caract\'{e}risation
des mat\'{e}riaux}\kwd{Mod\`{e}le de champ de phase}\kwd{Propagation
des fissures}\kwd{Minimisation altern\'{e}e}\kwd{M\'{e}thode des
\'{e}l\'{e}ments finis}}

%\thanks{}

\begin{abstract} 
This study investigates crack propagation in concrete mortar under
quasi-static tensile loading with imposed displacement. Several
rectangular cubic concrete specimens containing a predefined crack and
hole were cast and tested. The load--displacement curves and crack paths
were obtained experimentally and analyzed, enabling us to assess the
influence of the specimen age parameter on fracture behavior. The
evolution of mechanical properties as a function of age, in particular
the modulus of elasticity $E$ and Poisson's ratio $\vartheta$
(characterized by ultrasonic testing), was monitored. The fracture
energy $G_{c}$ was determined by direct tensile testing. These
parameters were fed into a phase field model implemented by the finite
element method in a computational code developed in
MATLAB{\textregistered}. The solution is based on a variational system
that treats the displacement and damage fields in an alternating
minimization scheme. The numerical results were compared in two
distinct frameworks: one in a logical domain and the other in a
physical domain. Finally, the results obtained in each of these
frameworks were compared with experimental results, revealing
relatively good agreement.
\end{abstract} 

\begin{altabstract} 
Cette \'{e}tude examine la propagation des fissures dans le mortier de
b\'{e}ton sous une charge de traction quasi-statique avec
d\'{e}placement impos\'{e}. Plusieurs \'{e}chantillons de b\'{e}ton
cubiques rectangulaires contenant une fissure et un trou
pr\'{e}d\'{e}finis ont \'{e}t\'{e} moul\'{e}s et test\'{e}s. Les
courbes charge-d\'{e}placement et les trajectoires des fissures ont
\'{e}t\'{e} obtenues exp\'{e}rimentalement et analys\'{e}es, ce qui
nous a permis d'\'{e}valuer l'influence du param\`{e}tre d'\^{a}ge de
l'\'{e}chantillon sur le comportement \`{a} la rupture. L'\'{e}volution
des propri\'{e}t\'{e}s m\'{e}caniques en fonction de l'\^{a}ge, en
particulier le module d'\'{e}lasticit\'{e} $E$ et le coefficient de
Poisson $\vartheta$ (caract\'{e}ris\'{e}s par des essais
ultrasoniques), a \'{e}t\'{e} surveill\'{e}e. L'\'{e}nergie de rupture
$G_c$ a \'{e}t\'{e} d\'{e}termin\'{e}e par des essais de traction
directs. Ces param\`{e}tres ont \'{e}t\'{e} int\'{e}gr\'{e}s dans un
mod\`{e}le de champ de phase mis en {\oe}uvre par la m\'{e}thode des
\'{e}l\'{e}ments finis dans un code de calcul d\'{e}velopp\'{e} dans
MATLAB{\textregistered}. La solution est bas\'{e}e sur un syst\`{e}me
variationnel qui traite les champs de d\'{e}placement et
d'endommagement dans un sch\'{e}ma de minimisation altern\'{e}e. Les
r\'{e}sultats num\'{e}riques ont \'{e}t\'{e} compar\'{e}s dans deux
cadres distincts : l'un dans un domaine logique et l'autre dans un
domaine physique. Enfin, les r\'{e}sultats obtenus dans chacun de ces
cadres ont \'{e}t\'{e} compar\'{e}s aux r\'{e}sultats
exp\'{e}rimentaux, r\'{e}v\'{e}lant une concordance relativement bonne.
\end{altabstract} 

\maketitle
{\vspace*{-.6pc}}

\end{noXML}

\section{Introduction}

Concrete plays a key role in engineering structures due to its
exceptional compressive strength and remarkable durability under
various environmental conditions. Despite these advantages, which make
it a popular choice among engineers and architects, concrete has two
notable drawbacks that can impact its overall performance: its
relatively low tensile strength and inherent brittleness. These
properties lead to cracking, which is one of the primary failure
mechanisms observed in nearly all brittle cementitious
materials~\cite{1,2}.

In 1961, Kaplan~\cite{3} conducted a series of experimental fracture
tests on concrete beams, applying Griffith's theory~\cite{4} to better
understand the complex behavior of cracks and fractures in solid
materials. This work marked the systematic introduction and development
of fracture mechanics and damage mechanics, providing a rigorous
framework for analyzing the non-linear properties of concrete
materials.

Bruhwiler et~al.~\cite{5}, performed a series of wedge-splitting tests
on concrete specimens with predefined cracks. They observed that
fracture toughness values remained unaffected by crack size,
demonstrating that the fracture behavior of concrete remained
consistent despite variations in crack dimensions.

Zhang et~al.~\cite{6} evaluated concrete toughness through
wedge-splitting tests and applied several analytical techniques,
including the simplified Green's function method, the simplified
equivalent cohesive force method, the four-term weighting function
method, and the Gauss-Chebyshev integral method. Their findings
revealed a strong agreement between the results of the wedge-splitting
tests and those obtained from each of the four analytical approaches.
Xu et~al.~\cite{7} developed a numerical method to derive the
stress--strain curve of concrete, incorporating the fracture process
zone observed during three-point bending tests.

Landis~\cite{8} investigated micro--macro fracture relationships in
concrete using Accoustic Emission (AE) and found that AE activity
increases sharply at the onset of microcracking, and that the spatial
distribution of AE sources predicts the eventual macrocrack trajectory.
To get more information about crack initiation and propagation, Dai
et~al.~\cite{9}, Alam et~al.~\cite{10}, Rouchier et~al.~\cite{11},
Verbruggen et~al.~\cite{12}, Sharma et~al.~\cite{13} conducted an
extensive study on the fracture process zone (FPZ) in concrete
employing both digital image correlation (DIC) and acoustic emission
(AE). They showed that AE could detect and localize the initiation of
microcracks within the specimen volume, while DIC provided full-field
strain and displacement measurements on the surface, enabling accurate
tracking of FPZ development.

Zhang et~al.~\cite{14} conducted an extensive study to capture crack
initiation, using the acoustic emission (AE) technique as a key
analytical tool, particularly during three-point bending tests on
concrete specimens with predefined cracks. Fang et~al.~\cite{15}
evaluated the performance characteristics related to the early crack
resistance of concrete, focusing on mixtures containing ternary mineral
components, using a thermal stress testing machine.

Guan et~al.~\cite{16} conducted a series of wedge splitting tests on
solid concrete specimens and found a significant correlation between
the age of the concrete and its mechanical \mbox{properties},
specifically noting that toughness and fracture energy increased as the
concrete aged. Refs.~\cite{17,18} conducted wedge splitting tests on
solid concrete specimens with predefined cracks and found that
specimens cast during the summer months exhibited higher fracture
energy values than those cast in the winter months. In addition, Xu
et~al.~\cite{19} conducted a detailed examination and found that
toughness and fracture energy were largely unaffected by changes in the
predefined crack length. Meanwhile, Wittmann et~al.~\cite{20} conducted
rigorous evaluations to determine the fracture parameters of concrete
using the tensile test, Chen et~al.~\cite{21} investigated the tensile
properties of concrete under high strain rates using the split
Hopkinson pressure bar technique.

Due to the challenges and complexities of conducting direct
experimental tensile tests on concrete materials, researchers often
turn to indirect methods as practical alternatives to effectively study
the tensile behavior of concrete Zheng et~al., Mesbah
et~al.~\cite{22,23}. Among the various indirect testing methods
introduced in the field is the Brazilian disc test, developed by Rocco
et~al.~\cite{24}, four-point bending tests by Wu et~al.~\cite{7}, and
the three-point bending test by Lok et~al.~\cite{25}. In addition,
research conducted by Lok et~al., Ruiz et~al.~\cite{25,26} involved
performing splitting tensile tests using the split Hopkinson pressure
bar (SHPB).

A considerable body of research has focused on the tensile response of
concrete using direct tension tests, despite the inherent experimental
difficulties associated with such methods. Van Mier and Van
Vliet~\cite{27} provided a comprehensive review of direct tension
testing techniques, identifying critical issues such as eccentricity,
secondary bending, and platen misalignment, and demonstrated the
necessity of servo-controlled systems to obtain reliable measurements
of tensile strength and fracture energy. Earlier, Reinhardt and
Cornelissen~\cite{28} conducted deformation-controlled uniaxial tension
tests under both static and cyclic loading on lightweight and
normal-weight concretes. Their results emphasized the significance of
post-peak softening behavior and yielded essential data for
constitutive model development under monotonic and cyclic load
conditions. More recently, Nouailletas et~al.~\cite{29} performed
cyclic tension--compression tests on notched specimens, combining
traditional displacement sensors with digital image correlation (DIC).
Their findings revealed that crack closure under compression partially
restored stiffness, while frictional interactions along crack faces
substantially contributed to inelastic strain and energy dissipation.

To address persistent experimental limitations, several improved
testing protocols have been introduced. Alhussainy et~al.~\cite{30}
proposed an innovative specimen geometry incorporating a reduced
cross-section and universal grips, ensuring systematic midspan failure
without slippage and thereby enabling reproducible stress--strain
characterization. In a related study, Alhussainy et~al.~\cite{31}
examined the direct tensile behavior of self-consolidating concrete
(SCC) using specially designed metallic grips. Their experiments
confirmed consistent fracture initiation at the specimen midsection and
showed that direct tensile strength values were systematically lower
than those obtained from indirect methods such as Brazilian or flexural
testing.

The mechanical response of concrete is governed not only by its
inherent heterogeneity arising from aggregate distribution, cement
matrix properties, and interfacial transition zones, but also by its
geometric configuration. The presence of initial defects (cracks,
voids, notches) or geometric discontinuities (section reductions,
eccentricities) modifies stress distribution and apparent stiffness,
thereby influencing crack initiation, propagation, and overall
resistance. Moreover, the age of the material is a decisive parameter,
as the progressive evolution of mechanical properties directly impacts
structural performance under load. This complexity highlights the need
to accurately characterize the mechanical parameters governing the
stiffness and deformation capacity of concrete. The estimation of
Young's modulus and Poisson's ratio has been the subject of numerous
studies, using both destructive and non-destructive approaches.
Compression or bending tests remain the traditional method for
determining the static modulus Es, but~they depend on experimental
conditions and involve the destruction of test
specimens~\cite{32,33,34}. To overcome these limitations,
non-destructive methods, particularly ultrasound, have been developed.
As early as the 1990s, Boumiz et~al.~\cite{35} and Martin
et~al.~\cite{36} demonstrated the relevance of P and S wave velocities
for characterizing hydration and elastic properties, while more recent
work~\cite{35,37} have confirmed their effectiveness in monitoring
setting and correlation with strength. Finally, several
studies~\cite{33,34,38,39} have established relationships between
static and dynamic modules, confirming the complementary nature of the
approaches: static tests provide references, while ultrasound allows
for continuous, non-destructive monitoring.

While experimental testing provides invaluable insights, it is often
time-consuming, costly, and limited in scope. Consequently, numerical
modeling has become a cornerstone of fracture research. A predominant
feature of most numerical models is their tendency to bypass the
singularities that appear in displacement fields as a result of crack
formation in structures. Various numerical methods for simulating crack
initiation and propagation have been founded and developed in recent
years. These methods can be classified into two main types, called
discrete crack models and diffuse crack models. In discrete crack
models, it is assumed that damage within the material is concentrated
on an exceptionally short length scale, allowing localized analysis of
crack behavior. Many researchers have used discrete crack models for
crack growth and propagation in solids. Examples include~\cite{40,41}
in the case of the brutal fracture and~\cite{42,43,44,45,46} in the
case of fatigue.

Diffuse crack models conceptualize the crack phenomenon as a zone of
damage distributed over a finite volume. These models are widely
adopted due to their relatively simple implementation in a calculation
code and their ability to handle complex crack configurations. Rather
than considering discrete cracks, these models are based on a diffuse
approximation that incorporates the effects of crack formation. The
phase field model belongs to this class of diffuse cracking models; it
extends Griffith's theory by simultaneously determining the
displacement field and the crack trajectory by minimizing the total
potential energy in a variational framework~\cite{47,48}. Bourdin
et~al.~\cite{49,50} have established the fundamental principles for the
numerical implementation of this approach. A damage phase field is
introduced to transform a clean crack into a regularized damage band,
whose width is controlled by a characteristic scaling parameter.

A central feature of phase-field fracture models is the introduction of
an internal length scale~$\ell$, which governs the thickness of the
diffused crack and ensures the well-posedness of the variational
formulation. Several strategies have been proposed for its calibration.
Early works~\cite{51,52} suggested experimental identification or
calibration through size-effect tests, though Iacono et~al.~\cite{53}
showed that this procedure does not lead to a unique solution. Later
studies~\cite{54,55,56} interpreted $\ell$ as a material parameter,
intrinsically linked to the tensile strength of a homogeneous medium.
More recently, Nguyen et~al.~\cite{57}, and Li et~al.~\cite{58}
emphasized the need to ensure consistency between $\ell$ and the
microstructural heterogeneity of cementitious composites. In
particular, mesoscale analyses of concrete and mortar have shown that
$\ell$ should be related to the characteristic size of sand grains or
aggregates in order to retain physical representativity~\cite{58,59}.
At the same time, Miehe et~al.~\cite{60}, Farrell and Maurini~\cite{61}
highlighted that $\ell$ also influences numerical stability,
conditioning, and convergence of the finite element solution.
Therefore, the choice of $\ell$ must balance physical consistency with
computational robustness, a point that is especially critical when
modelling granular quasi-brittle materials such as concrete.

In addition to the phase-field approach, several other continuous
regularized methods have been proposed in the literature to address the
challenges related to crack initiation, propagation, and localization
in quasi-brittle materials. Nonlocal models, of gradient or integral
type~\cite{55,62,63}, introduce an internal length scale that makes it
possible to account for size effects observed in concrete structures.
The Thick Level Set (TLS) method~\cite{64,65,66} combines a local
damage law with an imposed localization profile, enabling a natural
representation of crack initiation, coalescence, and branching. More
recently, the Lip-field approach~\cite{67,68,69} has been proposed as a
mathematically rigorous regularization technique, ensuring boundedness
of the damage gradient through a Lipschitz continuity constraint. These
phase field models share a common point: the introduction of an
internal length $\ell$ that regularizes the damage field and provides a
physically consistent description of fracture in quasi-brittle media.
In comparison with these methods, the phase field model also takes into
account the size effect while offering a numerically inexpensive
formalism given its well-established variational form with remarkable
simplicity in the treatment of nucleation and crack branching. All of
this led us to choose the phase field approach as our modeling method.

The numerical implementation of phase field models is generally carried
out within the \mbox{framework} of the finite element method (FEM), due
to the fact that the equations governing these~models are partial
differential equations, which can be solved efficiently using this
approach. The spatial discretization of the domain using finite
elements can be implemented \mbox{relatively} simply. Numerous studies
have documented finite element implementations of phase field models,
including the work of Amor et~al.~\cite{54}, Bourdin et~al.~\cite{70},
Lee et~al.~\cite{71}, Miehe et~al.~\cite{72}, Sargado et~al.~\cite{73},
Wu~\cite{74}, Kriaa et~al.~\cite{75}, Piska et~al.~\cite{76}, Yuanfeng
et~al.~\cite{77},{\break} Hachi et~al.~\cite{78}.

In this context, the present study adopts an integrated experimental
framework combining destructive and non-destructive testing in order to
better understand the evolution of the \mbox{mechanical} properties and
tensile behavior of mortar with predefined geometric heterogeneities
(cracks and voids). The evolution of Young's modulus and Poisson's
ratio was first monitored by ultrasound at 7, 14, and 28 days, allowing
repeated measurements on the same specimens and reducing variability
due to sample differences. In parallel, a calibration was performed
using three-point bending tests to establish a link between dynamic and
static properties, leading to a linear empirical equation of the type
$E_{\mathrm{s}}=a\times E_{\mathrm{d}}+b$ adapted to the studied
mortar. Quasi-static direct tensile tests were then used to analyze the
structural response and crack propagation. The experimental
characterization results show a gradual but non-regular increase in
stiffness as the specimen ages, reflecting a continuous densification
of the microstructure. The tensile tests confirm this evolution, with a
slight increase in the critical load between 7 and 14 days, followed by
a much more marked increase between 14 and 28 days. The experimental
originality of this work lies in combining the monitoring of changes in
mechanical \mbox{properties} with the age of mortar and the tensile response
of mortars on test specimens with predefined geometric heterogeneities
(cracks and voids), thus giving a more representative and integrated
view of the tensile behavior of cementitious mortars at different ages.

Numerically, this study adopts a specific implementation of the phase
field approach to modeling the crack propagation in the mortar under
study. The variational formulation has been decomposed into two coupled
equations, one for displacement and the other for damage, solved by an
explicit alternating minimization strategy implemented in a
MATLAB{\textregistered} code developed for this work. This methodology
ensures enhanced numerical stability and reproducibility of the
protocol, consistent with the experimentally measured mechanical
parameters $(E, \vartheta, G_{c})$. The choice of the regularization
parameter $\ell$ was physically justified by the sand grain~size,
\mbox{establishing} a direct correspondence between the microstructural
scale and the numerical model. Finally, an original double-domain
strategy was introduced: in the logical domain, the complete
experimental configuration is reproduced, while in the physical domain,
the mechanical effect of the tubes is replaced by a simplified
condition $(\Phi=0)$ at the edges. This dual modeling approach allows
for cross-validation and highlights the consistency between numerical
and experimental results. The numerical originality of this work thus
lies in the adaptation of the phase field to the experimental
conditions studied, paving the way for a better understanding of
cracking in cementitious mortars.

This article is structured as follows: Section~\ref{sec2} presents the
experimental procedure and protocol adopted, including a granulometric
analysis of the sand used, as well as a detailed description of the
mixture prepared and the specimens produced. A characterization of the
material parameters, notably the modulus of elasticity $E_{\mathrm{s}}$ and the
fracture energy of the material $G_{c}$, is then carried out using a
non-destructive characterization (NDT) technique via an ultrasonic
device. The direct tensile test adopted is then described in detail.
Section~\ref{sec3} summarizes the experimental results obtained during
the direct tensile tests, together with a comparative and explanatory
analysis. In Section~\ref{sec4}, the theoretical aspects and numerical
implementation of the phase-field approach adopted are presented.
Section~\ref{sec5} presents the numerical results, with comments and
comparisons with experimental observations. Finally, the most pertinent
conclusions are drawn in Section~\ref{sec6}, bringing this paper to a
close.

\section{Experimental procedure and protocol} \label{sec2}
\subsection{Granulometric analysis of sand} \label{sec2.1}

The test involves classifying the various grains in the specimen using
a series of sieves, arranged in a stack with progressively smaller mesh
sizes from top to bottom.

The sand specimen is placed in the top sieve, and the grain
classification is obtained by vibrating the sieve stack.

For this test, we use a series of sieves with square mesh by their
aperture classified as follows:
0.08--0.125--0.315--0.5--0.63--0.8--1--2 (mm).

The granulometric curve of sand is shown in Figure~\ref{fig1} according
to the standard NF P18-431 [Afn, P~18-431-1990].

\begin{figure}
\vspace*{-2pt}
\includegraphics{fig01}
\vspace*{-2pt}
\caption{\label{fig1}The granulometric curve of the sand.}
\end{figure}

The granulometric curve shows that our sand has a well-distributed
particle size range, with the majority of particles falling between
0.1~mm and 1~mm, classifying it as medium sand. The grading curve is
relatively smooth, indicating a uniform particle size distribution and
confirming that the sand is well-graded.

\subsection{Preparation of the mixture} \label{sec2.2}

The material used was a standard mortar with the following proportions:
22\% cement (Portland cement PC40 mixed with high alumina cement in a
4:1 ratio, conforming to ASTM C595), 66\% sand (particle size ${<}$
1~mm), and 12\% water (clean, impurity-free water), resulting in a
water/cement ratio of~0.55.

The dosages of the constituents used in this study are presented in
Table~\ref{tab1}. These standard specifications were selected to
minimize uncertainties in the mix by excluding other influencing
factors, such as additives.

\begin{table}%tab1
\caption{\label{tab1}Mixture proportions use\vspace*{-2pt}}
\begin{tabular}{cccc}
\thead
Component & Cement & Sand & Water\\
\endthead
Percentage & 22\% & 66\% & 12\%\\
Dosage & 459.53~Kg/m$^{3}$ & 1378.59~Kg/m$^{3}$ & 250.65~Kg/m$^{3}$
\botline
\end{tabular}
\vspace*{-2pt}
\end{table}

By adopting a standardized mortar formulation, composed of well-graded
sand (particle size range between 0.1 and 1~mm) and free from
admixtures, we deliberately minimized external sources of variability
in the mechanical behavior of the material. This methodological choice
ensures that the fracture response observed during the direct tensile
tests (see Section~\ref{sec3}) is mainly governed by the predefined
discontinuities introduced in the specimens and by the influence of
specimen age, which represent the central focus of the present study.

\subsection{Preparation of the test specimen} \label{sec2.3}

The test specimen used in this study is a rectangular parallelepiped
with dimensions of $W = 91$~mm (width), $H = 168$~mm (height), and $e =
40$~mm (thickness). The specimen features several geometric
characteristics: a circular cutout with a diameter of 28~mm located
near the center of the plate ($x = 51.1$~mm, $y = 71.4$~mm), and two
smaller cutouts with a diameter of 14~mm near the top and bottom,
positioned at ($x = 28$~mm, $y = 28$~mm) and ($x = 28$~mm, $y =
140$~mm), respectively. Additionally, the plate has an initial crack of
length $a_{0} = 14$~mm on its left side, located at $y = 91$~mm (see
Figure~\ref{fig2}).

\begin{figure}
\includegraphics{fig02}
\caption{\label{fig2}Geometry of the test specimen used in this study.}
\end{figure}

The steps for preparing the test specimen are detailed in
Appendix~\ref{appA}.

\subsection{Material characterization of the test specimen}
\label{sec2.4}
\subsubsection{Estimation of Young's modulus and Poisson's ratio using
non-destructive tests (ultrasound)} \label{sec2.4.1}

The estimation of mechanical properties through destructive testing
provides reliable results but requires the use of different specimens
for each measurement age. Even when mix proportions and curing
conditions are identical, slight variations related to particle size
distribution or microstructure may arise from one specimen to another,
introducing variability that can obscure the actual effect of age. To
overcome this limitation, the ultrasonic method was adopted. Its major
advantage lies in allowing repeated measurements on the same specimen
without altering it, thereby providing a more representative and
accurate monitoring of the evolution of Young's modulus and Poisson's
ratio as a function of age, while reducing experimental dispersion
between test series. In this study, this non-destructive ultrasonic
approach was implemented to characterize mortar at different curing
stages, namely at 7, 14, and 28 days. The targeted elastic \mbox{properties,}
namely the dynamic Young's modulus $E_{\mathrm{d}}$ and Poisson's
ratio~$\vartheta$, were determined using this method. The selected
propagation mode was direct transmission, as illustrated in
Figure~\ref{fig3}.\looseness=-1

\begin{figure}
\includegraphics{fig03}
\caption{\label{fig3}The ultrasonic testing principle.}
\end{figure}

The ultrasonic technique enables the measurement of the propagation
times of longitudinal pulses ($t_{\mathrm{p}}$, $P$-waves) and shear
pulses ($t_{\mathrm{s}}$, $S$-waves) through the mortar. These
propagation times are then used to calculate the corresponding wave
velocities, $V_{\mathrm{p}}$ and $V_{\mathrm{s}}$, which are required
for the determination of the dynamic Young's modulus $E_{\mathrm{d}}$
and Poisson's ratio $\vartheta$, using the equations presented
in~(\ref{eq2}) and~(\ref{eq3})~\cite{36,37}.

The ultrasonic device consists of a generator connected to a
transmitter and a receiver of electrical pulses, as illustrated in
Figure~\ref{fig3}. Piezoelectric crystals excited by impact generate
the pulses, while similar crystals in the receiver detect them. The
receiver, coupled with a timing unit, records the propagation time
between the emission and detection of the pulse. Based on these
measurements, the propagation velocities of longitudinal waves $V_{\mathrm{p}}$
and shear waves $V_{\mathrm{s}}$ are obtained as follows~\cite{36,37}:
{\begin{equation} \label{eq1}
V_{\mathrm{p}}=\frac{L}{(2\times t_{\mathrm{p}})},\quad V_{\mathrm{s}}=\frac{L}{(2\times t_{\mathrm{s}})}
\end{equation}}\unskip
where $L$ is the distance between the transducers, $t_{\mathrm{p}}$ and $t_{\mathrm{s}}$
are the propagation times of the longitudinal and shear pulses,
respectively. These velocities are then used to determine the dynamic
Young's modulus $E_{\mathrm{d}}$ and Poisson's ratio $\vartheta$ of the mortar,
according to the following relations~\cite{36,37}:
{\begin{eqnarray}
&  \displaystyle
\label{eq2}
E_{\mathrm{d}} = V_{\mathrm{p}}^{2}\gamma \frac{(1+\vartheta)(1-2\vartheta)}
{(1-\vartheta)}\Seqnsplit
&  \displaystyle
\label{eq3}
\vartheta = \frac{V_{\mathrm{p}}^{2}-2V_{\mathrm{s}}^{2}}{2(V_{\mathrm{p}}^{2}-V_{\mathrm{s}}^{2})}
\end{eqnarray}}\unskip
where $\gamma$ is the density of the specimen (kg/m$^{3}$),
$V_{\mathrm{p}}$ is the longitudinal wave velocity (m/s),
$V_{\mathrm{s}}$ is the shear wave velocity (m/s), $E_{\mathrm{d}}$ is
the dynamic Young's modulus (GPa), and $\vartheta$ is Poisson's ratio.

Equations~(\ref{eq2}) and~(\ref{eq3}) are valid under the assumption
that the material behaves as an elastic, homogeneous, and isotropic
medium at the macroscopic scale, and that ultrasonic waves propagate in
the linear regime without significant attenuation. Although
cementitious materials exhibit a complex microstructure, these
relations have been widely employed and validated in numerous studies
on pastes, mortars, and concretes for the evaluation of their
mechanical properties~\cite{35,36,37,39,79}. The use of a standard
mortar (sand, cement, water), as described in Section~\ref{sec2.2},
also helps to minimize the influence of heterogeneities associated with
coarse \mbox{aggregates} and to ensure a globally homogeneous and isotropic
behavior at the considered scale. In this context, the application of
Equations~(\ref{eq2}) and~(\ref{eq3}) is justified and consistent with
commonly accepted practices in the literature.

Table~\ref{tab2} summarizes the average values obtained from
non-destructive tests performed on nine specimens at 7, 14, and 28 days
of curing. For detailed results corresponding to each specimen, the
reader is referred to Tables~\ref{tabA1}--\ref{tabA3} in
Appendix~\ref{appA}.

\begin{table}%tab2
\caption{\label{tab2}Results of non-destructive tests}
\begin{tabular}{ccccc}
\thead
& \parbox[t]{2.3cm}{\centering Longitudinal wave $V_{\mathrm{p}}$ (m/s)\raisebox{-.43pc}{}} &
\parbox[t]{2.3cm}{\centering Shear wave $V_{\mathrm{s}}$ (m/s)} &
\parbox[t]{2.3cm}{\centering Poison's ratio average $\vartheta$} &
\parbox[t]{3.85cm}{\centering Dynamic Young's modulus average $E_{\mathrm{d}}$ (GPa)}\\
\endthead
7 days & $4.8942\times 10^{3}$\raisebox{.85pc}{} & $2.9773\times 10^{3}$ & 0.2062 & 44.6653\\
14 days & $5.3061\times 10^{3}$ & $3.2233\times 10^{3}$ & 0.2076 & 52.4081\\
28 days & $6.2760\times 10^{3}$ & $3.8022\times 10^{3}$ & 0.2100 & 73.0781
\botline
\end{tabular}
\end{table}

Table~\ref{tab2} highlights a progressive and significant increase in
the dynamic Young's modulus $E_{\mathrm{d}}$ with mortar age, while the
Poisson's ratio $\vartheta$ remains relatively stable. This evolution
is directly related to the ultrasonic velocities $V_{\mathrm{p}}$ and
$V_{\mathrm{s}}$, used in Equations~(\ref{eq2}) and~(\ref{eq3}).
Indeed, the gain in material stiffness is reflected by an increase in
propagation velocities, particularly the longitudinal wave velocity
$V_{\mathrm{p}}$, which mechanically induces an elevation
of~$E_{\mathrm{d}}$.

At 7 days, the dynamic Young's modulus $E_{\mathrm{d}}$ reaches approximately
45~GPa, indicating a still limited stiffness of the material, which can
be explained by the relatively low measured wave velocities $V_{\mathrm{p}}$ and
$V_{\mathrm{s}}$ (see Tables~\ref{tabA1}--\ref{tabA3} in Appendix~\ref{appA}).
Between 7 and 14 days, the increase in $E_{\mathrm{d}}$ remains moderate
$({{\approx} {+}}15\%)$, corresponding to a slight variation in ultrasonic
velocities and reflecting a gradual consolidation of the microstructure
with a limited gain in compactness. In contrast, between 14 and 28
days, the increase in $E_{\mathrm{d}}$ is much more pronounced $({{\approx}
{+}}35\text{--}40\%)$ and is accompanied by a significant rise in both
$V_{\mathrm{p}}$ and $V_{\mathrm{s}}$ (see Tables~\ref{tabA1}--\ref{tabA3} in
Appendix~\ref{appA}). This evolution indicates a substantial reduction
in porosity and improved matrix continuity, directly enhancing load
transfer and thus the elastic properties.

Furthermore, the Poisson's ratio, determined according to
Equation~(\ref{eq3}), remains nearly constant throughout the entire
period studied, indicating that the increase in Young's modulus is
mainly driven by the joint growth of $V_{\mathrm{p}}$ and $V_{\mathrm{s}}$, without
significant variation in their~ratio. For~\mbox{detailed} information
on the individual values of $V_{\mathrm{p}}$, $V_{\mathrm{s}}$, $E_{\mathrm{d}}$, and $\vartheta$
for each specimen, the reader is referred to
Tables~\ref{tabA1}--\ref{tabA3} in Appendix~\ref{appA}.

Particular attention has been given to the comparison between static
Young's modulus $E_{\mathrm{s}}$ and dynamic Young's modulus $E_{\mathrm{d}}$ for
different types of rocks, with numerous empirical formulas reported in
the literature~\cite{34,38}. These formulas are either linear, as shown
in Equation~(\ref{eq4}), or nonlinear, as shown in
Equation~(\ref{eq5}).
{\begin{eqnarray}
&  \displaystyle
\label{eq4}
E_{\mathrm{s}}=a\times E_{\mathrm{d}}+b\Seqnsplit
& \displaystyle 
\label{eq5} 
E_{\mathrm{s}}=a\times (E_{\mathrm{d}})^{b}.
\end{eqnarray}}\unskip
The nonlinear relationship Equation~(\ref{eq5}) is generally considered
more representative for materials and rocks exhibiting significant
microcracking or relatively high porosity, cases in which $E_{\mathrm{s}}$
increases more slowly than $E_{\mathrm{d}}$ ($b<1$). In contrast, the linear
form Equation~(\ref{eq4}) remains applicable to more compact and
homogeneous materials. In the present study, and for a standard mortar
as described in Section~\ref{sec2.2}, the linear form
Equation~(\ref{eq4}) was adopted, as it is widely used in the
literature for both rocks and cementitious
materials~\cite{80,81,82,83}. In all cases, the empirical coefficients
$a$ and $b$ must be experimentally determined for each mortar
formulation to ensure the reliability and relevance of the correlation.

The application of Equation~(\ref{eq4}) requires the definition of the
empirical parameters $a$ and $b$, which are strongly dependent on the
type of material investigated. To establish a representative relation
for the mortar under consideration, a specific experimental calibration
was undertaken. For this purpose, the static Young's modulus $E_{\mathrm{s}}$
was determined from three-point bending tests performed on nine
prismatic specimens with dimensions of 160 ${\times}$ 40 ${\times}$
40~mm$^{3}$ at the age of 14 days, in accordance with the Russian State
Standard GOST 310.4-81~\cite{32,33}. The test was conducted under a
monotonically applied load of 2.5~kN, as illustrated in
Figure~\ref{fig4}.

\begin{figure}
\includegraphics{fig04}
\caption{\label{fig4}Experimental setup of the three-point bending test
performed on prismatic specimens with dimensions of 160 ${\times}$ 40
${\times}$ 40~mm$^{3}$ at the age of 14~days.}
\end{figure}

The choice of 14 days for this calibration is justified by its
representativeness in the microstructural evolution of the mortar.
Indeed, this intermediate age makes it possible, on the one hand, to
reduce the influence of the high mechanical variability generally
observed at early ages (7~days), and on the other hand, to avoid the
quasi-complete stabilization effects of mechanical properties typically
observed at 28~days.

The static Young's modulus was then calculated from the mean
force--deflection curves, using the classical relation~\cite{32,33}:
{\begin{equation} \label{eq6}
E_{\mathrm{s}}=\frac{l^{3}}{4bh^{3}}\frac{\Delta F}{\Delta \delta}
\end{equation}}\unskip
where $l$ is the span between supports (100~mm), $b$ and $h$ are the
dimensions of the cross-section (40~mm ${\times}$ 40~mm), $F$ is the
applied load (kN), and $\delta$ is the measured deflection (mm).

The average load and deflection results from nine specimens, taken at
selected points within the linear elastic range, are presented in
Table~\ref{tab3}.

\begin{table}%tab3
\caption{\label{tab3}Experimental results of the three-point bending
test at 14 days: average values of the applied load, the corresponding
deflection, and the static Young's modulus $E_{\mathrm{s}}^{(14)}$}
\begin{tabular}{cccc}
\thead
\parbox[t]{1.9cm}{\centering The load $F$~(kN)} &
\parbox[t]{2.6cm}{\centering The~deflection $\delta$~(mm)} &
\parbox[t]{2.7cm}{\centering The slope ${\Delta F}/{\Delta \delta}$~(KN/mm)} &
\parbox[t]{3cm}{\centering The static Young's modulus $E_{\mathrm{s}}^{(14)}$\raisebox{-.43pc}{}}\\
\endthead
0.229 & 0.03878 & \morerows{8}{5.905105724600306} &
\morerows{8}{0.5764891}\\
0.458 & 0.07756 &&\\
0.687 & 0.11634 &&\\
0.916 & 0.15512 &&\\
1.145 & 0.19390 &&\\
1.374 & 0.23268 &&\\
1.603 & 0.27146 &&\\
1.832 & 0.31024 &&\\
2.060 & 0.34902 &&
\botline
\end{tabular}
\end{table}

In parallel, the dynamic Young's modulus $E_{\mathrm{d}}$ was measured
on the same specimens using the ultrasonic method, which provided a
representative average value at this age. The combined introduction of
$E_{\mathrm{s}}$ and $E_{\mathrm{d}}$ values into the linear empirical
equation (Equation~(\ref{eq4})) then allowed the determination of the
parameters~$a$ and~$b$. Considering the physical condition that
$E_{\mathrm{s}}=0$ when $E_{\mathrm{d}}=0$, it was deemed appropriate
to set $b=0$, thereby simplifying the relation to a direct
proportionality between the two moduli: $E_{\mathrm{s}}=a\times
E_{\mathrm{d}}$.

The proportionality coefficient $a$ was experimentally determined to be
0.011 for the mortar under investigation. Accordingly, the specific
relation obtained is:
{\begin{equation} \label{eq7}
E_{\mathrm{s}}=0.011\times E_{\mathrm{d}}.
\end{equation}}\unskip
This calibrated empirical equation was subsequently employed to
estimate the static Young's modulus $E_{\mathrm{s}}$ at other ages (7 and
28~days) from the available dynamic measurements. It thus provides a
reliable tool to non-destructively monitor the evolution of the
mortar's elastic modulus as a function of age.

Figure~\ref{fig5} illustrates the evolution of the average static
Young's modulus $E_{\mathrm{s}}$ with curing age. For detailed $E_{\mathrm{s}}$ values
corresponding to each specimen, calculated from Equation~(\ref{eq7}),
the reader is referred to Tables~\ref{tabA1}--\ref{tabA3} in
Appendix~\ref{appA}.

\begin{figure}
\vspace*{1pt}
\includegraphics{fig05}
\vspace*{1pt}
\caption{\label{fig5}Evolution of the average static Young's modulus
with curing age (7, 14, and 28~days).}
\vspace*{1pt}
\end{figure}

The evolution of Young's modulus with age, illustrated in
Figure~\ref{fig5}, highlights the progressive increase in material
stiffness during curing. The static modulus $E_{\mathrm{s}}$ exhibits a trend
similar to $E_{\mathrm{d}}$, since it is linearly proportional to it. More
specifically, the value of Young's modulus increases from approximately
0.49~GPa at 7 days to 0.58~GPa at 14 days. Between 14 and 28~days, this
increase becomes more pronounced, reaching 0.80~GPa at 28~days. This
evolution reflects the progressive densification of the microstructure,
with improved compactness and matrix continuity, which directly
enhances material stiffness. The histogram, complemented by the trend
curve, clearly illustrates this behavior and emphasizes that the
development of mechanical properties is not linear but intensifies
after two weeks of curing.

\subsubsection{Estimation of fracture energy} \label{sec2.4.2}

The fracture Energy $G_{c}~(\mathrm{N}/\mathrm{mm})$ is a crucial
parameter in the phase field model. For normal concrete mortar,
$G_{c}~(\mathrm{N}/\mathrm{mm})$ is determined from the fracture
resistance $W_{\mathrm{F}}$ as follows (see~in~\mbox{\cite{84,85,86,87}}):
{\begin{equation} \label{eq8}
G_{c}=\frac{W_{\mathrm{F}}+M\times g\times u_{0}}{A}.
\end{equation}}\unskip
In this equation, we define $M$ as the weight of the specimen (kg); $g$
as the acceleration due to gravity 9.81 (m/s$^2$); $u_{0}$ as the
critical imposed displacement (m); $A$ as the cross-sectional area at
the initial crack ($\mathrm{m}^{2}$). The fracture resistance $W_{\mathrm{F}}$
represents the area under the load--displacement curve at the maximum
value of~$u_{0}$.

The expression adopted for the estimation of the fracture energy
$G_{c}$ (Equation~(\ref{eq8})) is derived from the work-of-fracture
method in the framework of quasi-static tensile tests (direct tension,
three-point bending, or splitting test). Its validity relies on several
assumptions that must be specified. First, the loading must be
monotonic and sufficiently slow so that dynamic and rate effects are
negligible. Second, crack propagation is assumed to occur in pure
mode~I (opening), without significant contribution from shear or
mixed-mode effects. The material is considered homogeneous,
well-compacted, and isotropic at the macroscopic scale, with a
well-defined notch or initial crack to ensure reproducible initiation.
Finally, this formulation is restricted to monotonic loading up to
crack propagation and does not account for cyclic, viscoelastic, or
size effects. Under these conditions, the equation provides a
consistent estimation of the fracture energy of cement-based materials,
as demonstrated in the seminal works of Hillerborg et~al.\ and RILEM TC
50-FMC~\cite{86,87}.

For the mortar investigated in this study, the fracture energy
calculated using Equation~(\ref{eq8}) is given by:
{\begin{equation} \label{eq9}
G_{c}=\frac{0.2+1.2\times 9.81\times 0.25\times 10^{-6}}{40\times 70
\times 10^{-6}}=72.48~(\mathrm{N}/\mathrm{m}).
\end{equation}}\unskip

\subsection{The principle of the test adopted} \label{sec2.5}

Direct tensile tests were performed using an MTS hydraulic testing
machine (capacity 100~kN), as illustrated in Figure~\ref{fig6}. The
specimens were fixed by means of rods inserted into dedicated~holes,
ensuring symmetric alignment between the upper and lower grips. Prior
to each test, a~slight preload was applied to eliminate any residual
play and to guarantee the absence of initial stresses. The tensile test
was then conducted under displacement control: a total displacement of
1.5~mm was progressively applied over a duration of 3000~s (loading
velocity 0.0005~mm/s), thereby ensuring quasi-static conditions
suitable for accurately capturing the brittle behavior of the specimens
up to failure. The testing machine was equipped with a load cell with a
precision of ${\pm}$0.5\% and a displacement transducer with a
resolution of 0.001~mm, ensuring sufficient sensitivity to detect fine
variations in load and displacement.

\begin{figure}
\includegraphics{fig06}
\caption{\label{fig6}(a)~The direct tensile test; (b)~crack propagation
in specimen.}
\end{figure}

\section{Experimental results} \label{sec3}

In this section, the experimental results obtained during the test
adopted in Figure~\ref{fig6} are presented as follows.

Tables~\ref{tab4}, \ref{tab5}, and~\ref{tab6} present the experimental
results of the tensile tests conducted at 7, 14, and 28~days,
respectively. The load--displacement curves show that specimens of the
same age exhibit a consistent evolution, with only slight variations in
the critical load and displacement values, as detailed in
Tables~\ref{tab4}--\ref{tab6}. Analysis of these results indicates that
the differences observed, although modest, exceed the limit of
instrumental sensitivity and are likely related to the intrinsic
heterogeneity of the mortar (component distribution, initial
microcracks, local porosity), as well as slight variations in specimen
positioning or alignment. Nevertheless, the overall consistency of the
values obtained confirms the good reproducibility of the experimental
method.

\begin{table}%tab4
\caption{\label{tab4}Experimental results obtained at 07 days of age
(a)~damage evolution in specimen~01; (b)~damage evolution in
specimen~02; (c)~load--displacement curves}
\begin{tabular}{c}
\tbody
\raisebox{14.5pc}{}\inlinefig{fx01}
\botline
\end{tabular}
{\vspace*{.8pc}}
\end{table}

\begin{table}%tab5
\caption{\label{tab5}Experimental results obtained at 14 days of age
(a)~damage evolution in specimen~03; (b)~damage evolution in
specimen~04; (c)~load--displacement curves}
\begin{tabular}{c}
\tbody
\raisebox{15.2pc}{}\inlinefig{fx02}
\botline
\end{tabular}
\end{table}

\begin{table}%tab6
\caption{\label{tab6}Experimental results obtained at 28 days of age
(a)~damage evolution in specimen~05; (b)~damage evolution in
specimen~06; (c)~load--displacement curves}
\begin{tabular}{c}
\tbody
\raisebox{15.2pc}{}\inlinefig{fx03}
\botline
\end{tabular}
\end{table}

The structure behaves linearly up to the maximum loading. After that, a
brutal drop can be identified when the crack begins to propagate.

Figure~\ref{fig7} compares the effort--displacement curve as a function
of age factor. The evolution of the maximum load $F_{\max}$ as a
function of age exhibits two distinct phases. Between 7 and 14~days,
the increase remains limited (${+}$12\%), reflecting a progressive but
still incomplete consolidation of the microstructure. In contrast,
between 14 and 28~days, the progression becomes much more pronounced
(${+}$45\%), indicating a significant reduction in porosity and an
improved continuity of the matrix, which directly enhances stress
transmission.

\begin{figure}
{\vspace*{-2pt}}
\includegraphics{fig07}
{\vspace*{-2pt}}
\caption{\label{fig7}Comparison of the load--displacement curves at
ages 7, 14, and 28~days.}
{\vspace*{-2pt}}
\end{figure}

This trend is consistent with the results obtained for the Young's
moduli (Section~\ref{sec2.4}). In particular, the dynamic modulus
$E_{\mathrm{d}}$ shows a moderate increase between 7 and 14~days
(${{\approx}{+}}15\%$), followed by a marked acceleration between 14 and
28~days (${{\approx}{+}}35$--40\%), in line with the increase in ultrasonic
velocities. The static modulus $E_{\mathrm{s}}$ evolves in the same manner.
Thus, the correlation between mechanical strength, static modulus, and
dynamic modulus confirms that the significant increase observed between
14 and 28~days does not result from an experimental artifact, but
rather reflects a genuine microstructural transformation.

\section{Phase field modeling} \label{sec4}

The phase-field model has proved to be a reliable and practical
approach for predicting crack initiation and propagation. In this
section, we present the main equations of the phase-field approach to
brittle fracture, as well as its implementation by the finite element
method under our calculation code developed on MATLAB{\textregistered}.
The modeling strategy followed to reproduce the experimental tensile
tests is also detailed and validated.

\subsection{Theoretical aspects} \label{sec4.1}

In the present study, we have chosen to adopt and use the theoretical
framework defined by Francfort and Marigo~\cite{47,88,89}, with the aim
of addressing and overcoming the shortcomings of Griffith's 1998
theory~\cite{4}. They express the phenomenon fracture for a linear
elastic domain $\Omega$ that may contain a set of cracks $\Gamma$ as a
competition between the surface energy required for crack propagation
$\Pi_{\mathrm{s}}$ and the elastic energy $\Pi_{\mathrm{e}}$ stored in
the material in the form of a variational theory as a function of the
displacement field~$u$.

The total potential energy of the $\Omega$ domain can be written as:
{\begin{eqnarray}
& \displaystyle
\label{eq10} 
\Pi_{\mathrm{p}}(u,\Gamma) = \Pi_{\mathrm{e}}(u,\Gamma)+\Pi_{\mathrm{s}}
(u,\Gamma)-\mathcal{W}^{\mathrm{ext}}(u)\Seqnsplit
&  \displaystyle
\label{eq11}
\Pi_{\mathrm{e}}(u,\Gamma) = \int_{\Omega}W_{0}(\varepsilon
(u),\Gamma)\,\mathrm{d}\Omega\Seqnsplit
&  \displaystyle
\label{eq12}
\Pi_{\mathrm{s}}(\Gamma) = \int_{\Gamma}G_{c}H^{n-1}(\Gamma)\,
\mathrm{d}\Gamma\Seqnsplit
&  \displaystyle
\label{eq13}
\mathcal{W}^{\mathrm{ext}}(u) = \int_{\Omega}f\cdot
u\,\mathrm{d}\Omega +\int_{\partial \Omega_{u}} t\cdot u\,\mathrm{d}
\Gamma\Seqnsplit
&  \displaystyle
\label{eq14}
\Pi_{\mathrm{p}}(u,\Gamma) = \int_{\Omega}W_{0}(\varepsilon
(u), \Gamma)\,\mathrm{d}\Omega -\int_{\Omega}f\cdot u\,\mathrm{d}
\Omega -\int_{\partial \Omega_{u}} t\cdot u\,\mathrm{d}\Gamma +
\int_{\Gamma}G_{c}H^{n-1}(\Gamma)\,\mathrm{d}\Gamma
\end{eqnarray}}\unskip
where $\varepsilon (u)={1}/{2}[\nabla \otimes u+u\otimes \nabla]$ is
the symmetrical tensor of small strains, $t=\sigma \cdot n$, $\sigma$
and $n$ being respectively the stress tensor and the exterior normal to
the boundary $\partial \Omega_{u}$, $H^{n-1}$ is the $(n-1)$
dimensional Hausdorff measure giving the surface associated with
$\Gamma$, $W_{0}(\varepsilon (u))$ is the elastic energy density.

In the present work, the formulation is restricted to the case of
displacement-controlled loading (Dirichlet conditions). In this
configuration, the external work term $\mathcal{W}^{\mathrm{ext}}(u)$
vanishes, since at each loading step the displacement gradient is zero.
Conversely, when force-controlled loading (Neumann conditions) is
considered, this term must be explicitly accounted for, which may
induce a loss of convexity of the functional and lead to numerical
convergence difficulties toward a global minimum in the alternate
minimization scheme, as noted by Amor~\cite{90}. This restriction
therefore corresponds to an implementation choice motivated both by
numerical stability considerations and by the aim of faithfully
reproducing the experimental conditions described in
Section~\ref{sec2.5}, which were conducted under
displacement-controlled loading.

The total energy function is expressed by:
{\begin{equation} \label{eq15}
\Pi_{\mathrm{p}}(u,\Gamma)=\int_{\Omega}W_{0}(\varepsilon (u),\Gamma)\,
\mathrm{d}\Omega+\int_{\Gamma}G_{c}H^{n-1}(\Gamma)\,\mathrm{d}\Gamma.
\end{equation}}\unskip
Considering all possible cracking states, the evolution of the crack
set $\Gamma$ as well as the unknown displacements $u$ is determined by
the global minimization of functional (\ref{eq15}), under the
irreversibility constraint, as follows:
{\begin{equation} \label{eq16}
\Pi_{\mathrm{p}}(u_{i},\Gamma_{i})=\min_{\Gamma_{t+\Delta t}\supseteq
\Gamma_{t}} \Pi_{\mathrm{p}}(u,\Gamma).
\end{equation}}\unskip
In the case of brittle fracture, the fracture energy $G_{c}$ can be
interpreted as Griffith's critical energy restitution rate. The
challenge of solving the surface integral directly arises from the fact
that the new $\Gamma(t)$ surfaces are unfamiliar at first sight.
Conversely, classical discretization methods fail to produce a
satisfactory resolution.

The main difficulty in minimizing the previous variational formulation,
given by functional (\ref{eq15}), lies in the need to anticipate all
possible cracking states. In other words, it is necessary to construct
discontinuous displacement fields on arbitrary discontinuity surfaces,
which is beyond the capabilities of conventional finite element
methods. To make the problem tractable, Bourdin et~al.~\cite{48}
adopted an alternative approach, where minimization is performed on a
regularized approximation of the Griffith functional. This
regularization introduces smoothed displacement discontinuities,
ensuring the continuity of the displacement field $u$ across the entire
domain. This approach is inspired by the work of Ambrosio and
Tortorelli~\cite{91}, who solved a similar problem in the context of
image segmentation, in connection with the function of Mumford and
Shah~\cite{92}.

From the mechanical point of view, this regularization approach makes
it possible to reformulate fracture models using a phase field model,
thus introducing a parameter $\phi(x,t)\in [0,1]$. It is based on
replacing the sharp discontinuity in the displacement field,
representing the crack, with a gradual transition within small zones
where the displacement gradient is high (damaged zones). The size of
these zones is regulated by a parameter~$\ell$, while their
characteristics are determined by the damage variable $\phi (x,t)\in
[0,1]$, as illustrated in Figure~\ref{fig8}.

\begin{figure}
\includegraphics{fig08}
\caption{\label{fig8}A solid body $\Omega$ with a conventional sharp
crack $\Gamma$ and diffuse crack band in phase-field $\Phi(x)$.}
\end{figure}

Since the crack is represented as a damaged zone, an exponential decay
function is used to approximate its irregular topology. In one
dimension, in order to represent a crack located at $x=0$, a specific
formulation of the second-order phase field is adopted~\cite{48,72}:
{\begin{equation} \label{eq17}
\Phi(x)=\exp \left(\frac{{-}|x|}{2\ell}\right).
\end{equation}}\unskip
In Equation~(\ref{eq17}), $\phi (x,t)\in [0,1]$, with 0 representing
undamaged material and 1 indicating complete material damage. This
equation represents the solution of a homogeneous ordinary differential
equation~\cite{55,72,93}:
{\begin{equation} \label{eq18}
\Phi''(x)-\frac{1}{2\ell}\phi(x)=0\quad\text{in }\Omega.
\end{equation}}\unskip
The differential equation (\ref{eq18}) is solved under Dirichlet-type
boundary conditions:
{\begin{equation} \label{eq19}
\left\{\begin{array}{@{}c@{}}
\Phi(0)=1,\quad \phi'(0)=0\vspace*{3pt}\\
\displaystyle\lim_{x\rightarrow +\infty} \phi(x) = \lim_{x\rightarrow
-\infty} \phi(x)=0\quad\text{and},\vspace*{3pt}\\
\displaystyle\lim_{x\rightarrow +\infty} \phi'(x) = \lim_{x\rightarrow
-\infty} \phi'(x)=0.
\end{array}\right.
\end{equation}}\unskip
The Equation~(\ref{eq18}) represents the strong formulation
corresponding to the variational problem:
{\begin{equation} \label{eq20}
\Phi (x)=\mathrm{Arg}(\inf I(\phi))
\end{equation}}\unskip
were
{$$
I(\phi)=\frac{1}{2}\int_{\Omega}(\phi^{2}+(2\ell)^{2}G_{c} \nabla
\phi^{2})\,\mathrm{d}\Omega\quad\text{were }\mathrm{d}\Omega =
\Gamma_{\ell}\,\mathrm{d}x.
$$}\unskip
Therefore, we introduce $\Gamma_{\ell}(\phi)$ to describe the
second-order phase-field approach, defined as follows:
{\begin{equation} \label{eq21}
\Gamma_{\ell}(\phi)=\frac{1}{2\ell}I(\phi)=\frac{1}{4\ell}
\int_{\Omega}(\phi^{2}+(2\ell)^{2}\nabla \phi^{2})\,\mathrm{d}\Omega.
\end{equation}}\unskip
The choice of the coefficients of $\phi(x)$ and its derivatives in
Equation~(\ref{eq18}) is selected so that:
{\begin{equation} \label{eq22}
\int_{-\infty}^{+\infty}G_{c} \Gamma_{\ell}(\phi)\,\mathrm{d}x 
\approx \int_{\Gamma}G_{c}\Gamma=G_{c}.
\end{equation}}\unskip
Thus, in the context of the second-order phase field, we obtain:
{\begin{equation} \label{eq23}
\Pi_{\mathrm{s}}(\phi)=\int_{\Gamma} G_{c}\,\mathrm{d}\Gamma =
\frac{\mathcal{G}_{c}}{4\ell}\int_{\Omega}(\phi^{2}+(2\ell)^{2}
\nabla \phi^{2})\,\mathrm{d}\Omega.
\end{equation}}\unskip
As proposed by Bourdin et~al.~\cite{31,32}, we approximate the surface
energy as follows:
{\begin{equation} \label{eq24}
\Pi_{\mathrm{s}}(\phi)=\int_{\Omega} G_{c}\left(\frac{\phi^{2}}{4\ell}+\ell
\nabla \phi^{2}\right)\,\mathrm{d}\Omega
\end{equation}}\unskip
were the choice of the internal length $\ell$ in this study was guided
by the need to ensure good reproducibility of the experimental results.
This parameter must remain consistent with the microstructural
heterogeneity of the mortar~\cite{58,59}, while also complying with
numerical discretization constraints. In particular, the condition
$h_{\mathrm{e}}\ll \ell$ must be satisfied in order to guarantee an accurate
resolution of the crack front and the stability of the
simulations~\cite{60,94}.

$G_{c}$ represents the fracture toughness, defined as the energy
required to generate one unit of cracked surface.

As a means of modeling the loss of material stiffness in the fracture
zone, i.e.\ approximating the phase field of the damaged zone, we adopt
the approach of Bourdin et~al.~\cite{48,49,50}, who proposed an
isotropic formulation to model the progressive degradation of strain
energy and stiffness during crack propagation. Elastic energy is then
defined by the following expression:
{\begin{equation} \label{eq25}
\Pi_{\mathrm{e}}(u,\phi)=\int_{\Omega}[(D(\phi)+\varepsilon_{\ell})W_{0}
(\varepsilon [u])]\,\mathrm{d}\Omega 
\end{equation}}\unskip
were
{\begin{equation} \label{eq26}
W_{0}(\varepsilon [u])=\tfrac{1}{2}\sigma : \varepsilon(u)=\tfrac{1}
{2}\varepsilon : \mathbb{C} : \varepsilon =\tfrac{1}{2}\lambda\,
\mathrm{tr}^{2}\varepsilon +\mu\,\mathrm{tr}(\varepsilon^2).
\end{equation}}\unskip
From the potential of the phase field model, the stress can be
calculated by:
{\begin{equation} \label{eq27}
\sigma(u,\phi)=D(\phi)\frac{\partial W_{0}(\varepsilon [u])}
{\partial \varepsilon [u]}=(1-\phi)^{2}\frac{\partial W_{0}
(\varepsilon [u])}{\partial \varepsilon [u]}
\end{equation}}\unskip
$W_{0}(\varepsilon [u])$ represents the elastic energy density
function, while $\mathbb{C}$ is the elasticity tensor, $\varepsilon$
and $\sigma$ designate the strain and stress tensors respectively. The
parameters $\lambda$ and $\mu$ are the Lame coefficients.

The initiation and propagation of cracks are fully described by a
problem of minimizing a regularized energy functional. The latter is
parameterized by $\ell$ and expressed according to the following
equation:
{\begin{equation} \label{eq28}
\Pi_{\mathrm{t}}(u,\phi)=\int_{\frac{\Omega}{r}}\frac{1}{2}(A(\phi)+
\varepsilon_{\ell})\varepsilon [u] : \mathbb{C} : \varepsilon
[u]\,\mathrm{d}\Omega + \int_{\Omega}\left(\frac{\omega_{c}
(\phi)}{\ell}+G_{c}\ell\nabla \phi \nabla \phi\right)\,
\mathrm{d}\Omega 
\end{equation}}\unskip
where:

$D(\phi)=(1-\phi)^{2}$: represents the stiffness degradation function.
It is a continuous, strictly decreasing function, which reaches zero
only for $\phi=1$, corresponding to cracking zones.

$\omega_{c}(\phi)=({\phi^{2}}/{4})G_{c}$ represents the dissipation
function. It is a strictly increasing function that takes a zero value
when $\phi=0$.

$A(\phi)=E_{0}D(\phi)=E_{\mathrm{s}}(1-\phi)^{2}$ represents the elasticity
tensor affected by damage.

$\varepsilon_{\ell}$ is a numerical regularization parameter that must
be significantly smaller than~$\ell$.

The phase field model has received various interpretations. Different
families of damage models have been proposed, each based on specific
choices for the stiffness degradation function $D(\phi)$ and the local
dissipation function $\omega_{c}(\phi)$. For further information, the
interested reader may refer to the work mentioned
in~\cite{48,54,55,72,93}.\vspace*{-.25pc}

\subsection{The minimization variational framework for the phase field
approach} \label{sec4.2}

The phase-field model presented above enables us to elucidate crack
formation and its spatio-temporal trajectories in the material, without
the need for prior specification. This variational framework guarantees
its applicability in an elastic medium subjected to
displacement-controlled loading, making the quasi-static evolution
problem one that seeks to minimize the total energy characterized by
the following function:
{\begin{equation} \label{eq29}
\min \Pi_{\mathrm{t}}(u,\phi)\Rightarrow \Pi_{\mathrm{t}}^{\prime(u,\phi)}=0
\end{equation}}\unskip
were
{$$
(\cdot)^{\prime(a,b)}=\left(\frac{\partial (\cdot)}{\partial a\,
\partial b}\right).
$$}\unskip
Due to the presence of the term $((1-\phi)^{2}+\varepsilon_{\ell})
w(\varepsilon)$, the functional $\Pi_{\mathrm{t}}$ is not convex with respect to
the pair $(u,\phi)$. However, it is observed that, in general,
$\Pi_{\mathrm{t}}$ is convex with respect to each variable separately. In other
words, when $u$ is fixed, the functional becomes convex with respect
to~$\phi$, and when $\phi$ is fixed, it is convex with respect to~$u$.
Consequently, by fixing one of the variables, finding the minimum with
respect to the other becomes relatively straightforward, as each
problem becomes linear after derivation with respect to the
corresponding variable. To this end, an alternative has been presented
which consists in decomposing the functional (\ref{eq29}) into two
variational problems which are equivalent to it. Indeed, the:
{\begin{eqnarray}
&  \displaystyle 
\label{eq30}
\Pi_{\mathrm{t}}^{\prime(u, \phi)(\chi, 0)}
\left\{\begin{array}{@{}c@{}}
\text{find } u;\; u \epsilon A_{\delta} \text{ such that},\,
\forall \nu \epsilon \mathcal{A}_{0} :\vspace*{3pt}\\
\displaystyle\int_{\Omega}((1-\phi)^{2}+\varepsilon_{\ell})
\varepsilon (u)\varepsilon(u)\,\mathrm{d}\Omega=0
\end{array}\right.\Seqnsplit
& \label{eq31} \Pi_{\mathrm{t}}^{\prime(u, \phi)(0, \Theta)}
\left\{\begin{array}{@{}c@{}}
\text{find }\phi;\; \phi \epsilon H^{1}(\Omega)\text{ such that},\,
\forall \Theta \epsilon H^{1}(\Omega) :\vspace*{3pt}\\
\displaystyle\int_{\Omega}\left(-(1-\phi)\Theta\;
\varepsilon_{u}[\mathbb{C}]\varepsilon_{u}+\frac{G_{c}}{2\ell}
\phi \Theta +2G_{c}\ell\nabla \phi \nabla \Theta\right)\,\mathrm{d}
\Omega =0.
\end{array}\right.
\end{eqnarray}}\unskip
The test functions $\chi$ and $\Theta $ are defined to cancel on the
boundaries of $\mathrm{d}\Omega$, where the Dirichlet conditions are
applied.

By solving Equation~(\ref{eq30}), we obtain the displacement field $u$
that minimizes energy for a fixed~$\phi$. This solution is then
inserted into Equation~(\ref{eq31}) to determine the damage $\phi$ that
minimizes energy for a fixed~$u$. This process is iterated until the
two variables $u$ and $\phi$ converge.\vspace*{-.25pc}

\subsection{Finite element method implementation} \label{sec4.3}

In this section, we will attempt to discretize the variational problems
set out in Equations~(\ref{eq30}), (\ref{eq31}). We have chosen a
triangular element with three nodes located at the vertices and adopted
a polynomial interpolation basis of degree~1, where $N_{\mathrm{s}}$
and $N_{\mathrm{t}}$ are respectively the number of nodes and the number of
triangles in the mesh.

The development of Equation~(\ref{eq30}) by applying the finite element
method (FEM) discretization gives:
{\begin{equation} \label{eq32}
\varepsilon_{\ell}\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\int_{T_{\mathrm{e}}}[B^{u}][\mathbb{C}]
[B^{\chi}]\,\overline{u}_{i}\,\overline{\chi}_{j}\,\mathrm{d}T_{\mathrm{e}}+
\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\int_{T_{\mathrm{e}}}[B^{u}][\mathbb{C}][B^{\chi}]\,
\overline{u}_{i}\,\overline{\chi}_{j}\,\mathrm{d}T_{\mathrm{e}}
\int_{T_{\mathrm{e}}}(1-\phi)^{2}\,\mathrm{d}T_{\mathrm{e}}=0.
\end{equation}}\unskip
The application of Dirichlet-type boundary conditions on the
displacement field cancels out the arbitrary auxiliary test nodal
solution $\overline{\chi}_{j}$, Equation~(\ref{eq32}) becomes:
{\begin{equation} \label{eq33}
\varepsilon_{\ell}\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\int_{T_{\mathrm{e}}}[B^{u}]
[\mathbb{C}][B^{\chi}]\,\overline{u}_{i}\,\mathrm{d}T_{\mathrm{e}}+
\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\int_{T_{\mathrm{e}}}[B^{u}][\mathbb{C}][B^{\chi}]\,
\overline{u}_{i}\,\mathrm{d}T_{\mathrm{e}}\int_{T_{\mathrm{e}}}(1-\phi)^{2}\,
\mathrm{d}T_{\mathrm{e}}=0.
\end{equation}}\unskip
The Equation~(\ref{eq33}) gives the following linear system:
{\begin{equation} \label{eq34}
(\varepsilon_{\ell}[R_{ij}^{u}]+[R_{i,j}^{u,\phi}])
\{\overline{u}_{i}\}=0
\end{equation}}\unskip
were
{\begin{equation} \label{eq35}
\left\{\begin{array}{@{}c@{}}
\displaystyle [R_{ij}^{u}]=\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\int_{T_{\mathrm{e}}}
[B^{u}]^{\mathrm{T}}[\mathbb{C}][B^{u}]\,\mathrm{d}T_{\mathrm{e}}
\vspace*{3pt}\\
\displaystyle {}[R_{i,j}^{u,\phi}]=\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\int_{T_{\mathrm{e}}}
[B^{u}]^{\mathrm{T}}[\mathbb{C}][B^{u}]D(\phi)\,\mathrm{d}T_{\mathrm{e}}
\end{array}\right.
\end{equation}}\unskip
with
{\begin{equation} \label{eq36}
[D_{r,s}]=\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\int_{T_{\mathrm{e}}}(1-\phi)^{2}\,\mathrm{d}
T_{\mathrm{e}}=\frac{\text{Aire}(T_{\mathrm{e}})}{6P}\sum_{r,s=i,j,l}(1-\phi)_{r}
(1-\phi)_{\mathrm{s}},\quad\text{such as }P=\begin{cases}
1, & \mathrm{si}\,r=s\\
2, & \mathrm{si}\,r\neq s
\end{cases}
\end{equation}}\unskip
where $i$, $j$, $l$ are the indices corresponding to the node number
per element.
{\begin{equation} \label{eq37}
[D_{r,s}]=\frac{\text{Aire}(T_{\mathrm{e}})}{6}[(1-\phi_1)^{2\!}+(1-\phi_2)^{2\!}
+(1-\phi_{3})^{2\!}+(1-\phi_{1})(1-\phi_{2})+(1-\phi_{1})(1-\phi_{3})+
(1-\phi_{2})(1-\phi_{3})]
\end{equation}}\unskip
were
{$$
\text{Aire}(T_{\mathrm{e}})=\tfrac{1}{2}[(S_{y}^{l}-S_{y}^{i})x+
(S_{x}^{j}-S_{x}^{i})y+(S_{x}^{i}-S_{x}^{j})-(S_{y}^{j}-S_{y}^{i})]. 
$$}\unskip
Similarly, the development of Equation~(\ref{eq31}) gives:
\vspace*{-2pt}
{\advance\jot by -3pt\begin{eqnarray}
&  \displaystyle
\label{eq38}
\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}} \int_{T_{\mathrm{e}}}-(1-\phi)\Theta
\varepsilon[\mathbb{C}]\varepsilon\,\mathrm{d}T_{\mathrm{e}}+\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}
\int_{T_{\mathrm{e}}}\frac{G_{c}}{2\ell} \phi \Theta\, \mathrm{d}T_{\mathrm{e}}+
\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\int_{T_{\mathrm{e}}}2G_{c}\ell\nabla \phi \nabla\Theta\,
\mathrm{d}T_{\mathrm{e}}=0\Seqnsplit
&  \displaystyle
\label{eq39}
\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}} \int_{T_{\mathrm{e}}}\phi \Theta \varepsilon
[\mathbb{C}]\varepsilon\, \mathrm{d}T_{\mathrm{e}}+\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}
\int_{T_{\mathrm{e}}}\frac{G_{c}}{2\ell} \phi \Theta\, \mathrm{d}T_{\mathrm{e}}+
\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\int_{T_{\mathrm{e}}}2G_{c}\ell\nabla \phi \nabla \Theta\,
\mathrm{d}T_{\mathrm{e}} = \int_{T_{\mathrm{e}}} \sum_{\mathrm{e}=0}^{N_{\mathrm{t}}}\Theta \varepsilon
[\mathbb{C}]\varepsilon\, \mathrm{d}T_{\mathrm{e}}\Seqnsplit
& \displaystyle 
\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}} \int_{T_{\mathrm{e}}}[N^{\phi}][N^{\Theta}]
\overline{\phi}_{i}\,\overline{\Theta}_{j}\,\varepsilon [\mathbb{C}]
\varepsilon\, \mathrm{d}T_{\mathrm{e}}+\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\frac{G_{c}}{2\ell}
[N^{\phi}][N^{\Theta}]\overline{\phi}_{i}\,\overline{\Theta}_{j}\,
\mathrm{d}T_{\mathrm{e}}+\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\int_{T_{\mathrm{e}}}2G_{c}\ell[B^{\phi}]
[B^{\Theta}]\overline{\phi}_{i}\,\overline{\Theta}_{j}\,\mathrm{d}
T_{\mathrm{e}}
\nonumber\\ & \displaystyle 
= \sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\int_{T_{\mathrm{e}}}[N^{\Theta}]
\overline{\Theta}_{j}\varepsilon [\mathbb{C}]\varepsilon\,
\mathrm{d}T_{\mathrm{e}}. \label{eq40}
\end{eqnarray}}\unskip
The application of Dirichlet-type boundary conditions on the damage
field cancels out the arbitrary auxiliary test nodal solution
$\overline{\Theta}_{j}$, Equation~(\ref{eq40}) becomes:
{\begin{eqnarray}
& \displaystyle
\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}} \int_{T_{\mathrm{e}}}[N^{\phi}][N^{\Theta}]
\overline{\phi}_{i}\varepsilon [\mathbb{C}]\varepsilon\,\mathrm{d}
T_{\mathrm{e}}+\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\frac{G_{c}}{2\ell}[N^{\phi}][N^{\Theta}]
\overline{\phi}_{i}\,\mathrm{d}T_{\mathrm{e}}+\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}
\int_{T_{\mathrm{e}}}2G_{c}\ell[B^{\phi}][B^{\Theta}]\overline{\phi}_{i}\,
\mathrm{d}T_{\mathrm{e}}
\nonumber\\& \displaystyle 
=\sum_{\mathrm{e}=1}^{N_{\mathrm{t}}}\int_{T_{\mathrm{e}}}[N^{\Theta}]
\varepsilon [\mathbb{C}]\varepsilon\, \mathrm{d}T_{\mathrm{e}}. \label{eq41}
\end{eqnarray}}\unskip
The Equation~(\ref{eq41}) gives the following linear system:
{\begin{equation} \label{eq42}
(2G_{c}\ell[R_{ij}^{\phi}]+[M_{ij}])\{\overline{\phi}_{i}\}=
\{F_{j}\}
\end{equation}}\unskip
were
{\begin{equation} \label{eq43}
\left\{\begin{array}{@{}c@{}}
\displaystyle [R_{ij}^{\phi}]=\int_{T_{\mathrm{e}}}[B^{\Theta}][B^{\phi}]\,
\mathrm{d}T_{\mathrm{e}}\vspace*{3pt}\\
\displaystyle {}[M_{ij}]=\left(\varepsilon[\mathbb{C}]\varepsilon +
\frac{G_{c}}{2\ell}\right)\int_{T_{\mathrm{e}}}[N^{\phi}]^{\mathrm{T}}
[N^{\phi}]\,\mathrm{d}T_{\mathrm{e}}\vspace*{3pt}\\
\displaystyle \{F_{j}\}=\varepsilon [\mathbb{C}]\varepsilon \int_{T_{k}}
[N^{\Theta}]\,\mathrm{d}T_{\mathrm{e}}.
\end{array}\right.
\end{equation}}\unskip
The chronological sequence of computational steps associated with all
these numerical developments can be summarized by following the
framework of the alternating minimization scenario detailed above in
the context of the phase field model applied via the finite element
method (FEM) according to the following algorithmic structure:

\beqimg\begin{eqimg*}
\begin{algorithm}
\caption{\label{alg1}Solution algorithm alternate minimization for
phase field models} 
\small
\textbf{1:}\quad Set the parameters of this model: $\ell$, the internal
length of the damage model; $\varepsilon_{\ell}$, a rest of rigidity
(parameter that ensures the existence and uniqueness of the problem in
motion); $h$, the size of the mesh; $\delta$, the convergence tolerance
of the damage model; $\mathrm{d}u$ discretization of the imposed
displacement;

\textbf{2:}\quad\, \textbf{for} the moment: $t_{I}$ : \textbf{do}\\
Let $k = 0$ and $\phi^{0} := \phi_{t_{i-1}}$ ($k$ is the discretization
of the displacement).

\begin{inftab}
\begin{tabular}{@{}ll}
\xtbody
\textbf{3:} & \textbf{As long as} the convergence on $\mathrm{d}$ is
not achieved \textbf{do}\\
& Conditions at the limits of the Dirichlet type at the moment
$\mathrm{I}$ applied for the displacement\\
& \textbf{for} a damage field set in $\phi^{k}$, by using $(34)$
calculate the displacement $u^{k+1}$\\
& Conditions at the limits of the Dirichlet type at the moment $j$
applied for the damage\\
& \textbf{for} a displacement field set in $u^{k+1}$, by using $(42)$
calculate the damage $\phi^{k+1}$\\
& under reserve of the irreversibility condition $0\leq\phi_{t_{i-1}}
\leq \phi^{k+1}\leq 1$\\
& $k := k+1$\\
& \textbf{end as long as}\\
& $u_{t_{i}}:=u^{k+1}$;\, $\phi_{t_{i}} := \phi^{k+1}$\\
\textbf{4:} & \textbf{Post-processing calculation}: Evaluation of
energies ($\Pi_{\mathrm{ei}}$, $\Pi_{\mathrm{si}}$), reactions
$F_{\mathrm{xi}}\;R_{\mathrm{yi}}, \ldots$ etc. $i := i+1$\\
& \textbf{end~for}\\
& ${u}_{{i}}={u}_{{i}-{1}}+ \mathrm{d}{u}_{{i}}$ 
\end{tabular}
\end{inftab}
\end{algorithm}
\eeqimg\end{eqimg*}

\section{Numerical results} \label{sec5}

This section is devoted to validating the results obtained with our
numerical contribution. To do this, we will proceed in two main stages.
Firstly, we will compare our results with analytical solutions using a
reference example in simple tension, which will enable us to assess the
accuracy and precision of the model employed. At the same time, we will
carry out a convergence analysis to ensure the stability and
reliability of the calculation code we have developed.

Secondly, we will implement numerical modeling of the experimental
tensile tests described in Section~\ref{sec2.5}. The aim is to verify
experimentally the ability of our calculation code to reproduce the
actual behavior of concrete mortar. To this end, we will consider two
types of domains: the logical domain, which represents the actual
configuration of the problem, and the physical domain, which includes a
boundary condition on the damage field aimed at preventing excessive
localization of damage at imposed boundaries. A comparative analysis
between these two domains will be carried out in order to assess their
respective impact on the quality and accuracy of results, and to
determine which one proves more appropriate on the basis of
experimental results.

\subsection{Problem of a rectangular plate under simple tension
loading} \label{sec5.1}

The case studied concerns a single-axis tensile test applied to a
rectangular plate, illustrated in Figure~\ref{fig9}, corresponding to
the domain $\Omega=(0;L)*(0;H)$ with $L=91$~mm and $H=168$~mm. In
addition, a zero-displacement condition $u(y=0)=0$ is imposed at the
lower end, while a displacement-controlled imposed loading is applied
at the upper end $u(y=L)$. The other parts of the domain are free.
Furthermore, damage is assumed to be zero at the edges where
displacement-controlled loading is applied.

\begin{figure}
\includegraphics{fig09}
\caption{\label{fig9}Problem of a rectangular plate under tension
loading.}
\end{figure}

The material used is a concrete mortar with the same mechanical
parameters found experimentally; Young's modulus $E_{\mathrm{s}}=0.8$~GPa
surface energy density $G_{c}=72.48$~N/m, Poisson's ratio $\nu=0.21$.
The numerical regularization parameter is $\varepsilon_{\ell}=10^{-5}$
and $\ell=5\times 10^{-3}$~m and a loading controlled by imposed
displacement $u=N\cdot\check{u}=4\times 10^{-4}$~mm with
$\check{u}=4\times 10^{-7}$~mm and $N=1000$ imposed displacement steps.

The mesh consists of 32~122 triangular elements and 16~239 nodes,
with elements of length $h_{\mathrm{e}}=0.001$.

First, we discuss the analytical results of energy variations, such as
the elastic energy, surface energy and total energy, as a function of
the loadings~$u$. We then present the critical load that causes damage
initiation.

By using the analytical equations~\cite{94,95}, the energies are
presented by:
\begin{enumerate}[(1)]
\item Elastic energy before rupture:
{\begin{equation} \label{eq44}
\Pi_{\mathrm{e}}=\int_{0}^{L}\frac{1}{2}E_{\mathrm{s}}\frac{\mathbf{u}_{c}^{2}}
{H^{2}}\,\mathrm{d}x\,\mathrm{d}y=\frac{1}{2}\frac{E_{0}
\mathbf{u}_{c}^{2}}{H}L=6.5884\text{ Joule}.
\end{equation}}\unskip
\item Surface energy:
{\begin{equation} \label{eq45}
\Pi_{\mathrm{s}}=\int_{0}^{H} G_{c}\,\mathrm{d}x = G_{c}H=6.5884\text{ Joule}.
\end{equation}}\unskip
\end{enumerate}

We use energy balance to determine the critical load that causes the
structure to fail. This is expressed by the following relationship:
{\begin{equation} \label{eq46}
u_{c}^{th}=\sqrt{\frac{2H\,G_{c}}{E_{0}}}1.743\times 10^{-4}~\text{mm}.
\end{equation}}\unskip

From Figure~\ref{fig10}(a), we can see that as long as the loading
remains below~1, the structure behaves in a purely elastic manner.
However, as soon as this limit is exceeded, a progressive degradation
of the material begins, evolving in proportion to the increase in
loading. This degradation continues gradually until a critical value of
2.387 is reached, beyond which the damage becomes significant and
concentrates in the center of the structure (Figure~\ref{fig10}(b)),
leading to complete fracture. The evolution of energies confirms this
behavior. Elastic energy increases progressively until it reaches a
maximum of 0.8979. Once the critical loading value is reached, it
decreases sharply until it is cancelled out, and then remains stable at
zero. On the other hand, surface energy continues to rise, reaching a
maximum of 1.132, where it stabilizes. This represents the total
rupture of the structure, confirming the complete separation of the
material.\looseness=-1

\begin{figure}
\includegraphics{fig10}
\caption{\label{fig10}(a)~Evolution of energies as a function of
imposed displacement; (b)~damage distribution.}
\end{figure}

Table~\ref{tab7} compares the analytical solutions with the
numerical solutions in terms of the energies-imposed displacement
response.

\begin{figure}
\includegraphics{fig11}
\caption{\label{fig11}(a)~Relative error as a function of number of
elements; (b)~log(error) as a function of number of elements.}
\end{figure}

\begin{table}%tab7
\caption{\label{tab7}Comparison of the numerical results with the
analytical results}
\begin{tabular}{ccc}
\thead
& Analytical results & Numerical results\\
\endthead
Critical displacement $u_{c}$ (mm) & $1.743\times 10^{-4}$\raisebox{.85pc}{} & $2.407\times 10^{-4}$\\
Surface energy $\Pi_{\mathrm{s}}$ (Joule) & 6.5884 & 7.379\\
Elastic energy $\Pi_{\mathrm{e}}$ (Joule) & 6.5884 & 6.043\\
Total energy $\Pi_{\mathrm{t}}$ (Joule) & 6.5884 & 7.491
\botline
\end{tabular}
\end{table}

We find that the critical value of the imposed displacement,
$2.407\times 10^{-4}$~mm, has an estimated error of 38\% compared with
the analytical value of $1.743\times 10^{-4}$~mm. As for surface
energy, with a value of 7.379~Joule, the error is 11.99\% compared with
the analytical value of 6.5884~Joule. This latter error is considerably
smaller than that observed for critical loading. For the elastic
energy, whose value is 6.043~Joule, the error is 8.278\% compared with
the analytical value of 6.5884~Joule. Finally, for the total energy,
whose value is 7.491~Joule, the error is 13.699\% compared with the
analytical value of 6.5884~Joule.

From the results obtained, we note a relative correspondence between
analytical and numerical values, confirming the consistency of the
model used.

\subsubsection{Convergence study} 

A convergence analysis is proposed on the basis of the analytical
solution, in order to demonstrate the accuracy and stability of the
present approach, as well as the efficiency of our developed
calculation code. To this end, the relative error in energy norms is
evaluated, considering the analytical solution as the exact solution:
{\begin{equation} \label{eq47}
\mathrm{Er\%}=\frac{\|\Pi_{\mathrm{t}}^{h_{\mathrm{e}}}-\Pi_{\mathrm{t}}^{\mathrm{EX}}\|}
{\|\Pi_{\mathrm{t}}^{\mathrm{EX}}\|}.
\end{equation}}\unskip
The numerical solution $\Pi_{\mathrm{t}}^{h_{\mathrm{e}}}$ is obtained
by running a series of simulations on the current example of a
rectangular plate under tension, by varying the mesh
size~$h_{\mathrm{e}}$.

The evolution of the error (Er\%) as a function of the number of mesh
elements is shown in Figure~\ref{fig11}(a). It is observed that the
error decreases significantly with mesh refinement up to a certain
threshold (3.95 ${\times}$ 104 elements), after which it stabilizes,
confirming the stability of our developed calculation code.

The rate of convergence is shown in Figure~\ref{fig11}(b), where the
value of the convergence rate is 0.94 for all configurations, close
to~1. This represents near-ideal convergence for a mesh using linear
elements, and testifies to the efficiency of our developed calculation
code.

\subsection{Problem of a rectangular plate with a crack and a void}
\label{sec5.2}

Through our calculation code, we aim to numerically reproduce the
tensile test carried out experimentally in (Section~\ref{sec3}), using
the phase field model implemented by the finite element method (FEM) in
(Section~\ref{sec4}).

The material used is that experimentally employed in
(Section~\ref{sec2}), i.e.\ a concrete mortar with the following
mechanical characteristics: Young's modulus $E_{\mathrm{s}}=0.8$~GPa, Poisson's
ratio $\nu= 0.21$, fracture energy $G_{c}=72.48$~N/m.

The choice of the internal length $\ell=0.5\times 10^{-3}$~m was guided
by the granulometric analysis of the sand (Section~\ref{sec2.1},
Figure~\ref{fig1}), which showed that most grains fall within the
0.1--1~mm range. This value ensures consistency between the mortar's
microstructural scale and the regularization parameter of the
phase-field model, while guaranteeing both physical representativity
and satisfactory numerical stability.

The numerical regularization parameter $\varepsilon_{\ell}=10^{-5}$, a
loading $u=N\cdot\check{u}=1.5\times 10^{-3}$~mm with $N=1000$ step
displacement imposed.

To this end, two simulations have been carried out, one on the logical
domain and the other on the physical domain. The definition of each
domain is described below.

\subsubsection{Logical domain}

Is the reference domain, the configuration is similar to the
experimental test carried out in (Section~\ref{sec2.5}, see
Figure~\ref{fig6}), is a modification of the simple tensile test
designed to ensure a homogeneous deformation field in the concrete
mortar. The idea is to glue copper tubes to the holes where the
displacement-controlled loading is applied. Due to copper's
sufficiently high stiffness, this will prevent damage from localizing
at the edges of the holes where the loading is applied. The mechanical
parameters of the copper used are: Young's modulus $E_{\mathrm{s}2}=128$~GPa,
Poisson's ratio $\nu_{2}=0.33$, fracture energy $G_{c2}=50.096\times
10^{3}$~N/m (see Figure~\ref{fig12}(a)).

\begin{figure}
\includegraphics{fig12}
\vspace*{-.35pc}
\caption{\label{fig12}The geometry of the specimen studied (a)~logical
domain; (b)~physical domain.}
\vspace*{-.15pc}
\end{figure}

\subsubsection{Physical domain}

Instead of bonding copper tubes, we apply zero damage $(\phi=0)$ to the
edges of the holes where the loading is applied. This condition has the
same effect as the previous one, since it avoids localizing damage on
these edges (see Figure~\ref{fig12}(b)).

To ensure reliable numerical reproducibility of the experimental tests
in both configurations investigated (logical domain and physical
domain), it is essential to first carry out a mesh convergence study.
This step makes it possible to determine the optimal number of elements
that guarantees numerical stability, accuracy of the results, and
compliance with the regularization conditions imposed by the
phase-field model, while minimizing computational cost. To this end,
several discretizations of the physical domain were tested by varying
the total number of finite elements. The objective is to identify a
mesh sufficiently refined to accurately capture crack propagation,
while avoiding an unnecessarily high computational cost.

Figure~\ref{fig13} illustrates, on the one hand, the discretized
structures corresponding to the different meshes considered, and on the
other hand, the evolution of the associated damage field for each mesh,
shown in the upper right corner of each illustration. This dual
representation enables simultaneous visualization of the influence of
mesh refinement on both element distribution and the accuracy of damage
modeling.

\begin{figure}
\includegraphics{fig13}
\caption{\label{fig13}Mesh convergence study: (a)~17~858 elements
(TRI3); (b)~31~112 elements (TRI3); (c)~47~917 elements (TRI3);
(d)~680~025 elements (TRI3); (e)~728~129 elements (TRI3);
(f)~764~438 elements (TRI3).}
\end{figure}

Figure~\ref{fig14} shows the load--displacement curves obtained for
each of the considered meshes. This comparison highlights the influence
of mesh density on the simulated mechanical response, particularly in
terms of initial stiffness, peak load, and post-peak behavior.

\begin{figure}
\includegraphics{fig14}
\caption{\label{fig14}Comparison of the effect of the number of
elements on the load--displacement behavior.}
{\vspace*{.25pc}}
\end{figure}

The convergence study highlights the decisive influence of mesh
refinement on the accuracy of the numerical results. With coarse
meshes, crack trajectories appear diffuse and may deviate slightly from
the expected path, while the load--displacement curves reveal an
overestimation of the initial stiffness and an inaccurate prediction of
the peak load. Moreover, it should be noted that coarse meshes do not
satisfy the condition $h_{\mathrm{e}}\ll \ell$. Specifically, for 17~858
elements, the characteristic element size is $h_{\mathrm{e}}=1.49\times
10^{-3}~\mathrm{m}>\ell$; for 31~112 elements, it is 
\mbox{$h_{\mathrm{e}}=1.12\times
10^{-3}~\mathrm{m}>\ell$;} and for 47~917 elements, it remains
$h_{\mathrm{e}}=0.9\times 10^{-3}~\mathrm{m}>\ell$.

From a sufficiently refined mesh (680~025 elements) and beyond, crack
propagation becomes stable and well-defined, with no significant
variation across different refinement levels. The load--displacement
curves then overlap satisfactorily, indicating robust numerical
convergence. Accordingly, the adopted mesh (680~025 TRI3 elements) was
identified as optimal, ensuring both the reliability and
reproducibility of the results.

This choice provides remarkable stability while satisfying the
condition $h_{\mathrm{e}}=0.22\times 10^{-3}~\mathrm{m}<\ell$.

Figure~\ref{fig15} presents a comparison between the load--displacement
curve obtained numerically using the phase-field model with the adopted
mesh of 680~025 elements in both the physical and logical domains, and
that derived from the experimental results.

\begin{figure}
\includegraphics{fig15}
\caption{\label{fig15}Comparison of the experimental results with the
numerical results.}
\end{figure}

The load--displacement response illustrated in Figure~\ref{fig15} shows
identical behavior across all curves prior to fracture. The structure
first exhibits a linear phase up to the first load drop. Numerically,
an elastic rebound corresponding to a load increase (zone~A) is then
observed, followed by the re-initiation of the crack on the right side
of the void and a subsequent abrupt drop. This phenomenon, however,
only appears in the numerical results. Experimentally, the~load
decreases continuously as displacement increases until reaching zero,
which corresponds to the complete failure of the specimen.

The divergence observed in zone A can be explained by several
mechanisms related to the heterogeneous and relatively low-stiffness
nature of mortar. On the one hand, the actual material develops an
extended fracture process zone (FPZ) ahead of the crack tip,
characterized by microcracking and paste--sand decohesion. Once the
crack initiates, the stress field quickly interacts with the large void
due to the short distance separating them. A progressive decohesion of
the interfaces surrounding the void occurs even before the main crack
reaches it. Consequently, a significant portion of the mechanical
energy is dissipated, such that when the crack eventually reaches the
void, it is already ``prepared'' and can be crossed at a reduced
energetic cost. The global response thus manifests as a continuous
softening, without the appearance of a distinct second peak.

On the other hand, the numerical model assumes an ideal geometry and
treats the mortar as a homogeneous and isotropic material, with
dissipation concentrated essentially along the main crack. It does not
capture the diffuse decohesion around the void nor the local
heterogeneities (porosity, weak paste--sand interfaces, notch
irregularities). Moreover, the displacement is imposed ideally in the
simulation, allowing the model to follow unstable branches of the
energy path. By contrast, even under quasi-static conditions, the
specimen and testing setup dissipate this instability, leading to a
more gradual and continuous load decrease.

Finally, these results confirm that the phase-field model is capable of
capturing highly brittle fracture phenomena, including elastic rebound
behavior, even though certain features observed experimentally remain
attenuated in the simulation.

The damage evolution in the structure is shown in Figure~\ref{fig16}.

\begin{figure}
{\vspace*{-.25pc}}
\includegraphics{fig16}
{\vspace*{-.25pc}}
\caption{\label{fig16}The evolution of the damage field in the
structure (a)~logical domain; (b)~physical domain; (c)~experimental
result.}
{\vspace*{-.25pc}}
\end{figure}

We can see that the crack starts in a straight line from the notch. It
then continues along a curved path before joining the hole. At the
imposed displacement $u_{i}=0.73$~mm, a secondary crack is initiated
from the hole and propagates in a straight line to the right-hand edge
of the structure (total fracture of the specimen). We can see that the
results of the numerical simulations in terms of cracking trajectories
are in line with the results observed experimentally.

From a comparative point of view, there is a high degree of consistency
between the results obtained in the logical and physical domains.

In order to better understand the effect of reinforcement on the holes
where the loading is applied, whether in the logical domain through the
bonding of copper tubes to these holes, or in the physical domain
through the application of zero damage $(\phi=0)$ to the edges of these
holes, a simulation was carried out where we added no reinforcement to
these holes, and applied no conditions on the damage field at the edge
of the holes where the loading is applied, i.e.\ the damage field is
free throughout the structure.

The numerical and experimental results obtained are shown in
Figure~\ref{fig17}.

\begin{figure}
{\vspace*{-.4pc}}
\includegraphics{fig17}
{\vspace*{-.4pc}}
\caption{\label{fig17}The evolution of the damage field in the
structure (a)~numerical result; (b)~experimental result.}
{\vspace*{-.4pc}}
\end{figure}

The results obtained in Figure~\ref{fig17}(a) and (b) show qualitative
agreement between numerical and experimental results. The structure was
damaged on the holes where the loading was applied, confirming the
validity and reliability of the experimental protocol adopted in the
previous experiments in (Section~\ref{sec2}) on the one hand, and the
effectiveness of the numerical adaptation adopted to solve the problem
of localizing damage at the edges of the loading zones in concrete
mortar on the other.

{\vspace*{-.25pc}}

\section{Conclusion} \label{sec6}

{\vspace*{-.25pc}}

The present work presents an experimental study and numerical modeling
by the phase field model of the behavior of concrete mortar that
contains predefined discontinuities, several points on have been
addressed:
\begin{itemize}
\item The material parameters associated with the studied mortar were
characterized in an integrated manner. The dynamic modulus of
elasticity $E_{\mathrm{d}}$ and Poisson's ratio $\vartheta$ were monitored using
ultrasonic testing at the ages of 7, 14, and 28~days, allowing repeated
measurements on the same specimens. A three-point bending calibration
enabled the establishment of a linear empirical relationship between
the static and dynamic moduli ($E_{\mathrm{s}}=a\times E_{\mathrm{d}}+b$), specifically
adapted to the investigated mortar. The fracture energy $G_{c}$ was
estimated from the load--displacement curves obtained in direct tensile
tests and subsequently employed for numerical calibration.
\item Quasi-static direct tensile tests were conducted on prismatic
specimens containing a predefined crack and hole, using a
servo-hydraulic testing machine. Three categories of specimens (7, 14,
and 28~days) were tested to analyze the influence of age on the
structural response.
\item The phase-field method was implemented numerically within the
framework of brittle fracture by the finite element method (FEM)
through a MATLAB{\textregistered}-based computational code. The
variational formulation was decomposed into two coupled equations
(displacement and damage), solved using an explicit alternate
minimization strategy. This implementation ensured enhanced numerical
stability and reproducibility. The regularization parameter $\ell$ was
physically justified by the sand granulometry, thereby establishing a
direct link between the microstructure and the numerical model.
\item The computational code was validated on a reference uniaxial
tensile case and then applied to the experimental problem through the
introduction of an original double-domain strategy: the logical domain,
in which the complete experimental configuration was reproduced, and
the physical domain, where the mechanical effect of the tubes was
replaced by a simplified boundary condition $(\phi=0)$. This approach
allowed numerical simplification while preserving physical
representativeness. The phase-field method thus enabled the prediction
of both crack propagation paths and load--displacement curves. A
comparison was carried out between experimental and numerical results.
\end{itemize}

The following conclusions can be drawn from the results obtained:
\begin{itemize}
\item The mechanical properties of the mortar (elastic modulus,
Poisson's ratio, stiffness) evolve with age, as evidenced by the
characterization tests. This evolution is directly reflected in the
direct tensile tests, where both strength and stiffness increase
consistently with the characterization results.
\item The experimental results highlight a progressive but non-linear
increase in stiffness with age: a modest rise in the critical load
between 7 and 14~days, followed by a pronounced increase between 14 and
28~days. These observations reflect a progressive densification of the
microstructure.
\item The phase-field damage model accurately reproduces the initiation
and propagation of cracks in cement-based mortar, in overall agreement
with the experimental observations. The numerical load--displacement
curve, however, exhibits a rebound associated with the local resistance
at the vicinity of the hole a phenomenon not observed experimentally,
which can be attributed to the quasi-heterogeneous and highly brittle
behavior of the mortar, characterized by paste--sand decohesion leading
to a smoother global response. Numerically, dissipation nevertheless
remains concentrated along the main crack.
\item Instead of creating a logical domain, it would be sufficient to
apply zero damage to the edges of the structure where displacement
loading applies. This condition has the same effect as the previous
one, since it avoids the localization (note that if the structure
wanted to damage itself on these edges, this condition would not
prevent it, since it would happen anyway) of the damage on these edges.
The latter is easier to implement and is equivalent to the logical
domain, which can become a difficult task for structures with more
complex geometries.
\item Due to the intrinsic nature of concrete as a heterogeneous
composite material, characterized by its non-uniform distribution of
various constituents and the complex interfaces that exist between
these components, it becomes clear that these aspects play a crucial
role in provoking not only the location of cracks, but also the
subsequent paths through which these cracks propagate, as well as the
ultimate mode of failure that may ensue.
\end{itemize}

Finally, this work opens up several avenues for future research:
\begin{itemize}
\item From an experimental perspective, an interesting direction would
be the integration of local crack measurements (e.g., CMOD or DIC) in
future experimental campaigns, in order to complement the global
load--displacement measurements and achieve a more refined
characterization of the direct tensile behavior.
\item From a numerical perspective, a promising extension would be the
implementation of an advanced arc-length type procedure, to better
capture the snap-back phenomenon and analyze the unstable response of
the material beyond the limitations of displacement control adopted in
the present study.
\end{itemize}

\section*{Acknowledgements}

The first author would like to express his sincere gratitude to
Professor Hachi Brahim Elkhalil for his valuable guidance on the
phase-field damage model implemented in this work. He also wishes to
thank Professor Zouari Bassem for his kind invitation to his laboratory
at the National School of Engineers of Sfax and for his insightful
advice on the experimental aspects. Finally, he extends his heartfelt
thanks to Professor Mohamed Haboussi for his invitation to the LSPM
Laboratory at Sorbonne Paris Nord University, France, and for his
valuable recommendations regarding the structuring of this article.

\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.

\back{}

\appendix{}

\gdef\appendixlabel{Appendix A}

\setcounter{table}{0}
\renewcommand\thetable{A\arabic{table}}

\section{Specimen preparation and detailed characterization of
properties by ultrasonic tests} \label{appA}
\subsection{The steps for preparation of specimen} \label{appA1}

The specimens were prepared in the following steps:
\begin{itemize}
\item We prepared a wooden mold with the dimensions of the test
specimen mentioned above.
\item We ensured the proper clamping of the mold and eliminated any
voids between the formwork walls using silicone. Then, we oiled the
mold to facilitate de-molding.
\item We mixed the sand and cement dry for 30~s, then added half of the
water and mixed for 60~s. In the third step, we added the remaining
water and mixed for another 90~s.
\item The initial notch ($a_{0}=14$~mm) was directly created during
casting by inserting a metallic plate of the corresponding dimensions
into a dedicated slot in the mold. After filling the mixture and the
beginning of the setting, the plate was carefully removed, and the
mortar surface was leveled. This procedure ensured a clean and
reproducible notch geometry, while avoiding any post-hardening cutting
operations that could induce undesirable prestresses.
\item We placed copper tubes on the mold so that they would stick to
the inside surface of the small holes in the specimen after de-molding.
The sufficiently high stiffness of the copper prevents the specimen
from fracturing in the holes where the displacement loading is~applied.
This test is similar to the experimental PIED test proposed by Bazant
et~al.~\cite{96} and taken up by Ramtani et~al.~\cite{97}.
\item We introduced the mixture into the formwork with vibration to
distribute the grains evenly and ensure a compact fill.
\item Once the mold has been filled, the surface is levelled using a
ruler.
\item The de-molding is carried out after 24~h.
\item The test specimens were conserved at ambient temperature and
relative moisture in the Laboratory of development in mechanic and
materials (LDMM).{\vspace*{-.25pc}}
\end{itemize}

\begin{table}[b!]%tabA1
\footnotesize\tabcolsep=4pt
\caption{\label{tabA1}Detailed results of ultrasonic velocities
${V}_{\mathrm{p}}$, ${V}_{\mathrm{s}}$, Poisson's ratio
${\vartheta}$, and the moduli ${E}_{\mathrm{d}}$ and
${E}_{\mathrm{s}}$ for specimens at 7~days}
\begin{tabular}{ccccccccc}
\thead
& \parbox[t]{1.4cm}{\centering Longitudinal wave $V_{\mathrm{p}}$ $(\mathrm{m/s})\times 10^{3}$} &
\parbox[t]{1.25cm}{\centering Shear wave~$V_{\mathrm{s}}$ $(\mathrm{m/s}) \times 10^{3}$} &
\parbox[t]{.9cm}{\centering Poison's ratio $\vartheta$} &
\parbox[t]{.95cm}{\centering Average Poison's ratio $\vartheta$} &
\parbox[t]{1.2cm}{\centering Dynamic modulus of elasticity $E_{\mathrm{d}}$ (GPa)\raisebox{-.43pc}{}} &
\parbox[t]{1.35cm}{\centering Average dynamic modulus of elasticity $E_{\mathrm{s}}$ (GPa)} &
\parbox[t]{1.35cm}{\centering Average static modulus of elasticity $E_{\mathrm{s}}$ (GPa)} &
\parbox[t]{1.35cm}{\centering Average static modulus of elasticity $E_{\mathrm{s}}$ (GPa)}\\
\endthead
U-SP-07-D-01 & 4.9256 & 2.9851 & 0.2097 &\morerows{8}{0.2062}& 
45.0308 &\morerows{8}{44.6653}& 0.4953388 &\morerows{8}{0.4913188}\\
U-SP-07-D-02 & 4.9057 & 2.9695 & 0.2108 && 44.6001 && 0.4906011 &\\
U-SP-07-D-03 & 4.8598 & 2.9617 & 0.2046 && 44.1372 && 0.4855092 &\\
U-SP-07-D-04 & 4.8925 & 2.9931 & 0.2010 && 44.9404 && 0.4943444 &\\
U-SP-07-D-05 & 48533 & 2.9773 & 0.1983 &  & 44.3705 & & 0.4880755 & \\
U-SP-07-D-06 & 4.9389 & 3.0009 & 0.2074 && 45.4202 && 0.4996222 &\\
U-SP-07-D-07 & 4.8404 & 2.9539 & 0.2033 && 43.8596 && 0.4824556 &\\
U-SP-07-D-08 & 4.9189 & 2.9773 & 0.2109 && 44.8375 && 0.4932125 &\\
U-SP-07-D-09 & 4.9123 & 2.9773 & 0.2097 && 44.7918 && 0.4927098 &
\botline
\end{tabular}
{\vspace*{.5pc}}
\end{table}

\subsection{Detailed results of ultrasonic tests} \label{appA2}

{\vspace*{.2pc}}

This section of Annex~A compiles the complete set of results obtained
from the ultrasonic characterization of mortar specimens at the ages of
7, 14, and 28~days. The longitudinal wave velocities
${V}_{\mathrm{p}}$ and shear wave velocities
${V}_{\mathrm{s}}$, the Poisson's ratio
${\vartheta}$, as well as the dynamic modulus
${E}_{\mathrm{d}}$ and the static modulus ${E}_{\mathrm{s}}$ were
determined for each specimen. Tables~\ref{tabA1}, \ref{tabA2},
and~\ref{tabA3} present the results corresponding to~7, 14, and
28~days, respectively.

\begin{table}%tabA2
\footnotesize\tabcolsep=4pt
\caption{\label{tabA2}Detailed results of ultrasonic velocities
${V}_{\mathrm{p}}$, ${V}_{\mathrm{s}}$, Poisson's ratio
${\vartheta}$, and the moduli ${E}_{\mathrm{d}}$ and
${E}_{\mathrm{s}}$ for specimens at 14~days}
\begin{tabular}{ccccccccc}
\thead
& \parbox[t]{1.4cm}{\centering Longitudinal wave $V_{\mathrm{p}}$ $(\mathrm{m/s})\times 10^{3}$} &
\parbox[t]{1.25cm}{\centering Shear wave~$V_{\mathrm{s}}$ \mbox{$(\mathrm{m/s}) \times 10^{3}$}} &
\parbox[t]{.9cm}{\centering Poison's ratio $\vartheta$} &
\parbox[t]{.95cm}{\centering Average Poison's ratio $\vartheta$} &
\parbox[t]{1.2cm}{\centering Dynamic modulus of elasticity $E_{\mathrm{d}}$ (GPa)\raisebox{-.43pc}{}} &
\parbox[t]{1.35cm}{\centering Average dynamic modulus of elasticity $E_{\mathrm{s}}$ (GPa)} &
\parbox[t]{1.35cm}{\centering Average static modulus of elasticity $E_{\mathrm{s}}$ (GPa)} &
\parbox[t]{1.35cm}{\centering Average static modulus of elasticity $E_{\mathrm{s}}$ (GPa)}\\
\endthead
U-SP-14-D-01 & 5.3139 & 3.2233 & 0.2089 &\morerows{8}{0.2076}& 
52.4669 &\morerows{8}{52.4081}& 0.5771359 &\morerows{8}{0.5706923}\\
U-SP-14-D-02 & 5.2984 & 3.2325 & 0.2036 && 52.5324 && 0.5778564 &\\
U-SP-14-D-03 & 5.3216 & 3.2141 & 0.2129 && 52.3388 && 0.5757268 &\\
U-SP-14-D-04 & 5.2907 & 3.2051 & 0.2101 && 51.9258 && 0.5711838 &\\
U-SP-14-D-05 & 5.3061 & 3.2325 & 0.2049 &  & 52.5922 &  & 0.5785142 & \\
U-SP-14-D-06 & 5.3092 & 3.2233 & 0.2081 && 52.4319 && 0.5767509 &\\
U-SP-14-D-07 & 5.3031 & 3.2417 & 0.2017 && 52.7511 && 0.5802621 &\\
U-SP-14-D-08 & 5.2984 & 3.2141 & 0.2089 && 52.1665 && 0.5216650 &\\
U-SP-14-D-09 & 5.3139 & 3.2233 & 0.2089 && 52.4669 && 0.5771359 &
\botline
\end{tabular}
{\vspace*{.5pc}}
\end{table}

\begin{table}%tabA3
\footnotesize\tabcolsep=4pt
\caption{\label{tabA3}Detailed results of ultrasonic velocities
${V}_{\mathrm{p}}$, ${V}_{\mathrm{s}}$, Poisson's ratio
${\vartheta}$, and the moduli ${E}_{\mathrm{d}}$ and
${E}_{\mathrm{s}}$ for specimens at 28~days}
\begin{tabular}{ccccccccc}
\thead
& \parbox[t]{1.4cm}{\centering Longitudinal wave $V_{\mathrm{p}}$\ $(\mathrm{m/s})\times 10^{3}$} &
\parbox[t]{1.25cm}{\centering Shear wave~$V_{\mathrm{s}}$ $(\mathrm{m/s}) \times 10^{3}$} &
\parbox[t]{.9cm}{\centering Poison's ratio $\vartheta$} &
\parbox[t]{.95cm}{\centering Average Poison's ratio $\vartheta$} &
\parbox[t]{1.2cm}{\centering Dynamic modulus of elasticity $E_{\mathrm{d}}$ (GPa)\raisebox{-.43pc}{}} &
\parbox[t]{1.35cm}{\centering Average dynamic modulus of elasticity $E_{\mathrm{s}}$ (GPa)} &
\parbox[t]{1.35cm}{\centering Average static modulus of elasticity $E_{\mathrm{s}}$ (GPa)} &
\parbox[t]{1.35cm}{\centering Average static modulus of elasticity $E_{\mathrm{s}}$ (GPa)}\\
\endthead
U-SP-28-D-01 & 6.2867 & 3.7965 & 0.2129 &\morerows{8}{0.2100}& 
73.0323 &\morerows{8}{73.0781}& 0.8033553 &\morerows{8}{0.8038076}\\
U-SP-28-D-02 & 6.2651 & 3.8093 & 0.2067 && 73.1464 && 0.8046104 &\\
U-SP-28-D-03 & 6.2543 & 3.7965 & 0.2083 && 72.7479 && 0.8002269 &\\
U-SP-28-D-04 & 6.2976 & 3.8093 & 0.2115 && 73.4352 && 0.8077872 &\\
U-SP-28-D-05 & 6.2436 & 3.7965 & 0.2067 &  & 72.6521 &  & 0.7991731 & \\
U-SP-28-D-06 & 6.2759 & 3.8093 & 0.2083 && 73.2432 && 0.8056752 &\\
U-SP-28-D-07 & 6.2651 & 3.7965 & 0.2098 && 72.8433 && 0.8012763 &\\
U-SP-28-D-08 & 6.2543 & 3.7965 & 0.2083 && 72.7479 && 0.8002269 &\\
U-SP-28-D-09 & 6.3415 & 3.8093 & 0.2177 && 73.8125 && 0.8119375 &
\botline
\end{tabular}
\end{table}

\printbibliography
\refinput{crmeca20250348-reference.tex}

\end{document}
