# Automated formulation and resolution of limit analysis problems

Jeremy Bleyer<EMAIL_ADDRESS>Ghazi Hassen<EMAIL_ADDRESS>Laboratoire Navier, Ecole des Ponts ParisTech, Univ Gustave Eiffel, CNRS
6-8 av Blaise Pascal, Cité Descartes, 77455 Champs-sur-Marne, FRANCE

###### Abstract

The present manuscript presents a framework for automating the formulation and
resolution of limit analysis problems in a very general manner. This framework
relies on FEniCS domain-specific language and the representation of material
strength criteria and their corresponding support function in the conic
programming setting. Various choices of finite element discretization,
including discontinuous Galerkin interpolations, are offered by FEniCS,
enabling to formulate lower bound equilibrium elements or upper bound elements
including discontinuities for instance. The numerical resolution of the
corresponding optimization problem is carried out by the interior-point solver
Mosek which takes advantage of the conic representation for yield criteria.
Through various illustrative examples ranging from classical continuum limit
analysis problems to generalized mechanical models such as plates, shells,
strain gradient or Cosserat continua, we show that limit analysis problems can
be formulated using only a few lines of code, discretized in a very simple
manner and solved extremely efficiently. This paper is accompanied by a FEniCS
toolbox implementing the above-mentioned framework.

###### keywords:

limit analysis , yield design , convex optimization , conic programming ,
FEniCS , generalized continua

##  1 Introduction

Limit analysis [26], or more generally yield design theory [53, 54], is an
efficient method for computing the ultimate load, or bearing capacity, of a
structure based on the sole knowledge of a given local strength criterion (or
plasticity criterion for limit analysis) and applied external loads. Its main
advantage is that it is a direct method i.e. it avoids the need of performing
a step-by-step non-linear analysis until complete failure, which is estimated
poorly in general since it is associated with non-convergence of the
resolution algorithm. More precisely, limit analysis theory can be formulated
as a convex optimization problem and therefore benefits from variational
approximations on the primal and associated dual problem yielding the so-
called lower bound static and upper bound kinematic approaches. The exact
collapse load can therefore be bracketed by the bounding status of the static
and kinematic solutions. This method has found tremendous applications in
mechanical and civil engineering problems since analytical upper bounds can be
obtained very efficiently by considering simple collapse mechanisms but also
because lower bounds provide a safe approximation to the exact collapse load.
Typical fields of application include soil slope stability, footing bearing
capacity or other geotechnical problems [17], rigid-block masonry structures,
design of reinforced-concrete structures [16, 49], especially through strut-
and-tie methods [57] or yield-line analysis [28], collapse loads of frame,
plate or shell structures [55, 56], etc.

In all the above-mentioned applications, hand-based solutions are quite easy
to compute and can be found in different design codes. However, limit analysis
techniques have been somehow limited to hand-based solutions for quite a long
time because of the difficulties encountered in the past when automating their
resolution in a finite-element discrete setting for instance. Indeed, the
corresponding optimization problems are inherently non-smooth and large-scale
which makes them challenging to solve. Standard gradient-based optimizers for
instance are not a good candidate for their resolution due to the highly non-
smooth aspect. The major breakthrough in this field is associated with the
development of efficient optimization algorithms particularly well-suited for
this class of optimization problems. The first progress has been made with the
simplex algorithm [19] for solving linear programming (LP) problems. However,
limit analysis problems fall into the LP category only when the strength
criterion is a polytope i.e. when it can be written as a collection of linear
inequalities [1, 58], whereas the majority of strength/plasticity criteria are
non-linear. The second progress has been achieved with the advent of so-called
interior-point (IP) solvers which improved the complexity of solving LP
problems [30], the simplex algorithm has exponential complexity compared to
polynomial (and, in practice, quasi-linear) complexity for interior-point
solvers. More importantly, the interior-point method has been extended to more
complex convex programs [48] such as second-order cone programs (SOCP) or
semi-definite programs (SDP). Combining the weak complexity of IP solvers and
the fact that most strength criteria can be expressed using second-order cone
or semi-definite constraints [3, 40] now enables to solve complex and large-
scale limit analysis problems [37, 43, 51, 61]. Obviously, a few other
alternative methods have been proposed but, to the authors’ opinion, none of
them have been shown to exhibit a performance similar to state-of-the-art IP
solvers.

As regards numerical discretization techniques for limit analysis problem, the
vast majority of works relies on the finite-element method [39, 1, 50]. The
specificity of limit analysis problems, compared to more standard nonlinear
computations with displacement-based FE discretizations, lies in the use of
static equilibrium-based finite-elements for obtaining true lower bound [58,
31] (and therefore safe) estimates of the limit load but also in the use of
discontinuous finite-elements for the kinematic upper bound approach [59, 42,
41]. Indeed, most hand-based upper bound solutions have been obtained
considering rigid-block mechanisms, thus involving no deformation but only
displacement jumps in the plastic dissipation computation. Despite the higher
computational cost compared to equivalent continuous interpolations,
discontinuous interpolations provide more accurate limit load estimates [50,
59, 36, 42], especially if finite-element edges are well-oriented. They are
also more robust for certain problems since they do not suffer from locking
issues, see for instance [46] for volumetric locking in pressure-insensitive
materials or [9] for shear-locking in thin plates.

Due to the specific nature of the optimization problems, formulating a
discretized version of a limit analysis (either static or kinematic) approach
requires forming matrices representing, for instance, equilibrium, continuity
or boundary conditions but also other linear relations coupling mechanical
variables (such as stress or strain) with auxiliary variables used to express
the strength conditions in a LP/SOCP or SDP format. Besides, depending on the
specific choice made for the optimization solver, the standard input format of
the problem may differ. As a result, discrete limit analysis problems require
access to matrices that are not readily available from standard displacement-
based FE solver and must therefore be implemented in an external program
before calling the optimizer. Combining this aspect with the various types of
FE discretizations and mechanical models makes the automation of limit
analysis problem a challenge. As a result, limit analysis codes are usually
limited to specific situations, sometimes with specific strength criteria.

In the present manuscript, we describe a general framework for the formulation
of limit analysis problems for different mechanical models (2D/3D continua,
plates/shells or generalized continua). Relying on the FEniCS finite-element
library and symbolic representation of operators and code-generation
capabilities, different FE discretization schemes (including discontinuous or
equilibrium elements) can be easily formulated and generalized to more
advanced mechanical models. The proposed framework therefore offers four
levels of generality in the problem formulation:

  * 1.

choice of a mechanical model: limit analysis problems possess the same
structure and can be formulated in a symbolic fashion through generalized
stress/strain definitions (section 2).

  * 2.

choice of a strength criterion: the formulation of its conic representation at
a local point suffices to completely characterize the strength criterion, its
translation to the global optimization problem being automatically performed
(section 3).

  * 3.

choice of a FE discretization: including element type, interpolation degree or
quadrature rule, all compatible for a variable number of degrees of freedom
related to the choice of the mechanical model (section 4).

  * 4.

choice of the optimization solver: although the accompanying FEniCS toolbox
relies extensively on the Mosek optimization solver [44], once formulated in a
standard conic programming form, the problem can then be written in a specific
file format appropriate for another solver.

This versatility is further illustrated by considering homogenization theory
in a limit analysis setting (section 5), plate and shell problems (section 6)
and generalized continua such as Cosserat or strain gradient models (section
7).

As regards numerical implementation, the present paper is accompanied by a
Python module fenics_optim.limit_analysis implemented as a submodule of the
FEniCS convex optimization package fenics_optim [6] described in [5].

##### Notations

$\boldsymbol{A}:\boldsymbol{B}=A_{ij}B_{ji}$,
$\boldsymbol{\mathcal{A}}\mathbin{{\mathchoice{\textstyle\vbox{\offinterlineskip\vskip-0.43057pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}}}{\textstyle\vbox{\offinterlineskip\vskip-0.43057pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}}}{\scriptstyle\scalebox{0.8}{$\vbox{\offinterlineskip\vskip-0.3014pt\hbox{$.$}\vskip-2.1097pt\hbox{$.$}\vskip-2.1097pt\hbox{$.$}}$}}{\scriptscriptstyle\scalebox{0.6}{$\vbox{\offinterlineskip\vskip-0.21529pt\hbox{$.$}\vskip-1.50693pt\hbox{$.$}\vskip-1.50693pt\hbox{$.$}}$}}}}\boldsymbol{\mathcal{B}}=\mathcal{A}_{ijk}\mathcal{B}_{kji}$.

##  2 A general framework for limit analysis problems

In this section, we consider a material domain $\Omega\subset\mathbb{R}^{d}$
(with $d=1,2,3$) associated with a specific mechanical model. In the
subsequent applications, we will consider classical continuum theories such as
2D or 3D Cauchy continua or Reissner-Mindlin plate models for instance but
also generalized continuum models encompassing higher-grade or higher-order
theories. For this reason, the subsequent presentation will be made in a
generalized continuum framework in which the mechanical stress or strain
measures, equilibrium or continuity equations and boundary conditions will be
written in an abstract fashion, their precise expression remaining to be
specified for each particular mechanical theory. In particular, the
presentation will be in the line of Germain’s construction through the virtual
work principle [23, 22].

###  2.1 Virtual work principle for generalized continua

Let us therefore consider a generalized virtual velocity field
$\boldsymbol{u}(\boldsymbol{x})$ of dimension $n$ and a set of strain measures
$\boldsymbol{D}\boldsymbol{u}$ of dimension $m$ with $\boldsymbol{D}$ being a
generalized strain operator. Following [23, 22], such strain measures must be
objective i.e. null for any rigid body motion. The virtual power of internal
forces is assumed to be given by an internal force density depending linearly
upon the strain measures:

$\mathcal{P}^{(i)}(\boldsymbol{u})=-\int_{\Omega}\boldsymbol{\Sigma}\cdot\boldsymbol{D}\boldsymbol{u}\mathop{}\\!\mathrm{dx}$
(1)

in which $\boldsymbol{\Sigma}$ denotes the generalized stress measure
associated with $\boldsymbol{D}\boldsymbol{u}$ by duality. The above
expression must in fact be understood in the sense of distributions i.e.
$\boldsymbol{u}$ may exhibit discontinuities $\boldsymbol{J}\boldsymbol{u}$
(consistent with the definition of operator $\boldsymbol{D}$) across some
internal surface $\Gamma$. The power of internal forces therefore writes more
explicitly as:

$\mathcal{P}^{(i)}(\boldsymbol{u})=-\left(\int_{\Omega\setminus\Gamma}\boldsymbol{\Sigma}\cdot\boldsymbol{D}\boldsymbol{u}\mathop{}\\!\mathrm{dx}+\int_{\Gamma}\boldsymbol{\Sigma}\cdot\boldsymbol{J}\boldsymbol{u}\mathop{}\\!\mathrm{dS}\right)$
(2)

The power of external forces is assumed to consist of two contributions: long-
range interaction forces described by a volume density $\boldsymbol{F}$ and
boundary contact forces described by a surface density $\boldsymbol{T}$ acting
on the exterior boundary $\partial\Omega$. Each power density depends linearly
upon the generalized velocity so that the total power is given by:

$\mathcal{P}^{(e)}(\boldsymbol{u})=\int_{\Omega}\boldsymbol{F}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dx}+\int_{\partial\Omega}\boldsymbol{T}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dS}$
(3)

According to the virtual power principle, the system is in equilibrium if and
only if the sum of the internal and external virtual powers is zero for any
virtual velocity field:

$\mathcal{P}^{(i)}(\boldsymbol{u})+\mathcal{P}^{(e)}(\boldsymbol{u})=0\quad\forall\boldsymbol{u}$
(equilibrium)

###  2.2 General formulation of a limit analysis problem

Limit analysis (or yield design) theory amounts to finding the maximum loading
a system can sustain considering only equilibrium and strength conditions for
its constitutive material. The latter can be generally described by the fact
that the generalized stresses $\boldsymbol{\Sigma}(\boldsymbol{x})$ must
belong to a strength domain $G(\boldsymbol{x})$ for all point
$\boldsymbol{x}\in\Omega$:

$\boldsymbol{\Sigma}(\boldsymbol{x})\in
G(\boldsymbol{x})\quad\forall\boldsymbol{x}\in\Omega$ (strength condition)

The strength domain $G\subseteq\mathbb{R}^{m}$ is assumed to be a convex set
(it may be unbounded and non-smooth) which usually contains the origin $0\in
G$.

Finding the maximum loading will be achieved with respect to a given loading
direction i.e. by assuming that both volume and surface forces depend upon a
single load factor $\lambda$ in an affine manner:

$\displaystyle\boldsymbol{F}(\lambda)$
$\displaystyle=\lambda\boldsymbol{f}+\boldsymbol{f}_{0}$ (4)
$\displaystyle\boldsymbol{T}(\lambda)$
$\displaystyle=\lambda\boldsymbol{t}+\boldsymbol{t}_{0}$ (5)

where $\boldsymbol{f},\boldsymbol{t}$ are given loading directions and
$\boldsymbol{f}_{0},\boldsymbol{t}_{0}$ are fixed reference loads. In such a
case, the power of external forces can also be written in an affine manner
with respect to $\lambda$:

$\mathcal{P}^{(e)}_{\lambda}(\boldsymbol{u})=\lambda\mathcal{P}^{(e)}(\boldsymbol{u})+\mathcal{P}^{(e)}_{0}(\boldsymbol{u})$
(6)

with obvious notations. Let us mention that if one wants to describe the set
of ultimate load defined by multiple loading parameters (see for instance
section 5), the loading directions must be varied in order to describe the
boundary of the ultimate load domain.

The limit analysis problem can finally be formulated as finding the maximum
load factor $\lambda$ such that there exists a generalized stress field
$\boldsymbol{\Sigma}(\boldsymbol{x})$ in equilibrium with
$(\boldsymbol{F}(\lambda),\boldsymbol{T}(\lambda))$ and complying with the
material strength properties i.e. satisfying both (equilibrium) and (strength
condition) which can also be written as:

$\begin{array}[]{rl}\displaystyle{\lambda^{+}=\sup_{\lambda,\boldsymbol{\Sigma}}}&\lambda\\\
\text{s.t.}&\mathcal{P}^{(i)}(\boldsymbol{u})+\lambda\mathcal{P}^{(e)}(\boldsymbol{u})+\mathcal{P}^{(e)}_{0}(\boldsymbol{u})=0\quad\forall\boldsymbol{u}\\\
&\boldsymbol{\Sigma}(\boldsymbol{x})\in
G(\boldsymbol{x})\quad\forall\boldsymbol{x}\in\Omega\end{array}$ (7)

Let us mention that for the infinite-dimensional convex problem (7) to have a
solution, the fixed loading $(\boldsymbol{f}_{0},\boldsymbol{t}_{0})$ must be
a sustainable loading i.e. there must exist a stress field in equilibrium with
$(\boldsymbol{f}_{0},\boldsymbol{t}_{0})$ and satisfying (strength condition).

Formulation (7) will be the basis of the mixed finite-element formulation
discussed in section 4.3 when choosing proper interpolation spaces for
$\boldsymbol{\Sigma}$ and $\boldsymbol{u}$. We now turn to the general
formulation of the static and kinematic approaches.

####  2.2.1 Static approach

Starting from the weak formulation of equilibrium given by (equilibrium),
strong balance equations, continuity conditions and boundary conditions can be
obtained for the generalized stresses $\boldsymbol{\Sigma}$. These will
generally take the following form:

$\displaystyle\boldsymbol{\mathcal{E}}\boldsymbol{\Sigma}+\lambda\boldsymbol{f}+\boldsymbol{f}_{0}=0$
$\displaystyle\quad\text{in }\Omega$ (8)
$\displaystyle\boldsymbol{\mathcal{C}}\boldsymbol{\Sigma}=0$
$\displaystyle\quad\text{on }\Gamma$ (9)
$\displaystyle\boldsymbol{\mathcal{T}}\boldsymbol{\Sigma}=\lambda\boldsymbol{t}+\boldsymbol{t}_{0}$
$\displaystyle\quad\text{on }\partial\Omega$ (10)

where $\boldsymbol{\mathcal{E}}$ is an equilibrium operator (adjoint to
$\boldsymbol{\mathcal{D}}$) and $\boldsymbol{\mathcal{C}}$ and
$\boldsymbol{\mathcal{T}}$ are some continuity and trace operators related to
$\boldsymbol{\mathcal{E}}$. A generalized stress field
$\boldsymbol{\Sigma}(\boldsymbol{x})$ satisfying these conditions will be
termed as statically admissible with a given loading
$(\lambda\boldsymbol{f}+\boldsymbol{f}_{0},\lambda\boldsymbol{t}+\boldsymbol{t}_{0})$.

The pure static formulation can therefore be generally written as:

$\begin{array}[]{rl}\displaystyle{\lambda^{+}=\sup_{\lambda,\boldsymbol{\Sigma}}}&\lambda\\\
\text{s.t.}&\boldsymbol{\mathcal{E}}\boldsymbol{\Sigma}+\lambda\boldsymbol{f}+\boldsymbol{f}_{0}=0\quad\text{in
}\Omega\\\ &\boldsymbol{\mathcal{C}}\boldsymbol{\Sigma}=0\quad\text{on
}\Gamma\\\
&\boldsymbol{\mathcal{T}}\boldsymbol{\Sigma}=\lambda\boldsymbol{t}+\boldsymbol{t}_{0}\quad\text{on
}\partial\Omega\\\ &\boldsymbol{\Sigma}(\boldsymbol{x})\in
G(\boldsymbol{x})\quad\forall\boldsymbol{x}\in\Omega\end{array}$ (SA)

Obviously $\boldsymbol{\Sigma}$ must belong to an appropriate functional space
$\mathcal{W}$ consistent with the nature of the above operators. If one
restricts to a (finite-dimensional) subset $\mathcal{W}_{h}\subset\mathcal{W}$
such that all constraints of (SA) can be satisfied exactly, the corresponding
solution $\lambda_{s}$ of the corresponding (finite) convex optimization
problem will therefore be a lower bound to the exact limit load:
$\lambda_{s}\leq\lambda^{+}$.

####  2.2.2 Kinematic approach

The kinematic formulation of a limit analysis problem can be obtained from an
equivalent formulation of (7):

$\lambda\leq\lambda^{+}\Longleftrightarrow\exists\boldsymbol{\Sigma}\text{
s.t.
}\begin{cases}\lambda\mathcal{P}^{(e)}(\boldsymbol{u})+\mathcal{P}^{(e)}_{0}(\boldsymbol{u})=-\mathcal{P}^{(i)}(\boldsymbol{u})\quad\forall\boldsymbol{u}\\\
\boldsymbol{\Sigma}(\boldsymbol{x})\in
G(\boldsymbol{x})\quad\forall\boldsymbol{x}\in\Omega\end{cases}$ (11)

One therefore has that:

$\lambda\mathcal{P}^{(e)}(\boldsymbol{u})+\mathcal{P}^{(e)}_{0}(\boldsymbol{u})\leq\sup_{\boldsymbol{\Sigma}(\boldsymbol{x})\in
G(\boldsymbol{x})}\\{-\mathcal{P}^{(i)}(\boldsymbol{u})\\}=\mathcal{P}^{(mr)}(\boldsymbol{u})\quad\forall\boldsymbol{u}$
(12)

where we introduced the maximum resisting power as:

$\mathcal{P}^{(mr)}(\boldsymbol{u})=\int_{\Omega\setminus\Gamma}\pi_{G}(\boldsymbol{D}\boldsymbol{u})\mathop{}\\!\mathrm{dx}+\int_{\Gamma}\pi_{G}(\boldsymbol{J}\boldsymbol{u})\mathop{}\\!\mathrm{dS}$
(13)

in which $\pi_{G}$ is the support function of the convex set $G$:

$\pi_{G}(\boldsymbol{d})=\sup_{\boldsymbol{\Sigma}\in
G}\\{\boldsymbol{\Sigma}\cdot\boldsymbol{d}\\}$ (14)

One can conclude that:

$\lambda\leq\lambda^{+}\Longrightarrow\lambda\leq\dfrac{\mathcal{P}^{(mr)}(\boldsymbol{u})-\mathcal{P}^{(e)}_{0}(\boldsymbol{u})}{\mathcal{P}^{(e)}(\boldsymbol{u})}\quad\forall\boldsymbol{u}$
(15)

Minimizing the right-hand side of the above relation therefore gives an upper
bound to the exact maximal load $\lambda^{+}$. Under appropriate mathematical
assumptions [47, 21] (non-restrictive in practice), it can be shown that the
minimum is in fact $\lambda^{+}$ so that one has:

$\lambda^{+}=\inf_{\boldsymbol{u}}\dfrac{\mathcal{P}^{(mr)}(\boldsymbol{u})-\mathcal{P}^{(e)}_{0}(\boldsymbol{u})}{\mathcal{P}^{(e)}(\boldsymbol{u})}$
(16)

Observing that the above quotient is invariant when rescaling $\boldsymbol{u}$
by a positive factor, a normalization constraint can in fact be considered to
remove the denominator so that, one finally has for the kinematic approach:

$\begin{array}[]{rl}\displaystyle{\lambda^{+}=\inf_{\boldsymbol{u}}}&\mathcal{P}^{(mr)}(\boldsymbol{u})-\mathcal{P}^{(e)}_{0}(\boldsymbol{u})\\\
\text{s.t.}&\mathcal{P}^{(e)}(\boldsymbol{u})=1\end{array}$ (KA)

Similarly to the static approach (SA), $\boldsymbol{u}$ must belong to an
appropriate functional space $\mathcal{V}$. If one restricts to a (finite-
dimensional) subset $\mathcal{V}_{h}\subset\mathcal{V}$, the corresponding
solution $\lambda_{u}$ of the corresponding (finite) convex optimization
problem will therefore be an upper bound to the exact limit load:
$\lambda^{+}\geq\lambda_{u}$.

##  3  Conic representation of strength criteria

This section deals with the representation of material strength criteria or
associated support functions in the framework of conic programming. The first
subsection describes the standard conic programming format used by the Mosek
solver although other solvers (e.g. CVXOPT, Sedumi, SDPT3) use a format which
is quite similar. We then discuss how the conic programming framework is used
for limit analysis problems.

###  3.1 Conic programming

Optimization problems entering the conic programming framework can be written
as:

$\begin{array}[]{rl}\displaystyle{\min_{\mathbf{x}}}&\mathbf{c}^{\text{T}}\mathbf{x}\\\
\text{s.t.}&\mathbf{b}_{l}\leq\mathbf{Ax}\leq\mathbf{b}_{u}\\\
&\mathbf{x}\in\mathcal{K}\end{array}$ (17)

where vector $\mathbf{c}$ defines a linear objective functional, matrix
$\mathbf{A}$ and vectors $\mathbf{b}_{u},\mathbf{b}_{l}$ define linear
inequality (or equality if $\mathbf{b}_{u}=\mathbf{b}_{l}$) constraints and
where
$\mathcal{K}=\mathcal{K}^{1}\times\mathcal{K}^{2}\times\ldots\times\mathcal{K}^{p}$
is a product of cones $\mathcal{K}^{i}\subset\mathbb{R}^{d_{i}}$ so that
$\mathbf{x}\in\mathcal{K}\Leftrightarrow\mathbf{x}^{i}\in\mathcal{K}^{i}\>\forall
i=1,\ldots,p$ where
$\mathbf{x}=(\mathbf{x}^{1},\mathbf{x}^{2},\ldots,\mathbf{x}^{p})$. These
cones can be of different kinds:

  * 1.

$\mathcal{K}^{i}=\mathbb{R}^{d_{i}}$ i.e. no constraint on $\mathbf{x}^{i}$

  * 2.

$\mathcal{K}^{i}=(\mathbb{R}^{+})^{d_{i}}$ is the positive orthant i.e.
$\mathbf{x}^{i}\geq 0$

  * 3.

$\mathcal{K}^{i}=\mathcal{Q}_{d_{i}}$ the quadratic Lorentz cone defined as:

$\mathcal{Q}_{d_{i}}=\\{\mathbf{z}\in\mathbb{R}^{d_{i}}\text{ s.t.
}\mathbf{z}=(z_{0},\bar{\mathbf{z}})\text{ and
}z_{0}\geq\|\bar{\mathbf{z}}\|_{2}\\}$ (18)

  * 4.

$\mathcal{K}^{i}=\mathcal{Q}_{d_{i}}^{r}$ the rotated quadratic Lorentz cone
defined as:

$\mathcal{Q}_{d_{i}}^{r}=\\{\mathbf{z}\in\mathbb{R}^{d_{i}}\text{ s.t.
}\mathbf{z}=(z_{0},z_{1},\bar{\mathbf{z}})\text{ and
}2z_{0}z_{1}\geq\|\bar{\mathbf{z}}\|_{2}^{2}\\}$ (19)

  * 5.

$\mathcal{K}^{i}=\mathbb{S}_{n_{i}}^{+}$ is the cone of semi-definite positive
matrices of dimension $n_{i}$.

If $\mathcal{K}$ contains only cones of the first two kinds, then the
resulting optimization problem (17) belongs to the class of Linear Programming
(LP) problems. If, in addition, $\mathcal{K}$ contains quadratic cones
$\mathcal{Q}_{d_{i}}$ or $\mathcal{Q}_{d_{i}}^{r}$, then the problem belongs
to the class of Second-Order Cone Programming (SOCP) problems. Finally, when
cones of the type $\mathbb{S}_{n_{i}}^{+}$ are present, the problem belongs to
the class of Semi-Definite Programming (SDP) problems.

###  3.2 Conic-representable functions and sets

In order to use conic programming solvers for our application, limit analysis
problems must therefore be formulated following format (17). Inspecting the
structure of (SA), it can be seen that all constraints and the objective
functions are linear, except for the strength condition (strength condition).
As a result, in order to fit format (17), only (strength condition) must be
reformulated in a conic sense. Similarly for problem (KA), only the support
function $\pi_{G}$ must be expressed in conic form to fit the standard format
(17). This reformulation step therefore depends on the specific choice of a
strength criterion. To do so, we define a generic form of conic-representable
functions and sets. The fenics_optim package [6] relies on this specific
notion for automating the formulation of generic convex problems. The
fenics_optim.limit_analysis module accompanying this paper uses these notions
and particularizes them for limit analysis problems.

Conic-representable functions are defined as the class of convex functions
which can be expressed as follows:

$\begin{array}[]{rl}\displaystyle{F(\mathbf{x})=\min_{\mathbf{y}}}&\mathbf{c}_{x}^{\text{T}}\mathbf{x}+\mathbf{c}_{y}^{\text{T}}\mathbf{y}\\\
\text{s.t.}&\mathbf{b}_{l}\leq\mathbf{Ax}+\mathbf{By}\leq\mathbf{b}_{u}\\\
&\mathbf{y}\in\mathcal{K}\end{array}$ (20)

with $\mathbf{x}\in\mathbb{R}^{n}$ and in which $\mathcal{K}$ is again a
product of cones of the kinds detailed in section 3.1. As a by-product of the
previous definition, conic-representable convex sets correspond to sets for
which the indicator function is conic-representable. If $\mathcal{K}$ contains
only second-order cones (SOC), then we will speak about a SOC-representable
function. A SDP-representable function corresponds to the case when
$\mathcal{K}$ contains SDP cones, whereas linear-representable functions
correspond to the case when the conic constraint $\mathbf{y}\in\mathcal{K}$ is
absent or contains only positive constraints $y_{i}\geq 0$.

It is easy to see that if $F$ is SOC-representable (resp. SDP-representable,
resp. linear-representable), then its Legendre-Fenchel conjugate $F^{*}$ is
also SOC-representable (resp. SDP-representable, resp. linear-representable).
For limit analysis applications, this means that if a strength criterion
(expressed as a convex set) is SOC-representable (resp. SDP-representable,
resp. linear-representable), then so will be its support function, and vice-
versa.

As already mentioned in [3, 34, 40], a large class of classical strength
criteria can be formulated in terms of second-order cone (SOC) constraints or
semi-definite positive (SDP) matrix constraints so that their expression or
their support function expressions can be expressed in the form (20). Below,
we give an example on how conic-representation is used for the case of a 2D
plane-strain Mohr-Coulomb material. Many generic conic reformulations can be
found in [38, 2, 12] and especially in the Mosek Modeling Cookbook [45].

###  3.3 Example of the plane-strain Mohr-Coulomb criterion

In plane-strain conditions, the Mohr-Coulomb criterion with cohesion $c$ and
friction angle $\phi$ writes as:

$\boldsymbol{\sigma}\in\text{MC}_{2D}(c,\phi)\Longleftrightarrow\sqrt{(\sigma_{xx}-\sigma_{yy})^{2}+4\sigma_{xy}^{2}}\leq
2c\cos\phi-(\sigma_{xx}+\sigma_{yy})\sin\phi$ (21)

which can be also written as:

$\displaystyle\sqrt{y_{1}^{2}+y_{2}^{2}}\leq
y_{0}\Leftrightarrow\mathbf{y}\in\mathcal{Q}_{3}$ (22)
$\displaystyle\text{with }\begin{Bmatrix}y_{0}\\\ y_{1}\\\ y_{2}\end{Bmatrix}$
$\displaystyle=\begin{bmatrix}-\sin\phi&-\sin\phi&0\\\ 1&-1&0\\\
0&0&2\end{bmatrix}\begin{Bmatrix}\sigma_{xx}\\\ \sigma_{yy}\\\
\sigma_{xy}\end{Bmatrix}+\begin{Bmatrix}2c\cos\phi\\\ 0\\\ 0\end{Bmatrix}$
(23)

