Jeremy Bleyer
email: jeremy.bleyer@enpc.fr

Patrick De Buhan

A numerical approach to the yield strength of shell structures

Keywords: yield design, limit analysis, shells, generalized strength criteria, finite element method, second-order cone programming

This work investigates the formulation of lower and upper bound finite elements for the yield design (or limit analysis) of shell structures. The shell geometry is first discretized into triangular planar facets so that previously developed lower bound equilibrium and upper bound kinematic plate finite elements can be coupled to membrane elements. The other main novelty of this paper relies on the formulation of generalized strength criteria for shells in membrane-bending interaction via an implicit upscaling procedure. This formulation provides a natural strategy for constructing lower and upper bound approximations of the exact shell strength criterion and are particularly well suited for a numerical implementation using second-order cone programming tools. By combining these approximate strength criteria to the previously mentioned finite elements, rigorous lower and upper bound ultimate load estimates for shell structures can be computed very efficiently. Different numerical examples illustrate the accuracy as well as the generality and versatility of the proposed approach.

Introduction

Computational direct methods to estimate the ultimate load of a various range of structures using yield design (or limit analysis in the context of an elastic perfectly plastic behavior) theory have gained increasing attention in the last decades thanks to the development of efficient optimization solvers, in particular interior-points algorithm for conic programming problems. The static approach, which consists in maximizing the load multiplier over a set of statically admissible stress fields, while statisfying the local strength at each point of the structure, enables to obtain a lower bound estimate for the ultimate load. Conversely, the dual approach, namely the kinematic approach, which consists in minimizing the maximum resisting work (or plastic dissipation in the context of limit analysis) over a set of kinematically admissible virtual velocity fields, enables to obtain upper bound estimates for the ultimate load. The finite element method can be implemented in the context of both approaches using equilibrium elements for the static approach or kinematic finite elements (which may include potential discontinuities) for the kinematic approach. The subsequent optimization problems can then be solved using dedicated convex programming solvers, depending on the mathematical expression of the strength criterion. The general framework of the yield design theory [START_REF] Salençon | Calcul à la rupture et analyse limite[END_REF][START_REF] Salençon | Yield Design[END_REF] can be applied to a wide variety of mechanical models : 2D/3D continuous media, beam structures, plates and slabs, shells, etc. As regards the latter structures, one key issue is then to formulate a strength criterion in terms of the generalized internal forces of the considered model. From a numerical point of view, attention has also been devoted to the formulation of such criteria using conic (in particular second-order cone) constraints [START_REF] Bisbos | Second-order cone and semidefinite representations of material failure criteria[END_REF][START_REF] Makrodimopoulos | Remarks on some properties of conic yield restrictions in limit analysis[END_REF] in order to employ efficient second-order cone programming (SOCP) solvers, such as the Mosek software package [START_REF] Mosek | The Mosek optimization software[END_REF]. Focusing more specifically on the limit analysis/yield design of shell structures, early works have been mainly dedicated to deriving analytical lower and upper bound estimates for simple structures [START_REF] Hodge | The rigid-plastic analysis of symmetrically loaded cylindrical shells[END_REF][START_REF] Hodge | Plastic analysis of structures[END_REF][START_REF] Prager | Problems in Continuum Mechanics[END_REF][START_REF] Hodge | Limit analysis of rotationally symmetric plates and shells[END_REF][START_REF] Save | Plastic limit analysis of plates, shells, and disks[END_REF]. An important number of such solutions is compiled in [START_REF] Save | Atlas of limit loads of metal plates, shells and disks[END_REF]. In this context, the possible use of a generalized strength criterion for shells in membrane-bending interaction has been discussed [START_REF] Ilyushin | Plasticité[END_REF], and various approximate criteria have been proposed to simplify the analysis [START_REF] Robinson | A comparison of yield surfaces for thin shells[END_REF]. From a computational point of view, papers dedicated to a numerical implementation of limit analysis applied to such structures remain quite scarce. Limit anal-ysis of axi-symmetric shells using non-linear programming has been proposed in [START_REF] Hung | Bornes quasiinferieures et bornes superieures de la pression de ruine des coques de revolution par la methode des elements finis et par la programmation non-lineaire[END_REF]. Shakedown analysis (extension of limit analysis to cyclic loadings) of axisymmetric shells using a linearized yield surface has been investigated in (Franco and Ponter, 1997a,b) and extended in (Franco et al., 1997) including error analysis and adaptive remeshing. Sequential limit analysis has been used to investigate the post-collapse behavior of a structure by updating the geometry from using collapse mechanisms obtained from the kinematic approach. Such analyses have been performed in (Corradi andPanzeri, 2003, 2004) in the case of a von Mises criterion and in [START_REF] Raithatha | Rigid plastic model of incremental sheet deformation using second-order cone programming[END_REF] using Ilyushin approximate criterion (although some concerns with respect to the upper bound formulation have been raised in [START_REF] Makrodimopoulos | Comments on rigid plastic model of incremental sheet deformation using second-order cone programming[END_REF]). A continuum-based shell element has been developed to perform shakedown analysis on von Mises shells in [START_REF] Martins | A continuum-based mixed shell element for shakedown analysis[END_REF], but due to the mixed formulation, only approximate values of the limit loads can be obtained. Upper bound limit and shakedown analysis based on the exact Ilysuhin strength criterion has been considered in [START_REF] Tran | Upper bound limit and shakedown analysis of shells using the exact Ilyushin yield surface[END_REF], using a particular optimization procedure. Finally, one can mention the work of Bisbos and Papaioannou [START_REF] Bisbos | Shakedown analysis of FEMdiscretized steel shell structures under the Ilyushin yield criterion[END_REF], who used a Morley shell element and SOCP to solve problems involving Ilyushin approximate criterion. From all these previous works, it can first be observed that the derivation of lower bound shell elements is almost nonexistent, although it can be highly beneficial to obtain a reliable bracketing of the true ultimate load. Secondly, almost all existing works are focused on homogeneous metal shells, which can be modeled using the von Mises strength criterion. To the authors' opinion, this particular focus is partly due to the difficulty of deriving appropriate general strength criteria for other types of constitutive materials (reinforced concrete for example). Therefore, the present work aims at contributing to the derivation of both lower and upper bound yield design shell finite elements as well as proposing a general formulation of generalized shell strength criteria which can be used in conjunction with state-of-the-art conic programming solvers. The paper is organized as follows : in Section 2, the main features of the yield design analysis of shell structures are recalled, relying on the formulation of a strength condition expressed in terms of generalized stress variables, such as membrane forces and bending moments. Section 3 advocates the use of a piecewise linearization of the shell geometry by discretizing its curved surface into triangular facets, that is by replacing the initial shell by a continuous assemblage of triangular plates, for which lower and upper bound finite element formulations have been previously developed. Section 4 is then devoted to the key issue of formulating a generalized stress-based criterion for a multilayered shell, which can be obtained from solving a yield design problem, where the distribution of material local strength properties across the shell thickness is known. The numerical performance of such an upscal-ing procedure is favorably compared with classical existing solutions or approximations. Finally, Section 5 presents some illustrative applications of the whole numerical procedure, where the combination of the plate finite element formulation and strength criterion approximation, in the context of the lower or upper bound approach, leads to a single global SOCP optimization problem.

2. Yield design of shells : a brief outline 2.1. General formulation using the static approach

Referring to a Cartesian orthonormal frame (O; e x , e y , e z ), the shell occupies a two-dimensional manifold Ω. The shell geometry can be described locally by a unit normal ν and a tangent plane spanned by two unit vectors a 1 and a 2 (Figure 1).

The generalized internal forces of the shell model are described by a symmetric tensor N = N ij a i ⊗ a j of membrane forces, a symmetric tensor M = M ij a i ⊗ a j of bending moments and a vector V = V i a i of shear forces (i, j = 1, 2). Let then G(x) be the generalized strength domain of the shell at x ∈ Ω, which is a convex set in the (N , M , V ) space (of dimension 8).

Assuming that the shell external loading depends upon several loading parameters collected in a vector Q, the domain K of potentially safe loads Q, as introduced in the yield design theory [START_REF] Salençon | Calcul à la rupture et analyse limite[END_REF][START_REF] Salençon | Yield Design[END_REF], is defined as follows:

K = Q ; ∃(N , M , V ) SA with Q and ∀x ∈ Ω N (x), M (x), V (x) ∈ G(x) (1)
where a distribution of generalized internal forces (N , M , V ) must satisfy all local equilibrium, continuity and boundary conditions in order to be statically admissible (SA) with a given loading Q.

Dual formulation using the kinematic approach

A dual definition of the domain K of potentially safe loads using the kinematic approach is classically obtained by means of the virtual work principle. In the context of a shell model, the virtual kinematics of the shell is characterized at any point x ∈ Ω by (Figure 1) :

• a virtual velocity U = U x e x + U y e y + U z e z of the particle attached to the point. It will be convenient to express U in terms of its in-plane u = u i a i and out-of-plane wν components, so that : U = u + wν.

• a virtual angular velocity Θ = Θ x e x + Θ y e y + Θ z e z , representing the rotation of the microstructure attached to the same point. Again, it will be useful to express it in terms of in-plane and out-of-plane components : Θ = θ + ϑν with θ = θ i a i . Hence, the virtual work principle reads as :

(N , M , V ) S.A. with Q ⇐⇒ ∀(U , Θ) K.A. with q, P (e) (U , Θ) = Q • q = -P (i) (U , Θ)
where P (e) (U , Θ) = Q • q represents the virtual work of external loading ( q being the generalized kinematic parameters associated by duality to the loading parameters Q) and P (i) (U , Θ) the virtual work of internal forces which can be expressed as follows :

P (i) (U , Θ) = - Ω (N : + M : χ + V • γ)dΩ (2)
where are the virtual membrane strain rates, χ the virtual curvature strain rates and γ the virtual shear strain rates.

In the case when the virtual velocity fields U and Θ are discontinuous across a line Γ, the above expression (2) should be completed by the following additional line integral:

P [[ ]] (i) (U , Θ) = - Γ (N • n) • [[u]] + (M • n) • [[β]] + (V • n)[[w]] dΓ (3)
where β = ν ∧ θ and [[ * ]] represents the jump of a variable * when crossing Γ, of unit tangent vector t, along its unit normal n = ν ∧ t.

The so-called maximum resisting work (also called maximum plastic dissipation in the context of limit analysis) is now introduced. It is defined as :

P rm (U , Θ) = Ω π , χ, γ dΩ + Γ Π n; [[u]], [[β]], [[w]] dΓ (4)
where the support functions are defined as :

π , χ, γ = sup (N ,M ,V )∈G N : + M : χ + V • γ (5) Π n; [[u]], [[β]], [[w]] = sup (N ,M ,V )∈G (N • n) • [[u]] + (M • n) • [[β]] + (V • n)[[w]] (6)
From these expressions, it is worth noting that the support function Π of velocity discontinuities can be obtained from the support function π of strain rates through the following relation :

Π n; [[u]], [[β]], [[w]] = π [[u]] s ⊗ n, [[β]] s ⊗ n, [[w]]n (7)
where

a s ⊗ b = 1 2 (a ⊗ b + b ⊗ a) is the symmetrized tensorial product.
The kinematic characterization of the domain K of potentially safe loads is then given by :

Q ∈ K ⇐⇒ ∀(U , Θ) K.A. with q, P (e) (U , Θ) = Q • q ≤ P rm (U , Θ) (8)

General comments

In the previous formulation, it has been implicitly assumed that the shell presents an infinite strength towards the drilling rotation about the shell normal ν. Hence, the drilling strain and discontinuities associated to the out-of-plane rotation ϑ should always be zero.

Despite the generality of the approach proposed here, the remainder of this article will be restricted to the case of thin shells, i.e. presenting an infinite strength towards shear forces compared to membrane forces and bending moments. As intensively discussed in (Bleyer et al., 2015a), this assumption imposes that γ = 0 and [[w]] = 0 as kinematic constraints, reducing to the classical Love-Kirchhoff conditions in the case of a plate.

Finally, in the previous presentation, the (rather complex) expressions of local equilibrium equations and compatibility conditions between generalized strain rates and velocity fields have not been specified. Indeed, these expressions involve differential geometry quantities of the shell and induce a coupling between membrane on the one hand and bending/shear on the other hand. In the next Section, a discretization of the shell into planar facets will be considered so that membrane and bending effects are decoupled in each facet.

Lower and upper bound finite elements

Discretized approximation of the shell geometry

In the following, the shell will be discretized as an assembly of N E triangular planar facets (or plates). As a result of this discretization, the decoupling between membrane and bending equilibrium and strain compatibility equations in each facet makes it possible to use a combination of membrane and plate bending finite elements which have been developed in [START_REF] Bleyer | Lower bound static approach for the yield design of thick plates[END_REF]) (lower bound) and (Bleyer et al., 2015a) (upper bound). The practical implementation being very similar to the one described in these papers, only particular aspects related to the shell model will be highlighted for the sake of concision. While membrane and bending equilibrium and strain compatibility equations are decoupled in each planar facet, it is important to keep in mind that the membrane forces and bending moments will still be coupled by the generalized strength criterion in which the membrane strength is affected by the intensity of the bending moment and vice-versa, as is the case with beam section interaction diagrams.

Lower bound finite element for the static approach

The numerical resolution of the static approach will be performed by considering piecewise continuous membrane forces, bending moment and shear force fields, separated by statically admissible discontinuities across adjoining elements. For the membrane part, a linear interpolation of the membrane force tensor N is assumed in plate each element. As regards the plate bending part, the lower bound finite elements developed in [START_REF] Bleyer | Lower bound static approach for the yield design of thick plates[END_REF] is considered. This element assumes a quadratic variation of the bending moment tensor M and a linear variation of the shear force vector V .

The local equilibrium equations inside a plane facet read as (in the absence of distributed bending couples) :

div N + p = 0 (9) div V + q = 0 (10) div M + V = 0 ( 11
)
where p and q represent an in-plane and out-of-plane distributed force respectively. These equations can be exactly enforced in each element when considering piecewise constant distributed forces in each element (as in [START_REF] Bleyer | Lower bound static approach for the yield design of thick plates[END_REF], equation ( 11) is enforced at the three vertices of each triangular facet).

As regards jump equations, let us consider two adjacent facets (generally non coplanar) sharing a common edge Γ of unit tangent vector t (Figure 2). The normal to the two facets will be denoted by ν + and ν -and the following in-plane normals to the interface Γ are introduced : n ± = ν ± ∧ t. In the case when there are no external line forces applied along Γ, the jump equations for resulting forces and bending moments are given by :

[[R]] = R + -R -= 0 (12) [[M • n]] = M + • n + -M -• n -= 0 ( 13
)
where

R ± = N ± • n ± + (V ± • n ± )ν ± .
From the resulting forces balance equation ( 12), it can be seen that bending/membrane coupling occurs at the edge between two non coplanar adjacent facets. In the case when these facets are coplanar, the jump equations in terms of in-plane and out-of-plane reactions are decoupled, reducing to :

[[N • n]] = 0 and [[V • n]] = 0.
As in [START_REF] Bleyer | Lower bound static approach for the yield design of thick plates[END_REF], due to the linear variation of N and V , it is sufficient to enforce (12) at the two end nodes of a given edge, whereas, due to the quadratic variation of M , ( 13) is enforced at both end nodes plus the mid-side node of the edge. Finally, static boundary conditions are imposed by prescribing the value of the components of R and M • n on the appropriate boundary edges.

All these elementary contributions of local equilibrium, jump equations, and boundary conditions are then assembled into a global linear relationship between the vector of static unknowns Σ = . . . , N e ij , M e ij , V e i , . . . T and a generalized external force vector F of the form HΣ + F = 0. The discretized fields verifying this equation are then statically admissible with the loading corresponding to F. Assuming that the loading mode depends upon one parameter only (the global force vector being proportional to a reference load F 0 : F = λF 0 ), the formulation of the yield design static approach amounts to solving the following maximization problem :

λ l = max λ,Σ λ s.t. HΣ + λF 0 = 0 (N c ij , M c ij ) ∈ G ∀c = 1, . . . , N c • N E (14)
where the superscript c spans all the points where the strength criterion has to be verified as discussed in [START_REF] Bleyer | Lower bound static approach for the yield design of thick plates[END_REF] (N c denotes the number of check points used inside each element). Keeping in mind the remarks made in [START_REF] Bleyer | Lower bound static approach for the yield design of thick plates[END_REF]) about a potential violation of the strength criterion between these points, the solution of problem ( 14) yields a lower bound λ l ≤ λ + to the exact ultimate load factor λ + that the shell can sustain on account of adopted strength criterion.

Obviously, the resolution of problem ( 14) is highly dependent on the form of the strength criterion constraint expressed as (N ij , M ij ) ∈ G. For example, if this constraint can be reformulated using linear equality/inequality constraints only, the corresponding problem will belong to the class of linear programming (LP) problems. Recently, important progresses have been achieved in the development of efficient solvers for a more general class of optimization problems, namely second-order cone programming problems. Most strength criteria involved in yield design/limit analysis problems can be reformulated using second-order cone constraints, so that the corresponding problems can be dealt with using such efficient solvers. The case of shell strength criteria has not received much attention so far and this aspect will be the purpose of Section 4.

Upper bound finite element for the kinematic approach

As regards the upper bound kinematic approach, the w6-d plate bending finite element presented in (Bleyer et al., 2015a) is combined with an upper bound membrane element. Hence, a quadratic interpolation of the virtual velocity U and a linear interpolation of the virtual rotation velocity Θ are considered in each element. As in (Bleyer et al., 2015a), the considered interpolation may be discontinous from one element to another, the degrees of freedom being attached to the element itself and not to a geometrical node shared by several adjacent elements. Including potential discontinuities of the virtual velocity fields in upper bound kinematic approaches have already proved to be not only more efficient in terms of accuracy of the produced estimates, but also to avoid in a natural way the shear locking phenomenon in the thin plate limit (Bleyer et al., 2015a). For these reasons, this strategy has been retained here.

As regards virtual strain rates, membrane strains , bending curvatures χ, shear strains γ and drilling strains associated with the rotation around the facet normal are computed in each element using the following strain compatibility relations in a planar facet :

= 1 2 (∇u + T ∇u) (15) χ = 1 2 (∇β + T ∇β) (16) γ = ∇w -β (17) = a 2 • 1 2 (∇u -T ∇u) • a 1 -ϑ (18)
where ∇ denotes the gradient operator associated with the plane of the facet and a 1 • a 2 = 0. As mentioned before, in the thin shell case, the Love-Kirchhoff kinematics is prescribed inside each facet as an additional constraint by imposing γ = 0, that is β = ∇w, whereas the infinite strength with respect to the drilling rotation imposes = 0.

With the retained degrees of interpolation, the membrane strain rate varies linearly inside each element, whereas the bending curvature χ is uniform per element. The contribution of the strain term to the maximum resisting work (first integral in ( 4)) is performed as described in (Bleyer et al., 2015a).

As regards the contribution of virtual velocity discontinuities to the maximum resisting work (second integral in ( 4)) through each edge j = 1, . . . , N D of the mesh, the discontinuities of U are computed at both end nodes plus the mid-side node wheareas, due to its linear variation, the discontinuities of Θ are computed only at both end nodes. It has already been pointed out that the thin shell hypothesis requires the out-of-plane virtual velocity w to be continuous for the upper bound approach to be relevant. Hence, an average normal ν is defined along each edge shared by two non-coplanar facets (Figure 2). The out-of-plane velocity jump is then defined as

[[w]] = (U + -U -)• ν and is imposed to be zero. The same is done for the out-of-plane rotation [[ϑ]] = (Θ + -Θ -)• ν = 0.
Thus, considering a loading mode with one parameter only for a reference loading F 0 , the implementation of the upper bound kinematic approach amounts to solving the following minimization problem on the generalized vector of kinematic unknowns V = . . . , U e i , Θ e i , . . . T :

λ u = min V N E e=1 g
ω e,g π( e,g , χ e,g ) + ( 19)

N D j=1 g ω j,g Π(n j ; [[u]] j,g , [[β]] j,g ) s.t. F 0 T • V = 1 (20) B V = (21) B χ V = χ (22) B γ V = 0 (23) B V = 0 (24) D u V = [[u]] (25) D β V = [[β]] (26) D w V = 0 (27) D ϑ V = 0 (28)
where ω e,g (resp. ω j,g ) are integration weights for the contribution of virtual strain rates (resp. velocity discontinuities) to the maximum resisting work. Constraint (20) represents a normalization of the work of external loads. Constraints ( 21)-( 24) are the discretized form of the compatibility equations ( 15)-( 18) whereas constraints ( 25)-( 28) represent the computation of different components of velocity jumps across edges of the mesh.

In the case when the quadrature rules used to compute the discretized maximum resisting work ensure an approximation by excess of the true maximum resisting work, the solution of problem ( 19) yield an upper bound λ u ≥ λ + to the exact ultimate load factor λ + . Concerning the way of preserving the upper bound status, more details can be found in (Bleyer et al., 2015a). As before, the resolution of problem ( 19) is highly dependent on the form of the support functions, which itself depend on the form of the generalized strength criterion G. If G can be expressed using second-order cone constraints, then it will also be the case for its support functions. The corresponding optimization problem will then belong to the class of SOCP problems and can be solved with dedicated solvers.

Formulation of shell generalized strength criteria

Approximate shell strength criteria for a von Mises plate

Obtaining an exact analytical expression of a generalized strength criterion for a plate or shell in membrane-bending interaction is quite a difficult task. In the case of a homogeneous shell made of a von Mises material, Ilyushin [START_REF] Ilyushin | Plasticité[END_REF]) obtained a parametrisation of the exact 6-dimensional hypersurface using 3 parameters only. Unfortunately, the expressions involved in this parametrisation are very complex and cannot be formulated using conic constraints, thus precluding the use of efficient conic programming solvers to tackle the associated optimization problems ( 14) and ( 19). Various approximate expressions have then been proposed by several authors in order to derive analytical approximate limit analysis solutions or to obtain more tractable expressions from a numerical point of view. We will review some of these approximate strength criteria and discuss their representation, as well as the expression of their corresponding support function, using conic constraints.

Considering a homogeneous shell of thickness t made of a von Mises material of uniaxial traction/compression strength σ 0 , the following notations are introduced :

n = 1 σ 0 t N 11 N 22 N 12 T = 1 σ 0 t N (29) m = 4 σ 0 t 2 M 11 M 22 M 12 T = 4 σ 0 t 2 M (30) q n = n T Qn (31) q m = m T Qm (32) q nm = m T Qn = n T Qm (33) Q =   1 -1/2 0 -1/2 1 0 0 0 3   (34)
One can immediately observe that the von Mises membrane strength criterion (M = 0) is simply given by √ q n ≤ 1, whereas the von Mises bending strength criterion (N = 0) is given by √ q m ≤ 1. As mentioned before, several authors have proposed approximate expressions for the generalized strength criteria of a von Mises shell :

• Ilyushin [START_REF] Ilyushin | Plasticité[END_REF] :

q n + q m + 1 √ 3 |q nm | ≤ 1 (35)
• Hodge [START_REF] Hodge | Plastic analysis of structures[END_REF] :

√ q m + q n ≤ 1 (36)
• Prager ("sandwich" approximation) [START_REF] Prager | Problems in Continuum Mechanics[END_REF] :

q n + q m + 2|q nm | ≤ 1 (37)
• Robinson [START_REF] Robinson | A comparison of yield surfaces for thin shells[END_REF] :

√ q m + q n ≤ 1 (38)
It is interesting to note that all these approximations are exact for pure membrane (i.e. M ij = 0, hence q m = q nm = 0) sollicitations as well as for pure bending (i.e. N ij = 0, hence q n = q nm = 0) sollicitations, being equivalent, respectively, to the membrane strength criterion √ q n ≤ 1 and the bending strength criterion √ q m ≤ 1.

As a result, all of these approximations only differ in the way they describe the membrane-bending interaction. For a more complete presentation of these approximations, we refer to [START_REF] Robinson | A comparison of yield surfaces for thin shells[END_REF]. In this paper, the maximum distance of these approximate surfaces to the exact one is also computed. Results are reported in Table 1.

Let us also mention that it is possible to consider other very crude approximations such as :

• the "square"-shaped strength criterion which assumes no interaction between membrane and bending. It consists in the cartesian product of the pure membrane and the pure bending strength criteria : [START_REF] Robinson | A comparison of yield surfaces for thin shells[END_REF]))

√ q n ≤ 1 √ q m ≤ 1 (39)
• the "diamond"-shaped strength criterion which assumes a linear interaction between membrane and bending. It consists in the convex hull of the pure membrane and the pure bending strength criteria :

√ q n + √ q m ≤ 1 (40)
• the "circle"-shaped strength criterion which assumes a non-linear interaction between membrane and bending. It consists in a 6-dimensional ellipsoid generated by the pure membrane and the pure bending strength criteria and is equivalent to (38).

These three approximations rely only on the pure membrane and bending strength criteria and can thus easily be generalized to other constitutive materials by replacing √ q n (resp. √ q m ) by the corresponding expression of the pure membrane (resp. bending) criterion expressed, for instance, as g memb (N ij ) ≤ 1 (resp. g bend (M ij ) ≤ 1).

All of the approximations mentioned so far can be expressed using conic constraints as shown in Appendix A, where expressions of the corresponding support function are also given. Hence, it is possible to use any of these approximate strength criteria in conjunction with the lower and upper bound finite elements presented before and conic programming solvers. As an illustration, all these different strength criteria are represented in the plane (N 11 , M 11 ) (Figure 3a) and (N 22 , M 11 ) (Figure 3b).

Extension to other material strength criteria

In the case when the shell constitutive material obeys another strength criterion than von Mises (e.g. Tresca), other approximations are necessary. For example, it is possible to assimilate the homogeneous shell to two lower and upper membranes ("sandwich" approach) as suggested in [START_REF] Hodge | The rigid-plastic analysis of symmetrically loaded cylindrical shells[END_REF][START_REF] Lubliner | Plasticity Theory[END_REF]. It is also possible to use the previous approximations obtained in the von Mises case since the Tresca strength criterion can be circumscribed by the von Mises ellipsoid and inscribed by the same ellipsoid reduced by a factor √ 3/2. Finally, as mentioned before, the "square", "diamond" and "circle" criteria can be obtained from the pure membrane and bending strength criteria. However, these different strategies produce, in general, approximations of a relatively poor quality, especially for [START_REF] Prager | Problems in Continuum Mechanics[END_REF],

"diamond" "circle"

"square"

(a) (N11, M11)-plane "diamond" "circle" [START_REF] Prager | Problems in Continuum Mechanics[END_REF] "square" complex constitutive strength criteria or in the case when the constituent material is heterogeneous through the shell thickness (multilayered shell for example).

(b) (N22, M11)-plane

Formulation of generalized strength criteria using an

implicit up-scaling procedure In order to overcome the previously mentioned drawbacks and to obtain approximate strength criteria with a desired level of accuracy, we now propose a general approach which will be well suited for a numerical implementation using conic programming tools.

If the shell is made of a constituent material with homogeneous strength properties in the in-plane directions of the shell (case of multilayered shells), then it can be shown that the generalized membrane-bending strength criterion G can be expressed from the local constitutive strength criterion G(ξ) by considering a distribution of plane stress states depending on the thickness coordinate ξ only [START_REF] Save | Plastic limit analysis of plates, shells, and disks[END_REF][START_REF] Dallot | Limit analysis of multi-layered plates. Part I: the homogenized Love-Kirchhoff model[END_REF]

1 : (N , M ) ∈ G ⇔                  ∀ξ ∈ [-t/2; t/2] and i, j = 1, 2 ∃ σ(ξ) = σ ij (ξ)a i ⊗ a j ∈ G(ξ) N ij = t/2 -t/2 σ ij (ξ)dξ M ij = t/2 -t/2 (-ξ)σ ij (ξ)dξ (41)
This formulation results from an upscaling approach and makes it possible to express the generalized strength criterion as a function of the local constitutive strength criterion G, although in an implicit way through the solution of the auxiliary yield design problem (41).

It is also possible to obtain a similar determination of the generalized support function π( , χ). Indeed, recalling the definition

π( , χ) = sup (N ,M )∈G N : + M : χ (42)
and using (41), we obtain :

π( , χ) = sup σ(ξ)∈G(ξ) t/2 -t/2 σ(ξ)dξ : + t/2 -t/2 (-ξ)σ(ξ)dξ : χ = sup σ(ξ)∈G(ξ) t/2 -t/2 σ(ξ) : ( -ξχ)dξ = t/2 -t/2
π ps loc ( -ξχ)dξ (43)

1 Strictly speaking, (41) represents a lower bound approximation to G making use of plane stress distributions. Actually, it may be proved (from using for instance the upper bound kinematic approach) that it leads to the exact determination of G.

where π ps loc (d) corresponds to the support function of the local plane stress strength criterion. The generalized support function is then obtained by integrating through the thickness the local support function for a virtual strain rate d = -ξχ corresponding to a Love-Kirchhoff kinematics.

Numerical lower bound approximation

Nevertheless, the continuous formulation (41) cannot be used as such for a numerical implementation. Hence, a lower bound approximation of the true generalized strength criterion G will first be proposed. To do so, a piecewise constant distribution of plane stress states in n different layers [t k-1 , t k ] with k = 1, . . . , n and t 0 = -t/2, t n = t/2 is considered. Assuming that the local strength criterion G(ξ) = G k is uniform in each layer k (multilayered approximation), we define the following lower bound approximation G lb of G :

(N , M ) ∈ G lb ⊆ G ⇔                        ∃ σ k ij such that k = 1, . . . , n and i, j = 1, 2 σ k ij a i ⊗ a j ∈ G k N ij = n k=1 (t k -t k-1 )σ k ij M ij = n k=1 t 2 k-1 -t 2 k 2 σ k ij
(44) It can be immediately seen that this definition leads, for any value of n, to a lower bound estimate G lb ⊆ G in the case of a piecewise constant local strength criterion. Hence, G lb will converge to G, from the inside, at the limit n → ∞. From a numerical point of view, G lb can be determined by introducing 3n auxiliary optimization variables σ k ij . Besides, due to the linear relationship between the auxiliary variables σ k ij and the membrane forces N ij and bending moments M ij , G lb can be expressed using conic constraints if the local strength criterion G k can be expressed using such constraints, which is the case for an important number of classical plane stress criteria e.g. Rankine, von Mises, Tresca, Mohr-Coulomb, Drucker-Prager, etc...

Example of a homogeneous von Mises shell :

As an illustration, let us consider again the case of the homogeneous von Mises shell. Introducing

σ k = σ k 11 σ k 22 σ k 12
T for k = 1, . . . , n, the local plane stress strength criterion may be written as :

σ k ij a i ⊗ a j ∈ G k ⇐⇒ σ k T Qσ k ≤ σ 0 ⇐⇒ P T σ k ≤ σ 0 with P = 1 √ 3   1 0 0 -1/2 √ 3/2 0 0 0 √ 3   ( 45
)
where P is the Cholesky factor of Q = PP T , the local strength criterion P T σ k ≤ σ 0 being here written in the form of a second-order conic constraint.

In the case of a crude approximation using n = 2 layers of thickness t/2, G lb can be expressed as follows :

(N , M ) ∈ G lb ⇔              P T σ 1 ≤ σ 0 P T σ 2 ≤ σ 0 N = t 2 (σ 1 + σ 2 ) M = t 2 8 (σ 2 -σ 1 ) (46)
leading, after some computations, to :

q n + q m -2q nm = 1 σ 0 P T σ 1 ≤ 1 (47) q n + q m + 2q nm = 1 σ 0 P T σ 2 ≤ 1 (48)
Therefore, we see that for n = 2, (N , M ) ∈ G lb ⇔ q n + q m + 2|q nm | ≤ 1. As a conclusion, the lower bound multilayered approximation using 2 layers of equal thickness for a von Mises material is equivalent to approximation (37). This should not be a surprise since approximation (37) has been obtained using the "sandwich" approach mentioned before. However, as it will be seen later, increasing the number n of layers will produce more accurate approximations than those mentionned earlier.

Numerical upper bound approximation

The previous multilayered strategy makes it possible to construct lower bound approximations for the generalized strength criteria, which can then be used in conjunction with a lower bound finite element yield design of the shell, so as to keep the lower bound status of the so-computed ultimate load estimate. It would then be interesting to obtain an equivalent upper bound approximation to G, to be used in conjunction with upper bound finite elements for the same reasons.

To do so, let us recall expression (43) which relates the generalized support function to that of the local strength criterion :

π( , χ) = t/2 -t/2 π ps loc ( -ξχ)dξ (49)
Since π ps loc (d) is a convex function of its argument, which is itself an affine function of the thickness variable (d = -ξχ), the previous integral can then be estimated by excess using a trapezoidal quadrature rule with n ≥ 2 points :

π( , χ) ≤ π ub ( , χ) = n k=1 ω k π ps loc ( -ξ k χ) (50)
where ξ k are the quadrature points with ξ 1 = -t/2 and ξ n = t/2 and ω k = (ξ k+1 -ξ k-1 )/2 for 2 ≤ k ≤ n -1, ω 1 = (ξ 2 -ξ 1 )/2 and ω n = (ξ n -ξ n-1 )/2. For the previous estimate to be a true upper bound, the local strength criterion must be uniform across the shell thickness. Nevertheless, in the case of a piecewise constant strength criterion, an upper bound estimate can be obtained by splitting the integral on the different domains where the strength criterion is uniform. The trapezoidal quadrature is then adapted to take into account the discontinuities between the different domains.

The convex set defined as follows :

G ub = (N , M ) s.t. N : + M : χ ≤ π ub ( , χ) ∀ , χ (51)
is clearly an upper bound approximation to G for all values of n. G ub will converge to G, from the outside, at the limit n → ∞. Similarly, if the local strength criterion G can be expressed using conic constraints, it will also be the case for its support function by a duality argument. Therefore, the upper bound generalized support function π ub ( , χ) can also be expressed using such constraints and this approximation can efficiently be used in conjunction with conic programming solvers. Using duality, it can also be easily shown that G ub can be expressed in a way similar to (44) as follows :

(N , M ) ∈ G ub ⊇ G ⇔                        ∃ σ k ij such that k = 1, . . . , n and i, j = 1, 2 σ k ij a i ⊗ a j ∈ G k N ij = n k=1 ω k σ k ij M ij = n k=1 (-ω k ξ k )σ k ij
(52) This expression also involves 3n auxiliary optimization variables but, contrary to (44), G ub cannot be interpreted as a discretization of the shell thickness into n nonoverlapping layers with uniform plane stress states in each layer.

Other possible approximations

The two previous strategies enabled to build a lower and an upper bound approximation of G respectively. This bounding property is indispensable to preserve the bounding status of a numerical ultimate load estimate for the shell when used with lower or upper bound finite elements. In some cases, keeping the bounding status of the approximation may not be necessary, e.g. when using mixed finite elements for which the bounding status of the computed estimate is lost, even when using an exact strength criterion. For this reason, it is also possible to obtain nonbounding approximations of the exact generalized strength criterion G. A simple idea then consists in applying the same strategy as for the upper bound approximation 4.3.2 but using a quadrature rule of integral (43) which does not necessarily yield an estimate by excess :

