Jeremy Bleyer
email: jeremy.bleyer@enpc.fr

Patrick De Buhan

A computational homogenization approach for the yield design of periodic thin plates Part I : Construction of the macroscopic strength criterion

Keywords: yield design, limit analysis, homogenization theory, thin plate model, second-order cone programming, finite element method

The purpose of this paper is to propose numerical methods to determine the macroscopic bending strength criterion of periodically heterogeneous thin plates in the framework of yield design (or limit analysis) theory. The macroscopic strength criterion of the heterogeneous plate is obtained by solving an auxiliary yield design problem formulated on the unit cell, that is the elementary domain reproducing the plate strength properties by periodicity. In the present work, it is assumed that the plate thickness is small compared to the unit cell characteristic length, so that the unit cell can still be considered as a thin plate itself. Yield design static and kinematic approaches for solving the auxiliary problem are, therefore, formulated with a Love-Kirchhoff plate model. Finite elements consistent with this model are proposed to solve both approaches and it is shown that the corresponding optimization problems belong to the class of second-order cone programming (SOCP), for which very efficient solvers are available. Macroscopic strength criteria are computed for different type of heterogeneous plates (reinforced, perforated plates,...) by comparing the results of the static and the kinematic approaches. Information on the unit cell failure modes can also be obtained by representing the optimal failure mechanisms. In a companion paper, the so-obtained homogenized strength criteria will be used to compute ultimate loads of global plate structures.

Introduction

Heterogeneous periodic thin plates are frequently encountered in civil engineering applications and the assessment of their bending strength capacities is of great importance for engineers. The computation of the ultimate load of a structure can be performed using two different class of methods. The first one, called the incremental approach, relies on a step-by-step elasto-plastic computation of the whole loading path until failure. This approach is time consuming, especially for complex structures, and poses convergence issues when approaching the collapse of the structure. The second class concerns direct methods using the theory of limit analysis or, in a more general manner, yield design theory [START_REF] Salençon | Yield Design[END_REF] to bracket the ultimate load using two theorems : the lower bound static approach and the upper bound kinematic approach. The efficiency of direct methods is that they require only the verification of equilibrium equations and the fulfillment of the yield criterion at each point of the structure without any knowledge of the mechanical behavior (apart from the strength criterion) or the whole loading path. The resolution of the static and kinematic approaches requires solving convex nonlinear optimization problems. Numerical methods dedicated to yield design have gained recent attention due to the development of mathematical programming techniques. A traditional approach involves the linearization of the yield criteria so that the corresponding optimization problems can be formulated within linear programming [START_REF] Faccioli | A finite element, linear programming methods for the limit analysis of thin plates[END_REF][START_REF] Munro | Yield line method by finite elements and linear programming[END_REF][START_REF] Sloan | Lower bound limit analysis using finite elements and linear programming[END_REF][START_REF] Sloan | Upper bound limit analysis using finite elements and linear programming[END_REF], for which powerful softwares based on interior point algorithms are available. These algorithms have also been developed for a broader class of optimization problems called second order cone programming (SOCP) [START_REF] Andersen | Computing limit loads by minimizing a sum of norms[END_REF] and implemented in commercial codes such as the Mosek software package [START_REF] Mosek | The Mosek optimization toolbox for MATLAB manual[END_REF]. Remarkably, a large number of traditional yield criteria can be expressed using conic constraints so that limit analysis problems can be formulated within SOCP [START_REF] Makrodimopoulos | Remarks on some properties of conic yield restrictions in limit analysis[END_REF]. Recent works applied this method to 2D plane strain problems [START_REF] Makrodimopoulos | Upper bound limit analysis using simplex strain elements and secondorder cone programming[END_REF][START_REF] Ciria | Mesh adaptive computation of upper and lower bounds in limit analysis[END_REF], frame structures or thin plates in bending [START_REF] Le | Upper and lower bound limit analysis of plates using FEM and second-order cone programming[END_REF][START_REF] Bleyer | On the performance of non-conforming finite elements for the upper bound limit analysis of plates[END_REF]. The obtained results seem very promising in terms of computational time saving (problems with a large number of optimization variables can be solved within seconds) as well as accuracy. For these reasons, these numerical methods will be used in this work.

Despite the efficiency of the previously mentioned numerical procedures, yield design of periodic media can be very difficult to perform due to the presence of rapidly and strongly varying material properties on a large scale structure. Numerical computations on the heterogeneous structure are, therefore, out of reach due to the high degree of local refinement needed to correctly capture the properties of the inhomogeneities. Homogenization theory in yield design has therefore been developed to characterize the strength domain of an equivalent homogeneous media with the idea that the corresponding homogenized yield design problem would, then, be much easier to solve. Founding works are due to Suquet [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF] and de Buhan [START_REF] De Buhan | A fundamental approach to the yield design of reinforced soil structures[END_REF]. The determination of the homogenized or macroscopic strength properties are quite similar to homogenization theory in elasticity. The macroscopic strength domain is, indeed, determined by solving an auxiliary yield design problem formulated on the unit cell of the periodic plate. Homogenization of elastic in-plane periodic plates has been widely studied by different authors. One main feature of this problem is that a plate model is valid in the limit of a small thickness h compared to the typical length L of the plate structure whereas the homogenization procedure is valid in the limit of a small in-plane typical length a of the unit cell compared to the same length L. Therefore, different homogenization procedures have to be considered depending on the relative values of a and h. For instance, if a and h are of the same order, the unit cell has to be modeled as a 3D medium [START_REF] Bourgeois | Modélisation numérique des panneaux structuraux légers[END_REF][START_REF] Sab | Yield design of thin periodic plates by a homogenization technique and an application to masonry walls[END_REF][START_REF] Dallot | Limit analysis of multi-layered plates. Part I: the homogenized Love-Kirchhoff model[END_REF]. On the contrary, if h is small compared to a, it is possible to replace the original 3D heterogeneous body by a 2D heterogeneous Love-Kirchhoff plate, which is homogenized in a second step (see [START_REF] Duvaut | Homogénéisation d'une plaque mince en flexion de structure périodique et symétrique[END_REF] for elastic plates). Our work will be focused on this specific case within the framework of yield design theory and associated numerical methods. This paper is organized as follows : in section 2, the homogenization theory in yield design will be briefly described within the framework of thin plates in bending and the equilibrium equations of the associated Love-Kirchhoff plate model will be recalled. Section 3 is devoted to the formulation of the auxiliary yield design problem, either by the static approach or by the kinematic approach, and the definition of the macroscopic strength criterion G hom is given. Numerical methods to solve the static approach are then presented in section 4 and the corresponding optimization problem is formulated as a SOCP problem. Section 5 deals with the case of the kinematic approach in the same manner. Finally, different numerical examples are studied in section 6, in order to assess the performance of both numerical procedures.

Yield design of periodic thin plates : a homogenization approach

The heterogeneous problem

We consider a heterogeneous thin plate occupying a domain Ω in the (x, y)-plane. Internal forces of the plate are the tensor of membrane forces N , the tensor of bending moments M and the vector of shear forces V . Generally speaking, the set of admissible internal forces with respect to the local strength of the plate at a point x ∈ Ω can be represented as a bounded convex set G(x) :

(N, M , V ) ∈ G(x)

In the special case of thin plates in bending, it is generally assumed that the plate is infinitely resistant to both membrane and shear forces such that the local strength criterion depends on the bending moment only :

M ∈ G(x)
Now, assuming that the plate loading depends upon several loading parameters Q, the domain Λ of potentially safe loads Q is defined as the set of loads such that there exists a statically admissible (S.A.) bending moment field M(x) (i.e. which equilibrates the loading Q), satisfying the strength criterion at each point of the plate (see [START_REF] Salençon | Yield Design[END_REF]) :

Λ = Q | ∃M (x) S.A. with Q, ∀x ∈ Ω M (x) ∈ G(x)
Making use of the virtual work principle, one can obtain a kinematic definition of Λ, dual to the previous static one. In the case of thin plates in bending, the hypothesis of infinite membrane and shear strength imposes that the plate kinematics obey the Love-Kirchhoff condition. Let u be the virtual transversal velocity of the plate and q the generalized kinematic parameters defined by duality in the expression of the work of external forces, such that for all u kinematically admissible (K.A.) with q (i.e. piecewise continuous and differentiable satisfying the kinematic boundary conditions), the virtual work of external load is given by P ext ( u) = Q • q. We then introduce π( χ; x), the support function of G(x) defined as π( χ; x) = sup

M ∈G(x)
M : χ and the associated maximum resisting work P rm ( u) as follows 1 :

P rm ( u) = Ω π( χ; x)dΩ
where χ = ∇ s ∇ u(x) is the curvature tensor associated with the virtual velocity field u. The following kinematic definition of Λ is then obtained :

Q ∈ Λ =⇒ ∀ u K.A. with q, P ext ( u) ≤ P rm ( u)

The homogenized problem

Now, the special case of plates which are periodic in their in-plane direction will be considered. Therefore, there exist two vectors a 1 , a 2 such that G(x) can be reproduced by periodicity along a 1 and a 2 :

G(x + n 1 a 1 + n 2 a 2 ) = G(x) ∀x ∈ Ω, ∀n 1 , n 2 ∈ Z
The two vectors a 1 and a 2 define the unit cell of the periodic plate. In the case when the typical size a of the unit cell is small in comparison to the plate typical length L (a ≪ L), the natural idea of homogenization theory is to substitute the local heterogeneous strength criterion G(x) by a homogenized or macroscopic strength criterion G hom , as illustrated in figure 1.

Using the same definitions as before, we introduce

Λ hom = Q | ∃M (x) S.A. with Q, ∀x ∈ Ω M (x) ∈ G hom
which also admits the following kinematic definition : where Π hom ( χ) is the support function associated with G hom . Similarly to the case of 3D periodic heterogeneous media [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF][START_REF] De Buhan | A fundamental approach to the yield design of reinforced soil structures[END_REF], the following result of homogenization theory in yield design holds true :

Q ∈ Λ hom =⇒ ∀ u K.A. with q P ext ( u) ≤ P hom rm ( u) = Ω Π hom ( χ)dΩ
Λ -→ Λ hom when a/L → 0

Let mention that Γ-convergence results for the homogenization of rigid perfectly plastic periodic thin plates in the case a ≪ L are given in [START_REF] Lewiński | Plates, laminates, and shells: asymptotic analysis and homogenization[END_REF].

Love-Kirchhoff plate model

Since the plate is supposed to be infinitely resistant to shear forces, the virtual motions of the plate must satisfy the Love-Kirchhoff condition.

In this case, the weak form of the Love-Kirchhoff plate equations is obtained by writing the virtual work of internal forces for a continuously differentiable test function u(x) (transversal velocity) and integrating twice by parts :

P int (u) = - Ω M : ∇ s ∇udΩ = Ω div M • ∇udΩ - ∂Ω (M • n) • ∇udl + Γ [[M • n]] • ∇udl P int (u) = - Ω div (div M )udΩ + ∂Ω (div M • n)udl - Γ [[div M • n]]udl - ∂Ω (M • n) • ∇udl + Γ [[M • n]] • ∇udl
where Γ is a line of potential discontinuities in Ω, n a unit normal vector and where [[•]] represents the discontinuity of the quantity (•) through Γ following its normal n. In the case when the loading consists only of a distributed transversal pressure p(x), the virtual work of external forces is given by P ext (u) = -Ω p(x)u(x)dΩ. Then, the virtual work principle reads : 0 = P int (u) + P ext (u) ∀u K.A.

The weak formulation yields the following equilibrium equation : div (div M ) + p = 0 in Ω While the boundary conditions and jump equations require a specific attention because of the fact that, on a given line, u and ∇u•n are independent fields whereas u and ∇u•t are not (t being a unit tangent vector). For example, on the part of the boundary where u = ∇u•n = 0, it is impossible to choose an arbitrary value of ∇u•t to obtain the boundary condition M nt = 0. This specific aspect has been widely reported in the literature and one method to overcome this problematic aspect is to use the following equality

M nt ∂ t u = ∂ t (M nt u) -(∂ t M nt )u
in the variational form, here ∂ t denotes the derivative along the tangential direction. Introducing

K = (div M ) • n + ∂M nt ∂t
as the equivalent Kirchhoff shear force, one obtains :

P int (u) = - Ω div (div M)udΩ + ∂Ω Kudl - Γ [[K]]udl - ∂Ω M nn ∂ n udl + Γ [[M nn ]]∂ n udl - ∂Ω ∂ t (M nt u)dl + Γ [[∂ t (M nt u)]]dl
If ∂Ω ∪ Γ is regular except on a finite number of angular points c j , the last two terms can be rewritten as :

- ∂Ω ∂ t (M nt u)dl + Γ [[∂ t (M nt u)]]dl = - j R j u(c j )
where

R j = [[M nt (c j )]
] is the discontinuity of M nt at a corner c j . Finally, the variational formulation reduces to : ∀uK.A.

0 = Ω (div (div M) + p)udΩ + ∂Ω (M nn ∂ n u -Ku)dl - Γ ([[M nn ]]∂ n u -[[K]]u)dl + j R j u(c j )
The variational formulation is now written in terms of independent fields only and one can obtain the following boundary conditions :

K = 0 M nn = 0 on ∂Ω 1 K = 0 ∂ n u = 0 on ∂Ω 2 u = 0 M nn = 0 on ∂Ω 3 u = 0 ∂ n u = 0 on ∂Ω 4 ∂Ω = i ∂Ω i , ∂Ω i ∩ ∂Ω j = ∅ for i = j
and the following jump equations (figure 2):

[[M nn ]] = 0 [[K]] = 0 through Γ R j = [[M nt (c j )]] = 0 at c j

Determination of G hom by solving an auxiliary yield design problem

The process of determining the macroscopic strength criterion G hom requires to solve a specific auxiliary yield design problem formulated on the unit cell, allowing to relate macroscopic quantities to the microscopic ones. Obviously, it is to be recalled that the thin plate model is only an idealization of a 3D structure in the limit of a small thickness h compared to the characteristic length of the structure L. As already mentioned in the introduction, the present paper is interested in the case when the thickness h can still be considered as small when compared to the typical size of the unit cell : h ≪ a ≪ L. In this specific case, the original 3D unit cell problem can still be modeled as a plate problem.

Space of statically admissible moment fields and definition of G hom

Let A be such a unit cell of the periodic plate, ξ a point in A and Γ a potential discontinuity line in A. M is the macroscopic bending tensor which will play the role of loading parameters for the auxiliary problem. The appropriate conditions to formulate a correct auxiliary problem in the framework of periodic homogenization theory are the following :

• zero distributed forces;

• periodic conditions on the boundary of the unit cell (see figure 3);

• averaging relation relating the macroscopic to the microscopic bending moment tensors.

Hence, we introduce the set SA(M ) of local bending moment tensor m(ξ) defined on the unit cell which are statically admissible with M as follows :

div div m(ξ) = 0 in A (1) [[m nn ]] = 0 through Γ (2) [[K]] = 0 through Γ (3) m(ξ) ∈ SA(M ) ⇐⇒ R j = 0 at c j inside A (4) m nn A -periodic (5) K A -antiperiodic (6) R j A -antiperiodic (7) M = m(ξ) (8)
where

• = 1 |A| A (•)dA.
The macroscopic strength criterion G hom is then obtained as follows :

G hom = M | ∃m(ξ) ∈ SA(M ), ∀ξ ∈ A m(ξ) ∈ G(ξ) (9)
which states that the macroscopic strength domain consists of all the macroscopic bending moments for which there exists at least one microscopic bending moment distribution on the unit cell which is statically admissible with the macroscopic bending moment, while fulfilling the local strength criterion at each point of the unit cell.

Space of kinematically admissible velocity fields

Let KA(χ) bet the set of local transversal velocity fields u(ξ) defined on the unit cell, and kinematically admissible with a macroscopic curvature tensor χ as follows :

u(ξ) = 1 2 ξ • χ • ξ + v(ξ) in A (10) u(ξ) ∈ KA(χ) ⇐⇒ v(ξ) continuous in A (11) v(ξ) A -periodic (12) ∂v ∂n (ξ) A -antiperiodic (13)
This definition deserves some remarks :

• Due to the previous assumptions (Love-Kirchhoff kinematics), the following properties can be deduced :

u(ξ) is continuous;

the discontinuity of the gradient ∇u through a discontinuity line Γ is necessarily normal to this line by virtue of Hadamard's compatibility condition [START_REF] Salençon | Handbook of Continuum Mechanics: General Concepts, Thermoelasticity[END_REF]. Hence

: [[∇u]] = ∂v ∂n n -∇v(ξ) is A-periodic.
• Assumptions made on v(ξ) are consistent with Love-Kirchhoff kinematics. In particular, the local velocity field u can be seen as a Love-Kirchhoff velocity field obtained by the superposition of a homogeneous curvature term plus a periodic perturbation v(ξ).

• Rotation of the microstructure is characterized by the gradient of u : ∇u(ξ) = χ • ξ + ∇v(ξ) which is also obtained by the superposition of a homogeneous curvature term plus a periodic perturbation ∇v(ξ) of the rotation vector.

• The local curvature is κ(ξ) = ∇ s ∇u = χ + ∇ s ∇v(ξ). Using Green's theorem, its average is given by :

κ = χ + ∇ s ∇v = χ + 1 |A| ∂A 1 2 (∇v ⊗ n + n ⊗ ∇v)dl - Γ 1 2 ([[∇v]] ⊗ n + n ⊗ [[∇v]])dl
the first term vanishes because ∇v ⊗ n is anti-periodic. Thus, the following identity holds true :

χ = κ + 1 |A| Γ ∂v ∂n n ⊗ ndl (14)
Using the virtual work principle, it is quite easy to see that the two sets SA(M) and KA(χ) are in duality in the following sense :

∀m ∈ SA(M ), ∀u ∈ KA(χ) P ext (u) = |A| m : κ(ξ) = |A|M : χ i.e
. the average internal work on the unit cell is equal to the internal work of the average bending moment M and of the macroscopic curvature tensor χ. More precisely, owing to [START_REF] De Buhan | A fundamental approach to the yield design of reinforced soil structures[END_REF] in the case when there are no rotation discontinuities of the velocity field on the unit cell, χ is equal to the average of the local curvature κ(ξ) such that the previous relation can now be written :

m : κ(ξ) = m(ξ) : κ(ξ) = M : χ
This expression can be interpreted as a generalized Hill's lemma which is the basis of all homogenization procedures. In our case, M clearly appears as a three-parameter (M xx , M yy , M xy ) loading mode, the associated dual kinematic parameters being the three components of the macroscopic curvature tensor (χ xx , χ yy , χ xy ).

Kinematic definition of G hom

Using the previous notations, it can be shown that the macroscopic support function Π hom is given by :

Π hom (χ) = sup M ∈G hom M : χ = inf u∈KA(χ) 1 |A| A π(κ; ξ)dA + Γ π([[∇u]]; n)dl or introducing KA 0 such that KA(χ) = 1 2 ξ • χ • ξ ⊕ KA 0 : Π hom (χ) = inf v∈KA 0 1 |A| A π(χ + ∇ s ∇v; ξ)dA + Γ π ∂v ∂n ; n dl (15)
Using this kinematic definition, G hom can then be computed as the convex hull of all hyperplanes of equation M : χ = Π hom (χ) for all macroscopic curvatures χ. Note that, since G has been supposed to be bounded, G hom is also bounded (see definition [START_REF] Makrodimopoulos | Upper bound limit analysis using simplex strain elements and secondorder cone programming[END_REF]). Therefore, all suprema and infinima in the previous relations are finite and there is no additional condition on the virtual curvatures to ensure finite upper bounds in the kinematic approach.

Conclusion

In this section, an appropriate definition for the macroscopic strength criterion G hom of a periodic plate has been given. This definition relies upon the solution of an auxiliary yield design problem defined over the unit cell, with specific periodic boundary conditions. This problem can be solved either in its static or its kinematic form. In practice, none of the two problems can be solved analytically and numerical procedures are required based for instance upon the discretization of both problems. Let SA d (M ) ⊆ SA(M) and KA 0d ⊆ KA 0 be subspaces consistent with the previous definitions. If ( 9) is solved using SA d (M ), one would obtain a lower bound approximation G hom stat to the macroscopic strength criterion G hom . Similarly, if [START_REF] Bourgeois | Modélisation numérique des panneaux structuraux légers[END_REF] is solved using KA 0d , one would obtain a larger value of the macroscopic support function Π kin hom ≥ Π hom , corresponding to an upper bound approximation G hom kin ⊃ G hom to the macroscopic strength criterion. Obvisouly, using the lower (resp. upper) approximation of G hom in the global computation of a homogenized structure will lead to a lower bound (resp. upper bound) of the exact homogenized limit load. Such numerical procedures based on the finite element method are discussed in the next two sections and it will be shown, on numerical examples, that the obtained approximations of G hom are very close estimates of the exact macroscopic strength criterion.

Numerical resolution of the static approach

In this section, the static problem ( 9) is discretized using the finite element method and the resulting optimization problem is solved using second-order cone programming. Lower bound limit analysis of thin plate problems have already been successfully solved using such methods [START_REF] Le | Upper and lower bound limit analysis of plates using FEM and second-order cone programming[END_REF]. Since the present auxiliary problem is very similar to yield design problems formulated on plate structures, the formulations will be simply adapted to take periodic boundary conditions and averaging relations into account.

Finite element discretization

In the following, the unit cell A is discretized into N E triangular finite elements assuming a constant varying moment field in each element [START_REF] Le | Upper and lower bound limit analysis of plates using FEM and second-order cone programming[END_REF][START_REF] Olsen | The influence of the linearisation of the yield surface on the load-bearing capacity of reinforced concrete slabs[END_REF]. Each such element, known as the equilibrium Morley element, exhibits three degrees of freedom {m e } = {m xx m yy m xy } T per element. Despite its apparent simplicity, some difficulties may arise for yield design problems involving for instance uniformly distributed pressures, since the Morley element does not result in an exact equilibrium relation for such loadings. In that case, some modifications are required, as mentioned in [START_REF] Le | Upper and lower bound limit analysis of plates using FEM and second-order cone programming[END_REF]. However, since in the specific case of periodic homogenization all distributed loads are supposed to be zero, such modifications are unnecessary.

With constant moments, equations (1), ( 3) and ( 6) are identically satisfied. Hence, the generalized forces consist of three internal corner forces Z 1 , Z 2 , Z 3 and three constant normal bending moments along edges m 12 , m 23 , m 31 (figure 4). These generalized forces (collected in a vector {G}) are related to the moment field in the element {m e } by a constant matrix 2) and ( 4) exactly. It should be mentioned that, in the global equilibrium matrix [C], the equilibrium of edges and nodes belonging to the unit cell boundary are not taken into account. The latter require, indeed, a special treatment to impose the periodic boundary conditions ( 5) and [START_REF] Mosek | The Mosek optimization toolbox for MATLAB manual[END_REF].

