Mixed 1D/2D modelling of a reinforced slope stability#
\(\newcommand{\bt}{\boldsymbol{t}} \newcommand{\bn}{\boldsymbol{n}} \newcommand{\bf}{\boldsymbol{f}} \newcommand{\dx}{\,\text{dx}} \newcommand{\dS}{\,\text{dS}} \newcommand{\bu}{\boldsymbol{u}} \newcommand{\bx}{\boldsymbol{x}} \newcommand{\bsig}{\boldsymbol{\sigma}} \newcommand{\jump}[1]{[\![#1]\!]}\)
In this demo, we illustrate how we can formulate a limit analysis problem involving 1D reinforcements in a 2D solid. In the present case, we consider the slope stability problem with three horizontal reinforcements. Note that the same approach could also be used for a 2D/3D mixed modelling. However, 1D/3D mixed modelling is much more involved and cannot be tackled using the present approach.
Mixed modelling approach#
We consider a bulk solid domain \(\Omega\) consisting of a material obeying the strength criterion \(G\), e.g. a plane strain Mohr-Coulomb criterion in the subsequent application. The domain also consists of a set of uniaxial reinforcements \(\Gamma\) of unit tangent vector \(\boldsymbol{t}\). The reinforcements are assumed to be perfectly flexible so that the internal forces consist only of a normal force \(N\). The latter is subjected to a uniaxial tension/compression criterion:
where \(N_0\) is the corresponding uniaxial strength.
Defining the different static variable including the bulk stress \(\bsig\) and the normal force \(N\) with fenics-optim is rather straightforward. It therefore essentially remains the question of imposing the equilibrium equations. We have inside the bulk \(\Omega\):
where \(\lambda\bf\) denotes the body force.
We also have on \(\Gamma\), the normal force equibrium:
where \(\jump{\bsig} = \bsig^\oplus - \bsig^\ominus\) denotes the stress jump through \(\Gamma\), \(\bn\) denotes the unit normal to \(\Gamma\) and \(s\) is the curvilinear coordinate along \(\Gamma\).
The above strong equilibrium equations can be enforced directly using a cell and a facet Lagrange multipliers respectively, as in a lower bound static approach. Alternatively, it can also be enforced weakly as in upper bound kinematic approaches written in a static form. In the latter case, considering a continuous virtual displacement \(\bu\), the virtual work principle reads:
and, since \(\dfrac{d\bu}{ds} = \dfrac{d\bu}{d\bx}\cdot\bt\), it also reads:
Implementation#
As mentioned above, we consider a vertical slope stability problem under self-weigth loading. The slope is reinforced by an array of uniaxial reinforcements of length \(L_\text{reinf}\) located at the heights \(H/4, H/2, 3H/4\). We first define the soil mesh and identify the regions corresponding to the fixed boundary conditions and the reinforcement locations. In particular, the reinforcements \(\Gamma\) will correspond to the facet tag 2.
[14]:
from dolfin import (
RectangleMesh,
Point,
FacetNormal,
near,
MeshFunction,
Measure,
dx,
AutoSubDomain,
Constant,
DirichletBC,
FunctionSpace,
VectorFunctionSpace,
plot,
)
from ufl import sym, dot, inner, avg, grad, pi, perp, outer
from fenics_optim import MosekProblem, to_mat
import fenics_optim.limit_analysis as la
import matplotlib.pyplot as plt
Nmesh = 100
L, H = 1.2, 1.0
L_reinf = 0.8
mesh = RectangleMesh(Point(0, 0), Point(L, H), Nmesh, Nmesh, "crossed")
def border(x, on_boundary):
return near(x[1], 0) or near(x[0], 0)
reinf_heights = [H / 4, H / 2, 3 * H / 4]
def reinforcement(x, on_boundary):
return any([near(x[1], Hi) for Hi in reinf_heights]) and (L - L_reinf <= x[0] <= L)
facets = MeshFunction("size_t", mesh, 1)
AutoSubDomain(border).mark(facets, 1)
AutoSubDomain(reinforcement).mark(facets, 2)
ds = Measure("ds", subdomain_data=facets, domain=mesh)
dS = Measure("dS", subdomain_data=facets, domain=mesh)
We then define the loading, the soil Mohr-Coulomb strength criterion and the reinforcement strength \(N_0\).
[15]:
f = Constant((0.0, -1.0))
c, phi = 1.0, 30 * pi / 180
mat = la.MohrCoulomb2D(c, phi)
N0 = 1.0
We then define the various function spaces which will be used to define the variables and constraint Lagrange multipliers. We use DG 0 spaces for the stress and the normal force and a CG 1 space for the displacement Lagrange multiplier. Note that since \(N\) lives on mesh facets only, we use a Discontinuous Lagrange Trace space. Although \(\Gamma\) corresponds only to a subspace of the mesh facets, the use of such a space will define \(N\) everywhere on the mesh facets.
As a result, a drawback of this approach is that dummy variables corresponding to \(N\) outside \(\Gamma\) are created. This is however not an issue as the above integral will be performed only on the measure dS(2) corresponding to \(\Gamma\) and because solvers such as Mosek are able to identify such unused variables and remove them from the final optimization problem. Finally, we use the variable lower bound lx and upper bound ux keyword to enforce the uniaxial strength
criterion \(|N|\leq N_0\).
[16]:
R = FunctionSpace(mesh, "R", 0)
Wsig = VectorFunctionSpace(mesh, "DG", 0, dim=3)
WN = FunctionSpace(mesh, "Discontinuous Lagrange Trace", 0)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
prob = MosekProblem("Limit analysis with reinforcements")
lamb, Sig, N = prob.add_var([R, Wsig, WN], lx=[None, None, -N0], ux=[None, None, N0])
We now express the above equilibrium weak form. Note that internal facet integrals must be restricted. In the present case, the integrand is only well defined on the facet itself so that the choice of the restriction operator has no influence. For the facet tagent vector, we use UFL’s perp function on the mesh FacetNormal.
[17]:
sig = to_mat(Sig)
n = FacetNormal(mesh)
tang = perp(n)
def equilibrium(v):
return (
inner(sig, sym(grad(v))) * dx
+ inner((N * outer(tang, tang))("+"), grad(v)("+")) * dS(2)
- lamb * dot(v, f) * dx
)
bc = DirichletBC(V, Constant((0.0, 0.0)), border)
prob.add_eq_constraint(V, A=equilibrium, name="Pseudo-mechanism", bc=bc)
We finish with defining the objective function \(\lambda\), applying the Mohr-Coulomb criterion to the bulk stress and calling the solver.
[18]:
prob.add_obj_func([1, None])
crit = mat.criterion(Sig)
prob.add_convex_term(crit)
prob.optimize(sense="maximize")
Matrix shape: (160402, 300201)
Number of cones: 40000
Problem
Name :
Objective sense : max
Type : CONIC (conic optimization problem)
Constraints : 160402
Cones : 40000
Scalar variables : 300201
Matrix variables : 0
Integer variables : 0
Optimizer started.
Presolve started.
Eliminator started.
Freed constraints in eliminator : 80000
Eliminator terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries : 2 time : 0.00
Lin. dep. - tries : 0 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.41
GP based matrix reordering started.
GP based matrix reordering terminated.
Problem
Name :
Objective sense : max
Type : CONIC (conic optimization problem)
Constraints : 160402
Cones : 40000
Scalar variables : 300201
Matrix variables : 0
Integer variables : 0
Optimizer - threads : 2
Optimizer - solved problem : the primal
Optimizer - Constraints : 40000
Optimizer - Cones : 40001
Optimizer - Scalar variables : 120200 conic : 120002
Optimizer - Semi-definite variables: 0 scalarized : 0
Factor - setup time : 0.54 dense det. time : 0.06
Factor - ML order time : 0.03 GP order time : 0.21
Factor - nonzeros before factor : 3.18e+05 after factor : 1.58e+06
Factor - dense dim. : 3 flops : 1.96e+08
Factor - GP saved nzs : 6.64e+04 GP saved flops : 5.87e+07
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 2.2e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 1.13
1 1.6e-02 7.1e-03 2.6e-03 1.60e-01 -1.509808514e-02 -1.403595471e-01 7.1e-03 1.32
2 1.3e-03 5.6e-04 6.0e-05 9.82e-01 4.918123109e-02 3.854400084e-02 5.6e-04 1.49
3 2.9e-04 1.3e-04 7.1e-06 8.72e-01 7.138467507e-01 7.109820659e-01 1.3e-04 1.64
4 1.8e-04 8.2e-05 4.7e-06 2.94e-01 2.684262616e+00 2.681212927e+00 8.2e-05 1.79
5 9.9e-05 4.4e-05 2.6e-06 2.62e-03 6.636029143e+00 6.632656524e+00 4.4e-05 1.92
6 5.0e-05 2.2e-05 1.1e-06 1.06e-01 1.022119237e+01 1.021893544e+01 2.2e-05 2.04
7 2.2e-05 9.9e-06 3.4e-07 4.14e-01 1.338418404e+01 1.338307805e+01 9.9e-06 2.16
8 1.2e-05 5.4e-06 1.3e-07 6.91e-01 1.478708499e+01 1.478651457e+01 5.4e-06 2.28
9 6.2e-06 2.8e-06 4.7e-08 8.32e-01 1.558059633e+01 1.558033238e+01 2.8e-06 2.40
10 1.9e-06 8.3e-07 7.1e-09 9.22e-01 1.610883868e+01 1.610877313e+01 8.3e-07 2.51
11 6.0e-07 2.7e-07 1.2e-09 9.79e-01 1.626576268e+01 1.626574448e+01 2.7e-07 2.66
12 2.1e-07 9.2e-08 2.3e-10 9.94e-01 1.631459476e+01 1.631458899e+01 9.2e-08 2.78
13 5.4e-08 2.4e-08 2.9e-11 9.98e-01 1.633333451e+01 1.633333317e+01 2.4e-08 2.90
14 1.8e-08 8.1e-09 5.8e-12 9.99e-01 1.633740925e+01 1.633740881e+01 8.1e-09 3.05
15 5.1e-09 2.3e-09 8.4e-13 1.00e+00 1.633890289e+01 1.633890278e+01 2.3e-09 3.19
Optimizer terminated. Time: 3.30
Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: 1.6338902894e+01 nrm: 4e+01 Viol. con: 2e-10 var: 0e+00 cones: 0e+00
Dual. obj: 1.6338902775e+01 nrm: 1e+03 Viol. con: 0e+00 var: 3e-08 cones: 0e+00
[18]:
16.338902894084487
The collapse mechanism \(\bu\) clearly shows the restraining effect of the reinforcements.
[21]:
u = prob.get_lagrange_multiplier("Pseudo-mechanism")
plot(-0.004 * u, mode="displacement")
plt.show()
Calling FFC just-in-time (JIT) compiler, this may take some time.
Ignoring precision in integral metadata compiled using quadrature representation. Not implemented.
[ ]: