\makeatletter
\@ifundefined{HCode}
{\documentclass[screen, CRMECA, Unicode, published]{cedram}
\newenvironment{noXML}{}{}%
\usepackage[figuresright]{rotating}
\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\xtbody{\noalign{\relax}}
\def\tabnote#1{\vskip4pt\parbox{.51\linewidth}{#1}}
\def\xxtabnote#1{\vskip4pt\parbox{.65\linewidth}{#1}}
\def\mathbi#1{\boldsymbol{#1}}
\def\xmathop#1{\mathop{\text{#1}}}
\def\rmDelta{\Delta}
\usepackage{etoolbox}
\usepackage{multirow}
\usepackage{upgreek}
\usepackage{amssymb}
\csdef{Seqnsplit}{\\}
\def\botline{\\\hline}
\newcounter{runlevel}
\usepackage{multirow} 
\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\xmorerows#1#2{#2}
\let\MakeYrStrItalic\relax
\def\jobid{crmeca20230403}
\def\morerows#1#2{#2}
\def\sfrac#1#2{{#1}/{#2}}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}}
\def\sofrac#1#2{({#1}/{#2})}
\def\0{\phantom{0}}
\def\rmmu{\upmu}
\def\mn{\phantom{$-$}}
\def\refinput#1{}
\def\back#1{}
\newcommand{\bnabla}{\boldsymbol{\nabla}}
\newcommand{\bgamma}{\boldsymbol{\gamma}}
\newcommand{\bdelta}{\boldsymbol{\delta}}
\def\balpha{\boldsymbol{\alpha}}
\def\bfbeta{\boldsymbol{\beta}}
\def\bvarepsilon{\boldsymbol{\varepsilon}}
\def\bsigma{\boldsymbol{\sigma}}
\newenvironment{inftab}{}{}
\def\Custref#1#2{\hyperlink{#2}{#1}}
\def\og{\guillemotleft}
\def\fg{\guillemotright}
\makeatletter
\g@addto@macro{\UrlBreaks}{\UrlOrds}
\gappto{\UrlBreaks}{\UrlOrds}
\makeatother
\DOI{10.5802/crmeca.227}
\datereceived{2023-04-24}
\daterevised{2023-08-02}
\datererevised{2023-09-18}
\dateaccepted{2023-09-18}
\ItHasTeXPublished
}
{\documentclass[crmeca]{article}
\usepackage[T1]{fontenc} 
\usepackage{upgreek}
\def\CDRdoi{10.5802/crmeca.227}
\def\newline{\unskip\break}
\def\xmorerows#1#2{\morerows{#1}{#2}}
\def\sfrac#1#2{{#1}/{#2}}
\def\xxtabnote#1{\tabnote{#1}}
\def\xmathop#1{\mathop{{#1}}}
\def\xtbody{\tbody}
\def\hypertarget#1#2{\label{#1}}
\usepackage{stmaryrd}
\newenvironment{sidewaystable}{\begin{table*}}{\end{table*}}
}

\makeatother

\dateposted{2023-12-04}
\begin{document}


\begin{noXML}

\TopicEF{Research article}{Article de recherche}

\title{Exact expressions of the uniaxial perfectly elasto-plastic
stress wave and induced mechanical fields in the case of a finite
impact: application to laser shock peening}

\alttitle{Expressions exactes de l'onde de contrainte \'elasto-plastique
uniaxiale et des champs m\'ecaniques induits dans le cas d'un impact fini
: application au grenaillage laser}