Let an edge belonging to the boundary of the unit cell be indexed by i, and let e be the index of its corresponding element. Let i ′ (resp. e ′ ) be the index of the edge (resp. of its corresponding element) on the opposite boundary2 . The two normal moments on i and i ′ can then be obtained as :

m (i) nn = c e {m e } m (i ′ ) nn = c e ′ {m e ′ }
The periodicity of m nn is then enforced by the following condition :

m (i) nn -m (i ′ ) nn = c e {m e } -c e ′ {m e ′ } = 0
The same kind of equation can be written concerning corner loads for all pair of opposite nodes. Both equations are then assembled into a global matrix [P ] taking into account all periodic boundary conditions :

[P ]{m} = {0}
Finally, if a e is the area of element number e and [I 3 ] the 3 × 3 identity matrix, the averaging condition (8) may be written as :

{M} =    M xx M yy M xy    = 1 |A| a 1 [I 3 ] a 2 [I 3 ] . . . a N E [I 3 ] {m} = [A]{m}
From now on, the unit cell will be considered of area |A| = 1.

Finally, it should be mentioned that all equations defining SA(M ) have been exactly satisfied by the proposed discretization. Therefore, if SA d (M ) is the set of vectors {m} satisfying all the previous discretized equations, the following inclusion SA d (M ) ⊂ SA(M ) holds true.

Formulation as a standard second-order cone programming problem

The bending moment field has to satisfy the local strength criterion G(ξ) at each point ξ ∈ A. In the following, it will be assumed that the unit cell is made of materials obeying the von Mises criterion :

G(ξ) = m(ξ) | m 2 xx + m 2 yy -m xx m yy + 3m 2 xy -m p (ξ) ≤ 0 ( 16
)
where m p (ξ) is the ultimate bending moment3 characterizing the strength of the plate at point ξ.

Since m(ξ) is constant over each element e, the yield criterion will be satisfied if :

{m e } T [P vM ]{m e } -m p,e ≤ 0 with [P vM ] =   1 -1/2 0 -1/2 1 0 0 0 3  
where m p,e = inf m p (ξ) for ξ in the element e. In the numerical examples, m p (ξ) will be taken equal to a constant m p,e in each element. Since [P vM ] is positive definite, the yield criterion can be formulated as the following conic constraint :

[J 1 ]{m e } ≤ m p,e where [J 1 ] =   1 -1/2 0 0 √ 3/2 0 0 0 √ 3   is the Cholesky factor of [P vM ].
Let us mention that, although only the von Mises criterion will be used in the following, our study can be easily extended to other types of bending plate criteria, provided they can be formulated using conic constraints [START_REF] Makrodimopoulos | Remarks on some properties of conic yield restrictions in limit analysis[END_REF][START_REF] Bleyer | On the performance of non-conforming finite elements for the upper bound limit analysis of plates[END_REF]. This is the case for the Johansen square criterion, the Tresca criterion, etc... Since G hom ⊂ R 3 , it will be numerically obtained by considering a specific unitary direction {M 0 } ∈ R 3 , {M 0 } = 1 and a loading parameter λ such that {M} = λ{M 0 }. Thus, the static approach reduces to the resolution of the following maximization problem :

λ s = max λ s.t.   0 [C] 0 [P ] {M 0 } -[A]   λ {m} =    0 0 0    [J 1 ]{m e } ≤ m p,e ∀e = 1, . . . , N E (17)
This maximization problem is a second-order cone programming (SOCP) problem involving linear equality constraints and N E conic constraints.

x y Finally, since SA d (M ) ⊂ SA(M) and since the strength criterion is exactly satisfied over the whole unit cell, it can be concluded that λ s {M 0 } ∈ G hom . Therefore, the convex hull of all λ s {M 0 } for different values of the direction {M 0 } defines a lower bound approximation G hom stat ⊂ G hom .

u i ∂ x u i ∂ y u i u g 1,2,

Numerical resolution of the kinematic approach

The development of efficient finite elements for the upper bound kinematic approach of thin plates in bending was the main contribution of [START_REF] Bleyer | On the performance of non-conforming finite elements for the upper bound limit analysis of plates[END_REF]. More precisely, the present authors proposed to use non-conforming finite elements which incorporate not only curvature-type deformations, but also rotation discontinuities across element edges. In particular, the cubic Hermite triangle (H3) seems to be a very efficient element in terms of numerical accuracy as well as computational time saving. This element will, therefore, be used to solve [START_REF] Bourgeois | Modélisation numérique des panneaux structuraux légers[END_REF]. Some slight modifications will be made to incorporate the specific features of (15); for more details on the other aspects of the formulation, we refer to [START_REF] Bleyer | On the performance of non-conforming finite elements for the upper bound limit analysis of plates[END_REF].

Finite element discretization

The cubic Hermite triangle (H3) is a non-conforming finite element assuming a cubic variation of the transversal displacement, ensuring C 0 -continuity. It has 10 degrees of freedom : the transversal displacement and its two derivative along x and y at corner nodes plus the transverse displacement at the centroid (figure 5).

In fact, in our case, the unknown will not be the total velocity u(ξ) but the periodic fluctuation v(ξ), as the implementation seemed simpler this way. Let {V } be the vector collecting all degrees of freedom. The curvature associated with v(ξ) at a point ξ in e is computed using a matrix [B(ξ)] involving second derivatives of shape functions :

   ∂ xx v ∂ yy v 2∂ xy v    = [B(ξ)]{V }
The total curvature is then obtained by adding the constant term {χ} due to the macroscopic curvature :

{κ} =    κ xx κ yy 2κ xy    = {χ} + [B(ξ)]{V }
In the case of a von Mises criterion, the local support function writes :

π(κ; ξ) = 2m p √ 3 (ξ) κ 2 xx + κ 2 yy + κ xx κ yy + κ 2 xy
As before, this expression can be written as :

π(κ; ξ) = m p (ξ) [J 2 ]{κ} where [J 2 ] =   2/ √ 3 1/ √ 3 0 0 1 0 0 0 1/ √ 3   .
The contribution of the curvature term in the expression of the maximum resisting work (P curv rm ) is then evaluated on the whole unit cell by a n = 3 points Gauss quadrature :

P curv rm = N E e=1 m p,e n g=1 c e,g [J 2 ]{χ} + [J 2 ][B e,g ]{V }
where c e,g are constant terms coming from the quadrature of the integral.

Since C 1 -continuity is not assured, we have to take into account the contribution of rotation discontinuities across element edges to the maximum resisting work.

The support function associated to such discontinuities of normal rotation ∆θ n is :

π(∆θ n ; n) = π(κ = ∆θ n n ⊗ n; ξ) = 2m p (ξ) √ 3 |∆θ n |
where the rotation discontinuity across edge j can be written :

∆θ n = ∆Θ j {V }
The contribution of rotation discontinuities in the expression of the maximum resisting work (P disc rm ) is also evaluated over all N D active 4 edges by a m = n = 3 points Gauss quadrature on the element edges :