This expression shows that the criterion is SOC-representable in the sense of
format (20).

Similarly, its support function is given by:

$\displaystyle\pi(\boldsymbol{d})$
$\displaystyle=\sup_{\boldsymbol{\sigma}\in\text{MC}_{2D}(c,\phi)}\sigma_{ij}d_{ij}$
(24)
$\displaystyle=\begin{cases}c\operatorname*{cotan}\phi\operatorname{tr}(\boldsymbol{d})&\text{if
}\operatorname{tr}(\boldsymbol{d})\geq\sin\phi\sqrt{(d_{xx}-d_{yy})^{2}+4d_{xy}^{2}}\\\
+\infty&\text{otherwise}\end{cases}$

which can be expressed as:

$\displaystyle\pi(\boldsymbol{d})=\min_{\mathbf{y}}\>\>$ $\displaystyle
c\operatorname*{cotan}\phi\operatorname{tr}(\boldsymbol{d})$ (25) s.t.
$\displaystyle\begin{bmatrix}1&1&0\\\ \sin\phi&-\sin\phi&0\\\
0&0&2\sin\phi\end{bmatrix}\begin{Bmatrix}d_{xx}\\\ d_{yy}\\\
d_{xy}\end{Bmatrix}=\mathbf{y}$ $\displaystyle\mathbf{y}\in\mathcal{Q}_{3}$

which also fits format (20). Note that the representation is by no means
unique since, for instance, we could have equivalently replaced the linear
objective term by $c\operatorname*{cotan}\phi\,y_{0}$ or we could also have
inverted the linear relation between $\boldsymbol{d}$ and $\mathbf{y}$ and
introduced $\widetilde{\mathbf{y}}=\mathbf{y}/\sin\phi$ so that:

$\displaystyle\pi(\boldsymbol{d})=\min_{\widetilde{\mathbf{y}}}\>\>$
$\displaystyle c\cos\phi\,\widetilde{y}_{0}$ (26) s.t.
$\displaystyle\begin{Bmatrix}d_{xx}\\\ d_{yy}\\\
d_{xy}\end{Bmatrix}=\begin{bmatrix}(\sin\phi)/2&1/2&0\\\
(\sin\phi)/2&-1/2&0\\\ 0&0&1/2\end{bmatrix}\widetilde{\mathbf{y}}$
$\displaystyle\widetilde{\mathbf{y}}\in\mathcal{Q}_{3}$

This expression has the advantage over (25) of being still well-defined when
$\phi=0$, enabling to recover the Tresca/von Mises case.

###  3.4 A gallery of conic-representable strength criteria

The fenics_optim.limit_analysis module provides access to a large class of
usual strength criteria through the conic formulation of the criterion
indicator and support functions, the latter being provided both in terms of
strains and velocity discontinuities, see section 4. Table 1 provides a list
of currently available material strength criteria.

Strength criterion | Mechanical model | Representation type
---|---|---
vonMises | 2D/3D | SOC
DruckerPrager | 2D/3D | SOC
Tresca2D | 2D | SOC
MohrCoulomb2D | 2D | SOC
Rankine2D | 2D | SOC
Tresca3D | 3D | SDP
MohrCoulomb3D | 3D | SDP
Rankine3D | 3D | SDP
TsaiWu | 3D | SOC
Table 1: List of available conic-representable strength criteria for 2D/3D
continua

##  4 Finite-element limit analysis of 2D and 3D continua

In this section, we consider the standard continuum model where
$\boldsymbol{\Sigma}=\boldsymbol{\sigma}$ is the classical Cauchy symmetric
stress tensor, $\boldsymbol{u}$ a 2D or 3D velocity field and the associated
strain is its symmetric gradient
$\boldsymbol{D}\boldsymbol{u}=\nabla^{s}\boldsymbol{u}$. For simplicity, we
will consider the case where there is no fixed loading
$\boldsymbol{f}_{0}=\boldsymbol{t}_{0}=0$.

###  4.1 Kinematic-based formulation

Let us now consider the finite-element discretization of the kinematic limit
analysis approach (KA) for a continuous velocity field and imposed tractions
on some part $\partial\Omega_{T}$ of the boundary:

$\begin{array}[]{rl}\lambda^{+}\leq\displaystyle{\inf_{\boldsymbol{u}\in\mathcal{V}_{h}}}&\displaystyle{\int_{\Omega}\pi_{G}(\nabla^{s}\boldsymbol{u})\mathop{}\\!\mathrm{dx}}\\\
\text{s.t.}&\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dx}+\int_{\partial\Omega_{T}}\boldsymbol{t}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dS}=1\end{array}$
(27)

where $\mathcal{V}_{h}$ is a finite-element subspace of Lagrange elements
based on a given mesh of typical mesh size $h$. In the above problem, the
computed objective function is an upper bound of the exact limit load factor
$\lambda^{+}$ only if the integral of the objective function term is evaluated
exactly. In general, this is not possible because of the non-linearity of
function $\pi_{G}$, except in the special case of $\boldsymbol{u}$ being
interpolated with $\mathbb{P}^{1}$-Lagrange elements so that the gradient is
cell-wise constant and the integral becomes trivial.

As explained in [41], keeping an exact upper-bound estimate of $\lambda^{+}$
requires this integral to be estimated by excess. This is possible for a mesh
consisting of simplex (straight edges) triangles and a
$\mathbb{P}^{2}$-Lagrange interpolation for $\boldsymbol{u}$ when using the
following so-called vertex quadrature scheme:

$\int_{T}F(\mathbf{r}(\boldsymbol{x}))\mathop{}\\!\mathrm{dx}\lesssim\dfrac{|T|}{d+1}\sum_{i=1}^{d+1}F(\mathbf{r}(\boldsymbol{x}_{i}))$
(28)

where $F$ is a convex function and $\mathbf{r}$ is an affinely-varying
function over the mesh cell $T$ (either a triangle in dimension $d=2$ or a
tetrahedron for $d=3$ of area/volume $|T|$) with
$\mathbf{r}(\boldsymbol{x}_{i})$ being its value at the $d+1$-vertices. In
this case, we have:

$\begin{array}[]{rl}\lambda^{+}\leq\lambda_{u}=\displaystyle{\inf_{\boldsymbol{u}\in\mathcal{V}_{h}}}&\displaystyle{\sum_{T\in\mathcal{T}_{h}}\dfrac{|T|}{d+1}\sum_{i=1}^{d+1}\pi_{G}(\nabla^{s}\boldsymbol{u}(\boldsymbol{x}_{i}))}\\\
\text{s.t.}&\displaystyle{\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dx}+\int_{\partial\Omega_{T}}\boldsymbol{t}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dS}=1}\end{array}$
(29)

As discussed in [5], the fenics_optim package enables to solve convex
variational problem of the form:

$\begin{array}[]{rl}\displaystyle{\inf_{\boldsymbol{u}\in\mathcal{V}}}&\displaystyle{\int_{\Omega}\left(j_{1}\circ\ell_{1}(\boldsymbol{u})+\ldots+j_{p}\circ\ell_{p}(\boldsymbol{u})\right)\mathop{}\\!\mathrm{dx}}\\\
\text{s.t.}&\boldsymbol{u}\in\mathcal{K}\end{array}$ (30)

where $j_{i}$ are conic-representable convex functions and $\ell_{i}$ are
linear operators which can be expressed using UFL symbolic operators. Each
individual term $j_{i}$ can be specified independently along with its conic
representation and the quadrature scheme used for the computation over
$\Omega$ so that (30) is in fact of the form:

$\begin{array}[]{rl}\displaystyle{\inf_{\boldsymbol{u}\in\mathcal{V}}}&\displaystyle{\sum_{i=1}^{p}\sum_{g_{i}=1}^{N_{g,i}}\omega_{g_{i}}j_{i}\left(\ell_{i}(\boldsymbol{u};\boldsymbol{x}_{g_{i}})\right)}\\\
\text{s.t.}&\boldsymbol{u}\in\mathcal{K}\end{array}$ (31)

where $\ell_{i}(\boldsymbol{u};\boldsymbol{x}_{g_{i}})$ denotes the evaluation
of $\ell_{i}(\boldsymbol{u})$ at a quadrature point $\boldsymbol{x}_{g_{i}}$.

Local auxiliary variables $\mathbf{y}$ of the conic representation (20) for
each $j_{i}$ will be added to the optimization problem for each quadrature
point $\boldsymbol{x}_{g_{i}}$. More details can be found in [5].

Obviously, problem (27) is a problem of the form (30) in which
$j_{1}=\pi_{G}$, $\ell_{1}=\nabla^{s}$ and $j_{2}$ is the indicator function
of the linear constraint
$\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dx}+\int_{\partial\Omega_{T}}\boldsymbol{t}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dS}=1$
with $\ell_{2}$ being the identity. Ignoring import statements and mesh
generation, we now give a few lines of Python script to illustrate how the
fenics_optim.limit_analysis module enables to formulate an upper bound limit
analysis problem for a 2D Mohr-Coulomb material using a
$\mathbb{P}^{2}$-Lagrange interpolation for
$\boldsymbol{u}\in\mathcal{V}_{h}$. First, the corresponding function space V
is defined and fixed boundary conditions are imposed on the part named
”border” of the boundary. A MosekProblem object is instantiated and a first
optimization field u belonging to function space V is added to the problem and
is constrained to satisfy the Dirichlet boundary conditions: V =
VectorFunctionSpace(mesh, ”CG”, 2) bc = DirichletBC(V, Constant((0.,0.)),
border)

prob = MosekProblem(”Upper bound limit analysis”) u = prob.add_var(V, bc=bc)
The external work normalization constraint is then added by defining the
function space for the Lagrange multiplier corresponding to the constraint
(here it is scalar so we use a "Real" function space) and passing the
corresponding constraint in its weak form as follows (here
$\boldsymbol{t}=0$): R = FunctionSpace(mesh, ”R”, 0) def Pext(lamb): return
[lamb*dot(f,u)*dx] prob.add_eq_constraint(R, A=Pext, b=1) Now, a Mohr-Coulomb
material is instantiated and provides access to its convex support function.
The input arguments are the strain $\nabla^{s}\boldsymbol{u}$ written in terms
of UFL operators as well as the choice for the quadrature scheme used for its
numerical evaluation. Here the vertex scheme (28) is chosen. This convex
function is then added to the problem before asking for its optimization by
Mosek. mat = MohrCoulomb2D(c, phi) strain = sym(grad(u)) pi =
mat.support_function(strain, quadrature_scheme=”vertex”)
prob.add_convex_term(pi)

prob.optimize()

###  4.2 Static-based formulation

Conversely, a lower bound estimate of the limit load factor $\lambda^{+}$ can
be obtained by considering a statically admissible discretization of the
stress field $\boldsymbol{\sigma}$. For this purpose, we consider the lower
bound element of [39, 58], which includes statically admissible
discontinuities between facets. The problem we aim at solving is the
following:

$\begin{array}[]{rl}\displaystyle{\lambda^{+}\geq\lambda_{s}=\sup_{\lambda\in\mathbb{R},\boldsymbol{\sigma}\in\mathcal{W}_{h}}}&\lambda\\\
\text{s.t.
}&\operatorname{div}\boldsymbol{\sigma}+\lambda\boldsymbol{f}=0\quad\text{in
}\Omega\\\
&[\\![\boldsymbol{\sigma}]\\!]\cdot\boldsymbol{n}=0\quad\text{through
}\Gamma\\\
&\boldsymbol{\sigma}\cdot\boldsymbol{n}=\lambda\boldsymbol{t}\quad\text{on
}\partial\Omega_{T}\\\ &\boldsymbol{\sigma}(\boldsymbol{x})\in
G\quad\forall\boldsymbol{x}\in\Omega\end{array}$ (32)

where $\Gamma$ is the set of internal facets and
$[\\![\boldsymbol{\sigma}]\\!]$ the stress discontinuity across this facet of
normal $\boldsymbol{n}$.

For this problem, we have as optimization variables one real $\lambda$ and a
discontinuous piecewise-linear field $\boldsymbol{\sigma}$ represented by a
vector of dimension 3 in 2D (6 in 3D), $\mathcal{W}_{h}$ being the associated
discontinuous $\mathbb{P}^{1}_{d}$ function space: prob = MosekProblem(”Lower
bound limit analysis”) R = FunctionSpace(mesh, ”R”, 0) W =
VectorFunctionSpace(mesh, ”DG”, 1, dim=3) lamb, Sig = prob.add_var([R, W]) sig
= as_matrix([[Sig[0], Sig[2]], [Sig[2], Sig[1]]])

Assuming piecewise constant body force over each cell, the first local
equilibrium equation can be equivalently written weakly using $\mathbb{P}^{0}$
velocity fields as Lagrange multipliers: V_eq = VectorFunctionSpace(mesh,
”DG”, 0) def equilibrium(u): return [dot(u,f)*lamb*dx, dot(u,div(sig))*dx]
prob.add_eq_constraint(V_eq, A=equilibrium) Note that in the equilibrium
function each block respectively corresponds to the optimization variable
$\lambda$ then $\boldsymbol{\sigma}$, in the order they have been initially
defined. See more details in [5] on the underlying block-structure of the
problem.

Similarly, discontinuous affine displacements $\boldsymbol{v}$ defined on the
mesh facets only111They are called Discontinuous Lagrange Trace elements in
FEniCS are used as Lagrange multipliers for the second and third constraints
(here $\boldsymbol{t}=0$): V_jump = VectorFunctionSpace(mesh, ”Discontinuous
Lagrange Trace”, 1) def continuity(v): return [None,
dot(avg(v),jump(sig,n))*dS \+ dot(v, dot(sig,n))*ds(0)]
prob.add_eq_constraint(V_jump, A=continuity) where the dS term corresponds to
the integral over all internal facets $\Gamma$ and ds(0) corresponds to the
integral over $\partial\Omega_{T}$ in this case.

Finally, the problem objective function is added to the problem as well as the
strength criterion condition. The latter is treated as a convex function
through its indicator function. A quadrature scheme is still needed since
quadrature points will correspond to points at which the strength condition
will be enforced. In the present case, the vertex scheme will enforce the
strength condition at the triangle vertices so that it will be satisfied
everywhere inside the cell by convexity. The maximization problem is then
solved by Mosek: prob.add_obj_func([1, None])

crit = mat.criterion(Sig, quadrature_scheme=”vertex”)
prob.add_convex_term(crit)

prob.optimize(sense=”maximize”)

