<figure><img src="image_1.png" /></figure>

HAL
open science

# A majoration method to obtain true upper bounds for the limit analysis of von Mises plates

Jérémy Bleyer, Patrick de Buhan

<figure><img src="image_2.png" /></figure>

## To cite this version:

Jérémy Bleyer, Patrick de Buhan. A majoration method to obtain true upper bounds for the limit analysis of von Mises plates. International Conference on Computational Mechanics, Mar 2013, Durham, United Kingdom. pp.16-CM13. hal-00807542

HAL Id: hal-00807542

https://enpc.hal.science/hal-00807542v1

Submitted on 8 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A majoration method to obtain true upper bounds for the limit analysis of von Mises plates

J. Bleyer<sup>a,*</sup>, P. de Buhan<sup>a</sup>

<sup>a</sup>Université Paris-Est, Laboratoire Navier, Ecole des Ponts ParisTech-IFSTTAR-CNRS (UMR 8205)
6-8 av Blaise Pascal, Cité Descartes, 77455 Champs-sur-Marne, FRANCE

## Abstract

Limit analysis has become an appealing method to compute the limit load of various kind of structures in engineering practice. Thin plates in bending are an important class of structures for which engineers are interested in obtaining such limit loads. Many methods like mixed element approach, meshless methods or even displacement-based finite element methods provide only approximate upper or lower bounds due to discretization or integration approximations. This paper considers the upper bound limit analysis approach for Love Kirchhoff plates obeying a von Mises criterion. With the use of nonconforming finite elements, we propose a novel approach based on a majoration of the plastic dissipation at the finite element level to obtain true upper bounds without any integration approximation. We obtain an already-integrated dissipation operator which avoids any use of numerical quadrature leading to a true upper bound. The problem is formulated as a standard second-order cone programming and solved with the Mosek software package. We show, through various numerical examples, that this method is more efficient than traditional formulation and gives very accurate upper bounds. By way of example, for the uniformly loaded square plate with clamped supports, we obtain for the limit load $\lambda = 44.196 \frac{m_p}{qL^2}$, which is the best upper bound estimate computed so far for this problem.

**Keywords:** limit analysis ; yield design ; upper bound ; thin plates in bending ; finite element method ; second order cone programming

## 1. Introduction

Thanks to important progress in mathematical programming and development of robust optimization methods, limit analysis problems have gained increasing attention in the last two decades. Various method involving finite element discretization have been proposed to solve lower bound and upper bound problems [1, 2, 3, 4]. Mixed method using both stress and displacement discretization have also been proposed [5, 6].

Plates in bending are an important class of structures for which engineers need to be able to compute limit loads. Yield line theory, based on rigid mechanisms separated by yield lines of rotation discontinuity, has been the first method to compute analytical or numerical upper bounds [7]. In order to improve the quality of these bounds, plastic dissipation through curvature has also been taken into account in conjunction with the use of finite element method [8, 9, 10]. However, few authors used finite elements combining plastic dissipation through curvature as well as rotation discontinuities. A recent paper under submission considered such non-conforming elements which proved to be much more efficient [11]. Traditional discretization procedure of the upper bound limit analysis approach with finite elements require to compute plastic dissipation terms as integrals of nonlinear functions with a Gauss-Legendre quadrature. Therefore, limit loads resulting from this approximation cannot be considered anymore as true upper bounds of the exact limit load, but only as approximations. Other methods like mixed element approach with particular choice of stress and displacement spaces or meshless methods can also lead to approximate limit loads instead of true bounds [12]. This can be very inconvenient for engineers since the exact limit load cannot be rigorously bracketed by a lower and an upper bound. Fortunately, numerical experiments tend to prove that the use of sufficient Gauss points at element level still provide accurate results.

In this work, we suggest to take advantage of the specific form of the dissipation in the case of a von Mises criterion to obtain a new limit analysis problem. The obtained problem can be easily integrated before any optimization procedure, thus avoiding any use of approximate quadrature and leading to a true upper bound. Moreover, this problem can still be cast as a second-order cone programming (SOCP) problem, for which efficient interior point solvers have been implemented in commercial codes such as the Mosek software package [13].

In the first section, the general formulation of upper bound limit analysis in the case of bending plates is briefly recalled.

<footnote>*Laboratoire Navier, Ecole des Ponts ParisTech Tel : +33 (0)1 64 15 36 59
Email address: jeremy.bleyer@enpc.fr (J. Bleyer)</footnote>
Finite element discretization procedure is addressed in section 2 and the proposed majoration method is presented in section 3. In section 4, numerical results on classical problems in bending plate limit analysis are investigated in order to assess the efficiency of the proposed method compared to the classical use of Gauss quadrature.

## 2. Upper bound kinematic formulation of bending plates

The purpose of the present section is to formulate the upper bound kinematic approach of yield design in the specific case of bending plates.

We consider a plate occupying a geometrical domain S in the Oxy-plane. Its transverse virtual velocity field will be noted u. In the following, only the case of bending plates will be investigated, that is the plate is assumed to be infinitely resistant to both shear and membrane forces. It can be shown that this assumption implies that the well-known Love-Kirchhoff condition must be satisfied by the velocity fields.

As a consequence, the implementation of the yield design upper bound kinematic approach for bending plates is based upon considering distributions of kinematically admissible transverse velocity fields u subject to the following properties:

* u is everywhere continuous;

* its second gradient (rate of curvature $\bar{\chi}$) can be calculated everywhere, except on a number of discontinuity lines $\Gamma$ (called yield lines) where jumps of angular velocities have to be taken into account.

We assume that the plate is loaded by a transversal loading distribution $q(x, y)e_z$ and that its boundary $\partial S$ can be decomposed into four parts $\partial S_i$ such that $\partial S = \cup \partial S_i \forall i = 0, \dots, 3$ with $\partial S_i \cap \partial S_j = \emptyset \forall i \neq j$ and:

* on $\partial S_1 : u = 0$ (simply supported)

* on $\partial S_2 : \frac{\partial u}{\partial n} = 0$ (fixed rotation, free displacement)

* on $\partial S_3 : u = 0$ and $\frac{\partial u}{\partial n} = 0$ (clamped boundary)

* on $\partial S_0 :$ free boundary

We note KA the set of kinematically admissible velocity fields. The work of external forces is here :

$$P_{(e)}(u) = \int_S q(x, y)u(x, y)dS \qquad (1)$$

While the maximum resisting work is $\forall u \in KA$:

$$P_{rm}(u) = P_{rm}^{curv}(u) + P_{rm}^{disc}(u) = \int_{S} \pi(\chi(x, y))dS + \int_{\Gamma} \pi(\delta\theta_{n}, \underline{n})dl \qquad (2)$$

where in the case of a von Mises criterion, we have:

$$\pi(\chi(x, y)) = \frac{2m_p}{\sqrt{3}} \sqrt{\chi_{xx}^2 + \chi_{yy}^2 + \chi_{xx}\chi_{yy} + \chi_{xy}^2} \\ \pi(\delta\theta_n, \underline{n}) = \frac{2m_p}{\sqrt{3}} |\delta\theta_n|$$

where $\delta\theta_n = [\frac{\partial u}{\partial n}]$ is the angular velocity jump through $\Gamma$ along its normal $\underline{n}$.

The maximum resisting work consists of two terms: $P_{rm}^{curv}$ comes from the contribution of the dissipation through curvature, while $P_{rm}^{disc}$ comes from the contribution of the dissipation through a discontinuity of the angular velocity $\delta\theta_n$.

The kinematic formulation for the ultimate load factor $\lambda^*$ is then given by the following optimization problem:

$$\lambda^* = \min_{\forall u \in KA} P_{rm}(u) \\ P_{(e)}(u) = 1 \qquad (3)$$
### 3. Finite element discretization

In this section, we will investigate the finite element discretization of problem (3). The way of computing the different contributions to the maximum resisting work will be particularly highlighted.

We consider the general case of a triangular finite element e with s degrees of freedom $q_i^{(e)}$ and given shape functions $N_i(\xi, \eta)$ $\forall i = 1, \dots, s$ defined in the $(\xi, \eta)$-plane of a reference element T through an affine transformation. We have therefore, for $\zeta = (\xi, \eta)$:

$$ \{\chi(\zeta)\} = \begin{cases} \partial_{xx}u(\zeta) \\ \partial_{yy}u(\zeta) \\ 2\partial_{xy}u(\zeta) \end{cases} = [b(\zeta)]\{q_e\} \quad \text{and} \quad \delta\theta_n = \langle\delta\Theta_j(\tau)\rangle\{q_j\} $$

where $\{q_e\}$ collects the dof of element e and $\{q_j\}$ collects the dof of the elements sharing a given edge j. $\delta\theta_n$ being, thus, the discontinuity of the angular velocity around edge j.