P disc rm = N D j=1 m p,j m g ′ =1 c ′ j,g ′ | ∆Θ j {V }|
where c ′ j,g ′ are constant terms coming from the quadrature of the integral and where m p,j is chosen as the minimum value of m p (ξ) over each side of edge j.

Finally, periodic boundary conditions ( 12) and ( 13) have to be satisfied. First, it is to be noticed that, since v admits a cubic interpolation over an element, it is not sufficient to impose that v is periodic at both end nodes of an edge for v to be periodic over the whole edge. Two supplementary conditions on the derivative of v along the edge have to be imposed. Hence, we must impose that ∂ t v is periodic too. Therefore, since v and ∂ t v must be periodic and ∂ n v must be antiperiodic, it is equivalent to impose that v, ∂ x v and ∂ y v are periodic. For all pair (i, i ′ ) of opposite nodes, this can be written as :

   v i ∂ x v i ∂ y v i    -    v i ′ ∂ x v i ′ ∂ y v i ′    = {0}
which leads to the condensed form [P ′ ]{V } = {0}.

Formulation as standard second-order cone programming

For a given value of χ, the approximate macroscopic support function Π kin hom (χ) is then obtained by solving the following optimization problem :

Π kin hom (χ) = min N E e=1 m p,e n g=1 c e,g [J 2 ]{χ} + [J 2 ][B e,g ]{V } + N D j=1 m p,j m g ′ =1 c ′ j,g ′ | ∆Θ j {V }| s.t. [P ′ ]{V } = {0}
Introducing various auxiliary variables, the preceding problem can be reformulated as :

Π kin hom (χ) = min N E •n k=1 m p,e c k t k + N D •m k ′ =1 m p,j c ′ k ′ s k ′ s.t. [P ′ ]{V } = {0} [J 2 ] -1 {r k } -[B e,g ]{V } = {χ} t k ≥ {r k } u k ′ = ∆Θ j {V } s k ′ ≥ |u k ′ | (18)
4 i.e. which actually contribute to P disc rm which amounts to a standard SOCP problem involving linear equalities and inequalities, as well as N E • n conic constraints.

Finally, since the curvature is exactly computed, and since conditions ( 10)-( 13) are exactly satisfied over the whole unit cell, the following inclusion KA 0h ⊂ KA 0 holds true where KA 0h is the space of periodic fluctuations obtained with this discretization. Therefore, if numerical errors due to the quadrature of the integral can be neglected, it can be concluded that Π hom (χ) ≤ Π kin hom (χ) with Π kin hom being the support function of an upper bound approximation G hom kin ⊃ G hom of the macroscopic strength criterion.

Numerical examples

In the following, different examples of periodic plates will be treated such as reinforced or perforated plates. Lower and upper approximations will be computed by solving the static and kinematic approach on the unit cell using the methods presented in sections 4 and 5. Unit cells are meshed using Gmsh, optimization problems ( 17) and ( 18) are formulated under the Matlab environment and solved using the dedicated SOCP solver Mosek [START_REF] Mosek | The Mosek optimization toolbox for MATLAB manual[END_REF]. For a given value of the macroscopic variables (M or χ), the resolution of the associated optimization problem derived from the static or the kinematic approach, took in general less than 1-2 seconds. Pre-processing including mesh generation and matrices assembling procedures took approximately 10-20 seconds.

Plate reinforced in one direction

The unit cell geometry represented in figure 6 corresponding to a plate made of a reference material of ultimate moment m p1 = 1 reinforced by a band of width e = 0.05 along the y direction made of a material of ultimate moment m p2 = 40 will first be considered.

A trivial static approach consisting of a homogeneous bending moment field m(ξ) = M shows that G hom contains the domain defined by the von Mises criterion of the weakest material m p1 = 1, denoted by G vM (m p1 ). Similarly, a trivial kinematic approach with v(ξ) = 0 (homogeneous curvature) shows that G hom is contained in the domain defined by the von Mises criterion of average ultimate moment m p = 0.95

• 1 + 0.05 • 40 = 2.95, denoted by G vM ( m p ).
First, the intersection of G hom with the plane M xy = 0, noted G hom (xx,yy) will be investigated. Due to the unit cell symmetry with respect to its center, it is easily shown that this intersection is also equal to the projection onto the plane M xy = 0, which is also obtained by computing Π hom (χ xx , χ yy , χ xy = 0). Therefore, G hom (xx,yy) is a convex domain of R 2 , the support function of which is given by Π hom (χ xx , χ yy , 0). The results of both approaches are represented in figure 7. It can be observed that the points obtained from the static approach correspond exactly to the convex hull of the tangent hyperplanes obtained from the kinematic approach, ensuring that G hom (xx,yy) is exactly computed. As expected, the reinforcement is ineffective concerning bending moments along M xx , only whereas for bending moments along M yy only, the average strength m p = 2.95 is recovered. Moreover, it is possible to show that G hom (xx,yy) is actually given by the following analytical relations :

(M xx , M yy ) ∈ G hom (xx,yy) ⇐⇒ ∃(X, Y ) such that    (M xx , X) ∈ G vM (m p1 ) (M xx , Y ) ∈ G vM ( m p ) M yy = (1 -η)X + ηY
where η is the fraction of the reinforcement, i.e. η = e/a. These relations mean that a point on ∂G hom (xx,yy) is obtained as the barycenter of two points on ∂G vM (m p1 ) and ∂G vM ( m p ) which share the same abscissa M xx , the weighting factor being given by 1η and η. This construction is very similar to the construction of the macroscopic strength criterion of a multilayer purely cohesive material under plane strain conditions [START_REF] De Buhan | A fundamental approach to the yield design of reinforced soil structures[END_REF]. In the limit of a small reinforcement fraction, G hom (xx,yy) can be obtained by translating the von Mises ellipse with m p = 1 (in red) by a vector (M xx = 0, M yy = ±1.95) which corresponds to the effect of the reinforcement along the direction y.

In figure 8, the whole domain G hom has been represented in the (M xx , M yy , M xy ) space. Black dots correspond to points obtained from the static approach, whereas the gray domain (xx,yy) for the reinforced plate obtained by the static approach (black dots). Blue tangent planes are obtained with Π hom (χ xx , χ yy 0) from the kinematic approach correspond to the convex hull of the obtained from the kinematic approach. Clearly, both approaches lead to the same convex set with an excellent precision. It is also worth noting that G hom has a relatively complex shape (which is not an ellipsoid) bounded by two planes of equation M xx = ± 2m p1 √ 3 .

Plate perforated by a rectangular hole

The unit cell geometry represented in figure 9, corresponding to a plate made of a reference material of ultimate moment m p1 = 1 perforated by a rectangular hole of width e = 0.1 and length l = 0.6, will now be considered. It is to be mentioned that this example involves a moment-free boundary (noted B) inside the unit cell, which was not considered in the description of the auxiliary problem in section 3. The only modifications for the static approach are that the relations m nn = 0 and K = 0 have to be imposed on B. Concerning the kinematic approach, the only modification is that edges belonging to B remain inactive, i.e. they do not contribute to P disc rm via rotation discontinuities.

In figure 10, isovalues of the different components of the local bending moments m(ξ) obtained by the static approach for a macroscopic bending moment of the form M = e x ⊗ e x Figure 8: Domain G hom for the reinforced plate obtained by the static approach (black dots). The grey domain is obtained as the convex hull of the hyperplanes obtained from the kinematic approach.

Figure 9: Geometry and mesh of the unit cell for the plate perforated by a rectangular hole have been represented. It can be observed that the components satisfy all the conditions defining SA(M ) (periodicity, boundary conditions on the free boundary, averaging relation,...). The complexity of the spatial variations of the local bending moment also underlines the need of numerical methods to solve the auxiliary problem, despite the simplicity of the macroscopic moment.

In figure 11, the failure mechanism (i.e. the field u(ξ) ∈ KA(χ)) obtained by the kinematic approach has been represented for different values of the macroscopic curvature : pure curvature in the x-direction (χ xx = 1, χ yy = χ xy = 0), pure curvature in the y-direction (χ xx = 0, χ yy = 1, χ xy = 0) and pure xy-torsion (χ xx = χ yy = 0, χ xy = 1). The mechanism obtained for χ xx = 1 in figure 10(a) corresponds to a yield line in the direction of the hole length. On the contrary for χ yy = 1 (figure 10(b)), a better mechanism than a simple yield line along the hole width is obtained via a complex local field with double curvature regions around the hole. Finally, for χ xy = 1 (figure 10(c)), the optimal mechanism seems to be a homogeneous torsion curvature field u(ξ) = 1 2 ξ • χ • ξ = ξ x ξ y obtained with a zero periodic fluctuation v(ξ) = 0. Those results also illustrate the complexity of the auxiliary problem, justifying the use of numerical methods to solve the static as well as the kinematic approach.

Finally, in figure 12 the intersection of G hom with the plane M xy = 0 has been represented for various values of the hole length l while keeping the hole width equal to e = 0.1. In this example, it seemed interesting to study the influence of pattern arrangement in the case of circular holes. Holes can, for example, be disposed in staggered rows or in straight rows. The ultimate bending strength can be estimated by solving the problem of an infinite plate solicited by a bending moment at infinity, or equivalently, by solving the auxiliary yield design problem on the unit cell with this bending moment as the macroscopic loading. Hence, in the following, the auxiliary problem will be solved to assess the ultimate strength in bending of those plates for different hole patterns.

Different unit cell geometries are considered in function of the pattern angle θ. The pattern geometries are obtained by translating circular holes of radii R along vectors a 1 = ae x and a 2 = a sin θe x + a cos θe y (see figure 13). The unit cell area is then |A| = a 2 cos θ. The opening area ratio (or porosity) is defined as η = πR 2 /|A|. Since we want to investigate the influence of the pattern angle on the macroscopic strength properties, a constant opening area ratio will be assumed for different values of θ. Hence, the radius R is set to R = R 0 √ cos θ where R 0 is the hole radius for θ = 0 such that η = πR 2 0 /a 2 is now independent of θ.