Note that Mosek solutions give access to dual variables (Lagrange multipliers)
so that a pseudo velocity field can be obtained from u for instance.

###  4.3 Mixed finite-element discretizations

The kinematic formulation (27) is sometimes difficult to use because the
support function expression may be cumbersome to derive. As a result, a
static-based formulation as in (32) is usually more attractive as it requires
only the expression of the strength criterion itself. Formulation (7) enables
to produce upper bounds equivalent to (27), provided a proper choice of
interpolation spaces for $\boldsymbol{u}$ and $\boldsymbol{\Sigma}$ and
quadrature rules. In the general case, such a formulation forms the basis of
mixed finite-element interpolations for which the bounding status is lost in
general [14, 1, 15, 18]. This formulation reads here as:

$\begin{array}[]{rl}\displaystyle{\lambda_{m}=\sup_{\lambda\in\mathbb{R},\boldsymbol{\sigma}\in\mathcal{W}_{h}}}&\lambda\\\
\text{s.t.
}&\displaystyle{\int_{\Omega}\boldsymbol{\sigma}:\nabla^{s}\boldsymbol{u}\mathop{}\\!\mathrm{dx}=\lambda\left(\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dx}+\int_{\partial\Omega_{T}}\boldsymbol{t}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dS}\right)}\quad\forall\boldsymbol{u}\in\mathcal{V}_{h}\\\
&\boldsymbol{\sigma}(\boldsymbol{x})\in
G\quad\forall\boldsymbol{x}\in\Omega\end{array}$ (33)

in which the static equilibrium conditions have been replaced by their weak
counterpart using the virtual work principle for a class of kinematically
admissible continuous velocity fields $\boldsymbol{u}\in\mathcal{V}_{h}$. As
mentioned in [5], if $\mathcal{V}_{h}$ corresponds to a continuous
$\mathbb{P}^{1}$-Lagrange interpolation, $\mathcal{W}_{h}$ to a discontinuous
$\mathbb{P}^{0}$ interpolation of the stress field and the strength criterion
is enforced at one point in each cell (this is enough since
$\boldsymbol{\sigma}$ is cell-wise constant), then (33) is equivalent, it is
even the dual problem, to (27) for the same $\mathcal{V}_{h}$.

For other cases, such as $\mathcal{V}_{h}$ being $\mathbb{P}^{2}$-Lagrange and
$\mathcal{W}_{h}$ discontinuous $\mathbb{P}^{1}_{d}$-Lagrange, quadrature
rules must be specified both for the equilibrium constraint as well as for the
strength criterion (criterion enforcement points). If both quadrature rules
are identical, then problem (33) is equivalent to problem $\eqref{kinematic-
LA}$ for the chosen quadrature rule. Upon choosing a vertex scheme, the
objective value $\lambda_{m}$ will be an upper-bound to $\lambda^{+}$. Other
choices such as quadrature points located at the mid-side points do not
produce rigorous upper bounds but are usually observed to converge from above
in practice. When both quadrature rules are different, one obtains a true
mixed-interpolation and again the bounding status is lost.

The FEniCS formulation for the first case $\mathbb{P}^{1}/\mathbb{P}^{0}$ for
$\mathcal{V}_{h}/\mathcal{W}_{h}$ would read as (here quadrature rules are
trivial one-point rules by default): prob = MosekProblem(”Upper bound from
static formulation”) R = FunctionSpace(mesh, ”R”, 0) W =
VectorFunctionSpace(mesh, ”DG”, 0, dim=3) lamb, Sig = prob.add_var([R, W])

V = VectorFunctionSpace(mesh, ”CG”, 1) bc = DirichletBC(V, Constant((0, 0)),
border)

sig = as_matrix([[Sig[0], Sig[2]], [Sig[2], Sig[1]]]) def equilibrium(u):
return [lamb*dot(u, f)*dx, -inner(sig, sym(grad(u)))*dx]
prob.add_eq_constraint(V, A=equilibrium, bc=bc)

prob.add_obj_func([1, None])

crit = mat.criterion(Sig) prob.add_convex_term(crit)

prob.optimize(sense=”maximize”)

###  4.4 Discontinuous velocity fields

As mentioned in the introduction, the use of discontinuous velocity fields is
interesting in a limit analysis context due to a higher accuracy and the
absence of volumetric or shear locking effects.

In case of discontinuous velocity fields across a set $\Gamma$ of
discontinuity surfaces, the kinematic limit analysis formulation (27) now
becomes:

$\begin{array}[]{rl}\lambda^{+}\leq\lambda_{u}=\displaystyle{\inf_{\boldsymbol{u}\in\mathcal{V}_{h}}}&\displaystyle{\int_{\Omega}\pi_{G}(\nabla^{s}\boldsymbol{u})\mathop{}\\!\mathrm{dx}+\int_{\Gamma}\Pi_{G}([\\![\boldsymbol{u}]\\!];\boldsymbol{n})\mathop{}\\!\mathrm{dS}}\\\
\text{s.t.}&\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dx}+\int_{\partial\Omega_{T}}\boldsymbol{t}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dS}=1\end{array}$
(34)

where the second term of the objective function denotes the dissipated power
contribution of the velocity jumps $[\\![\boldsymbol{u}]\\!]$ through a
surface $\Gamma$ of normal $\boldsymbol{n}$. It is computed from the
discontinuity support function:

$\Pi_{G}(\boldsymbol{v};\boldsymbol{n})=\sup_{\boldsymbol{\sigma}\in
G}\\{(\boldsymbol{\sigma}\cdot\boldsymbol{n})\cdot\boldsymbol{v}\\}=\pi_{G}\left(\boldsymbol{v}\overset{s}{\otimes}\boldsymbol{n}\right)$
(35)

i.e. the jump operator is here
$\boldsymbol{J}\boldsymbol{u}=[\\![\boldsymbol{u}]\\!]\overset{s}{\otimes}\boldsymbol{n}$.
As a result, the second term will also be conic-representable and will be
treated similarly. Despite this relation, the conic representation of
$\Pi_{G}(\boldsymbol{v};\boldsymbol{n})$ is usually implemented explicitly in
the material library due to the potential savings in terms of auxiliary
variables compared to calling directly
$\pi_{G}\left(\boldsymbol{v}\overset{s}{\otimes}\boldsymbol{n}\right)$. A
local projection of $[\\![\boldsymbol{u}]\\!]$ on the facet
$(\boldsymbol{n},\boldsymbol{t})$-plane is performed as $\Pi_{G}$ functions
are in general naturally expressed in this local frame. Finally, adding the
discontinuity term to the optimization problem is similar to the first term,
one must just specify the integration measure (the inner facets) and the
quadrature rule which must be used to perform the facet integration.

###  4.5 Mesh refinement

For improving the quality of the computed limit load estimates, a mesh
refinement procedure is also implemented based on the contribution of each
cell to maximum resisting work $\mathcal{P}^{(mr)}$ in the context of a
kinematic approach. The cell contributions are sorted in descending order and
we compute the cumulated contribution to $\mathcal{P}^{(mr)}$ until reaching a
user-specified threshold $\eta\mathcal{P}^{(mr)}$ (with a ratio $\eta$ of 0.5
typically) of the total dissipation. The first $k$ cells whose total
contribution is at least $\eta\mathcal{P}^{(mr)}$ are then marked for mesh
refinement. In the case of discontinuous elements, the facet contribution to
$\mathcal{P}^{(mr)}$ is computed for each facet, split evenly between the two
sharing cells and added to the cell contributions. For the lower bound
approach, a similar procedure is used except that we use the dual variable
associated with the local equilibrium equation to reconstruct a piecewise
linear velocity field from which we compute a cell-wise contribution to the
total dissipation.

###  4.6 Vertical cut-off stability

The different discretization choices are illustrated on the stability of a
vertical slope under its self-weight $\gamma$. The rectangular domain, of
height $H$, is clamped on both bottom and left boundaries and traction-free on
the remaining boundaries (see Figure 1-right). The soil is modelled as a Mohr-
Coulomb material with cohesion $c$ and friction angle $\phi=30^{\circ}$ under
plane strain conditions. The slope factor of safety is given by the maximum
value of the non-dimensional quantity $(\gamma H/c)^{+}$. Convergence of he
factor of safety estimates for various finite-element discretizations are
reported on Figure 1-left as a function of the total number of elements, the
concentration of the local dissipation $\pi_{G}(\nabla^{s}\boldsymbol{u})$
along the slip-line is represented on Figure 1-right. The comparison between
continuous and discontinuous upper bound finite elements is reported in Figure
2. Similarly, the comparison between uniform and adaptive mesh refinement is
reported in Figure 3 along with the final adapted mesh.

Figure 1: Left: Convergence of the vertical slope factor of safety for various
finite element discretization: UB1 (resp. UB2) are the $\mathbb{P}^{1}$ (resp.
$\mathbb{P}^{2}$) Lagrange upper bound elements, LB the lower bound element
and Disp1, Disp2, Mixed correspond to non-upper bound elements considered in
[34]. Analytical upper bound from Chen [17]. Right: local dissipation map
Figure 2: Convergence of the vertical slope factor of safety for various
finite element discretization for continuous and discontinuous upper bound
elements. Analytical upper bound from Chen [17].

Figure 3: Left: Convergence of the vertical slope factor of safety for adapted
meshes. Analytical upper bound from Chen [17]. Right: adapted mesh after 7
refinement steps

##  5 Computation of macroscopic strength properties through homogenization
theory

In this section, we show how to adapt the above-described framework to compute
macroscopic strength properties of heterogeneous materials. Indeed, limit
analysis can be combined to homogenization theory [60, 13] to compute
macroscopic strength criteria of heterogeneous materials by solving an
auxiliary limit analysis problem formulated on a representative elementary
volume or a periodic unit-cell for periodically heterogeneous materials. We
will consider here the latter case and denote $\mathcal{A}$ the unit cell
domain. The strength conditions are described by a spatially varying strength
criterion $G(\boldsymbol{x})$. The main goal of homogenization in limit
analysis is to compute the macroscopic strength domain $G^{hom}$ defined as
follows:

$\boldsymbol{\Sigma}\in
G^{hom}\Longleftrightarrow\begin{cases}\operatorname{div}\boldsymbol{\sigma}=0&\text{in
}\mathcal{A}\\\
[\\![\boldsymbol{\sigma}]\\!]\cdot\boldsymbol{n}=0&\text{through }\Gamma\\\
\boldsymbol{\sigma}\cdot\boldsymbol{n}&\text{antiperiodic}\\\
\boldsymbol{\sigma}(\boldsymbol{x})\in
G(\boldsymbol{x})&\forall\boldsymbol{x}\in\mathcal{A}\\\
\displaystyle{\dfrac{1}{|\mathcal{A}|}\int_{\mathcal{A}}\boldsymbol{\sigma}\mathop{}\\!\mathrm{dx}=\boldsymbol{\Sigma}}\end{cases}$
(36)

Upon choosing a given loading direction $\boldsymbol{\Sigma}_{0}$ of arbitrary
magnitude, one can look for an estimate to the maximum load factor
$\lambda^{+}$ such that $\lambda^{+}\boldsymbol{\Sigma}_{0}\in G^{hom}$,
yielding the following maximization problem:

$\begin{array}[]{rll}\displaystyle{\lambda^{+}\geq\lambda_{s}=\sup_{\lambda\in\mathbb{R},\boldsymbol{\sigma}\in\mathcal{W}_{h}}}&\lambda&\\\
\text{s.t. }&\operatorname{div}\boldsymbol{\sigma}=0&\text{in }\mathcal{A}\\\
&[\\![\boldsymbol{\sigma}]\\!]\cdot\boldsymbol{n}=0&\text{through }\Gamma\\\
&\boldsymbol{\sigma}\cdot\boldsymbol{n}&\text{antiperiodic}\\\
&\boldsymbol{\sigma}(\boldsymbol{x})\in
G(\boldsymbol{x})&\forall\boldsymbol{x}\in\mathcal{A}\\\
&\displaystyle{\dfrac{1}{|\mathcal{A}|}\int_{\mathcal{A}}\boldsymbol{\sigma}\mathop{}\\!\mathrm{dx}=\lambda\boldsymbol{\Sigma}_{0}}\end{array}$
(37)

Definition (37) correspond to the static formulation of a limit analysis
problem formulated on $\mathcal{A}$ with loading being parametrized by
$\boldsymbol{\Sigma}_{0}$. Its dual counter-part (kinematic formulation) can
be shown to be given by222This formulation can obviously be easily extended to
take into account discontinuous velocity fields.:

$\begin{array}[]{rl}\lambda^{+}\leq\lambda_{u}=\displaystyle{\inf_{\boldsymbol{D}\in\mathbb{R}^{6},\boldsymbol{u}\in\mathcal{V}_{h}}}&\displaystyle{\int_{\Omega}\pi(\boldsymbol{D}+\nabla^{s}\boldsymbol{u};\boldsymbol{x})\mathop{}\\!\mathrm{dx}}\\\
\text{s.t.}&|\mathcal{A}|\boldsymbol{\Sigma}_{0}:\boldsymbol{D}=1\\\
&\boldsymbol{u}\text{ periodic}\end{array}$ (38)

Let us remark that:

$\Pi_{hom}(\boldsymbol{D})=\lambda^{+}\boldsymbol{\Sigma}_{0}:\boldsymbol{D}\leq\dfrac{1}{|\mathcal{A}|}\int_{\Omega}\pi(\boldsymbol{D}+\nabla^{s}\boldsymbol{u};\boldsymbol{x})\mathop{}\\!\mathrm{dx}$
(39)

where $\displaystyle{\Pi_{hom}(\boldsymbol{D}):=\sup_{\boldsymbol{\Sigma}\in
G^{hom}}\\{\boldsymbol{\Sigma}:\boldsymbol{D}\\}}$ is the macroscopic support
function.

In (38), the macroscopic strain is considered as an additional optimization
variable since the loading direction $\boldsymbol{\Sigma}_{0}$ is fixed. It is
also possible to prescribe directly the macroscopic strain direction
$\boldsymbol{D}$, leaving free the loading direction, by removing the first
constraint in (38).

Formulation (38) is applied to a 3D periodic porous medium with pores
following a face-centered cubic system (see Figure 4(a)) made of a Drucker-
Prager material $(c=1,\phi=30^{\circ})$. The unit cell response
$\boldsymbol{U}(\boldsymbol{x})=\boldsymbol{D}\cdot\boldsymbol{x}+\boldsymbol{u}(\boldsymbol{x})$
to a pure shear loading $\boldsymbol{\Sigma}_{xy,0}=1$ is represented on
Figure 4(b). The macroscopic strength domain in the
$(\boldsymbol{\Sigma}_{xx},\boldsymbol{\Sigma}_{yy})$ plane and in the
$(\boldsymbol{\Sigma}_{xx},\boldsymbol{\Sigma}_{xy})$ with other
$\boldsymbol{\Sigma}_{ij}=0$ have been represented on Figure 5. It can be
observed how much the presence of the pores reduces the original Drucker-
Prager criterion of the skeleton.