#### 3.1. Contribution of the curvature to the dissipation

The contribution of the curvature term $P_{rm}^{curv}$ to the total dissipation.

$$ P_{rm}^{curv} = m_p \sum_{e=1}^{N_e} \det[J_e] \int_{T} \sqrt{\{\chi(\zeta)\}^T [Q] \{\chi(\zeta)\}} dS \quad \text{with} \quad [Q] = \frac{1}{3} \begin{bmatrix} 4 & 2 & 0 \\ 2 & 4 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad (4) $$

where $N_e$ is the total number of elements and $\det[J_e]$ the jacobian of the geometric transformation for e.

#### 3.2. Contribution of rotation discontinuities to the dissipation

The contribution of all $N_d$ active edges to $P_{rm}^{disc}$ is:

$$ P_{rm}^{disc} = \frac{2m_p}{\sqrt{3}} \sum_{j=1}^{N_d} \det[J_j] \int_{L} |\langle \delta\Theta_j(\tau) \rangle\{q_j\}| d\tau \quad (5) $$

where $\det[J_j]$ is the jacobian of the geometric transformation for j to a reference line L.

**Remark:** Active edges represent edges which can potentially contribute to $P_{rm}$ via a rotation discontinuity. This involves internal edges of the mesh and also edges belonging to the $\partial S_2 \cup \partial S_3$ part of the boundary where the normal rotation is set equal to zero.

#### 3.3. Discrete minimization problem

Finally, the discrete minimization problem can be expressed as:

$$ \lambda_h = \min_{\{U\}} m_p \sum_{e=1}^{N_e} \det[J_e] \int_{T} \sqrt{\{\chi(\zeta)\}^T [Q] \{\chi(\zeta)\}} dS + \frac{2m_p}{\sqrt{3}} \sum_{j=1}^{N_d} \det[J_j] \int_{L} |\langle \delta\Theta_j(\tau) \rangle\{q_j\}| d\tau \quad (6) \\ \text{subject to } \{F\}^T \cdot \{U\} = 1 $$

with appropriate boundary conditions which can be taken into account during the assembling process of the global vectors, $\{F\}$ and $\{U\}$ being the global load and displacement vectors.

Due to the finite element interpolation of the transverse displacement, $u$ belongs here to a smaller space included in the real optimization space KA. Hence, the corresponding load factor $\lambda_h$ is a strict upper bound of the ultimate load factor: $\lambda^* \le \lambda_h$.

However, elementary integrals still remain in the formulation and cannot be computed as such. A standard method to compute these integrals is to use a quadrature rule (Gauss-Legendre mostly) with a fixed number of integration points. Since integrands are non-polynomial, the computation is not exact and there is no way to tell beforehand if it leads to an approximation by excess or by default. Hence, letting:

$$ \lambda_G = \min_{\{U\}} m_p \sum_{e=1}^{N_e} \det[J_e] \sum_{g=1}^{n} \omega_g \sqrt{\{\chi(\zeta_g)\}^T [Q] \{\chi(\zeta_g)\}} + \frac{2m_p}{\sqrt{3}} \sum_{j=1}^{N_d} \det[J_j] \sum_{g'=1}^{m} \omega_{g'} |\langle \delta\Theta_j(\tau_{g'}) \rangle\{q_j\}| \quad (7) \\ \text{subject to } \{F\}^T \cdot \{U\} = 1 $$

be the load factor obtained with a Gauss quadrature with $n$ points on T and $m$ points on L, the previous upper bound result is now lost and we can only say that $\lambda_G$ is an approximation of the true limit load $\lambda^*$.
## 4. Majoration of the maximum resisting work

In this section, we propose a new method to compute an exact upper bound in the particular case of the von Mises criterion by approximating the previous integrals by excess.

Our idea is to use the Cauchy-Schwartz inequality on $P_{rm}^{curv}$:

$$P_{rm}^{curv} \leq m_p \sum_{e=1}^{N_e} \det[J_e] \left( |T| \int_T \{\chi(\zeta)\}^T [Q] \{\chi(\zeta)\} dS \right)^{1/2}$$

With $[B(\zeta)] = [b(\zeta)]^T [Q] [b(\zeta)]$, we have:

$$P_{rm}^{curv} \leq m_p \sum_{e=1}^{N_e} \det[J_e] \left( |T| \{q_e\}^T \int_T [B(\zeta)] dS \{q_e\} \right)^{1/2} = m_p \sqrt{\{q_e\}^T [\widehat{B}_e] \{q_e\}}$$

where $[\widehat{B}_e] = (\det[J_e])^2 |T| \int_T [B(\zeta)] dS$ is some kind of average operator on element $e$. Since $[B(\zeta)]$ is polynomial (shape functions and as a consequence their derivatives are polynomial), the main interest is that $[\widehat{B}_e]$ can be computed exactly without any quadrature approximation.

Following the same procedure for the discontinuity term $P_{rm}^{disc}$ we obtain the following new optimization problem:

$$\lambda_U = \min_{\{U\}} m_p \sum_{e=1}^{N_e} \sqrt{\{q_e\}^T [\widehat{B}_e] \{q_e\}} + \frac{2m_p}{\sqrt{3}} \sum_{j=1}^{N_d} \sqrt{\{q_j\}^T [\widehat{D}_j] \{q_j\}} \quad (8)$$

$$\text{subject to } \{F\}^T \cdot \{U\} = 1$$

where $[\widehat{D}_j] = (\det[J_j])^2 |L| \int_L \{\delta\Theta_j\} \{\delta\Theta_j\}^T dl$ can also be computed exactly.

Evidently, we have $\lambda^* \leq \lambda_h \leq \lambda_U$. $\lambda_U$ is, thus, a true upper bound of the ultimate load factor.

It is worth mentioning that problem (8) is of the same form as problem (7) i.e. a sum of $\| \cdot \|_2$-norms which can easily be cast as an SOCP optimization problem.

The second main interest of the proposed method is that it avoids the use of Gauss points in the optimization formulation since both average operators are a priori computed. Consequently, this tends to reduce the total number of variables required in the optimization process which is very interesting in terms of computational performances.

## 5. Numerical examples

In this section, we will provide some examples of the performance of this method for computing very accurate load factors despite the overestimation made through the majoration procedure.

### 5.1. Numerical choices

It has recently been shown that the most performing finite element for limit analysis of plates in bending seems to be the cubic Hermite triangle (H3) [11]. This element is a non-conforming triangular finite element with $s = 10$ degrees of freedom (the transverse displacement and the two rotations at corner nodes, plus the transverse displacement at centroid) with a cubic interpolation of the transverse displacement. We will compare the performance of the majoration method to the method with a Gauss-Legendre quadrature with $n = m = 3$. The H3 element will thus be adopted in both cases.

As already mentioned, both resulting optimization problems have the same form and can be rewritten as a standard second-order cone programming problem which can be solved very efficiently by industrial software packages like Mosek (version 6.0).

### 5.2. Circular plate with clamped boundaries

Our first example will be the classical problem of a uniformly loaded circular plate with clamped supports (Fig. 1(a)). The exact load factor for this problem is known to be $\lambda^* = \frac{q^* R^2}{m_p} = 12.5$ where $R$ is the radius of the plate [7]. Only a quarter of plate was modeled with appropriate boundary conditions to account for symmetry.

Unstructured meshes of different size were considered by varying the number $M$ of elements along each radius (from $M = 3$ to $M = 29$). Both methods were tested on the same meshes. In Figure 1(b), the relative error to exact load $\lambda^* = 12.5$ is plotted for both methods along the corresponding number of minimization variables $N$.
<figure><img src="image_4.png" /><figcaption>Figure 1: Circular clamped plate problem</figcaption></figure>

We can clearly observe that for coarse meshes (values of $N$ smaller than $10^3$), limit loads obtained with the majoration method are higher than those given by the Gauss quadrature but are still a very good estimation, the error being less than 2.5% for the smaller mesh with $M = 3$. With increasing mesh size, we observe that the majoration method gives roughly the same values as the Gauss method and seems to be even more efficient for fine meshes. We observed, indeed, some numerical difficulties in the optimizing process with the Gauss method leading to higher bounds with finer meshes whereas those difficulties were not encountered with the majoration method. It is also worth mentioning that the optimization stage has also been carried out much faster in general.

### 5.3. Square plate

The second example will treat the case of a $L \times L$ square plate with simply supported and clamped boundaries. For the simply supported case, the best numerical upper bound seems to be $\lambda^* = \frac{m_p}{qL^2} = 25.02$ computed in [9] with conforming BFS elements. For the clamped case, the best upper bound seems to be $\lambda^* = 44.29$ computed in [11] with the considered H3 elements.

In the following, computed limit loads were compared to these best values. Computations were performed on a quarter of plate with $M$ elements along each side. Relative error are represented in figure 2(a) for the simply supported case and in figure 2(b) for the clamped case.

The same observations can be made as for the circular plate examples. The majoration gives higher limit loads for smaller mesh but results clearly improve with mesh refinement. Numerical difficulties to improve the bounds were still observed for fine meshes with the Gauss method but not with the majoration method. With the finest mesh, it was even possible to obtain a relative error of 0.25% for the simply supported case compared to Capsoni and Corradi's bound whereas, for the clamped case, the best upper bound obtained was improved to the value of $\lambda = 44.196$. Moreover, this result is assured to be a true upper bound with this majoration method.

<figure><img src="image_5.png" /><figcaption>Figure 2: Square plate problem results</figcaption></figure>
## 6. Conclusion

A new method has been presented to compute true upper bounds for the limit analysis of bending plates in the particular case of a von Mises yield criterion. This method is based on the majoration of elementary dissipation terms at the finite element level which avoids any use of quadrature rule. Leading to an approximation by excess of the maximum resisting work, the computed limit loads are true upper bounds to the exact limit load.

Through numerical examples, it has been demonstrated that the use of non-conforming cubic Hermite triangles (H3), combined with the SOCP minimization procedure, produces excellent results in terms of numerical accuracy for both methods. Moreover, compared to the classical method involving Gauss quadrature, this majoration method presents some important advantages. First, the upper bound status is theoretically assured and computed upper bounds are quite close from bounds obtained with a Gauss quadrature with a quite similar convergence rate (even better in some cases). One can observe better stability towards the optimization process (no numerical difficulties have been observed to obtain an optimum). Finally, for a similar number of optimization variables, the optimization stage is faster. These different aspects suggest that the resulting optimization problem seems to be better posed, which enables to improve some upper bounds. For instance, this method allowed to improve the best known upper bound for the uniformly loaded clamped square plate.

## References

[1] S. Sloan, Lower bound limit analysis using finite elements and linear programming, International Journal for Numerical and Analytical Methods in Geomechanics 12 (1) (1988) 61–77.

[2] S. Sloan, Upper bound limit analysis using finite elements and linear programming, International Journal for Numerical and Analytical Methods in Geomechanics 13 (3) (1989) 263–282.

[3] A. Makrodimopoulos, C. Martin, Lower bound limit analysis of cohesive-frictional materials using second-order cone programming, International Journal for Numerical Methods in Engineering 66 (4) (2006) 604–634.

[4] A. Makrodimopoulos, C. Martin, Upper bound limit analysis using simplex strain elements and second-order cone programming, International Journal for Numerical and Analytical Methods in Geomechanics 31 (6) (2007) 835–865.

[5] J. Munoz, J. Bonet, A. Huerta, J. Peraire, Upper and lower bounds in limit analysis: adaptive meshing strategies and discontinuous loading, International Journal for Numerical Methods in Engineering 77 (4) (2008) 471–501.

[6] H. Ciria, J. Peraire, J. Bonet, Mesh adaptive computation of upper and lower bounds in limit analysis, International Journal for Numerical Methods in Engineering 75 (8) (2008) 899–944.

[7] M. Save, C. Massonnet, G. de Saxce, Plastic limit analysis of plates, shells, and disks, Vol. 43, North Holland, 1997.

[8] P. Hodge Jr, T. Belytschko, Numerical methods for the limit analysis of plates, Journal of Applied Mechanics 35 (1968) 796.

[9] A. Capsoni, L. Corradi, Limit analysis of plates- a finite element formulation, Structural Engineering and Mechanics 8 (4) (1999) 325–341.

[10] C. Le, H. Nguyen-Xuan, H. Nguyen-Dang, Upper and lower bound limit analysis of plates using fem and second-order cone programming, Computers & Structures 88 (1-2) (2010) 65–73.

[11] J. Bleyer, P. de Buhan, On the performance of non-conforming finite elements for the upper bound limit analysis of plates, International Journal for Numerical Methods in Engineering - accepted (2012).

[12] J. Muñoz, A. Huerta, J. Bonet, J. Peraire, A note on upper bound formulations in limit analysis, International Journal for Numerical Methods in Engineering.

[13] Mosek, The mosek optimization toolbox for matlab manual (2008).
URL http://www.mosek.com/