Due to the chosen geometry, the strength capacity which will be most influenced by the pattern angle is the maximum bending moment in the x-direction, denoted by M xx,0 and defined as : This quantity has been computed using the static approach for varying values of θ. In figure 14, the variation of M xx,0 as a function of pattern angle θ is represented. Three different values of the opening area ratio have been considered : η = 0.03, η = 0.13 and η = 0.28. Obviously, for a non-perforated plate (η = 0), we have M xx,0 = m p = 1. Results indicate that M xx,0 is decreasing with increasing values of η but the evolution of M xx,0 with respect to θ is non-monotonous. Besides, the influence of the pattern angle is increasing with η. In particular, for η = 0.28, the staggered 30 • pattern has a reinforcing effect of approximately 30% compared to the straight pattern (θ = 0 • ). This reinforcement effect is around 15% for η = 0.13 but only 2% for η = 0.03. Therefore, the pattern arrangement may have a non-negligible reinforcing effect on the homogenized strength capacities of the plate for high opening area ratios. Now, defining e α = cos αe x + sin αe y , the ultimate bending moment in the α-direction M αα,0 is defined as :

M αα,0 = max M =M e α ⊗e α ∈G hom

M

It is interesting to note that in the case of a straight pattern (θ = 0 • ), M αα,0 tends to increase with α from 0 • to 45 • and to decrease symmetrically (due to the orthogonal symmetry of the problem) from 45 • to 90 • (figure 15(a)). This effect can be explained because the smallest path between two holes in the orthogonal direction of the loading is equal to a for α = 0 • or 90 • and equal to √ 2a for α = 45 • . On the contrary, in the case of staggered pattern (θ = 30 • ), M αα,0 decreases with α from 0 • to 30 • and then decreases symmetrically (due to the hexagonal symmetry of the problem) from 30 • to 60 • (figure 15(b)). In this case, the smallest path is larger than a for α = 0 • or 60 • , whereas it is equal to a for α = 30 • . Hence, for the straight pattern the optimal loading direction seems to be around 45 • whereas for the staggered pattern the optimal loading direction is aligned with the holes.

Finally,the intersection of G hom with the plane M xy = 0 has been represented in figure 16 for η = 0.13 and two values of the pattern angle : θ = 0 • and θ = 30 • . The reinforcing effect of the 30 • pattern, compared to the straight pattern, can clearly be observed, not only for the particular value of M xx,0 , but also for the whole strength domain.

Conclusion

Yield design of thin periodic plates in bending has been addressed using a homogenization procedure to derive the macroscopic strength criterion of an equivalent homogeneous plate. The macroscopic strength criterion is obtained by solving an auxiliary yield design problem formulated on the unit cell characterizing the periodic microstructure. Static and kinematic approaches are formulated using a thin plate model, which is valid in the limit h ≪ a of a small plate thickness h compared to the unit cell characteristic length a. Numerical tools are introduced to solve both approaches efficiently in terms of computational time and accuracy. Indeed, the discrepancy between the evaluation of the macroscopic domain obtained from the two different approaches was negligible for all examples. The global efficiency and performance of the proposed method is based upon the following decisive ingredients :

• The discretization of the static approach is performed with finite elements satisfying exactly all equilibrium equations, periodic boundary conditions and local yield criterion.

• The discretization of the kinematic approach is performed with non-conforming finite elements, whose performance regarding upper bound limit analysis of thin plates had previously been demonstrated in [START_REF] Bleyer | On the performance of non-conforming finite elements for the upper bound limit analysis of plates[END_REF].

• The related optimization problems are formulated as standard SOCP problems for which a dedicated software package has already been successfully employed for solving yield design problems.

It has been shown on various examples that the macroscopic strength criteria are usually complicated anisotropic three-dimensional surfaces in the space of bending moments. Therefore, even if such criteria can be obtained by efficient numerical procedures as those presented in this paper, the main difficulty of homogenization theory in yield design relies on the implementation of such complex criteria in a yield design problem formulated on the global structure. Approximate representations of those yield criteria are, therefore, required to perform such computations. This specific aspect will be the purpose of the second part of this work.

Figure 1 :

 1 Figure 1: Illustration of the homogenization procedure for a periodic heterogeneous thin plate : initial heterogeneous (left) and equivalent homogeneous (right) yield design plate problems

Figure 2 :

 2 Figure 2: Jump equations of the Love-Kirchhoff plate model : through a line Γ (left) and at a corner c j (right)

Figure 3 :

 3 Figure 3: Periodic and antiperiodic boundary conditions on the unit cell for static and kinematic fields. ξ + and ξ -are the position vectors of two points on opposite sides of the unit cell and deduced from each other by periodicity.

Figure 4 :

 4 Figure 4: Equilibrium Morley finite element (from [11])

