\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2012 (2012), No. 167, pp. 1--9.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2012 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document}
\title[\hfilneg EJDE-2012/167\hfil Peak load evaluation of structural elements]
{An algorithmic approach for peak load evaluation of structural elements
obeying a Men\'etrey-Willam type yield criterion}

\author[A. A. Pisano \hfil EJDE-2012/167\hfilneg]
{Aurora Angela Pisano}  % in alphabetical order

\address{Aurora A. Pisano \newline
Department PAU, University Mediterranea of Reggio Calabria\\
Feo di Vito, 89100 Reggio Calabria, Italy}
\email{aurora.pisano@unirc.it}

\thanks{Submitted July 9, 2012. Published September 28, 2012.}
\subjclass[2000]{65D25, 74S05, 74C05}
\keywords{Numerical methods; finite element methods; rate-independent theory}

\begin{abstract}
 A numerical procedure for the evaluation of the peak load of a general
 class of structures is presented. Namely the addressed structures obey
 to a three stress invariants constitutive criterion of Men\'etrey-Willam-type
 endowed with cap in compression. The formulation, developed in the
 Haigh-Westergaard coordinates, allows an easy treatment of the mechanical
 problem in a 3D context as well as an interesting geometrical interpretation
 in the principal stress space.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{remark}[theorem]{Remark}
\allowdisplaybreaks

\section{Introduction}

 The broadest area of success of plasticity theory with concrete structures
 is the treatment of situations in which the material acts primarily in compression.
 The usual procedure is nowadays to apply plasticity theory in the compression zone
 and treat the zones in which at least one principal stress is tensile by one of several
 version of fracture and/or damage mechanics also with the use of non local internal
 variables, see e.g. \cite{JiBa02}. To guarantee a reliable design is essential to possess
 codes able to catch the behaviour of concrete structures even above the elastic limit;
 post elastic step-by-step analyses are the common tools available to this aim.
 However, the obtainable results,  in the case of concrete elements, are strongly
 influenced by the dependence of strength, stiffness and ductility of this material
 on the load path. A valid alternative can be given by direct methods, i.e. the methods
 that are able to predict the load bearing capacity of a structure, or a structural element,
 in terms of peak load value, the one corresponding to the exhaustion of the strength
 capabilities, so avoiding any evolutive analysis.
 
 A number of recent contributions are given in \cite{WePo09} and references therein.
The numerical procedure for limit analysis, here adopted and known as
\emph{Linear Matching Method} (LMM), is the one
originally proposed by Ponter and Carter  \cite{PoCa97} for Von Mises materials
and recently generalized in Pisano and Fuschi \cite{PiFu07,PiFu11},
 Pisano et al \cite{PiFuDD12}, with reference
to Tsai-Wu-type anisotropic, non associative composite plates. 
The method is here deeply reformulated in the Haigh-Westergaard coordinates 
with reference to a 3D-plasticity model for concrete
with a pressure-sensitive yield surface with curved meridians in the Rendulic section,
$J{_3}$ dependence, cap in compression and a non associated flow rule.
The yield surface turns into the failure surface proposed by Men\'etrey and Willam \cite{MeWi95},
calibrated from elementary strength data and including the standard strength hypotheses
of Huber-Mises, Drucker-Prager, Rankine, Mohr-Coulomb and Leon as special cases. The
assumed yield surface is given in terms of three independent stress invariants and,
in principal stress space it is convex and smooth. The LMM is applied to compute an upper bound
to the peak load multiplier for simple concrete elements and one academic example
is considered to estimate the gap between the computed upper bounds and the corresponding
theoretical peak load values.

\section{Constitutive model and theoretical background}

The adopted three dimensional constitutive model for concrete is due
 to Men\'etrey and Willam \cite{MeWi95}.
It provides a three parameter failure surface having the following  expression:
\begin{equation} \label{P_F1}
    f(\xi,\rho,\theta)=[\sqrt{1.5}\frac{\rho}{f_c'}]^2+
    m[\frac{\rho}{\sqrt{6}f_c'}\, r(\theta,e)+\frac{\xi}{\sqrt{3}f_c'}]-c=0,
\end{equation}
where
\begin{gather} \label{P_F2}
    r(\theta,e)=\frac{4(1-e^2)\cos^2\theta+(2e-1)^2}
   {2(1-e^2)\cos\theta+(2e-1)\Bigl[4(1-e^2)\cos^2\theta+5e^2-4e\Bigr]^{1/2}},
\\ \label{P_F3}
     m:=3\frac{{f_c'}^2-{f_t'}^2}{f_c'\,f_t'}\,\frac{e}{e+1}.
\end{gather}

Equation \eqref{P_F1} is expressed in terms of three stress invariants $\xi,\rho$ 
and $\theta$ known as the Haigh-Westergaard coordinates; $m$ is the friction 
parameter of the material depending on
the uniaxial compressive strength $f_c'$, on the uniaxial tensile strength 
$f_t'$ as well as on the eccentricity parameter $e$.
The eccentricity $e$ defines the smoothness of the Men\'etrey-Willam surface 
and its value strongly influences
the failure description either in biaxial tension or in compression.
In \eqref{P_F1}, $\xi$ is the hydrostatic stress invariant, 
$\rho$ is the deviatoric stress invariant, and $\theta$ is the deviatoric 
polar Lode angle; the Haigh-Westergaard (H-W) coordinates are given by
\begin{gather}
 \xi =\frac{1}{\sqrt{3}}I_{1} , \quad I_{1} =\sigma_{ii}\,, \label{P_F4} \\
 \rho =\sqrt{2\, J_{2}} , \quad J_{2} =\frac{1}{2}s_{ij}s_{ij}\,, \label{P_F5} \\
 \cos3\theta =\frac{3\sqrt{3}}{2}\frac{J_3}{J_{2}^{3/2}} ,
\quad  J_3=\frac{1}{3}s_{ij}s_{jk}s_{ki}\,, \label{P_F6}
\end{gather}
with $s_{ij}$ denoting the deviatoric stress components; i.e.,
$s_{ij}=\sigma_{ij}-\frac{1}{3}\sigma_{kk}\delta_{ij}$ being
$\delta_{ij}$ the Kronecker symbol. It is worth to note that,
for $0\leq\theta\leq\frac{\pi}{3}$, the following relations,
between principal stresses of $\sigma_{ij}$ and H-W coordinates,
hold:
\begin{equation}\label{P_F7}
\begin{Bmatrix}
\sigma_{1}\\
\sigma_{2}\\
\sigma_3
\end{Bmatrix}
= \frac{1}{\sqrt{3}}
\begin{Bmatrix}
\xi\\
\xi\\
\xi
\end{Bmatrix}
+\sqrt{\frac{2}{3}}\rho\,
\begin{Bmatrix}
\cos\theta\\
\cos(\theta-\frac{2\pi}{3})\\
\cos(\theta+\frac{2\pi}{3})
\end{Bmatrix}.
\end{equation}


\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.9\textwidth]{fig1}
\end{center}
\caption{Triaxial failure surface of Men\'etrey-Willam with cap:
(a) deviatoric sections at three generic values of hydrostatic pressure;
(b) compressive and tensile meridians in the Rendulic plane}
\label{fig1} \end{figure}


The sector $0\leq\theta\leq\frac{\pi}{3}$ confining the values
of the Lode angle $\theta$ given by \eqref{P_F6}$_1$ or, equivalently, by
$\cos\theta=(2\sigma_{1}-\sigma_{2}-\sigma_3)
/2\sqrt{3J_{2}}$, is consequence of a regular ordering
of the eigenvalues of the stress tensor $\sigma_{ij}$; i.e.,
$\sigma_{1}\geq\sigma_{2}\geq\sigma_3$, as assumed in the
following. Remembering that the Lode angle $\theta$ is the
deviatoric projection of the angle between the radius vector
of the current stress point in the principal stress space and
the axis $\sigma_{1}$, the arbitrariness in the ordering of the
eigenvalues of a tensor implies that six different points (i.e.
six $\theta$ values satisfying \eqref{P_F6}$_1$) correspond in the H-W
representation to a given stress tensor.
As a result, the M-W surface has to be symmetric with respect
to the projections of the principal axes on the deviatoric plane.
The surface represented by \eqref{P_F1}, defined for
$0\leq\theta\leq\frac{\pi}{3}$, so extends to all polar directions
$0\leq\theta\leq 2\pi$ using the three-fold symmetry shown in 
Figure \ref{fig1}(a).
The above condition is obviously met by the postulated isotropy
of the material for which the labels 1,2,3 attached to the coordinate
axes are arbitrary.

With the assumed ordering of the eigenvalues of the stress tensor there are
two extreme cases:
\begin{align}\label{P_F8}
\sigma_{1}=\sigma_{2}>\sigma_3; \quad \sigma_{1}>\sigma_{2}=\sigma_3;
\end{align}
corresponding to $\theta=\frac{\pi}{3}$ and $\theta=0$, respectively.
The meridian at $\theta=\frac{\pi}{3}$ is called \emph{compressive
meridian} (see also Figure \ref{fig1}(b)) in that \eqref{P_F8}$_1$ represents a stress state
corresponding to an hydrostatic stress state with a compressive stress
superimposed in one direction. The meridian determined by $\theta=0$
 (see again Figure \ref{fig1}(b)), represents an hydrostatic
stress state with a tensile stress superimposed in one direction and
is therefore called the \emph{tensile meridian}.
A realistic representation of concrete, viewed as a frictional
material, requires to take into account the dilatancy, i.e. the
volumetric expansion under compression. To this aim a \emph{non associated
flow rule} has to be \emph{postulated}. Moreover, to capture the concrete
behaviour near the hydrostatic axis, \emph{a cap is adopted} to ``close''
the surface given by \eqref{P_F1}. This cap, in the H-W coordinates has the shape
\begin{equation}\label{P_F9}
\begin{gathered}
\rho^{CAP}(\xi,\theta) =-\frac{\rho^{MW}(\xi_{a},\theta)}{(\xi_{a}-\xi_{b})^2}
[\xi^2-2\xi_{a}(\xi-\xi_{b})-\xi_{b}^2]\\
\textrm{with } \xi_{b}\leq\xi\leq\xi_{a},\quad 
0\leq\theta\leq\frac{\pi}{3} 
\end{gathered}
\end{equation}
where $\rho^{MW}(\xi,\theta)$ is the explicit form of the
parabolic meridian of the M-W surface that, looking at \eqref{P_F1},
can be given by the expression
\begin{equation}\label{P_F10}
\begin{gathered}
\rho^{MW}(\xi,\theta)=\frac{1}{2a}\{ -b(\theta)
+[b^2(\theta)-4a\,c(\xi)]^{1/2}\} \\
 \textrm{with } 
\xi_{a}\leq\xi\leq\xi_{v}, \quad 
0\leq\theta\leq\frac{\pi}{3}
\end{gathered}
\end{equation}
where:
\begin{equation}\label{P_F11}
a:=\frac{1.5}{(f_c')^2}, \quad 
b(\theta):=\frac{m}{\sqrt{6}f_c'}r(\theta,e), \quad 
c(\xi) :=\frac{m}{\sqrt{3}f_c'}\xi-1.
\end{equation}

The three-fold symmetry of the M-W surface is maintained and the cap surface,
given by \eqref{P_F9}, extends to all polar directions $0\leq\theta\leq 2\pi$
as sketched in Figure \ref{fig1}. The values $\xi_{a}$ and $\xi_{b}$ appearing in
\eqref{P_F9} locate the cap position that can be detected experimentally
as given by Li and Crouch \cite{LiCr10}.

 The Men\'etrey-Willam surface with a cap in compression  is here assumed
as yield criterion for the material,
while the lack of associativity is overcome by means of the non-standard
limit analysis approach \cite{Ra61}.

\section{The kinematic approach of limit analysis and the LMM}


\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.9\textwidth]{fig2}
\end{center}
\caption{Stress point at yield $\mathbb{P}^Y$:
(a) deviatoric section at $\xi=\xi^Y$;
(b) Rendulic section at $\theta=\theta^Y$}
\label{fig2} \end{figure}


To set the problem, let consider a body of volume $V$  and   external 
surface $\partial V$ both referred to a
Cartesian co-ordinate system ($x_{i}\,$, $i=1,2,3$).
The body is, for simplicity, subjected only to surface loads
$P\bar{\mathbf{p}}(\mathbf{x})$, where: $P$ is a scalar
load multiplier; $\bar{\mathbf{p}}(\mathbf{x})$ denotes
the reference load vector collecting all the surface force components,
$\bar{p}_{i}$, acting on points of a portion of the body surface,
namely $\partial V_t$. The remaining part of $\partial V$, namely
$\partial V_{u}=\partial V-\partial V_t$, is assumed to suffer
displacements $\mathbf{u}=\mathbf{0}$.
Following the directions of non-standard limit analysis,
reference can be made to a standard material "bounding from above"
the real non-standard one and obeying to the same yield criterion.
Rewriting the M-W-type yield criterion  in the abridged form
$F(\xi,\rho,\theta)=0$, the volumetric and deviatoric strain
rate components at collapse
can be expressed in the form
 \begin{equation}\label{S_F1}
\boldsymbol{\dot{\varepsilon}}_{v}^{c}=\dot{\lambda}
\frac{\partial F}{\partial \xi} ;\quad
\boldsymbol{\dot{\varepsilon}}_{d}^{c}=\dot{\lambda}
\frac{\partial F}{\partial \rho} ;
\end{equation}
with $\dot{\lambda}\geq 0$ scalar multiplier. The kinematic approach of the limit
analysis theory states that
\begin{equation}\label{S_F2}
P_{_{UB}}\,\int_{\partial V_t}\bar{p}_{i}\,\dot{u}_{i}^{c}\,\textrm{d} (\partial V)
=\int_{V}(\xi^{Y}\,\dot{\varepsilon}_{v}^{c}+
\rho_{x}^{Y}\,\dot{\varepsilon}_{d_{x}}^{c}+\rho_{y}^{Y}\,
\dot{\varepsilon}_{d_{y}}^{c})\,\textrm{d}V .
\end{equation}
In this equation, $P_{_{UB}}$ denotes the searched upper bound to the peak 
load multiplier; $\dot{u}_{i}^{c}$ are the displacement rates at collapse 
of the points
where surface loads, $\bar{p}_{i}$, act. Moreover, $\dot{u}_{i}^{c}$ satisfy
the condition $\dot{u}_{i}^{c}=0$ on $\partial V_{u}$ and
are compatible with the strain rates at collapse.
Finally, $\xi^{Y}$, $\rho_{x}^{Y}$,
$\rho_{y}^{Y}$ are the stress components at yield associated to the
volumetric and deviatoric strain rate components at collapse
$\dot{\varepsilon}_{v}^{c}$, $\dot{\varepsilon}_{d_{x}}^{c}$,
$\dot{\varepsilon}_{d_{y}}^{c}$, respectively.
Referring
to Figures \ref{fig2}(a)(b) an orthogonal cartesian system ($x ,y$), centered
at $\xi=\xi^{Y}$, associated to the cylindrical coordinates
$(\xi,\rho,\theta)$ has been introduced for simplicity,
$(\cdot)_{x}$, $(\cdot)_{y}$ denoting
components of $(\cdot)$ along $x$ and $y$ respectively.
The \emph{set} $(\dot{\varepsilon}_{v}^{c},
\dot{\varepsilon}_{d_{x}}^{c}, \dot{\varepsilon}_{d_{y}}^{c},
\dot{u}_{i}^{c})$, if given \emph{at all stress point
at yield} (such as $\mathbb{P}^{Y}$ in Figures \ref{fig2}(a)(b)), defines a collapse
mechanism for the analyzed structure.
The LMM in practice furnishes all the required quantities to compute a $P_{_{UB}}$.
The key idea is to build such quantities assuming the analysed structural element as made by a
\emph{linear viscous fictitious} material.
This material is fictitious in the sense that it
possesses \emph{spatially varying moduli} and is also subjected  
to \emph{imposed initial stresses}.
A Finite Element (FE) analysis is performed on such fictitious structure so obtaining a fictitious
solution. The fictitious moduli and initial stresses are then adjusted
so that the computed fictitious stresses are
brought on the yield surface at a fixed strain rate distribution.
 This allows one to define a collapse
mechanism (compatible strain and displacement rates
$(\dot{\varepsilon}_{v}^{c},
\dot{\varepsilon}_{d_{x}}^{c}, \dot{\varepsilon}_{d_{y}}^{c},
\dot{u}_{i}^{c})$,
plus related stresses at yield  $\xi^{Y}$, $\rho_{x}^{Y}$,
$\rho_{y}^{Y}$) and,
eventually, allows to evaluate a value of $P_{_{UB}}$. 
Such procedure has to be carried on iteratively
to satisfy the equilibrium conditions altered by the adjusting of the moduli
 and imposed stresses.
A theoretical  sufficient condition for convergence \cite{PoFuEn00} assures
 the reliability of the final $P_{_{UB}}$;
i.e., the one computed at last iteration.

\subsection{Geometrical interpretation of the LMM}

 It is worth to remind the formal analogy existing between a linear viscous problem
and a linear elastic problem on the basis of which
the complementary dissipation rate function pertaining to the linear viscous problem
plays the same role of the complementary energy potential function, say $W$,
 pertaining to the elastic problem.
As a consequence, the fictitious linear solution
(in terms of rate quantities) can be computed as fictitious elastic solution 
(in terms of finite quantities).
In the H-W representation  and for the chosen initial values
(hereafter denoted by the apex $(^{0})$) of material parameters and initial stresses,
$W$ can be given by the expression
\begin{equation}\label{S_F3}
W^{(0)}(\xi,\rho)=\frac{(\xi-\bar{\xi}^{(0)})^2}{6K^{(0)}}+
\frac{(\rho-\bar{\rho}^{(0)})^2}{4G^{(0)}}=\bar{W}^{(0)}.
\end{equation}
In equation, $\bar{W}^{(0)}$ is a constant value, $K^{(0)}$ and $G^{(0)}$ are 
the initial elastic bulk and shear moduli respectively, while
$\bar{\xi}^{(0)},\bar{\rho}_{x}^{(0)},\bar{\rho}_{y}^{(0)}$ are the initial stresses.
The volumetric and deviatoric strain rates of the linear problem can then be 
computed as:
\begin{align}\label{S_F4}
\boldsymbol{\dot{\varepsilon}}_{v}^{\ensuremath{\ell}}
 =\frac{\partial W^{(0)}}{\partial \xi} ;  \quad
\boldsymbol{\dot{\varepsilon}}_{d}^{\ensuremath{\ell}}
 =\frac{\partial W^{(0)}}{\partial \rho}.
\end{align}

The fictitious linear solution is also given in terms of
compatible displacement rates of the points at which surface
loads act, say $\dot{u}_{i}^{\ensuremath{\ell}}$, and
in terms of stress values $\xi^{\ensuremath{\ell}},
\rho_{x}^{\ensuremath{\ell}},\rho_{y}^{\ensuremath{\ell}}$,
associated to the strain rate components
$\dot{\varepsilon}_{v}^{\ensuremath{\ell}},
\dot{\varepsilon}_{d_{x}}^{\ensuremath{\ell}},
\dot{\varepsilon}_{d_{y}}^{\ensuremath{\ell}}$, respectively.
The latter simply being the fictitious strain rates \eqref{S_F4}
referred to the orthogonal cartesian system ($x,y,\xi$) associated
to the cylindrical H-W coordinates.
Grounding on the above analogy, the fictitious linear analyses, involved by the LMM
procedure, can be carried on at each iteration as fictitious elastic analyses
performed by any commercial FE code, rendering
the whole procedure easy.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.8\textwidth]{fig3}
\end{center}
\caption{Construction of the prolate spheroid $W^{(*)}=$ const. at matching:
(a) deviatoric section, at $\hat{\xi}=\xi=\xi_M$, of the M-W-type 
 yield surface and of the spheroid    $W^{(*)}=$ const.;
(b) section on the plane belonging to the sheaf of axis $\hat{\xi}$ 
and passing through $\mathbb{P}_M$}
\label{fig3} 
\end{figure}


 The above analogy leads also to an
interesting geometrical interpretation of the LMM.
First of all, as
demonstrated in Ponter et al \cite{PoFuEn00}, a sufficient condition 
\emph{to assure convergence}
of the iterative procedure relies on a geometrical requirement; i.e.,
the \emph{complementary energy}
equipotential \emph{surface} of the fictitious material
\emph{at matching} (circumstance hereafter denoted by the apex $(^{(*)})$)
has indeed to touch
the yield surface at the matching point $\mathbb{P}_M$ but
\emph{must} otherwise \emph{lie outside the yield surface}
in stress space. In the specific case,
with reference to \eqref{S_F3},
the complementary equipotential surface is a
\emph{prolate spheroid}; i.e., an \emph{ellipsoid of revolution}
obtained by rotating a generatrix ellipse about its major axis.
As sketched in Figure \ref{fig3}(b), in the principal stress space, the semi-axes
lengths of this spheroid at matching (say $W^{(*)}=$ const.) are related 
to the \emph{unknown} values
of the updated fictitious elastic moduli $K^{(*)}$ and $G^{(*)}$.
Its origin, $\mathbb{C}$, is located by the \emph{unknown} values of the
updated initial hydrostatic pressure and deviatoric stresses
$\bar{\xi}^{(*)},\bar{\rho}_{x}^{(*)},\bar{\rho}_{y}^{(*)}$.
In order to satisfy either the matching or the sufficient condition
for convergence the elastic moduli and the initial stresses must
be updated, at each iteration, according to the following requisites
to be met by the complementary energy surface $W^{(*)}=$ const.:
\begin{itemize}
\item[(i)]
it must contain the matching point $\mathbb{P}_M$;

\item[(i)] it has to be tangent in $\mathbb{P}_M$
 to the M-W-type yield surface admitting at
$\mathbb{P}_M$ outward normal equal to
$\boldsymbol{\dot{\varepsilon}}^{\ensuremath{\ell}} \equiv
\boldsymbol{\dot{\varepsilon}}^{c}$;

\item[(ii)]
it has otherwise lie outside the yield surface.
The latter requisite can be easily accomplished by imposing
that the ellipse passes through points (outer to the yield surface
but ``below'' the tangent through $\mathbb{P}_M$) like
$\mathbb{P}_{a}$ and $\mathbb{P}_{v}$ at $\hat{\xi}=\xi_{a}$
and $\hat{\xi}=\xi_{v}$, respectively, shown in Figure \ref{fig3}(b)
for a $\mathbb{P}_M$ on the M-W surface.
Analogously it can be done for a $\mathbb{P}_M$ on the cap.
\end{itemize}
All the above conditions yield to a non linear system of five equations to be solved
at each Gauss point of each finite element in the mesh.
In the next section the whole procedure is presented in a flow-chart style.

\subsection{Flow-chart of the iterative scheme of LMM} \quad

$\bullet$ \emph{Initialization}

\noindent Assign to all FEs an initial set of fictitious elastic parameters,
say $E^{(0)}$, $G^{(0)}$, $\nu^{(0)}$, and initial
stresses $\bar{\xi}^{(0)},\bar{\rho}_{x}^{(0)},\bar{\rho}_{y}^{(0)}$.
Set also: $k=1$,
$P_{_{UB}}^{(k-1)} = P_{_{UB}}^{(0)} =1$ (for $k=1$,
$P_{_{UB}}^{(0)}$ can be any arbitrary value) and compute
the bulk modulus $K^{(0)}=E^{(0)}$/ $3(1-2\nu^{(0)})$.

$\bullet$ \emph{Start iterations}

\noindent\textbf{Step 1:}
perform a fictitious elastic analysis with elastic
parameters $E^{(k-1)}$, $G^{(k-1)}$, $\nu^{(k-1)}$;
initial stresses $\bar{\xi}^{(k-1)},\bar{\rho}_{x}^{(k-1)},\bar{\rho}_{y}^{(k-1)}$;
loads $P_{_{UB}}^{(k-1)} \bar p_i$ (at the first iteration use the 
initialization quantities).
Compute a fictitious linear solution (at Gauss point level),
 namely: $ \dot{\varepsilon}_{v}^{{\ensuremath{\ell}}(k-1)}$,
$\dot{\varepsilon}_{d_{x}}^{{\ensuremath{\ell}}(k-1)}$,
$\dot{\varepsilon}_{d_{y}}^{{\ensuremath{\ell}}(k-1)}$,
$\dot{u}_{i}^{{\ensuremath{\ell}}(k-1)}$  together with the stress values
$\xi^{{\ensuremath{\ell}}(k-1)},\rho_{x}^{{\ensuremath{\ell}}(k-1)},
\rho_{y}^{{\ensuremath{\ell}}(k-1)}$.

\medskip
\noindent\textbf{Step 2:}
compute the constant value of the complementary
potential energy for later use:
$$
\bar{W}^{(k-1)} = {1 \over 2}  \xi^{{\ensuremath{\ell}}(k-1)}
                            \dot{\varepsilon}_{v}^{{\ell}(k-1)}+
                            \rho_{x}^{{\ell}(k-1)}
                            \dot{\varepsilon}_{d_{x}}^{{\ell}(k-1)}+
                            \rho_{y}^{{\ell}(k-1)}\
                            \dot{\varepsilon}_{d_{y}}^{{\ell}(k-1)}
$$

\noindent\textbf{Step 3:}
impose the matching conditions to update the fictitious quantities:
$\hat{\xi}_{\mathbb{C}}^{{(k)}}$, $\hat{\rho}_{\mathbb{C}}^{{(k)}}$,
$G^{{(k)}}$, $K^{{(k)}}$, $\Omega^{{(k)}}$ to be used, if necessary,
at next iteration (refer also to Figures \ref{fig3}(a)(b))
\begin{gather*}
 \frac{\hat{\xi}_M^{{(k-1)}}-\hat{\xi}_{\mathbb{C}}^{{(k)}}}{3K^{{(k)}}}
 =\dot{\varepsilon}_{v}^{\ensuremath{\ell}^{{(k-1)}}}\,, \\
 \frac{\hat{\rho}_M^{{(k-1)}}-\hat{\rho}_{\mathbb{C}}^{{(k)}}}{2G^{{(k)}}}
 =\dot{\varepsilon}_{d}^{\ensuremath{\ell}^{{(k-1)}}}\,,\\
 \frac{[\hat{\xi}_{\mathbb{P}_{i}}^{{(k-1)}}-\hat{\xi}_{\mathbb{C}}^{{(k)}}]^2}
{6K^{{(k)}}\Omega^{{(k)}}\bar{W}^{{(k-1)}}}
 +\frac{[\hat{\rho}_{\,\mathbb{P}_{i}}^{{(k-1)}}
 -\hat{\rho}_{\mathbb{C}}^{{(k)}}]^2}{4G^{{(k)}}\Omega^{{(k)}}\bar{W}^{{(k-1)}}}=1
  \quad \textrm{with }i=a,v,M.
\end{gather*}

\noindent\textbf{Step 4:}
evaluate the upper bound multiplier
$$
P_{_{UB}}^{(k)}={{\int_{\partial V_t}\bar{p}_{i}\,\dot{u}_{i}^{{c}{(k-1)}}\,\textrm{d} (\partial V)}\over
{\int_{V} \Big(\xi^{{Y}{(k-1)}}\,\dot{\varepsilon}_{v}^{{c}{(k-1)}}+ \,
\rho_{x}^{{Y}{(k-1)}}\, \dot{\varepsilon}_{d_{x}}^{{c}{(k-1)}}+
\rho_{y}^{{Y}{(k-1)}}\,\dot{\varepsilon}_{d_{y}}^{{c}{(k-1)}}\Big) } }
$$

\noindent\textbf{Step 5:} check for convergence
$$
 \vert P_{_{UB}}^{(k)} - P_{_{UB}}^{(k-1)} \vert \le
       \textrm{TOL} \quad
     \begin{cases}
          \textrm{YES} & \Rightarrow  \textrm{EXIT} \\
          \textrm{NOT} & \Rightarrow  \textrm{set $k = k - 1$ and GOTO step 1 }
\end{cases}
$$

$\bullet$ \emph{End iterations}


\section{Numerical application and concluding remarks}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig4}
\end{center}
\caption{Values of the upper bound $P_{_{UB}}$ to the peak load multiplier 
versus iteration numbers for a cubic  concrete sample subjected to a
 central line load} \label{fig4}
\end{figure}


The reliability of the discussed approach is verified by a simple numerical 
application of engineering interest concerning a cubic concrete sample 
of side $d=150\,$mm, with: $E=0.588\,$GPa; $\nu=0.2$.
Moreover, referring to Eqs.(2.1-3), the following material data have been assumed:
${f_c'}=30\,$MPa, ${f_t'}=3\,$MPa and ${e}=0.539$.
The sample is subjected to a central line reference load of global value 
equal to $100\,$kN.
The expected peak load is given by the formula 
$P_{_{U}}=\pi d^{2} f_{st} / 2 \,\bar{p}$,
after \cite{Roetal99}, with $f_{st}=f_t' / 0.9885$
so that, for the given data, the expected peak load multiplier is 
$P_{_{U}}=1.06$.
All the elastic analyses have been carried out using the FE code ADINA,
with meshes of 3D-Solid elements, while a Fortran main program has been 
developed
to drive the iterative procedure and the matching at each Gauss point.
Figure \ref{fig4} shows the numerical evaluated upper bound to the peak load multiplier,  
$P_{_{UB}}$, versus the iteration number. A monotonic and rapid convergence 
of the iterative  procedure is exhibited and the numerical $P_{_{UB}}$ value 
appears to be quite close to  the expected one.
Even if the considered example is simple, the obtained results seems to be 
quite encouraging, comparison with experimental findings on concrete prototypes 
of engineering interest exhibiting
a greater complexity are the object of an ongoing research.


\begin{thebibliography}{00}

\bibitem{JiBa02} M. Jir\'{a}sek, Z.P. Ba\v{z}ant;
\emph{Inelastic Analysis of Structures}, John Wiley \& Sons, LTD,
West Sussex, England, (2002).

\bibitem{LiCr10} T. Li, R. Crouch;
\emph{A $C_{2}$ plasticity model for structural concrete}.
Computers and Structures \textbf{88} (2010), 1322-1332.

\bibitem{MeWi95} P. Men\'etrey, K.J. Willam;
\emph{A triaxial failure criterion for concrete and its generalization}.
ACI Structural Journal \textbf{92}, 311-318.

\bibitem{PiFu07} A. A. Pisano, P. Fuschi;
 \emph{A numerical approach for limit analysis of orthotropic composite laminates}.
Int. J. Numerical Methods Eng \textbf{70} (2007), 71-93.

\bibitem{PiFu11} A. A. Pisano, P. Fuschi;
\emph{Mechanically fastened joints in composite laminates: Evaluation of 
load bearing capacity}.
Composites: Part B \textbf{42} (2011), 949-961.

\bibitem{PiFuDD12} A. A. Pisano, P. Fuschi, D. De Domenico;
\emph{A layered limit analysis of pinned-joints composite laminates:
Numerical versus experimental findings}.
Composites: Part B \textbf{43} (2012), 940-952.

\bibitem {PoCa97} A. R. S. Ponter, K. F. Carter;
\emph{Limit state solutions, based upon linear elastic solutions with
spatially varying elastic modulus}.
Comput Methods Appl Mech Eng \textbf{140} (1997), 237-258.

\bibitem{PoFuEn00} A. R. S. Ponter, P. Fuschi, M. Engelhardt;
\emph{Limit analysis for a general class of yield conditions}.
Eur. J. Mechanics/A Solids  \textbf{19} (2000), 401-421.

\bibitem{Ra61} D. Radenkovic;
\emph{Th\'eor\`emes limites pour un materiau de Coulomb \`a dilatation 
non standardis\'ee}.
CR Acad Sci Paris \textbf{252} (1961), 4103-4104.

\bibitem{Roetal99} C. Rocco, G.V. Guinea, J. Planas, M. Elices;
\emph{Size effect and boundary conditions in the Brazilian test: 
theoretical analysis}. Materials and Structures \textbf{32} (1999), 437-444.

\bibitem{WePo09}  D. Weichert and A. R. S. Ponter, Eds.;
 \emph{Limit States of Materials and Structures -  Direct Methods}, 
Springer Verlag, Wien (2009).

\end{thebibliography}

\end{document}