π( , χ) ≈ π approx ( , χ) = n k=1 ω k π ps loc ( -ξ k χ) (53)
where ξ k are the quadrature points in [-t/2; t/2] and ω k the quadrature weights. The choice of different quadrature rules (e.g. Simpson's rule, Gauss-Legendre quadrature, ...) will therefore produce different approximations.

As before, the convex set defined as follows :

G approx = (N , M ) s.t. N : + M : χ ≤ π approx ( , χ) ∀ , χ (54) is a non-bounding approximation of G which converges to G as n → ∞.
G approx can also be expressed using implicit auxiliary optimization variables as in ( 52) by simply replacing ω k and ξ k by ω k and ξ k .

Assessing the quality of the different approximations

Two testing cases

To assess the accuracy of the different proposed approximations, two testing cases of different complexity will be investigated.

The first one corresponds to a homogeneous von Mises shell, for which existing approximations have already been discussed in 4.1. The local strength criterion being isotropic, it can directly be concluded that the macroscopic strength criterion G is also isotropic. This property will also be true for all the considered approximations. Other remarkable properties can be deduced such as :

• G is symmetric with respect to the origin : (N , M ) ∈ G ⇔ -(N , M ) ∈ G. This is easily proven using (41) and noting that the von Mises criterion itself is symmetric with respect to the origin.

• G is symmetric with respect to the membrane subspace : (N , M ) ∈ G ⇔ (N , -M ) ∈ G. This is easily proven using ( 41) and observing that the shell is symmetric with respect to ξ = 0 (homogeneous shell).

• From the two previous properties, G is also symmetric with respect to the bending subspace :

(N , M ) ∈ G ⇔ (-N , M ) ∈ G.
The first property is also true for all approximations, whereas the second will be true if the discretization is symmetric with respect to the shell mid-surface ξ = 0 (even (resp. odd) value of n for the lower (resp. upper) bound approximation with a uniform discretization).

The second example consists in two layers of thickness t/2 made of a material obeying a generalized orthotropic plane stress Rankine strength criterion (used for example to decribe the strength properties of a concrete). The two orthogonal directions of orthotropy for each layer are denoted by e α = cos αa 1 + sin αa 2 and e β = -sin αa 1 + cos αa 2 . In the following, we assume that α = 0 • in the lower layer [-t/2; 0] and α = 45 • in the upper layer. The orthotropic generalization of the plane stress Rankine strength criterion reads as :

σ ∈ G ⇐⇒            (σ - 0α + σ αα )(σ - 0β + σ ββ ) ≥ σ 2 αβ (σ + 0α -σ αα )(σ + 0β -σ ββ ) ≥ σ 2 αβ -σ - 0α ≤ σ αα ≤ σ + 0α -σ - 0β ≤ σ ββ ≤ σ + 0β ( 55
)
where σ ± 0α,β are the tensile/compressive uniaxial strength in the directions of orthotropy. In the following, we will assume that σ + 0α = 5σ 0 and σ - 0α = σ + 0β = σ - 0β = σ 0 for both layers. The conic formulation of this strength criterion as well as of his support function are given in Appendix B. Contrary to the simple case of the homogeneous von Mises shell, the generalized strength criterion G is not isotropic (not even orthotropic) and does not possess the symmetries with respect to the origin (the local strength criterion being non-symmetric with respect to the origin) or with respect to the pure bending/membrane sollicitations (the shell being non-symmetric with respect to its mid-plane).

Projections on different planes

To investigate the accuracy of the proposed lower and upper bound approximations for the von Mises shell, the corresponding strength criteria have been represented in the (N 11 , M 11 ) plane (Figure 4a) and in the (N 22 , M 11 ) plane (Figure 4b) for different values of the discretization parameter n. Uniform discretization for both lower and upper bound approximations have been chosen. Only even (resp. odd) values of n have been chosen for the lower (resp. upper) bound approximations so that approximations are exact for pure bending sollicitations. One can observe the increasing quality with n of the approximation from the inside and from the outside. As a comparison, approximate criteria obtained with a Gauss-Legendre quadrature have been represented in Fig- ure 5. The accuracy of the approximation seems to be higher for similar values of the discretization parameter n compared to the previous lower and upper bound approximations. But, it can clearly be seen that this approximation has no bounding status and that the pure bending strength, in particular, is never exact.

As regards the bilayered orthotropic Rankine shell, projections of the different approximate strength criteria on the (N 11 , M 11 ) and (N 22 , M 11 ) planes have been reported in Figure 6a and 6b for the lower and upper bound approximations and in Figure 7a and 7b for the Gauss-Legendre approximation. Despite the relative complexity of this example, the proposed approximations are of excellent quality, especially for n ≥ 4 for both the lower and upper bound approximations. Furthermore, in spite of the lack upper bound lower bound exact of symmetry, it is also interesting to observe that both approximations always yield the exact value of the maximum strength in membrane and bending. On the contrary, this is not the case for the Gauss-Legendre quadrature.

Maximum approximation errors

In order to quantify more precisely the relative errors induced by the different approximations, the boundary of the different yield surfaces have been captured by sampling the six-dimensional space R 6 with 10,000 points. The maximum distance between the approximated criteria and the exact one (identified with either of the two approximations using n = 100) has then been computed. Results in the case of the bilayered shell have been represented in Figure 8a for the lower and upper bound approximations ( 44) and ( 50). First, it can be verified that the lower (resp. upper) bound approximation always induces an approximation of the exact strength criterion by default (resp. by excess). Besides, such approximations are of relatively good quality for n = 4, 5 (roughly ±10% maximum error). The convergence is also relatively fast, since it can be estimated to decrease proportionally to n -2 . In Figure 8b, the evolution of the maximum relative error for a Gauss-Legendre quadrature has been represented. The previous results are also recalled in dashed lines. It can be observed that the Gauss-Legendre quadrature strategy always produce by default and by excess errors, which seem to be approximately of the same order. It is also interesting to note that these maximum errors are very close to those produced by the lower and upper bound strategy. Therefore, it can be concluded that the Gauss-Legendre strategy does not necessarily yield better estimates than the lower and upper bound approximations. The latter should, therefore, be preferred, due to their rigorous bounding status and equivalent computational cost.

Illustrative applications to the yield design of shell structures

The last part of this paper is devoted to applying the above described numerical procedure to evaluating the ultimate yield strength of various shell structures, ranging from classical to more complex problems.

Simply supported spherical cap under uniform pressure

The first example considers the problem of a spherical cap shell of radius R and thickness t made of a consitutent material obeying the Tresca strength criterion of uniaxial strength σ 0 . The shell is subjected to a uniformly distributed external radial pressure p and is simply supported on its boundary (see Figure 9). The cap opening angle is denoted by α and the following non-dimensional slenderness ratio is introduced k = t/4R. In the following, estimates for the ultimate pressure will be given in terms of the non-dimensional pressure p * = pR/(σ 0 t). This classical problem has already been considered by various authors who provided different analytical lower and upper bound estimates for the ultimate pressure p . In particular, by considering estimates derived by Hodge in [START_REF] Hodge | Limit analysis of rotationally symmetric plates and shells[END_REF] based on plate and membrane solutions, as well as different approximate generalized strength criteria, the following analytical estimates for the ultimate pressure can be obtained [START_REF] Save | Plastic limit analysis of plates, shells, and disks[END_REF] :

p * LB = max{p * 2 ; 0.618p * 1 } ≤ p * ≤ min{p * 1 ; 1.25p * 3 } = p * U B
(56) where :

p * 1 = 2 + 2k sin α ln 1 + sin α cos α -sin α -1 (57) p * 2 = max 2; 2k (1 + k)(1 -α cot α) (58) p * 3 =      2k 1 -α cot α if cos α ≥ 1 -k 2 (sin α -ϕ cos α -(1 -k)(sin α -sin ϕ)) sin α -α cos α if cos α < 1 -k (59) with cos ϕ = cos α/(1 -k)
Numerical lower (resp. upper) bound estimates for the ultimate load p * have been computed using the static lower (resp. kinematic upper) bound finite elements presented in Section 3, in conjunction with the approximate lower (resp. upper) bound formulation (44) (resp. ( 50)) of the generalized strength criterion. Typical meshes for different values of the opening angle α consisted of roughly 700-800 triangular facets (Figure 9), while the discretization parameter n of the shell thickness has been varied from n = 2 to n = 6. The so-obtained SOCP problems have been solved using the Mosek software package [START_REF] Mosek | The Mosek optimization software[END_REF]. Typical optimization times needed to solve such problems ranged from 10s to less than a minute, for increasing values of n, on a standard desktop computer.

In Figure 10, lower bound estimates using n = 2, 4, 6 and upper bound estimates using n = 2, 3, 5 have been reported along with the analytic lower (p LB ) and upper (p U B ) bounds (56) for a shell with a moderate value (k = 0.1) of the slenderness ratio. First, it can be observed that the ultimate pressure is estimated with an excellent accuracy in the regime of high values of the opening angle (α ≥ 45 • ) for all values of n. This can be attributed to the fact that the solution is close to a membrane state and that the approximate formulations of the generalized strength criterion are exact for such states. Secondly, the quality of the bracketing is improved when increasing the value of the discretization parameter n. It is to be noted that, even for a not so refined mesh and moderate values of n, the quality of the computed estimates is relatively good, the maximum relative gap between the static approach with n = 6 and the kinematic approach with n = 5 being less than 10% for this mesh density. It is to be noted that these estimates fall inside the domain delimited by the analytic lower and upper bounds. Finally, the poor quality of the upper bound estimates obtained with n = 2 can be explained by the fact that the corresponding approximate generalized strength criterion is not exact for pure bending states, which deteriorates the quality of the corresponding estimate. Hence, using a Gauss-Legendre quadrature for the approximate formulation as suggested in (53) would lead to similar results. In Figure 11, lower bound estimates using n = 6 and upper bound estimates using n = 5 have been reported along with the analytic lower (p LB ) and upper (p U B ) bounds of (56) for a shell with a small value k = 0.005 of the slenderness ratio. It can be observed that the quality of the estimate is similar, the maximum relative gap between the two approaches being less than 8%. As an illustration, the numerical estimates obtained with the Ilysuhin approximate strength criterion (35) have also been represented. These bounds have been computed using lower and upper bound finite elements and appropriate coefficients have then been applied to the so-obtained ultimate loads in order to bound the Tresca generalized strength criterion by Ilyushin approximate surface as discussed in Section 4.2. The quality of the so-obtained estimates is relatively bad, with a maximum relative gap of 40%. From a computational point of view, computing times with Ilyushin surface took approximately 15s. Finally, it can be noted that the solution approaches the limit p * = 2 for a value of the opening angle much smaller (α ≈ 15 • ) than in the case k = 0.1 (α ≈ 60 • ). One can observe that the lower bound analytic estimate p * LB seems to be very close to the exact solution in this case.

cap opening angle normalized ultimate load

Clamped cylindrical shell under vertical distributed loading

The second example investigates the ultimate load of a cylindrical shell of radius R, thickness t and length 2L, clamped at both ends x = 0 and x = 2L. This shell is subjected to a uniformly distributed vertical loading -f e z and only a quarter of the shell has been modeled by taking appropriate symmetry boundary conditions into account (Figure 12). In the following, the thickness to radius ratio will be kept constant to t/R = 0.01 while the structural slenderness ratio 2L/R will be varied. The shell constitutive material is assumed to obey the von Mises strength criterion of uniaxial strength σ 0 . In the limit of high slenderness ratios 2L/R 1, this shell problem is expected to be appropriately modeled by a one-dimensional bending beam. The corresponding ultimate load of such a beam model can be computed exactly: the corresponding collapse mechanism consisting of two rigidly rotating regions separated by three hinges located at both clamped extremities x = 0, 2L and at the mid-span x = L. Hence, this ultimate load is given by q + = 16M 0,beam /(2L) 2 where q corresponds to the equivalent lineic load (i.e. q = 2πRf ) while M 0,beam = 4t(R 2 + t 2 12 )σ 0 ≈ 4tR 2 σ 0 is the ultimate bending moment of a beam with a tubular section of radius R and thickness t. Finally, the ultimate load obtained with the bending beam model is given by :

f + beam = 32 π σ 0 t R 2L 2 (60)
Lower and upper bound finite element computations, using respectively n = 6 and n = 5 for the approximate von Mises generalized strength criterion, have been performed on structured meshes with 20 elements along the half-circumference and elements of length R/4 along the x direction. The so-obtained lower and upper bound estimates of the ultimate load f + have been normalized with respect to the beam model ultimate load f + beam and represented in Figure 13 for varying values of the structural slenderness ratio 2L/R. As a comparison, ultimate load estimates obtained with Ilyushin approximate strength criterion, multiplied by appropriate scaling factors (see Table 1) to ensure the bounding status of Ilyushin surface, have been represented. Finally, ultimate load estimates obtained by elastoplastic shell finite element computations using Cast3m [START_REF] Verpeaux | Castem 2000: une approche moderne du calcul des structures[END_REF] software have also been reported.

First, it can be observed that elastoplastic computation estimates fall between the lower bounds of the static approach and the upper bounds of the kinematic approach, thus comforting our results. As expected, an evolution of the ultimate load towards a value close to the beam model for increasing values of the slenderness ratio can clearly be observed. For moderate values of 2L/R, the computed ultimate load is much smaller than the one predicted by the beam model. These aspects are further confirmed when inspecting the collapse mechanisms obtained with the kinematic approach for different slenderness ratios (Figure 14). In spite of the infinite shear strength of the shell itself, the structural collapse mechanism for low values of the slenderness ratio is similar to a global shear failure of the beam. In particular for 2L/R = 5, the circular section of the cylinder exhibits a collapse mode involving three hinge (or yield) lines along the x longitudinal direction (Figure 14a). For high values of the slenderness ratio, a bending beam-like collapse mechanism involving localized deformations near the clamped supports and at mid-span, separated by rigidly rotating regions, is retrieved, recalling the three hinges involded in the one dimensional beam collapse mechanism (Figure 14c-14d). Finally, let us mention that such collapse mechanisms are very similar to those obtained by the last increments of the elastoplastic computations using Cast3m. As regards the bounds obtained with Ilyushin surface, the relative gap is larger (15% in average) than with our approach (8% in average). By using the same strength criterion for both the static and the kinematic approach, it is possible to estimate the contribution of the mesh discretization to this relative gap. It is interesting to note that the relative gap due to the mesh discretization only is about 4%, the remainder (4% for our approach, 11% for Ilyushin surface) being due to the distance between the lower and upper bound approximate strength criteria.

Hyperbolic paraboloid shell with orthotropic strength criterion

This last short example aims at illustrating the capability of the proposed approach to deal with more complex local strength criteria. A hyperbolic parabolic shell of equation z = -4hxy/L 2 for x, y ∈ [-L/2; L/2] (with L = 10 and h = 2 in the following) is clamped on four square supports of width c = 0.4 at its four corners, the whole surface being subjected to a uniformly distributed vertical load -f e z while the shell thickness is t = 0.1 (see Figure 15). The shell constitutive material is supposed to obey the orthotropic generalized Rankine strength criterion ( 55) with e α = cos αe x + sin αe y and e β = -sin αe x + cos αe y as directions of orthotropy and σ + 0α = σ + 0β = σ - 0β = σ 0 and σ - 0α = 100σ 0 (direction e α being, thus, highly resistant in compression).

Collapse mechanisms for different values of the maximum resisting direction angle α have been represented in Figure 16. First, it is to be observed that, for α = 0 • , the geometrical reflection symmetry with respect to both diagonals y = x and y = -x is broken by the orthotropy of the strength criterion. The collapse mechanism reduces here to a yield line along y = 0 due to a weakest bending strength in the e y = e β direction. One can observe exactly the same for α = 90 • along x = 0. For α = 45 • and α = 135 • the diagonals symmetry is restored since the orthotropy directions are aligned with the latter. However, a slight difference between both mechanisms can be observed. This is due to the difference between tensile and compressive strength when the strongest direction is aligned either with the first diagonal x = y, mainly sollicited in compression, or with the second one y = -x, mainly sollicited in traction. The same differences can be observed between α = 60 • and α = 120 • . Finally, all of these remarks are corroborated by the rather erratic evolution of the so-computed ultimate load estimates represented in Figure 17.

Conclusions and future work

A numerical procedure for performing lower and upper bound yield design analysis of shells with general strength criteria has been presented. Using a faceted discretization of the shell surface, lower and upper bound plate finite elements are combined to membrane elements in a simple manner. Apart from the implementation of lower bound equilibrium and upper bound kinematic shell finite elements, the main novelty of this paper relies in the formulation of generalized shell strength criteria in membrane-bending interaction. Implicit through-thethickness formulations are proposed to derive lower and upper bound approximations of the exact generalized strength criterion for any local constitutive material. Such a formulation is particularly suited for a numerical implementation of the lower and upper bound yield design shell problems using second-order cone programming. Different illustrative applications have been considered to illustrate the potentials of the proposed approach.

From such numerical results, it can be concluded that the proposed procedure combines the following aspects :

• accuracy : the ultimate load estimates enable to bound the exact ultimate load with a relatively good accuracy, which can be controled by mesh refinement as well as by increasing the discretization parameter n. Despite the fact that the optimization problem size scales with n, numerical tests showed that the computational cost remained reasonable (e.g. for n = 6, depending on the local strength criterion, between 2 to 4 times more expensive than when using Ilyushin approximate surface).

• generality : the proposed approximate formulation of the generalized strength criteria combined to conic programming enables to treat a wide range of problems with different local strength criteria, provided they can be formulated using conic constraints (Rankine, Tresca, von Mises, Drucker-Prager, etc.).

In particular, it can tackle anisotropy as well as through-the-thickness heterogeneity of the local constitutive material (Bleyer et al., 2015b). Besides, although this paper has only considered shells with infinite shear strength, a finite resistance to shear be easily incorporated in the model (see (Bleyer et al., 2015a;[START_REF] Bleyer | Lower bound static approach for the yield design of thick plates[END_REF] in the case of homogeneous plates).

• bounding status : the proposed methodology ensures, at each of its steps (finite element discretization and strength criteria approximation), that the bounding status of the static and kinematic yield design approaches is preserved2 . This particular aspect allows to obtain reliable estimates of the ultimate load.

Perspectives on the basis of this work are quite numerous. First, such shell finite elements could be used for shakedown analysis in the case of cyclic loading for pipes, tanks and vessels. In the case of reinforced concrete shell structures, the proposed methodology can also be very easily extended to take into account the presence of reinforcing steel bars (Bleyer et al., 2015b). It would also be interesting to investigate the performance of the proposed method on multilayered structures, in which shear effects may become particularly important. Finally, the development of curved shell finite element would certainly increase the accuracy of the numerical estimates with respect to the mesh size. When interpreting results from yield design/limit analysis on slender structures, one should keep in mind that the influence of geometrical non-linearities is not taken into account. For such structures, local or global buckling can occur at a lower load than the ultimate load predicted by a yield design analysis and will therefore influence the design of the structure. In some other cases, geometrical changes can occur slightly before or during collapse and can possibly reduce the computed ultimate load. Previous works [START_REF] Save | Plastic limit analysis of plates, shells, and disks[END_REF]Corradi andPanzeri, 2003, 2004;[START_REF] Raithatha | Rigid plastic model of incremental sheet deformation using second-order cone programming[END_REF]Bleyer et al., 2015b) have already attempted to provide a means for assessing the influence of geometrical changes but this question still need further investigation in order to be included in applications on real structures.

  Figure 1: Shell kinematics

Figure 2 :

 2 Figure 2: Two adjacent triangular plates of the discretized shell

Figure 3 :

 3 Figure 3: Representation of exact and approximate generalized strength criteria for a von Mises shell

Figure 4 :

 4 Figure 4: Representation of lower and upper bound approximate strength criteria for a homogeneous von Mises shell

Figure 5 :

 5 Figure 5: Representation of approximate strength criteria obtained with a Gauss-Legendre quadrature for a von Mises shell

Figure 6 :

 6 Figure 6: Representation of lower and upper bound approximate strength criteria for the bilayered Rankine shell. The values of n corresponding to different line styles are the same as those of Figure 4.

Figure 7 :Figure 8 :

 78 Figure 7: Representation of approximate strength criteria obtained with a Gauss-Legendre quadrature for the bilayered shell. The values of n corresponding to different line styles are the same as in Figure 5.

Figure 9 :

 9 Figure 9: Simply supported spherical cap problem and typical mesh

Figure 10 :

 10 Figure 10: Ultimate load estimates for the simply supported spherical cap in the case k = 0.1 for different discretization parameters

Figure 11 :

 11 Figure 11: Ultimate load estimates for the simply supported spherical cap in the case k = 0.005

Figure 12 :

 12 Figure 12: Clamped cylinder problem and typical mesh

Figure 13 :

 13 Figure 13: Normalized ultimate load estimates for the clamped cylinder problem as a function of the slenderness ratio, using present lower and upper bound approximate formulation of the generalized strength criterion, Ilysuhin approximate criterion and elastoplastic (EP) computations.

Figure 14 :Figure 15 :

 1415 Figure 14: Collapse mechanism and local dissipation for different slenderness ratios

Figure 17 :

 17 Figure 17: Evolution of the normalized ultimate load estimate f * = f + /(σ 0 t 2 4 ) as a function of the orthotropy angle α

Table 1 :

 1 Maximum under-and over-estimation relative errors of different approximate generalized strength criteria with respect to the exact criterion of a homogeneous von Mises shell (from

	Approximation under-estimation over-estimation
	(35)	7.1%	3.4%
	(36)	16.7%	0%
	(37)	20%	0%
	(38)	4.5%	15.5%

at least with respect to the faceted geometry. Indeed, the facetization of the true shell geometry introduces inevitably a geometrical discretization error which has no bounding status.

Acknowledgments

The authors would like to thank the Mosek software package team for their free release policy for academic re-

search.

Appendix A. Conic formulation of generalized strength criteria for shells

This appendix aims at formulating different generalized strength criteria, as well as their support functions, using second-order conic constraints.

Appendix A.1. Second-order cone programming Second-order cone programming (SOCP) is a particular class of convex programming which consists in the minimization/maximization of a linear function, under linear equality and/or inequality constraints and particular non-linear constraints, namely second-order cone (SOC) constraints. A SOCP is, therefore, a convex optimization problem which can be cast under the following form [START_REF] Boyd | Convex optimization[END_REF] :

and each K i is one of the following second-order cone of dimension q i (with

. Pure membrane and bending von Mises strength criteria

Recalling matrix P introduced in (45) as the Cholesky factor of Q = PP T , the pure membrane von Mises strength criterion can be written as :

It is, thus, expressed as a SOC constraint on a Lorentz cone of dimension 4 : (1, P T n) ∈ L 4 or, equivalently, (σ 0 t, P T N) ∈ L 4 . It follows immediately that the pure bending criterion is given by (1, P T m) ∈ L 4 or, equivalently, (σ 0

Appendix A.3. "Square","circle" and "diamond" interaction criteria

Using the previous results, the "square" strength criterion (39) is given by :

The "circle" strength criterion (38) can be rewritten as :

where Σ = n m T . Hence, the "circle" strength criterion is expressed as a SOC constraint on a Lorentz cone of dimension 7 : (1, P T Σ) ∈ L 7 .

The "diamond" strength criterion ( 40) is equivalent to :

Hence, the "diamond" strength criterion is expressed as a two SOC constraints on Lorentz cones of dimension 4, (s 1 , P T n) ∈ L 4 and (s 2 , P T m) ∈ L 4 , plus one linear inequality constraint s 1 + s 2 ≤ 1.

Appendix A.4. "Square","circle" and "diamond" support functions

First, let us recall the following classical results of convex analysis [START_REF] Boyd | Convex optimization[END_REF]) :

• the support function of a convex set defined by E = {x s.t. Ax ≤ 1} (equation of an ellipsoid) is given by :

• the support function of the convex hull CH(A; B) of two convex sets A and B (with support functions π A and π B ) is given by :

• the support function of the intersection A ∩ B of two convex sets A and B (with support functions π A and π B ) is given by :

The support function of the "square" strength criterion follows directly from the absence of interaction between membrane and bending :

2 xy T and χ = χ xx χ yy 2χ xy T . In order to obtain a conic formulation of the corresponding kinematic approach, the expression of the support function has to be reformulated using conic constraints. In the case of the "square" criterion, this reformulation reads as :

The support function of the "circle" strength criterion follows from expression (A.4) and result (A.6) : .11) where d = χ T . Its conic reformulation is immediate :

Finally, as regards the "diamond" strength criterion, one has to recall that it is the convex hull of the pure membrane strength criterion : (n, m) such that P T n ≤ 1 and m = 0 and the pure bending strength criterion : (n, m) such that n = 0 and P T m ≤ 1. Using result (A.7), its support function is given by :

Ilyushin and Prager strength criteria and support functions

Let us now consider the Ilysuhin (35) and the Prager ("sandwich") (37) strength criteria written as :

where the Ilyushin strength criterion is obtained for α = 1 2 √ 3 and the Prager "sandwich" criterion for α = 1. The Cholesky factorization of matrix Q α = P α P T α is given by :

The conic formulation of these criteria is, therefore, given by :

corresponding to the intersection of two Lorentz cone constraints of dimension 7 (this dimension reduces to 4 for the Prager criterion with α = 1 due to the null second column of P α ). Now, using result (A.8), the support function of such criteria is given by :

where

Hodge strength criterion and support function

Hodge strength criterion (36) can be rewritten as follows :

Hence, Hodge strength criterion can be expressed as a two SOC constraints on a Lorentz cone of dimension 4, (s 1 , P T m) ∈ L 4 , and on a rotated Lorentz cone of dimension 5, (s 2 , s 3 , P T n) ∈ L 5 r , plus one linear inequality constraint s 1 +s 2 ≤ 1 and one equality constraint s 3 = 1/2.

After some computations, it is possible to show that the support function of the Hodge strength criterion can be expressed as :

Appendix B. Conic formulation of the generalized Rankine strength criterion

We recall expression (55) :

Introducing the following auxiliary variables :

) is therefore equivalent to :

This local strength criterion is, therefore, expressed using two rotated Lorentz cone constraints of dimension 3.

Its support function can be readily computed as :