\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRMECA,Unicode,published]{cedram}
\newenvironment{noXML}{}{}%
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.935\linewidth}{#1}}
\usepackage{etoolbox}
\usepackage{upgreek}
\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}}
\let\MakeYrStrItalic\relax
\def\jobid{crmeca20221189}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
\skip\footins=10pt
\def\0{\phantom{0}}
\def\ndash{\text{--}}
\def\rmmu{\upmu}
\def\mn{\phantom{$-$}}
\def\refinput#1{}
\def\back#1{}
\def\og{\guillemotleft}
\def\fg{\guillemotright}
\DOI{10.5802/crmeca.230}
\datereceived{2022-12-05}
\daterevised{2023-01-27}
\datererevised{2023-06-12}
\dateaccepted{2023-10-09}
\ItHasTeXPublished
\setcounter{secnumdepth}{4}
}
{\documentclass[crmeca]{article}
\usepackage[T1]{fontenc} 
\usepackage{upgreek}
\def\CDRdoi{10.5802/crmeca.230}
\def\newline{\unskip\break}
}
\makeatother
\begin{DefTralics}
\newcommand{\tsup}[1]{#1}
\end{DefTralics}

\usepackage{upgreek}

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

\begin{noXML}

\RubricEF{History of Sciences and Ideas}{Histoire des sciences et des id\'ees}

\title{Implementation of the c-phi reduction procedure in Cast3M code
for calculating the stability of retaining walls in the layered
backfill with strength parameters reduction by elasto-plastic finite
element analysis using fields data}

\alttitle{Mise en {\oe}uvre de la proc\'{e}dure c-phi de r\'{e}duction
des param\`{e}tres de r\'{e}sistance dans le code Cast3M pour le calcul
de la stabilit\'{e} des murs de sout\`{e}nement dans le remblai
stratifi\'{e} par analyse elements finis elasto-plastique en utilisant
les donn\'{e}es de terrain}

\author{\firstname{Zoa} \lastname{Ambassa}\CDRorcid{0000-0002-0359-695X}\IsCorresp}
\address{Laboratory of Energy Modeling Materials and Methods (E3M),
National Higher Polytechnic School of Douala, University of Douala,
P.O. Box: 2701 Douala, Cameroon}
\address{Department of Civil Engineering, National Higher Polytechnic
School of Douala, University of Douala, Cameroon}
\email[Z. Ambassa]{zoa.amabassa@enspd-udo.cm}

\author{\firstname{Jean Chills} \lastname{Amba}\CDRorcid{0000-0002-5801-1349}}
\addressSameAs{1}{Laboratory of Energy Modeling Materials and Methods
(E3M), National Higher Polytechnic School of Douala, University of
Douala, P.O. Box: 2701 Douala, Cameroon}
\addressSameAs{2}{Department of Civil Engineering, National Higher
Polytechnic School of Douala, University of Douala, Cameroon}
\email[J. C. Amba]{jc.amba@enspd-udo.cm}

\author{\firstname{Nandor} \lastname{Tamaskovics}}
\address{TU Bergakademie Freiberg, Geotechnical Institute, Chair of
Soil Mechanics and Ground Engineering, Gustav-Zeuner-Stra{\ss}e 1/Room
014/Rom 016, Germany}
\email[N. Tamaskovics]{Nandor.Tamaskovics@ifgt.tu-freiberg.de}

\shortrunauthors

\keywords{\kwd{Implementation}\kwd{C-phi reduction procedure}\kwd{Code
Cast3M}\kwd{Stability factor}\kwd{Geotechnical structure}\kwd{Finite
element calculations}}

\altkeywords{\kwd{Mise en {\oe}uvre}\kwd{Proc\'{e}dure de r\'{e}duction
c-phi}\kwd{Code Cast3M}\kwd{Facteur de stabilit\'{e}}\kwd{Structure
g\'{e}otechnique}\kwd{Calculs par \'{e}l\'{e}ments finis}}

\begin{abstract}
This paper presents an implementation of the c-phi reduction procedure
in the finite element code Cast3M. This procedure is first validated on
a simple example of bearing capacity of a strip footing and then used
to evaluate the stability factor of a geotechnical structure composed
of two levels of backfill retained by two retaining walls, with a road
built on the top of the second backfill. The results of the c-phi
reduction procedure for the construction of first and second backfill
are also compared to other analysis methods classically used in
geotechnics and to the results of finite element calculations.
\end{abstract}

\begin{altabstract}
Cet article pr\'{e}sente une impl\'{e}mentation de la proc\'{e}dure de
r\'{e}duction c-phi dans le code aux \'{e}l\'{e}ments finis Cast3M.
Cette proc\'{e}dure est d'abord valid\'{e}e sur un exemple simple de
capacit\'{e} portante d'une semelle filante, puis utilis\'{e}e pour
\'{e}valuer le facteur de stabilit\'{e} d'une structure
g\'{e}otechnique compos\'{e}e de deux niveaux de remblai retenus par
deux murs de sout\`{e}nement, avec une route construite sur le dessus
du second remblai. Les r\'{e}sultats de la proc\'{e}dure de
r\'{e}duction c-phi pour la construction du premier et du second
remblai sont \'{e}galement compar\'{e}s \`{a} d'autres m\'{e}thodes
d'analyse classiquement utilis\'{e}es en g\'{e}otechnique et aux
r\'{e}sultats des calculs par \'{e}l\'{e}ments finis.
\end{altabstract}

\maketitle
\vspace*{-.5pc}

\end{noXML}

\section{Introduction}

For decades, the majority of slope stability analyses performed in
practice still use traditional limit equilibrium approaches involving
methods of slices that have remained essentially unchanged. This was
not the outcome envisaged when Whitman and Bailey~\cite{1} set criteria
for the then emerging methods to become readily accessible to all
engineers. The Finite Element Method represents a powerful alternative
approach for slope stability and retained wall analysis which is
accurate, versatile and requires fewer {\it a priori} assumptions,
especially, regarding the failure mechanism. Slope failure in the
Finite Element model occurs ``naturally'' through the zones in which
the shear strength of the soil is insufficient to resist the shear
stresses. Elasto-plastic analysis of geotechnical problems using the
Finite Element (FE) Method has been widely accepted in the research
arena for many years; however, its routine use in geotechnical practice
for slope and retaining stability analysis still remains limited. The
reason for this lack of acceptance is not entirely clear; however,
advocates of FE techniques in academe must take some responsibility.
Practising engineers are often sceptical of the need for such
complexity, especially in view of the poor quality of soil property
data often available from routine site investigations~\cite{2,3,4}.
Although this scepticism is often warranted, there are certain types of
geotechnical problem for which the FE approach offers real benefits.
The challenge for an experienced engineer is to know which kind of
problem would benefit from a FE treatment and which would not. In
general, linear problems such as the prediction of settlements and
strains, the calculation of flow quantities due to steady seepage or
the study of transient effects due to consolidation are all highly
amenable to solution by Finite Elements. Traditional approaches
involving charts, tables or graphical methods will often be adequate
for routine problems but the Finite Elements approach may be valuable
if awkward geometries or material variations are encountered which are
not covered by traditional chart solutions. 

The use of nonlinear analysis in routine geotechnical practice is
harder to justify, because there is usually a significant increase in
complexity which is more likely to require the help of a modelling
specialist. Nonlinear analyses are inherently iterative in nature,
because the material properties and/or the geometry of the problem are
themselves a function of the ``solution''. Objections to nonlinear
analyses on the grounds that they require excessive computational
power, however, have been largely overtaken by developments in, and
falling costs of, computer hardware. Slope stability and retained wall
represents an area of geotechnical analysis in which a nonlinear FE
approach offers real benefits over existing methods. As this paper will
show, Slope stability and retained wall analysis by elastoplastic
Finite Elements is accurate, robust and simple enough for routine use
by practising engineers. Perception of the FE method as complex and
potentially misleading is unwarranted and ignores the real possibility
that misleading results can be obtained with conventional ``slip
circle'' approaches. The graphical capabilities of FE programs also
allow better understanding of the mechanisms of failure, simplifying
the output from reams of paper to manageable graphs and plots of
displacements.

The objective of this work is the development of an elasto-plastic
calculation method for the prediction of the stability for the large
geotechnical structures in interaction by Finite Element technique
thanks to with methodological and sophisticated programming from Cast3M
code. This work does not only concern the slope stability calculation,
but the stability calculation of bearing capacity, overturning and
sliding of the layered backfill and the retaining walls by the
\mbox{c-phi} parameters reduction because it concerns several
structures in interaction that encompass the failure mechanisms of the
majority of geotechnical structures. The paper describes several
standards in geotechnical practice of Finite Element retained wall and
slope stability analysis with comparison against other solution
methods, including the influence of a free surface on slope and layered
backfill stability. Graphical output is included to illustrate
displacements, deformations, stresses and mechanisms of failure.

\section{Traditional methods of slope stability analysis} \label{sec2}

Most textbooks on soil mechanics or geotechnical engineering will
include reference to several alternative methods of slope stability
analysis. In a survey of equilibrium methods of slope stability
analysis reported by Duncan~\cite{5}, the characteristics of a large
number of methods were summarized, including the ordinary method of
slices~\cite{6}, Bishop's Modified Method~\cite{7}, force equilibrium
methods (e.g. Lowe and Karafiath)~\cite{8}, Janbu's generalized
procedure of slices~\cite{9}, Morgenstern and Price's
method~\cite{10,11} and Spencer's method~\cite{12}. Although there
seems to be some consensus that Spencer's method is one of the most
reliable, textbooks continue to describe the others in some detail, and
the wide selection of available methods is at best confusing to the
potential user. For example, the controversy was revisited by Lambe and
Silva~\cite{13}, who maintained that the ordinary method of slices had
an undeservedly bad reputation. A difficulty with all the equilibrium
methods is that they are based on the assumption that the failing soil
mass can be divided into slices. This in turn necessitates futher
assumptions relating to side force directions between slices, with
consequent implications for equilibrium. The assumption made about the
side forces is one of the main characteristics that distinguishes one
limit equilibrium method from another, and yet is itself an entirely
artificial distinction.

\section{Finite elements method for slope and retaining walls stability
analysis} \label{sec3}

Duncan's review of FE analysis of slopes concentrated mainly on
deformation rather than stability analysis of slopes; however,
attention was drawn to some important early papers in which
elasto-plastic soil models were used to assess stability. Smith and
Hobbs~\cite{14} reported results of $\varphi_u = 0$ slopes and obtained
reasonable agreement with Taylor's~\cite{15} charts. Zienkiewicz
{\etal}~\cite{16} considered a $c'$, $\varphi'$ slope and obtained good
agreement with slip circle solutions. Griffiths~\cite{17,18} extended
this work to show reliable slope stability results over a wide range of
soil properties and geometries as compared with charts of Bishop and
Morgenstern~\cite{19}. Subsequent use of the FE method in slope
stability analysis has added further confidence in the
method~\cite{20,21,22,23,24,25,26,27,28,29,30}. Griffiths~\cite{20}
investigate the probability of failure of a cohesive slope using both
simple and more advanced probabilistic analysis tools. The influence of
local averaging on the probability of failure of a test problem has
been investigated. The advanced method, called the random finite
element method (RFEM), uses elastoplasticity combined with random field
theory. For the authors, the RFEM method is shown to offer many
advantages over traditional probabilistic slope stability techniques,
because it enables slope failure to develop naturally by seeking out
the most critical mechanism. Griffiths~\cite{20} reported that
simplified probabilistic analysis, in which spatial variability is
ignored by assuming perfect correlation, lead to unconservative
estimates of the probability of failure. To avoid human error in the
analyses criteria and find a unified mechanical explanation for slope
instability, a criterion has been proposed by Chengya~\cite{21}. For
the author, this criterion, of the stability is determined by
identifying the positive or negative of the variational value. The
Factor of safety of two and three dimensional slopes calculated by this
criterion are compared with those calculated by other criteria. The
effect of mesh size on the accuracy of the results have been analyzed.
The study made by Anindita~\cite{22} proposes a comparison of Limit
Equilibrium and Finite Element methods by numerical approach. The
method developed by Dakshith~\cite{23} presents a numerical technique
combining the scaled boundary finite element method and image-based
meshing for slope stability analysis. Quad tree decomposition has been
applied to simultaneously generate meshes and consider the spatial
variation of material properties directly from the images through a
mapping technique. The shear strength technique has been applied to
evaluate the shear reduction factor of safety of the slope. Coal slopes
at Yallourn open-pit mine, Victoria, Australia was considered, forming
the basis of a case study to demonstrate the applicability of the
method. Huang~\cite{24} extended the use of strength reduction FEM to
include the effects of unsaturated transient seepage and some primary
numerical results concerning the stability of an earth dam under rapid
drawdown have been presented. A stability analysis performed on the
Yashigou earth dam in China was conducted by Huy~\cite{25}. Factor of
safety results were considered for vacuum dam, steady state water level
and water level drawdown using finite element stresses method analysis
to compare the Morgenstern--Price method. Jiang~\cite{26} has developed
a slope stability analysis considering spatially variable shear
strength parameters using a non-intrusive stochastic finite element
method. Zoa~\cite{27} conducted a finite element analysis of a soft
soil reinforced by an industrial nailing technique using a soft soil
behavior model. He determined the main influential parameters. The
probabilistic stability analysis of an earth dam by stochastic finite
element method based on field data has been performed by
Mouyeaux~\cite{28} for realistic prediction of the failure mechanism
and seepages under this dam. Recently, a hybrid stacking ensemble
approach Navid~\cite{29} for enhancing the prediction of slope
stability. In this hybrid stacking ensemble, he used an artificial bee
colony algorithm to find out the best combination of base classifiers
(level~0) and determined a suitable meta-classifier (level~1) from a
pool of 11 individual optimized machine learning algorithm. Finite
element analysis was conducted in order to form the synthetic database
for the training stage (150 cases) of proposed model while 107 real
field slope cases were used for the testing stage. The results by the
hybrid stacking ensemble were compared with that obtained by the 11
optimized machine learning. The comparisons showed that a significant
improvement in the prediction ability of slope stability has been
achieved by the hybrid stacking ensemble. A numerical model termed as
ISR3DNMM-GPS method has been developed by Yantao~\cite{30} for factor
of safety calculation and 3D failure surface determination involved in
a 3D slope stability analysis. In this proposed numerical model, an
improved strength reduction technique (ISRT) which can avoid
unreasonable plastic zones appearing in the deep area of a slope has
been implemented into 3D numerical manifold method, and improved
strength reduction based 3DNNM is developed. Duncan mentions the
potential for improved graphical results and reporting utilizing FE,
but cautions against artificial accuracy being assumed when the input
parameters themselves are so variable. Wong~\cite{31} gives a useful
summary of potential sources of error in the FE modelling of slope
stability. Other important contributions in this area come from the
work of Sloan~\cite{32} who published the paper concerning the advances
in stability analysis that combine the limit theorems of classical
plasticity with finite elements to give rigorous upper and lower bounds
on the failure load.{\break} He proposed a new development, which
incorporates pore water pressures in finite-element limit analysis and
proposed of the new techniques for stability solutions including
foundations, anchors, slopes, excavations and tunnels.
Tschuchnigg~\cite{33} carried out a comparative study concerning the
strength reduction method and rigorous limit analyses which are based
on collapse theorems of plasticity.


\smallskip
\bigskip\noindent
\textit{Advantages of the finite element method}

\medskip\noindent
The advantages of a FE approach to slope stability analysis over
traditional limit equilibrium methods can be summarized as follows:
\begin{enumerate}[(a)]
\item No assumption needs to be made in advance about the shape or
location of the failure surface. Failure occurs ``naturally'' through
the zones within the soil mass in which the soil shear strength is
unable to sustain the applied shear stresses.

\item Since there is no concept of slices in the FE approach, there is
no need for assumptions about slice side forces. The FE method
preserves global equilibrium until ``failure'' is reached.

\item If realistic soil compressibility data are available, the FE
solutions will give information about deformations at working stress
levels.

\item The FE method is able to monitor progressive failure up to and
including overall shear failure.
\end{enumerate}

The developments which follow in this paper, are inspired by these
important contributions on this topic.

\section{Methodogogy} \label{sec4}

After this literature review related to the problematic, the
methodology adopted in this paper consists first of all in describing
the site location of the geotechnical structures design on the
geological and geomorphological aspects. Then the numerical modeling of
the structures will be approached by first validating the program for
calculating the bearing capacity and settlement of a strip footing
placed on an elastoplastic soil whose analytical solutions are known.
Then the numerical modeling of the problem will be compared with that
other software for the same problem. The different calculation phases
obeying the staged construction of these geotechnical structures will
be respected in our numerical modeling. Settlements will be analyzed at
the serviceability limit state and all other failure criteria at the
ultimate limit state.

\subsection{Presentation of the construction site of the geotechnical
structures} \label{sec4.1}

The project site is located in the Douala city in Cameroon. It is very
rugged with a steep incline. The objective is to build a series of
geotechnical structures to guarantee the stability of the road whose
projection from the level of the platform is at the top of slope.
Figure~\ref{fig1} illustrates the shadings of the contour lines of the
site, the position and the direction in cross-section of the
geotechnical structures under investigation.

\begin{figure}
\includegraphics{fig01}
{\vspace*{.4pc}}
\caption{\label{fig1}Shadings of the contour lines of the site, the
position and the direction in cross-section of the geotechnical
structures.}
\end{figure}

\subsection{Geology of site project} \label{sec4.2}

In the Littoral region of Cameroon, the sedimentary rocks form low
lying and gently undulating hills along the western side of the Dibamba
river. The strata have experienced tropical weathering, which has
resulted in the development of a lateritic residual soil profile. This
type of rock mass tropical leads to weathering profiles as presented by
Fookes~\cite{34}. The material viewed at site comprises Weathering
Grade~6, as seen in the trial pit photograph from site, with possibly
weathering Grade~5 material encountered towards the base of deeper
rotary boreholes~\cite{35,36}. Below this material, the stratum may
then transition more rapidly to weathered parent rock material. This
transition was not identified in the deepest exploratory boreholes that
went 30~m below ground level. In the exploratory holes observed there
was no clear evidence for overlying transported material, so material
may be weathered {\it in-situ}.

\subsection{Behaviour laws of subsoil, basement and backfill materials}
\label{sec4.3}

The first step of the proposed approach consists in characterizing on
the one hand the behaviour of the natural materials constituting the
structures whose stability is sought to be assessed by choosing an
elastoplastic law to represent them, and on the other hand in
estimating in a relevant way a value deterministic for each of the
parameters of the constitutive law, which associated with other
variables, will constitute the input data of the Finite Element model.
The choice of a constitutive law depends on the mechanism to be
modeled, the data available for the structures, but above all on the
precision sought and numerical considerations. We are interested in the
mechanisms that can lead structures under extreme stresses to failure.
However, to calculate the safety factor relating to the behavior of
geotechnical structures\vadjust{\pagebreak} by the Finite Element
Method, the evaluation of the elastic deformations does not matter
because it is the modeling of the appearance and the evolution of the
unstabilized plastic deformations which is important. Among the laws at
our disposal and described in Mestat~\cite{37}, we present here two
classic behavior laws of the elastoplastic type: those of Mohr--Coulomb
and Drucker--Prager. These two laws have the advantage of depending on
few parameters directly from triaxial and oedometric tests. They allow
a certain ``economy'' of data acquisition and management which are, in
particular those concerning the mechanical characterization of
materials, relatively few in number and not very varied on geotechnical
projects.

The Mohr--Coulomb criterion is the first that was proposed for
soils~\cite{37}. It is used for long-term for frictional and cohesive
soils. It is characterized by the relations:
{\begin{eqnarray}
&& \label{eq1} F(\sigma_{ij})= (\sigma_{1} - \sigma_{3}) - (\sigma_{1}
+ \sigma_{3}) \sin \varphi - 2 c \cos \varphi \leq 0\Seqnsplit
&& \label{eq2} G(\sigma_{ij}) = (\sigma_{1} - \sigma_{3}) - (\sigma_{1}
+ \sigma_{3}) \sin \psi + c^{\mathrm{ste}}
\end{eqnarray}}\noindent
where $c$ is the cohesion of the material, $\varphi$ is the friction
angle, and $\psi$ is the dilatancy angle. $\sigma_{1}$ and $\sigma_{3}$
represent the principal stresses ($\sigma_{1} \geq \sigma_{2} \geq
\sigma_{3}$). In the principal stress space, the surface defined by the
Mohr--Coulomb criterion is a hexagonal pyramid with axis ($\sigma_{1} =
\sigma_{2} = \sigma_{3}$) (Figure~\ref{fig2}b).

\begin{figure}
\includegraphics{fig02}
\caption{\label{fig2}Mohr--Coulomb and Drucker--Prager yield criteria.}
\end{figure}

Analogies are possible between the Mohr--Coulomb and Drucker--Prager
criteria (Figure~\ref{fig2}). Relations can be established between the
parameters ($\alpha, \gamma, k$) and (${\varphi}, {\psi}, c$) in
certain situations~\cite{37}. The Drucker--Prager criterion is written:
{\begin{eqnarray}
&& \label{eq3}f(\sigma) = \sqrt{J_{2}}+\alpha I_{1}-K\leq 0\Seqnsplit
&& \label{eq4}f(\sigma) = \sigma_{e}+\sin \alpha I_{1}-K'\leq 0.
\end{eqnarray}}\noindent

The correspondence between the Mohr--Coulomb criteria and those of
Drucker--Prager are as follows~\cite{37}:
{\begin{equation}\label{eq5}
\text{ALFA} = \text{ETA} = \frac{2\sqrt{3}\sin \varphi}
{9-\sin^{2}\varphi},\quad \text{GAMMA} = \frac{2\sqrt{3}\sin\psi}
{9-\sin^{2}\psi},\quad K = KL = \frac{6\sqrt{3}c\cdot\cos\varphi}
{9-\sin^{2}\varphi}.
\end{equation}}\noindent

Moreover, the choice of a criterion must also be guided by the possible
numerical difficulties of implementing the criterion. Some constitutive
laws, such as that of Mohr--Coulomb for example, have a failure
criterion comprising singular edges
(Figure~\ref{fig2}b),\vadjust{\pagebreak} which can numerically result
in convergence difficulties~\cite{28,38}. In view of the above, the
Drucker--Prager model implemented in the Cast3M Finite Elements code
will be used in this paper to describe the elastoplastic behavior of
untreated materials.

\subsection{Conceptual basis of the c-phi reduction procedure in Cast3M}
\label{sec4.4}

The c-phi reduction procedure developed is used to determine the factor
of safety in this part and part~\ref{sec5}. The concept of the strength
reduction method is based:

{Strength reduction method---Global safety factor concept: FOS}
{\begin{eqnarray}
&\displaystyle \label{eq6}
R_{n}=\|f-K(\varphi_{k,\mathrm{mob}},c_{k,\mathrm{mob}},
\psi_{k,\mathrm{mob}},\sigma_{t,\mathrm{mob}})\cdot u_{n}\|<\varepsilon
\Seqnsplit
&\displaystyle \label{eq7} \tan
\varphi_{k,\mathrm{mob}}=\frac{\tan\varphi_{k}} {\text{FOS}},\quad
c_{k,\mathrm{mob}}=\frac{c_{k}}{\text{FOS}},\quad \tan
\psi_{k,\mathrm{mob}}=\frac{\tan\psi_{k}}{\text{FOS}},\quad
\sigma_{t,\mathrm{mob}}=\frac{\sigma_{t,k}}{\text{FOS}}.
\end{eqnarray}}\noindent

{Strength reduction method---Partial safety factor concept: FCUT,}
(${\gamma}$)
{\begin{eqnarray}
&& R_{n}=\|f-K(\varphi_{d,\mathrm{mob}}, c_{d,\mathrm{mob}},
\psi_{d,\mathrm{mob}},\sigma_{t,d,\mathrm{mob}})\cdot u_{n}\|<
\varepsilon\nonumber\\
&& \tan\varphi_{d,\mathrm{mob}}=\text{FCUT}\cdot\frac{\tan\varphi_{k}}
{\gamma_{\varphi}},\quad c_{d,\mathrm{mob}}=\text{FCUT}\cdot\frac{c_{k}}
{\gamma_{c}},\nonumber\\
&& \tan\psi_{d,\mathrm{mob}}=\text{FCUT}\cdot\frac{\tan
\psi_{k}}{\gamma_{\psi}},\quad \sigma_{t,d,\mathrm{mob}}=\text{FCUT}
\cdot\frac{\sigma_{t,k}}{\gamma_{\sigma,t}} \label{eq8}
\end{eqnarray}}\noindent
where $R$ represents the strength reduction method; $f$ represents the
total force mobilized during shear failure of model stability; $K$ is
the matrix of design strength parameters as a function of friction
angle, dilatance angle and cohesion; $U_n$ is the displacement obtained
during strength failure as a function of design strength parameters and
the total force mobilized. $\varepsilon$ represents the relative error
of the method.

For the analysis of strength reduction method in ULS on Cast3M code, we
have used a variable FCUT (FaCtor UTilisation): this factor is the
unknown in the stability computation step. In order to use the partial
safety concept with Cast3M, one has to reduce the strength parameters
with a prescribed partial safety factor to their design value. The
model must be stable with this design strength parameters of friction
angle and cohesion. The design values of the friction angle and
cohesion are multiplied with FCUT and the model is recalculated. As
long as there is convergence, the value of FCUT can be reduced. For
example, from 1.0 in steps of 0.01 and latter is the accuracy of the
solution. If FCUT gets very low, a certain part of the model plastifies
and the convergence gets lost, leading to the solution for the factor
of utilization. We have using this brute force strategy but if there
was an response from PASAPAS or UNPAS in some of its variables on the
convergence quality, a more efficient and intelligent servo mechanism
could be used for a bracketed search with a lower and upper limit value
focusing in to the final solution in few steps. In our brute force
algorithm, we have to count down and if the solution of FCUT is 0.32
for example, we need 68 computational runs starting from 1.0 with step
size 0.01 with successful convergence, before we arrive at the
solution. The failure and shadings are strongly influenced by the
value{\break} of FCUT. The closer we get to FCUT which makes the
solution non-convergent, the clearer the failure curve appears.

\subsection{Validation of the elastoplastic calculation program on
Cast3M: calculation of bearing capacity and settlement of a strip
footing} \label{sec4.5}

Before modeling the problem investigated in this paper, we propose to
first validate the code written on Cast3M for the elastoplastic
calculation of the\vadjust{\pagebreak} stability of geotechnical
structures.{\break} This is the calculation of the bearing capacity of
a strip footing whose solution is known. Recall that shortly after
isotropic or anisotropic linear elasticity had been introduced into a
computer code by the Finite Element Method around 1965, computer
scientists proposed algorithms to deal with nonlinear behaviour laws
(nonlinear elasticity of Duncan or the perfectly plastic
laws)~\cite{14,39,40}. These laws are developed by rheologists from
tests carried out in the laboratory (mainly triaxial and oedometric
tests), and they make it possible to better represent the behaviour of
soils (appearance of irreversible deformations and rupture).

The plasticity (or failure) criteria are the basis for calculating the
limit loads that can be borne by soil masses~\cite{35,41,42,43,44,45}.
The Finite Element Method (the unknowns are the displacements and the
stresses) makes it possible, by carrying out an incremental
calculation, to highlight the value of the limit
loads~\cite{46,47,48,49,50,51,52}. We cite here some contributions to
the Finite Element study of foundations on elastoplastic materials in
plane strain~\cite{16,17,18,46,47,48,49,52,53,54,55,56,57,58,59,60,61}.

Figure~\ref{fig3} shows the view of Cast3M 2-D mesh of a strip footing
with $B = 2$~m on a 0.25~m thick concrete on a frictional soil. The
thickness of the soil layer is taken 4~m and the material behaviour is
represented by Drucker--Prager model for Cast3M~\cite{62} calculation
and Mohr--Coulomb model for Plaxis calculation (for comparison). The
characteristics of material are presented in Table~\ref{tab1}. The
phreatic level is located in the base of the strip footing.

\begin{figure}
\includegraphics{fig03}
\caption{\label{fig3}Presentation of the geometry and the view of 2-D
mesh of a strip footing.}
{\vspace*{.8pc}}
\end{figure}

\begin{table}[t!] %tab1
\caption{\label{tab1}The characteristics of materials for calculation
of the footing}
\begin{tabular}{ccccc}
\thead
Parameter & Symbol & Sandy clay & Thick concrete & Unit\\
\endthead
Material model & Model & Mohr--Coulomb & Linear elastic & -\\
Type of behaviour & Type & Drained & Non-porous & -\\
Weight above phreatic level & $\gamma_{\mathrm{unsat}}$ & 16 & 24 & kN/$\text{m}^{3}$\\
Weight below phreatic level & $\gamma_{\mathrm{sat}}$ & 18 & - & kN/$\text{m}^{3}$\\
Young's modulus & $E_{\mathrm{ref}}$ & $5 \times 10^{3}$ & $2 \times 10^{7}$ & kN/$\text{m}^{2}$\\
Poisson's ratio & ${\nu}$ & 0.35 & 0.15 & -\\
Cohesion & $c$ & 5 & - & kN/$\text{m}^{2}$\\
Friction angle & ${\varphi}$ & 20 & - & {\textdegree}\\
Dilatancy angle & ${\psi}$ & 0 & - & {\textdegree}
\botline
\end{tabular}
\end{table}

The Cast3M elasto-plastic calculation is performed in 2-D plane
strains. Two calculation phases were carried out, the first phase
concerns the simple elastoplastic calculation which makes it possible
to obtain the elastoplastic settlement according to the applied load as
long as convergence is ensured. The second phase concerns the
calculation, the presentation of the shear failure mechanism, the
factor of safety as well as the bearing capacity of the strip footing
by strength parameters reduction c-phi method until failure
(non-converging calculation).

Figures~\ref{fig4}--\ref{fig9} present the displacements and vertical
stresses (in phase~1) of the strip footing submitted at vertical
pressure of 100~kPa from elastoplastic analysis. Figure~\ref{fig4}
clearly shows the asymmetry of the horizontal displacements in the
ground. The vertical displacements isobars are higher under the
foundation and gradually decrease until they cancel out as one move
away from the point of application of the vertical load
(Figures~\ref{fig5} and~\ref{fig6}). Figure~\ref{fig7} shows the load
displacement curve for the footing at top of thick concrete. The
deformed mesh of the model{\break} through \mbox{Figure}~\ref{fig8},
represents the incremental displacements. The footing sinks following
the application of the load causing a heaving of the ground around the
footing. The vectors of incremental displacements indeed provide the
form of the kinematics failure; this is in agreement with the
theoretical results (rigid corner under the foundation
(Figure~\ref{fig11})). The total displacement has calculated by the
formula: $|U_{\mathrm{total}}| =\sqrt{U_{x}^{2}+U_{y}^{2}}$.

\begin{figure}
\includegraphics{fig04}
{\vspace*{-.1pc}}
\caption{\label{fig4}Horizontal displacements (m) at 5th incremental
time in SLS of the strip footing submitted at vertical pressure of
100~kPa.}
{\vspace*{-.1pc}}
\end{figure}

\begin{figure}
\includegraphics{fig05}
{\vspace*{-.1pc}}
\caption{\label{fig5}Vertical displacements (m) at 5th incremental time
in SLS of the strip footing submitted at vertical pressure of 100~kPa.}
{\vspace*{-.1pc}}
\end{figure}

\begin{figure}
\includegraphics{fig06}
\caption{\label{fig6}Shadings of total displacements (m) in SLS of the
strip footing submitted at vertical pressure of 100~kPa: (a)~at 5th
incremental time; (b)~at full incremental time.}
\end{figure}

\begin{figure}
{\vspace*{-.25pc}}
\includegraphics{fig07}
\caption{\label{fig7}Load displacement curve for the footing.}
{\vspace*{-.4pc}}
\end{figure}

\begin{figure}
{\vspace*{-.25pc}}
\includegraphics{fig08}
\caption{\label{fig8}Deformed incremental total displacements (m) at
full time in SLS of the strip footing submitted at vertical pressure of
100~kPa (displacements scaled by factor of~2).}
\end{figure}

\begin{figure}
\includegraphics{fig09}
{\vspace*{.15pc}}
\caption{\label{fig9}Vertical stresses (Pa) in ULS of the strip footing
submitted at vertical pressure of 100~kPa.}
{\vspace*{.4pc}}
\end{figure}

In the calculation program, each load is applied in 20 increments. For
100~kPa, applied, the pressure corresponding at the 5th increment is
25~kPa (Figures~\ref{fig5}--\ref{fig7}).

The input value of pressure is 100 kPa and the Factor of Safety reaches
approximately 1.0101 (for a calculation carried out with very fine mesh
of 1192 elements of the type TRI3 or TRI6). Therefore the failure
stress is equal to 100~kPa $\times$ 1.0101 $=$ 101.01~kPa. This bearing
capacity value should be compared with the reference analytical
solutions of Vesic~\cite{63}, Brinch-Hansen~\cite{64,65},
Meyerhof~\cite{66,67} and the solution obtained by the Plaxis
software~\cite{68}. Given the formula for analytical bearing capacity
of a strip footing (Figure~\ref{fig11}):

Figure~\ref{fig10} presents the calculation result of phase~2. The
shear failure mechanism is clearly represented through the shadings of
displacements under the foundation. The calculation stops automatically
as soon as the solution is no longer convergent (plastification of the
model).
{\begin{eqnarray}
&& \label{eq9}\frac{Q_{f}}{B}=cN_{c}+\frac{1}{2}\gamma B\cdot
N_{\gamma}\Seqnsplit
&& N_{q}=\mathrm{e}^{\uppi\tan\varphi}\cdot\tan^{2}\left(45+\tfrac{1}
{2}\varphi\right)\nonumber\\
&& \label{eq10} N_{c}=(N_{q}-1)\cot \varphi\Seqnsplit
&& \label{eq11} N_{\gamma}=\begin{cases}
2(N_{q}+1)\tan \varphi => (\text{Vesic})\\
1.5(N_{q}-1)\tan \varphi => (\text{Brinch-Hansen})\\
(N_{q}-1)\tan (1.4\varphi) => (\text{Meyerhof}).
\end{cases} 
\end{eqnarray}}\noindent

\begin{figure}
\includegraphics{fig10}
{\vspace*{.2pc}}
\caption{\label{fig10}Shadings of shear displacements in ULS (m) and
FOS (FOS $=$ 1.0101) of the strip footing submitted at vertical
pressure of 100~kPa.}
\end{figure}

\begin{figure}
\includegraphics{fig11}
\caption{\label{fig11}Shear failure mechanism of a shallow foundation
for theoretical results.}
{\vspace*{.5pc}}
\end{figure}

Filling in given soil data:
{\begin{eqnarray*}
&& N_{q}=\mathrm{e}^{\uppi\tan 20}\cdot\tan^{2}(55)=6.4\\
&& N_{c}=(6.4-1)\cot(20)=14.83\\
&& N_{\gamma}=\begin{cases}
2(6.4+1)\tan (20)=5.39 => (\text{Vesic})\\
1.5(6.4-1)\tan(20)=2.95 => (\text{Brinch-Hansen})\\
(6.4-1)\tan (28)=2.87 => (\text{Meyerhof}).
\end{cases} 
\end{eqnarray*}}\noindent

The effective weight of the soil: $\gamma' = \gamma_{\mathrm{sat}} -
\gamma_{w} = 18-10 = 8~\text{kN}/\text{m}^{2}$
{\[
\frac{Q_{f}}{B}=c\cdot N_{c}+0.5\gamma'B\cdot N_{\gamma}=\begin{cases}
5*14.83+0.5*8*2*5.39\approx 117~\text{kN/m}^{2} => (\text{Vesic})\\
5*14.83+0.5*8*2*2.95\approx 98~\text{kN/m}^{2} => (\text{Brinch-Hansen})\\
5*14.83+0.5*8*2*2.87\approx 97~\text{kN/m}^{2} => (\text{Meyerhof}).
\end{cases} 
\]}\noindent

For the results in PLAXIS modeling, in addition to the mesh used in
these calculations were performed using a very coarse mesh with a local
refinement at the bottom of the footing and a very fine mesh. Fine
meshes will normally give more accurate results than coarse meshes.
Instead of refining the whole mesh, it is generally better to refine
the most important parts of the mesh, in order to reduce computing
time. Here we see that the differences are small (when considering
15-noded elements), which means that we are close to the exact
solution. The accuracy of the \mbox{15-noded} element is superior to
the 6-noded element, especially for the calculation of failure loads.
The results of fine/coarse and 6-noded/15-noded analyses are given
below. Table~\ref{tab2} presents the results for the maximum load
reached on a strip footing on the drained sub-soil for different 2D and
3D meshes.

\begin{table}[t!] %tab2
\fontsize{9}{11}\selectfont
\caption{\label{tab2}Results for the maximum load reached on a strip
footing for 2D (for Cast3M) and different 2D and 3D meshes (for PLAXIS)
and analytical solutions}
\tabcolsep 3pt
\begin{tabular}{cccccc}
\thead
Solver & Mesh size & \parbox[t]{1cm}{\centering Element type} &
\parbox[t]{1.5cm}{\centering Number of elements} &
\parbox[t]{1.65cm}{\centering Maximum load (kN/m)} &
\parbox[t]{1.6cm}{\centering Failure load
(kN/$\text{m}^{2}$)}\vspace*{2pt}\\
\endthead
\morerows{5}{PLAXIS Code} & \parbox[t]{4.5cm}{\raggedright\hangindent 1pc Very coarse
mesh with local refinements under footing}\vspace*{2.5pt} & 6-noded &
\0\079 & 281 & 146\\
& Coarse mesh & 6-noded & \0121 & 270 & 141\\
& Very fine mesh & 6-noded & 1090 & 229 & 121\\
& \parbox[t]{4.5cm}{\raggedright\hangindent 1pc Very coarse mesh with local refinements
under the footing}\vspace*{2.5pt} & 15-noded & \0\079 & 236 & 124\\
& Coarse mesh & 15-noded & \0121 & 248 & 130\\
& Very fine mesh & 15-noded & 1090 & 220 & \textbf{116}\vspace*{10pt}\\
Cast3M code & Very fine mesh & 3-noded & 1192 & - &
\textbf{101}\vspace*{10pt}\\
\morerows{2}{\parbox[t]{1.55cm}{\centering Analytical results}} & Vesic
&&&& 117\\
& Brinch-Hansen &&&& \0\textbf{98}\\
& Meyerhof &&&& \0\textbf{97}\vspace*{2pt}
\botline
\end{tabular}
\end{table}

In this Table~\ref{tab2}, the failure load for PLAXIS calculation has
been calculated as:
{\begin{equation}\label{eq12}
\frac{Q_f}{B}=\frac{\text{Maximumforce}}{B}+\gamma_{\mathrm{concrete}}
* d = \frac{\text{Maximumforce}}{2}+24*0.25.
\end{equation}}\noindent

From the above results of PLAXIS, it is clear that fine FE meshes give
more accurate results. On the other hand the performance of the
15-noded elements is superior over the performance of the lower order
6-noded elements. Needless to say that computation times are also
influenced by the number and type of elements. The Cast3M results are
generally in good agreement with the analytical results of the collapse
load. The bearing capacity obtained by our elastoplastic modeling on
Cast3M is 101~kPa very close to the Brinch-Hansen and Meyerhof
reference solutions, validating our calculation program very
satisfactorily. For very fine mesh of the numerical model of the
footing, the difference between Plaxis and Cast3M results is 7\%
(Table~\ref{tab2}). This gap comes the difference between
Drucker--Prager and Mohr--Coulomb materials. Drucker--Prager is the
advanced ideal elasto-plasticity model by the number of parameters to
determined. He is adequate to describe the sophisticated and realistic
behaviour of geotechnical structure submitted by any load by correcting
the shortcomings of Mohr--Coulomb model which only includes five
parameters ($E_{\mathrm{ref}}$, $\nu$, $c$, $\varphi$, $\psi$)
determined directly from conventional triaxial tests.

The analysis on the mesh density and element type (TRI3 and TRI6) in
Cast3M (Table~\ref{tab3}) show that, the difference between element
type TRI3 and TRI6 results varies between 30 (for the very coarse mesh
with local refinements under footing) to 0\% (for the very fine mesh).
In the developed procedure in this paper, when the mesh is very fine,
the adjustment of the factor FCUT (this factor is the unknown in the
stability computation step. See part~\ref{sec4.4}) makes it possible to
regulate numerical gap between the very coarse mesh effect to very fine
mesh of the model. Table~\ref{tab3} clearly shows the interest of this
factor. For the very fine mesh with TRI3 element, it is reduced
progressively by step of 0.01 from 1 to 0.65 (35~steps) to obtain the
same result with a TRI6 element whose FCUT is reduced progressively
from 1 to 0.91 (9~steps).

\begin{table}[t!] %tab3
\caption{\label{tab3}Results for the maximum load reached on a strip
footing for different 2D meshes in Cast3M (partial safety factors:
$\gamma_{\sigma t}$) and analytical solutions}
\fontsize{8.5}{10.5}\selectfont
\tabcolsep3pt
\begin{tabular}{ccccccccc}
\thead
Solver & Mesh size & \parbox[t]{.85cm}{\centering Element type} &
\parbox[t]{1.2cm}{\centering Number of elements} &
\parbox[t]{.9cm}{\centering Failure load (kN/$\text{m}^{2}$)} & FCUT
& $\gamma_{\sigma t}$ & FOS & \parbox[t]{1.6cm}{\centering
($((F_{3\mbox{-}\mathrm{noded}}) -
(F_{6\mbox{-}\mathrm{noded}}))$/{\unskip\break}$(F_{6\mbox{-}\mathrm{noded}})$)
$\times$ 100}\vspace*{2pt}\\
\endthead
\morerows{9}{\parbox[t]{1.05cm}{\centering Cast3M code}} &
\parbox[t]{4.3cm}{\raggedright\hangindent 1pc Very coarse mesh with local refinements
under footing}\vspace*{1.5pt} & 3-noded & \0\076 & 166.66 & 0.65 & 1.08
& 1.67 & 29.93\\
& Coarse mesh & 3-noded & \0120 & 138.89 & 0.65 & 0.90 & 1.39 & 14.59\\
& Coarse mesh & 3-noded & \0242 & 123.46 & 0.65 & 0.80 & 1.23 & 13.89\\
& Fine mesh & 3-noded & \0370 & 111.11 & 0.65 & 0.72 & 1.11 & \06.66\\
& Very fine mesh & 3-noded & 1192 & \textbf{101.11} & 0.65 &
\textbf{0.66} & \textbf{1.01} & \0\textbf{0.00}\\
& \parbox[t]{4.3cm}{\raggedright\hangindent 1pc Very coarse mesh with local refinements
under the footing}\vspace*{1.5pt} & 6-noded & \0\076 & 128.28 & 0.91 &
1.17 & 1.28 & -\\
& Coarse mesh & 6-noded & \0120 & 121.21 & 0.91 & 1.10 & 1.21 &\\
& Coarse mesh & 6-noded & \0242 & 108.40 & 0.91 & 0.99 & 1.08 &\\
& Fine mesh & 6-noded & \0370 & 104.17 & 0.91 & 0.95 & 1.04 &\\
& Very fine mesh & 6-noded & 1192 & \textbf{101.10} & 0.91 &
\textbf{0.92} & \textbf{1.01} &{\vspace*{7pt}}\\
\morerows{2}{\parbox[t]{1.05cm}{\centering Analytical results}} & Vesic &&&
\textbf{117} &&&& -\\
& Brinch-Hansen &&& \0\textbf{98} &&&&\\
& Meyerhof &&& \0\textbf{97} &&&&
\botline
\end{tabular}
\end{table}

The code programming algorithm is made so that at the starting of the
calculation a value of FCUT $=$ 1 is introduced. As the calculation
evolves, this value is gradually reduced from 1 by steps of 0.01 until
a value which allows to have the solution (end of convergence). The
global safety factor is therefore obtained automatically at the end of
this calculation. It is not necessary to introduce the partial safety
factors. After calculation the global safety factor; the partial safety
factors are obtained automatically from FCUT adjusted during the
calculation and the FOS obtained according to Equation~(\ref{eq8}). The
FOS is not equal to 1/FCUT, because the partial safety coefficients are
different from~1. It is also possible to introduce any value $n$ of
FCUT less than 1 at the start of the calculation. If this calculation
converges, the calculation time which was necessary to reduce FCUT
progressively up to $n$ is gained and the calculation continues until
the final solution. The values of the partial safety coefficients have
been inserted in Table~\ref{tab3}.

\section{Numerical analysis of retaining walls--slope stability and
backfill} \label{sec5}

The materials used for the backfills come from two different quarries
located a few kilometers from the project area. Backfill~1 consits of
greyish pozzolan with a strong granular predominance with a low fines
content. Backfill~2 is constituted by the Reddish clay gravel at
gravelly sand with a spread out grain size.

The numerical modeling of the retaining structures is made in 2D, plane
deformations from the Finite Elements Cast3M calculation code with
triangular elements with 3 nodes (TRI3). \mbox{Figure}~\ref{fig12}
shows the 2D mesh view of the model and the borehole log with the
altimetry of the \mbox{different} backfills and the name of the
sedimentary rocks involved. The boundary conditions of the model are
standard, i.e. blocking of all displacements at the bottom of the model
\mbox{($Ux = Uy = 0$)} and blocking of horizontal displacements on the
left and right bound of the model ($Ux = 0$).

\begin{figure}
\includegraphics{fig12}
\caption{\label{fig12}Detailed view of the 2D mesh and soil clusters of
retaining walls and backfills.}
\end{figure}

In this paper, our analyze is performed at the Serviceability Limit
State (SLS) for deformation-based approach and the Ultimate Limit State
analysis (ULS) for stress and strength-based approach. The forces
generated by the self-weight of the soil are computed using a standard
gravity ``turnon'' procedure involving integrals over each element of
the form:
{\begin{equation} \label{eq13}
p^{(e)}=\gamma \int_{{V^{e}}}N^{T}\,\mathrm{d}(\text{vol})
\end{equation}}\noindent
where $N$ values are the shape functions of the element and the
superscript $e$ refers to the element number. This integral evaluates
the area of each element, multiplies by the total unit weight of the
soil and distributes the net vertical force consistently to all the
nodes. These element forces are assembled into a global gravity force
vector that is applied to the FE mesh in order to generate the initial
stress state of the problem. 

\subsection{Materials model} \label{sec5.1}

Drucker--Prager model~\cite{69}, ideal elasto-plasticity without
hardening or softening, they constant stiffness parameters are adjusted
of strength parameters to the average of Mohr--Coulomb compression and
tension. Apart from unit weight gravity~$\gamma$, Young Modulus~$E$ and
Poisson's \mbox{ratio},{\break} this~model requires 8 additional
parameters, some of which are deduced from the Mohr--Coulomb model for
its use in the $F\cdot E$ Cast3M code. Its other parameters are as
follows:
{\begin{eqnarray}
&& \label{eq14} \text{ALFA} = \text{ETA} =\frac{2\sqrt{3}\sin\varphi}
{9-\sin^{2}\varphi}\Seqnsplit
&& \label{eq15} \text{BETA} = \text{MU} = \text{DELTA} = \sqrt{\frac{2}
{3}}\Seqnsplit
&& \label{eq16} K = KL = \frac{6\sqrt{3}c\cdot\cos \varphi}{9-\sin^{2}
\varphi}\Seqnsplit
&& \label{eq17} \text{GAMMA} = \frac{2\sqrt{3}\sin\psi}{9-\sin^2\psi}.
\end{eqnarray}}\noindent

In our computations on Cast3M, we have used a variable FCYS (FaCtor
Yield Surface) for the built in Drucker--Prager plasticity model. In
order to adjust the Mohr--Coulomb strength parameters to the circular
cone of the Drucker--Prager failure surface (see Figure~\ref{fig2}),
different strategies are known. The factor FCYS regulates the
adjustment~\cite{70}. The value of 0.0 is tension, 1.0 is compression,
0.5 is their average. Other values would mean a weighted result and
values above 1.0 use the surface equality approach. Values below 0.0
default to 1.0 and the compression adjustment. In this paper FCYS is
fixed at 0.5. Table~\ref{tab4} presents the materials properties of the
model. $ \text{BETA} = \text{MU} = \text{DELTA} = \sqrt{{2}/{3}}$ is
implementation specific for Drucker--Prager model in Cast3M.

\begin{table}[t!] %tab4
\fontsize{7.8}{9.8}\selectfont
\tabcolsep2.75pt
\caption{\label{tab4}Materials properties of the model}
\begin{tabular}{cccccccc}
\thead
Parameter & Symbol & \parbox[t]{1.7cm}{\centering Subsoil: thick Sandy
clay} & Basement & \parbox[t]{1.1cm}{\centering Concrete retaining wall
1 \& 2}\vspace*{2pt} & Backfill 1 & Backfill 2 & Unit\\
\endthead
Material model & Model & Mohr--Coulomb & Mohr--Coulomb & Linear elastic &
Mohr--Coulomb & Mohr--Coulomb & -\vspace*{2pt}\\
Type of behaviour & Type & Drained & Drained & Non-porous & Drained &
Drained & -\vspace*{2pt}\\
\parbox[t]{1.8cm}{\centering Weight above phreatic level}\vspace*{2pt}
& $\gamma_{\mathrm{unsat}}$ & 20 & 20 & 25 & 20 & 20 &
kN/$\text{m}^{3}$\\
Young's modulus & $E_{\mathrm{ref}}$ & $40 \times 10^{3}$ & $40 \times
10^{3}$ & $4 \times 10^{7}$ & $40 \times 10^{3}$ & $40 \times 10^{3}$ &
kN/$\text{m}^{2}$\vspace*{2pt}\\
Poisson's ratio & ${\nu}$ & 0.3 & 0.3 & 0.2 & 0.3 & 0.3 &
-\vspace*{2pt}\\
Cohesion & $c'$ & 10 & 0.1 & - & 0.1 & 1 &
kN/$\text{m}^{2}$\vspace*{2pt}\\
Friction angle & $\varphi'$ & 20 & 35 & - & 35 & 35 &
{\textdegree}\vspace*{2pt}\\
Dilatancy angle & ${\psi}$ & 0 & 5 & - & 0 & 0 & {\textdegree}\
\botline
\end{tabular}
\end{table}

\subsection{Analysis of problem phase by phase: staged construction}
\label{sec5.2}

To know the contribution of each phase to the global stability, the
displacements are calculated at the Serviceability Limit State (SLS)
and the Ultimate Limit State analysis (ULS) for stress and strength.
The solution strategy is based by treatment of the problem in 2-D plane
strain (geometry and mesh), the definition of soils and retaining walls
and the definition of all external loads (self-weight and charges) and
the implementation of the construction phases by stepwise imposition of
loads. The following staged constructions are analyzed in Cast3M: the
first phase is based on the application of geostatic stress state in
the subsoil and basement (stress adjustment or gravity loading). In
phase~2, the retaining wall construction~1 is realized. The layered
backfill~1 is realized in phase~3; the retaining wall construction~2 in
phase~4; the 2nd layered backfill in phase~5; the application of
surface load for road construction following the GTR
recommendation~\cite{71} (regarded as variable) in phase~6 and the
calculation of general stability of the construction in phase~7. 

The ``turn-on'' procedure is described clearly according to the
following calculation steps{\break} \cite{72,73}: after generation of
the mesh of the global model, the calculation of the elastic phase is
carried out in order to obtain the stresses and the elastic
displacements in the structure. The elastoplastic calculation is
carried out PASAPAS in the code cast3M. Displacements are initialized
to zero at the beginning of each phase. Each construction stage of the
structure corresponds to a PASAPAS calculation phase in the cast3M
code. In each phase, the parameters of the entire model are inserted
and the calculation results are logged in a table. The phase number
corresponds to the number of calculation tables. The results of each
phase are recorded in the tables. From these tables, the next phase
begins with the previous phase.

\subsection{Phase 0: elastic modeling in one phase} \label{sec5.3}

The purpose of the elastic calculation carried out is to obtain the
stresses and displacements throughout the model. The entire
geotechnical structure is built in a single phase. The gravity load are
applied on the model. The objective of this part is to show the limits
and insufficiencies of a modeling of a geotechnical structure by using
a model with elastic behavior. The results of this part are to be
compared with those of part~\ref{sec5.4} whose behavior model
realistically describes the behavior of geotechnical structures.
Figures~\ref{fig13}--\ref{fig16} present the displacements (in phase~0)
of the model from elastic analysis.

\begin{figure}
\includegraphics{fig13}
\caption{\label{fig13}Horizontal displacements: $Ux$ (m) in SLS of the
model.}
\end{figure}

\begin{figure}
\includegraphics{fig14}
\caption{\label{fig14}Vertical displacements: $Uy$ (m) in SLS of the
model.}
{\vspace*{.5pc}}
\end{figure}

\begin{figure}
\includegraphics{fig15}
{\vspace*{-.6pc}}
\caption{\label{fig15}Total displacements (m) in SLS of the model.}
\end{figure}

\begin{figure}
\includegraphics{fig16}
{\vspace*{-.15pc}}
\caption{\label{fig16}Deformed of the total displacements (m) in SLS of
the model (displacements scaled by factor of 20).}
\end{figure}

The interest of the computation of the elastic phase is to obtain
displacements and stresses in the model. The total displacement in the
elastic calculation phase is 21~cm.

\subsection{Drucker--Prager elasto-plastic modeling in staged
construction} \label{sec5.4}

In this part, the stresses and displacements are calculated for each
construction stage up to the Serviceability phase.

\subsubsection{Phase~1: elasto-plastic modeling} \label{sec5.4.1}

This phase is based on the application of geostatic stress state in the
subsoil and basement (gravity loading). The installation of the subsoil
and basement causes the appearance of the vertical, horizontal and
total displacements of 3.22~cm due to gravity loading.
Figures~\ref{fig17} and~\ref{fig18} present the total displacements and
vertical stress (in phase~1) of the model from elastoplastic analysis.

\begin{figure}
\includegraphics{fig17}
{\vspace*{-.45pc}}
\caption{\label{fig17}Total displacements (m) in SLS of the subsoil and
basement in phase~1.}
\end{figure}

\begin{figure}
\includegraphics{fig18}
{\vspace*{-.45pc}}
\caption{\label{fig18}Vertical stress (Pa) of the model in ULS
(phase~1).}
\end{figure}

\subsubsection{Phase~2: elasto-plastic modeling} \label{sec5.4.2}

In this phase~2, the retaining wall 1 is built. This phase begging
after phase~1, the displacements are initialized to zero in the model
and is based on the application of geostatic stress of the retaining
wall~1 on drained sub-soil (gravity loading). Although displacements on
the basement and subsoil are controlled with the realization of this
one with extra embankment{\break} (over-elevation) whose thickness
corresponds to the vertical displacement, the construction of the
retaining wall~1, causes an appearance of settlements of 8.83~cm.
Figures~\ref{fig19}--\ref{fig21} present the total displacements,
vertical plastic strains and vertical stress (in phase~2) of the model
from elastoplastic analysis.

\begin{figure}
\includegraphics{fig19}
\caption{\label{fig19}Deformed of the total incremental displacements
(m) in SLS of the retaining wall~1 on drained sub-soil (displacements
scaled by factor of 20) in phase~2.}
\end{figure}

\begin{figure}
\includegraphics{fig20}
{\vspace*{-.15pc}}
\caption{\label{fig20}Vertical plastic strains in SLS of the retaining
wall~1 on drained sub-soil in phase~2.}
{\vspace*{.3pc}}
\end{figure}

\begin{figure}
\includegraphics{fig21}
{\vspace*{-.3pc}}
\caption{\label{fig21}Vertical stress (Pa) of the model in ULS
(phase~2).}
{\vspace*{.5pc}}
\end{figure}

\subsubsection{Phase 3: elasto-plastic modeling and the strength
reduction method} \label{sec5.4.3}
\paragraph{Phase 3.1: elasto-plastic modeling.} \label{sec5.4.3.1}

In this phase~3.1, the layered backfill~1 is realized. This phase
starts after phase~2, the displacements are initialized to zero in the
model and is based on the application of geostatic stress of the
layered backfill~1 on drained sub-soil (gravity loading).
Figures~\ref{fig22}--\ref{fig24} present the total displacements,
plastic strains and vertical stress (in phase~3.1) of the model from
elastoplastic analysis. Although the displacements are reset to zero on
the subsoil, basement and on the retaining wall~1, the construction of
layered backfill~1 generates a horizontal displacement of the wall and
the total displacements of 15~cm due to its own weight.

\begin{figure}
\includegraphics{fig22}
{\vspace*{-.8pc}}
\caption{\label{fig22}Total displacements (m) in SLS of the model in
phase~3.1.}
\end{figure}

\begin{figure}
\includegraphics{fig23}
{\vspace*{-.25pc}}
\caption{\label{fig23}Vertical plastic strains in SLS in phase~3.1.}
{\vspace*{.5pc}}
\end{figure}

\begin{figure}
\includegraphics{fig24}
{\vspace*{-.4pc}}
\caption{\label{fig24}Vertical stress (Pa) of the model in ULS
(phase~3.1).}
\end{figure}

\paragraph{Phase 3.2: strength reduction method.} \label{sec5.4.3.2}

This phase starts after phase~3.1, the displacements are initialized to
zero in the model. The analysis is based on the limit equilibrium of
external and internal forces in the model (gradual reduction of the
mobilized strength). In the Cast3M calculation code, the advanced
solution is performed on recursive reduction of the mobilized strength
using a search strategy based on a convergence criterion. The
part~\ref{sec4.4} describe the concept of the strength reduction method
developed in this paper. Phi-c reduction phase cannot be used as a
starting condition for another calculation phase because it ends in a
state of failure. Therefore it is advised to define all safety analyses
at the end of the list of calculation phase and to use the start from
phase parameter as a reference to the calculation phase for which a
safety factor is calculated.

The Factor of Safety calculation procedure is based on the calculation
of a stable initial stress and corresponding deformation state
recursive reduction of the shear and tensile strength until a
convergence to a prescribed tolerance is still obtainable. The
restrictions on suitable material behaviour of Drucker--Prager, which is
the ideal elasto-plasticity are considered. This model is adequate to
describe the sophisticated and realistic behaviour of geotechnical
structure submitted by any load. Figure~\ref{fig25} presents the total
shear displacements (in phase~3.2) of the model from c-phi reduction
strength method (ULS, GEO-3). In this phase the Factor of Safety is FOS
$=$ 1.5625. The duration of the calculation in this phase is 40~min.

\begin{figure}
\includegraphics{fig25}
\caption{\label{fig25}(a)~Shadings of shear total incremental
displacements (m); (b)~shadings of deformed of the total shear
incremental displacements (m) in ULS of the model (displacements scaled
by factor of~0.25) in phase~3.2.}
{\vspace*{-.2pc}}
\end{figure}

\looseness=-1
The benefit of this intermediate calculation phase is to verify the
stability of the structure after the installation of the thick backfill
before continuing with the other construction phases. According to the
result of this phase, the sliding FOS is equal to 1.56 guaranteeing the
stability of the structure despite the shear displacements of 71~cm
which are canceled by the extra embankment and the abutment of the
retaining wall~1. In order to validate our elastoplastic calculation
program for the c-phi parameters reduction on Cast3M, the results of
this part were compared with those of other geotechnical software
commonly used in industrial environments. For these modeling with these
other softwares, the behaviour of materials obeys the criterion of
rupture of Mohr--Coulomb. Several software or methods were used for the
comparison of the results.

These include: ULS-Geo-3, failure kinematics with the code Flac2D
(FDM); ULS-Geo-3, failure kinematics with the code Tochnog (FEM); LEM,
Limit Equilibrium Method; KEM, Kinematic Element Method; Method of
slices with circular slip surface (Bishop method); Method of slices
with polygonal slip surface (Janbu method) and Strength Reduction
Finite Element Limit Analysis (SR-FELA)~\cite{33} implemented in OPTUM
G2~\cite{74}. The results in Figure~\ref{fig26} show the shear failure
mechanisms of the structure. The comparative FOS values obtained for
each method~\cite{9,62,68,74,75,76,77} are summarized in
Table~\ref{tab5}. The FOS obtained in our calculation is~1.56. It is
the same value that is obtained for the other Finite Element
calculation codes presented (TOCHNOG and PLAXIS) in Table~\ref{tab5}.
This calculation phase~3.2 is very useful, because it makes it possible
to know whether the structure is stable before to continue
construction. With this FOS $=$ 1.56, the safety of this part of the
structure is guaranteed.

\begin{figure}
\includegraphics{fig26}
\caption{\label{fig26}(a)~Total displacements in Ultimate Limit State
with the code Flac3D (FDM); (b)~total displacements in ULS with the
code Tochnog (FEM); (c)~ULS analysis with the Limit Equilibrium Method
(LEM); (d)~ULS analysis with the Kinematic Element Method (KEM);
(e)~ULS analysis with the method of slices (Bishop method); (f)~ULS
analysis with the method of slices (e.g, Janbu method).}
\end{figure}

\begin{table}[t!] %tab5
\caption{\label{tab5}The comparison the results of different
method~\cite{9,62,68,74,75,76,77}}
\begin{tabular}{cccc}
\thead
Solver & Method & FOS & Deviation (\%)\\
\endthead
\morerows{1}{OPTUM G2 (2022)} & \morerows{1}{SR-FELA} & 1.5--1.6 (lower
and upper bound) & $-$3.22 to $+$3.22\\
&& \textbf{1.55} ${\pm}$ \textbf{0.032 (exact solution)} & -\\
PLAXIS (2018) & FEM & \textbf{1.56} & $+$0.65\\
TOCHNOG (2022) & FEM & \textbf{1.56} & $+$0.65\\
FLAC (2019) & FDM & 1.63 & $+$5.16\\
KEM & KEM & 1.58 & $+$1.94\\
Bishop (Talren 2020) & Slices & 1.89 & $+$21.94\\
Janbu & Slices & 1.44 & $-$7.10\\
Cast3M & FEM & \textbf{1.56} & $+$\textbf{0.65}
\botline
\end{tabular}
\end{table}

The limit equilibrium methods for assessing the stability of earth
structures are now used routinely in practice. In spite of this
extensive use, the fundamentals of the methods are often not that well
understood and expectations exceed what the methods can provide. The
fact and implications that limit equilibrium formulations are based on
nothing more than equations of statics with a single, constant factor
of safety is often not recognized. A full appreciation of the
implications reveals that the method has serious limitations. The main
limitation is the need to guess the general form of the failure surface
in advance, with poor choices giving poor estimates of the failure
load. In practice, the correct form of the failure surface is often not
intuitively obvious, especially for problems with an irregular
geometry, complex loading, or complicated stratigraphy~\cite{32,33}.
There are other shortcomings of the technique, as follows.
\begin{enumerate}[(a)]
\item The resulting stresses do not satisfy equilibrium at every point
in the domain. 
\item There is no simple means of checking the accuracy of the
solution. 
\item It is hard to incorporate anisotropy, discontiniuty and
inhomogeneity. 
\item The analysis does not take into account the dilatancy of the
materials and does not directly obtain the stresses and displacements
in the structure. 
\item The failing soil mass is divided into slices for the calculation.
The accuracy of the result depends on the width of each slice. A thin
slice width allows to have a valid result. 
\item It is difficult to generalise the procedure from two to three
dimensions.
\item It is now possible to deal highly irregular porewater pressure
conditions, a variety of linear and nonlinear shear strength models,
virtually any kind of slip surface shape, concentrated loads, and
structural reinforcement.
\end{enumerate}

However, owing to the inherent assumptions of limit equilibrium
analyses, this method does not always furnish a unique factor of safety
and it is therefore unsuitable for generating a reference solution for
assessing the accuracy of alternative methods. To use limit equilibrium
methods effectively, it is important to understand and comprehend the
inherent limitations. The c-phi reduction method has been used in
Plaxis and Tochnog software. c-phi reduction is an option available in
Plaxis~\cite{68} to compute safety factors. In the phi-c reduction
approach the strength parameters $\tan\varphi$ and $c$ of the soil are
successively reduced until failure of the structure occurs. The total
multiplier $\Sigma Msf$ is used to define the value of the soil
strength parameters at a given stage in the analysis:
{\begin{eqnarray}
&& \label{eq18} \Sigma Msf = \tan(\varphi'_{\mathrm{input}}) / \tan
(\varphi'_{\mathrm{reduced}}) = c'_{\mathrm{input}} /
c'_{\mathrm{reduced}}.\Seqnsplit
&& \label{eq19} \text{SF} = \text{available strength/strength at
failure} = \text{value of}\ \Sigma Msf\ \text{at failure}.
\end{eqnarray}}\noindent

A feature to OptumG2~\cite{74} is the ability to compute rigorous upper
and lower bounds on the factor of safety. This allows the user to
bracket the exact solution from above and below. More precisely, if we
denote the exact solution by $E$ and the lower and upper bounds by $L$
and $U$ respectively, the following inequalities hold: $L \leq E \leq
U$. Moreover, if we denote the mean value of the upper and lower bounds
by $M = (L + U)/2$, we may define an absolute error by
$\varepsilon_{\mathrm{abs}} = (M - L) = (U - M)$. 

So that the exact solution may be expressed as: 
{\begin{equation}\label{eq20}
E = M \pm \varepsilon_{\mathrm{abs}}\ [\%].
\end{equation}}\noindent

The comparison of Cast3M result with different methods commonly used in
geotechnics and exact solution (gap: 0.65\%) shows that the
implementation of c-phi reduction method is successful.

\subsubsection{Phase 4: elasto-plastic modeling} \label{sec5.4.4}

In this phase~4, the retaining wall~2 is built on first layered
backfill. This phase starts after phase~3.1, the displacements are
initialized to zero in the model and is based on the application of
geostatic stress (gravity loading)\vadjust{\pagebreak} and the stresses
add up from one phase to the following.{\break}
Figures~\ref{fig27}--\ref{fig29} present the total displacements,
vertical plastic strains and vertical stress (in phase~4) of the model
from elastoplastic analysis. The construction of the retaining wall~2
above the first layered backfill generates minor settlements on it
without prejudice to the stability of the structure.

\begin{figure}
\includegraphics{fig27}
{\vspace*{-.25pc}}
\caption{\label{fig27}Total displacements (m) in SLS of the model in
phase~4.}
{\vspace*{.35pc}}
\end{figure}

\begin{figure}
\includegraphics{fig28}
{\vspace*{-.25pc}}
\caption{\label{fig28}Vertical plastic strains in SLS in phase~4.}
\end{figure}

\begin{figure}
\includegraphics{fig29}
{\vspace*{-.45pc}}
\caption{\label{fig29}Shadings of the vertical stress (Pa) of the model
in ULS (phase~4).}
{\vspace*{.4pc}}
\end{figure}

\subsubsection{Phase 5: elasto-plastic modeling} \label{sec5.4.5}

In this phase~5, the 2rd layered backfill is realized on first layered
backfill. This phase starts after phase~4, the displacements are
initialized to zero in the model and is based on the application of
geostatic stress (gravity loading). Figures~\ref{fig30} and~\ref{fig31}
present the total displacements and vertical stress (in phase~5) of the
model from elastoplastic analysis. The construction of the layered
backfill~2 above the first layered backfill, generates a negligible
horizontal displacement of the wall~2 and minor settlements on the
latter without prejudice to the stability of the structure. The~thrust
of the wall is controlled by the construction of an additional stop
corresponding to the projected horizontal displacement due to the
thrust of the retaining wall~2.

\begin{figure}
\includegraphics{fig30}
{\vspace*{-.25pc}}
\caption{\label{fig30}Deformed of the total incremental displacements
(m) in SLS of the model (displacements scaled by factor of~320.0) in
phase~5.}
{\vspace*{.2pc}}
\end{figure}

\begin{figure}
\includegraphics{fig31}
{\vspace*{-.25pc}}
\caption{\label{fig31}Shadings of the vertical stress (Pa) of the model
in ULS (phase~5).}
\end{figure}

\subsubsection{Phase 6: elasto-plastic modeling} \label{sec5.4.6}

In this phase~6, the application of surface load for road construction
following the GTR recommendation~\cite{71} (20~kN/$\text{m}^{2}$,
regarded as variable in Figure~\ref{fig12}) is realized on 2rd layered
backfill. This phase starts after phase~5, the displacements are
initialized to zero in the model and is based on the application of
surface load and the stresses add up from one phase to the following.
Figures~\ref{fig32}--\ref{fig35} present the total displacements,
plastic strains and vertical stress (in phase~6) of the model from
elastoplastic analysis. The total displacements obtained are small
(0.56~cm), these are quickly canceled by the creation of extra
embankment whose thickness is equal to the provisional displacement of
this phase.

\begin{figure}
\includegraphics{fig32}
\caption{\label{fig32}Shadings of the deformed of the total incremental
displacements (m) in SLS of the model (displacements scaled by factor
of 320.0) in phase~6.}
\end{figure}

\begin{figure}
\includegraphics{fig33}
{\vspace*{-.25pc}}
\caption{\label{fig33}Horizontal plastic strains in SLS in phase~6.}
{\vspace*{.5pc}}
\end{figure}

\begin{figure}
\includegraphics{fig34}
{\vspace*{-.25pc}}
\caption{\label{fig34}Vertical plastic strains in SLS in phase~6.}
\end{figure}

\begin{figure}
\includegraphics{fig35}
\caption{\label{fig35}Shadings of the vertical stress (Pa) of the model
in ULS (phase~6).}
\end{figure}

\subsubsection{Phase 7: elasto-plastic modeling and strength reduction
method} \label{sec5.4.7}
\paragraph{Phase 7.1: elasto-plastic modeling.} \label{sec5.4.7.1}

In this phase~7.1, the entire geotechnical structure is built in a
single phase, there are no previous phases. The gravity load and
surface load are applied on the model. Figure~\ref{fig36} presents the
total displacements (in phase~7.1) of the model from elastoplastic
analysis. The total displacement in a single phase is 49~cm. This
calculation is complementary in order to have the difference in term of
displacements if the structure were built in one phase. Which does not
represent the practical reality. On the other hand, it is to be
compared with that of phase~0 (part~\ref{sec5.3}), carried out in
linear elastic behavior.

\begin{figure}
\includegraphics{fig36}
\caption{\label{fig36}Shadings of the total displacements (m) in SLS of
the model in phase~7.1.}
\end{figure}

\paragraph{Phase 7.2: strength reduction method.} \label{sec5.4.7.2}

In this phase~7.2, the gravity load and surface load are applied on the
model. This phase starts from phase~7.1. Figures~\ref{fig37}
and~\ref{fig38} present the horizontal displacements (in phase~7.2) and
the shadings of shear total incremental displacements of the model from
c-phi strength reduction method (ULS, GEO-3) analysis and displayed
adequate failure kinematics. The factor of safety for this phase is
equal to 1.72 guaranteeing the stability of the geotechnical structure
despite the shear displacements of 73~cm.

\begin{figure}
\includegraphics{fig37}
{\vspace*{-.15pc}}
\caption{\label{fig37}Shadings of horizontal displacements (m) in ULS
of the model in phase~7.2.}
{\vspace*{.8pc}}
\end{figure}

\begin{figure}
\includegraphics{fig38}
\caption{\label{fig38}(a)~Shadings of shear total incremental
displacements (m); (b)~shadings of deformed of the total shear
incremental displacements (m) in ULS of the model (displacements scaled
by factor of~2) in phase~7.2.}
{\vspace*{-.3pc}}
\end{figure}

In this phase, Plaxis software has been used (using Mohr--Coulomb model
for soils) for the comparison of the results (Figure~\ref{fig39}). The
factor of safety for this phase obtained by Plaxis software is equal
to~1.61. The Cast3M results are generally in good agreement with the
Plaxis results of the factor of safety. The difference between the two
safety factor values is 6.83\%. This~gap can be justified by the fact
that the very fine mesh of Plaxis only has 818 elements at 15-noded
(not possible to refine beyond) whereas the global model of Cast3M has
1800 elements TRI3.

\begin{figure}
\includegraphics{fig39}
\caption{\label{fig39}Shadings of shear total incremental displacements
(m) calculated by Plaxis.}
\end{figure}

Table~\ref{tab6} presents the results of the displacements and the
Factor of Safety on the model phase by phase. The displacements of the
elasto-plastic calculation are higher compared to the elastic
calculation (Table~\ref{tab6}). In this paper, our computations were
carried out respecting the practical phasing of construction of the
geotechnical structures. The advantage of staged construction is to
know the displacements generated after the construction of each phase.
This technique makes it possible to build each phase by providing on
the one hand the extra embankment whose extra elevation corresponds to
the provisional settlements of this phase and on the other hand an
additional buttress if necessary so the thickness corresponds to the
provisional horizontal displacement caused on the retaining wall.

\begin{table}[t!] %tab6
\fontsize{9.2}{11.2}\selectfont
\caption{\label{tab6}Summary of the results of the displacements by
phase}
\begin{tabular}{ccccc}
\thead
\parbox[t]{1.5cm}{\centering Calculation phase} & Description of phase
& \parbox[t]{1.5cm}{\centering Calculation type} &
\parbox[t]{2cm}{\centering Total displacement~(cm)} &
\parbox[t]{1.15cm}{\centering Factor of safety}\vspace*{2pt}\\
\endthead
Phase 0 & Initialization of initial stresses & Linear elastic &
\textbf{21} & -\vspace*{3pt}\\
Phase 1 & \parbox[t]{5.3cm}{\centering Application of geostatic stress in
the Subsoil and basement}\vspace*{5pt} & Elasto-plastic & 3.22 & -\\
Phase 2 & \parbox[t]{5.3cm}{\centering Application of geostatic stress of
the retaining wall~1}\vspace*{3.5pt} & Elasto-plastic & 8.83 & -\\
Phase 3.1 & \parbox[t]{5.3cm}{\centering Application of geostatic stress
of the layered backfill~1}\vspace*{3.5pt} & Elasto-plastic & 15 & -\\
Phase 3.2 & \parbox[t]{5.3cm}{\centering Gradual reduction of the
mobilized strength}\vspace*{3.5pt} & c-phi reduction & 71 &
\textbf{1.56}\\
Phase 4 & \parbox[t]{5.3cm}{\centering Application of geostatic stress of
the retaining wall~2}\vspace*{3.5pt} & Elasto-plastic & 0.24 & -\\
Phase 5 & \parbox[t]{5.3cm}{\centering Application of geostatic stress of
layered backfill~2}\vspace*{3.5pt} & Elasto-plastic & 0.087 & -\\
Phase 6 & \parbox[t]{5.3cm}{\centering Application of the surface
load}\vspace*{3.5pt} & Elasto-plastic & 0.56 & -\\
Phase 7.1 & \parbox[t]{5.3cm}{\centering The gravity load and surface
load are applied on the model in single phase}\vspace*{4pt} &
Elasto-plastic & \textbf{49} & -\\
Phase 7.2. & \parbox[t]{5.3cm}{\centering Gradual reduction of the
mobilized strength on the model}\vspace*{2pt} & c-phi reduction & 73 &
\textbf{1.72}
\botline
\end{tabular}
\end{table}

After construction of the phase, the displacements are obtained and the
structure returns to its initial position of stability provided during
the design. From this staged construction modeling, it becomes easy to
control and limit settlement according to the functional requirements
of the geotechnical structure. 

The safety factor obtained after the construction of the entire
geotechnical structure (phase~7) is 1.72, therefore higher than the
1.56 obtained after the construction of backfill~1 (phase~3.2). This
increase is due to the fact that, backfill~2 has a cohesion value 10
times greater than that used on backfill~1. These materials come from
two different quarries. Backfill~1 consits of materials with a strong
granular predominance (greyish pozzolan) with a low fines content.
Backfill~2 is constituted by the materials with a spread out grain size
(Reddish clay gravel at gravelly sand).

\section{Discussion} \label{sec6}

The implementation of a c-phi reduction procedure in Cast3M shows that
this approach, which is often used as a black box in geotechnical
softwares, can be implemented in rigorous way. This procedure is first
validated on a simple example of bearing capacity of a strip footing
and then used to evaluate the stability factor of a geotechnical
structure composed of two levels of backfill retained by two retaining
walls, with a road built on the top of the second backfill.

\subsection{Footing results} \label{sec6.1}

Concerning the modeling results of the strip footing, the difference in
the theoretical results obtained (formulas~(\ref{eq9})
and~(\ref{eq11})) is in agreement with the assumptions made for each
method for the parameter $N\gamma$. The ``exact'' theoretical
expression of $N\gamma$ is not known. Only the upper and lower bounds
have been determined in particular by Vesic~\cite{63},
Brinch-Hansen~\cite{64,65} and Meyerhof~\cite{66,67}.
Brinch-Hansen~\cite{64,65,66} indicates its value of $N\gamma$ is a
lower bound and that the addition of the three bearing capacity terms,
which does not correspond exactly to the real failure mechanism,
corresponds to an underestimation of the bearing capacity which is
generally less than 20\%. It is also this lower bound modified, which
is used in Meyerhof's formula. For the Vesic formula, the upper bound
of $N\gamma$ is used. It indicates that the addition of the three
bearing capacity terms gives an approximative estimate with less than
10\% error over the interval $15 < \varphi < 45\mbox{\textdegree}$.

The bearing capacity value of 101~kPa obtained by Cast3M modelling with
a very fine mesh is regarded as the exact solution because included
between the theoretical values obtained respectively by the lower and
upper bounds. For very fine mesh of the numerical model of the footing
(validation on a simple example of bearing capacity of a strip
footing), the difference between Plaxis and Cast3M results is 7\%
(Table~\ref{tab2}). This gap comes the difference between
Drucker--Prager and Mohr--Coulomb materials. The analysis on the mesh
density and element type (TRI3 and TRI6) in Cast3M (Table~\ref{tab3})
shows that, the difference between element type TRI3 and TRI6 results
varies between 30 (for the very coarse mesh with local refinements
under footing) to 0\% (for the very fine mesh). The results of
Table~\ref{tab3} show clearly that the element type, the mesh density
and the convergence tolerances have a pronounced influence on the
factor of safety obtained from the displacement finite-element method.
Therefore, the influence of these parameters must be minimised by the
use of high-order elements, very fine meshes and stringent tolerances
in order to obtain a realistic results.

\subsection{Retained walls--slope stability results} \label{sec6.2}

The Factor of Safety calculation procedure is based on the calculation
of a stable initial stress and corresponding deformation state
recursive reduction of the shear and tensile strength until a
convergence to a prescribed tolerance is still obtainable. With regard
to the retained walls--slope stability and backfill results at the
phase~3.2 (part~\ref{sec5.4.3.2}), the form of the failure surface
obtained from Bishop and Janbu methods is very different from the other
methods (Figure~\ref{fig26}). The same is true for values of resulting
safety factor (Table~\ref{tab5}) which deviate from the value which
taken for reference~1.55, against 1.89 for Bishop and 1.44 for Janbu.
The Factor of Safety obtained by the FE methods (Cast3M, PLAXIS and
TOCHNOG) is identical (1.56) and between the upper and lower bounds of
OPTUM G2 (1.5 and 1.6). On the other hand, the value of 1.58 obtained
by the Kinematic Element Method is an upper bounds and corresponds well
to $1.55+0.03$ obtained by OPTUM G2 (Table~\ref{tab5}). In view of
these results, the Janbu method seems the most conservative. The
difference in the results obtained is in agreement with the initial
assumptions made for each method as follows.
\begin{enumerate}[(a)]
\item For the limit equilibrium analysis, the need to guess the
general form of the failure surface in advance makes it impossible to
obtain a rigorous solutions. The poor choices giving poor estimates of
the failure load. In practice, the correct form of the failure surface
is often not intuitively obvious, especially for problems with an
irregular geometry, complex loading, or complicated
stratigraphy~\cite{32}. The resulting stresses do not satisfy
equilibrium at every point in the domain. The analysis does not take
into account the dilatancy of the materials and does not directly
obtain the stresses and displacements in the structure. The failing
soil mass is divided into slices for the calculation. The accuracy of
the result depends on the width of each slice. 

\item For the Kinematic Element Method, the theorem used is the
approach of limit loads by upper bounds. It consists in constructing
kinematics respecting boundaries conditions on displacements and in
finding for each of them forces that are too great for the
\mbox{resistance} of the soil~\cite{51}. This method overcomes the
limits of the slice method in \mbox{rigorous} framework. It allows
taking into account any geometries and stratigraohy, all flow
conditions, point or distributed overloads in any direction and any
type of \mbox{reinforcement}.

\item The finite-difference analysis is based in the explicit
Lagrangian calculation scheme. The timestep of the calculation, the
choice of the elements of the mesh and the numerical convergence have a
pronounced influence on the factor of safety obtained from this method.
Therefore, the influence of these parameters must be minimised by the
use very small timestep of the calculation, convenable elements and
stringent tolerances to obtain a realistic results.

\item The finite-element analysis is routinely applied for assessing
displacements, stresses and strains for working load conditions, this
technique is increasingly being used to calculate ultimate limit states
and, consequently, factors of safety. In the finite-element codes the
factor of safety is obtained by means of the strength reduction method;
that is, an analysis is performed with mobilised strength properties
for the friction angle and the cohesion~$c$, followed by an incremental
decrease of $\tan\varphi'$ and $c'$ (assuming a Mohr--Coulomb failure
criterion). This gives stress states that violate the strength criteria
which are resolved in an iterative manner using the same stress point
algorithm employed for a standard elasto-plastic analysis in
finite-element codes, leading to a stress redistribution in the system
until equilibrium can no longer be established and failure is
reached~\cite{32}. However, close inspection of the developed failure
mechanism and displacements of appropriate control points is required
in order to avoid misinterpretation. It is well known that the element
type, the mesh discretisation and the convergence tolerances have a
pronounced influence on the factor of safety obtained from the
displacement finite-element method. Therefore, the influence of these
parameters must be minimised by the use of high-order elements, fine
meshes and stringent tolerances in order to obtain a realistic results.
Our numerical modelling were made using a very fine mesh in Cast3M
models.

\item The new emerging method so-called finite-element limit analysis
(FELA) is particularly powerful when both upper and lower bound
estimates are calculated so that the true collapse load is bracketed
from above and below. The difference between the two bounds then
provides an exact measure of the error in the solution, and can be used
to refine the meshes until a suitably accurate estimate of the collapse
load is found~\cite{32}.
\end{enumerate}

The originality of our computation method is based on the adjustment of
the factor FCUT allowing to regulate numerical gap between the very
coarse mesh effect to very fine mesh of the model.

\section{Conclusion} \label{sec7}

For geotechnical problems: Serviceability Limit State (SLS),
computations are based on the deformation analysis. Ultimate Limit
State (ULS): computations for stability analysis Strength reduction
method are a powerful numerical technique implementing different safety
concepts. It is argued that the Finite Element method of retaining wall
and slope stability analysis is a more powerful alternative to
traditional limit equilibrium methods and its widespread use should now
be standard in geotechnical practice. Cast3M is an FE code dedicated
mainly to the calculation of structures. At present, from a methodical,
rigorous and sophisticated computation, we manage to solve complex
geotechnical problems both at Serviceability Limit State (SLS) and the
Ultimate Limit State analysis (ULS) by respecting the functional
requirements and the different stages of construction allowing to have
the overall stability which guarantees the durability of the
geotechnical structure in comfortable safety conditions. The results
presented in this paper were compared with those obtained by commercial
industrial software dedicated to geotechnical problems and by reference
solutions testifying to the reliability of our rigorous computation
from the Cast3M code. A unique feature to Cast3M is the ability to
compute rigorous factor of safety and analyses deliver adequate failure
kinematics. The main contributions of this paper are the following:
\begin{itemize}
\item The implementation of a c-phi reduction procedure in Cast3M shows
that this approach, which is often used as a black box in geotechnical
softwares, can be implemented in rigorous way.
\item The comparison with different methods commonly used in
geotechnics shows that the implementation of c-phi reduction method is
successful.
\item The originality of the method is based on the adjustment of the
factor FCUT allowing to regulate numerical gap between the very coarse
mesh effect to very fine mesh of the model.
\end{itemize}

\section*{Declaration of interests}

The authors confirm that there are no conflicts of interest associated
with this publication and there has been no significant financial
support for this work that could have influenced its outcome.

\vspace*{-.3pc}

\back{}

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

\end{document}