(a) Porous 3D medium unit cell

(b) Unit-cell response to a pure-shear loading

Figure 4: Homogenization of a 3D porous medium

(a) $(\boldsymbol{\Sigma}_{xx},\boldsymbol{\Sigma}_{yy})$-plane

(b) $(\boldsymbol{\Sigma}_{xx},\boldsymbol{\Sigma}_{xy})$-plane

Figure 5: Homogenized strength criterion

##  6 Plates and shells

###  6.1 Thin and thick plates

In this section, we show how limit analysis of thin and thick plates can be
tackled by following exactly the general format described in section 2 and
implemented easily by taking advantage of the high-level symbolic
representation of operators in FEniCS. First, let us recall that the kinematic
limit analysis of thin plates obeying a Love-Kirchhoff kinematics corresponds
to the following minimization problem:

$\begin{array}[]{rl}\displaystyle{\inf_{w}}&\displaystyle{\int_{\Omega}\pi_{G}(\nabla^{2}w)\mathop{}\\!\mathrm{dx}+\int_{\Gamma}\Pi_{G}([\\![\partial_{n}w]\\!]\boldsymbol{n};\boldsymbol{n})\mathop{}\\!\mathrm{dS}}\\\
\text{s.t.}&\displaystyle{\int_{\Omega}fw\mathop{}\\!\mathrm{dx}=1}\end{array}$
(40)

for a distributed loading $f$. It can be seen that (40) complies with (KA)
where the generalized velocity is only the virtual deflection $w$, the strain
operator is the curvature $\boldsymbol{D}\boldsymbol{u}=\nabla^{2}w$ and the
velocity jump is given by a normal rotation jump
$\boldsymbol{J}\boldsymbol{u}=[\\![\partial_{n}w]\\!]\boldsymbol{n}\otimes\boldsymbol{n}$.
In the above, $\pi_{G}$ and $\Pi_{G}$ are the support functions of the
corresponding thin plate strength criterion $G_{\text{bend}}$, expressed
solely on the bending moment tensor
$\boldsymbol{M}=\begin{bmatrix}M_{xx}&M_{xy}\\\ M_{xy}&M_{yy}\end{bmatrix}$
with
$\Pi_{G}([\\![\partial_{n}w]\\!]\boldsymbol{n};\boldsymbol{n})=\pi_{G}([\\![\partial_{n}w]\\!]\boldsymbol{n}\otimes\boldsymbol{n})$.
As a result, thin plate strength criteria and their support functions are
treated exactly as 2D continua with the only difference coming from the
definition of the strain and discontinuity operators. We refer to [5] for an
implementation example of thin plate limit analysis using the fenics_optim
package. Let us also remark that implementing the corresponding lower bound
static approach is more involved due to the more complicated continuity
conditions involving equivalent Kirchhoff shear forces (see [31, 7]).

As regards thick plates involving shear and bending strength conditions and
following Reissner-Mindlin kinematics, the main difference with respect to
2D/3D continua or thin plates is that two unknown fields must be considered
instead of one: bending moments $\boldsymbol{M}$ and shear forces
$\boldsymbol{Q}=(Q_{x},Q_{y})$ for the static approach and out-of-plane
deflection $w$ and plate rotations
$\boldsymbol{\theta}=(\theta_{x},\theta_{y})$ for the kinematic approach. We
now discuss only the latter since the former will be discussed in the more
general case of shells in the next subsection. Following [10], kinematic limit
analysis of thick plates can be written as:

$\begin{array}[]{rl}\displaystyle{\inf_{w,\boldsymbol{\theta}}}&\displaystyle{\int_{\Omega}\pi_{G}((\nabla\boldsymbol{\theta},\nabla
w-\boldsymbol{\theta}))\mathop{}\\!\mathrm{dx}+\int_{\Gamma}\Pi_{G}(([\\![\boldsymbol{\theta}]\\!],[\\![w]\\!]);\boldsymbol{n})\mathop{}\\!\mathrm{dS}}\\\
\text{s.t.}&\displaystyle{\int_{\Omega}fw\mathop{}\\!\mathrm{dx}=1}\end{array}$
(41)

Again, the structure is the same as before with
$\boldsymbol{u}=(w,\boldsymbol{\theta})$, the (generalized) strain
$\boldsymbol{D}\boldsymbol{u}=(\boldsymbol{\chi},\boldsymbol{\gamma})$
consisting of the curvature $\boldsymbol{\chi}=\nabla\boldsymbol{\theta}$ and
the shear strain $\boldsymbol{\gamma}=\nabla w-\boldsymbol{\theta}$. In the
above, the support functions $\pi_{G}$ and $\Pi_{G}$ are defined with respect
to a thick plate strength criterion involving both the bending moment tensor
$\boldsymbol{M}$ and the shear force vector $\boldsymbol{Q}$ i.e.
$\boldsymbol{\Sigma}=(\boldsymbol{M},\boldsymbol{Q})$. In [8, 10], different
choices of thick plate criteria are discussed, especially regarding the
bending/shear interaction. For the sake of simplicity, we consider here only
the case of no interaction between bending and shear, meaning that the thick
plate strength criterion is in fact decoupled between bending and shear,
taking the following form:

$(\boldsymbol{M},\boldsymbol{Q})\in G_{\text{thick
plate}}\Longleftrightarrow\begin{cases}\boldsymbol{M}\in G_{\text{bend}}\\\
\boldsymbol{Q}\in G_{\text{shear}}\end{cases}$ (42)

where typically $G_{\text{shear}}=\\{\boldsymbol{Q}\text{ s.t.
}\|\boldsymbol{Q}\|_{2}\leq Q_{0}\\}$ with $Q_{0}$ being the plate pure shear
strength. Since typical thin plate criteria are SOC-representable, so will be
the thick plate criterion (42). We also have the following expression for the
corresponding support function:

$\displaystyle\pi_{G_{\text{thick
plate}}}\left((\boldsymbol{\chi},\boldsymbol{\gamma})\right)$
$\displaystyle=\sup_{(\boldsymbol{M},\boldsymbol{Q})\in G_{\text{thick
plate}}}\\{\boldsymbol{M}:\boldsymbol{\chi}+\boldsymbol{Q}\cdot\boldsymbol{\gamma}\\}$
$\displaystyle=\sup_{\boldsymbol{M}\in
G_{\text{bend}}}\\{\boldsymbol{M}:\boldsymbol{\chi}\\}+\sup_{\boldsymbol{Q}\in
G_{\text{shear}}}\\{\boldsymbol{Q}\cdot\boldsymbol{\gamma}\\}$
$\displaystyle=\pi_{G_{\text{bend}}}(\boldsymbol{\chi})+\pi_{G_{\text{shear}}}(\boldsymbol{\gamma})$
(43)

with
$\pi_{G_{\text{shear}}}(\boldsymbol{\gamma})=Q_{0}\|\boldsymbol{\gamma}\|_{2}$
and $\pi_{G_{\text{bend}}}$ depending on the choice of the bending strength
criterion. Finally, the (generalized) velocity jump consisting of the rotation
and velocity jumps
$\boldsymbol{J}\boldsymbol{u}=([\\![\boldsymbol{\theta}]\\!]\overset{s}{\otimes}\boldsymbol{n},[\\![w]\\!]\boldsymbol{n})$
so that we also have:

$\Pi_{G}\left(([\\![\boldsymbol{\theta}]\\!],[\\![w]\\!]);\boldsymbol{n}\right)=\pi_{G}\left(\left([\\![\boldsymbol{\theta}]\\!]\overset{s}{\otimes}\boldsymbol{n},[\\![w]\\!]\boldsymbol{n}\right)\right)$
(44)

As discussed in [10], finite-element discretization for the upper bound limit
analysis of thick plates must be chosen with care. Indeed, continuous Lagrange
interpolations for both the deflection $w$ and the rotation field
$\boldsymbol{\theta}$ will exhibit shear locking in the thin plate limit.
Reference [10] considered fully discontinuous interpolation of both fields and
showed extremely good performances. To illustrate the versatility of the
proposed framework, we consider here a small variant, namely a continuous
$\mathbb{P}^{2}$-Lagrange interpolation for the deflection $w$ and a
discontinuous $\mathbb{P}^{1}_{d}$-Lagrange interpolation for the rotation
$\boldsymbol{\theta}$. This interpolation choice is expected to yield similar
performances, although slightly higher, as the fully discontinuous
$\mathbb{P}^{2}_{d}/\mathbb{P}^{1}_{d}$ interpolation of [10]. The problem of
a uniformly distributed (intensity $q$) clamped square plate of side length
$L$ and thickness $h$ is considered with a von Mises bending strength
criterion (bending strength $M_{0}=\sigma_{0}h^{2}/4$). The shear strength is
taken as $Q_{0}=\sigma_{0}h/\sqrt{3}$. As already observed in [8, 10], for the
case of no bending/shear interaction, the limit load for this problem is well
approximated by:

$q^{+}=\min\\{q_{\text{shear}}^{+};q_{\text{bending}}^{+}\\}$ (45)

where $q_{\text{shear}}^{+}=(2+\sqrt{\pi})Q_{0}/L$ is the pure shear solution
for a square plate (see the Cheeger set example of [5]) and
$q_{\text{bending}}^{+}\approx 44.2M_{0}/L^{2}$ is the pure bending thin plate
solution. This is indeed what is also observed for the present interpolation
with a 50$\times$50 mesh, see Figure 6. The correct pure shear (Figure 7) and
pure bending (Figure 8) mechanisms are also retrieved.

Figure 6: Normalized limit load $q^{+}L^{2}/M_{0}$ of a square clamped plate
as function of plate slenderness $L/h$: comparison with pure shear and pure
bending failure mechanisms

(a) Collapse mechanism

(b) Norm of the shear strain $\|\boldsymbol{\gamma}\|$

Figure 7: Pure shear collapse solution for a thick plate ($h/L=0.5$)

(a) Collapse mechanism

(b) Norm of the curvature $\|\boldsymbol{\chi}\|$

Figure 8: Pure bending collapse solution for a thin plate ($h/L=0.01$)

###  6.2 Shells and multilayered plates

Let us now briefly consider the case of shells, for which we will discuss here
a lower bound implementation only. We refer to [9] and references therein for
more details concerning limit analysis of shells, especially the upper bound
kinematic approach. The shell geometry will be approximated by a plane facet
discretization into triangular elements and will be described locally by a
unit normal $\boldsymbol{\nu}$ and a tangent plane spanned by two unit vectors
$\boldsymbol{a}_{1}$ and $\boldsymbol{a}_{2}$. This local frame is therefore
constant element-wise. The generalized internal forces of the shell model are
described by a symmetric membrane force tensor
$\boldsymbol{N}=N_{ij}\boldsymbol{a}_{i}\otimes\boldsymbol{a}_{j}$, a
symmetric bending moment tensor
$\boldsymbol{M}=M_{ij}\boldsymbol{a}_{i}\otimes\boldsymbol{a}_{j}$ and a shear
force vector $\boldsymbol{Q}=Q_{i}\boldsymbol{a}_{i}$ ($i,j=1,2$), the
components of which are expressed in the local tangent plane. We will consider
only thin shells, meaning that the shell strength criterion $G_{\text{shell}}$
is a convex set in the 6-dimensional $(\boldsymbol{N},\boldsymbol{M})$ space
(infinite shear strength assumption). Introducing
$\boldsymbol{T}=\boldsymbol{N}+\boldsymbol{\nu}\otimes\boldsymbol{Q}$, the
local equilibrium equations in a plane facet are given by:

$\displaystyle\operatorname{div}_{T}\boldsymbol{T}+\lambda\boldsymbol{f}$
$\displaystyle=0$ (46)
$\displaystyle\operatorname{div}_{T}\boldsymbol{M}+\boldsymbol{Q}$
$\displaystyle=0$ (47)

where $\operatorname{div}_{T}$ is the tangent plane divergence operator and
$\lambda\boldsymbol{f}$ a distributed loading with an amplification factor
$\lambda$. In addition to local equilibrium, continuity equations of the force
resultant $\boldsymbol{R}=\boldsymbol{T}\cdot\boldsymbol{n}$ and the normal
bending moment
$\boldsymbol{\mathcal{M}}=\boldsymbol{M}\cdot\boldsymbol{n}\times\boldsymbol{\nu}$
must be satisfied where $\boldsymbol{n}$ is the in-plane normal to a facet
edge:

$\displaystyle[\\![\boldsymbol{R}]\\!]$ $\displaystyle=0$ (48)
$\displaystyle[\\![\boldsymbol{\mathcal{M}}]\\!]$ $\displaystyle=0$ (49)

Finite-element discretization relies on a discontinuous $\mathbb{P}^{1}_{d}$
interpolation for both membrane forces $\boldsymbol{N}$ and shear forces
$\boldsymbol{Q}$ and a discontinuous $\mathbb{P}^{2}_{d}$ interpolation for
the bending moments $\boldsymbol{M}$. Considering cell-wise uniform
distributed loadings, local force equilibrium (46) is uniform and can
therefore be satisfied exactly using a cell-wise constant Lagrange multiplier
$\boldsymbol{u}$. Local moment equilibrium (47) can be satisfied exactly using
a cell-wise linear Lagrange multiplier $\boldsymbol{\theta}$. Force and moment
continuity equations (48)-(49) are respectively satisfied by considering
facet-wise linear and quadratic Lagrange multipliers $\boldsymbol{v}$ and
$\boldsymbol{\vartheta}$. We again give below the main part of the
corresponding Python code for formulating this rather complex problem:

prob = MosekProblem(”Shell lower bound limit analysis”) R =
FunctionSpace(mesh, ”R”, 0) Ne = VectorElement(”DG”, mesh.ufl_cell(), 1,
dim=3) Me = VectorElement(”DG”, mesh.ufl_cell(), 2, dim=3) Qe =
VectorElement(”DG”, mesh.ufl_cell(), 1, dim=2) W = FunctionSpace(mesh,
MixedElement([Ne, Me, Qe]))

lamb, Sig = prob.add_var([R, W]) prob.add_obj_func([1, None])

(N, M, Q) = split(Sig) T = as_matrix([[N[0], N[2]], [N[2], N[1]], [Q[0],
Q[1]]]) M = to_mat(M) where the different unknown fields
$\boldsymbol{\Sigma}=(\boldsymbol{N},\boldsymbol{M},\boldsymbol{Q})$ have been
defined and collected into a global vector. Equilibrium and continuity
equations are then defined as in section 4.2 for continua by specifying the
Lagrange multiplier function space and writing the weak form of the
constraint: V_f_eq = VectorFunctionSpace(mesh, ”DG”, 0, dim=3) def
force_equilibrium(u): u_loc = dot(Ploc, u) return [dot(u, f)*lamb*dx,
dot(u_loc, divT(T))*dx] prob.add_eq_constraint(V_f_eq, A=force_equilibrium)