[

  C e ] depending on the geometry of the element : {G} = [C e ]{m e } If {m} collects the degrees of freedom of all N E elements, the global equilibrium of the unit cell is obtained by assembling all contributions into the following relation : [C]{m} = {0} thereby satisfying (

Figure 5 :

 5 Figure 5: H3 finite element (from [12])

Figure 6 :

 6 Figure 6: Geometry and mesh of the unit cell for the plate reinforced by periodic bands in the y-direction

Figure 7 :

 7 Figure 7: Domain G hom (xx,yy) for the reinforced plate obtained by the static approach (black dots). Blue tangent planes are obtained with Π hom (χ xx , χ yy 0) from the kinematic approach

Figure 10 : 1 Figure 11 :Figure 12 : 6 . 3 .

 101111263 Figure 10: Mapping of the local bending moment m(ξ) for a macroscopic bending moment of the form M = e x ⊗ e x

Figure 13 :

 13 Figure 13: Perforated plate geometry for a θ = 30 • pattern angle

Figure 14 :Figure 15 :Figure 16 :

 141516 Figure 14: Variation of M xx,0 as a function of pattern angle θ for different values of the opening area ratio

This expression assumes that the rotation vector associated with the gradient of the transversal velocity field u is everywhere continuous. If this is not the case, another term taking into account the contribution of angular jumps has to be considered in the expression of the maximum resisting work. For more details, we refer to[START_REF] Bleyer | On the performance of non-conforming finite elements for the upper bound limit analysis of plates[END_REF] and section 5.

Obviously, two opposite boundaries must have exactly the same mesh so that periodic boundary conditions can be imposed.

m p (ξ) = σ 0 (ξ)h 2 /4 if σ 0 is the uniaxial yield strength of the plate material


Jeremy Bleyer
email: jeremy.bleyer@enpc.fr

Patrick De Buhan

A computational homogenization approach for the yield design of periodic thin plates Part II : Upper bound yield design calculation of the homogenized structure

Keywords: yield design, limit analysis, periodic thin plates, yield surface approximation, second-order cone programming, finite element method

In the first part of this work [Bleyer, J., de Buhan, P., 2013, A computational homogenization approach for the yield design of periodic thin plates -Part I : Construction of the macroscopic strength criterion], the determination of the macroscopic strength criterion of periodic thin plates has been addressed by means of the yield design homogenization theory and its associated numerical procedures. The present paper aims at using such numerically computed homogenized strength criteria in order to evaluate limit load estimates of global plate structures. The yield line method being a common kinematic approach for the yield design of plates, which enables to obtain upper bound estimates quite efficiently, it is first shown that its extension to the case of complex strength criteria as those calculated from the homogenization method, necessitates the computation of a function depending on one single parameter. A simple analytical example on a reinforced rectangular plate illustrates the simplicity of the method. The case of numerical yield line method being also rapidly mentioned, a more refined finite element-based upper bound approach is also proposed, taking dissipation through curvature as well as angular jumps into account. In this case, an approximation procedure is proposed to treat the curvature term, based upon an algorithm approximating the original macroscopic strength criterion by a convex hull of ellipsoids. Numerical examples are presented to assess the efficiency of the different methods.

Introduction

In this joint work, the yield design of periodic thin plate structures is investigated. The first part of this work has been dedicated to the determination of the homogenized strength properties of different periodic plates through the numerical computation of a macroscopic strength criterion by means of finite elements and mathematical programming. Homogenization theory in the framework of yield design (or limit analysis) of periodic structures has first been proposed in the work of Suquet [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF] and de Buhan [START_REF] De Buhan | A fundamental approach to the yield design of reinforced soil structures[END_REF], where a proper definition of the macroscopic strength criterion involving the resolution of an auxiliary yield design problem formulated on the unit cell has been given. An analytical determination of the macroscopic strength criterion is very rare (e.g. the case of the multilayered soil under plane strain [START_REF] De Buhan | A fundamental approach to the yield design of reinforced soil structures[END_REF]) and often restricted to symmetric unit cell geometries and simple macroscopic loading [START_REF] Maghous | Détermination du critère de résistance macroscopique d'un matériau hétérogène à structure périodique. Approche numérique[END_REF]. Therefore, numerical methods are required, notably to conveniently capture the anisotropy of the homogenized material when preferential reinforcing directions are involved. The numerical resolution of the auxiliary problem can be tackled using incremental elasto-plastic approaches [START_REF] Marigo | Plasticité et homogénéisation: un exemple de prévision des charges limites d'une structure hétérogène périodique[END_REF] but a more natural method is to perform numerical limit analysis computations directly. This method, in conjunction with a finite element discretization, has been widely applied to different type of structures like porous media [START_REF] Pastor | Approche numérique des charges limites pour un matériau orthotrope de révolution en déformation plane[END_REF][START_REF] Turgeman | Comparaison des charges limites d'une structure hétérogène et homogénéisée[END_REF], periodic plates solicited in their own plane [START_REF] Maghous | Détermination du critère de résistance macroscopique d'un matériau hétérogène à structure périodique. Approche numérique[END_REF][START_REF] Francescato | Résistance de plaques multiperforées: comparaison calcul-expérience[END_REF], masonry walls [START_REF] Sab | Yield design of thin periodic plates by a homogenization technique and an application to masonry walls[END_REF], stone columns reinforced soils [START_REF] Hassen | A homogenization approach for assessing the yield strength properties of stone column reinforced soils[END_REF] whereas the first part of this work [START_REF] Bleyer | A computational homogenization approach for the yield design of periodic thin plates -Part I : Construction of the macroscopic strength criterion[END_REF] deals with thin periodic plates in bending. Different numerical techniques have also been used to solve the corresponding optimization problem. In particular, linear programming (LP) associated to a piecewise linearization of the original local strength criterion has been very attractive due to the efficiency of interior point algorithms to solve LP problems. The extension of these algorithms to a wider class of convex programming problems, namely second-order cone programming (SOCP) enables today to solve limit analysis problems with their original nonlinear criterion very efficiently. SOCP has been notably used in the first part of this work to solve the static as well as the kinematic approach of the auxiliary problem. Although an important amount of work has been dedicated to the numerical determination of macroscopic strength criteria of heterogeneous media, only a few papers, to the authors' knowledge, have been aimed at estimating the ultimate load of global structures made of such macroscopic strength criteria derived from the homogenization procedure. One can mention the work of Milani [START_REF] Milani | Homogenization approach for the limit analysis of out-of-plane loaded masonry walls[END_REF] concerning brick masonry or the following papers [START_REF] Hassen | A homogenization approach for assessing the yield strength properties of stone column reinforced soils[END_REF][START_REF] Maghous | Failure design of jointed rock structures by means of a homogenization approach[END_REF] on geotechnical problems. The small amount of work dedicated to this specific aspect is certainly due to the absence of closed-form expressions for the homogenized yield surfaces, which restrains their use to simple failure mechanisms for analytical applications. As regards numerical approaches, the homogenized yield surface has to be approximated, usually by piecewise linearization, to be dealt with. However, the number of hyperplanes describing the surface with a sufficient accuracy can be quite important, which is not desirable for efficient computations. As regards the specific case of thin plates in bending, which is the scope of the present work, the strength criterion depends on the bending moment only. The yield line method, originally proposed by Johansen [START_REF] Johansen | Yield-line theory[END_REF], is an efficient upper bound kinematic method which considers only rigid mechanisms separated by yield lines where angular rotation discontinuities occur. Analytical upper bound estimates are, therefore, easily available and a numerical implementation using linear triangular finite element is also possible [START_REF] Munro | Yield line method by finite elements and linear programming[END_REF], although there are some inherent difficulties due to mesh dependency [START_REF] Johnson | Mechanism determination by automated yield-line analysis[END_REF][START_REF] Jennings | On the identification of yield-line collapse mechanisms[END_REF]. However, to obtain tight upper bound estimates of plates in bending, dissipation through curvature has also to be taken into account. This requires to use a quadratic interpolation, at least, of the plate velocity field. Some authors proposed to use C 1 -continuous high order finite elements to perform the upper bound approach [START_REF] Capsoni | Limit analysis of plates -a finite element formulation[END_REF][START_REF] Le | Upper and lower bound limit analysis of plates using FEM and second-order cone programming[END_REF] but better estimates have been obtained by the present authors using only C 0 -continuous element [START_REF] Bleyer | On the performance of non-conforming finite elements for the upper bound limit analysis of plates[END_REF], dissipation being produced by curvature as well as angular rotation discontinuities. The aim of this work is, thus, to perform yield design computations of plate structures by adapting these numerical methods to the case of complex anisotropic strength criteria computed from homogenization. It will be shown that the yield line method can easily be extended to these criteria without much difficulties, whereas a specific approximation procedure will be required for a complete finite element upper bound approach. It should be noticed that the proposed methods aim at taking advantage of the efficiency of SOCP solvers, which enable to manipulate nonlinear strength criteria, so that a more efficient approximation procedure than piecewise linearization is possible.

Section 2 will first be devoted to the extension of the yield line method to complex strength criteria and an analytical example on a simply supported rectangular plate made of a reinforced material will be presented. Section 3 will present a procedure to approximate a numerically computed three-dimensional yield surface by a convex hull of ellipsoids so that a finite element kinematic approach can be formulated and treated by SOCP solvers. Finally, numerical examples making use of some macroscopic strength criteria previously computed in Part I [START_REF] Bleyer | A computational homogenization approach for the yield design of periodic thin plates -Part I : Construction of the macroscopic strength criterion[END_REF] will be investigated.

2.

A first attempt at evaluating the bearing capacity of heterogeneous thin plates by the yield line method

Yield line method for a numerically computed strength criterion

The yield line method is a simple upper bound approach for the yield design of plates in bending which considers only rigid mechanisms separated by yield lines, where jumps of angular velocity have to be taken into account in the expression of the maximum resisting work. The formulation of the corresponding upper bound yield design problem reads as :

Q ∈ Λ =⇒ ∀ u K.A.
Π([[θ n ]]; n) = sup M ∈G M nn [[θ n ]] = Π(χ = [[θ n ]]n ⊗ n)
We now consider that the local strength criterion of the plate is a macroscopic strength criterion G hom obtained from a homogenization procedure, as described in the first part of this work. In particular, it is described by its support function Π hom (χ). Using the previous remarks, we have :

Π([[θ n ]]; n) = Π hom ([[θ n ]]n ⊗ n) = |[[θ n ]]|Π hom (n ⊗ n)
since support functions are positively 1-homogeneous. Let (e x , e y ) be the orthonormal frame attached to the periodic unit cell, then,

Π([[θ n ]]; n) = |[[θ n ]]|Π hom (cos 2 α(e x ⊗ e x ) + sin 2 α(e y ⊗ e y ) + sin 2α(e x s ⊗ e y )) = |[[θ n ]]|Π 0 (α)
where α is such that n = cos αe x + sin αe y . Therefore, the support function of rotation discontinuities is entirely described by function Π 0 (α) depending on the sole normal orientation angle α. This function can be determined by solving a series of auxiliary yield design problems attached to the unit cell with a macroscopic curvature of the form : χ xx = cos 2 α, χ yy = sin 2 α and χ xy = sin 2α/2 for different values of α. Figure 1 represents such a function Π 0 (α) corresponding to the reinforced plate example presented in the first part of this work is represented. The dependance of this function with respect to α is characteristic of the reinforced plate anisotropy.

Finally, once the geometry of the yield line mechanism has been fixed, Π 0 (α) can be computed for all yield line normal orientations α by interpolating the numerically computed discrete values, deduced from the auxiliary problem.

Application to a reinforced rectangular plate under uniform loading

Consider a rectangular plate of length b in the e X direction and a in the e Y direction. The plate strength criterion corresponds to the macroscopic criterion G hom of the reinforced plate example considered in Part I. The reinforcements, oriented along the unit cell direction e y , form an angle θ with the e Y direction of the plate (see figure 2). The plate is simply supported and loaded by a uniform transversal load -qe Z . The objective of this subsection is to determine an upper bound of the ultimate load q * of such a reinforced plate by considering two different simple yield line mechanisms. For more details on the yield line method, the reader may refer to [START_REF] Save | Plastic limit analysis of plates, shells, and disks[END_REF].

Mechanism 1. The first mechanism (figure 3) consists of two triangles BF C and AED (numbered 2 and 4, respectively) parametrized by the same angle β and two trapezes CDEF and ABF E (numbered 1 and 3). The central yield line EF is animated by a uniform vertical velocity -f e Z with f > 0. The following upper bound estimate is obtained for the reinforced plate limit load (cf. Appendix A):

q * ≤ q 1 = min β∈[0;arctan(b/a)] 4 sin 2β (Π 0 (β -θ) + Π 0 (π -β -θ)) + 4(b/a -tan β)Π 0 (π/2 -θ) ab( 1 2 -a 6b tan β)
Mechanism 2. The second mechanism (figure 4) is similar to the first one, but with the triangles along the plate length and the trapezes along its width. The expression of the upper bound is very similar to the expression obtained for the first mechanism if one exchanges the role of a and b and replaces β by π/2 -β ′ (with a special care for the contribution of the yield line E ′ F ′ in the expression of P rm ). The following upper bound is then obtained :

q * ≤ q 2 = min β ′ ∈[0;arctan(a/b)] 4 sin 2β ′ (Π 0 (π/2 -β ′ -θ) + Π 0 (π/2 + β ′ -θ)) + 4(a/b -tan β ′ )Π 0 (π -θ) ab( 1 2 -b 6a tan β ′ )
The best upper bound is given by the minimum of q 1 and q 2 : q * ≤ q u = min(q 1 , q 2 ) Both minimization problems are numerically solved for different values of the reinforcement orientation θ for a square plate a = b = 1 as well as for a rectangular plate with a = 1 The obtained upper bounds are normalized by a numerical estimate of the corresponding unreinforced plate problem limit load, namely q 0 = 25.02 for the square plate and q 0 = 17.69 for the rectangular plate.

Results for the reinforced square plate are presented in figure 5. These results clearly illustrate the anisotropy effect on the ultimate load, induced by the privileged direction of the reinforcement. Contrary to a isotropic homogeneous square plate for which the optimal mechanism correspond to yield lines along both diagonals, optimal yield line layouts and associated upper bound estimates vary with respect to the orientation angle. For instance, the upper bound estimate between the case θ = 0 • and θ = 45 • is decreased by 18%. It is also to be noted that, for angles comprised in the range 0 • ≤ θ ≤ 45 • , optimal mechanisms are of type 1, whereas for 45 • ≤ θ ≤ 90 • , optimal mechanisms are of type 2, such that the curve of ultimate load estimates is symmetric with respect to θ = 45 • .

As regards the rectangular plate problem (figure 6), the anisotropy effect is still more pronounced, since the upper bound estimate is decreased by 31% between θ = 0 • and θ = 56 • for which the reinforcement orientation yields the minimum limit load estimate. On this example, the competition between mechanisms 1 and 2 is more complex since optimal mechanisms of type 2 are obtained for values in the range 56 • ≤ θ ≤ 74 • .

Finite element formulation of the yield line method

Yield line method for complex strength criteria can be also easily implemented in a finite element framework. Consider, indeed, a triangular mesh and a linear interpolation of the velocity in each element. These elements are pure yield line elements since dissipation can be only produced by rotation discontinuities through element edges. Let {U} denote the global vector of nodal velocities and {F } the global loading vector. The normalization of the work of external loads reads as : P ext (U) = F {U} = 1. Let N D denote the number of active edges i.e. all external and internal edges except free edges and simply supported edges (which do not contribute to the maximum resisting work). For a given active edge j, its length l j and its normal of orientation angle α j are computed. The value Π 0j = Π 0 (α j ) is computed for all edges and the rotation discontinuity is given by [[θ nj ]] = ∆Θ j {U} such that the maximum resisting work is given by :

P rm (U) = N D j=1 l j Π 0j | ∆Θ j {U}|
Finally, the associated linear programming problem reads as :

q Y L = min U,t j N D j=1 l j Π 0j t j s.t. F {U} = 1 | ∆Θ j {U}| ≤ t j j = 1, ..., N D (1)
Although the yield line approach may be very attractive due to its simplicity of implementation, it suffers from many drawbacks. The first one is that it is highly dependent on the mesh orientation because rotation discontinuities can occur along the element edges only. Limit load upper bounds produced by a pure yield line finite element computation can be, therefore, highly dependent on the mesh orientation. Different approaches have been proposed in the literature to overcome this drawback [START_REF] Johnson | Mechanism determination by automated yield-line analysis[END_REF][START_REF] Jennings | On the identification of yield-line collapse mechanisms[END_REF][START_REF] Askes | Adaptive analysis of yield line patterns in plates with the arbitrary lagrangian-eulerian method[END_REF]. Unfortunately, it has been pointed out by Braestrup [22] that, even if exact solutions can be obtained for some problems, yield line theory fails, in general, to predict the exact limit load even when assuming a very complex mechanism or with an infinitely refined mesh. Indeed, for the simply supported square plate under uniform pressure with a von Mises strength criterion, the yield line theory cannot do better than predict an upper bound which is 10% higher than the exact limit load. This is due to the fact that the exact solution exhibits areas where dissipation is produced by curvature instead of rotation discontinuities only.

3. Upper bound yield design approach with a complex strength criterion

Numerical challenges

In the previous section, it has been highlighted that the yield line method can be easily extended to the case when the plate strength criterion is complex, provided that function Π 0 (α) can be computed for all potential yield lines orientations. However, dissipation through curvature has also to be taken into account to obtain better upper bound estimates. Even if a preferred orientation for the rotation discontinuities will still be present due to the fixed positions of the finite element edges, it will be compensated by a localized zone of curvature field so that a finite element approach combining curvature and rotation discontinuities, besides being much more accurate, is also much less sensitive to the mesh layout than a pure yield line approach.

Unfortunately, the yield design problem may become very difficult when dealing with strength criteria of general shape such as those obtained from a homogenization procedure. For example, if the numerical value of the support function is known only on a finite set of curvature directions, the maximum resisting power cannot be analytically expressed in terms of the unknowns of the problem. The main challenge is, thus, to obtain easy to handle semi-analytical expressions for the support function, tractable by a numerical solver.

Piecewise linear approximation of the strength criterion is one possible approach. Although it was mainly used to approximate classical non-linear strength criterion in the early developments of numerical yield design, one can imagine using this technique to approximate a complex strength criterion obtained from a homogenization procedure and solve the associated linear programming (LP) problem. However, the number of planes required to approximate the strength criterion with a sufficient accuracy may become rapidly large.

The approximation procedure proposed in this work is a further step in this direction, which combines two important ingredients : first, a natural idea is to use primitives of higher order than planes so that a smaller number of these will be required in the approximation process; the second ingredient is to obtain an approximate criterion which can be expressed using conic constraints, in order to take maximum advantage of the efficiency of SOCP optimization solvers.

Based on these two arguments, our goal is to approximate a convex strength criterion using a convex hull of ellipsoids, since the support function of an ellipsoid E i can, indeed, be written as :

π E i (χ) = J i • χ + µ i • χ
where χ = {χ xx χ yy 2χ xy } T , J i is a 3 × 3 positive semi-definite upper triangular matrix (defined by the ellipsoid axes length and orientation) and µ i is the coordinate vector of the ellipsoid center. It is to be recalled that the support function of the convex hull of two sets A and B is equal to the maximum between the support function of A and that of B. Hence, the support function of the convex hull CH(E i ) of r ellipsoids E i for i = 1, ..., r parametrized by J i and µ i is given by :

π CH(E i ) (χ) = max i=1,...,r { J i • χ + µ i • χ}
which shows that the support function of the convex hull of r ellipsoids can be easily expressed by means of conic constraints.

Approximation procedure of the strength criterion

In a recent work, the authors proposed a relatively simple algorithm to approximate numerically computed strength criteria using a convex hull of ellipsoids. The general principle of the algorithm will be briefly recalled here, more details may be found in [? ].

Assuming that the support function of the strength criterion is known for M directions, uniformly distributed on the unit sphere in the space of curvature, the following iterative algorithm constructs, at each iteration, an ellipsoid which is a local approximation of the original convex set.

1. Let χ 0 be the direction whose tangent plane is the furthest away from the current approximation G app,r ; 2. a second-order Taylor expansion of Π(χ) is performed around χ 0 ; 3. the local radii of curvature of G at χ 0 are computed from the Taylor expansion; 4. a series of ellipsoid with the same curvature and tangent to G at χ 0 are considered; 5. the best ellipsoid is the one which minimizes the root mean squared error between the support function of G and the one of the ellipsoid in the region around χ 0 ; 6. the optimal ellipsoid E n is added to the convex hull of the current approximation G app,r+1 = CU(G app,r , E r ); 7. go to step 1.

This approximation procedure has been applied to the macroscopic strength criterion of the reinforced plate problem. In figure 7, the original strength criterion and approximate strength criteria obtained with different number of ellipsoids are represented to assess the performance of the proposed procedure. More precisely, figure 8 represents the evolution of the root mean square (RMS) and maximal error made by the approximation procedure as a function of the number r of ellipsoids. It can be observed that the error is rapidly decreasing. Maximum error is of 17% for 10 ellipsoids, 5% for 30 ellipsoids and 2.8% for 50 ellipsoids. The determination of the 50 optimal ellipsoids took less than 2 minutes. Finally, it is worth noting that, although inner approximations are initially produced by the algorithm, it is always possible to expand them by an appropriate scaling factor so as to obtain outer approximations of the initial convex set. This was done for the approximating criterion involving 50 ellipsoids.

Finite element and SOCP formulation

In this subsection, the discrete upper bound kinematic approach on a plate structure, the strength criterion of which has been previously approximated by a convex hull of r ellipsoids, will be presented. The plate is discretized into N E cubic Hermite triangles (H3) which were introduced in [START_REF] Bleyer | On the performance of non-conforming finite elements for the upper bound limit analysis of plates[END_REF] and successfully used to solve the kinematic approach of the auxiliary problem in Part I.

The curvature at a given Gauss point g in a given element e can be expressed as {χ} = [B e,g ]{U}. The support function at this point is then :

π e,g CH(E i ) = max i=1,...,r { [J i ][B e,g ]{U} + µ i [B e,g ]{U}}
The contribution of the curvature term to the maximum resisting work is then obtained after integration over all Gauss points of the structure : where c e,g are constant terms coming from the n-Gauss points quadrature.

P curv rm = N E e=1 n g=1 c e,g max i=1,...,r { [J i ][B e,g ]{U} + µ i [B e,g ]{U}}
As regards the contribution of inter-element rotation discontinuities to the maximum resisting work, the quantities Π 0j are directly evaluated as a function of the approximating ellipsoids parameters. Indeed, introducing χ disc,j = {cos 2 α j sin 2 α j sin 2α j } T :

Π 0j = π CH(E i ) (χ disc,j ) = max i=1,...,r { J i • χ disc,j + µ i • χ disc,j } so that P disc rm = N D j=1 Π 0j m g ′ =1 c ′ j,g ′ | ∆Θ j {U}|
where c ′ j,g ′ are constant terms coming from the m-Gauss points quadrature on edges.

Finally, introducing different auxiliary variables, the upper bound kinematic problem can be formulated as :

q u = min N E •n k=1 c k t k,0 + N D •m k ′ =1 Π 0j c ′ k ′ s k ′ s.t. F {U} = 1 {r k,i } = [J i ][B e,g ]{U} t k,i ≥ {r k,i } t k,i + µ i [B e,g ]{U} ≤ t k,0 u k ′ = ∆Θ j {U} s k ′ ≥ |u k ′ | i = 1, ..., r k = 1, ..., n • N E k ′ = 1, ..., m • N D (2)
which is a standard SOCP problem involving r • n • N E conic constraints, r • n • N E + m • N D linear inequality constraints and 1 + 3r • n • N E + m • N D equality constraints.
It is worth noting that, even with a relatively simple mesh of N E = 500 elements, a n = 3 Gauss points rule and r = 50 ellipsoids, r • n • N E = 75 000 so that (2) already constitutes a large scale SOCP problem.

Numerical examples

In this section, numerical examples are presented to assess the performance of the different proposed methods. In the first example, problem (2) associated to the reinforced plate problem is solved using the Mosek software package and the influence of the number r of approximating ellipsoids is investigated. The value of r is, indeed, the main factor which determines the difficulty of the optimization problem. It is to be noted that approximations with r = 10 and r = 30 are not outer approximations of the macroscopic criterion G hom so that a global computation with such domains will not produce an upper bound of the limit load. Therefore, problem (2) has been solved for these values of r, then, in a second step, the so-obtained optimal velocity field {U} has been used to evaluate the maximum resisting work with r = 50 ellipsoids, the convex hull of which has been considered to be sufficiently close to the original strength criterion. Such a post-processing procedure has been applied to all SOCP upper bounds presented in this section. The numerical yield line problem (1) has also been solved using Mosek for the reinforced plate problem. The second example will treat the case of a perforated circular plate and compare the influence of the hole sizes and shapes on the limit load.

Reinforced plate under uniform loading

Finite element meshes used for the SOCP computations are represented in figure 9(a) for the square plate problem and in figure 9(b) for the rectangular plate problem. As regards Nb. of ellipsoids Square plate Rectangular plate r = 10 2.5s 4s r = 30 6.5s 13s r = 50 12s 24s

Table 1: Typical optimization times for the SOCP problems on the square and rectangular plate problem.

(performed on a Intel-P4 2.4 GHz running Linux 32-bits with Mosek v7.0) numerical yield line computations, the same type of structured mesh has been used with 1568 elements for the square plate and 3248 elements for the rectangular plate.

In figure 10 for the square plate and in figure 11 for the rectangular plate, the SOCP upper bounds obtained with r = 10, 30 and r = 50 ellipsoids have been compared to the upper bound obtained with the analytical yield line method of subsection 2.2 as well as to the upper bound obtained with a numerical implementation of the yield line method presented in subsection 2.3. It is first to be noted that numerical yield line upper bound does not necessarily improve the analytical yield line upper bound. This is mainly due to the considered finite element mesh which restrains the potential set of yield lines. In this case too, the numerical yield line method also proved to be highly sensitive to the mesh layout.

Secondly, SOCP computations yield better estimates, by almost 10% in some cases. One can also observe that upper bounds obtained with r = 30 or 50 ellipsoids are very close to those obtained with r = 10 ellipsoids, differing by less than 2%. Therefore, despite the fact that the original strength criterion is not particularly well approximated with r = 10 ellipsoids, it seems that, due to the post-processing procedure, it is sufficient to obtain a good optimal velocity field and related upper bound estimate. Besides, typical computation times for the optimization procedure have been reported in table 1. It can be observed that the computational cost on a standard desktop computer is very reasonable even for the highest number of ellipsoids used to describe the yield surface.

Finally, failure mechanisms of the square plate problem have also been represented in figure 12 for the homogeneous plate and the reinforced plate with θ = 0 • and θ = 45 • . The anisotropy induced by the presence of the reinforcements can be clearly observed for these two cases by comparison to the homogeneous case. For example with θ = 45 • , the optimal mechanism seems exhibit a kind of yield line along the diagonal parallel to the reinforcement direction and a region with a distributed curvature along the other diagonal.

Perforated circular plate under uniform loading

In this last example, the case of a simply supported circular plate of radius R subject to a uniform loading of intensity q is considered. The plate is made of a von Mises material of ultimate bending moment m p perforated by a series of circular or square holes arranged as depicted in figure 13. Different square unit cells are considered by varying the non-dimensional hole size parameter λ which was respectively taken as λ = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8 for the circular hole as well as for the square hole. Obviously, for λ = 0, the solution relative to the homogeneous circular plate problem obeying a von Mises criterion is recovered. For simple supports, the limit load is known to be q 0 = 6.52m p /R 2 [START_REF] Save | Plastic limit analysis of plates, shells, and disks[END_REF]. Besides, for λ = 1, all holes are connected and the plate bearing capacity vanishes which gives then q = 0 as the limit load.

The macroscopic strength criterion G hom has been determined by solving the kinematic upper bound auxiliary problem for the 2 × 7 different unit cell geometries. For one unit cell geometry, 10 000 values of Π hom have been computed in around 1 hour. Simultaneous computations have been made to reduce the computational cost. Then, the approximation procedure has been performed to approximate the strength criterion by 50 ellipsoids. This also took around 2 minutes for each yield surface. Finally, computations on the homogenized problem have been performed on a quarter of circular plate (see figure 14) with 10 ellipsoids and the previously mentioned post-processing step. Each optimization procedure was performed by Mosek in 5-10 seconds.

In figure 15, the ratio between the upper bounds of each homogenized problems and the solution for the homogeneous plate have been represented with respect to the plate porosity defined by η = πλ 2 /4 in the case of circular holes and η = λ 2 in the case of square holes. The black dotted line represents the Voigt upper bound which is obtained by considering that the exact macroscopic strength criterion is contained in the von Mises criterion with a reduced moment capacity m p,red = (1 -η)m p , which leads to an upper bound equal to q V oigt = (1 -η)q 0 . One can observe that the Voigt upper bound gives a very poor estimate of the limit load for all types of hole shape. Moreover, for a given value in the small porosity range, the limit load does not seem to depend on the hole shape. Finally, considering the evolution of the limit load estimates with respect to the hole size λ, a maximum discrepancy of 6.3% can be observed for λ = 0.5 between the circular and the square shape.

Conclusions and perspectives

In the second part of this work on yield design homogenization method applied to periodic thin plates, the computation of limit load estimates of homogenized plate structure problems derived from upper bound kinematic approaches has been addressed, and associated numerical tools have been proposed. It has been shown, in particular, that the yield line method, which is commonly used in civil engineering computations, can be easily extended to the case of complex strength criteria obtained from homogenization procedures. Indeed, taking advantage of the 1-homogeneity property of support functions and provided that the function Π 0 (α) depending on the sole normal orientation angle α is determined from the homogenization procedure, the yield line method can be applied, in a straightforward manner, analytically as well as numerically, as shown on the example of the simply supported reinforced rectangular plate. A more sophisticated and performing upper bound kinematic approach has also been proposed based on the use of finite elements taking into account dissipation through curvature as well as angular jumps. For the curvature term, the support function has no simple ex- pression and an approximation procedure is, thus, required to implement a semi-analytical expression of the support function in the optimization solver. An approximation algorithm based on the construction of a convex hull of ellipsoids is proposed leading to a simple description of the macroscopic strength criterion support function which is compatible with a SOCP formulation. Despite a large number of resulting conic constraints in the global optimization problem, accurate solutions can be obtained in extremely reasonable computing times on a standard desktop. Results indicate that good optimal velocity fields and accurate upper bounds can still be found with a small number (e.g. 10) ellipsoids, a post-processing procedure with a more accurate description of the strength criterion allowing to obtain a more precise upper bound estimate. Finally, the numerical results indicate that refined kinematic approaches with approximate criteria and SOCP provide better upper bound by at most 15% compared to the yield line method. Therefore, as regards civil engineering applications, the yield line method seems to be a good compromise between implementation simplicity and predictability, although mesh dependency can be an important drawback.

Possible extensions of this work are multiple. A first direction of research may consist in applying the general idea of approximating complex strength criteria obtained from homogenization procedures to other mechanical models (plane strain, soils,...). The development of more performing approximation procedures which would require even less primitives to describe the yield surface is also possible. As regards periodic plates, a current work aims at treating the case of the plate thickness being of the same order as the unit cell characteristic length, by solving the auxiliary problems defined on a three-dimensional unit cell with the same numerical tools. Finally, further work can also deal with the case of plates having a finite resistance to shear forces or the more complicated case of shells.

Figure 1 :

 1 Figure1: Function Π 0 (α) for the reinforced plate example[START_REF] Bleyer | A computational homogenization approach for the yield design of periodic thin plates -Part I : Construction of the macroscopic strength criterion[END_REF]. The circles correspond to numerically computed values obtained from the resolution of auxiliary yield design problems on the unit cell. The horizontal dotted line corresponds to the value of Π 0 (α) in the case of an unreinforced von Mises plate with ultimate bending moment m p1 = 1.

Figure 2 :Figure 3 :

 23 Figure 2: Homogenization procedure for a rectangular reinforced plate under uniform loading

6

 6

Figure 4 :

 4 Figure 4: Mechanism 2 for the rectangular plate under uniform loading

Figure 5 :

 5 Figure 5: Yield line upper bound for the reinforced square plate as a function of the reinforcement orientation θ. Optimal yield line mechanisms are represented for θ = 10 • , θ = 45 • and θ = 80 •

Figure 6 :

 6 Figure 6: Yield line upper bound for the reinforced rectangular plate as a function of the reinforcement orientation θ. Optimal yield line mechanisms are represented for θ = 45 • , θ = 60 • and θ = 90 •

Figure 7 :

 7 Figure 7: Approximation procedure of the macroscopic strength domain G hom for the reinforced plate

Figure 8 :

 8 Figure 8: Evolution of the relative error with the number of ellipsoids. Errors are computed by considering the difference between both support functions.

Figure 9 :

 9 Figure 9: Finite element meshes for SOCP computations

Figure 10 :

 10 Figure 10: Comparison of the upper bound estimates for the reinforced rectangular plate obtained with various methods : analytical yield line (YL) method, numerical yield line method, SOCP with a convex hull of r = 10, 30 and 50 ellipsoids

Figure 11 :

 11 Figure 11: Comparison of the upper bound estimates for the reinforced rectangular plate obtained with various methods : analytical yield line (YL) method, numerical yield line method, SOCP with a convex hull of r = 10, 30 and 50 ellipsoids

Figure 12 :

 12 Figure 12: Comparison of different failure mechanisms for the simply supported square plate

Figure 13 :

 13 Figure 13: Circular plate perforated by a series of circular or square holes

Figure 14 :

 14 Figure 14: 413 elements mesh used for the global computation

Figure 15 :

 15 Figure 15: Evolution of the obtained upper bounds with respect to the plate porosity

  with q, P ext ( u) ≤ P rm ( u)

	where P ext ( u) is the work of external loads in the kinematically admissible velocity field u and P rm ( u) = Γ Π([[θ n ]]; n)dl is the maximum resisting work associated with a set of yield lines Γ of unit normal n and angular velocity jumps [[θ n ]] across Γ following the normal n. Hence, Π([[θ n ]]; n), which corresponds to a particular value of the support function of the strength criterion associated with rotation discontinuities, is defined as :

Appendix A. Analytical yield line computations

The work of external loads developed in mechanism 1 is given by :

The angular velocity of regions 1 and 3 are given by θ 01 = θ 03 = 2 f a , whereas the angular velocity of regions 2 and 4 are given by θ 02 = θ 04 = 2 f a tan β . The hodograph of angular velocities represented in figure A.16 is used to compute the discontinuity of angular velocity θ ij between regions i and j. We have :

The orientation angles α ij of the yield line normals n ij are given by : associated with this mechanism is given by :

Finally, observing that Π 0 (α) is π-periodic and symmetric with respect to π/2 due to the unit cell symmetry, the following upper bound estimate is obtained for the reinforced plate limit load :