%~Mouliné par MaN_auto v.0.36.0 (4a84613e) 2025-10-13 16:19:48
\documentclass[CRMECA,Unicode,biblatex,published]{cedram}

\addbibresource{CRMECA_Kamagate_20250374.bib}

\usepackage{subcaption}
\usepackage{siunitx}
\usepackage{mathtools}
\usepackage[ruled]{algorithm2e}

\newcommand{\dd}{\mathrm{d}}
\newcommand{\nd}{\mathrm{nd}}
\newcommand{\cal}{\mathrm{cal}}
\newcommand{\err}{\mathrm{err}}
\newcommand{\ext}{\mathrm{ext}}
\newcommand{\tol}{\mathrm{tol}}
\newcommand{\dam}{\mathrm{dam}}
\newcommand{\diss}{\mathrm{d}}
\newcommand{\intr}{\mathrm{int}}
\newcommand{\therm}{\mathrm{th}}
\newcommand{\internal}{\mathrm{int}}

\DeclareMathOperator{\Arg}{Arg}
\DeclareMathOperator{\Div}{div}
\DeclareMathOperator{\grad}{grad}
\DeclareMathOperator{\stat}{stat}

\DeclareMathOperator{\divop}{div}
\newcommand{\divergence}{\divop}

\newcommand{\tq}[1]{\mathbb{#1}}
\newcommand{\bm}{\boldsymbol}

\newcommand{\bmnabla}{\bm{\nabla}}

\newcommand{\Twx}{T_{\omega}(\underline{x})}

\DeclareMathSymbol{\R}{\mathbin}{AMSb}{"52}

\makeatletter
\renewcommand\p@subfigure{}
\makeatother

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

%%% BAR & OVERLINE

%%% Better widebar macro
\makeatletter
\let\save@mathaccent\mathaccent
\newcommand*\if@single[3]{%
  \setbox0\hbox{${\mathaccent"0362{#1}}^H$}%
  \setbox2\hbox{${\mathaccent"0362{\kern0pt#1}}^H$}%
  \ifdim\ht0=\ht2 #3\else #2\fi
 }
%The bar will be moved to the right by a half of \macc@kerna, which is computed by amsmath:
\newcommand*\rel@kern[1]{\kern#1\dimexpr\macc@kerna}
%If there's a superscript following the bar, then no negative kern may follow the bar;
%an additional {} makes sure that the superscript is high enough in this case:
\newcommand*\widebar[1]{\@ifnextchar^{{\wide@bar{#1}{0}}}{\wide@bar{#1}{1}}}
%Use a separate algorithm for single symbols:
\newcommand*\wide@bar[2]{\if@single{#1}{\wide@bar@{#1}{#2}{1}}{\wide@bar@{#1}{#2}{2}}}
\newcommand*\wide@bar@[3]{%
  \begingroup
  \def\mathaccent##1##2{%
%Enable nesting of accents:
    \let\mathaccent\save@mathaccent
%If there's more than a single symbol, use the first character instead (see below):
    \if#32 \let\macc@nucleus\first@char \fi
%Determine the italic correction:
    \setbox\z@\hbox{$\macc@style{\macc@nucleus}_{}$}%
    \setbox\tw@\hbox{$\macc@style{\macc@nucleus}{}_{}$}%
    \dimen@\wd\tw@
    \advance\dimen@-\wd\z@
%Now \dimen@ is the italic correction of the symbol.
    \divide\dimen@ 3
    \@tempdima\wd\tw@
    \advance\@tempdima-\scriptspace
%Now \@tempdima is the width of the symbol.
    \divide\@tempdima 10
    \advance\dimen@-\@tempdima
%Now \dimen@ = (italic correction / 3) - (Breite / 10)
    \ifdim\dimen@>\z@ \dimen@0pt\fi
%The bar will be shortened in the case \dimen@<0 !
    \rel@kern{0.6}\kern-\dimen@
    \if#31
      \overline{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@}%
      \advance\dimen@0.4\dimexpr\macc@kerna
%Place the combined final kern (-\dimen@) if it is >0 or if a superscript follows:
      \let\final@kern#2%
      \ifdim\dimen@<\z@ \let\final@kern1\fi
      \if\final@kern1 \kern-\dimen@\fi
    \else
      \overline{\rel@kern{-0.6}\kern\dimen@#1}%
    \fi
 }%
  \macc@depth\@ne
  \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar
  \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax}%
  \macc@set@skewchar\relax
  \let\mathaccentV\macc@nested@a
%The following initialises \macc@kerna and calls \mathaccent:
  \if#31
    \macc@nested@a\relax111{#1}%
  \else
%If the argument consists of more than one symbol, and if the first token is
%a letter, use that letter for the computations:
    \def\gobble@till@marker##1\endmarker{}%
    \futurelet\first@char\gobble@till@marker#1\endmarker
    \ifcat\noexpand\first@char A\else
      \def\first@char{}%
    \fi
    \macc@nested@a\relax111{\first@char}%
  \fi
  \endgroup
}
\makeatother

\let\oldbar\bar
\renewcommand*{\bar}[1]{{\mathchoice{\widebar{#1}}{\widebar{#1}}{\widebar{#1}}{\oldbar{#1}}}}

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

%%% DELIMITERS

%\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\braces}{\{}{\}}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert}
\DeclarePairedDelimiter{\bracks}{[}{]}

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


\graphicspath{{./figures/}}

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

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

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}
\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}


\title{An incremental variational method to the coupling between gradient damage, thermoelasticity and heat conduction}
\alttitle{Une méthode variationnelle incrémentale pour le couplage entre gradient d’endommagement, thermoélasticité et conduction thermique}

\author{\firstname{Banouho} \lastname{Kamagaté}}
\address{Institut Jean Le Rond d'Alembert, CNRS UMR7190, Sorbonne Université, 4 place Jussieu, 75005 Paris, France}
\address{Laboratoire de Mécanique et Informatique, UFR Mathématique et Informatique, Université Félix Houphouët-Boigny, rue des écoles, Abidjan, Côte d'Ivoire}

\author{\firstname{Long} \lastname{Cheng}}
\address{GeoRessources Laboratory, Universit{\'e} de Lorraine (ENSG), CNRS UMR7359, 54518 Vandoeuvre-{l\`e}s-Nancy, France}

\author{\firstname{Radhi} \lastname{Abdelmoula}}
\address{Laboratoire des Sciences des Procédés et des Matériaux, Institut Galilée, Université Sorbonne Paris Nord, CNRS UPR3407, 93430 Villetaneuse, France}

\author{\firstname{Emile} \lastname{Danho}}
\address[2]{Laboratoire de Mécanique et Informatique, UFR Mathématique et Informatique, Université Félix Houphouët-Boigny, rue des écoles, Abidjan, Côte d'Ivoire}

\author{\firstname{Djimédo} \lastname{Kondo}\IsCorresp}
\address[1]{Institut Jean Le Rond d'Alembert, CNRS UMR7190, Sorbonne Université, 4 place Jussieu, 75005 Paris, France}
\email{djimedo.kondo@sorbonne-universite.fr}

\begin{abstract}
In this work, we propose an incremental variational approach to study the coupling between gradient damage, thermoelasticity and heat conduction phenomena. To this end, we first extend the thermodynamics of linear thermoelasticity to incorporate gradient damage phenomena. {After carefully} introducing the concept of \emph{kinetic entropy} to describe the interaction between thermoelasticity and heat conduction, this extension is implemented to establish a four-field {incremental} energy minimization procedure. By considering a suitable \emph{kinetic entropy} approximation, the latter is then consistently reduced to a three-field (displacement, damage, and absolute temperature) dependency, numerically implemented by means of a staggered optimization algorithm. Applications consisting in a study of the {cracking} of a plate under thermal shocks are considered. The approach is shown to deliver reliable predictions, based on comparison to available experimental observations which is also provided.
\end{abstract}

\begin{altabstract}
Dans ce travail, nous proposons une approche variationnelle incrémentale pour étudier le couplage entre les phénomènes d’endommagement non local, thermoélasticité et conduction thermique. À cette fin, nous étendons d’abord le cadre thermodynamique de la thermoélasticité linéaire afin d’y intégrer les phénomènes d’endommagement non local. Après avoir soigneusement introduit le concept d’\emph{entropie cinétique} pour décrire l’interaction entre thermoélasticité et conduction thermique, cette extension est mise en œuvre pour établir une procédure de minimisation de l’énergie incrémentale à quatre champs. En considérant une approximation adaptée de l’\emph{entropie cinétique}, cette procédure est ensuite réduite de manière cohérente à une dépendance en trois champs (déplacement, endommagement et température absolue), mise en œuvre numériquement au moyen d’un algorithme d’optimisation alternée. Des applications sont proposées, consistant en l’étude de la fissuration d’une plaque soumise à des chocs thermiques. L’approche montre sa capacité à fournir des prédictions fiables, une comparaison avec les observations expérimentales disponibles étant également présentée.
\end{altabstract}

\keywords{\kwd{Gradient damage model} \kwd{thermoelasticity} \kwd{heat conduction} \kwd{thermodynamics-based formulation} \kwd{incremental variational approach}}
\altkeywords{\kwd{Modèle à gradient d’endommagement} \kwd{thermoélasticité} \kwd{conduction thermique} \kwd{formulation fondée sur la thermodynamique} \kwd{approche variationnelle incrémentale}}

\editornote{Article submitted by invitation}
\alteditornote{Article soumis sur invitation}

\dateposted{2025-10-20}
\begin{document}

\maketitle

\vspace{-.5\baselineskip}

\section{Introduction}\label{sec:introduction}

Thermoelasticity theory describes the interaction between deformations and temperature variations in materials and constitutes a fundamental framework in the physics of continuous media. The first variational formulation of coupled linearized thermoelasticity can be traced back to the pioneering work of~\cite{Biot1956} (see also~\cite{Biot1957}), which was based on the thermodynamics of irreversible processes and employed Fourier's law for heat conduction. This approach has been later followed and extended by several works, such as~\cite{boley1962, BEN-AMOZ1965, Herrmann1963, NOWACKI1977, Nowacki1986, Batra1989}. These studies contributed to a systematic treatment of the coupled equations governing mechanics and heat conduction. However, most of these approaches do not provide a variational framework; that is, they do not allow for the construction of an energy functional whose minimization allows to rigorously yield the field equations of linear thermoelasticity. As a result, and strictly speaking, models proposed in this framework do not possess a fully variational structure; rather, they fall within the scope of so-called quasi-variational formulations (see for instance~\cite{finlayson1967,lebon1980}). An alternative formulation of thermoelasticity has been proposed by~\cite{green1993} who distinguished three types of coupled thermoelasticity models among which Type~II which is of particular interest for variational formulation due to the fact that its heat conduction law is dissipationless (at the difference of Fourier law). A systematic analysis of these three Green--Naghdi types of models has been done by~\cite{bargmann2006}. These studies highlight the lack of a truly well-established variational formulation for coupled linear thermoelasticity with Fourier type heat conduction, and they also lack an explicit recognition of the need for thermodynamics-based variational formulations capable of accounting for the coupling between thermoelasticity, heat conduction, and dissipative phenomena.

To address the issue of establishing variational formulations in presence of thermomechanical couplings, \cite{Stainier2006} proposed an incremental variational approach based on thermodynamics of irreversible processes. Their approach consists in introducing an integration factor that is based on a time rescaling of the dissipation potential (see also~\cite{Stainier2013}). This factor, expressed as the ratio of the external temperature to the equilibrium temperature, must remain equal to one throughout the system at equilibrium. Essentially, the equality of these two temperatures acts as an internal constraint in the variational analysis, which relaxes at equilibrium and enables the incorporation of strong coupling between damage and thermal mechanisms during entropy production. Independently, in the same year, \cite{peigney2006} proposed a time-discretization scheme, in a study on shape memory alloys, to obtain solutions to the incremental problem by optimizing the corresponding functional $\mathcal{F}$. Following these works, a series of studies have been proposed in the incremental variational setting. Among them, particular attention can be paid to~\cite{peigney2013,stainier2010, Stainier2013} and more recently~\cite{Mosler2015,romero2021,teichtmeister2022}. Note also the study by~\cite{bargmann2006} who also attempted to provide an incremental variational formulation for the three types of Green and Naghdi's thermoelasticity models.

Despite their interest, these approaches still need to be extended in order to account for non-local dissipative mechanisms such as in plasticity or damage processes. The present study aims to extend the incremental variational formulation to non local damage models for which there is no equivalent publication, except for~\cite{teichtmeister2022}, in which a four-field mixed variational approach was recently proposed without numerical implementation nor simulation.

Coming to the coupling between gradient damage with thermoelasticity, it has attracted increasing attention, particularly for simulating thermomechanical degradation phenomena such as thermal shock-induced cracking. This has led to a variety of modeling strategies, mainly based on phase-field methods, where the temperature field is often assumed to be known, typically obtained from an independent solution to the heat equation. This simplification allows for a variational formulation restricted to the displacement and damage fields. In this context, one may mention the works of~\cite{bourdin2014, sicsic2014, Mandal2021, dhahri2023}, in which the evolution of the temperature field is solved independently of the mechanical problem. The resulting thermal field is then used as a prescribed input in the phase-field model. Another class of models, such as those proposed in~\cite{chu2017, zheng2020, Zhang2023, li2023}, lacks a variational structure and does not rigorously integrate damage within a thermoelastic coupling framework.

In the present study, by taking advantage of thermodynamics of irreversible processes, we aim at developing variational gradient damage models with coupled thermoelasticity. Moreover, based on some experimental data (see for instance~\cite{Chrysochoos2003}) and following~\cite{peigney2006}, the intrinsic contribution of the damage in the heat production could be neglected. In this perspective, the above mentioned integration factor will not be considered in the corresponding variational analysis. Consequently, the variational structure which is preserved and addressed herein corresponds to a so-called moderate coupling between the damage and thermoelasticity. This is of course an approximation but will lead to easier calculations which are expected to yield sufficiently accurate predictions. {It is important to emphasize that this moderate coupling framework still accounts for the influence of thermoelasticity on entropy evolution, which remains affected by the presence of damage.}

The paper is organized as follows. Based on the framework of thermodynamics of irreversible processes, we develop in Section~\ref{sec: general variational principle} an incremental variational principle for the coupled problem, which involves four fields: displacement, damage, temperature, and the so-called kinetic entropy.

The latter is then carefully approximated in each increment that allows to reduce and express the incremental variational principle as based on functional of the three former fields. This formulation is further applied in Section~\ref{sec:variational model} to obtain a variational model for the concerned coupled system.

Finally, by making use of a staggered optimization algorithm, we apply the proposed model to simulate cracking of a plate under thermal shocks loadings. A comparison to available experimental observations of~\cite{shao2011} is also provided.

\section{General formulation of incremental variational principle for the coupled problem}\label{sec: general variational principle}

This section deals with a thermodynamics based formulation of an incremental variational framework for gradient damage process coupled with thermoelasticity as well as heat conduction. We assume that the material is subject to infinitesimal strains and small temperature variations.

\subsection{Thermodynamics of the coupling between linearized thermoelasticity and gradient damage}\label{sec:thermodynamics framework}

Let us consider a continuum medium, denoted by $\Omega$, for which the global form of the first and second principles of thermodynamics are classically expressed as:
\begin{equation}\label{eq:global form thermo}
\dot{\mathcal{E}} = - \mathcal{P}_{\internal} + \mathcal{P}_{\cal} ;
\qquad
\dot{{\mathcal{S}}} \geq \int_{\Omega} \frac{{r}}{T} \, \dd \Omega - \int_{\partial \Omega} \frac{\bm{{q}} \cdot {\bm{n}}}{T} \, \dd S
\end{equation}
where $\mathcal{E}$ is the global internal energy; $\mathcal{P}_{\internal}$ accounts for the power of internal work including that of possible microforces; $\mathcal{P}_{\cal}$ symbolizes the calorific power received by the system, which is furnished by the external heat supplies per unit volume, denoted by $r$, and per unit surface, denoted by $\bm{q}$; $\mathcal{S}$ gives the global entropy; $\bm{n}$ is the unit outward normal on the boundary of $\partial \Omega$ and $T$ the local absolute temperature.

In order to incorporate the gradient damage mechanism, we mainly rely on the work of~\cite{fremond1996damage} by introducing the internal variables such as the damage, denoted by $\alpha$, and its gradient $\bm{\nabla} \alpha$, with the latter characterizing the non-local effect of damage. In this context, the power of interior forces not only requires the classical term ${\bm{\sigma}} : \dot{{\bm{\varepsilon}}}$, where $\bm{\sigma}$ denotes the Cauchy stress tensor; $\dot{{\bm{\varepsilon}}}$ is the rate of the linearized strain tensor $\bm{\varepsilon}= \bm{\nabla}_s \bm{u}$ (with $\bm{u}$ the displacement vector and $\bm{\nabla}_s$ the symmetric gradient operator), but also incorporates the microscopic forces associated to damage $\alpha$ and its gradient $\bmnabla \alpha$, which are respectively denoted by ${Y}_{\alpha}$ and $ {\bm{Y}}_{\bmnabla {\alpha}} $. One has then:
\begin{equation}\label{eq:P_int}
\mathcal{P}_{\internal} = -
\int_{\Omega} {\bm{\sigma}} : \dot{{\bm{\varepsilon}}} \, \dd \Omega -
\int_{\Omega} \parens[\big]{{Y}_{\alpha} \dot{\alpha} + {\bm{Y}}_{\bmnabla {\alpha}} \cdot \bmnabla \dot{{\alpha}}
} \, \dd \Omega.
\end{equation}

As classically, let us introduce the Helmholtz free energy per unit volume, denoted by $w$, a function of the state variables, such that:
\begin{equation}\label{eq:Helmholz bis}
w \coloneqq w (\bm{\varepsilon}, T, \alpha, \bm{\nabla} \alpha).
\end{equation}
Note that $w$ is separately convex with respect to $\bm{\varepsilon}$, $\alpha$ and $\bmnabla \alpha$, but concave with respect to $T$.

The state laws which furnish the reversible forces are then derived as (see Appendix~\ref{sec:thermodynamics appendix}):
\begin{equation}\label{eq:state eqs with w bis}
\bm{\sigma}^{\nd} = \frac{\partial w}{\partial \bm{\varepsilon}};
\qquad s = - \frac{\partial w}{\partial {T}};
\qquad Y_{\alpha}^{\nd} = \frac{\partial w}{\partial {\alpha}};
\qquad \bm{Y}_{\bm{\nabla}{\alpha}}^{\nd} = \frac{\partial w}{\partial {\bmnabla \alpha}}
\end{equation}
with $s$ the entropy per unit volume of the solid. The superscripts ``nd'' denotes the non-dissipative part of the related thermodynamic forces.

Also, the Clausius--Duhem inequality (positivity of the dissipation) can be derived as:
\begin{equation}\label{eq:C-D inequality in w}
{\mathcal{D}} =
\underbrace{\vphantom{\frac{T}{T}}\bm{\sigma} : \dot{\bm{\varepsilon}} - \dot{w} - s \dot{T} + {Y}_{\alpha} \dot{\alpha} + {\bm{Y}}_{\bmnabla {\alpha}} \cdot \bmnabla \dot{{\alpha}}}_{{\mathcal{D}}_{\intr}}
\underbrace{- \frac{\bm{\nabla} T}{T} \cdot \bm{q}}_{{\mathcal{D}}_{\therm}}
\geq 0.
\end{equation}
This local dissipation $\mathcal{D}$ can be classically decomposed into two parts: the intrinsic dissipation ${\mathcal{D}}_{\intr}$ and the thermal dissipation due to heat conduction ${\mathcal{D}}_{\therm}$.

Since the strain field $\bm{\varepsilon}$ is non-dissipative in thermoelasticity, it follows that the dissipative component of the associated thermodynamic force vanish, i.e.\ $\bm{\sigma}^{\diss} = \bm{\sigma} - \bm{\sigma}^{\nd}= 0$. Consequently, eq.~\eqref{eq:C-D inequality in w} reduces then to:
\begin{equation}\label{eq:C-D inequality in w simplified}
{\mathcal{D}} =
\underbrace{\vphantom{\frac{T}{T}}\bm{Y}_{{\alpha}}^{\diss} \dot{\alpha} + \bm{Y}_{\bm{\nabla}{\alpha}}^{\diss} \cdot \bmnabla \dot{\alpha}}_{{\mathcal{D}}_{\intr}}
\underbrace{- \frac{\bm{\nabla} T}{T} \cdot \bm{q}}_{{\mathcal{D}}_{\therm}}
\geq 0
\end{equation}
where $\bm{Y}_{{\alpha}}^{\diss}$ and $\bm{Y}_{\bm{\nabla}{\alpha}}^{\diss}$ are the dissipative thermodynamic forces of $\alpha$ and $\bmnabla \alpha$, satisfying $ Y_{\alpha}^{\diss} = Y_{\alpha} - Y_{\alpha}^{\nd}$ and $
\bm{Y}_{\bm{\nabla}{\alpha}}^{\diss} = \bm{Y}_{\bm{\nabla}{\alpha}} - \bm{Y}_{\bm{\nabla}{\alpha}}^{\nd}$, respectively.

To characterize the evolution of dissipative systems in the context of Generalized Standard Materials (GSM) framework~\cite{Nguyen1975}, one can resort to the dissipation potential, denoted by $\varphi$. Following~\cite{Stainier2006}, and in the case of a coupling between thermoelasticity and gradient damage, a joint form can be adopted:
\begin{equation}\label{eq:dissipation potential 1}
\varphi\parens[\Bigg]{\dot{\alpha}, \bm{\nabla} \dot{\alpha}, - \frac{\bmnabla T}{T} ; \alpha, \bm{\nabla} \alpha, T} = \varphi^{M} (\dot{\alpha}, \bm{\nabla} \dot{\alpha}; \alpha, \bm{\nabla} \alpha) - \varphi^{T}\parens[\Bigg]{- \frac{\bmnabla T}{T}; T}
\end{equation}
where $\varphi^{M}$ and $\varphi^{T} $ respectively correspond to the damage dissipation and that of heat conduction. The latter owns the positiveness of the dissipation associated with heat conduction. Note that both $\varphi^{M}$ and $\varphi^{T}$ are convex with respect to their respective arguments. It follows from eq.~\eqref{eq:dissipation potential 1} that $\varphi$ is convex with respect to $\dot{\alpha}$ and $\bmnabla \dot{\alpha}$, but concave with respect to $- {\bmnabla T} / {T}$. One has:\footnote{In the case where the potential $\varphi$ is not differentiable, one should refer to the notion of subdifferential $\partial \varphi$, e.g.\ the evolution law for $\alpha$ then taking the form: $\mathcal{A}_{\alpha}\,\in\,\partial_{\dot{\alpha}} \varphi$. The irreversible force $\mathcal{A}_{\alpha}$ is said to belong to the sub-gradient of $\varphi$ at the considered point.}
\begin{equation}\label{eq:thermo forces d}
{Y}^{\diss}_{{\alpha}} = \frac{\partial \varphi}{\partial \dot{\alpha}}; \qquad
\bm{Y}^{\diss}_{\bm{\nabla} {\alpha}}= \frac{\partial \varphi}{\partial \bmnabla \dot{\alpha}};
\qquad \bm{q}= \frac{\partial \varphi^{T}}{\partial \parens[\big]{- \frac{\bmnabla T}{T}}}
\end{equation}
Concerning the dissipative processes, we recall that the corresponding dissipation potential is assumed to take a joint form (see eq.~\eqref{eq:dissipation potential 1}), wherein the damage dissipation potential $\varphi^M$ and that of the thermal conduction $\varphi^T$ are combined but remain decoupled. Since this study aims to provide a first attempt at formulating an incremental variational framework for modeling gradient damage mechanisms coupled with thermoelasticity and heat conduction, a \emph{moderate} coupling framework, as mentioned in Section~\ref{sec:introduction}, is adopted, in which the damage process affects the entropy evolution only through the thermoelastic potential.

For more details on the thermodynamic formulations presented above, readers are referred to Appendix~\ref{sec:thermodynamics appendix}. In the following part of this section, inspired by~\cite{Stainier2006} (see also~\cite{Stainier2013}), a rate form of the total energy functional which involves the two above potentials will be proposed. Its integral over the time increment will allow to construct an incremental energy functional. Stationarity of the latter will allow to establish the solution of the concerned coupled problem.

\subsection{Incremental variational principle}\label{sec:variational principle}

Let us consider a system subjected to some mechanical and/or {external thermal actions}. Specifically, the mechanical loads could be composed of surface forces $\bm{T}$ applied on the Neumann boundary $\partial_{\bm{T}} \Omega$ and an external displacement $\bm{\bar{u}}$ on the Dirichlet part $\partial_{\bm{u}} \Omega$, where $\partial_{\bm{T}} \Omega \, \cup \, \partial_{\bm{u}} \Omega = \partial \Omega$ and $\partial_{\bm{T}} \Omega \, \cap \, \partial_{\bm{u}} \Omega = \emptyset$. Concerning the external thermal charge, the boundary of the medium could be subjected to a heat supply ${Q}$ on its Neumann surface $\partial_{Q} \Omega$ and prescribed temperature $\bar{T}$ on the Dirichlet border $\partial_{T} \Omega$, for which we have $\partial_{Q} \Omega \, \cup \, \partial_{T} \Omega = \partial \Omega$ and $\partial_{Q} \Omega \, \cap \, \partial_{T} \Omega = \emptyset$. Furthermore, a body force per unit volume, denoted by $\bm{F}$, can be considered. It follows that the external power functional reads (see also~\cite{Stainier2006}):
\begin{equation}\label{eq:external work}
P_{\ext} \coloneqq 
\int_{\Omega} \bm{F} \cdot \dot{\boldsymbol{u}} \, \dd \Omega +
\int_{\partial_{\bm{T}} \Omega} \bm{T} \cdot \dot{\boldsymbol{u}} \, \dd S +\int_{\partial_{Q} \Omega} Q \, \log \frac{T}{T_0} \, \dd S.
\end{equation}

Since we consider the Helmholtz free energy, it may be more appropriate to consider the temperature field to be in equilibrium at any given time $t$. Yet, due to thermal diffusion, heat tends to enter the material point and alter its entropy, which corresponds to some kinetic behaviors. To describe this process, we assume that the entropy change occurs at a non-equilibrium state, denoted as time $t^{-}$. In this perspective and inspired by recent work by~\cite{Zhang2024JMPS} on the coupling between gradient damage, poroelasticity, and fluid flow in porous media, we distinguish between equilibrium entropy $s$ and a so-called \emph{kinetic entropy}, denoted by $\hat{s}$, which must be equivalent at equilibrium states at time $t$. Furthermore, the reversible process of the concerned open system is described by the power $\dot{w} + \hat{s} \dot{T}$. Following~\cite{Stainier2006, Stainier2013}, and inspired by~\cite{Ortiz1999} for elasto(visco)-plasticity, \cite{han1997,nguyen2016quasi} for gradient plasticity, or~\cite{petryk2003}, etc., let us introduce the rate form of total energy functional for damage coupled with thermoelastic behavior:
\begin{equation}\label{eq:total energy functional Pi 0}
\Pi(\dot{\bm{u}}, \dot{T}, {s}, \dot{\alpha}; T) = \int_{\Omega} \dot{w} (\dot{\boldsymbol{\varepsilon}}, \dot{T}, \dot{\alpha}, \bmnabla \dot{\alpha}) + \hat{s} \dot{T} + \varphi \parens[\Bigg]{\dot{\alpha}, \bm{\nabla} \dot{\alpha}, - \frac{\bmnabla {T}}{T} ; \alpha, \bm{\nabla} \alpha, T} \, \dd \Omega - {P}_{\ext}.
\end{equation}

\begin{rema}
Regarding the concept of \emph{kinetic entropy}, note that the overall system is thermodynamically open due to the presence of external heat sources and/or fluxes. At the local scale, each material point is assumed to be in a thermodynamically non-equilibrium state at time $t^{-}$ due to the thermal diffusion, whereby heat tends to enter the point and modify its entropy. The same material point is then considered to reach equilibrium at time $t$. The characterization of this kinetic entropy will be done in the following by considering eqs.~\eqref{eq:Euler--Lagrange t_n+1} and~\eqref{eq:Euler--Lagrange interval} in the variational point of view.
\end{rema}

We now proceed to establish an incremental variational problem in order to reduce the time-dependent nature such that the solution is expected to be given by the optimization of an incremental total energy functional. Specifically, the total time path $[t_0, t_{f}]$ is discretized into a sequence of time steps: $t_0 = 0, \dots, t_n, t_{n+1}, \dots, t_N = t_f$. The state at each step is determined by consistently optimizing within each time interval $[t_n, t_{n+1}]$. Given the known state $(\bm{\varepsilon}_n, T_n, \hat{s}_n, \alpha_n)$ at time $t_n$, we seek to find the corresponding unknown state $(\bm{\varepsilon}_{n+1}, T_{n+1}, \hat{s}_{n+1}, \alpha_{n+1})$ at $t_{n+1}$. In this context, we define the following functional to characterize the incremental evolution:
\begin{equation}\label{eq:total energy functional I_0}
I_n(\bm{\varepsilon}_{n+1}, T_{n+1}, \hat{s}_{n+1}, \alpha_{n+1}) = \int_{t_n}^{t_{n+1}} \Pi(\dot{\bm{u}}, \dot{T}, \hat{s}, \dot{\alpha}; T) \, \dd t.
\end{equation}

Taking variation of eq.~\eqref{eq:total energy functional I_0} respectively respect to $\bm{u}$, $T$ and $\alpha$, and enforcing the stationarities $t_{n+1}$ yield (see Appendix~\ref{sec:variational calculation} for details):
\begin{equation}\label{eq:Euler--Lagrange t_n+1}
\begin{cases}
	\divergence{\bm{\sigma}_{n+1}} +\bm{F}_{n+1}= \bm{0}; & \forall \bm{x} \in \Omega 
\\	\bm{\sigma}_{n+1} \cdot \boldsymbol{n}=\bm{T}_{n+1}; &\forall \bm{x} \in \partial_{\bm{T}} \Omega 
\\	\hat{s}_{n+1} = {s}_{n+1} ; & \forall \bm{x} \in \Omega 
\\	Y_{\alpha} \rvert_{n+1} - \divergence{\bm{Y}_{\bmnabla \alpha} \rvert_{n+1}} =0; &\forall \bm{x} \in {\Omega}
\\	\bm{Y}_{\bmnabla \alpha} \rvert_{n+1} \cdot \bm{n}=0; & \forall \bm{x} \in {\partial \Omega}
\end{cases}
\end{equation}
and in the interval $[t_n, t_{n+1}]$:
\begin{equation}\label{eq:Euler--Lagrange interval}
\begin{dcases}
	\dot{\hat{s}} + \divergence{\frac{\bm{q}}{T}} = {0}; & \forall \bm{x} \in \Omega 
\\	\bm{q} \cdot \bm{n} = Q; & \forall \bm{x} \in \partial_{Q} \Omega 
\\	\bracks[\Bigg]{\frac{\partial \varphi}{\partial \alpha} - \frac{\dd}{\dd t} \parens[\Bigg]{\frac{\partial \varphi}{\partial \dot{\alpha}}}} -\divop\bracks[\Bigg]{\frac{\partial \varphi}{\partial \bmnabla {\alpha}} - \frac{\dd}{\dd t} \parens[\Bigg]{\frac{\partial \varphi}{\partial \bmnabla \dot{\alpha}}}
} =0 ; & \forall \bm{x} \in \Omega
\\	\bracks[\Bigg]{\frac{\partial \varphi}{\partial \bmnabla {\alpha}} - \frac{\dd}{\dd t} \parens[\Bigg]{\frac{\partial \varphi}{\partial \bmnabla \dot{\alpha}}}} \cdot \bm{n}=0 ; & \forall \bm{x} \in {\partial \Omega}
\\	\dot{T} = 0 ; & \forall \bm{x} \in {\Omega}
\end{dcases}
\end{equation}
eqs.~\eqref{eq:Euler--Lagrange t_n+1} furnish the Euler--Lagrange equations of $I_n$ and the Neumann boundary conditions at $t_{n+1}$, while eqs.~\eqref{eq:Euler--Lagrange interval} allow to determine the optimization path. Their derivation is presented in detail in Appendix~\ref{sec:variational calculation}. All of them are expressed thanks to the thermodynamics potentials allowing to take into account the non-local damage effects. Essentially, these equations make $I_n$ to be stationary in $[t_n, t_{n+1}]$. Let us suppose the thermodynamics potentials as classically being quadratic, and recall their convexity as discussed in Section~\ref{sec:thermodynamics framework}. It follows that $I_n$ must be convex with respect to $u_{n+1}$ and $\alpha_{n+1}$ but concave with respect to $T_{n+1}$. We can hence establish the following incremental variational principle:
\begin{equation}\label{eq:incremental variational principle}
\parens{\boldsymbol{u}_{n+1}, T_{n+1}, \hat{s}_{n+1}, \alpha_{n+1}} = \Arg \inf_{\boldsymbol{u}'_{n+1}, \alpha'_{n+1}}
\stat{{T}{'}_{n+1}, \hat{s}'_{n+1}}
\braces[\big]{I_n\parens{\bm{u}{'}_{n+1}, T{'}_{n+1}, \hat{s}'_{n+1}, \alpha{'}_{n+1}}
}
\end{equation}
with $\parens{\boldsymbol{u}_{n+1}, T_{n+1}, \hat{s}_{n+1}, \alpha_{n+1}}$ henceforth the estimate of the expected solution and $(\boldsymbol{u}{'}_{n+1}$, $T{'}_{n+1}, \hat{s}'_{n+1} $, $\alpha{'}_{n+1})$ denoting the admissible fields of the displacement, the absolute temperature, the \emph{kinetic entropy} as well as that of the damage, respectively.

\subsection{A reduced form of the variational scheme}\label{sec:different principle}

As demonstrated above, the application of the incremental variational principle~\eqref{eq:incremental variational principle} necessitates the calculation of the optimization path, whose numerical implementation might be performed via a fully implicit algorithm. However, as discussed in~\cite{Ortiz1989} (see also~\cite{Stainier2006, Zhang2024JMPS}), this path calculation could be treated as avoidable except for some particular case. Instead, a \emph{weak} quantification of $I_n$ that aligns with the field equations at $t_{n+1}$ could be sufficient. From this perspective, we employ backward Euler approximations for the following rate fields:
\begin{equation}\label{eq:numerical velocities - 0}
\dot{\bm{u}} \simeq \frac{\bm{u}_{n+1} - \bm{u}_{n}}{\Delta t};
\qquad
\dot{\alpha} \simeq \frac{\alpha_{n+1} - \alpha_{n}}{\Delta t};
\qquad
\dot{T} \simeq \frac{T_{n+1} - T_{n}}{\Delta t};
\end{equation}
where $\Delta t = t_{n+1} - t_n$.

Recall that, as aforementioned, the entropy simultaneously presents some \emph{equilibrium} and \emph{kinetic} behaviors at any given time.\footnote{Specifically, the entropy is kinetic at the non-equilibrium state $t^{-}$ and equilibrium at $t$.} In this perspective and in the context of incremental variational framework, we assume that the entropy is considered to be in a non-equilibrium (i.e., kinetic) state at time $t_n$, and it relaxes to an equilibrium state at $t_{n+1}$. The continuity of $\hat{s}$ in time dimension makes:
\begin{equation}\label{eq:s_n}
\hat{s} \simeq s_n = - \frac{\partial w}{\partial T} (\bm{u}_n, T_n, \alpha_n, \bmnabla \alpha_n).
\end{equation}
This expression actually corresponds to an optimization of the kinetic entropy over the previous increment. Indeed, as shown for instance in~\cite{Heuze2022}, the kinetic entropy $\hat{s}$ can be approximated at any time within the interval $[t_n, t_{n+1})$, except at the final time $t_{n+1}$, which is considered here as the stationary state to be determined. In summary, this type of approximation, as will be demonstrated in the following, ensures a consistent representation of the coupled constitutive relations and guarantees satisfactory convergence as $\Delta t \to 0$. By considering eqs.~\eqref{eq:Helmholz bis}, \eqref{eq:dissipation potential 1}, \eqref{eq:numerical velocities - 0} and~\eqref{eq:s_n} and making use of the rectangular rule for the time integration of $\varphi$ as well as the term of thermal boundary condition, the incremental energy $I_n$ can be calculated as:
\begin{multline}\label{eq:I^* delta t}
I_n(\bm{u}_{n+1}, T_{n+1}, \alpha_{n+1}) \simeq \int_{\Omega} \parens{w_{n+1} - w_{n}} + {s}_n \parens{T_{n+1} - T_{n}} \, \dd \Omega + \Delta t \int_{\Omega} \parens{\varphi^{M}_{n+1} - \varphi^{T}_{n+1}} \, \dd \Omega
\\	- \int_{\Omega} \bm{F}_{n+1} \cdot \parens{\bm{u}_{n+1} - \bm{u}_{n}} \, \dd \Omega - \int_{\partial_{\bm{T}} \Omega} \bm{T}_{n+1} \cdot \parens{\bm{u}_{n+1} - \bm{u}_{n}} \, \dd S - \Delta t \int_{\partial_{Q} \Omega} Q_{n+1} \, \log \frac{T_{n+1}}{T_n} \, \dd S.
\end{multline}
Let us compute the variation of eq.~\eqref{eq:I^* delta t} with respect to $\bm{u}_{n+1}, T_{n+1}, \alpha_{n+1} $; the stationarity gives:
\begin{equation}\label{eq:Euler--Lagrange delta t}
\begin{dcases}
	\divergence{\frac{\partial w_{n+1}}{\partial \bm{\varepsilon}_{n+1}}} + \bm{F}_{n+1} = 0; & \forall \bm{x} \in \partial \Omega
\\	\frac{\partial w_{n+1}}{\partial \bm{\varepsilon}_{n+1}} \cdot \bm{n} - \bm{T}_{n+1} = 0 ; & \forall \bm{x} \in \partial_{\bm{T}} \Omega
\\	{T_{n+1}} \frac{- \frac{\partial w_{n+1}}{\partial T_{n+1}} - {s}_n}{\Delta t} + \divergence{\frac{\partial \varphi^{T}_{n+1}}{\partial \parens[\big]{- \frac{\bmnabla {T_{n+1}}}{T_{n+1}}}}} = 0; & \forall \bm{x} \in \Omega
\\	\frac{\partial \varphi^{T}_{n+1}}{\partial \parens[\big]{- \frac{\bmnabla {T_{n+1}}}{T_{n+1}}}} \cdot \bm{n} - Q_{n+1}= 0; & \forall \bm{x} \in \partial_{Q} \Omega
\\	\parens[\Bigg]{\frac{\partial w_{n+1}}{\partial \alpha_{n+1}} + \frac{\partial \varphi^{M}_{n+1}}{\partial \dot{\alpha}_{n+1}}} - \divergence{\frac{\partial w_{n+1}}{\partial \bmnabla \alpha_{n+1}} +
\frac{\partial \varphi^{M}_{n+1}}{\partial \bmnabla \dot{\alpha}_{n+1}}} = \mathcal{O} (\Delta t); & \forall \bm{x} \in {\Omega}
\\	\parens[\Bigg]{\frac{\partial w_{n+1}}{\partial \bmnabla \alpha_{n+1}} + \frac{\partial \varphi^{M}_{n+1}}{\partial \bmnabla \dot{\alpha}_{n+1}}} \cdot \bm{n} = \mathcal{O} (\Delta t); & \forall \bm{x} \in {\partial \Omega}.
\end{dcases}
\end{equation}
Furthermore, by considering eqs.~\eqref{eq:state eqs with w bis} and~\eqref{eq:thermo forces d}, \eqref{eq:Euler--Lagrange delta t} can be readily recast into:
\begin{equation}\label{eq:Euler--Lagrange delta t bis}
\begin{dcases}
	\divergence{\bm{\sigma}_{n+1}} + \bm{F}_{n+1} = 0; & \forall \bm{x} \in \Omega
\\	\bm{\sigma}_{n+1} \cdot \bm{n} - \bm{T}_{n+1} = 0 ; & \forall \bm{x} \in \partial_{\bm{T}} \Omega
\\	T_{n+1} \frac{s_{n+1} - {s}_n}{\Delta t} + \divergence{\bm{q}_{n+1}} = 0; &\forall \bm{x} \in \Omega
\\	\bm{q}_{n+1}\cdot \bm{n} - Q_{n+1} = 0; & \forall \bm{x} \in \partial_{Q} \Omega
\\	{Y}_{\alpha} \rvert_{n+1} - \divergence{\bm{Y}_{\bm{\nabla} \alpha} \rvert_{n+1}} = \mathcal{O} (\Delta t); & \forall \bm{x} \in {\Omega}
\\	\bm{Y}_{\bm{\nabla} \alpha} \rvert_{n+1} \cdot \bm{n} = \mathcal{O} (\Delta t); & \forall \bm{x} \in {\partial \Omega}
\end{dcases}
\end{equation}
which are the governing equations (equilibrium, heat equations, damage criterion) and Neumann boundary conditions associated to $\bm{u}_{n+1}, T_{n+1}, \alpha_{n+1} $, respectively.

Recalling the convexity of $I_n(\bm{u}_{n+1}, T_{n+1}, \alpha_{n+1})$ with respect to $\bm{u}_{n+1}$ and $\alpha_{n+1}$ and the concavity with respect to $T_{n+1}$, the incremental variational principle of eq.~\eqref{eq:incremental variational principle} could be reduced into the following form as three fields dependent, such that:
\begin{equation}\label{eq:incremental variational principle - *}
\parens{\boldsymbol{u}_{n+1}, T_{n+1}, \alpha_{n+1}} = \Arg \inf_{\boldsymbol{u}_{n+1}{'}, \alpha_{n+1}{'}} \sup_{\smash{T_{n+1}{'}}}
\braces[\big]{{I}_n\parens{\boldsymbol{u}{'}_{n+1}, T{'}_{n+1}, \alpha{'}_{n+1} ; \hat{s}}
}
\end{equation}
$\hat{s}$ being aforehand identified in the current increment.

\section{A variational model for coupled thermoelasticity with gradient damage}\label{sec:variational model}

This section aims to apply the established incremental variational principle to a gradient damage model coupled with thermoelasticity and heat conduction. As mentioned before, we assume the material undergoes infinitesimal deformations and small temperature variations. Accordingly, the formulation of the variational model will rely on the choice of a suitable Helmholtz free energy, a damage dissipation potential, and the Fourier potential for heat conduction.

\subsection{Helmholtz free energy and dissipation potentials}

We propose as:
\begin{description}
\item[Helmholtz free energy for linear thermoelasticity]
\begin{equation}\label{eq:G0}
{w} (\bm{\varepsilon}, T, \alpha) =
\frac{1}{2} \bm{\varepsilon} : \tq{C} (\alpha) : \bm{\varepsilon} - s_0 T - (T-T_0) \bm{\beta}(\alpha) : \bm{\varepsilon} - \frac{c(\alpha)}{2 T_0} (T-T_0)^2
\end{equation}
where $\tq{C} (\alpha) $, $\bm{\beta}(\alpha)$, $c(\alpha)$ are respectively the stiffness tensor, the second order tensor of linear thermal expansion and the heat capacity. A priori, all these material parameters are affected by the damage.

Next, the reversible thermodynamic forces obtained from the state laws read:
\begin{equation}\label{eq:nd forces model}
\begin{aligned}
	& \bm{\sigma} = \tq{C} (\alpha) : \bm{\varepsilon} - (T-T_0) \bm{\beta}(\alpha);
\\	& s - s_0 = \bm{\beta}(\alpha) : \bm{\varepsilon} + \frac{c(\alpha)}{T_0} (T-T_0);
\\	& Y_{\alpha}^{\nd} = \frac{\partial {w}}{\partial \alpha};
\\	& \bm{Y}_{\bmnabla \alpha}^{\nd} = \bm{0}.
\end{aligned}
\end{equation}
\item[Dissipation potential] Let us recall that the dissipation potential is considered as in a joint form of the damage dissipation potential $\varphi^M$ and that for thermal conduction $\varphi^T$ (see eq.~\eqref{eq:dissipation potential 1}). Following~\cite{KPO2021, Zhang2024JMPS}, the dissipation potential $\varphi^M$ is chosen:
\begin{equation}\label{eq:Ds0-general}
{\varphi}^{M}(\dot{\alpha}, \bmnabla \dot{\alpha}; \alpha, \bmnabla \alpha) = Y_c(\alpha) \dot{\alpha} +2 l_{0}^2 w_1 \bmnabla \alpha \cdot \bmnabla \dot{\alpha}
\end{equation}
where $Y_c(\alpha)$ denotes the critical damage energy at current state, $l_{0}$ is the material internal length and $w_1$ the first threshold in energy release rate that can be identified from different damage descriptions (e.g.~\cite{pham2011gradient} for AT1 model and~\cite{bourdin2000numerical} for AT2 model).

In addition to satisfying the classical requirements (being a positive scalar-valued function, convex with respect to its arguments, and minimal when the arguments vanish), $\varphi^{M}$ is positively homogeneous of degree~1 and have the following remarkable property. Indeed, it corresponds to the so-called
\emph{simple dissipative systems},\footnote{See~\cite{Ehrlacher1989} for an account of simple dissipative systems.} for which the total dissipated energy until the current time is function only of the current value of internal variables.

Owing to this property, it readily follows that:
\begin{equation}\label{eq:dissipation energy Phi 0}
\int_{t_n}^{t_{n+1}} \varphi^{M} (\dot{\alpha}, \bmnabla \dot{\alpha}; \alpha, \bmnabla \alpha) \, \dd t = w_1 l_{0} \bracks[\Bigg]{\frac{\omega (\alpha_{n+1}) - \omega (\alpha_{n})}{l_{0}} + l_{0}
\parens[\big]{\parallel \bmnabla \alpha_{n+1} \parallel^2 - \parallel \bmnabla \alpha_{n} \parallel^2
}}.
\end{equation}

It is worthy to note here that eq.~\eqref{eq:dissipation energy Phi 0} allows to obtain a closed-form expression of the related incremental dissipation due to the simple dissipative process without the approximation made in eq.~\eqref{eq:I^* delta t}. Nevertheless, the incremental variational method proposed in Section~\ref{sec: general variational principle} remains applicable to more general cases.

Coming now to the dissipation potential $\varphi^T$ corresponding to thermal dissipation, we assume the validity of the Fourier potential to describe heat conduction:
\begin{equation}\label{eq:explicit fourier}
\varphi^{T} \parens[\Bigg]{- \frac{\bmnabla T}{T}; T} = \frac{T}{2} \parens[\Bigg]{- \frac{\bmnabla T}{T}} \cdot \bm{k} \cdot \parens[\Bigg]{- \frac{\bmnabla T}{T}} = \frac{\bmnabla T \cdot \bm{k} \cdot \bmnabla T}{2 T}
\end{equation}
where $\bm{k}$ is the isotropic thermal conductivity tensor that is assumed here as constant.

Since small temperature variations are considered, the force term in the Fourier potential can be approximated as (see for instance~\cite{Stainier2013}):
\begin{equation}\label{eq:appro fourier force}
-\frac{\bmnabla T}{T} \simeq - \frac{\bmnabla T}{T_r}
\end{equation}
with $T_r$ a reference temperature that is supposed to be the initial temperature $T_0$. Additionally, we consider that $\varphi^T$ is parametrically affected by $T_0$. Then, eq.~\eqref{eq:explicit fourier} can be readily reduced as:
\begin{equation}\label{eq:explicit fourier reduced}
\varphi^{T} \parens[\Bigg]{- \frac{\bmnabla T}{T_0}; T_0} = \frac{T_0}{2} \parens[\Bigg]{- \frac{\bmnabla T}{T_0}} \cdot \bm{k} \cdot \parens[\Bigg]{- \frac{\bmnabla T}{T_0}} = \frac{\bmnabla T \cdot \bm{k} \cdot \bmnabla T}{2 T_0}.
\end{equation}
\end{description}

\subsection{Derivation of the incremental functional \texorpdfstring{$I_n$}{I\_n} for the proposed variational model}

As discussed in Section~\ref{sec: general variational principle}, the non-dissipative process of the concerned system is described by $\dot{w} + \hat{s} \dot{T}$, in which the kinetic entropy was approximated as $\hat{s} \simeq s_n$. From this assumption, one has:
\[
\hat{s} \simeq - \frac{\partial w}{\partial T} ({\bm{\varepsilon}_n}, {T}_n, {\alpha}_{n}) = s_0 + \bm{\beta} (\alpha_{n}) : \bm{\varepsilon}_n + \frac{c(\alpha_n)}{T_0} (T_n-T_0)
\]
which can be reported as $s_n$ in the incremental energy functional~\eqref{eq:I^* delta t}.

As mentioned before, this choice $\hat{s} \simeq s_n$ is one of the possible approximations for describing the kinetic entropy. In order to well incorporate the effect of damage on the thermoelastic degradation during the variational calculation, we assume henceforward that the dependence of $\beta({\alpha})$ and $c(\alpha)$ on the damage field $\alpha$ in the current increment is incorporated through a semi-implicit scheme, such that:
\begin{equation}\label{eq: model equi porosity 0}
\hat{s} \simeq s_0 + \bm{\beta} (\alpha_{n+1}) : \bm{\varepsilon}_n + \frac{c(\alpha_{n+1})}{T_0} (T_n-T_0)
\end{equation}
which also corresponds to
\begin{equation}\label{eq:w at n calibrated}
{w}_n \simeq
\frac{1}{2} \bm{\varepsilon}_{n} : \tq{C} (\alpha_{n}) : \bm{\varepsilon}_{n} - s_0 T_{n} - (T_{n} -T_0) \bm{\beta}(\alpha_{n+1}) : \bm{\varepsilon}_n - \frac{c(\alpha_{n+1})}{2 T_0} (T_n-T_0)^2.
\end{equation}
It is important to note here that this approximation is variationally consistent.

In the increment $[t_n, t_{n+1}]$, these two equations lead to:
\begin{equation}\label{eq: model equi porosity}
\hat{s} \dot{T} \simeq
\hat{s} \frac{T_{n+1} - T_{n}}{\Delta t} = \frac{T_{n+1} - T_{n}}{\Delta t}
\bracks[\Bigg]{s_0 + \bm{\beta} (\alpha_{n+1}) : \bm{\varepsilon}_n + \frac{c(\alpha_{n+1})}{T_0} (T_n-T_0)
}
\end{equation}
and
\begin{multline}\label{eq: model phi - phi_n}
{w}_{n+1} - {w}_{n} \simeq
\frac{1}{2} \bm{\varepsilon}_{n+1} : \tq{C} (\alpha_{n+1}) : \bm{\varepsilon}_{n+1} - \frac{1}{2} \bm{\varepsilon}_{n} : \tq{C} (\alpha_{n}) : \bm{\varepsilon}_{n} - \parens{T_{n+1} - T_{n}} s_0
\\	- \bm{\beta} (\alpha_{n+1}) :
\bracks[\big]{
\parens{T_{n+1} - T_{0}} \bm{\varepsilon}_{n+1} - \parens{T_{n} - T_{0}} \bm{\varepsilon}_{n}
} - \frac{c(\alpha_{n+1})}{2 T_0}
\parens{T_{n+1} - T_{n}} \parens{T_{n+1} + T_{n} - 2 T_0}
\end{multline}

Still in the context of small temperature variations, the external power functional expressed in eq.~\eqref{eq:external work} can be reduced into:
\begin{equation}\label{eq:external work linear}
P_{\ext} \coloneqq 
\int_{\Omega} \bm{F} \cdot \dot{\boldsymbol{u}} \, \dd \Omega +
\int_{\partial_{\bm{T}} \Omega} \bm{T} \cdot \dot{\boldsymbol{u}} \, \dd S +\int_{\partial_{Q} \Omega} Q \, {\frac{T - T_0}{T_0}} \, \dd S.
\end{equation}

Gathering eqs.~\eqref{eq:Ds0-general}, \eqref{eq:dissipation energy Phi 0}, \eqref{eq:explicit fourier reduced}, \eqref{eq: model equi porosity}, \eqref{eq: model phi - phi_n} and~\eqref{eq:external work linear}, the incremental energy functional given by~\eqref{eq:I^* delta t} can be recast in the form:
\begin{equation}\label{eq:I_n Model}
\begin{split}
I_n (\bm{\varepsilon}_{n+1}, \alpha_{n+1}, T_{n+1})
	& = \int_{\Omega} \bracks[\Bigg]{\frac{1}{2} \bm{\varepsilon}_{n+1} : \tq{C} (\alpha_{n+1}) : \bm{\varepsilon}_{n+1} - \parens{T_{n+1} - T_0} \bm{\beta} (\alpha_{n+1}): \parens{\bm{\varepsilon}_{n+1} - \bm{\varepsilon}_{n}}
\\	& \quad - \frac{c(\alpha_{n+1})}{2 T_0} (T_{n+1}-T_n)^2
} \, \dd \Omega
\\	& \quad + \int_{\Omega} w_1 l_{0} \bracks[\Bigg]{\frac{\omega (\alpha_{n+1}) - \omega (\alpha_{n})}{l_{0}} + l_{0}
\parens[\big]{\parallel \bmnabla \alpha_{n+1} \parallel^2 - \parallel \bmnabla \alpha_{n} \parallel^2
}} \, \dd \Omega
\\	& \quad - \Delta t \int_{\Omega} \frac{\bmnabla T_{n+1} \cdot \bm{k} \cdot\bmnabla T_{n+1}}{2T_0} \, \dd \Omega - \int_{\Omega} \bm{F}_{n+1} \cdot \parens{\bm{u}_{n+1} - \bm{u}_{n}} \, \dd \Omega
\\	& \quad - \int_{\partial_{\bm{T}} \Omega} \bm{T}_{n+1}\cdot \parens{\bm{u}_{n+1} - \bm{u}_{n}} \, \dd S - \Delta t \int_{\partial_{Q} \Omega} Q_{n+1} {\frac{T_{n+1} - T_0}{T_0}} \, \dd S.
\end{split}
\end{equation}
Note that in the context of variational optimization, the quantities that depend only on the state at $t_n$ can be treated as constants in the above equation. Moreover, due to the form of the state and dissipation potentials expressed by~\eqref{eq:G0}, \eqref{eq:Ds0-general} and~\eqref{eq:explicit fourier reduced}, it can be readily understood that this incremental energy functional is quadratic in each of its arguments and convex with respect to $\bm{u}_{n+1}$ and $\alpha_{n+1}$, but concave with respect to $T_{n+1}$. Hence, the resolution of the concerned problem can be achieved by applying the incremental variational principle~\eqref{eq:incremental variational principle - *}. We provide below the governing equations and corresponding boundary conditions, which are obtained based on the variational calculus of eq.~\eqref{eq:I_n Model} with respect to $\bm{u}_{n+1}$, $T_{n+1}$ and $\alpha_{n+1}$:
\begin{equation}\label{eq:model eqs in tn+1}
\begin{dcases}
	\divop \bracks[\big]{\tq{C}(\alpha_{n+1}) : \bm{\varepsilon_{n+1}} - (T_{n+1} - T_0) \bm{\beta} (\alpha_{n+1})} + \bm{F}_{n+1} = \bm{0}; & \forall \bm{x} \in \Omega
\\	\bracks[\big]{
\tq{C}(\alpha_{n+1}) : \bm{\varepsilon}_{n+1} - (T_{n+1} - {T_{0}})\bm{\beta} (\alpha_{n+1})
} \cdot \bm{n} - \bm{T}_{n+1} = \bm{0}; & \forall \bm{x} \in \partial_{\bm{T}} \Omega
\\	{c(\alpha_{n+1})}
\frac{T_{n+1} - T_n}{\Delta t} + {T_0} \bm{\beta} (\alpha_{n+1}):
\frac{\bm{\varepsilon}_{n+1} - \bm{\varepsilon}_{n}}{\Delta t} - \divergence{\bm{k} \cdot \bmnabla T_{n+1}} = 0; &\forall \bm{x} \in \Omega
\\	\parens{- {\bm{k}} \cdot \bmnabla T_{n+1}} \cdot \bm{n} + Q_{n+1} = 0; & \forall \bm{x} \in \partial_{Q} \Omega
\\	\frac{\partial {w}_{n+1}}{\partial \alpha_{n+1}} + w_1 \frac{\partial \omega (\alpha_{n+1})}{\partial \alpha_{n+1}} - 2 w_1 l_{0}^2 {{\bmnabla}^2} \alpha_{n+1} = 0; & \forall \bm{x} \in {\Omega}
\\	\parens{2 w_1 l_{0}^2 {{\bmnabla}} \alpha_{n+1}} \cdot\bm{n} = 0; & \forall \bm{x} \in {\partial \Omega}
\end{dcases}
\end{equation}
with {${\bmnabla}^2$} the Laplacian operator.

\section{Numerical implementation and application}\label{sec:Numerical}

In this section, we deal with a numerical implementation procedure for the proposed variational model coupling gradient damage, thermoelasticity and heat conduction. Since its total energy functional~\eqref{eq:I_n Model} includes a quadratic term of the product between $\bmnabla{u}$ and $\alpha$, $I_n$ is hence separately but not globally convex with respect to them. This makes the corresponding global minimization unachievable. However, this issue does not arise between the displacement $\bm{u}$ and the temperature field $T$. For this reason, the variational optimization is performed through a semi-staggered algorithm, which is implemented in this work with the open-source computing platform FEniCS~\cite{AlnaesEtal2015}. This approach is then applied to the study of plate cracking under thermal shocks in quasi-static conditions.

\subsection{Numerical algorithm}\label{sec:algorithm}

In order to make the above mentioned semi-staggered optimization in the current increment $[t_n, t_{n+1}]$, we first switch off the damage variation, such that $\alpha \coloneqq \alpha_n$, and solving the following stationary problem where the thermoelasticity and heat conduction process are coupled:
\begin{equation}\label{eq:staggered 1}
\parens{\bm{u}_{n+1}, T_{n+1}} =
\Arg \stat{\bm{u}{'}_{n+1}, T{'}_{n+1}}
\braces[\big]{I^{\therm}_n\parens{\bm{u}{'}_{n+1}, T{'}_{n+1} ; \alpha_{n}}
}.
\end{equation}
By considering eq.~\eqref{eq:I_n Model} and neglecting the constant terms, $I^{\therm}_n$ can be explicitly expressed as:
\begin{multline}\label{eq:I^th bis}
I_n^{\therm} (\bm{\varepsilon}_{n+1}, T_{n+1}; \alpha_{n}) =
\int_{\Omega} \bracks[\Bigg]{\frac{1}{2} \bm{\varepsilon}_{n+1} : \tq{C} (\alpha_{n}) : \bm{\varepsilon}_{n+1} - \parens{T_{n+1} - T_0} \bm{\beta} (\alpha_{n}): \parens{\bm{\varepsilon}_{n+1} - \bm{\varepsilon}_{n}} 
\\	-  \frac{c(\alpha_{n})}{2 T_0} (T_{n+1}-T_n)^2} \, \dd \Omega - \Delta t \int_{\Omega} \frac{\bmnabla T_{n+1} \cdot \bm{k} \cdot\bmnabla T_{n+1}}{2T_0} \, \dd \Omega- \int_{\Omega} \bm{F}_{n+1} \cdot \parens{\bm{u}_{n+1} - \bm{u}_{n}} \, \dd \Omega
\\	- \int_{\partial_{\bm{T}} \Omega} \bm{T}_{n+1}\cdot \parens{\bm{u}_{n+1} - \bm{u}_{n}} \, \dd S - \Delta t \int_{\partial_{Q} \Omega} Q_{n+1} {\frac{T_{n+1} - T_0}{T_0}} \, \dd S.
\end{multline}
Next, we estimate the damage solution $\alpha_{n+1}$ from the following minimizer:
\begin{equation}\label{eq:staggered 2}
\alpha_{n+1} =
\Arg \min{\alpha{'}_{n+1}}
\braces[\big]{I^{\dam}_n\parens{\alpha{'}_{n+1}}
}
\end{equation}
where
\begin{multline}\label{eq:I^dam}
I^{\dam}_n (\alpha_{n+1}) =
\int_{\Omega} \bracks[\Bigg]{\frac{1}{2} \bm{\bar{\varepsilon}}_{n+1} : \tq{C} (\alpha_{n+1}) : \bm{\bar{\varepsilon}}_{n+1} - \parens{\bar{T}_{n+1} - T_0} \bm{\beta} (\alpha_{n+1}): \parens{\bm{\varepsilon}_{n+1} - \bm{\varepsilon}_{n}} 
\\	 - \frac{c(\alpha_{n+1})}{2 T_0} (\bar{T}_{n+1}-T_n)^2 + w_1 l_{0} \parens[\Bigg]{\frac{\omega (\alpha_{n+1})}{l_{0}} + l_{0}
\parallel \bmnabla \alpha_{n+1} \parallel^2}} \, \dd \Omega
\end{multline}
where the bar superscript denotes the optimization of $\bm{u}_{n+1}$ and $T_{n+1}$ computed from eq.~\eqref{eq:staggered 1}.

Specifically, eq.~\eqref{eq:staggered 1} is solved by using the \emph{NonlinearVariationalSolver} module in the open-source DOLFIN~\cite{LoggEtal2012}, which allows an automatic iteration in each time increment with user-defined tolerances of absolute and relative errors as well as the maximal number of iterations, etc. Moreover, mixed finite elements of Lagrange type, i.e.\ CG1 and CG2 elements, are adopted to discretize the corresponding mixed function space that includes the temperature and the displacement fields. For the damage minimization problem~\eqref{eq:staggered 2}, linear Lagrange finite elements (i.e.\ CG1 elements) and the \emph{PETScTAOSolver} of DOLFIN are utilized. Note that the latter allows to numerically enforce the damage irreversibility condition, i.e.\ $\dot{\alpha} \geq 0$, by uploading its lower bound at each endpoint of increment, such that $\alpha \in [\alpha_{\min}, 1]$ with $\alpha_{\min} \leftarrow \alpha_{n+1}$. The whole numerical procedure is shown in Algorithm~\ref{algo:algorithm total energy minimization}.

\begin{algorithm}[h!]
\caption{Semi-staggered algorithm proposed and implemented for the model in $\bracks{t_0, t_{f}}$}\label{algo:algorithm total energy minimization}
\KwIn{$\braces{\boldsymbol{u}_{0}, \, {\alpha}_{0}, \, T_{0}}$ at $t_0$}
\KwOut{$\braces{\boldsymbol{u}_{n}, {\alpha}_{n}, T_{n}}$ at $t=t_1, \dots, t_{n}, t_{n+1}, \dots, t_N=t_{f}$} Set initial and boundary conditions;

Initialize $t = t_0$; set final time $t_{f}$, absolute tolerance ${\err}_{\tol}$ and calculate the constant time increment $\Delta t = \frac{t_f - t_0}{N}$; 

Initialize the increment counter with $n =0$;

\While{t $\leq t_{f}$}
{
$t = t_{n+1} = t_0 + (n+1) \Delta t$;

Initialize an iteration counter with $i =0$; set $(\bm{u}_{n+1}^0, \alpha_{n+1}^0, T_{n+1}^0) \coloneqq (\bm{u}_n, \alpha_n, T_{n})$;

\Repeat{$\norm{\boldsymbol{u}_{n+1}^{i+1}-\boldsymbol{u}_{n+1}^{i}}_{2} < {\err}_{\tol}$ \textbf{and} $\norm{T_{n+1}^{i+1}-T_{n+1}^{i}}_{2} < {\err}_{\tol}$ \textbf{and} $\norm{\boldsymbol{\alpha}_{n+1}^{i+1}-\boldsymbol{\alpha}_{n+1}^{i}}_{2} < {\err}_{\tol}$}
{
Compute $\parens{\boldsymbol{u}_{n+1}^{i+1}, \, T_{n+1}^{i+1}}$ with $\parens{\boldsymbol{u}_{n+1}^{i}, \, T_{n+1}^{i}}$ by solving eq.~\eqref{eq:staggered 1}

Compute ${\alpha}_{n+1}^{i+1}$ with ${\alpha}_{n+1}^{i}$ from minimization of eq.~\eqref{eq:staggered 2}

$i=i+1$;
}

update $\parens{\bm{u}_n, \alpha_n, T_n} \gets \parens{\bm{u}_{n+1}, \alpha_{n+1}, T_{n+1}} $;

update the lower bound of phase field $\alpha_{\min} = \alpha_{n+1}$.
}
\end{algorithm}

\subsection{Application to plate cracking under thermal shocks}

This section is devoted to the numerical simulation of plate cracking under thermal shock conditions. For application purpose, we consider, as classically $\tq{C} (\alpha) = g(\alpha) \tq{C}^0$ with $g(\alpha)$ the degradation function that is taken as $g(\alpha) = (1-\alpha)^2$ and $\tq{C}^0$ the stiffness tensor of the sound solid. {Having in hand this expression, the dilatation tensor can be computed as $\bm{\beta} (\alpha) = \tq{C} (\alpha) : \bm{a}$ with $\bm{a}$ the constant tensor of linear thermal expansion (thermal strain tensor). The heat capacity $c$ is also taken constant.} As studied in~\cite{li2023} (see also Figure~\ref{thermal_plaque}), we consider a rectangular plate with a length $L=\qty{25}{mm}$ and a height $H=\qty{9.8}{mm}$. The plate is initially at a temperature, denoted by $T_0$, whose value depends on the specific test conditions. The right edge of the plate is constrained with zero horizontal displacement to account for the symmetry of the problem, while the other boundaries remain free of stress. A temperature $T_B=\qty{300}{K}$ is applied to these free boundaries to induce the thermal shock process with a temperature drop, such that $\Delta T = T_0 - T_B$.

The material parameters are summarized in Table~\ref{table:parameters} which can also be found in~\cite{li2023}.

\begin{figure}[h!]
\includegraphics[width=0.6\textwidth]{Plaque_thermique_second}
\caption{Geometry and boundary conditions of the quenching test on ceramics \cite{shao2011}.}\label{thermal_plaque}
\end{figure}

\begin{table}[htbp]
\caption{Material parameters for numerical quenching test reported by~\cite{li2023}.}\label{table:parameters}
\begin{tabular}{clcc}
\hline
\\[-2.7ex]
\hline
	Parameter 	& Name 							& Value				& Unit\\
\hline
\\[-2.7ex]
\hline
	$E$ 			& Young's modulus 				& $340$ 				& \unit{GPa}
\\ 	$\nu$ 		& Poisson coefficient 			& $0,22$ 			& --
\\	$G_c$ 		& Critical fracture energy 		& $42.47$ 			& \unit{J / m^2}
\\	$\kappa$		& Thermal conductivity 			& 300 				& \unit{W(m~K)}
\\	$c$			& Specific heat capacity 		& $0.775$ 			& \unit{J/(kg~K)}
\\	$a$ 			& Thermal expansion coefficient 	& $8 \cdot 10^{-6}$	& \unit{K^{-1}}
\\	$l_{0}$ 		& Internal length of damage 		& $0.092$			& \unit{mm}
\\
\hline
\end{tabular}
\end{table}

\subsubsection{AT1-type gradient damage model applied to the quenching test}

Numerical simulations are carried out by employing the AT1-type damage model. This model is defined by a linear damage function $\omega(\alpha) = \alpha$, and an elastic energy density given by $w_1 = \frac{3G_c}{8 l_0}$. {To investigate the impact of thermal loading on crack nucleation and propagation,} we perform numerical simulations for two initial temperatures $T_0 = \qty{550}{K}$ and $\qty{880}{K}$, corresponding to $\Delta T = \qty{250}{K}$ and $\Delta T = \qty{580}{K}$ on the stress free boundaries, respectively.

Figure~\ref{fig:contour 250} shows the temperature field $T$ as well as damage field $\alpha$ at times $t=\qty{1}{\mu s}$ and $t=\qty{10}{\mu s}$ for $\Delta T=\qty{250}{K}$. It is observed that the temperature field evolves more rapidly near the boundaries of the plate due to a high thermal gradient. This phenomenon leads to a progressive cooling of the plate, starting at the thermal shock boundaries and propagating inward. Concerning cracking (represented by strongly localized damage zone with $\alpha = 1$), it occurs at the boundaries subjected to thermal shock and logically propagates inward. Smaller cracks appear between larger ones, and their propagation is hindered by the presence of the latter.

Figure~\ref{fig:contour 580} illustrates the different fields at $t=\qty{10}{\mu s}$ for another temperature variation $\Delta T = \qty{580}{K}$.

By comparing with Figure~\ref{fig:contour 250}, the localized damage zone is shown to be more pronounced. It is also observed that cracks, particularly smaller cracks, become more numerous when the temperature variation $\Delta T = \qty{580}{K}$. This is due to a more intense thermal effect at this time $t = \qty{10}{\mu s}$ in the plate (see Figure~\ref{fig:contour 580}\eqref{t_10_temp_580}). Moreover, we show in Figure~\ref{fig:contour u_x} the contour of the horizontal displacement, denoted by $u_x$, at $t = \qty{10}{\mu s}$ respectively for $\Delta T = \qty{250}{K}$ and $\qty{580}{K}$. These horizontal displacement fields show jump at the location of vertical cracks. Note that the amplitude of the displacement is different for the two thermal shock conditions.

Finally, it is worth noting that at a sufficiently early time (e.g., $t = \qty{0.01}{\mu s}$), no crack appears under the thermal shock load of $\Delta T = \qty{250}{K}$. In contrast, at the same time, cracks are already nucleated under the higher thermal load $\Delta T = \qty{580}{K}$, which induce displacement discontinuities on the two sides of them. Readers are referred to Appendix~\ref{sec:complement figures} for more details about this effect.

\begin{figure}[h!]
\begin{subfigure}{0.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_TEMPERATURE_DELTA_T_250_1_us}
\caption{$T$ at $t=\qty{1}{\mu s}$}\label{t_1_temp_250}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_TEMPERATURE_DELTA_T_250_10_us}
\caption{$T$ at $t=\qty{10}{\mu s}$}\label{t_10_temp_250}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{0.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_FISSURATION_DELTA_T_250_1_us}
\caption{$\alpha$ at $t=\qty{1}{\mu s}$}\label{d_1_temp_250}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_FISSURATION_DELTA_T_250_10_us}
\caption{$\alpha$ at $t=\qty{10}{\mu s}$}\label{d_10_temp_250}
\end{subfigure}
\caption{Contours of temperature $T$ and damage fields $\alpha$ at $t=\qty{1}{\mu s}$ and $\qty{10}{\mu s}$ with $\Delta T = \qty{250}{K}$.}\label{fig:contour 250}
\end{figure}

\begin{figure}[h!]
\begin{subfigure}{0.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_TEMPERATURE_DELTA_T_580_10_us}
\caption{$T$ at $t=\qty{10}{\mu s}$}\label{t_10_temp_580}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_FISSURATION_DELTA_T_580_10_us}
\caption{$\alpha$ at $t=\qty{10}{\mu s}$}\label{d_10_temp_580}
\end{subfigure}
\caption{Contours of temperature $T$ and damage fields $\alpha$ at $\qty{10}{\mu s}$ with $\Delta T = \qty{580}{K}$.}\label{fig:contour 580}
\end{figure}

\begin{figure}[h!]
\begin{subfigure}{0.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_horizontal_diplacement_DELTA_T_250_10_us}
\caption{$u_x$ at $t=\qty{10}{\mu s}$ with $\Delta T = \qty{250}{K}$}\label{ux_temp_250}
\end{subfigure}
\hfill
\begin{subfigure}{0.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_horizontal_diplacement_DELTA_T_580_10_us}
\caption{$u_x$ at $t=\qty{10}{\mu s}$ with $\Delta T = \qty{580}{K}$}\label{ux_temp_580}
\end{subfigure}
\caption{Contours of horizontal displacement (denoted by $u_x$) at $\qty{10}{\mu s}$ with $\Delta T = \qty{250}{K}$ and $\qty{580}{K}$, respectively.}\label{fig:contour u_x}
\end{figure}

\subsubsection{Comparison to experimental observations}

In this section, for comparison purpose, we present numerical predictions by the developed model to experimental data reported by~\cite{shao2011}. To this end, we adopt the AT2 model (dissipation function $\omega(\alpha) = \alpha^2$ and $\omega_1 = \frac{G_c}{2\ell_0}$). This choice is justified only by the fact that the thermoelastic parameters considered for the simulations are taken from~\cite{li2023}, who also used an AT2-type model to compare their numerical results with the same experimental observations.

The configuration chosen is the same as that previously considered with the AT1 model. The temperature imposed on the edges is: $T_{B} = \qty{300}{K}$. The initial temperatures are set to $T_0 = \qty{550}{K} $, $\qty{680}{K}$ and $\qty{880}{K}$. These values respectively correspond to temperature variations of $\Delta T = \qty{250}{K}$, $\qty{380}{K}$ and $\qty{580}{K}$ on the hole boundary subjected to thermal shock.

Figures~\ref{QS_plaque_thermique_comp1}, \ref{QS_plaque_thermique_comp2}, and~\ref{QS_plaque_thermique_comp3} illustrate the evolution of the damage pattern $\alpha$ at $t = \qty{10}{\mu s}$ for different thermal shock amplitudes, as well as the comparison with the experimental results reported by~\cite{shao2011}. The contours of the predicted damage field show good qualitative agreement with the experimental observations. This confirms the ability of the proposed model to realistically reproduce the cracking phenomena induced by transient thermal stress. As previously observed for the AT1 model, cracks develop nucleate from the boundary exposed to thermal shock and propagate towards the interior of the plate. A gradual increase in crack density is also observed as the amplitude of the thermal shock is high, which is in agreement with the experimental results. During the propagation process, numerous short and parallel cracks appear, among which some stop while others continue to propagate. Finally, it can be observed from Figure~\ref{fig:contour_alpha_AT2} that increasing $\Delta T$ results in a larger number of cracks.

\begin{figure}[!h]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[width=\linewidth, height=4.5cm]{New_AT2_FISSURATION_DELTA_T_250_10_us}
\caption{$\alpha$ at $t=\qty{10}{\mu s}$}\label{QS_plaque_thermique_comp1_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[width=\linewidth, height=4.5cm]{AT2_experience_delta_T_250}
\caption{$\alpha$ at $t=\qty{10}{\mu s}$ (experimental result)}\label{QS_plaque_thermique_comp1_b}
\end{subfigure}
\caption{Predicted damage field $\alpha$~\eqref{QS_plaque_thermique_comp1_a} and the corresponding experimental result~\eqref{QS_plaque_thermique_comp1_b} at $\qty{10}{\mu s}$ with $\Delta T = \qty{250}{K}$.}\label{QS_plaque_thermique_comp1}
\end{figure}

\begin{figure}[!h]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[width=\linewidth, height=4.5cm]{New_AT2_FISSURATION_DELTA_T_380_10_us}
\caption{$\alpha$ at $t=\qty{10}{\mu s}$}\label{QS_plaque_thermique_comp2_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[width=\linewidth, height=4.5cm]{AT2_experience_delta_T_380}
\caption{$\alpha$ at $t=\qty{10}{\mu s}$ (experimental result)}\label{QS_plaque_thermique_comp2_b}
\end{subfigure}
\caption{Predicted damage field $\alpha$~\eqref{QS_plaque_thermique_comp2_a} and the corresponding experimental result~\eqref{QS_plaque_thermique_comp2_b} at $\qty{10}{\mu s}$ with $\Delta T = \qty{380}{K}$.}\label{QS_plaque_thermique_comp2}
\end{figure}

\begin{figure}[!h]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[width=\linewidth, height=4.5cm]{New_AT2_FISSURATION_DELTA_T_580_10_us}
\caption{$\alpha$ at $t=\qty{10}{\mu s}$}\label{QS_plaque_thermique_comp3_a}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\includegraphics[width=\linewidth, height=4.5cm]{AT2_experience_delta_T_580}
\caption{$\alpha$ at $t=\qty{10}{\mu s}$ (experimental result)}\label{QS_plaque_thermique_comp3_b}
\end{subfigure}
\caption{Predicted damage field $\alpha$~\eqref{QS_plaque_thermique_comp3_a} and the corresponding experimental result~\eqref{QS_plaque_thermique_comp3_b} at $\qty{10}{\mu s}$ with $\Delta T = \qty{580}{K}$.}\label{QS_plaque_thermique_comp3}
\end{figure}

\begin{figure}[h!]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=4.5cm]{AT2_FISSURATION_DELTA_T_250_1_us}
\caption{$\alpha$ at $t=\qty{1}{\mu s}$ with $\Delta T = \qty{250}{K}$}\label{alpha_temp_250}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\includegraphics[width=1\textwidth,height=4.5cm]{AT2_FISSURATION_DELTA_T_580_1_us}
\caption{$\alpha$ at $t=\qty{1}{\mu s}$ with $\Delta T = \qty{580}{K}$}\label{alpha_temp_580}
\end{subfigure}
\caption{Predicted damage field $\alpha$ for AT2 at $\qty{1}{\mu s}$ with $\Delta T = \qty{250}{K}$ and $\qty{580}{K}$, respectively.}\label{fig:contour_alpha_AT2}
\end{figure}

\section{Conclusion}

In this work, we developed an incremental variational approach for gradient damage coupled with thermoelasticity and heat conduction. This was formulated by relying on the Generalized Standard Materials framework with a particular account of the non local damage. The proposed incremental variational principle consists in the minimization of a four-fields incremental energy functional, in which the consideration of a \emph{kinetic entropy} is required for an appropriate description of the coupling between the heat conduction and the thermoelastic deformations. More specifically, the irreversible process is quantified via a dissipation potential that is expressed in a joint form between the damage dissipation and the heat conduction. By simultaneously considering a suitable dissipation potential together with a Helmholtz free energy we succeed to establish a total energy functional. The incremental form of the latter is dependent on the displacement, the damage, the absolute temperature fields as well as on the \emph{kinetic entropy}. It was then demonstrated that the corresponding weak forms lead to the appropriate Euler--Lagrange equations and boundary conditions for the concerned coupled system. An incremental variational principle was hence proposed. With a suitable kinetic entropy approximation, the variational procedure is reduced to a three-field dependency: displacement, damage, and temperature fields.

This has led to the variational formulation of a complete thermoelastic-gradient damage model, whose numerical implementation has been carried out through a semi-staggered algorithm by using the open-source platform FEniCS. For illustration purpose, the model has been applied to the simulate the fracturing process of a plate under thermal shocks. A comparison with available experimental data of~\cite{shao2011} has shown qualitatively good agreements.

Finally, it should be noted that some extensions/improvements of the proposed variational approach are possible. We are particularly interested by:
\begin{itemize}
\item consideration of a strong coupling between the thermoelasticity and the gradient damage: this could be achieved by introducing an integration factor (proposed by~\cite{Stainier2006}) in the dissipation potential of gradient damage in order to well incorporate its contribution to entropy production;
\item an extension of the present study to the context thermo-poroelastic couplings for which a thermodynamic basis is already available (see for instance~\cite{Coussy2004}): note that this type of coupling has been recently investigated by~\cite{liu2024} in the context of phase-fields methods.
\end{itemize}
Another interesting point could consist in comparing predictions of fracture nucleation under thermal shocks to that which can be established by considering the coupled criterion proposed by D.~Leguillon and coauthors~\cite{leguillon2013, ricardo2020} (see several papers in the present issue).

\appendix

\section{Thermodynamics formulation for the coupling between gradient damage and thermoelasticity}\label{sec:thermodynamics appendix}

Thermodynamics provides a theoretical and essential framework allowing to rigorously describe the evolution of a multiphysical system depending on all forms of energy. In this section, we formulate a suitable thermodynamic description of thermoelasticity with heat conduction coupled to gradient damage. Given the continuum medium $\Omega$, the global internal energy $\mathcal{E}$, the global entropy $\mathcal{S}$ and the calorific power are classically expressed as:
\begin{equation}\label{eq: global form E and S}
\mathcal{E} = \int_{\Omega} e \, \dd \Omega;
\qquad
\mathcal{S} = \int_{\Omega} s \, \dd \Omega;
\qquad
\mathcal{P}_{\cal} =
\int_{\Omega} {r} \, \dd \Omega - \int_{\partial \Omega} {\bm{q}} \cdot {\bm{n}} \, \dd S
\end{equation}
where $e$ is the density of internal energy per unit volume.

Introducing eqs.~\eqref{eq: global form E and S} and~\eqref{eq:P_int} into~\eqref{eq:global form thermo} yields the following local energy balance equation and the Clausius--Duhem inequality for a thermomechanical medium with gradient damage effects:
\begin{equation}\label{eq:local_form_damage}
\begin{aligned}
	& \dot{e} = \bm{\sigma} : \dot{\bm{\varepsilon}} + r - \divergence{\bm{q}} + {Y}_{\alpha} \dot{\alpha} + {\bm{Y}}_{\bmnabla {\alpha}} \cdot \bmnabla \dot{{\alpha}},
\\	& {\mathcal{D}} = \bm{\sigma} : \dot{\bm{\varepsilon}} - \dot{e} + {T} \dot{s} + {Y}_{\alpha} \dot{\alpha} + {\bm{Y}}_{\bmnabla {\alpha}} \cdot \bmnabla \dot{{\alpha}} - \frac{\bm{\nabla} T}{T} \cdot \bm{q} \geq 0.
\end{aligned}
\end{equation}
Yet, since the entropy $s$ is less ``controllable'', it would be more appropriate to adopt its conjugate variable, i.e.\ the absolute temperature $T$, as an external state variable. In this context, it can be preferable to consider the Helmholtz free energy per unit volume $w$, which satisfies the following partial Legendre transform:
\begin{equation}\label{eq:legendre w}
-w = Ts-e.
\end{equation}
It is convenient to emphasize here that $e$ is convex with respect to all of its arguments. It follows by duality that $w$ is a concave function of the temperature field $T$, while its convexity with respect to other variables still holds. Substituting eq.~\eqref{eq:legendre w} into~\eqref{eq:local_form_damage} yields eq.~\eqref{eq:C-D inequality in w}, namely:
\[
{\mathcal{D}} =
\underbrace{\vphantom{\frac{T}{T}}\bm{\sigma} : \dot{\bm{\varepsilon}} - \dot{w} - s \dot{T} + {Y}_{\alpha} \dot{\alpha} + {\bm{Y}}_{\bmnabla {\alpha}} \cdot \bmnabla \dot{{\alpha}}}_{{\mathcal{D}}_{\intr}}
\underbrace{- \frac{\bm{\nabla} T}{T} \cdot \bm{q}}_{{\mathcal{D}}_{\therm}}
\geq 0
\]
where ${\mathcal{D}}_{\intr}$ is the intrinsic dissipation and ${\mathcal{D}}_{\therm}$ the dissipation of the heat conduction.

Moreover, through a straightforward rearrangement of the expression of ${\mathcal{D}}_{\intr}$, the following equality can be readily obtained:
\begin{equation}\label{eq:discussion w}
\dot{w} =
\bracks[\big]{
\parens{\bm{\sigma} : \dot{\bm{\varepsilon}} + Y_{\alpha} \dot{\alpha} + \bm{Y}_{\bm{\nabla}{\alpha}} \cdot \bmnabla \dot{\alpha}
} - \mathcal{D}_{\intr}
} - s \dot{T}.
\end{equation}
From a similar perspective as presented in~\cite{Maitournam2013} in the case without gradient damage, the above equation reveals that the variation of $w$ can be divided into two parts: the first one in the bracket corresponds to the isothermal and non-dissipative mechanical work ; the other, i.e.\ $- {s} \dot{T}$, accounts for variations of $w$ when all other parameters hold constant and is obviously associated with the temperature variation. Hence, one can define $w$ as a function of $\bm{\varepsilon}$, $\alpha$, $\bmnabla \alpha$ and $T$, which has been given by eq.~\eqref{eq:Helmholz bis}, namely:
\[%\label{sec:state variables and laws}
w \coloneqq w (\bm{\varepsilon}, T, \alpha, \bm{\nabla} \alpha).
\]
Since the state variables are independent, eq.~\eqref{eq:discussion w} can be reformulated as:
\begin{equation}\label{eq:identification forces rev01}
\dot{w} =
\frac{\partial w}{\partial \bm{\varepsilon}} : \dot{\bm{\varepsilon}} + \frac{\partial w}{\partial {\alpha}} \dot{\alpha} + \frac{\partial w}{\partial {\bmnabla \alpha}} \cdot \bmnabla \dot{\alpha} + \frac{\partial w}{\partial T} \dot{T} =
\bracks[\big]{
\parens{\bm{\sigma} : \dot{\bm{\varepsilon}} + Y_{\alpha} \dot{\alpha} + \bm{Y}_{\bm{\nabla}{\alpha}} \cdot \bmnabla \dot{\alpha}
} - \mathcal{D}_{\intr}
} - s \dot{T}
\end{equation}
which allows to identify the following non-dissipative thermodynamic forces (denoted with superscript ``nd'') as the state laws given by eq.~\eqref{eq:state eqs with w bis}. Inserting the latter into~\eqref{eq:C-D inequality in w} yields the local dissipation ${\mathcal{D}}$ described by eq.~\eqref{eq:C-D inequality in w simplified}.

Finally, as mentioned in Section~\ref{sec: general variational principle}, the evolution laws are described through a joint dissipation potential $\varphi$, i.e.\ eq.~\eqref{eq:dissipation potential 1}. It is worthy to point out that the in the dissipation potential, $\varphi^{M}$ is supposed as a function of the rate of the internal state variables (i.e., $\dot{\alpha}$ and $\bmnabla \dot{\alpha}$) and parametrically depends on its current state variables. For $\varphi^{T}$, and in link to the definition of the thermal dissipation ${\mathcal{D}}_{\therm}$ (see eq.~\eqref{eq:C-D inequality in w simplified}), the simplest form could be the Fourier potential, which is defined as function of an \emph{additional} variable $- {\bmnabla T} / {T} $ and parametrically depends on $T$.

\section{Variational analysis of the incremental total energy functional \texorpdfstring{$I_n$}{I\_n}}\label{sec:variational calculation}

Having in hand eq.~\eqref{eq:total energy functional I_0}, along with the integral of the exact differential of $\dot{w}$ and the substitution eqs.~\eqref{eq:Helmholz bis}, \eqref{eq:dissipation potential 1}, and~\eqref{eq:external work} into~\eqref{eq:total energy functional I_0}, we obtain:
\begin{multline}\label{eq:total energy functional I_1}
I_n(\boldsymbol{\bm{u}}_{n+1}, T_{n+1}, \hat{s}_{n+1}, \alpha_{n+1}) =
\int_{\Omega} {w}_{n+1} (\bm{u}_{n+1}, {T}_{n+1}, {\alpha}_{n+1}, \bmnabla {\alpha}_{n+1}) - {w}_{n} (\bm{u}_{n}, {T}_{n}, {\alpha}_{n}, \bmnabla {\alpha}_{n}) \, \dd \Omega
\\	+ \int_{t_n}^{t_{n+1}} \int_{\Omega} \hat{s} \dot{T} + {\varphi}\parens[\Bigg]{\dot{\alpha}, \bm{\nabla} \dot{\alpha}, - \frac{\bmnabla {T}}{T}; \alpha, \bm{\nabla} \alpha} \, \dd \Omega \, \dd t - \int_{t_n}^{t_{n+1}} \int_{\partial_{Q} \Omega} Q \, \log \frac{T}{T_0} \, \dd S \, \dd t
\\	- \int_{\Omega}
\bm{F}_{n+1} \cdot \parens{\bm{u}_{n+1} - \bm{u}_{n}} \, \dd \Omega - \int_{\partial_{\bm{T}} \Omega} \bm{T}_{n+1} \cdot \parens{\bm{u}_{n+1} - \bm{u}_{n}} \, \dd S \, \dd t.
\end{multline}
Taking variation of eq.~\eqref{eq:total energy functional I_1} with respect to $T$ and simultaneously considering~\eqref{eq:state eqs with w bis} and~\eqref{eq:thermo forces d} gives the related weak form:
\begin{equation}\label{eq:first variation p 1}
- \int_{\Omega} s_{n+1} \, \delta {T}_{n+1} \dd \Omega +
\int_{t_{n}}^{t_{n+1}} \int_{\Omega}
\parens[\Bigg]{\hat{s} \, \delta \dot{T} + \frac{\bm{q}}{T} \cdot \bmnabla\delta T} \, \dd \Omega \, \dd t -\int_{t_{n}}^{t_{n+1}}\int_{\partial_{Q} \Omega} Q \frac{\delta {T}}{T} \, \dd S \, \dd t = 0
\end{equation}
whose integration by parts with respect to time and space yields:
\begin{multline}\label{eq:first variation p 2}
\int_{\Omega}
\parens{- {s}_{n+1} + \hat{s}_{n+1}} \delta T_{n+1} \dd \Omega - \int_{t_n}^{t_{n+1}}\int_{\Omega}
\bracks[\Bigg]{
\dot{\hat{s}} + \divergence{\frac{\bm{q}}{T}}
} \delta T \, \dd \Omega \, \dd t
\\	+ \int_{t_{n}}^{t_{n+1}} \int_{\partial_{Q} \Omega}
\parens{{\bm{q}} \cdot \bm{n} - {Q}
}
\frac{\delta T}{T}
\dd S \, \dd t =0.
\end{multline}
Ensuring stationarity yields at $t = t_{n+1}$:
\begin{equation}\label{eq: case phi hat = phi n first method}
\hat{s}_{n+1} = {s}_{n+1} ; \qquad \forall \bm{x} \in \Omega
\end{equation}
and in the time interval $[t_n, t_{n+1}]$
\begin{align}
	& \dot{\hat{s}} + \divergence{\frac{\bm{q}}{T}} = {0}; & \quad	& \forall \bm{x} \in \Omega, \label{eq:balance equation for T incr}
\\	\hspace{4cm} & \bm{q} \cdot \bm{n} = Q; &	\quad & \forall \bm{x} \in \partial_{Q} \Omega. \hspace{4cm} \label{eq:boundary condition of T incr}
\end{align}
It follows from eq.~\eqref{eq: case phi hat = phi n first method} that the \emph{kinetic entropy} $\hat{s}_{n+1}$ is optimized at time $t_{n+1}$ as equal to the \emph{equilibrium entropy} ${s}_{n+1}$, which indicates that the material point reaches the local thermodynamic equilibrium at $t_{n+1}$. Also, eq.~\eqref{eq:balance equation for T incr} reveals that the evolution of \emph{kinetic entropy} satisfies the heat equation in $[t_n, t_{n+1}]$ with the corresponding Neumann boundary condition on $ \partial_{Q} \Omega$ given by~\eqref{eq:boundary condition of T incr}.

Similarly, by taking variations of eq.~\eqref{eq:total energy functional I_1} with respect to the other variables and considering~\eqref{eq:state eqs with w bis} and~\eqref{eq:thermo forces d}, the corresponding weak forms for $\bm{u}$, $\alpha$ and $\hat{s}$ can be derived as follows:
\begin{align}
& -\int_{\Omega}
\bracks[\big]{
\divergence{\bm{\sigma}_{n+1}} +
\bm{F}_{n+1}
} \cdot \delta\boldsymbol{u}_{n+1} \, \dd \Omega + \int_{\partial_{\bm{T}} \Omega}
\parens[\big]{
\bm{\sigma}_{n+1} \cdot \boldsymbol{n} - \bm{T}_{n+1}
}
\cdot \delta\boldsymbol{u}_{n+1} \, \dd S =0, \label{eq: first variation u 2}
\\
& \begin{aligned}\label{eq:variation alpha}
	& \int_{\Omega}
\bracks[\big]{Y_{\alpha} \rvert_{n+1} - \divergence{\bm{Y}_{\bmnabla \alpha} \rvert_{n+1}}}
\delta \alpha_{n+1} \, \dd \Omega +\int_{\partial\Omega}
\parens[\big]{
\bm{Y}_{\bmnabla \alpha} \rvert_{n+1} \cdot \bm{n}} \delta \alpha_{n+1} \, \dd S
\\	& \quad +\int^{t_{n+1}}_{t_{n}} \int_{\Omega}
\braces[\Bigg]{
\bracks[\Bigg]{
\frac{\partial \varphi}{\partial \alpha} - \frac{\dd}{\dd t} \parens[\Bigg]{\frac{\partial \varphi}{\partial \dot{\alpha}}}
} -\divop
\bracks[\Bigg]{
\frac{\partial \varphi}{\partial \bmnabla {\alpha}} - \frac{\dd}{\dd t} \parens[\Bigg]{\frac{\partial \varphi}{\partial \bmnabla \dot{\alpha}}}
}} \delta \alpha \, \dd \Omega \, \dd t
\\	& \quad + \int^{t_{n+1}}_{t_{n}} \int_{\partial \Omega}
\braces[\Bigg]{
\bracks[\Bigg]{
\frac{\partial \varphi}{\partial \bmnabla {\alpha}} - \frac{\dd}{\dd t} \parens[\Bigg]{\frac{\partial \varphi}{\partial \bmnabla \dot{\alpha}}}
} \cdot \bm{n}} \delta \alpha \, \dd S \, \dd t = 0,
\end{aligned}
\\
& \int^{t_{n+1}}_{t_{n}} \int_{\Omega} \dot{T} \, \delta \hat{s} \, \dd \Omega \, \dd t = 0, \label{eq:variation phi hat}
\end{align}
where $\cdot\,\rvert_{n+1}$ symbolizes the estimate of considered quantity at time $t_{n+1}$.

\section{Damage and displacement fields at \texorpdfstring{$t = 0.0\qty{1}{\mu s}$}{t = 0.01 mu s} and \texorpdfstring{$0.\qty{1}{\mu s}$}{0.1 mu s} with \texorpdfstring{$\Delta T = \qty{250}{K}$}{Delta T = 250 K}}\label{sec:complement figures}

In this section, we display the damage field (i.e.\ $\alpha$) as well as the horizontal and vertical displacements (i.e.\ $u_x$ and $u_y$) fields respectively at $t = 0.0\qty{1}{\mu s}$ and $0.\qty{1}{\mu s}$ under the thermal shock loads $\Delta T = \qty{250}{K}$. The primary objective is to show the transition from damage processes at early stage to occurence cracks. It can be observed in Figure~\ref{fig:contour annexe}\eqref{d_001_temp_250} that at $t = 0.0\qty{1}{\mu s}$, damage starts on the thermal loading boundaries without any strongly localized field (i.e.\ smeared cracks) such that $\alpha < 1$. While at $t= 0.\qty{1}{\mu s}$ (see Figure~\ref{fig:contour annexe}\eqref{d_01_temp_250}), the damage has already been localized leading to the propagation of the diffused cracks. In addition, it is very interesting to note that displacement jumps are observed (in Figures~\ref{fig:contour annexe}\eqref{ux_001_temp_250}--\eqref{uy_01_temp_250}) at $t = 0.\qty{1}{\mu s}$ for horizontal (resp.\ vertical) displacements on the top and bottom (resp.\ left) boundaries, which somehow represent the opening of these smeared cracks during the damage evolution.

\begin{figure}[h!]
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_FISSURATION_DELTA_T_250_0.01_us}
\caption{$\alpha$ at $t= 0.0\qty{1}{\mu s}$}\label{d_001_temp_250}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_FISSURATION_DELTA_T_250_0.1_us}
\caption{$\alpha$ at $t= 0.\qty{1}{\mu s}$}\label{d_01_temp_250}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_horizontal_displacement_DELTA_T_250_0.01_us}
\caption{$u_x$ at $t= 0.0\qty{1}{\mu s}$}\label{ux_001_temp_250}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_horizontal_displacement_DELTA_T_250_0.1_us}
\caption{$u_x$ at $t= 0.\qty{1}{\mu s}$}\label{ux_01_temp_250}
\end{subfigure}
\vskip .5\baselineskip
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_vertical_displacement_DELTA_T_250_0.01_us}
\caption{$u_y$ at $t= 0.0\qty{1}{\mu s}$}\label{uy_001_temp_250}
\end{subfigure}
\hfill
\begin{subfigure}{.48\linewidth}
\centering
\includegraphics[width=1\textwidth,height=5cm]{AT1_vertical_displacement_DELTA_T_250_0.1_us}
\caption{$u_y$ at $t= 0\qty{1}{\mu s}$}\label{uy_01_temp_250}
\end{subfigure}
\caption{Contours of damage $\alpha$, horizontal displacement $u_x$ and vertical displacement $u_y$ and at $t=\qty{0.011}{\mu s}$ and $\qty{0.1}{\mu s}$ with $\Delta T = \qty{250}{K}$.}\label{fig:contour annexe}
\end{figure}

\section*{Declaration of interests}

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

\printbibliography

\end{document}