V_m_eq = VectorFunctionSpace(mesh, ”DG”, 1, dim=2) def
moment_equilibrium(theta): return [None, dot(theta, divT(M)+Q)*dx]
prob.add_eq_constraint(V_m_eq, A=moment_equilibrium)

V_f_jump = VectorFunctionSpace(mesh, ”Discontinuous Lagrange Trace”, 1, dim=3)
Tglob = dot(Ploc.T, T) def force_continuity(v): return [None,
dot(avg(v),jump(Tglob, n_plan))*dS] prob.add_eq_constraint(V_f_jump,
A=force_continuity)

V_m_jump = VectorFunctionSpace(mesh, ”Discontinuous Lagrange Trace”, 2, dim=3)
Mglob = dot(Ploc_plane.T, M) def moment_continuity(vtheta): return [None,
dot(avg(vtheta), cross(jump(Mglob, n_plan), avg(nu)))*dS]
prob.add_eq_constraint(V_m_jump, A=moment_continuity) where the Ploc (resp.
Ploc_plane) variable corresponds to a rotation matrix transforming fields
expressed in the global
$(\boldsymbol{e}_{x},\boldsymbol{e}_{y},\boldsymbol{e}_{z})$ into the local
$(\boldsymbol{a}_{1},\boldsymbol{a}_{2},\boldsymbol{\nu})$ (resp.
$(\boldsymbol{a}_{1},\boldsymbol{a}_{2})$) frame. Note that we approximated
the bending moment continuity equation by
$[\\![\boldsymbol{M}\cdot\boldsymbol{n}]\\!]\times\widehat{\boldsymbol{\nu}}$
using an average normal vector
$\widehat{\boldsymbol{\nu}}=(\boldsymbol{\nu}^{+}+\boldsymbol{\nu}^{-})/2$.

If the previous snippets illustrate the efficiency of FEniCS high-level
symbolic formulations for this kind of complex problem, the conic
representation framework will also be extremely beneficial for formulating the
shell strength criterion. Indeed, as discussed in length in [9], even for the
simple case of a homogeneous von Mises thin shell, the strength condition
expressed in terms of $(\boldsymbol{N},\boldsymbol{M})$ stress-resultant
becomes extremely complicated [27] so that simple SOC-representable
approximations have been proposed in the past for the von Mises shell [52]. In
[9], we proposed a general way of formulating an
$(\boldsymbol{N},\boldsymbol{M})$ shell strength criterion for a general
multilayered shell through an up-scaling procedure. It is given by:

$(\boldsymbol{N},\boldsymbol{M})\in
G_{\text{shell}}\Longleftrightarrow\begin{cases}\exists\>\boldsymbol{\sigma}(z)\in\mathcal{G}_{\text{ps}}\quad\forall
z\in[-h/2;h/2]\text{ and s.t.}\\\
\displaystyle{\boldsymbol{N}=\int_{-h/2}^{h/2}\boldsymbol{\sigma}(z)\operatorname*{dz}}\\\
\displaystyle{\boldsymbol{M}=\int_{-h/2}^{h/2}(-z)\boldsymbol{\sigma}(z)\operatorname*{dz}}\end{cases}$
(50)

where $h$ is the shell thickness and $\mathcal{G}_{ps}$ is the material local
plane stress criterion, which may potentially depend on coordinate $z$ for a
multilayered shell. To make formulation (50) usable in practice, the local
plane stress distribution $\boldsymbol{\sigma}(z)$ is replaced by a discrete
set of plane stress states
$\boldsymbol{\sigma}_{g}=\boldsymbol{\sigma}(z_{g})$ expressed at quadrature
points $z_{g}$ which are used to approximate the two integrals:

$(\boldsymbol{N},\boldsymbol{M})\in
G_{\text{shell}}^{\text{approx}}\Longleftrightarrow\begin{cases}\exists\>\boldsymbol{\sigma}_{g}\in\mathcal{G}_{\text{ps}}\quad\forall
g=1,\ldots,n_{z}\text{ and s.t.}\\\
\displaystyle{\boldsymbol{N}=\sum_{g=1}^{n_{z}}\omega_{g}\boldsymbol{\sigma}_{g}}\\\
\displaystyle{\boldsymbol{M}=\sum_{g=1}^{n_{z}}(-z_{g})\omega_{g}\boldsymbol{\sigma}_{g}}\end{cases}$
(51)

where $\omega_{g}$ are the corresponding quadrature weights of the
$n_{z}$-points quadrature rule. As discussed in [9], the precise choice of the
quadrature rule leads to different kinds of approximations to
$G_{\text{shell}}$: e.g. an upper bound approximation is obtained with a
trapezoidal quadrature rule, a Gauss-Legendre quadrature leads to an
approximation with no lower or upper bound status, a rectangular rule will
give a lower bound approximation. In the following, we choose the latter to be
consistent with the lower bound status of the static approach.

From (51) it can readily be seen that if $\mathcal{G}_{\text{ps}}$ is SOC-
representable, so will be $G_{\text{shell}}^{\text{approx}}$. For instance, in
the case of a plane stress von Mises criterion of uniaxial strength
$\sigma_{0}$, it can be shown that:

$\boldsymbol{\sigma}\in\mathcal{G}_{\text{ps von
Mises}}\Leftrightarrow\begin{cases}\boldsymbol{\sigma}=\mathbf{J}\overline{\mathbf{y}}\\\
y_{0}=\sigma_{0}\\\ \|\overline{\mathbf{y}}\|\leq
y_{0}\end{cases}\quad\text{where }\mathbf{J}=\begin{bmatrix}1&1/\sqrt{3}&0\\\
0&2/\sqrt{3}&0\\\ 0&0&1/\sqrt{3}\end{bmatrix}$ (52)

in which the last constraint is a $\mathcal{Q}_{4}$ quadratic cone. Choosing a
rectangular quadrature rule with $\omega_{g}=h/n_{z}$ and
$z_{g}=-\frac{h}{2}+\frac{h}{n_{z}}(g-\frac{1}{2})$ for $g=1,\ldots,n_{z}$, we
therefore have:

$(\boldsymbol{N},\boldsymbol{M})\in
G_{\text{shell}}^{\text{approx}}\Longleftrightarrow\begin{cases}\exists\>\mathbf{y}_{g}=(y_{0g},\overline{\mathbf{y}}_{g})\in\mathcal{Q}_{4}\quad\forall
g=1,\ldots,n_{z}\text{ and s.t.}\\\ y_{0g}=\sigma_{0}\\\
\displaystyle{\boldsymbol{N}=\sum_{g=1}^{n_{z}}\frac{h}{n_{z}}\mathbf{J}\overline{\mathbf{y}}_{g}}\\\
\displaystyle{\boldsymbol{M}=\sum_{g=1}^{n_{z}}(-z_{g})\frac{h}{n_{z}}\mathbf{J}\overline{\mathbf{y}}_{g}}\end{cases}$
(53)

which obviously fits format (20). Let us finally remark that the approximation
will converge to the shell criterion $G_{\text{shell}}$ when increasing
$n_{z}$. In the following we took $n_{z}=6$.

(a) Limit load results: reference LB and UB limit loads from [9]

(b) Cylindrical shell pseudo-collapse mechanisms and normalized normal force
magnitude $\|\boldsymbol{N}\|/N_{0}$. Top: slenderness $2L/R=10$, bottom:
slenderness $2L/R=30$.

Figure 9: Cylindrical shell under self-weight

As an illustrative application, we consider the problem of a cylindrical shell
of length $2L$, radius $R$ and thickness $h=0.01R$, clamped at both
extremities and loaded by a self-weight uniform vertical loading
$\boldsymbol{f}=-q\boldsymbol{e}_{z}$ (see Figure 9(a)). The shape of the
collapse mechanism varies depending on the cylinder slenderness $2L/R$. For
sufficiently long cylinders, the computed limit load $q^{+}$ is well described
by the one obtained when representing the cylinder as a 1D beam
$q_{\text{beam}}^{+}=\frac{32}{\pi}N_{0}\left(\dfrac{R}{2L}\right)^{2}$ with
$N_{0}=\sigma_{0}h$ being the membrane uniaxial strength. The obtained limit
loads agree very well with the lower bound result of [9]. Having access to the
Lagrange multiplier $\boldsymbol{u}$, we reconstructed a pseudo-collapse
mechanism by performing a projection of $\boldsymbol{u}$ on a continuous
$\mathbb{P}^{1}$ space. The obtained deformations have been represented in
Figure 9(b) along with the normal force magnitude distribution. It can be seen
that the mechanisms agree well with those obtained from an upper bound
kinematic approach in [9], with a beam-like mechanism involving plastic hinges
at the clamped supports and mid-span for the case $2L/R=30$.

##  7 Generalized continua

In this last section, we further illustrate the proposed framework on two
generalized continuum models, namely a strain gradient and a Cosserat
continuum. We would like to point out that the numerical implementation of
limit analysis for these kinds of model is almost non-existent whereas we will
show that it can now be easily formulated with the proposed framework.

###  7.1 Strain-gradient material

In this subsection, we consider the extension of limit analysis to strain
gradient materials. We do not attempt at providing physical justifications for
using this kind of model but let us just mention that it can, for instance, be
obtained when considering the elastically rigid version of a strain-gradient
plasticity model. We therefore consider the following strain-gradient
generalization of the kinematic limit analysis theorem:

$\begin{array}[]{rl}\displaystyle{\inf_{\boldsymbol{u}\in\mathcal{V}_{h}}}&\displaystyle{\int_{\Omega}\pi_{G_{\text{SG}}}((\nabla^{s}\boldsymbol{u},\nabla^{2}\boldsymbol{u}))\mathop{}\\!\mathrm{dx}}\\\
\text{s.t.}&\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dx}+\int_{\partial\Omega_{T}}\boldsymbol{t}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dS}=1\end{array}$
(54)

where we considered only classical loadings (body forces or surface tractions)
and where the strain operator
$\boldsymbol{D}\boldsymbol{u}=(\nabla^{s}\boldsymbol{u},\nabla^{2}\boldsymbol{u})$
now includes both the first and second displacement gradient of
$\boldsymbol{u}$ with $(\nabla^{2}\boldsymbol{u})_{ijk}=u_{i,jk}$. We do not
consider here the equivalent static formulation but let us just point out that
the generalized stress
$\boldsymbol{\Sigma}=(\boldsymbol{\sigma},\boldsymbol{\tau})$ includes both
the standard Cauchy stress $\boldsymbol{\sigma}$ and the third-rank couple
stress tensor $\boldsymbol{\tau}=(\tau_{ijk})$ which is associated by duality
with the second gradient $\nabla^{2}\boldsymbol{u}$. The generalized strength
criterion therefore depends both on $\boldsymbol{\sigma}$ and
$\boldsymbol{\tau}$. For simplicity, we will consider the following extended
von Mises criterion:

$(\boldsymbol{\sigma},\boldsymbol{\tau})\in
G_{\text{SG}}\Leftrightarrow\sqrt{\frac{1}{2}(\boldsymbol{s}:\boldsymbol{s}+\ell^{-2}\boldsymbol{\tau}^{\text{T}}\mathbin{{\mathchoice{\textstyle\vbox{\offinterlineskip\vskip-0.43057pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}}}{\textstyle\vbox{\offinterlineskip\vskip-0.43057pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}}}{\scriptstyle\scalebox{0.8}{$\vbox{\offinterlineskip\vskip-0.3014pt\hbox{$.$}\vskip-2.1097pt\hbox{$.$}\vskip-2.1097pt\hbox{$.$}}$}}{\scriptscriptstyle\scalebox{0.6}{$\vbox{\offinterlineskip\vskip-0.21529pt\hbox{$.$}\vskip-1.50693pt\hbox{$.$}\vskip-1.50693pt\hbox{$.$}}$}}}}\boldsymbol{\tau})}\leq
k$ (55)

where $\boldsymbol{s}=\operatorname*{dev}\boldsymbol{\sigma}$,
$\boldsymbol{\tau}^{\text{T}}\mathbin{{\mathchoice{\textstyle\vbox{\offinterlineskip\vskip-0.43057pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}}}{\textstyle\vbox{\offinterlineskip\vskip-0.43057pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}}}{\scriptstyle\scalebox{0.8}{$\vbox{\offinterlineskip\vskip-0.3014pt\hbox{$.$}\vskip-2.1097pt\hbox{$.$}\vskip-2.1097pt\hbox{$.$}}$}}{\scriptscriptstyle\scalebox{0.6}{$\vbox{\offinterlineskip\vskip-0.21529pt\hbox{$.$}\vskip-1.50693pt\hbox{$.$}\vskip-1.50693pt\hbox{$.$}}$}}}}\boldsymbol{\tau}=\tau_{ijk}\tau_{ijk}$
and $\ell$ is an internal length scale. The associated support function is:

$\pi_{G_{\text{SG}}}((\boldsymbol{d},\boldsymbol{\eta}))=\sup_{(\boldsymbol{\sigma},\boldsymbol{\tau})\in
G_{\text{SG}}}\\{\boldsymbol{\sigma}:\boldsymbol{d}+\boldsymbol{\tau}^{\text{T}}\mathbin{{\mathchoice{\textstyle\vbox{\offinterlineskip\vskip-0.43057pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}}}{\textstyle\vbox{\offinterlineskip\vskip-0.43057pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}}}{\scriptstyle\scalebox{0.8}{$\vbox{\offinterlineskip\vskip-0.3014pt\hbox{$.$}\vskip-2.1097pt\hbox{$.$}\vskip-2.1097pt\hbox{$.$}}$}}{\scriptscriptstyle\scalebox{0.6}{$\vbox{\offinterlineskip\vskip-0.21529pt\hbox{$.$}\vskip-1.50693pt\hbox{$.$}\vskip-1.50693pt\hbox{$.$}}$}}}}\boldsymbol{\eta}\\}=\begin{cases}k\sqrt{2(\boldsymbol{d}:\boldsymbol{d}+\ell^{2}\boldsymbol{\eta}^{\text{T}}\mathbin{{\mathchoice{\textstyle\vbox{\offinterlineskip\vskip-0.43057pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}}}{\textstyle\vbox{\offinterlineskip\vskip-0.43057pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}}}{\scriptstyle\scalebox{0.8}{$\vbox{\offinterlineskip\vskip-0.3014pt\hbox{$.$}\vskip-2.1097pt\hbox{$.$}\vskip-2.1097pt\hbox{$.$}}$}}{\scriptscriptstyle\scalebox{0.6}{$\vbox{\offinterlineskip\vskip-0.21529pt\hbox{$.$}\vskip-1.50693pt\hbox{$.$}\vskip-1.50693pt\hbox{$.$}}$}}}}\boldsymbol{\eta})}&\text{if
}\operatorname{tr}\boldsymbol{d}=0\\\ +\infty&\text{otherwise}\end{cases}$
(56)