\author{\firstname{Lucas} \lastname{Lapostolle}\CDRorcid{0000-0002-1511-9421}\IsCorresp}
\address{PIMM, Arts et Metiers Institute of Technology, CNRS, Cnam,
HESAM University, 151~boulevard de l'Hopital, 75013, Paris, France}
\email{lucas.lapostolle@ensam.eu}
\email[L. Lapostolle]{lucas.lapostolle@ensta-bretagne.org} 
                   
\author{\firstname{L\'eo} \lastname{Morin}\CDRorcid{0000-0002-6694-1212}}
\address{University of Bordeaux, CNRS, Arts et metiers Institute of
Technology, Bordeaux INP, I2M, Bordeaux, F-33400, Talence, France}
\email[L. Morin]{leo.morin@ensam.eu}

\author{\firstname{Katell} \lastname{Derrien}\CDRorcid{0000-0002-3564-2586}}
\address[1]{PIMM, Arts et Metiers Institute of Technology, CNRS, Cnam,
HESAM University, 151~boulevard de l'Hopital, 75013, Paris, France}
\email[K. Derrien]{katell.derrien@ensam.eu}

\author{\firstname{Laurent} \lastname{Berthe}\CDRorcid{0000-0002-9978-9843}}
\address[1]{PIMM, Arts et Metiers Institute of Technology, CNRS, Cnam,
HESAM University, 151~boulevard de l'Hopital, 75013, Paris, France}
\email[L. Berthe]{laurent.berthe@cnrs.fr}

\author{\firstname{Olivier} \lastname{Castelnau}\CDRorcid{0000-0001-7422-294X}}
\address[1]{PIMM, Arts et Metiers Institute of Technology, CNRS, Cnam,
HESAM University, 151~boulevard de l'Hopital, 75013, Paris, France}
\email[O. Castelnau]{olivier@castelnau@ensam.eu}

\keywords{\kwd{Laser shock peening}
\kwd{Elastic-plastic wave propagation}
\kwd{Numerical simulation}
\kwd{Analytical solution}
\kwd{Plastic strains}
\kwd{Residual stresses}}

\altkeywords{\kwd{Grenaillage laser}
\kwd{Propagation d'ondes \'elasto-plastiques}
\kwd{Simulation num\'erique}
\kwd{Solution analytique}
\kwd{D\'eformations plastiques}
\kwd{Contraintes r\'esiduelles}}

\begin{abstract} 
This paper aims at providing exact expressions for the mechanical
fields induced by Laser Shock Peening and comparing them to their
numerical estimations. We use a uniaxial strain field hypothesis with
an elastic perfectly plastic behavior to derive the stress wave
equation. An exact solution to this equation is given using the method
of characteristics for a step time profile for the pressure loading,
and numerically using finite differences schemes adapted for this
hyperbolic equation. An additional residual stress modeling is used,
providing the residual stress distribution assuming a planar infinite
plate with a finite thickness. Results are presented for three loading
pressures, each one corresponding to a different structure in the exact
solution. The exact and numerical results present a good match,
allowing either the use of the exact solution for an initial estimation
of the mechanical fields, or to test the accuracy of other numerical
methods.
\end{abstract}


\begin{altabstract} 
Cet article vise \`a fournir des expressions exactes pour les champs
m\'ecaniques induits par le grenaillage laser et \`a
les comparer \`a leurs estimations num\'eriques. Nous utilisons une
hypoth\`ese de champ de d\'eformation uniaxial avec un comportement
\'elastique parfaitement plastique pour obtenir l'\'equation
de l'onde de contrainte. Une solution exacte de cette \'equation est donn\'ee en
utilisant la m\'ethode des caract\'eristiques pour un profil temporel de
la charge de pression, et num\'eriquement en utilisant des sch\'emas de
diff\'erences finies adapt\'es \`a cette \'equation hyperbolique. Une
mod\'elisation suppl\'ementaire des contraintes r\'esiduelles est
utilis\'ee, fournissant la distribution des contraintes r\'esiduelles en
supposant une plaque infinie plane avec une \'epaisseur finie. Les
r\'esultats sont pr\'esent\'es pour trois pressions de chargement, chacune
correspondant \`a une structure diff\'erente dans la solution exacte. Les
r\'esultats exacts et num\'eriques pr\'esentent une bonne concordance, ce
qui permet soit d'utiliser la solution exacte pour une estimation
initiale des champs m\'ecaniques, soit de tester la pr\'ecision d'autres
m\'ethodes num\'eriques.
\end{altabstract}

\shortrunauthors

\thanks{Carnot Institute ARTS, French Agence Nationale de la Recherche
(ANR), research project ForgeLaser (Grant agreement no.
ANR-18-CE08-0026)}

\maketitle

{\vspace*{-6pt}}

\end{noXML}

\section{Introduction}\label{sec1}

\def\tn#1{\text{#1}}
\def\ds{\displaystyle}
\def\eqref#1{(\ref{#1})}
\def\aj{}
\def\sofrac#1#2{({#1}/{#2})}

This work tackles aspects of the modeling of the Laser Shock Peening
(LSP) process. For this process, a high-energy laser is used to impact
the surface of a specimen in order to introduce compressive residual
stresses, eventually improving its fatigue behavior. The energy used
(${\sim}$10~J) is applied on a laser spot with a diameter of a few
millimeters, with an impulsion of a dozen of nanoseconds, leading to
high power densities at the surface (${\sim}$10~GW{${\cdot}$}cm$^{-2}$) and
pressure (a few GPa). These pressures are enough to induce
heterogeneous plastic strain fields in the target, leading to residual
stresses. Their compressive nature have been related to improved
fatigue performances \cite{ding_laser_2006,peyre_laser_1996}.

Efforts have been made to model LSP for more than 30 years, starting
with the works of Ballard~\cite{ballard_contraintes_1991} (see also
\cite{ballard_residual_1991}) and Braisted~\cite{braisted_finite_1999} who
pioneered the modeling and simulation of LSP. Since then, numerous
works were done to evaluate the effects of the different parameters of
LSP including size and shape of the laser spot, pressure amplitude,
pulse duration, overlapping, etc.
\cite{julan_simulation_2013,peyre_fem_2007,peyre_fem_2003,xiang_numerical_2021,brockman_prediction_2012} 
using finite elements modeling. The effect of the sample microstructure
on the wave propagation and the residual stress field has been
investigated for 1D laminate specimens \cite{lapostolle_modeling_2022}
and 2D polycrystalline aggregates \cite{lapostolle_modeling_2023}, and
a numerical method for fast residual stress estimation has been
proposed recently using an elasto-viscoplastic relation
\cite{lapostolle_fast_2023}.

Beside the numerical simulations of LSP, an important aspect of the
work of Ballard~\cite{ballard_contraintes_1991} is to provide an exact
solution for the stress wave propagation under certain hypotheses.
Indeed, this work shows that for a uniaxial state of strain and an
elastic perfectly plastic material, it is possible to predict the
stress and velocity wave caused by a pressure step loading. Exact
solutions are particularly useful, because they can notably provide a
baseline to assess the quality of a numerical method. This was
exploited by Heuz\'e~\cite{heuze_laxwendroff_2017} who considers the case of a
finite impact (i.e. a loading with a given duration after which it
becomes zero), but using a bar modeling, i.e. a single  {non} vanishing
stress component (see also \cite{wang_foundations_2007}), which
provides different results than with a uniaxial strain modeling, as
considered here and better suited for  modelling the mechanical fields
at the center of the LSP impact.  {The structure of the exact solution
in the uniaxial strain case is introduced in the work of
Wang~\cite{wang_unloading_2004} for a specific range of pressure amplitudes
leading to complex interaction between elastic and plastic waves.}  In
addition, semi-analytical solutions also exist for the residual stress
profile, based on a heterogeneous eigenstrain distribution. The 
{eigenstrain} incompatibilities caused by this distribution generate
elastic strains for the structure to be at equilibrium, leading to
residual stresses. Such an approach was already widely used in the
literature
\cite{korsunsky_modelling_2005,gelineau_etude_2018,johns_thermal_1965,lapostolle_fast_2023,ahdad_contraintes_1996},
and particularly in the case of Almen strips deflection analysis
\cite{glaser_evaluation_2022}.

Several numerical methods can be used to simulate LSP. The most widely
used method is the finite elements method (see for example
\cite{braisted_finite_1999,brockman_prediction_2012,peyre_fem_2007}).
However other methods can be used for the problem of a stress wave
propagation, using notably a finite differences approach
\cite{leveque_finite_2002,heuze_laxwendroff_2017,wang_conservative_1993}. 
{The numerical procedures of these methods do not concern the same
formulation, and the wave equation with the displacement field as
unknown used in finite element solvers can be recast as a set of
conservation laws on the stress and velocity fields, for which
finite-differences based numerical methods can be used.} 

The aim of this work is to extend the work of
Ballard~\cite{ballard_contraintes_1991}, as well as reporting some aspects of
the exact solutions that were previously overlooked or not
investigated. In particular, we provide a detailed exact solution for
the elasto-plastic uniaxial stress/velocity wave and the various
reflection phenomenon caused by plasticity, but also for the plastic
strain distribution caused by the propagation, which has not been
previously reported in such detail. We use these results to give an
in-depth exact solution for the residual stress profile, and an exact
expression of the surface displacement caused by the laser shock. 
{Such an exact solution dedicated to laser shock applications can be
used to provide an initial sensitivity study of some laser parameters,
and thus helps optimize the number of experiments needed to observe
specific behaviors.} Moreover, by comparing the exact solution to the
various numerical methods, we also provide some insights as to which
ones are better suited to tackle this problem. The article is organized
as follows: Section~\ref{sec:mech_model} details the mechanical
modeling of the problem, the stress wave equation and its exact and
numerical solutions. Section~\ref{sec:residual_stresses} presents the
residual stresses modeling and the underlying hypotheses. 
Section~\ref{sec:results} presents results, both exact and numerical, for two
loading amplitudes, each leading to a different exact solution.
Finally, Section~\ref{sec:discussion} discusses the results, provides
some additional information for their interpretation  {and further
clarifies the position of the current work with respect to the cited
literature in light of the presented results}.


\vspace*{-2pt}
\section{Mechanical modeling of the stress wave
propagation}\label{sec:mech_model}
\vspace*{-2pt}

We are concerned with the propagation of elastic-plastic waves produced
by a laser shock. An illustration of the problem considered is
represented in Figure~\ref{fig:illustration}.


\begin{figure}[b!]
\includegraphics{fig01}
\vspace*{-3pt}
\caption{Illustration of the problem.\label{fig:illustration}}
\vspace*{-3pt}
\end{figure}

\subsection{Elasto-plastic stress wave equation}

\subsubsection{Material behavior} 
The problem of the elasto-plastic
shockwave propagation in small strains can be modeled using the
following set of equations:
{\begin{equation}
    \left\{
    \begin{array}{lr}
        {\tn{div}}\left( {\bsigma} \right)=\ds \rho \frac{\partial^2 \mathbi{u}}{\partial t^2} & \tn{Dynamic equilibrium} \vspace*{2pt} \\
        {\bsigma}=\mathbb{C}:{\bvarepsilon}_{e}=\mathbb{C}:({\bvarepsilon}-{\bvarepsilon}_{p}) & \tn{Elastic behavior} \vspace*{2pt}\\
        {\bvarepsilon}=\ds \frac{1}{2}\left(\nabla \mathbi{u}+\nabla^{T}\mathbi{u} \right) & \tn{Strain-displacement relation} \vspace*{2pt}\\
        f=\sigma_{\mathrm{vM}}\left({\bsigma}-\mathbi{X} \right)-R(p) \leq 0 & \tn{Plasticity criterion} \vspace*{2pt}\\
        \left\{
        \begin{array}{ll}
            \dot{{\bvarepsilon}}_{p}=\ds \dot{\lambda}\frac{\partial f}{\partial {\bsigma}} \vspace*{2pt}\\
            \dot{\lambda} \geq 0 \vspace*{2pt}\\
             {\dot{\lambda} f=0} 
        \end{array}
        \right. & \tn{Flow rule and consistency conditions} \vspace*{2pt}\\
        \left\{
        \begin{array}{ll}
            \sigma_{11}(x=0,t)=-P\left( H(t)-H(t-t_0) \right) \ \forall t \vspace*{2pt}\\
            \sigma_{11}(x=L,t)=0 \ \forall t \vspace*{2pt}\\
            \sigma_{11}(x,t=0)=0 \ \forall x.
        \end{array} 
        \right. & \tn{Boundary and initial conditions}
    \end{array}
    \right.
\label{eq:dynamics_eq}
\end{equation}}\noindent

In Equation \eqref{eq:dynamics_eq}, 
${\bsigma}$ is the Cauchy stress
tensor, $\mathbi{u}$ the displacement field, ${\bvarepsilon}$ the total
strain tensor, ${\bvarepsilon}_e$ the elastic strain tensor and
${\bvarepsilon}_p$ the plastic strain tensor. The material properties
are the stiffness tensor $\mathbb{C}$ and the density $\rho$. $f$
denotes the plasticity criterion, in which $\mathbi{X}$ is the backstress
used for kinematic hardening, $R$ the isotropic hardening function, $p$
the accumulated plastic strain and  {$\dot{\lambda}$} the plastic
multiplier. $\sigma_{\mathrm{vM}}$ is the von Mises equivalent stress, defined
using $J_2$, the second stress tensor invariant:
{\begin{equation}
\sigma_{\mathrm{vM}}=\sqrt{3J_2} \quad \tn{and} \quad J_2({\bsigma})=\ds 
\tfrac{1}{2}{\bsigma}_d : {\bsigma}_d,
\end{equation}}\noindent
where ${\bsigma}_d$ is the deviatoric part of the stress tensor
defined as 
${\bsigma}_d={\bsigma}-\sofrac{1}{3}\mathrm{Tr}({\bsigma})\mathbi{I}$,
$\mathbi{I}$ being the second-order identity tensor. The $\tn{Tr}(\cdot)$
notation corresponds to the Trace operator. In the boundary conditions
of Equation \eqref{eq:dynamics_eq}$_6$, the $H$ function is the
Heaviside step function, $P$ the strictly positive plasma pressure, and
$t_0$ a strictly positive impulsion time.  The material is at rest
initially, and not previously loaded, so that $\mathbi{v}(x,t=0)=\mathbf{0}$,
$\mathbi{u}(x,t=0)=\mathbf{0}$ ($\mathbi{v}=\partial \mathbi{u}/\partial t$) and
${\bvarepsilon}_p(x,t=0)=\mathbf{0}$, for all~$x$.

Thereafter, we consider the following hypotheses. Following several
works
\cite{ballard_contraintes_1991,heuze_laxwendroff_2017,lapostolle_fast_2023,lapostolle_modeling_2022} 
and with in mind the derivation of an exact solution, we assume a
uniaxial total strain field, such that
{\begin{equation}
{\bvarepsilon}=\varepsilon_{11}(x,t) \mathbi{e}_1 \otimes \mathbi{e}_1,
\end{equation}}\noindent
$\mathbi{e}_1$ being a unit vector in the direction of the propagation.
Experimentally, such a situation may be reached when the diameter of
the laser spot is much greater than the specimen thickness
\cite{ayad_modeling_2022,lapostolle_fast_2023}. In this paper, we
consider a LSP-type situation, where $L$ is very large (a few cm, as
opposed to the thickness of specimens used for backface velocity VISAR
measurements which are thin
\cite{seddik_identification_2022,ayad_modeling_2022}). We also make
the (standard) hypothesis of incompressible plasticity, so that we have
{\begin{equation}
{\bvarepsilon}_p=\ds \varepsilon_p(x,t)\left(\mathbi{e}_1 \otimes \mathbi{e}_1 - 
\tfrac{1}{2}\mathbi{e}_2 \otimes \mathbi{e}_2 - \tfrac{1}{2}\mathbi{e}_3 \otimes \mathbi{e}_3 \right).
\label{eq:plastic_strain}
\end{equation}}\noindent

For simplicity, the developments will be made in the case of
homogeneous isotropic material (even though they could be extended to
more complex material symmetries). Thus, the stiffness tensor is
expressed by:
{\begin{equation}
    \mathbb{C}=2\mu \mathbb{K}+3\kappa \mathbb{J},
    \label{eq:C_isotropic}
\end{equation}}\noindent
with 
\vspace*{2pt}
{\begin{equation}
\mathbb{J}=\ds \tfrac{1}{3}\mathbi{I} \otimes \mathbi{I} \quad 
\tn{and} \quad \mathbb{K}=\mathbb{I}-\mathbb{J},
\end{equation}}\noindent
$\mathbb{I}$ being the fourth-order identity tensor. In Equation
\eqref{eq:C_isotropic}, $\mu$ is the shear modulus and $\kappa$ the
bulk modulus. In these conditions, the stress tensor takes the
following form:
\vspace*{2pt}
{\begin{equation}
{\bsigma}=\sigma_{11}(x,t)\mathbi{e}_1 \otimes \mathbi{e}_1+\sigma_{22}
\left( \mathbi{e}_2 \otimes \mathbi{e}_2+\mathbi{e}_3 \otimes \mathbi{e}_3\right).
\end{equation}}\noindent
Finally, we consider the case of an elastic perfectly plastic behavior,
i.e. with no kinematic nor isotropic hardening: $\mathbi{X}=\mathbf{0}$ and
$R(p)=\sigma_{y}$, $\sigma_{y}$ being the yield strength of the
material. The plasticity criterion of Equation
\eqref{eq:dynamics_eq}$_4$ reduces to:
\vspace*{2pt}
{\begin{equation}
    f=\left\lvert \sigma_{11}-\sigma_{22} \right\rvert-\sigma_y.
    \label{eq:plast_crit}
\end{equation}}\noindent

It can be noted that in this case, the Tresca and von Mises yield
criteria adopt the same expression, given in 
Equation~\eqref{eq:plast_crit}. As pointed out by 
Heuz\'e~\cite{heuze_laxwendroff_2017}, it is useful to rewrite 
Equation~\eqref{eq:plast_crit} in the following form:
\vspace*{2pt}
{\begin{equation}
f=\ds \left\lvert \frac{6\mu}{3\kappa+4\mu}\sigma_{11}-\frac{9\kappa \mu}{3\kappa+4\mu}
\varepsilon_p \right\rvert-\sigma_y.
    \label{eq:f_alternate}
\end{equation}}\noindent

Notably, Equation \eqref{eq:f_alternate} allows the determination of
the axial stress value corresponding to the onset of plasticity in a
material with no residual plastic strain: 
\vspace*{2pt}
{\begin{equation}
\sigma_{\mathrm{hel}}=\left(\frac{2}{3}+\frac{\kappa}{2\mu} \right)\sigma_y,
\end{equation}}\noindent
where $\sigma_{\mathrm{hel}}$ is called the Hugoniot Elastic Limit (HEL).

\vspace*{2pt}

\subsubsection{Stress wave equation}

\vspace*{2pt}

We derive here the stress wave
propagation equation, which we put in the form of a hyperbolic system,
following the same developments as in previous works
\cite{heuze_laxwendroff_2017,lapostolle_modeling_2022,lapostolle_fast_2023}. 
Equation \eqref{eq:dynamics_eq}$_1$ combined with Equation
\eqref{eq:dynamics_eq}$_2$ derived with respect to time yields the
following system:
\vspace*{2pt}
{\begin{equation}
\left\{
\begin{array}{ll}
\ds \frac{\partial \sigma_{11}}{\partial x}=\rho \frac{\partial^2 u_1}{\partial t^2}=
\rho \frac{\partial v_1}{\partial t} \vspace*{0.3cm}\\
\ds \frac{\partial \sigma_{11}}{\partial t}=\left(\kappa+\frac{4}{3}\mu \right)
\frac{\partial \varepsilon_{11}}{\partial t}-2\mu \frac{\partial \varepsilon_p}{\partial t},
    \end{array}
    \right.
    \label{eq:system_init}
\end{equation}}\noindent
where $v_1$ is the axial velocity component. Under a small strain
hypothesis, we can write
\vspace*{2pt}
{\begin{equation}
\ds \frac{\partial \varepsilon_{11}}{\partial t}=\frac{\partial}{\partial t}
\left(\frac{\partial u_1}{\partial x} \right)=\frac{\partial}{\partial x}
\left(\frac{\partial u_1}{\partial t} \right)=\frac{\partial v_1}{\partial x}.
\end{equation}}\noindent

Thus the system of Equation \eqref{eq:system_init} can be rewritten as
\vspace*{2pt}
{\begin{equation}
    \left\{
    \begin{array}{ll}
\ds \frac{\partial \sigma_{11}}{\partial x}=\rho \frac{\partial v_1}{\partial t} \vspace*{0.3cm}\\
\ds \frac{\partial \sigma_{11}}{\partial t}=\left(\kappa+\frac{4}{3}\mu \right)
\frac{\partial v_{1}}{\partial x}-2\mu \frac{\partial \varepsilon_p}{\partial t}.
    \end{array}
    \right.
    \label{eq:system_inter}
\end{equation}}\noindent

Firstly we consider the case of a purely elastic propagation, i.e.
$\partial \varepsilon_p / \partial t=0$, Equation
\eqref{eq:system_inter} thus becomes 
\vspace*{2pt}
{\begin{equation}
\ds \frac{\partial \mathbi{U}}{\partial t}+\mathbi{A}_{\mathrm{el}}\frac{\partial \mathbi{U}}{\partial x}=\mathbf{0},
\label{eq:adv_el}
\end{equation}}\par\pagebreak\noindent
where
\vspace*{2pt}
{\begin{equation}
    \mathbi{U}=
    \begin{pmatrix}
    \sigma_{11} \\
    v_1
    \end{pmatrix} \quad \tn{and} \quad \mathbi{A}_{\mathrm{el}}=
    \begin{pmatrix}
    0 & \ds -\left( \kappa+\frac{4}{3}\mu \right) \\
    \ds -\frac{1}{\rho} & 0 
    \end{pmatrix}.
\end{equation}}\noindent

Equation \eqref{eq:adv_el} is the so-called advection equation, and can
be solved analytically and using dedicated numerical schemes to compute
the elastic stress wave propagation. 

We consider then the case where the propagation is elasto-plastic, i.e.
$\partial \varepsilon_p / \partial t \neq 0$. Here, by taking the time
derivative of the plasticity criterion (Equation
\eqref{eq:dynamics_eq}$_4$), and expressing the stress components using
the elastic behavior (Equation \eqref{eq:dynamics_eq}$_2$), we obtain
the following relation:
{\begin{equation}
\ds \frac{\partial \varepsilon_p}{\partial t}=\frac{2}{3}\frac{\partial \varepsilon_{11}}{\partial t}.
\label{eq:epsp_eps}
\end{equation}}\noindent
By injecting Equation \eqref{eq:epsp_eps} into Equation \eqref{eq:system_inter}, we obtain:
{\begin{equation}
\ds \frac{\partial \mathbi{U}}{\partial t}+\mathbi{A}_{\mathrm{pl}}\frac{\partial 
\mathbi{U}}{\partial x}=\mathbf{0},
\label{eq:adv_pl}
\end{equation}}\noindent
where the matrix $\mathbi{A}_{\mathrm{pl}}$ reads 
{\begin{equation}
    \mathbi{A}_{\mathrm{pl}}=
    \begin{pmatrix}
    0 & \ds -\kappa \\
    \ds -\frac{1}{\rho} & 0 
    \end{pmatrix}.
\end{equation}}\noindent

It can be noted that the eigenvalues (in absolute value) of matrices
$\mathbi{A}_{\mathrm{el}}$ and $\mathbi{A}_{\mathrm{pl}}$ provide the stress wave
velocities when the wave is elastic, noted $c_{\mathrm{el}}$ and when the
wave is plastic, noted  $c_{\mathrm{pl}}$:
{\begin{equation}
c_{\mathrm{el}}=\ds \sqrt{\ds \frac{\kappa+\frac{4}{3}\mu}{\rho}} \quad \tn{and} 
\quad c_{\mathrm{pl}}=\ds \sqrt{\frac{\kappa}{\rho}}.
\end{equation}}\noindent

We thus have $c_{\mathrm{el}}>c_{\mathrm{pl}}$, which means that a stress wave
loses speed when it induces plastic strains. Finally, the displacement
field can be recovered with:
{\begin{equation}
    u_1(x,t)=\int^t_{t=0}v(x,t){\rm d}t+u_1(x,t=0).
    \label{eq:u_1}
\end{equation}}\noindent

This quantity is interesting because the displacement of the front face
after the shock (i.e. at large $t$) can be experimentally measured by
profilometry
\cite{hu_investigation_2009,kamkarrad_fem_2016,achintha_eigenstrain_2011,xiang_numerical_2021}.

{In summary, Equation \eqref{eq:adv_el} is solved as an elastic
predictor step at each time increment. If the stress state complies
with the yield condition, then it is adopted as the final state, and
the simulation goes on to the next time increment. If not, then
Equation \eqref{eq:adv_pl} is solved, and the solution replaces the one
computed with the elastic predictor step.}  {The fields of
$\sigma_{22}$ and $\varepsilon_p$ can be recovered by using Equation
\eqref{eq:plast_crit}, and Equations \eqref{eq:dynamics_eq}$_2$ and
\eqref{eq:epsp_eps} respectively.}

\subsection{Exact solution}
Equations \eqref{eq:adv_el} and \eqref{eq:adv_pl} are hyperbolic
systems, and they can be solved analytically, using the methods of
characteristics, whose description can be found in several works
\cite{leveque_finite_2002,wang_foundations_2007,ballard_contraintes_1991,heuze_laxwendroff_2017}. 
The solution $\mathbi{U}$ of the advection equation is constant along lines
in the $(x,t)$ plane called characteristic lines, defined by the
following equations:
{\begin{equation}
\begin{cases}
 {\rm d}x=\pm c \, {\rm d}t & \tn{characteristic line equation} \\
 {\llbracket \sigma \rrbracket =} \aj{\mp}  {\rho c \llbracket v \rrbracket} & 
 {\tn{momentum conservation condition},}
\end{cases}
\label{eq:charac_lines}
\end{equation}}\noindent
where $c$ is a generic notation for the stress wave velocity, whether
it is elastic or plastic. From Equation \eqref{eq:charac_lines}$_1$, we
deduce the general equation for the characteristic lines: $t=x/c+t_0$.
These lines delimitate zones in the $(x,t)$ plane in which $\mathbi{U}$,
and thus $\sigma_{11}$ and $v_1$, are constant. There are as many
characteristic lines as there are wave fronts, whether they are loading
or unloading the material, or propagating following $\mathbi{e}_1$ or along
the reversed direction. Equation \eqref{eq:charac_lines}$_1$, by
convention, is written with a ``$+$'' sign when the wave is propagating
following $\mathbi{e}_1$, and with a ``$-$'' sign in the opposite direction.
Equation \eqref{eq:charac_lines}$_2$ gives a relation between the
stress and velocity variations from one side of the  {wave front to the
other \cite{wang_foundations_2007}}.  {The notation $\llbracket \cdot
\rrbracket$ corresponds to the difference of a given quantity between
each side of the wave front.} Thus, knowing the stress values from
either sides of a  {wave front}, and one value of velocity, one can
deduce the values of velocities at all instants of the propagation.
Here, following  {\cite{wang_foundations_2007,wang_unloading_2004}},
Equation \eqref{eq:charac_lines}$_2$ is used with a ``$-$'' sign when the
wave is propagating following $\mathbi{e}_1$,  {but the opposite sign must
be applied} when it is propagating in the opposite direction. One can
thus use these relations to plot a so-called $x$--$t$ diagram, displaying
the various characteristic lines, which can be used to compute the
stress and velocity states for any given time or location.

\subsection{Numerical implementation of the stress wave equation}\label{sec:num_methods}

We present here the numerical implementation of the stress wave
equation. As a hyperbolic system, the stress wave equation can be
solved using an explicit finite differences approach based on several
numerical schemes. We follow the methodology of
Leveque~\cite{leveque_finite_2002}, and introduce the space and time
increments, $\Delta x$ and $\Delta t$ respectively, such that
$x_{i+1}=x_i+\Delta x$ and $t^{n+1}=t^n+\Delta t$. The notation
$(\cdot)_i$ refers to quantities evaluated at $x=x_i$, and $(\cdot)^n$
at $t=t^n$. The explicit procedures that can be used for this problem
must however abide the so-called Courant--Friedrichs--Lewy (CFL)
condition:
{\begin{equation}
    \tn{CFL}=c_{\mathrm{el}}\frac{\Delta t}{\Delta x}\leq 1,
    \label{eq:cfl}
\end{equation}}\noindent
which implies that the lower $\Delta x$ is, the lower $\Delta t$ must
be for the CFL condition to be verified, which increases the number of
steps of the simulation. The elastic velocity is used in 
Equation~\eqref{eq:cfl}
so that the CFL condition is valid even when the wave is
propagating at speed $c_{\mathrm{pl}}$. The numerical solution of Equations
\eqref{eq:adv_el} and \eqref{eq:adv_pl} has the following form:
{\begin{equation}
    \mathbi{U}_i^{n+1}=\mathbi{U}_i^n-\frac{\Delta t}{\Delta x}\left(\mathbi{F}^n_{i+1/2}-\mathbi{F}^n_{i-1/2} \right).
    \label{eq:adv_flux}
\end{equation}}\noindent

The quantity $\mathbi{F}$ is called a flux, and its expression depends on
the chosen numerical scheme. In this work we will compare several
schemes, namely the Lax--Friedriechs (LF), Lax--Wendroff (LW), Godunov
and Godunov with High Resolution (GHR) schemes. Because of their
different expressions, these schemes introduce various degrees of
artificial viscosity or oscillations. As an example, the GHR scheme
actively tries to limit these phenomenon by adding or reducing the
artificial viscosity when the rates of the unknowns are high or low. A
comprehensive description of these schemes can be found in
Leveque~\cite{leveque_finite_2002} and a short summary is provided in
Appendix~\ref{app:num_schemes}. Additionally, this problem can be
solved using a finite element approach (and in particular the
commercial software ABAQUS). In this case, an axisymetric model is used
with an standard/explicit procedure, and with the pressure applied over
a very large spot (radius 1 cm) so as to mimic a uniaxial case. In 3D,
stress wave edge effects originate from the edges of the laser spot and
propagate inward, eventually introducing a perturbation at the center
of the impact \cite{ayad_modeling_2022,ecault_numerical_2016}. However
the larger the spot size, the longer this perturbation takes to arrive
at the center, which is where the analytical expressions are valid.

For all the numerical schemes, the expression of $\mathbi{F}$ depends on
matrices $\mathbi{A}_{\mathrm{el}}$ or $\mathbi{A}_{\mathrm{pl}}$ of 
Equations~\eqref{eq:adv_el} and \eqref{eq:adv_pl} depending on the elastic or
plastic nature of the propagation. The numerical resolution will thus
adopt the following strategy: a predictor step will be computed by
solving Equation \eqref{eq:adv_el}. \aj{In the case of the GHR scheme,
the flux limiter (see Equation \eqref{eq:superbee}) is also applied to
plastic waves during the correction.} Then the plasticity criterion
will be checked, and Equation~\eqref{eq:adv_pl} will be solved at the
points were the predictor steps violates the plasticity criterion. At
those points, the increment of plastic strain can be computed by using
Equation \eqref{eq:epsp_eps} in Equation \eqref{eq:system_init}$_2$:
{\begin{equation}
\Delta \varepsilon_p=\frac{2 \Delta \sigma_{11}}{3 \kappa},
\label{eq:delt_epsp_sig}
\end{equation}}\noindent
with $\Delta \sigma_{11}$ known from the resolution of the advection
equation. \aj{We specify that the notation $\Delta \left( \cdot
\right)$ applied to $\sigma_{11}$ or $\varepsilon_p$ refers to a time
increment: $\Delta \sigma_{11}=\sigma_{11}^{n+1}-\sigma_{11}^n$.} This
approach has the \unskip\break advantage of preserving numerical accuracy during the
plastic correction, since the same numerical method is used. Another
possible approach would be to compute the plastic correction not with
the stress wave equation but with a radial return procedure (see for
instance \cite{lapostolle_fast_2023}). Both approaches are almost
equivalent as shown in Appendix A of Lapostolle~\cite{lapostolle_etude_2022}. 
\aj{Using a radial return procedure was found to increase the numerical
viscosity while avoiding plastic strain overshoot at the boundary}.

\section{Residual stresses modeling}\label{sec:residual_stresses}

In this section we detail the hypotheses for the residual stresses
modeling, following several works of the literature
\cite{lapostolle_fast_2023,cellard_laser_2012,ahdad_contraintes_1996,johns_thermal_1965,gelineau_etude_2018,korsunsky_modelling_2005}. 
The modeling assumes a finite plate of thickness $L$ with infinite
dimensions in the plane directions, with a heterogeneous plastic strain
field along the thickness. The plate is allowed to bend. This plastic
strain field corresponds to the one computed by the propagation of the
stress wave in Section~\ref{sec:mech_model}, which is axisymmetric and
only varies in the $\mathbi{e}_1$ direction. The total strain field is also
axisymmetric and has the following expression:
{\begin{equation}
{\bvarepsilon}(x)=\varepsilon_{11}(x)\mathbi{e}_1+\varepsilon_{22}(x)
\left( \mathbi{e}_2 \otimes \mathbi{e}_2 + \mathbi{e}_3 \otimes \mathbi{e}_3 \right).
\label{eq:total_strain}
\end{equation}}\noindent

The resolution of the compatibility equations imposes the form of the
radial strain components:
{\begin{equation}
    \varepsilon_{22}(x)=\varepsilon_{33}(x)=\alpha x + \beta,
    \label{eq:eps_alph_beta}
\end{equation}}\noindent
where $\alpha$ and $\beta$ are constants to be computed so that the
residual stresses are at equilibrium along the thickness in the absence
of external forces. Because the gradient of the mechanical fields is
only non zero in the $\mathbi{e}_1$ direction, and since the axial stress
is 0 at the surface, the static equilibrium yields
{\begin{equation}
    \sigma_{11}(x)=0 \ \forall x.
    \label{eq:s11_null}
\end{equation}}\noindent

Using Hooke's law Equation \eqref{eq:dynamics_eq}$_2$ with Equations
\eqref{eq:total_strain}, \eqref{eq:plastic_strain} and
\eqref{eq:s11_null}, we have 
{\begin{equation}
\ds \sigma_{22}(x)=\sigma_{33}(x)=\frac{18\kappa\mu}{3\kappa+4\mu }
\left( \alpha x + \beta +\frac{\varepsilon_p(x)}{2} \right).
\label{eq:RS_expression}
\end{equation}}\noindent

The conditions to compute the constants $\alpha$ and $\beta$ of
Equation \eqref{eq:eps_alph_beta} are the enforcement of zero normal
force and bending moment:
{\begin{equation}
\begin{cases}
        \ds \int^L_0 \sigma_{22}(x){\rm d}x=0 & \tn{normal force}  \vspace*{0.35cm}\\
        \ds \int^L_0 \sigma_{22}(x)x{\rm d}x=0 & \tn{bending moment} ,
    \end{cases}
\end{equation}}\par\pagebreak\noindent
leading to 
{\begin{equation}
    \left\{
    \begin{array}{ll}
        \ds L\left(\frac{\alpha L}{2}+\beta \right)
        =-\int^L_0 \frac{\varepsilon_p(x)}{2}\,{\rm d}x=\Gamma_1  \vspace*{0.35cm}\\
        \ds L^2\left(\frac{\alpha L}{3}+\frac{\beta}{2} \right) 
        = -\int^L_0\frac{\varepsilon_p(x)}{2}x\,{\rm d}x=\Gamma_2,
    \end{array}
    \right.
    \label{eq:res_sys_integ}
\end{equation}}\noindent
which resolution leads to
{\begin{equation}
    \left\{
    \begin{array}{ll}
        \ds \alpha=\frac{12}{L^3}\left(\Gamma_2-\frac{L\Gamma_1}{2} \right) \vspace*{0.35cm}\\
        \ds \beta=\frac{2}{L}\left( 2\Gamma_1-\frac{3\Gamma_2}{L} \right).
    \end{array}
    \right.
\end{equation}}\noindent

Knowing the distribution of $\varepsilon_p(x)$ from the propagation
simulation, one is thus able to compute the residual stresses using
Equation \eqref{eq:RS_expression}. It is also possible to compute the
curvature radius of the plate \cite{korsunsky_modelling_2005}:
{\begin{equation}
R=\frac{L^3}{6\left(2\Gamma_2-L\Gamma_1 \right)}.
\end{equation}}\noindent

As we will see with the results of Section~\ref{sec:results}, the
plastic strain field has step-wise constant in-depth profile, which
will allow a fully analytical expression of $\Gamma_1$ and $\Gamma_2$
of Equation \eqref{eq:res_sys_integ}, and thus of the residual
stresses.

\section{Results}\label{sec:results}

In this section, we present results for the stress wave propagation,
both analytically and numerically. We distinguish three cases: $P\leq
\sigma_{\mathrm{hel}}$, $\sigma_{\mathrm{hel}}<P \leq 2\sigma_{\mathrm{hel}}$ and 
{$P>(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}>2\sigma_{\mathrm{hel}}$}, $P$ being the
positive plasma pressure.  {The intermediary subcase of
$P>2\sigma_{\mathrm{hel}}$, namely $2\sigma_{\mathrm{hel}}<
P\leq(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$ is detailed below}. Indeed when
a material is loaded such that $\sigma_{11}=-P<-\sigma_{\mathrm{hel}}$, the
induced plastic strain is
{\begin{equation}
\varepsilon_p=\ds -\frac{2}{3}\frac{P-\sigma_{\mathrm{hel}}}{\kappa},
\end{equation}}\noindent
according to the plasticity criterion of Equation
\eqref{eq:f_alternate}. The value of axial stress inducing plastic
strain during the unloading, given the already induced plastic strain,
is $P-2\sigma_{\mathrm{hel}}$. Given that the value of the axial stress cannot
change sign during the unloading, then no plastic strains are induced
during the unloading if $P \leq 2\sigma_{\mathrm{hel}}$. On the contrary,
plastic strain will be induced during the unloading if
$P>2\sigma_{\mathrm{hel}}$.  {A subcase of this last condition is
$2\sigma_{\mathrm{hel}}< P\leq(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$, whose exact
solution structure is presented in Wang~\cite{wang_unloading_2004}. The
main difference is that though plastic strains are induced during the
unloading in both subcases (same exact solution in the initial
instants), the case with $P>(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$ induces
plastic strains over a greater depth, changing the subsequent layout of
the exact solution compared to the other subcase, since it gives rise
to different wave interactions. See Figure~\ref{fig:diag_lag_5GPa} for
more details.}

In the following, we consider a steel-like material with $\kappa=175$
GPa, $\mu=80.8$ GPa, $\rho=7800$~kg${\cdot}$m$^{-3}$, and with a yield
strength of $\sigma_{y}=870$ MPa, leading to
$\sigma_{\mathrm{hel}}=1522.5$ MPa. The loading remains a step function
of amplitude $-P$ and of duration $t_0=80$ ns. The boundary conditions
are specified in Equation \eqref{eq:dynamics_eq}$_6$. For the whole
section, numerical results will be presented, using the following
parameters: the length of the domain is $L=14$ mm, the duration of the
simulation is $T=800$ ns, the spatial increment is $\Delta
x=3.5~\upmu$m, and the CFL number if $0.95$.\pagebreak

Hereafter, in order to shorten the notations, the axial stress
component $\sigma_{11}$ will be simply noted $\sigma$, and the axial
velocity component $v_1$ simply $v$.

\subsection{Case 1: $P\leq  \sigma_{\mathrm{hel}}$}

This trivial case of pure elastic response allows to illustrate the
$x$--$t$ diagram for a simple solution. The corresponding $x$--$t$ diagram is
presented in Figure~\ref{fig:diag_lag_elas}, in which two parallel
characteristics lines are plotted, one for loading front of the wave,
the other for the unloading front, the latter being shifted by $t_0$.


\begin{figure}[h!]
\includegraphics{fig02}
\caption{$x$--$t$ diagram in the case of a purely elastic propagation. The
letters identify areas in the diagram. \label{fig:diag_lag_elas}}
\end{figure}

For each zones of the $x$--$t$ diagram of Figure~\ref{fig:diag_lag_elas},
it is straightforward to attribute an axial stress value. Using the
compatibility equation of Equation \eqref{eq:charac_lines}$_2$, one can
then compute the velocity values of the $x$--$t$ diagram. Their
expressions are given in Table~\ref{tab:diag_elas}.

%tab1
\begin{table}[h!]
\caption{Values of axial stress and velocities of
the $x$--$t$ diagram of Figure~\ref{fig:diag_lag_elas}\label{tab:diag_elas}}
\begin{tabular}{ccc}
\thead
Zone $i$ of the $x$--$t$ diagram  & $\sigma_i$ & $v_i$ \\
\endthead
    a & 0 & 0 \\
    b & $-P$ & $\ds \frac{P}{\rho c_e}$ \\
    c & 0 & 0 
\botline
    \end{tabular}
\end{table}


As no plastic deformation occurs in this case, no residual stress is
generated by the shock wave. This methodology will be applied to the
more complex wave propagations, where plastic strain is induced.

\subsection{Case 2: $\sigma_{\mathrm{hel}}< P\leq 2\sigma_{\mathrm{hel}}$}

We now consider the case where $\sigma_{\mathrm{hel}}<P\leq 2\sigma_{\mathrm{hel}}$, with
$P=2$ GPa. Figure~\ref{fig:diag_lag_2GPa} presents the corresponding
exact $x$--$t$ diagram, as well a numerical $x$--$t$ diagram computed using
the Godunov High Resolution method presented in
Section~\ref{sec:num_methods}. 

\begin{figure}
\includegraphics{fig03}
\caption{$x$--$t$ diagrams in the case where $\sigma_{\mathrm{hel}}<P\leq
2\sigma_{\mathrm{hel}}$. (a) Exact $x$--$t$ diagram. The letters identify areas
within the diagram, the number the wave front which induces plasticity.
(b)~Comparison between the exact and numerical $x$--$t$ diagrams for the
axial stress. \label{fig:diag_lag_2GPa}}
\end{figure}

\pagebreak
Figure~\ref{fig:diag_lag_2GPa}(a) presents characteristic lines with
different slopes. The red line, identified by the number (1)
corresponds to the plastic loading front, which propagates at the speed
$c_p$. Since plastic strain does not occur during the unloading in this
case, no characteristic line of slope $1/c_p$ originates at $t_0$. The
values of stress and velocity for the $x$--$t$ diagram of
Figure~\ref{fig:diag_lag_2GPa}(a) are summarized in
Table~\ref{tab:diag_lag_2GPa}.

%tab2
\begin{table}[h!]
\caption{Values of axial stress and velocity in the zones of the $x$--$t$
diagram of Figure~\ref{fig:diag_lag_2GPa}\label{tab:diag_lag_2GPa}}
\begin{tabular}{ccc}
\thead
Zone $i$ of the $x$--$t$ diagram & $\sigma_{i}$ & $v_{i}$ \\
\endthead
a & 0 & 0 \\
b & $-\sigma_{\mathrm{hel}}$ & $\ds \frac{\sigma_{\mathrm{hel}}}{\rho c_e}$ \vspace*{4pt}\\
c & $-P$ & $\ds \frac{Pc_e+\sigma_{\mathrm{hel}}(c_p-c_e)}{\rho c_e c_p}$\vspace*{4pt}\\
d & 0 & $\ds \frac{\left(P-\sigma_{\mathrm{hel}}\right)(c_e-c_p) }{\rho c_e c_p}$\vspace*{4pt}\\
e & $\ds \frac{\left(P-\sigma_{\mathrm{hel}} \right)\left(c_p-c_e\right)}{2c_p}$ & $\ds \frac{\left(P-\sigma_{\mathrm{hel}}\right)(c_e-c_p) }{2\rho c_e c_p}$\vspace*{4pt} \\
f & 0 & 0 
\botline
\end{tabular}
\end{table}

Table~\ref{tab:epsp_2GPa} presents the exact value of plastic strain
induced by the plastic loading front in the $x$--$t$ diagram of
Figure~\ref{fig:diag_lag_2GPa}. It can be computed using the plasticity
criterion of Equation \eqref{eq:dynamics_eq}$_4$.

%tab3
\begin{table}[h!]
\caption{Values of the axial plastic strain
induced by the plastic front of the $x$--$t$ diagram of
Figure~\ref{fig:diag_lag_2GPa}\label{tab:epsp_2GPa}}
\begin{tabular}{cc}
\thead
Plastic front $i$ of the $x$--$t$ diagram & $\varepsilon_{p,i}$ \vspace*{1pt}\\
\endthead
1 & $\ds -\frac{2}{3}\frac{P-\sigma_{\mathrm{hel}}}{\kappa}$ \raisebox{1.5pc}{}\vspace*{2pt}
\botline
\end{tabular}
\end{table}

\pagebreak
For each of the zones (identified with letters from (a) to (f)) of the
$x$--$t$ diagram of Figure~\ref{fig:diag_lag_2GPa}, the axial stress and
velocity can be computed. During the loading front, the stress in zone
(b) reaches the lowest value possible without inducing plasticity, i.e.
$-\sigma_{\mathrm{hel}}$. Once plasticity is induced, the stress reaches the
prescribed value $-P$ in zone (c), before coming back to 0 once the
wave has passed in zone (d). For each of these zones, the velocity can
be deduced using Equation \eqref{eq:charac_lines}$_2$. Since the
elastic and plastic fronts do not propagate at the same speed, the
elastic unloading front eventually catches up the plastic loading
front. Because of the different speeds, the fronts interact with each
other, in the similar way of an interface between two materials with
different impedance. Transmitted and reflected waves are created, thus
delimiting here zone (e). To compute the values of stress and velocity
in this zone, one can take advantage of the loading fronts that zone
(e) shares with zones for which the stresses and velocities are already
known. In this instance, one must solve the following linear system:
{\begin{equation}
    \left\{
    \begin{array}{ll}
        \sigma_{(e)}-\sigma_{(d)}=\rho c_{(e)} (v_{(e)}-v_{(d)} ) \vspace*{3pt}\\
        \sigma_{(e)}-\sigma_{(b)}=-\rho c_{(e)} (v_{(e)}-v_{(b)}).
    \end{array}
    \right.
\label{eq:sys_s_v}
\end{equation}}\noindent

Eventually, the plastic loading front is overtaken by the elastic
unloading front, and the stress wave stops inducing plastic strains. At
this point, the stress is attenuated (in absolute value) to the highest
value of stress that does not generate plastic strains, i.e.
$\sigma_{\mathrm{hel}}$. 

We can now use the different values of Tables~\ref{tab:diag_lag_2GPa}
and~\ref{tab:epsp_2GPa} with the $x$--$t$ diagram of
Figure~\ref{fig:diag_lag_2GPa} to compute analytically the stress,
velocity and plastic strain distributions at any given time $t$. We can
see in Figure~\ref{fig:diag_lag_2GPa}(b) that the numerical solution
corresponds very well to the exact $x$--$t$ diagram, as zones of different
stress values are clearly identifiable and delimited by the
characteristic lines. These numerical results obtained using various
numerical schemes can be compared to the exact solution, which is
presented in Figure~\ref{fig:s_v_ep_2GPa}.

\begin{figure}
\vspace*{-2pt}
\includegraphics{fig04}
\vspace*{-2pt}
\caption{Comparisons of the mechanical fields obtained with
the exact solution and with the various numerical methods. (a) Axial
stress field at $t=100$ ns. (b) Axial stress field at $t=490$ ns. (c)
Axial velocity field at $t=100$ ns. (d) Axial velocity field at $t=490$
ns. (e) Axial plastic strain field at $t=100$ ns. (f) Axial plastic
strain field at $t=490$ ns.\label{fig:s_v_ep_2GPa}}
\vspace*{-2pt}
\end{figure}

The results of Figure~\ref{fig:s_v_ep_2GPa} show that all numerical
methods used are in a globally good agreement with the exact solution,
either for the stress, velocity or plastic strain.  {It must be noted
that, for the sake of visibility, the points in the plots of
Figure~\ref{fig:s_v_ep_2GPa} do not correspond to each node, but rather
a fixed interval of nodes. This also applies for future figures.} Each
scheme however has different numerical properties, which can be
observed. For instance, the Lax--Wendroff scheme, though precise,
introduce some small amounts of artificial oscillations, which can be
observed after the unloading front in the stress and velocity profiles.
The Lax--Friedrichs scheme is characterized by a large amount of
artificial viscosity, which can clearly be seen in the plastic strain
profile of Figure~\ref{fig:s_v_ep_2GPa}(f). In all instances, the
Godunov with High Resolution scheme proves to be the most precise, as
it follows the more closely the exact solution. A quantification of the
errors of the schemes with respect to the exact solution will be
discussed in Section~\ref{sec:discussion}. For the plastic strain, the
first points ($x\simeq 0$ mm) are far from the exact solution. This is
caused by  {a large stress increment} at the first cell ($x=0$)  {due
to the instantaneous nature of the loading}. As a result, the computed
plastic strains exhibits this locally high values.  {This behavior was
also observed for the finite elements solution in the work of
Heuz\'e~\cite{heuze_laxwendroff_2017}, and is discussed in more details 
in\break
Section~\ref{sec:discussion}.} 

The plastic strain profile of Figure~\ref{fig:s_v_ep_2GPa}(f)
corresponds to the residual plastic strain field (after propagation)
and can thus be used to compute the residual stresses presented in
Section~\ref{sec:residual_stresses}. Moreover, the value of the induced
plastic strain, and the depth along which it is induced are known
analytically, which allows the exact expression of the integrals in
Equation \eqref{eq:res_sys_integ}. We thus have:
{\begin{equation}
\left\{
\begin{array}{ll}
\ds \alpha=\frac{2c_e c_p t_0 \left(P-\sigma_{\mathrm{hel}} \right)\left(L(c_p-c_e)+c_e c_p t_0 \right)}
{L^3\kappa (c_e-c_p)^2} \quad \\
\ds \beta=\frac{c_e c_p t_0\left(P-\sigma_{\mathrm{hel}} \right)
\left(4L(c_e-c_p)-3c_e c_p t_0 \right)}
{3L^2\kappa(c_e-c_p)^2}.
\end{array}
\right.
\end{equation}}\noindent

The comparison between the exact and numerical residual stresses is
presented in Figure~\ref{fig:RS_2GPa}. For the numerical residual
stresses, the integrals of Equation \eqref{eq:res_sys_integ} are
computed numerically using a trapezoidal rule based of the numerical
plastic strain profiles.

\begin{figure}
\includegraphics{fig05}
\caption{Comparison of the residual stress field obtained
with the exact solution and the various numerical 
methods.\label{fig:RS_2GPa}}
\end{figure}

The analytical (exact) and numerical residual stress profiles are in
very good agreement which is expected since the plastic strain profiles
were already close between the different methods. The residual stresses
are in compression at the surface, and become in tension further in
depth (at $x \approx 1.9$ mm) due to the equilibrium. The surface
values of the numerical residual stress profiles are much more in
compression than the rest, which is due to the numerical errors at the
beginning of the plastic strain profiles.

In a similar manner, we can take advantage of the piece-wise constant
nature of the velocity profiles at any given time $t$ to give an exact
expression of the displacement field. More precisely, we will take an
interest in the surface displacement field, since it can be
experimentally measured. To illustrate this, we compute the time
evolution of the displacement at the surface of the specimen ($x=0$)
using Equation \eqref{eq:u_1} with the exact solution and the numerical
methods. This is presented in Figure~\ref{fig:ux_2GPa}. The agreement
between the exact solution and the numerical methods is relatively
good, with some errors introduced by the differences in the velocity
predictions, which lead to slightly different slopes in the
displacement profiles. The GHR scheme is the closest to the exact
solution.

\begin{figure}
\includegraphics{fig06}
\caption{Comparison of the temporal profile of the displacement at the
surface of the specimen ($x=0$) obtained with the exact solution and
the various numerical methods. \label{fig:ux_2GPa}}
\end{figure}

The final value of the displacement at the surface (i.e. once the
stress wave has completely passed) can be computed, and reads:
{\begin{equation}
u_{1,\mathrm{final}}=\frac{t_0\left(P(c_e+2c_p)-\sigma_{\mathrm{hel}}(c_e+c_p) \right)}{\rho c_e c_p}.
\label{eq:ux_2GPa}
\end{equation}}\noindent

\subsection{Case 3: $P>(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$}
We finally consider here the more complex case where the loading
amplitude $P$ is such that  {$P>(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$}. We
use the value of $P=5$ GPa to illustrate this case. The main
particularity in this case is that plastic strains will be induced
during the unloading front. This will also create a plastic
characteristic line during the unloading, creating more complex
interactions between the elastic and plastic parts of the wave. This is
illustrated in\break Figure~\ref{fig:diag_lag_5GPa}.

\begin{figure}[h!]
\vspace*{6pt}
\includegraphics{fig07}
\vspace*{6pt}
\caption{$x$--$t$ diagrams in the case where  {$P>
(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$}. (a) Exact $x$--$t$ diagram. The
letters identify areas within the diagram, the number the wave front
which induces plasticity. (b) Enlarged view of the $x$--$t$ diagram. (c)
Comparison between the exact and numerical $x$--$t$ diagrams for the axial
stress. \label{fig:diag_lag_5GPa}}
\vspace*{6pt}
\end{figure}

Figure~\ref{fig:diag_lag_5GPa}(a) presents a much more complex $x$--$t$
diagram than the one in Figure~\ref{fig:diag_lag_2GPa}. Notably, a
plastic wave front is created during the unloading front. The elastic
unloading wave front is thus trapped between two plastic wave fronts of
lower velocity, creating a complex pattern of transmitted and reflected
wave. The stress and velocity values associated with each zone of
Figure~\ref{fig:diag_lag_5GPa}(a) are summarized in
Table~\ref{tab:diag_lag_5GPa}, and the values of plastic strains
induced by the plastic wave fronts in Table~\ref{tab:epsp_5GPa}. These
values of stress and velocity are determined either by mechanical
analysis (for instance, it is straightforward to associate the stress
value of zone b with ${-}\sigma_{\mathrm{hel}}$ and the one of zone c
with ${-}P$) when possible, or either by solving a linear system
similar to  Equation~\eqref{eq:sys_s_v} when a zone is delimited by two
wave fronts that can be\break exploited. 

%tab4
\begin{table}
\caption{Values of axial stress and velocity
in the zones of the $x$--$t$ diagram of Figure~\ref{fig:diag_lag_5GPa}
\label{tab:diag_lag_5GPa}\vspace*{4pt}}
\fontsize{8}{11}\selectfont\tabcolsep=2pt
\begin{tabular}{ccc}
\thead
$i$ & $\sigma_{i}$ & $v_{i}$ \\
\endthead
a & 0 & 0 \\
b & {$-\sigma_{\mathrm{hel}}$} & {$ \frac{\sigma_{\mathrm{hel}}}{\rho c_{e}}$}\vspace*{3pt} \\
c & {$-P$} & {$ \frac{P c_e+\sigma_{\mathrm{hel}}(c_p-c_e)}{\rho c_e c_p}$} \vspace*{3pt}\\
d & {$-P+2\sigma_{\mathrm{hel}}$} & {$ \frac{P c_e-\sigma_{\mathrm{hel}}(c_e+c_p)}{\rho c_e c_p}$} \vspace*{3pt}\\
e & 0 & {$ \frac{\sigma_{\mathrm{hel}} (c_{e} - c_{p})}{\rho c_{e} c_{p}}$} \vspace*{3pt}\\
f & {$ -\frac{P(c_e+c_p)-4\sigma_{\mathrm{hel}}c_p}{c_e+c_p}$} & {$ \frac{Pc_e(c_e+ c_p)+\sigma_{\mathrm{hel}}(c_p^{2}-c_e^{2}-4c_e c_p)}{\rho c_e c_p(c_e+c_p)}$} \vspace*{3pt}\\
g & {$ -\frac{4\sigma_{\mathrm{hel}}c_p(c_e-c_p)}{(c_e+c_p)^{2}}$} & {$ \frac{\sigma_{\mathrm{hel}} (c_{e} - c_{p})^{2} (c_{e} + 3 c_{p})}{\rho c_{e} c_{p} (c_{e} + c_{p})^{2}}$} \vspace*{3pt}\\
h & {$ \frac{(c_e-c_p)\left(P(c_e+c_p)-2\sigma_{\mathrm{hel}}(c_e+3c_p) \right)}{2c_p(c_e+c_p)}$} & {$ \frac{(c_{e} - c_{p}) \left(P (c_{e} +  c_{p}) - 4 \sigma_{\mathrm{hel}} c_{p}\right)}{2 \rho c_{e} c_{p} (c_{e} + c_{p})}$} \vspace*{3pt}\\
i & {0} & {$ \frac{\sigma_{\mathrm{hel}} (c_{e} - c_{p}) \left(c_{e}^2 + 2c_e c_p - 7 c_{p}^2\right)}{\rho c_{e} c_{p} (c_{e} + c_{p})^2}$} \vspace*{3pt}\\
j & {$ -\frac{(c_e-c_p)\left(P(c_e^3+c_e^2c_p)-c_ec_p^2-c_p^3\right)+\sigma_{\mathrm{hel}}(-c_e^3-3c_e^2c_p+13c_e^2c_p-c_p^3)}{2 c_p(c_e+c_p)^3}$} & {$ \frac{(c_e-c_p)\left(P(c_e^3+c_e^2c_p-c_ec_p^2-c_p^3)+\sigma_{\mathrm{hel}}(-c_e^3-3c_e^2c_p+13c_ec_p^2-c_p^3) \right)}{2\rho c_ec_p(c_e+c_p)^3}$}\vspace*{3pt}\\
k & {$ \frac{(c_e-c_p)\left(P(c_e+c_p)^2+2\sigma_{\mathrm{hel}}(-c_e^2-4c_ec_p+c_p^2)\right)}{2c_e(c_e+c_p)^2}$} & {$ \frac{(c_e-c_p)\left(P(c_e+c_p)^2-4\sigma_{\mathrm{hel}}(c_ec_p+3c_p^2) \right)}{2\rho c_ec_p(c_e+c_p)^2}$} \vspace*{3pt}\\
l & {$- P + 2 \sigma_{\mathrm{hel}}$} & {$ \frac{Pc_e(c_e+c_p)+\sigma_{\mathrm{hel}}(3c_p^2-c_e^2-6c_e c_p)}{\rho c_e c_p (c_e+c_p)}$}  \vspace*{3pt}\\
m & {$ -\frac{P(c_e+c_p)^{2}-8\sigma_{\mathrm{hel}}c_e c_p}{(c_e+c_p)^{2}}$} & {$ \frac{Pc_e(c_e+c_p)^{2}+\sigma_{\mathrm{hel}}(c_p^{3}+c_e c_p^{2}-9c_e^{2}c_p-c_e^{3})}{\rho c_e c_p (c_e+c_p)^{2}}$}\vspace*{3pt} \\
n & {0} & {$ \frac{(c_e-c_p)\left( P(c_e+c_p)-\sigma_{\mathrm{hel}}(c_e+5c_p)\right)}{\rho c_e c_p (c_e+c_p)}$} \\
o & {$ \frac{(c_e-c_p)\left(P(c_e+c_p)^{2}-2\sigma_{\mathrm{hel}}(c_e^{2}+6c_e c_p+c_p^{2}) \right)}{2c_p(c_e+c_p)^{2}}$} & {$ \frac{(c_e-c_p)\left(P(c_e+c_p)^{2}-8\sigma_{\mathrm{hel}}c_e c_p \right)}{2\rho c_e c_p (c_e+c_p)^{2}}$} \vspace*{3pt}\\
p & 0 & {$ \frac{(c_e-c_p)\left(P(c_e+c_p)^3+2\sigma_{\mathrm{hel}}c_p(-3c_e^2-10c_ec_p+c_p^2) \right)}{2\rho c_e^2 c_p (c_e+c_p)^2}$} \vspace*{3pt}\\
q & {$ \frac{(c_e-c_p)\left(P(c_e+c_p)^2-2\sigma_{\mathrm{hel}}(c_e^2+6c_ec_p+c_p^2) \right)}{(c_e+c_p)^3}$} & {$ \frac{(c_e-c_p)\left(Pc_e(c_e+c_p)^2+\sigma_{\mathrm{hel}}(-c_e^3-5c_e^2c_p+c_ec_p^2-3c_p^3) \right)}{\rho c_ec_p(c_e+c_p)^3}$} \vspace*{3pt}\\
r & {$ \frac{(c_e-c_p)\left(P(c_e+c_p)^2-4\sigma_{\mathrm{hel}}(c_e^2+3c_ec_p) \right)}{(c_e+c_p)^3}$} & {$ \frac{(c_e-c_p)\left(Pc_e(c_e+c_p)^2+\sigma_{\mathrm{hel}}(-c_e^3-7c_e^2c_p+c_ec_p^2-c_p^3) \right)}{\rho c_ec_p(c_e+c_p)^3}$} \vspace*{3pt}\\
s & {$ -\frac{2\sigma_{\mathrm{hel}}(c_e-c_p)^2}{(c_e+c_p)^2}$} & {$ \frac{(c_e-c_p)\left(P(c_e+c_p)^2-\sigma_{\mathrm{hel}}(c_e^2+8c_e c_p+3c_p^2) \right)}{\rho c_e c_p (c_e+c_p)^2}$} \vspace*{3pt}\\
t & 0 & {$ \frac{(c_e-c_p)\left( P(c_e+c_p)^2-\sigma_{\mathrm{hel}}(c_e^2+10c_e c_p+c_p^2) \right)}{\rho c_e c_p (c_e + c_p)^2}$} \vspace*{3pt}\\
u & {$ \frac{(c_e-c_p)\left(P(c_e+c_p)^3+2\sigma_{\mathrm{hel}}(-c_e^3-9c_e^2c_p+c_ec_p^2+c_p^3) \right)}{2c_e(c_e+c_p)^3}$} & {$ \frac{(c_e-c_p)\left(P(c_e+c_p)^3-4\sigma_{\mathrm{hel}}c_p(c_e^2+6c_ec_p+c_p^2) \right)}{2\rho c_e c_p(c_e+c_p)^3}$} \vspace*{3pt}\\
v & {$ -\frac{\sigma_{\mathrm{hel}}(c_e-c_p)}{2c_p}$} & {$ \frac{\sigma_{\mathrm{hel}}\left( c_e-c_p \right)}{2\rho c_e c_p}$}\vspace*{4pt}
\botline
\end{tabular}
\vspace*{8pt}
\end{table}

%tab5
\begin{table}
\caption{Values of the axial plastic strain
induced by the plastic front of the $x$--$t$ diagram of
Figure~\ref{fig:diag_lag_5GPa}\label{tab:epsp_5GPa}\vspace*{4pt}}
\begin{tabular}{cc}
\thead
\parbox[t]{3cm}{\centering Plastic front $i$ of the $x$--$t$ diagram} 
& $\varepsilon_{p,i}$\vspace*{2pt} \\
\endthead
1 & $\ds -\frac{2}{3}\frac{P-\sigma_{\mathrm{hel}}}{\kappa}$ \raisebox{1.5pc}{}\vspace*{5pt} \\
2 & $\ds -\frac{2}{3}\frac{\sigma_{\mathrm{hel}}}{\kappa}$  \vspace*{5pt}\\
3 & $\ds \frac{2\left(\sigma_{\mathrm{hel}}(c_e+5c_p)-P(c_e+c_p) \right)}{3\kappa (c_e+c_p)}$  \vspace*{3pt}\\
4 & $\ds -\frac{2\sigma_{\mathrm{hel}}(c_e^2+6c_e c_p-3c_p^2)}{3\kappa(c_e+c_p)^2}$  \vspace*{3pt}\\
5 & $\ds -\frac{2\left( P(c_e+c_p)^{2}-\sigma_{\mathrm{hel}}(c_e^{2}+c_p^{2}+10c_e c_p) \right)}{3\kappa(c_e+c_p)^{2}}$  \vspace*{3pt}\\
6 & $\ds \frac{P(c_e^3+c_e^2c_p-c_ec_p^2-c_p^3)+2\sigma_{\mathrm{hel}}(-2c_e^3-5c_e^2c_p+4c_ec_p^2-c_p^3)}{3\kappa c_e(c_e+c_p)^2}$  \vspace*{3pt}\\
7 & $\ds 2\frac{P(c_e^3+c_e^2c_p-c_ec_p^2-c_p^3)+\sigma_{\mathrm{hel}}(-3c_e^3-13c_e^2c_p+7c_ec_p^2+c_p^3)}{3\kappa(c_e+c_p)^3}$ \vspace*{3pt}\\
8 & $\ds \frac{P(c_e^3+c_e^2c_p-c_ec_p^2-c_p^3)+2\sigma_{\mathrm{hel}}(-2c_e^3-5c_e^2c_p+4c_ec_p^2-c_p^3)}{3\kappa c_e(c_e+c_p)^2}$ \vspace*{3pt}\\
9 & { $\ds \frac{P(c_e^4+2c_e^3c_p-2c_ec_p^3-c_p^4)+2\sigma_{\mathrm{hel}}(-2c_e^4-11c_e^3c_p+7c_e^2c_p^2-c_ec_p^3-c_p^4)}{3\kappa c_e(c_e+c_p)^3}$}\vspace*{2pt}
\botline
\end{tabular}
\vspace*{6pt}
\end{table}

\pagebreak
In particular, the stress in zone (f) is here high enough to induce
plastic strains. This is true only if
$\sigma_{(f)}>\sigma_{\mathrm{hel}}$, which yields
$P>(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$, hence the definition of
the loading case considered.\footnote{{An alternative way to derive
this condition is by making the hypothesis of the wave front (3) being
elastic, which yields a value of $\sigma_{(f)}$ different from the one
in Table~\ref{tab:diag_lag_5GPa}. However the condition
$\sigma_{(f)}<\sigma_{\mathrm{hel}}$ will still lead to the same
threshold of $(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$.}} The $x$--$t$
diagram for $2\sigma_{\mathrm{hel}}\leq
P<(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$ is presented in
Wang~\cite{wang_unloading_2004}. In effect, the results for zones (a) to (e)
are the same as the ones presented in Table~\ref{tab:diag_lag_5GPa} and
fronts (1) and (2) of Table~\ref{tab:epsp_5GPa}, but the subsequent
structure of the $x$--$t$ diagram, and thus of the exact solution, is
different.

When compared to the characteristic lines, the numerical results
(obtained with the GHR scheme) in Figure~\ref{fig:diag_lag_5GPa}(b)
correspond very well to the stress zones displayed. In
Figure~\ref{fig:s_v_ep_5GPa}, we compare results of stress, velocity
and plastic strain between the numerical and exact solutions for the
various numerical schemes at several time steps.

\begin{figure}
\vspace*{2pt}
\includegraphics{fig08}
\vspace*{2pt}
\caption{Comparisons of the mechanical fields obtained with the exact
solution and with the various numerical methods. (a) Axial stress field
at $t=440$ ns. (b) Axial stress field at $t=790$ ns. (c) Axial velocity
field at $t=440$ ns. (d) Axial velocity field at $t=790$ ns. (e) Axial
plastic strain field at $t=440$ ns. (f) Axial plastic strain field at
$t=790$ ns. \label{fig:s_v_ep_5GPa}}
\vspace*{4pt}
\end{figure}

All the numerical schemes achieve a relatively good approximation of
the exact solution, but with more pronounced differences than for the
$P=2$ GPa case. These differences correspond to different degrees of
artificial viscosities in the numerical schemes, leading some of them
to miss certain fine features of the exact solution. In fact, only the
GHR scheme is able to correctly reproduce the smallest variations of
stress, velocity and plastic strains (see for instance the last step of
plastic strain in Figure~\ref{fig:s_v_ep_5GPa}(f) at $x\approx3.5$ mm).
The other numerical schemes sometimes completely miss the fine features
of the exact solution (for instance the part of the stress wave in
tension at $x\approx0.5$ mm in Figure~\ref{fig:s_v_ep_5GPa}(b), which
only the LW and GHR schemes reproduce to some~degree). Here again, the
errors of the scheme with respect to the exact solution are discussed
in Section~\ref{sec:discussion}.

Similarly to Figure~\ref{fig:RS_2GPa}, it is possible to compute the
residual stresses for the current case. The plastic strain profiles
being stepwise constant, it is possible to analytically evaluate the
integrals from Equation \eqref{eq:res_sys_integ}. The resulting
coefficients $\alpha$ and $\beta$ are detailed in
Appendix~\ref{sec:app_alpha_beta} for the sake of conciseness.
Figure~\ref{fig:RS_5GPa} presents the final plastic strain profile once
the stress wave stops inducing them, and the corresponding residual
stress profile.

\begin{figure}
\includegraphics{fig09}
\caption{Comparisons of the mechanical fields obtained with the exact
solution and with the various numerical methods. (a) Stationary plastic
strain field at $t=1050$ ns. (b) Residual stress field.
\label{fig:RS_5GPa}}
\end{figure}

Here again, we can compute the temporal evolution of the displacement
of the specimen surface using Equation \eqref{eq:u_1}. The results are
presented in Figure~\ref{fig:ux_5GPa}. The agreement between the exact
and numerical solution is here again good, the GHR scheme proving the
best to approximate the exact solution. Contrary to the profile of
Figure~\ref{fig:ux_2GPa}, the one of Figure~\ref{fig:ux_5GPa} presents
a decreasing slope at $t\approx 740$ ns. This is due to the secondary
wave which has been reflected in tension at the specimen surface (see
Figure~\ref{fig:s_v_ep_5GPa}(b) at $x\approx 0.5$ mm).

\begin{figure}
\includegraphics{fig10}
\caption{Comparison of the temporal profile of the displacement at the
surface of the specimen ($x=0$) obtained with the exact solution and
the various numerical methods. \label{fig:ux_5GPa}}
\end{figure}

The exact value of the final displacement of the surface is given in
that case by:
{\begin{eqnarray}
u_{1,\tn{final}}&= & \frac{-t_0}{2\rho c_e^2 c_p^2(c_e+c_p)^3}
\left(P(c_e^6-9c_e^4c_p^2-18c_e^3c_p^3-15c_e^2c_p^4-6c_ec_p^5-c_p^6) \right.\nonumber \\
&&  \left.+\, 2\sigma_{\mathrm{hel}}(-c_e^6-4c_e^5c_p+c_e^4c_p^2+24c_e^3c_p^3+c_e^2c_p^4+12c_ec_p^5-c_p^6) \right) .
\label{eq:ux_5GPa}
\end{eqnarray}}\noindent

\section{Discussion}\label{sec:discussion}

For both $P=2$ GPa and $P=5$ GPa, since $c_e>c_p$, the initial elastic
loading front becomes more distant to the plastic loading front,
creating the so-called elastic precursor, which is a typical feature of
stress wave propagation induced by laser impact. This precursor can be
experimentally observed in backface velocity measurements
\cite{ayad_new_2022,seddik_identification_2022,ayad_modeling_2022}.

As a complement to Figures~\ref{fig:s_v_ep_2GPa}
and~\ref{fig:s_v_ep_5GPa}, Table~\ref{tab:errors} presents the average
relative errors between the numerical schemes and analytical
expressions. These errors are computed using the profiles at $t=490$ ns
and $t=790$ ns for the cases $\sigma_{\mathrm{hel}}<P\leq 2\sigma_{\mathrm{hel}}$ and 
{$P>(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$} respectively. The relative error
is obtained by normalizing the absolute error with the highest value
(in absolute value) of each corresponding profile.
Table~\ref{tab:errors} confirms that the GHR scheme is the most
precise. Note that, while the errors are of the same order of
magnitude, the results from the ABAQUS software give higher errors than
the other schemes in most cases.

%tab6
\begin{table}
\caption{Average relative errors between the exact solutions and the
numerical schemes, computed along the length $L$, at $t=490$ sn in the
case $\sigma_{\mathrm{hel}}<P\leq 2\sigma_{\mathrm{hel}}$ and at
$t=790$ ns in the case  $P>(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$
\label{tab:errors}}
\tabcolsep 4pt\begin{tabular}{ccccccc}
\thead
& \multicolumn{3}{c}{Case $\sigma_{\mathrm{hel}}<P\leq 2\sigma_{\mathrm{hel}}$} &
\multicolumn{3}{c}{Case  {$P>(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$}}\vspace*{2pt} \\\cline{2-4}\cline{5-7}
 & $\sigma$ & $v$ & $\varepsilon_p$ & $\sigma$ & $v$ & $\varepsilon_p$\vspace*{2pt} \\ 
\endthead
 Lax--Friedrichs & $4.77\times 10^{-3}$ & $4.86\times 10^{-3}$ & $2.16\times 10^{-2}$ & $3.75\times 10^{-1}$ & $3.77\times 10^{-1}$ & $2.48\times 10^{-2}$\raisebox{1pc}{} \\
 Lax--Wendroff & $3.11\times 10^{-3}$ & $4.52\times 10^{-3}$ & $2.09\times 10^{-2}$ & $3.74\times 10^{-1}$ & $3.73\times 10^{-1}$ & $1.42\times 10^{-2}$ \\
 Godunov & $3.40\times 10^{-3}$ & $3.41\times 10^{-3}$ & $1.51\times 10^{-2}$ & $3.74\times 10^{-1}$ & $3.74\times 10^{-1}$ & $1.88\times 10^{-2}$ \\
 GHR & $8.92\times 10^{-4}$ & $9.16\times 10^{-4}$ & $7.93\times 10^{-3}$ & $3.70\times 10^{-1}$ & $3.70\times 10^{-1}$ & $6.41\times 10^{-3}$ \\
 ABAQUS & $6.52\times 10^{-3}$ & $6.95\times 10^{-3}$ & $8.87\times 10^{-3}$ & $3.66\times 10^{-1}$ & $3.66\times 10^{-1}$ & $2.62\times 10^{-2}$ 
\botline
\end{tabular}
\vspace*{6pt}
\end{table}

Most of the exact results presented in this paper are in agreement with
work of Ballard~\cite{ballard_contraintes_1991}. Namely, in the
$\sigma_{\mathrm{hel}}<P\leq 2\sigma_{\mathrm{hel}}$ case, the value of the induced
plastic strain (see Table~\ref{tab:epsp_2GPa}) linearly depends on the
applied pressure $P$. In the  {$P>(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$}
case, the initial plastic loading induces a level of plastic strain
that is independent of the applied pressure, and is thus considered as
saturated (see line 2 of Table~\ref{tab:epsp_5GPa}).
Ballard~\cite{ballard_contraintes_1991} also presents exact results for the
saturated $l_{\mathrm{sat}}$ and affected depth $l_{\mathrm{plast}}$. The plastic strain
is considered saturated when its value does not depend on the applied
pressure, but only on the material parameters (see~line~2 
of~Table~\ref{tab:epsp_5GPa}). The saturated depth thus corresponds to the
depth along which the plastic strain is saturated. The affected depth
corresponds to the depth along which plastic strain is induced,
regardless of its saturated or unsaturated nature. The affected depth
is larger than the saturated depth. Their expressions are the
following, without any particular consideration on the value of~$P$:
{\begin{equation}
l_{\mathrm{sat}}=\frac{c_p c_e t_0}{c_e-c_p}\left\lfloor \frac{P}{2\sigma_{\mathrm{hel}}} \right\rfloor 
\quad \tn{and} \quad l_{\mathrm{plast}}=\frac{c_p c_e t_0}{c_e-c_p}\left\lfloor 
\frac{P+\sigma_{\mathrm{hel}}}{2\sigma_{\mathrm{hel}}} \right\rfloor,
\label{eq:l_sat_plast_ballard}
\end{equation}}\noindent
where the notation $\left\lfloor \cdot \right\rfloor$ denotes the floor
function. In our case we find, using the $x$--$t$ diagram of
Figure~\ref{fig:diag_lag_5GPa}:
{\begin{equation}
l_{\mathrm{sat}}=\frac{c_p c_e t_0}{c_e-c_p} \quad \tn{and} \quad l_{\mathrm{plast}}=2l_{\mathrm{sat}},
\end{equation}}\noindent
which is in agreement with Equation \eqref{eq:l_sat_plast_ballard}.
However we show that the plastic strain, at the surface, is saturated
only for a limited period of time. Indeed, the interactions between the
plastic and elastic fronts create a smaller secondary wave
(corresponding to zone (g) in Figure~\ref{fig:diag_lag_5GPa}(a)) that
is propagating in the opposite direction as the initial wave. This
secondary wave eventually reaches the left free surface boundary of the
domain ($x=0$), at which point it is reflected with an opposite sign.
The opposite sign makes the wave generate new plastic strain in a zone
where it was previously saturated. The plastic strain at the surface is
saturated for $t\in[t_1,t_2]$, with
{\begin{equation}
t_1=t_0 \quad \tn{and} \quad t_2=\frac{t_0\left(c_e^2+c_p^2 \right)}{c_p(c_e-c_p)}.
\end{equation}}\noindent

The new plastic strain taking place over the saturated plastic strain
is dependent on the applied pressure (see lines 6 and 8 of
Table~\ref{tab:epsp_5GPa}), and is thus in this sense not saturated.
These results are due to the perfectly plastic nature of the material
model and may change if one considers hardening.

The present paper also improves the work of
Ballard~\cite{ballard_contraintes_1991}, mainly concerning the structure of
the $x$--$t$ diagram of Figure~\ref{fig:diag_lag_5GPa}(a) and the plastic
strain distribution. We present here the full description of the stress
wave propagation for $P=5$ GPa, with complex reflection patterns,
leading to the plastic strain being modified in several zones after the
initial plastic loading, more specifically by the loading fronts
identified by the numbers 4, 5, 6, 7 and 8 in
Figure~\ref{fig:diag_lag_5GPa}(a). Even though the work of
Heuz\'e~\cite{heuze_laxwendroff_2017} tackles the same problem, it is modeled
using a bar model, which is different from our uniaxial strain
hypothesis, and does not result in the same stress wave propagation.
These particularities of the stress wave propagation are confirmed by
the numerical results, especially the GHR scheme (see for instance
Figure~\ref{fig:s_v_ep_5GPa}(e) and (f)). Even if these features do not
change the plastic affected depth, they change the plastic strain
distribution, which eventually influences the residual stress
distribution. Similar considerations in a general case on the stress
and velocity can be found in Wang~\cite{wang_foundations_2007} (see 
Figure~7.5 of their work), without however any quantitative information
about the plastic strain values.

The results of the case  {$P>(c_e+5c_p)/(c_e+c_p)\sigma_{\mathrm{hel}}$} are
however not valid for any $P$. If $P$ increases too much, the structure
of the corresponding $x$--$t$ diagram will change, and with it the whole
exact solution. The situation will be different if the wave front
between zones (f) and (l) becomes fully plastic. This will be true if
$\sigma_{(f)}>2\sigma_{\mathrm{hel}}$ (see Table~\ref{tab:diag_lag_5GPa}),
leading to:
{\begin{equation}
    P>2\sigma_{\mathrm{hel}}\frac{c_e+3c_p}{c_e+c_p}.
    \label{eq:cond_P}
\end{equation}}\noindent

For the sake of simplicity, we chose not to include the analysis of
this case in this paper, which already covers most of the pressure
range reachable in LSP for the chosen material. However, the change of
structure of the exact solution implied by Equation \eqref{eq:cond_P}
was verified numerically. Ballard~\cite{ballard_contraintes_1991} also
proposed a pressure limit above which the material starts to behave
hydrodynamically, where a simple elastic behavior is not enough to
describe the material, which needs to make use of an equation of state
\cite{peyre_analysis_2012,ayad_modeling_2022}. This condition reads:
{\begin{equation}
P\lesssim 0.1\kappa.
\end{equation}}\noindent

In our case, $0.1\kappa=17.5$ GPa, while the maximum applied pressure
was 5 GPa.

{The plastic strain profiles of Figures~\ref{fig:s_v_ep_2GPa}(e)
and~\ref{fig:s_v_ep_2GPa}(f) present large values close to $x=0$
diverging from the exact solution. These plastic strain values can be
linked to the sharp nature of the loading profile, which induces large
stress} \aj{time}  {increments in the very first steps of the
simulation, and thus large plastic strain} \aj{time}  increments.
Later in the simulation, \aj{the plastic
fronts become less sharp or smaller than the one in the initial
instants,}  {leading to smaller} \aj{stress time increments and thus to
smaller}  {plastic strain} \aj{time}  {increments} \aj{in accordance
with Equation \eqref{eq:delt_epsp_sig}.}  {To illustrate the influence
of the sharpness of the pressure time profile, we can compute the
plastic strain field induced by a pressure time profile of amplitude 2
GPa (to replicate the conditions of Figure~\ref{fig:s_v_ep_2GPa}),
whose edges have been smoothed. The results are computed for the
Lax--Wendroff and Godunov with High Resolution schemes, since they are
the ones with the least amount of artificial viscosity. The results are
presented in Figure~\ref{fig:ep_sharp_smooth}. For both schemes, the
plastic strains close to the edge of the domain are reduced once a
smoother profile is used, though not entirely eliminated. It could be
achieved by increasing the sampling of the time profile, which however
implies to reduce the time increment of the simulation in order to get
more point in the rise and decrease part of the profile. Here the
discretization of the time profiles in
Figure~\ref{fig:ep_sharp_smooth}(a) corresponds exactly to the time
increment $\Delta t$ used in the simulations. However the smoothed
version of the time profiles presents the disadvantage of producing
stress waves that are more quickly attenuated, resulting in lower
plastically affected depths. This makes the comparison with the exact
solution less relevant, since the latter assume a instantaneous loading
of the material. }

\begin{figure}
\includegraphics{fig11}
\caption{Comparisons of the plastic strain fields obtained with a sharp
and smooth pressure time profile. (a) Different pressure time profiles.
(b) Plastic strain fields computed with the Lax--Wendroff scheme. (c)
Plastic strain fields computed with the Godunov with High Resolution
scheme. \label{fig:ep_sharp_smooth}}
\end{figure}

The exact solution presented in this work relies on some necessary
hypotheses (uniaxial strain field, perfectly plastic material, pressure
loading as a step function) that may not be met in actual LSP
configuration. Nonetheless, this solution provides insight and
explanations to phenomenon that can be observed during LSP. The
solution is also useful from a numerical point of view, to assess the
quality of a numerical solution, before it is applied to more realistic
modeling hypotheses. An important hypothesis in this work that allowed
the computation of analytic stress and velocity fields is the step time
pressure loading (limiting the number of characteristic lines needed
for the $x$--$t$ diagrams). Actual experimental pressure time profiles are
smoother \cite{scius-bertrand_laser_2021}, and the exact fields are
likely to overestimate the actual fields. 

Finally, it is worth reminding that the expressions of the final
displacement of the surface of the specimen in Equations
\eqref{eq:ux_2GPa} and \eqref{eq:ux_5GPa} are valid in this uniaxial
case, but may not be directly comparable to experimental (e.g.
profilometry) or even 3D numerical results. Indeed, because of the edge
effects originating at the edge of the laser spot and propagating
inward \cite{ayad_modeling_2022,ecault_numerical_2016}, a secondary
velocity wave passes through the center of the impact, modifying the
displacement field. This phenomenon cannot be represented with the
current uniaxial model, and thus discrepancies can be expected when
comparing its results to experimental data.

\section{Conclusion}

This article tackled the exact solution for the stress and velocity
wave induced by laser shock peening, as well as the induced plastic
strain and residual stress fields. This paper extends the work of
Ballard~\cite{ballard_contraintes_1991}. To achieve this, the elasto-plastic
stress wave equation was first derived for uniaxial strains in an
elastic perfectly plastic material. It was then showed that the
hyperbolic nature of the stress wave equation allowed the determination
of an exact solution using the methods of characteristics. It also
allows the use of several dedicated numerical schemes to compute a
numerical solution, each with its own numerical properties. Comparisons
between the exact and numerical solutions were made for two different
loading amplitudes, each one leading to a different architecture of the
stress wave. It was shown that even though all numerical methods were
able to correctly reproduce the exact solution, only the Godunov with
High Resolution was able to capture the finest feature of the exact
results. The plastic strains results were used to estimate the residual
stress field, with an exact analytic formulation and with numerical
results, both being in agreement. This work can thus be used to compute
an initial estimation of the mechanical fields without any need of
numerical computation. It can also be used to verify the precision of
different numerical methods. An improvement of this work would be to
determine the next analytical configuration for the $x$--$t$ diagram in
case the pressure is much higher than the Hugoniot elastic limit, but
being still realistic for LSP application, and without stepping into
the hydrodynamic regime. Applying the method for smoother loading
profiles is also possible, but the latter must be discretized into
several step-wise constant portions.

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

\section*{Acknowledgments}

This research was partly funded by the Carnot Institute ARTS and by the
French Agence Nationale de la Recherche (ANR), research project
ForgeLaser (grant number: ANR-18-CE08-0026).

\CDRGrant[ANR]{ANR-18-CE08-0026}

\back{}

\appendix{}

\def\appendixlabel{Appendix \Alph{section}}

\renewcommand\theequation{A\arabic{equation}}

\section{Exact expressions of the parameters $\alpha$ and $\beta$ for
the evaluation of the residual stresses in the case
$P>\sigma_{\mathrm{hel}}$}\label{sec:app_alpha_beta}
{\begin{equation}
    \left\{
    \begin{array}{ll}
       \begin{aligned}
        \ds \alpha= & \frac{1}{L^3\kappa(c_e-c_p)^2(c_e+c_p)^7}\left[ c_{p} t_{0} \left( L \left[ P \left( - c_{e}^{9} - 7 c_{e}^{8} c_{p} - 20 c_{e}^{7} c_{p}^{2} - 28 c_{e}^{6} c_{p}^{3} - 14 c_{e}^{5} c_{p}^{4} +
        14 c_{e}^{4} c_{p}^{5}  \right. \right. \right. \right. \vspace*{0.35cm}\\
        &  \left. \left. \left.+\, 28 c_{e}^{3} c_{p}^{6} + 20 c_{e}^{2} c_{p}^{7} + 7 c_{e} c_{p}^{8} + c_{p}^{9} \right) + 2\sigma_{\mathrm{hel}} \left( 3 c_{e}^{9} - 15 c_{e}^{8} c_{p} + 20 c_{e}^{7} c_{p}^{2} + 60 c_{e}^{6} c_{p}^{3} + 2 c_{e}^{5} c_{p}^{4} - 10 c_{e}^{4} c_{p}^{5}  \right. \right. \right. \vspace*{0.35cm}\\
        & \left. \left.-\, 12 c_e^3c_p^6 - 36 c_{e}^{2} c_{p}^{7} - 13 c_{e} c_{p}^{8} + \left. \left. c_{p}^{9} \right) \right] + t_0 \left[ P \left( 5 c_{e}^{9} c_{p} + 36 c_{e}^{8} c_{p}^{2} + 112 c_{e}^{7} c_{p}^{3} + 196 c_{e}^{6} c_{p}^{4} + 210 c_{e}^{5} c_{p}^{5}  \right. \right. \right. \right. \vspace*{0.35cm}\\
        & \left. \left. \left.+\, 140 c_{e}^{4} c_{p}^{6} + 56 c_{e}^{3} c_{p}^{7} + 12 c_{e}^{2} c_{p}^{8} + c_{e} c_{p}^{9} \right) + 2 \sigma_{\mathrm{hel}}\left( - 4 c_{e}^{10} + 23 c_{e}^{9} c_{p} - 102 c_{e}^{8} c_{p}^{2} + 10 c_{e}^{7} c_{p}^{3} -  342 c_{e}^{6} c_{p}^{4}  \right. \right. \right. \vspace*{0.35cm}\\
        & \left. \left. \left. \left.-\, 280 c_{e}^{5} c_{p}^{5} - 146 c_{e}^{4} c_{p}^{6} - 138 c_{e}^{3} c_{p}^{7} - 46 c_{e}^{2} c_{p}^{8} + c_{e} c_{p}^{9} \right) \right] \right) \right] 
        \end{aligned}\vspace*{5pt}\\
        \begin{aligned}
        \ds \beta = & \frac{1}{6L^2\kappa(c_e-c_p)^2(c_e+c_p)^7}\left[ - c_{p} t_{0} \left( L \left[ 4P \left( - c_{e}^{9} - 7 c_{e}^{8} c_{p} - 20 c_{e}^{7} c_{p}^{2} - 28 c_{e}^{6} c_{p}^{3} - 14 c_{e}^{5} c_{p}^{4} + 14 c_{e}^{4} c_{p}^{5}  \right. \right. \right. \right. \vspace*{0.35cm}\\
        & \left. \left. \left.+\, 28 c_{e}^{3} c_{p}^{6} + 20 c_{e}^{2} c_{p}^{7} + 7 c_{e} c_{p}^{8} + c_{p}^{9} \right) + 8\sigma_{\mathrm{hel}} \left( 3 c_{e}^{9} - 15 c_{e}^{8} c_{p} + 20 c_{e}^{7} c_{p}^{2} + 60 c_{e}^{6} c_{p}^{3} + 2 c_{e}^{5} c_{p}^{4} - 10 c_{e}^{4} c_{p}^{5}  \right. \right. \right. \vspace*{0.35cm}\\ 
        & \left. \left.-\, 12 c_{e}^{3} c_{p}^{6} - 36 c_e^3c_p^7 - 13 c_{e} c_{p}^{8} + \left. \left. c_{p}^{9} \right) \right] + t_0 \left[ P \left(  15 c_{e}^{9} c_{p} + 108 c_{e}^{8} c_{p}^{2} + 336 c_{e}^{7} c_{p}^{3} + 588 c_{e}^{6} c_{p}^{4} + 630 c_{e}^{5} c_{p}^{5}    \right. \right. \right. \right. \vspace*{0.35cm}\\
        & \left. \left. \left.+\, 420 c_{e}^{4} c_{p}^{6} + 168 c_{e}^{3} c_{p}^{7} + 36 c_{e}^{2} c_{p}^{8} + 3 c_{e} c_{p}^{9} \right) + 6 \sigma_{\mathrm{hel}} \left( - 4 c_{e}^{10} + 23 c_{e}^{9} c_{p} - 102 c_{e}^{8} c_{p}^{2} + 10 c_{e}^{7} c_{p}^{3}   \right. \right. \right. \vspace*{0.35cm}\\ 
        & \left. \left. \left. \left. -\, 342 c_{e}^{6} c_{p}^{4}- 280 c_{e}^{5} c_{p}^{5} - 146 c_{e}^{4} c_{p}^{6} - 138 c_{e}^{3} c_{p}^{7} - 46 c_{e}^{2} c_{p}^{8} + c_{e} c_{p}^{9} \right) \right] \right) \right] .
        \end{aligned}
    \end{array}
    \right.
\end{equation}}\noindent

\section{Expression of the fluxes for the numerical
schemes}\label{app:num_schemes}

\renewcommand\theequation{B\arabic{equation}}

We detail here the expression of the fluxes of Equation
\eqref{eq:adv_flux}. They can be found in Leveque~\cite{leveque_finite_2002}.

\subsection{Lax--Friedrichs} 
The Lax--Friedrichs flux reads:
{\begin{equation}
{\mathbi{F}_{i-1/2}^{n}=\frac{1}{2}\mathbi{A}\left(\mathbi{U}^n_{i-1}+\mathbi{U}^n_{i} \right)-
\frac{\Delta x}{2\Delta t}\left( \mathbi{U}^n_i -\mathbi{U}^n_{i-1}\right)}
\end{equation}}\noindent

However to ensure stability, the complete scheme must be:
{\begin{equation}
\mathbi{U}_i^{n+1}=\frac{1}{2}\left(\mathbi{U}^n_{i+1}+\mathbi{U}^n_{i-1} \right)-\frac{\Delta t}{2\Delta x} 
{\mathbi{A}}\left(\mathbi{U}^n_{i+1}-\mathbi{U}^n_{i-1} \right).
\end{equation}}\noindent

\subsection{Lax--Wendroff} 
The Lax--Wendroff flux reads:
{\begin{equation}
\mathbi{F}_{i-1/2}^{n}=\frac{1}{2}\mathbi{A}\left(\mathbi{U}^n_{i-1}+\mathbi{U}^n_i \right)-\frac{1}{2}
\frac{\Delta t}{\Delta t}\mathbi{A}^2\left(\mathbi{U}^n_i-\mathbi{U}^n_{i-1} \right).
\end{equation}}\noindent

\subsection{Godunov} The expression for the Godunov flux is:
{\begin{equation}
\mathbi{F}^{n}_{i-1/2}=\mathbi{A}^+\mathbi{U}^n_{i-1}+\mathbi{A}^-\mathbi{U}^n_i.
\end{equation}}\noindent

Matrices $\mathbi{A}^+$ and $\mathbi{A}^-$ are obtained in the following manner:
{\begin{equation}
\mathbi{A}^+=\mathbi{P}\mathbi{D}^+\mathbi{P}^{-1} \quad \tn{and} \quad \mathbi{A}^-=\mathbi{P}\mathbi{D}^-\mathbi{P}^{-1},
\end{equation}}\noindent
with $\mathbi{D}^+$ and $\mathbi{D}^-$ being diagonal matrices such
that $(\mathbi{D}^+)_{ii}=\max( \gamma^{(i)},0)$ and
$(\mathbi{D}^-)_{ii}=\min( \gamma^{(i)},0)$. $\gamma^{(i)}$ are the
eigenvalues of $\mathbi{A}$ and $\mathbi{P}$ the matrix containing the
eigenvectors of $\mathbi{A}$.

\subsection{Godunov with high resolution}

The Godunov with High Resolution scheme is similar to the standard
Godunov scheme:
{\begin{equation}
\mathbi{F}^{n}_{i-1/2}=\mathbi{A}^+\mathbi{U}^n_{i-1}+\mathbi{A}^-\mathbi{U}^n_i+\widetilde{\mathbi{F}}^n_{i-1/2}.
\end{equation}}\noindent

The expression of $\widetilde{\mathbi{F}}^n_{i-1/2}$ is:
{\begin{equation}
\ds \widetilde{\mathbi{F}}^n_{i-1/2}=\frac{1}{2}\sum_p \left\lvert \gamma^{(p)} \right\rvert
\left( 1-\frac{\Delta t}{\Delta x}\left\lvert \gamma^{(p)} \right\rvert \right)
\widetilde{\alpha}^{(p)}_{i-1/2}\mathbi{r}^{(p)},
\end{equation}}\noindent
where $\mathbi{r}^{(p)}$ is the $p{\mathrm{th}}$ eigenvector of $\mathbi{A}$, and
${\widetilde{\alpha}}^{(p)}_{i-1/2}$ is given by:
{\begin{equation}
\widetilde{\alpha}^{(p)}_{i-1/2}=\alpha^{(p)}_{i-1/2}\phi ( \theta^{(p)}_{i-1/2}),
\end{equation}}\noindent
with $\alpha^{(p)}_{i-1/2}=(\mathbi{P}^{-1}(\mathbi{U}_i-\mathbi{U}_{i-1}
))_p$. In this equation, $\phi$ is a function called a flux
limiter. In this work we consider the Superbee flux limiter which
reads:
{\begin{equation}
\phi(\theta)=\max\left(0,\min (1,2\theta),\min(2,\theta) \right).
\label{eq:superbee}
\end{equation}}\noindent

Finally in our case:
{\begin{equation}
\theta^{(p)}_{i-1/2}=\frac{\alpha^{(p)}_{I-1/2}}{\alpha^{(p)}_{i-1/2}},
\end{equation}}\noindent
with
{\begin{equation}
    I=\left\{
    \begin{array}{lr}
        i-1 & \quad\text{if } \gamma^{(p)}>0 \\
        i+1 & \quad\text{if } \gamma^{(p)}<0 
    \end{array}.
    \right.
\end{equation}}\noindent

\bibliographystyle{crunsrt}
\bibliography{crmeca20230403}
\refinput{crmeca20230403-reference.tex} 

\end{document}