where $\boldsymbol{d}=\nabla^{s}\boldsymbol{u}$ and
$\boldsymbol{\eta}=\nabla^{2}\boldsymbol{u}$.

Restricting to a plane strain situation, one has
$\eta_{ij3}=\eta_{i3j}=\eta_{3ij}=0$ and:

$\displaystyle\eta_{111}=u_{1,11}$ $\displaystyle\qquad\eta_{211}=u_{2,11}$
$\displaystyle\eta_{122}=u_{1,22}$ $\displaystyle\qquad\eta_{222}=u_{2,22}$
(57) $\displaystyle\eta_{112}=\eta_{121}=u_{1,12}$
$\displaystyle\qquad\eta_{212}=\eta_{221}=u_{2,12}$

Introducing
$\mathbf{D}=(d_{11},d_{22},\sqrt{2}d_{12},\ell\eta_{111},\ell\eta_{122},\sqrt{2}\ell\eta_{112},\ell\eta_{211},\ell\eta_{222},\sqrt{2}\ell\eta_{212})$,
one has
$\mathbf{D}^{\text{T}}\mathbf{D}=\boldsymbol{d}:\boldsymbol{d}+\ell^{2}\boldsymbol{\eta}^{\text{T}}\mathbin{{\mathchoice{\textstyle\vbox{\offinterlineskip\vskip-0.43057pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}}}{\textstyle\vbox{\offinterlineskip\vskip-0.43057pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}\vskip-3.01385pt\hbox{$.$}}}{\scriptstyle\scalebox{0.8}{$\vbox{\offinterlineskip\vskip-0.3014pt\hbox{$.$}\vskip-2.1097pt\hbox{$.$}\vskip-2.1097pt\hbox{$.$}}$}}{\scriptscriptstyle\scalebox{0.6}{$\vbox{\offinterlineskip\vskip-0.21529pt\hbox{$.$}\vskip-1.50693pt\hbox{$.$}\vskip-1.50693pt\hbox{$.$}}$}}}}\boldsymbol{\eta}$
so that
$\pi_{G_{\text{SG}}}((\boldsymbol{d},\boldsymbol{\eta}))=k\sqrt{2}\|\mathbf{D}\|_{2}$.
Since $\pi_{G_{\text{SG}}}$ involves a $L_{2}$-norm on a 9-dimensional vector,
it can be represented using a 10-dimensional quadratic cone $\mathcal{Q}_{10}$
(see [5]).

As regards the finite-element discretization, we choose a
$\mathbb{P}^{2}$-Lagrange interpolation for $\boldsymbol{u}$. In formulation
(54), it is implicitly assumed that both $\boldsymbol{u}$ and
$\nabla\boldsymbol{u}$ are continuous. The latter condition will not be
achieved easily by a standard FE discretization so that we supplement (54) by
a discontinuity term for
$\partial_{n}\boldsymbol{u}=\nabla\boldsymbol{u}\cdot\boldsymbol{n}$,
similarly to thin plates:

$\begin{array}[]{rl}\displaystyle{\inf_{\boldsymbol{u}\in\mathcal{V}_{h}}}&\displaystyle{\int_{\Omega}\pi_{G_{\text{SG}}}((\nabla^{s}\boldsymbol{u},\nabla^{2}\boldsymbol{u}))\mathop{}\\!\mathrm{dx}+\int_{\Gamma}\pi_{G_{\text{SG}}}((0,[\\![\partial_{n}\boldsymbol{u}]\\!]\otimes\boldsymbol{n}\otimes\boldsymbol{n}))\mathop{}\\!\mathrm{dS}}\\\
\text{s.t.}&\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dx}+\int_{\partial\Omega_{T}}\boldsymbol{t}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dS}=1\end{array}$
(58)

with
$\pi_{G_{\text{SG}}}((0,\partial_{n}\boldsymbol{u}\otimes\boldsymbol{n}\otimes\boldsymbol{n}))=k\sqrt{2}\ell\|[\\![\partial_{n}\boldsymbol{u}]\\!]\|_{2}$
for (56). Again, this can be easily implemented in very few lines of code,
regarding the problem complexity: prob = MosekProblem(”Strain gradient limit
analysis”) u = prob.add_var(V, bc=bc)

D = as_vector([u[0].dx(0), u[1].dx(1), (u[0].dx(1)+u[1].dx(0))/sqrt(2),
l*u[0].dx(0).dx(0), l*u[0].dx(1).dx(1), sqrt(2)*l*u[0].dx(0).dx(1),
l*u[1].dx(0).dx(0), l*u[1].dx(1).dx(1), sqrt(2)*l*u[1].dx(0).dx(1)]) pi =
L2Norm(D, quadrature_scheme=”vertex”) prob.add_convex_term(k*sqrt(2)*pi)

isochoric = EqualityConstraint(div(u), quadrature_scheme=”vertex”)
prob.add_convex_term(isochoric)

n = FacetNormal(mesh) pi_d = L2Norm([jump(k*sqrt(2)*l*grad(u), n)],
on_facet=True) prob.add_convex_term(pi_d)

prob.optimize()

As an illustrative application, we consider a rectangular domain of dimensions
$L\times 1.5L$ perforated at its center by a circular hole of radius $R=0.2L$.
The bottom boundary is fully clamped and the top one is displaced vertically
$\boldsymbol{u}=(0,U)$. No other loading is applied and the computed objective
value of (58) will be $Q^{+}U$ where
$Q=\int_{y=H}\sigma_{yy}\mathop{}\\!\mathrm{dS}$ is the resultant force on the
top boundary. The evolution of the normalized uniaxial strength $Q^{+}/(kL)$
is plotted for various values of the internal length parameter $\ell/L$ in
Figure 10 for two different mesh sizes. As expected, the plate apparent
strength is size-dependent and exhibits a strengthening behaviour for larger
values of $\ell$ or, equivalently, smaller sample size $L$. The standard
continuum limit analysis results (dashed lines) are retrieved when $\ell/L\to
0$. Collapse mechanisms along with dissipation fields $\pi_{\text{SG}}$ are
represented in Figure 11. Broadening of the plastic dissipation slip zones can
clearly be observed for increasing values of $\ell/L$.

Figure 10: Evolution of the normalized uniaxial strength for the strain-
gradient perforated traction as function of the internal length scale $\ell/L$
for two mesh sizes (dashed lines correspond to standard continuum limit
analysis results $\ell=0$).

(a) $\ell/L=0$ (standard continuum)

(b) $\ell/L=0.001$

(c) $\ell/L=0.01$

(d) $\ell/L=0.1$

Figure 11: Collapse mechanism and plastic dissipation as function of internal
length scale ratio $\ell/L$

###  7.2 A Cosserat-continuum model for jointed rocks

We further illustrate the ability of the proposed framework to tackle
generalized continua by considering a Cosserat (or micropolar continuum) model
for jointed rocks, initially proposed in [20]. The governing equations of the
model, in plane strain conditions, involve a non-symmetric stress tensor
$\boldsymbol{\Sigma}=\begin{bmatrix}\Sigma_{11}&\Sigma_{12}\\\
\Sigma_{21}&\Sigma_{22}\end{bmatrix}$ and a couple stress vector
$\boldsymbol{H}=(H_{1},H_{2})$ both expressed in the local reference frame
$(\boldsymbol{e}_{1},\boldsymbol{e}_{2})$ of the jointed rock mass (see Figure
12). The corresponding equilibrium equations read as:

$\displaystyle\operatorname{div}\boldsymbol{\Sigma}+\boldsymbol{f}$
$\displaystyle=0$ (59)
$\displaystyle\operatorname{div}\boldsymbol{H}+\Sigma_{21}-\Sigma_{12}$
$\displaystyle=0$ (60)

the corresponding weak form obtained from the virtual work principle being:

$\int_{\Omega}\left(\boldsymbol{\Sigma}^{\text{T}}:(\nabla\boldsymbol{u}-\operatorname{skew}\omega)+\boldsymbol{H}\cdot\nabla\omega\right)\mathop{}\\!\mathrm{dx}=\int_{\Omega}\boldsymbol{f}\cdot\boldsymbol{u}\mathop{}\\!\mathrm{dx}$
(61)

for any continuous test function $\boldsymbol{u}$ and $\omega$ with
$\operatorname{skew}\omega=\omega(\boldsymbol{e}_{2}\otimes\boldsymbol{e}_{1}-\boldsymbol{e}_{1}\otimes\boldsymbol{e}_{2})$
and where we considered only body forces as loading parameters.

As regards strength properties, the rock mass is assumed to obey a Mohr-
Coulomb criterion of cohesion $c_{m}$ and friction angle $\phi_{m}$. The
joints are represented as an orthogonal array, spaced by a length $\ell$ and
making an angle $\theta$ with the horizontal axis. They are assumed to also
obey a Mohr-Coulomb condition with parameters $(c_{j},\phi_{j})$. The
generalized strength condition for a jointed rock mass modelled as a Cosserat
continuum is expressed as [20]:

$(\boldsymbol{\Sigma},\boldsymbol{H})\in
G_{\text{Cosserat}}\Longleftrightarrow\begin{cases}\Sigma_{11}\tan\phi_{j}+|\Sigma_{21}|\leq
c_{j}\\\ \Sigma_{22}\tan\phi_{j}+|\Sigma_{12}|\leq c_{j}\\\
\dfrac{\ell}{2}\Sigma_{11}+|H_{1}|\leq\dfrac{\ell}{2}\dfrac{c_{j}}{\tan\phi_{j}}\\\
\dfrac{\ell}{2}\Sigma_{22}+|H_{2}|\leq\dfrac{\ell}{2}\dfrac{c_{j}}{\tan\phi_{j}}\\\
|H_{1}|\leq\dfrac{\ell}{2}\dfrac{c_{j}}{\cos\phi_{j}}\\\
|H_{2}|\leq\dfrac{\ell}{2}\dfrac{c_{j}}{\cos\phi_{j}}\\\
\operatorname*{sym}\boldsymbol{\Sigma}\in
G_{\text{MC,2D}}(c_{m},\phi_{m})\end{cases}$ (62)

where the last condition expresses the rock mass Mohr-Coulomb criterion on
$\operatorname*{sym}\boldsymbol{\Sigma}=(\boldsymbol{\Sigma}+\boldsymbol{\Sigma}^{\text{T}})/2$
and where all other conditions involve the joints resistance. Let us point out
that the case $\ell=0$ induces $H_{i}=0$ and thus
$\boldsymbol{\Sigma}=\boldsymbol{\Sigma}^{\text{T}}$ due to (60), one
therefore retrieves a Cauchy model with a strength criterion described by the
first, second and last conditions of (62). Finally, $G_{\text{Cosserat}}$
involves only linear inequality constraints in addition to the Mohr-Coulomb
criterion $G_{\text{MC,2D}}$. It is, therefore, SOC-representable, the part
involving joints only being linear-representable.

A mixed approach for this model has been implemented in the spirit of (33)
which avoids the need to compute the support function expression associated
with (62). Continuous $\mathbb{P}^{2}$ (resp. $\mathbb{P}^{1}$) Lagrange
elements have been used for $\boldsymbol{u}$ (resp. $\omega$) and
discontinuous $\mathbb{P}^{1}_{d}$-Lagrange elements for both
$\boldsymbol{\Sigma}$ and $\boldsymbol{H}$. The strength conditions have been
imposed at the vertices of each element. We considered the stability problem
of an excavation of height $H$, making a $25^{\circ}$ angle with the vertical
and subjected to its self-weight of intensity $\gamma$. The problem amounts to
find the maximum value of the non-dimensional stability factor
$K^{+}=\left(\dfrac{\gamma H}{c_{m}}\right)^{+}$. For numerical applications,
we took $c_{j}=0.5c_{m}$, $\phi_{j}=20^{\circ}$, $\phi_{m}=40^{\circ}$,
$\theta=10^{\circ}$ and varied the joint spacing $\ell$. The evolution of the
stability factor estimates as a function of $\ell/H$ has been represented in
Figure 13 for two different mesh sizes. As for the strain gradient model,
strengthening is observed for increasing $\ell/H$ ratios. Interestingly, size-
effects are much stronger for this problem than those of Figure 10. The
obtained value in the standard Cauchy ($\ell=0$) case is quite close to the
analytical upper bound of $K^{+}\leq 1.47$ derived for the same problem in
[20]. Finally, collapse mechanisms and a measure of the pure Cosserat
contribution
$(\Sigma_{21}-\Sigma_{12})(u_{2,1}-u_{1,2}-\omega)+\boldsymbol{H}\cdot\nabla\omega$
to the total dissipation have been represented in Figure 14. It can be
observed that the shape of the collapse mechanism and the location of
”shearing” zones involving Cosserat effects is quite dependent on the joint
spacing. For $\ell=0$, a triangular sliding block with a concentrated slip
zone is obtained, approximately corresponding to the merging of the two slip
bands of Figure 14(a).

Figure 12: Stability of a jointed rock excavation Figure 13: Evolution of the
stability factor as function of the joints spacing $\ell/H$ for two mesh sizes
(dashed lines correspond to standard continuum limit analysis results
$\ell=0$).

(a) $\ell/H=0.001$

(b) $\ell/H=0.01$

(c) $\ell/H=0.1$

Figure 14: Collapse mechanism of the jointed rock excavation and pure Cosserat
contribution to the total dissipation

##  8 Conclusions

This paper proposed a way to easily formulate and solve limit analysis
problems by taking advantage of three essential ingredients:

  * 1.

the FEniCS finite element library with its high-level domain specific language
and large choice of finite-element interpolations

  * 2.

the representation of limit analysis criteria and associated support functions
in a conic programming format

  * 3.

the resolution of the corresponding conic optimization problems by a dedicated
and efficient interior-point solver implemented in the Mosek software package

The first two ingredients offer an extremely large versatility in the problem
definition, giving access to an extremely large range of applications. In
particular, various finite-element interpolations can be easily defined,
giving access to traditional continuous displacement-based upper-bound
formulations but also their discontinuous counterpart through Discontinuous
Galerkin function spaces. Equilibrium lower-bound elements are therefore also
accessible through such spaces as well as mixed formulations through different
choices of quadrature rules. We also showed that plates and shells problems
could also be discretized without effort. The conic representation format (20)
offers a unified way of defining strength conditions and associated support
functions for different mechanical models, ranging from classical 2D/3D
continuum mechanics to plate bending criteria including potential shear
conditions, shell criteria with membrane/bending interaction or even
generalized continua such as Cosserat or strain gradient models. If the conic
representation format is large enough to encompass many strength criteria, it
is also sufficiently disciplined to yield optimization problems of the conic
programming class for which dedicated solvers like Mosek have been designed.
Mosek is indeed known to be a state-of-the-art optimizer for this class of
problems and therefore offers efficiency and robustness of the solution
procedure.

Obviously, the present work could still improve upon some aspects, in
particular regarding computational efficiency. For instance, many additional
auxiliary variables are usually introduced when complying with the conic
programming format. Some of them may be handled and eliminated by Mosek during
its pre-processing phase, although this is not entirely clear since it is used
as a black-box. Devising an interior-point solver specific to limit analysis
problems can take advantage of the problem structure and may be more
efficient. Besides, there are two points which prevent solving extremely
large-scale 3D problems. The first one is related to the use of direct solvers
in the interior-point inner iterations which requires large memory capacities
for 3D problems. The use of iterative solvers is still an active research
topic due to the difficulty of efficiently preconditioning the interior-point
linear systems. The second one concerns the need to solve SDP problems when
SDP-representable criteria like Mohr-Coulomb, Tresca or Rankine are used in
3D. Even though interior-point solvers efficiency has greatly improved for SDP
problems over the last decade, it is still more difficult to solve than an
SOCP problem of similar size. Improving even more their efficiency or finding
alternate strategies would be a great benefit for such 3D problems.

If mathematical programming (and more particularly conic programming) tools
for solving limit analysis has now emerged as the state-of-the art method,
some extensions have already been proposed in the literature to apply them
also to closely related problems. One can mention, for instance,
elastoplasticity [35, 34], viscoplasticity for yield stress fluids [11, 4],
contact in granular materials [33, 62], limit analysis-based topology
optimization [29, 25], etc. The present framework is sufficiently general to
also extend to these related problems, see for instance the application to
viscoplastic fluids in [5]. Other situations appear however more difficult to
include such as non-associative behaviours or geometrical non-linearities
since such problems cannot be formulated as convex optimization problems
anymore. Nonetheless, some works have already proposed some iterative
strategies for tackling non-associativity [24, 32, 51], it would therefore be
interesting to pursue in this direction.

## References

  * Anderheggen and Knöpfel [1972] Anderheggen, E., Knöpfel, H., 1972\.  Finite element limit analysis using linear programming.  International Journal of Solids and Structures 8, 1413–1431.
  * Ben-Tal and Nemirovski [2001] Ben-Tal, A., Nemirovski, A., 2001\.  Lectures on modern convex optimization: analysis, algorithms, and engineering applications. volume 2.  Siam.
  * Bisbos and Pardalos [2007] Bisbos, C., Pardalos, P., 2007\.  Second-order cone and semidefinite representations of material failure criteria.  Journal of Optimization Theory and Applications 134, 275–301.
  * Bleyer [2018] Bleyer, J., 2018.  Advances in the simulation of viscoplastic fluid flows using interior-point methods.  Computer Methods in Applied Mechanics and Engineering 330, 368–394.
  * Bleyer [2019] Bleyer, J., 2019.  Automating the formulation and resolution of convex variational problems: applications from image processing to computational mechanics.  arXiv preprint arXiv:1911.13185 .
  * Bleyer [2020] Bleyer, J., 2020.  fenics-optim – Convex optimization interface in FEniCS.  URL: https://doi.org/10.5281/zenodo.3778848, doi:10.5281/zenodo.3778848.
  * Bleyer and de Buhan [2014a] Bleyer, J., de Buhan, P., 2014a.  A computational homogenization approach for the yield design of periodic thin plates. Part I: Construction of the macroscopic strength criterion.  International Journal of Solids and Structures 51, 2448–2459.
  * Bleyer and de Buhan [2014b] Bleyer, J., de Buhan, P., 2014b.  Lower bound static approach for the yield design of thick plates.  International Journal for Numerical Methods in Engineering 100, 814–833.
  * Bleyer and de Buhan [2016] Bleyer, J., de Buhan, P., 2016\.  A numerical approach to the yield strength of shell structures.  European Journal of Mechanics-A/Solids 59, 178–194.
  * Bleyer et al. [2015a] Bleyer, J., Le, C.V., de Buhan, P., 2015a.  Locking-free discontinuous finite elements for the upper bound yield design of thick plates.  International Journal for Numerical Methods in Engineering 103, 894–913.
  * Bleyer et al. [2015b] Bleyer, J., Maillard, M., De Buhan, P., Coussot, P., 2015b.  Efficient numerical computations of yield stress fluid flows using second-order cone programming.  Computer Methods in Applied Mechanics and Engineering 283, 599–614.
  * Boyd and Vandenberghe [2004] Boyd, S., Vandenberghe, L., 2004\.  Convex optimization.  Cambridge university press.
  * de Buhan [1986] de Buhan, P., 1986.  A fundamental approach to the yield design of reinforced soil structures.  Ph.D. thesis. Thèse d’Etat, Paris VI.
  * Capurso [1971] Capurso, M., 1971.  Limit analysis of continuous media with piecewise linear yeld condition.  Meccanica 6, 53–58.
  * Casciaro and Cascini [1982] Casciaro, R., Cascini, L., 1982\.  A mixed formulation and mixed finite elements for limit analysis.  International Journal for Numerical Methods in Engineering 18, 211–243.
  * Chen [2007] Chen, W.F., 2007.  Plasticity in reinforced concrete.  J. Ross Publishing.
  * Chen [2013] Chen, W.F., 2013.  Limit analysis and soil plasticity.  Elsevier.
  * Christiansen and Andersen [1999] Christiansen, E., Andersen, K., 1999\.  Computation of collapse states with von Mises type yield condition.  International Journal for Numerical Methods in Engineering 46, 1185–1202.
  * Dantzig et al. [1955] Dantzig, G.B., Orden, A., Wolfe, P., et al., 1955.  The generalized simplex method for minimizing a linear form under linear inequality restraints.  Pacific Journal of Mathematics 5, 183–195.
  * De Buhan et al. [2002] De Buhan, P., Freard, J., Garnier, D., Maghous, S., 2002\.  Failure properties of fractured rock masses as anisotropic homogenized media.  Journal of Engineering Mechanics 128, 869–875.
  * Frémond and Friaa [1978] Frémond, M., Friaa, A., 1978\.  Analyse limite. comparaison des méthodes statique et cinématique.  Comptes Rendus de l’Academie des Sciences, Paris 286, 107–110.
  * Germain [1973a] Germain, P., 1973a.  La méthode des puissances virtuelles en mécanique des milieux continus, premiere partie: théorie du second gradient.  Journal de mécanique 12, 235–274.
  * Germain [1973b] Germain, P., 1973b.  The method of virtual power in continuum mechanics. part 2: Microstructure.  SIAM Journal on Applied Mathematics 25, 556–575.
  * Gilbert et al. [2006] Gilbert, M., Casapulla, C., Ahmed, H., 2006.  Limit analysis of masonry block structures with non-associative frictional joints using linear programming.  Computers & structures 84, 873–887.
  * Herfelt et al. [2019] Herfelt, M.A., Poulsen, P.N., Hoang, L.C., 2019.  Strength-based topology optimisation of plastic isotropic von mises materials.  Structural and Multidisciplinary Optimization 59, 893–906.
  * Hill [1950] Hill, R., 1950.  The mathematical theory of plasticity.  Clarendon Press, Oxford.
  * Ilyushin [1956] Ilyushin, A., 1956.  Plasticité.  Eyrolles, Paris.
  * Johansen [1962] Johansen, K., 1962.  Yield-line theory.  Cement and Concrete Association London.
  * Kammoun and Smaoui [2014] Kammoun, Z., Smaoui, H., 2014\.  A direct approach for continuous topology optimization subject to admissible loading.  Comptes Rendus Mécanique 342, 520–531.
  * Karmarkar [1984] Karmarkar, N., 1984.  A new polynomial-time algorithm for linear programming, in: Proceedings of the sixteenth annual ACM symposium on Theory of computing, ACM. pp. 302–311.
  * Krabbenhoft and Damkilde [2002] Krabbenhoft, K., Damkilde, L., 2002\.  Lower bound limit analysis of slabs with nonlinear yield criteria.  Computers & structures 80, 2043–2057.
  * Krabbenhoft et al. [2012a] Krabbenhoft, K., Karim, M., Lyamin, A., Sloan, S., 2012a.  Associated computational plasticity schemes for nonassociated frictional materials.  International Journal for Numerical Methods in Engineering 90, 1089–1117.
  * Krabbenhoft et al. [2012b] Krabbenhoft, K., Lyamin, A., Huang, J., da Silva, M.V., 2012b.  Granular contact dynamics using mathematical programming methods.  Computers and Geotechnics 43, 165–176.
  * Krabbenhøft et al. [2007] Krabbenhøft, K., Lyamin, A., Sloan, S., 2007.  Formulation and solution of some plasticity problems as conic programs.  International Journal of Solids and Structures 44, 1533–1549.
  * Krabbenhoft et al. [2007] Krabbenhoft, K., Lyamin, A., Sloan, S., Wriggers, P., 2007\.  An interior-point algorithm for elastoplasticity.  International Journal for Numerical Methods in Engineering 69, 592–626.
  * Krabbenhoft et al. [2005] Krabbenhoft, K., Lyamin, A.V., Hjiaj, M., Sloan, S.W., 2005\.  A new discontinuous upper bound limit analysis formulation.  International Journal for Numerical Methods in Engineering 63, 1069–1088.
  * Krabbenhøft et al. [2008] Krabbenhøft, K., Lyamin, A.V., Sloan, S.W., 2008.  Three-dimensional mohr–coulomb limit analysis using semidefinite programming.  Communications in Numerical Methods in Engineering 24, 1107–1119.
  * Lobo et al. [1998] Lobo, M.S., Vandenberghe, L., Boyd, S., Lebret, H., 1998\.  Applications of second-order cone programming.  Linear algebra and its applications 284, 193–228.
  * Lysmer [1970] Lysmer, J., 1970.  Limit analysis of plane problems in soil mechanics.  Journal of Soil Mechanics & Foundations Div 96, 1311–1334.
  * Makrodimopoulos [2010] Makrodimopoulos, A., 2010.  Remarks on some properties of conic yield restrictions in limit analysis.  International Journal for Numerical Methods in Biomedical Engineering 26, 1449–1461.
  * Makrodimopoulos and Martin [2007] Makrodimopoulos, A., Martin, C., 2007\.  Upper bound limit analysis using simplex strain elements and second-order cone programming.  International journal for numerical and analytical methods in geomechanics 31, 835–865.
  * Makrodimopoulos and Martin [2008] Makrodimopoulos, A., Martin, C., 2008\.  Upper bound limit analysis using discontinuous quadratic displacement fields.  Communications in Numerical Methods in Engineering 24, 911–927.
  * Martin and Makrodimopoulos [2008] Martin, C.M., Makrodimopoulos, A., 2008\.  Finite-element limit analysis of Mohr–Coulomb materials in 3d using semidefinite programming.  Journal of Engineering Mechanics 134, 339–347.
  * MOSEK ApS [2018] MOSEK ApS, ., 2018.  The MOSEK optimization API for Python 8.1.0.  URL: http://docs.mosek.com/8.1/pythonapi/index.htm.
  * MOSEK ApS [2019] MOSEK ApS, ., 2019.  MOSEK Modeling Cookbook.  https://docs.mosek.com/modeling-cookbook/index.html.
  * Nagtegaal et al. [1974] Nagtegaal, J.C., Parks, D.M., Rice, J., 1974.  On numerically accurate finite element solutions in the fully plastic range.  Computer methods in applied mechanics and engineering 4, 153–177.
  * Nayroles [1970] Nayroles, B., 1970.  Essai de théorie fonctionnelle des structures rigides plastiques parfaites.  Journal de Mécanique 9, 491–506.
  * Nesterov et al. [1994] Nesterov, Y., Nemirovskii, A., Ye, Y., 1994.  Interior-point polynomial algorithms in convex programming. volume 13.  SIAM.
  * Nielsen and Hoang [2016] Nielsen, M.P., Hoang, L.C., 2016\.  Limit analysis and concrete plasticity.  CRC press.
  * Pastor and Turgeman [1976] Pastor, J., Turgeman, S., 1976\.  Mise en œuvre numérique des méthodes de l’analyse limite pour les matériaux de von mises et de coulomb standards en déformation plane.  Mechanics Research Communications 3, 469–474.
  * Portioli et al. [2014] Portioli, F., Casapulla, C., Gilbert, M., Cascini, L., 2014\.  Limit analysis of 3d masonry block structures with non-associative frictional joints using cone programming.  Computers & Structures 143, 108–121.
  * Robinson [1971] Robinson, M., 1971.  A comparison of yield surfaces for thin shells.  International Journal of Mechanical Sciences 13, 345–354.
  * Salençon [1983] Salençon, J., 1983.  Calcul à la rupture et analyse limite.  Presses de l’Ecole Nationale des Ponts et Chaussées.
  * Salençon [2013] Salençon, J., 2013.  Yield Design.  London, Hoboken : ISTE Ltd., John Wiley & Sons, Inc.
  * Save [1995] Save, M., 1995.  Atlas of limit loads of metal plates, shells and disks.  North-Holland Series in Applied Mathematics & Mechanics, Elsevier Science Ltd.
  * Save et al. [1997] Save, M.A., Massonnet, C.E., de Saxce, G., 1997.  Plastic limit analysis of plates, shells, and disks.  North Holland.
  * Schlaich et al. [1987] Schlaich, J., Schäfer, K., Jennewein, M., 1987.  Toward a consistent design of structural concrete.  PCI journal 32, 74–150.
  * Sloan [1988] Sloan, S.W., 1988.  Lower bound limit analysis using finite elements and linear programming.  International Journal for Numerical and Analytical Methods in Geomechanics 12, 61–77.
  * Sloan and Kleeman [1995] Sloan, S.W., Kleeman, P.W., 1995\.  Upper bound limit analysis using discontinuous velocity fields.  Computer Methods in Applied Mechanics and Engineering 127, 293 – 314.
  * Suquet [1985] Suquet, P., 1985.  Elements of homogenization for inelastic solid mechanics.  Homogenization techniques for composite media 272, 193–278.
  * Vincent et al. [2018] Vincent, H., Arquier, M., Bleyer, J., de Buhan, P., 2018\.  Yield design-based numerical analysis of three-dimensional reinforced concrete structures.  International Journal for Numerical and Analytical Methods in Geomechanics 42, 2177–2192.
  * Zhang et al. [2014] Zhang, X., Krabbenhoft, K., Sheng, D., 2014.  Particle finite element analysis of the granular column collapse problem.  Granular Matter 16, 609–619.