Michele Godio
email: michele.godio@enpc.fr

Ioannis Stefanou

Karam Sab

Jean Sulem

Multisurface plasticity for Cosserat materials: plate element implementation and validation

Keywords: Multisurface plasticity, Cosserat continuum, finite elements, micropolar plates, strain localisation

The macroscopic behaviour of materials is affected by their inner micro-structure. Elementary considerations based on the arrangement, and the physical and mechanical features of the micro-structure may lead to the formulation of elastoplastic constitutive laws, involving hardening/softening mechanisms and nonassociative properties. In order to model the non-linear behaviour of micro-structured materials, the classical theory of time-independent multisurface plasticity is herein extended to Cosserat continua. The account for plastic relative strains and curvatures is made by means of a robust quadratic-convergent projection algorithm, specifically formulated for non-associative and hardening/softening plasticity. Some important limitations of the classical implementation of the algorithm for multisurface plasticity prevent its application for any plastic surfaces and loading conditions. These limitations are addressed in this paper, and a robust solution strategy based on the Singular Value Decomposition technique is proposed. The projection algorithm is then implemented into a finite element formulation for Cosserat continua. A specific finite element is considered, developed for micropolar plates. The element is validated through illustrative examples and applications, showing able performance.

INTRODUCTION

Trabecular bones, wood, soils, metals and many other materials existing in nature are provided with an internal micro-structure which is apparent at specific observation scales. The macroscopic response of these materials is therefore affected by the geometrical (characteristic length) and mechanical internal properties of their micro-structure. This holds true for man-made materials and structures such as composites, grid-works, masonry, etc.

An efficient and elegant way to model the overall response of such materials without the need to resort to non-local theories is represented by the Cosserat continuum. With respect to the classical (Cauchy) continuum, the Cosserat (or micropolar) continuum has independent rotational degrees of freedom attached to the material particle. This confers an internal length to the medium, and makes it suitable for the description of micro-structured materials. The additional degrees of freedom also provide a larger set of boundary conditions for the medium. It is well known how the presence of an internal length and the use of additional boundary conditions may affect the material response, especially in the non-linear range. Boundary layers and regions with concentrated irreversible deformation such as shear bands are controlled by internal material lengths [START_REF] Muhlhaus | The thickness of shear bands in granular materials[END_REF][START_REF] Vardoulakis | Bifurcation analysis in geomechanics[END_REF][START_REF] Sulem | Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure[END_REF]. The use of the Cosserat continuum is hence crucial to model complex test configurations of micro-structured materials.

The use of Cosserat continua for representing the macroscopic behaviour of micro-structured materials is well assessed nowadays, at least for what concerns the elastic behaviour [START_REF] Anderson | Size effects due to Cosserat elasticity and surface damage iin closed-cell polymethacrylimide foam[END_REF][START_REF] Forest | Cosserat overall modeling of heterogeneous materials[END_REF][START_REF] Stefanou | Three-dimensional Cosserat homogenization of masonry structures: Elasticity[END_REF]. Most of the models proposed in the literature are devoted to the formulation of homogenization schemes, i.e. mathematical or numerical procedures through which one can retrieve the micropolar properties of the material by starting from those of its micro-structure. However, there is still a lack of knowledge concerning the modeling of the inelastic behaviour of these materials. In the non-linear range, irreversible, dissipative, time-dependent phenomena relative to the evolution and damage of the material's inner structure occur. In this context, the formulation of models based on micropolar continua seems complicated, if not unfeasible. However, based on elementary considerations and often motivated by micromechanical observations (such as friction and interlocking phenomena in soils, or creep phenomena in cement materials), researchers may be led to the formulation of extended constitutive models for inelasticity. It has been the case, for instance, for the generalised J2-plasticity model for soils [START_REF] Muhlhaus | The thickness of shear bands in granular materials[END_REF], which has made possible the study and the understanding of complex phenomena such as the localisation of deformation in microstructured materials [START_REF] Steinmann | Localization within the Framework of Micropolar Elasto-Plasticity[END_REF][START_REF] De Borst | Localisation in a Cosserat continuum under static and dynamic loading conditions[END_REF].

In this paper, a numerical solution procedure for the modeling of the elastoplastic behaviour of micro-structured materials is presented in the framework of the Cosserat continuum. Notice that problems related to the non-linear material behaviour always require the use of numerical tools. Concerning plasticity, some have gained recent attention. For instance, one may mention the use of mathematical programming techniques for the solution of incremental elastoplastic problems [START_REF] Krabbenhø Ft | Formulation and solution of some plasticity problems as conic programs[END_REF][START_REF] Krabbenhø Ft | An interior-point algorithm for elastoplasticity[END_REF]. Those techniques have been also employed for the solution of limit analysis problems [START_REF] Makrodimopoulos | Lower bound limit analysis of cohesive-frictional materials using 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 numerical procedure adopted in this paper is based on a projection algorithm, originally formulated for classical hardening-plastic materials by [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF]. The projection algorithm is a robust quadratic-convergent algorithm which is widely used in finite elements for the solution of elastoplastic boundary value problems. The algorithm is based on a return map, which solves the incremental elastoplastic problem at every increment of the load-path (and point of the model), by assuring the respect of the consistency requirement. The solution will be then found at the projection point of the given increment, falling on the boundary of the elastic domain of the material [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF]. Herein, we made use of the time-independent version of the projection algorithm. In other words, viscoplastic effects are not considered. In this setting, the most general version of the algorithm for the multisurface plasticity theory is adopted [START_REF] Simo | Non-smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms[END_REF].

The theory of time-independent multisurface strain hardening plasticity is very general [START_REF] Simo | Computational inelasticity[END_REF]. It allows to define multiple sets of plastic surfaces, non-associated plastic potentials, and softening/hardening laws of every kind. When interested in computing not only the ultimate (or limit) load, but the whole material response until failure, this theory gives great flexibility to the user in: a) prescribing failure criteria and b) choosing appropriate plastic evolution laws for the material. Nevertheless, the numerical implementation of the multisurface plasticity theory in finite elements suffers from important limitations, which prevent its application for any plastic surfaces and loading conditions. These limitations are explored in this paper, and a solution strategy based on the use of the Singular Value Decomposition (SVD) technique is proposed. This represent an important improvement in the numerical implementation of the algorithm for multisurface plasticity that can be used not only in the framework of the Cosserat continuum theory, but also in more classical finite element formulations.

In this paper, the theory of multisurface plasticity is adapted for Cosserat continua. The onset and the evolution of irreversible deformations within the micro-structured material is here detected in terms of plastic relative strains and curvatures. It is important to note that the use of projection algorithms for plasticity has been already encountered in the frame of the Cosserat continuum theory. Most of the existing works were devoted to the study of strain localisation. As known, the Cosserat continuum provides a natural regularization effect versus the pathological mesh-sizedependency intrinsic to classical numerical formulations. To this purpose, one can see the works of de Borst and co-workers [START_REF] De Borst | Localisation in a Cosserat continuum under static and dynamic loading conditions[END_REF][START_REF] De Borst | Continuum models for discontinuous media[END_REF][START_REF] De Borst | Simulation of strain localisation: a reappraisal of the Cosserat continuum[END_REF][START_REF] De Borst | Fundamental issues in finite element analyses of localization of deformation[END_REF], but also [START_REF] Papanastasiou | Numerical treatment of progressive localization in relation to borehole stability[END_REF][START_REF] Li | A consistent return mapping algorithm for pressure-dependent elastoplastic Cosserat continua and modelling of strain localisation[END_REF][START_REF] Neff | a Numerical Solution Method for an Infinitesimal Elasto-Plastic Cosserat Model[END_REF], among others. Always in the frame of the Cosserat continuuum theory, [START_REF] Addessi | A 2D Cosserat finite element based on a damage-plastic model for brittle materials[END_REF] proposed a coupled damage-plastic model for microstructured materials showing brittle macroscopic behaviour. In that case, a projection algorithm was used to integrate simultaneously the evolution of the plastic and damage variables, see also [START_REF] Steinmann | Theory and numerics of ductile micropolar elastoplastic damage[END_REF]. Another use of the projection algorithm is found in [START_REF] Khoei | An enriched-FEM model for simulation of localization phenomenon in Cosserat continuum theory[END_REF], where an elastoplastic Cosserat model was incorporated into an extended finite element (X-FEM) formulation for the propagation of cracks. Nevertheless, in the aforementioned works it has always been considered a single plastic surface. In the present paper, the novelty is to consider a set of plastic surfaces which are multiple in number, of general form and intersecting. One example of Cosserat multisurface plasticity has been given by [START_REF] Godio | Cosserat Elastoplastic Finite Elements for Masonry Structures[END_REF], where an application to masonry structures is proposed as an extension to more classical approaches based on the Cauchy continuum [START_REF] Stefanou | Three dimensional homogenization of masonry structures with building blocks of finite strength: A closed form strength domain[END_REF].

The multisurface plasticity theory is here implemented into a Cosserat finite element and therefore incorporated within an implicit time-integration scheme for the non-linear analysis of microstructured materials. A Newton-Raphson iterative method is employed for the solution of the elastoplastic problem both at the local level, i.e. on the integration point of the element, and at the global (element) level. At the local level, the Newton's method is used in conjunction with the SVD matrix inversion technique. This enables the projection algorithm to attain all the features of the multisurface plasticity theory, i.e. the definition of an unlimited number of plastic surfaces and potentials of general form. At the global level, the method requires the implementation of the algorithmic version of the elastoplastic stiffness tangent matrix in order to assess the full rate of convergence. The finite element that we use is a 8-node quadrilateral, that accommodates small displacements and small rotations at every node [START_REF] Godio | Dynamic finite element formulation for Cosserat elastic plates[END_REF]. The element is formulated for a Cosserat (or micropolar) plate model [START_REF] Eringen | Theory of micropolar plates[END_REF][START_REF] Green | Theory of an elastic-plastic Cosserat surface[END_REF]. Such model enables the analysis of materials and structures with an inner micro-structure mainly developed in two directions. Although this may certainly represent a simplification with respect to the three-dimensional case, a plate model allows the representation of complex spatial configurations in a quite refined way, with the advantage of being much less computational demanding. To this purpose, notice that the plate element is stress resultant, which means that no integration over the thickness is made. Moreover, small deflections and small angle approximations are involved in the formulation. In this sense, we are not dealing with a geometrically exact plate formulation.

The paper is structured as follows. Section 2 provides basic background in Cosserat (or micropolar) plate models. In Section 3 we present the multisurface plasticity theory for microstructured materials within the frame of the Cosserat continuum. In this case, we will make use of general notations. In Section 4 the projection algorithm is described, and details on the implementation of the plasticity theory in finite elements are provided in Section 5. Finally in Section 6 we focus more specifically on the micropolar plate element. We present numerical tests for the assessment of the projection algorithm accuracy and benchmarks for the validation and the actual use of the element. These benchmarks concern both the in-plane (or membrane) and the outof-plane (or flexural) behaviour of the element. It is worth noticing that numerical benchmarks for Cosserat (multisurface) plasticity are nearly absent in literature. Those proposed in Section 6 prove then useful for the validation of all kind of algorithms and finite elements formulated for micropolar plasticity. The benchmarks for the in-plane behaviour, in particular, apply for 2D Cosserat continua.

Matrix notation is adopted throughout the paper. We use symbol d[ ] to indicate an increment between two successive load steps, and [ ] to denote time derivative. Symbol ∂ V [ ] is used for vector derivative.

THE COSSERAT (OR MICROPOLAR) PLATE MODEL

The Cosserat continuum (or micropolar) theory differs from the classical theory of the Cauchy continuum in two aspects [START_REF] Eringen | Microcontinuum field theories. I: Foundations and solids[END_REF]. First, the couple-stress is completely absent in the Cauchy continuum. Second, in the Cosserat continuum the stress tensor is not symmetric, as consequence of the couple-stresses. Figure 1 shows the stresses and couple-stresses featuring in a micropolar plate model [START_REF] Godio | Dynamic finite element formulation for Cosserat elastic plates[END_REF][START_REF] Eringen | Theory of micropolar plates[END_REF]. The former are assembled in the vector τ and the latter are contained into the polar (

As far as it concerns the membrane (or in-plane) behaviour (Figure 1-left), one recognizes: the inplane tractions (τ 11 , τ 22 ), the membrane shears (τ 12 , τ 21 ) and the in-plane couple-stresses (µ 31 , µ 32 ). These stress measures are those of 2D Cosserat continuum. In addition one has the components controlling the flexural (or out-of-plane) behaviour of the plate (Figure 1-right), namely: the torsions (µ 11 , µ 22 ), the out-of-plane flexions (µ 12 , µ 21 ), the transverse shears (τ 31 , τ 32 ), and the longitudinal shears (τ 13 , τ 23 ). Due to the presence of transverse shears, the micropolar plate model can be regarded as a thick plate theory. Its convergence towards the thin plate solution has been numerically assessed by [START_REF] Godio | Dynamic finite element formulation for Cosserat elastic plates[END_REF].

τ
In particular one has: the in-plane extensions/compressions (γ 11 , γ 22 ), the membrane shear strains (γ 12 , γ 21 ) and the in-plane curvatures (κ 31 , κ 32 ). Moreover, one has to take into account the out-ofplane curvatures, as the torsion deformations (κ 11 , κ 22 ) and the bending deformations (κ 12 , κ 21 ), together with the transverse and longitudinal shear strains, respectively (γ 31 , γ 32 ) and (γ 13 , γ 23 ).

The deformation measures are derived from the kinematics of the plate, which consists of 3 displacements, indicated with u, and 3 rotations about the orthogonal axis, called Cosserat rotations ω c :

u = u 1 u 2 u 3 t ω c = ω c 1 ω c 2 ω c 3 t . (3)
For plates under small deflections and rotations, one defines the relative strain vector γ as [START_REF] Steinmann | Localization within the Framework of Micropolar Elasto-Plasticity[END_REF]:

γ = ∇u + e • ω c , γ ij = u i,j + e ijk ω c k (4)
and the curvature vector as:

κ = ∇ω c , κ ij = ω c i,j . (5)
Dealing with finite elements, it is useful to express the balance equations in weak form. Herein, we make use of the virtual work principle on a region of the plate A and its boundary ∂A. Following the vector form introduced above and indicating a virtual variation by δ[ ], the principle reads:

R -P -Q = 0, (6)
where:

R = A τ t δγ + µ t δκ dA (7)
is the work of the stresses and couple-stresses,

P = A b t δu + l t δω c dA + ∂A T t δu + M t δω c dA (8)
is the work done respectively by the body forces b and couples l and the traction T and moment M on the boundary, and

Q = A ρü t δu + ρω ct Iδω c dA (9)
is the work due to the inertial mass ρ and second order moment I.

MULTISURFACE PLASTICITY THEORY FOR COSSERAT MATERIALS

The time-independent theory of non-smooth multisurface hardening plasticity is presented in this section. This theory, which was formerly formulated for simple (or Cauchy) materials in its original version [START_REF] Vermeer | Non-associated plasticity for soils, concrete and rock[END_REF] and then extended for the computation of multiple yield surfaces [START_REF] Simo | Non-smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms[END_REF][START_REF] Simo | Computational inelasticity[END_REF], is here used in the framework of the Cosserat continuum theory. Proceeding in general fashion and adopting incremental notation, relative strains (4) and curvatures (5) can be assembled into a general total strain vector dε, defined as follows:

dε = dγ dκ t . (10)
Similarly, the stresses and the couple-stresses can be assembled into a general stress vector dσ:

dσ = dτ dµ t . (11)
The formulation (and the implementation) of the plasticity theory for Cosserat materials is obtained in a way similar to that followed for Cauchy materials [START_REF] De Borst | Simulation of strain localisation: a reappraisal of the Cosserat continuum[END_REF]. By means of the additive decomposition, the vector of the total (small) strains is divided into the elastic [ ] e and plastic [ ] p parts:

dε = dε e + dε p . ( 12
)
The elastic response of the material is governed by the general linear constitutive law, relating stresses and couple-stresses to the elastic strain increment:

dσ = C [dε -dε p ] , ( 13
)
where C is the matrix of the elastic moduli. Dealing with an incremental problem, it is also useful to express the stresses as a function of the total strains:

dσ = C ep dε ( 14
)
where C ep is now the matrix of the elastoplastic tangent moduli.

The hardening/softening variables are introduced in similar way. These are defined by α, and its dual q. The former may be considered as a deformation measure, while the latter as a stress measure. Their relation is based on the definition of D, the matrix of the hardening moduli: dq = -Ddα. [START_REF] Green | Theory of an elastic-plastic Cosserat surface[END_REF] It follows the definition of generalised strains and stresses, useful for the sequel, respectively given by the couples (ε; α) and (σ; q). The salient feature of the multisurface plasticity theory lies on the definition of the elastic domain. Denoted with E σ , the elastic domain is defined as the convex region of the generalised stress space bounded by multiple N F plastic surfaces F β (σ; q) [START_REF] Simo | Non-smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms[END_REF]:

E σ = (σ; q) |F β (σ; q) ≤ 0, ∀β ∈ [1, . . . , N F ] (16)
These surfaces represent the prescribed material yield loci (or yield criteria). They are smooth, but intersect in non-smooth manner (Figure 2). The boundary of E σ is:

∂E σ = (σ; q) |F β (σ; q) = 0, ∀β ∈ [1, . . . , N F ] . ( 17
)
τ μ F β (τ,μ;q) ≤ 0 According to the multisurface plasticity theory, it is a convex region in the generalised stress space (τ , µ; q) defined by multiple yield criteria.

Non-associative plasticity requires the definition of plastic potentials G β (σ; q), and hardening/softening functions H β (σ; q). Through these, one formulates the equations of evolution. On one hand, the flow rule defines the plastic strain increment:

dε p = N F β=1 λβ ∂ σ G β (σ; q) . (18)
On the other hand, the hardening law controls the increment of the hardening variables:

dα = N F β=1 λβ ∂ q H β (σ; q) . ( 19
)
The plastic multipliers λβ give the magnitude of the plastic strain increment and, by extension, that of the hardening variables. They must respect the following restrictions, namely a) the complementary conditions:

F β (σ; q) ≤ 0, λβ ≥ 0, λβ F β (σ; q) = 0, (20)
and b) the consistency requirement:

λβ dF β (σ; q) = 0. (21)
The above restrictions are written for the N act activated surfaces, with N act ≤ N F . They allow one to determine whether the material response to loading will be elastic or elastoplastic. In the presence of multiple surfaces, this is a central task in the formulation of the algorithm for plasticity. The unloading response is considered elastic, and not affected by damage (cf. [START_REF] Addessi | A 2D Cosserat finite element based on a damage-plastic model for brittle materials[END_REF]).

Further notations

The multisurface plasticity theory requires, especially in its algorithmic version, a large number of variables and indices. When used in conjunction with the Cosserat continuum theory, notation may become quite complicated. Vector notation has already been introduced. From now on, it is convenient to introduce the following additional notations.

The set of plastic surfaces, plastic potential and hardening/softening functions will be indicated in the following vector form:

F =    F 1
. . .

F N F    , G =    G 1 . . . G N F    , H =    H 1 . . . H N F    . (22)
In turn, functions F , G and H will have gradients given by:

∇ σ F = ∂ σ F 1 . . . ∂ σ F N F =    ∂ τ11 F 1 . . . ∂ τ11 F N F . . . . . . ∂ µ32 F 1 . . . ∂ µ32 F N F    (23)
∇ σ G = ∂ σ G 1 . . . ∂ σ G N F =    ∂ τ11 G 1 . . . ∂ τ11 G N F . . . . . . ∂ µ32 G 1 . . . ∂ µ32 G N F    , (24)
and:

∇ q F = ∂ q F 1 . . . ∂ q F N F =    ∂ q1 F 1 . . . ∂ q1 F N F . . . . . . ∂ q Nq F 1 . . . ∂ q Nq F N F    (25)
∇ q H = ∂ q H 1 . . . ∂ q H N F =    ∂ q1 H 1 . . . ∂ q1 H N F . . . . . . ∂ q Nq H 1 . . . ∂ q Nq H N F    . (26)
Second derivatives of G and H will be the third-order tensors:

∇ 2 σσ G = ∇ 2 σσ G 1 . . . ∇ 2 σσ G N F ∇ 2 qq H = ∇ 2 qq H 1 . . . ∇ 2 qq H N F . ( 27
)
Similarly, the plastic multipliers will be assembled into the following vector:

Λ =    λ1 . . . λN F    . (28)
With a notation of this type, expressions for multisurfaces plasticity can be recast in a simplified form. For example, the flow rule (Eq.( 18)) reads:

dε p = ∇ σ G • Λ.

SOLUTION ALGORITHM

Solution algorithm for multisurface plasticity

The elastoplastic problem is incremental, since the elastoplastic solution is in general stresspath dependent. Therefore an iterative solution scheme must be introduced. Herein, an implicit (backward-Euler) Closest-Point-Projection (CPP) solution algorithm is adopted. The CPP algorithm [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF] is probably the most employed method in elastoplasticity. Its use for the multisurface plasticity [START_REF] Simo | Non-smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms[END_REF] is based on a return map, which solves the incremental elastoplastic problem (given by Eqs.( 12),( 13),( 15),( 18),( 19)) and fulfills the restrictions ( 20)- [START_REF] Makrodimopoulos | Lower bound limit analysis of cohesive-frictional materials using second-order cone programming[END_REF]. In this algorithm, the plastic multipliers play an essential role, since they control the loading/unloading conditions and allow to determine the set of activated plastic surfaces by means of a backward procedure. According to [START_REF] Simo | Non-smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms[END_REF], two procedures may be used for the determination of the set of active surfaces (see also [START_REF] Simo | Computational inelasticity[END_REF], Ch.V, p.212). Procedure 1 consists in solving the incremental elastoplastic problem by holding fixed the set of active surfaces during the iterative process and checking Eq.( 20) then. In Procedure 2, the set of active surfaces is updated during the iterative process, by enforcing Eq.( 20) at every iteration. Herein we make use of the Procedure 1. It is in fact more robust than the Procedure 2, even though it is slightly more demanding from a computational point of view [START_REF] Simo | On a Stress Resultant Geometrically Exact Shell Model. Part V: Nonlinear Plasticity: Formulation and Integration Algorithms[END_REF][START_REF] Dujc | Stress resultant plasticity for shells revisited[END_REF].

The algorithm consists of the following steps (Figure 3):

(1) Increments of the displacement du n+1 and Cosserat rotation fields dω c n+1 are given for the generic load step (n + 1). These fields are contained in a single vector ds n+1 .

(2) Total strain increments dε n+1 (relative strains and curvatures) are computed through Eq.( 4)-( 5). ( 3) An elastic solution is predicted as trial [ ] T solution:

                   ε e T n+1 = ε e n + dε n+1 ε p T n+1 = ε p n α T n+1 = α n σ T n+1 = Cε e T n+1 q T n+1 = -Dα T n+1 F T n+1 = F σ T n+1 ; q T n+1 . (29)
In this phase all the N F plastic surfaces F β can be potentially activated (N act = N F ). ( 4) Plastic surfaces are checked (yield criteria). ( 5) If the trial state σ T n+1 ; q T n+1 falls inside the elastic domain E σ , i.e.:

F β ≤ 0, ∀β ∈ [1, . . . , N F ] , (30)
then the trial solution is retained as the solution for the step (elastic increment). ( 6) If at least one of the surfaces is activated, i.e.:

F β > 0, ∀β ∈ [1, . . . , N act ] , (31)
then an elastoplastic increment occurs, and the trial state needs to be corrected. A notable feature of the multisurface plasticity theory is that, when several surfaces are activated, condition [START_REF] Simo | On a Stress Resultant Geometrically Exact Shell Model. Part V: Nonlinear Plasticity: Formulation and Integration Algorithms[END_REF] does not necessarily imply that λβ > 0 for all the activated surfaces (Eq.( 20)-2). Procedure 1 consists then in a) solving the incremental elastoplastic problem, here reported as:

             σ n+1 = C ε n+1 -ε p n+1 q n+1 = -Dα n+1 ε p n+1 = ε p n + ∇ σ G n+1 • Λn+1 α n+1 = α n + ∇ q H n+1 • Λn+1 F n+1 = 0, (32)
Copyright c 0000 John Wiley & Sons, Ltd.

Int. J. Numer. Meth. Engng (0000) Prepared using nmeauth.cls DOI: 10.1002/nme by holding fixed the set of the activated surfaces, and b) checking the solution by testing condition (20)-2. (7) Surfaces for which condition (20)-2 is not satisfied (non-positive plastic multiplier) are dropped from the set of the activated plastic surfaces and a new trial solution is demanded:

(3) ← (7), with N act ≤ N F . (8) If condition (20)-2 is satisfied for all the surfaces, the algorithm has converged to the elastoplastic solution.

Total strain and curvature increment dγ, dκ

Trial solution dτ T , dμ T ; dq T

Displacement and rotation increment du, dω c

Definition of the activated surfaces F β (τ,μ;q) Check plastic surfaces F β (τ,μ;q) ≤ 0 for all activated surfaces F β (τ,μ;q) > 0 for at least one surface

Elastic increment (solution)

Elastoplastic problem Elastoplastic increment (solution) Check consistency λ β > 0 λ β ≤ 0 for at least one surface

λ β > 0 for all activated surfaces (1) (2)
(3)

(4) (5) (6) (7) (8)

Procedure 1 of Simo

Update for the next increment

Implementation of the Newton-Raphson method in the CPP algorithm

The elastoplastic problem ( 32) is solved at every increment (n) by employing a classical (or full) Newton-Raphson iterative method. Generally speaking, the method consists in, given a set of equations of the type r (ζ) -p = 0, solving the equations iteratively, by replacing them at every increment with the linear approximation:

In the above, r (k) = r(ζ (k) ) is the residual vector and ∂ ζ r (k) is the Jacobian matrix of r (k) . Both quantities are known at the increment (k) and, since p is constant, the solution for the increment (k + 1) is simply given by:

ζ (k+1) = ζ (k) -∂ ζ r (k) -1 r (k) -p . ( 34
)
Solution of the starting equations is then found when the norm of the updated residual becomes sufficiently small: r (k+1) -p < T OL N R . In general, the method converges after a number of iterations which depends on the type of equations being considered.

Referring to the CPP algorithm presented in Section 4.1, the residual vector for the incremental elastoplastic problem [START_REF] Stefanou | Three-dimensional Cosserat homogenization of masonry structures: Elasticity[END_REF] writes:

r (k) =    C -1 σ (k) n+1 + ∇ σ G (k) n+1 • Λ(k) n+1 F (k) n+1 D -1 q (k) n+1 + ∇ q H (k) n+1 • Λ(k) n+1    , (35)
where

G (k) n+1 = G σ (k) n+1 ; q (k)
n+1 , while the constant vector p is given by:

p =   C -1 σ T n+1 0 D -1 q T n+1   . ( 36
)
Notice that the equation r (k) -p = 0 corresponds exactly with Eqs. [START_REF] Stefanou | Three-dimensional Cosserat homogenization of masonry structures: Elasticity[END_REF]. In particular, the first and the third term of Eq.( 35)-( 36) are obtained by replacing the set of generalised deformations (ε n+1 ; α n+1 ) in the Eqs. [START_REF] Stefanou | Three-dimensional Cosserat homogenization of masonry structures: Elasticity[END_REF] with the trial state of the generalised stresses σ T n+1 ; q T n+1 (Eq.( 29)). The Jacobian matrix for the incremental elastoplastic problem then writes:

∂ ζ r (k) =     C -1 + ∇ 2 σσ G (k) n+1 • Λ(k) n+1 ∇ σ G (k) n+1 ∇ 2 σq G (k) n+1 • Λ(k) n+1 ∇ σ F (k)t n+1 0 ∇ q F (k)t n+1 ∇ 2 qσ H (k) n+1 • Λ(k) n+1 ∇ q H (k) n+1 D -1 + ∇ 2 qq H (k) n+1 • Λ(k) n+1     . ( 37
)
It is important to note that in case of perfect plasticity one can to consider the submatrix formed by taking the block of the first two entries of the matrix [START_REF] Vardoulakis | Bifurcation analysis in geomechanics[END_REF]. Similarly, the vectors r (k) and p are reduced to their first two terms. The solution for the (k)-th iteration is given by:

ζ (k+1) =    σ (k+1) n+1 Λ(k+1) n+1 q (k+1) n+1    . (38)
The Newton-Raphson method is the core of the CPP algorithm: it assures the existence and the accuracy of the elastoplastic solution at each increment (n). The Newton's scheme is initialized at k = 0 by injecting the trial state (29) into the linearized system [START_REF] Stefanou | Three dimensional homogenization of masonry structures with building blocks of finite strength: A closed form strength domain[END_REF], namely:

ζ (0) =   σ T n+1 0 q T n+1   , (39)
and:

r (0) -p =   0 F (0) n+1 0   , (40)
with

F (0)
n+1 > 0 the vector of the activated surfaces. After this first iteration, the iterative process starts. At each iteration (k), the solution vector for the (k + 1)-th iteration is computed according to Eq. [START_REF] Steinmann | Theory and numerics of ductile micropolar elastoplastic damage[END_REF]. The residual vector [START_REF] Steinmann | Localization within the Framework of Micropolar Elasto-Plasticity[END_REF] is next updated, and its norm is finally checked. Usual conditions for assessing method's convergence towards the elastoplastic solution are [START_REF] Simo | Non-smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms[END_REF]:

F (k+1) n+1 < T OL 1 r (k+1) -p < T OL 2 .
(41)

At that point, the solution of the elastoplastic problem (32) at the increment (n) is presumably found. According to the Procedure 1 (Figure 3), the CPP algorithm requires the plastic multipliers to be all positive. If this condition is respected, the elastoplastic solution given by the Newton's method is retained and the iterative scheme is updated for the (n + 1)-th successive increment. If not, the procedure is reinitialized within the same increment, but with a new, updated trial state (Figure 3). It is to note that, generally, the convergence of this method to the solution is not always assured but if the first-iteration (or trial) vector ζ (0) n+1 does not differ in norm from the expected solution. This could suggests the user to subdivide the analysis in a number of sufficiently small steps. 4.2.1. Limitations of the Newton-Raphson method in the CPP algorithm. When used within the CPP algorithm, the Newton-Rapshon method has some limitations. The limitations derive from the definition of the elastic domain (Eq.( 16)). According to [START_REF] Simo | Computational inelasticity[END_REF], this domain has finite dimension:

dim {E σ } = N σ + N q < ∞. ( 42
)
It follows that each point of the generalised stress space (σ; q) belonging to the boundary of the elastic domain ∂E σ can be intercepted by at most (N σ + N q ) independent plastic surfaces, i.e. surfaces leading to non-redundant constraints [START_REF] Simo | Computational inelasticity[END_REF]. In fact, if one of these surfaces was not independent with the others, the normal vectors ∂ σ F β and ∂ q F β would become linearly dependent. The consequence thereof is that matrices ∇ σ F and ∇ q F and the Jacobian matrix (37) would become singular. In this case, the use of standard matrix inversion techniques in Eq.( 34) would be ineffective and convergence of the Newton's iterative process would not be achieved. Figure 4 shows two distinct configurations leading to singular Jacobian matrix. For illustration purposes, perfect plasticity is chosen and only two components of the stress tensor are kept, i.e. (σ 11 , σ 22 ). In this case, N σ = 2, N q = 0. In the first example (Figure 4-left), 3 surfaces enclose the elastic domain. Therefore N F > (N σ + N q ). Consequently, the projection of the trial stress dσ T at the point of intersection would cause Jacobian matrix singularity at the first Newton's iteration. In the second example (Figure 4-right), only 2 surfaces are defined (N F = 2), but the normals to the surfaces are pointwise coincident. In this case the criteria are clearly redundant. The Jacobian matrix would be then rank deficient, and the Newton's process would newly stop.

σ 11 σ 22 F 1 F 2 F 3 dσ T Redundant plastic surfaces F 1 = F 2 dσ T σ 11 σ 22
Coincident plastic surfaces Both examples are admissible from a physical point of view, and the CPP algorithm should not fail in cases like these. In order to overcome this numerical inefficiency intrinsic to the Newton-Raphson method, and to enable the projection algorithm to fully attain all the features of the multisurface plasticity theory, the Singular Value Decomposition (SVD) technique is adopted in the present formulation for the inversion of the Jacobian matrix in Eq.( 34). This technique is general and robust. In this context, it allows the definition of a theoretically unmlimited number of yield criteria, plastic potentials and hardening/softening laws of every form. The next Section briefly outlines the SVD technique.

Implementation of the Singular Value

Decomposition technique in the Netwton-Raphson method for the CCP algorithm. The Singular Value Decomposition (or SVD) is a technique allowing to compute the pseudo-inverse of a rectangular matrix. Here we use the algorithm based on this technique, presented in [START_REF] Press | Numerical recipes in Fortran 77: the art of scientific computing[END_REF].

The method consists in factorizing a generic rectangular [M × N ] matrix A as the inner product between three matrices, namely U , V and Λ, as follows:

A = U ΛV t . ( 43
)
U and V are two orthogonal [M × M ] and [N × N ] square matrices, respectively. The columns of U are computed as the eigenvectors of the inner product AA t , whereas the columns of V are computed as the eigenvectors of the product A t A. Λ is a [M × N ] diagonal, positive definite matrix containing a i , the non-negative singular values of A (the singular values are the square root of the eigenvalues). Under this form, the computation of the inverse of A reduces to:

A -1 = V Λ -1 U t . ( 44
)
When A is non singular, i.e. its singular values are all positive, the SVD technique allows to compute the inverse matrix A -1 as other methods do (notice however that A is rectangular). When A is singular (or badly conditioned), some of its singular values a 0 i are null (or very close to zero). This would make the inversion of Λ (and A) impossible. In that case, the SVD technique consists in constructing the inverse matrix Λ-1 , by imposing the following condition:

if a 0 i < T OL SV D , then 1/a 0 i = 0. ( 45
)
The resultant matrix A + , indicated as the pseudo-inverse of A, represents the the closest matrix approximation of A -1 :

A + = V Λ-1 U t . ( 46
)
The SVD technique has been implemented within the Newton's iterative process for the inversion of the Jacobian matrix (Eq.( 34)). Cases as those depicted in Figure 4 can be properly treated through this method. When the Jacobian matrix is singular, the SVD gives its pseudo-inverse. When no singularity occurs, the SVD technique always provides an accurate alternative to other classical inversion methods. The effect of the application of the SVD technique on elastoplastic problems is shown in Section 6, by means of specific benchmarks.

FINITE ELEMENT FORMULATION

Element implementation

The theory of multisurface plasticity for Cosserat materials is herein incorporated into a displacement and rotation-based finite element formulation. For further details on the formulation of Cosserat finite elements of this type, the reader is referred to [START_REF] Providas | Finite element method in plane Cosserat elasticity[END_REF][START_REF] Zervos | Finite elements for elasticity with microstructure and gradient elasticity[END_REF][START_REF] Godio | Dynamic finite element formulation for Cosserat elastic plates[END_REF]. The account for plastic strains is made through the adoption of the implicit (backward-Euler) Closest-Point-Projection algorithm and the related solution strategies presented in Section 4.

Dealing with irreducible finite elements, approximations of the kinematic fields are first needed. In the case of Cosserat element, displacement and rotational fields are interpolated with respect to nodal translations U and Cosserat nodal rotations Ω c . This can be made, for instance, by means of specific shape functions N and Φ, associated respectively to the displacements and the rotations [START_REF] Godio | Dynamic finite element formulation for Cosserat elastic plates[END_REF]. However, in order to keep the formulation simple and general, and being in analogy to other more conventional theories, here we express the element interpolation as:

ds = N dS, ( 47
)
where s contains both the displacement u and the rotational fields ω c , and S is the vector of the nodal degrees-of-freedom. Its increment writes:

dS = dU dΩ c t . (48)
Similarly, the increment of the total strain vector ε (Eq.( 11)) of the element is:

dε = BdS, (49)
with B the element compliant matrix. With this notation, the principle of virtual work (Eq.( 6)) on the element of area A e holds at the increment (n + 1):

[R (S n+1 ) -P n+1 ] t δS = 0, (50)
where:

R (S n+1 ) = A e B t σ (S n+1 ) dA, P n+1 = A e N t f n+1 dA (51)
are respectively the resultant vectors of the internal and the external forces f = b l t acting on the element. The dynamic part of Eq.( 6) and the work done by the tractions and moments on the boundary are omitted. Solution to Eq.( 50) can be obtained for every δS and at every increment (n) through an iterative solution procedure. This procedure is controlled at the global level, i.e. for the whole element discretisation, and encloses that adopted for the material (Section 4.2), which is now solved at the local level, i.e. at the quadrature points of the element. The majority of the commercial codes implementing the finite element method employ for this purpose a classical Newton-Rapshon scheme. Therefore, from the knowledge of the iteration (h), by using Eq.( 34) one obtains the solution vector S (h+1) n+1 at the iteration (h + 1):

S (h+1) n+1 = S (h) n+1 -∂ S R (h) n+1 -1 R (h) n+1 -P n+1 . ( 52
)
Once the solution vector S

n+1 is known, the residual vector

R (h)
n+1 is also updated to (h + 1)-th iteration, and (after the assembly operation) convergence is checked at the global level. In order to avoid confusion, in Eq.(52) the index (h) is used to distinguish the global Newton's iteration from the local one, denoted with (k) (Section 4.2), made at the quadrature point.

By using the chain rule along with Eqs.( 14),( 49) and (51), the Jacobian matrix in Eq.( 52) takes the following form:

∂ S R (h) n+1 = ∂ σn+1 R (h) n+1 • ∂ εn+1 σ (h) n+1 • ∂ Sn+1 ε (h) n+1 = A e B t C (h) n+1 B dA = K (h) n+1 , (53)
where

C (h)
n+1 is the matrix of the elastoplastic tangent moduli (Eq.( 14)), computed at the time step (n + 1) and updated at the h-th global iteration. K (h) n+1 is the corresponding element's elastoplastic tangent stiffness matrix . Its expression proves crucial in preserving the rate of quadratic convergence that distinguishes the full Newton-Raphson method from its modified version [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF][START_REF] Simo | Computational inelasticity[END_REF]. The explicit expression of the elastoplastic element stiffness matrix is derived in the next Section, in the framework of the Cosserat multisurface plasticity theory.

Derivation of the elastoplastic tangent stiffness matrix for multisurface plasticity

The peculiar aspect of finite elements based on Cosserat continuum is that both the terms associated with the nodal translations and the nodal rotations confer rigidity to the element. However, in order to avoid further complexity in notation, here the derivation of the stiffness matrix K (h) n+1 is carried out in general form. This can make the derivation similar to that made in other elastoplastic theories, as in the classical continuum theory [START_REF] Simo | Computational inelasticity[END_REF] or in the generalised theories of beams and plates [START_REF] Zienkiewicz | The finite element method for solid and structural mechanics[END_REF]. However, one must have in mind the complex structure of the Cosserat continuum in the context of the multisurface plasticity theory. For instance, due to plasticity, coupling between the terms related to the stresses and the couple-stresses may occur, even when dealing with centro-symmetric materials.

For completeness, we derive below two versions of the tangent moduli

C (h)
n+1 entering Eq.( 53): the continuum and the algorithmic. The continuum version refers directly to the theory of multisurface plasticity developped in Section 3. The algorithmic version refers specifically to the iterative solution scheme presented in Section 4. It is important to note that only the use of the algorithmic moduli in the expression of the element stiffness matrix (53) preserves, whatever the increment size, the quadratic convergence of the Newton's method. The continuum moduli reaches this rate only in the limit of very small step increments [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF][START_REF] Simo | Computational inelasticity[END_REF].

Continuum tangent modulus.

Eqs.( 12),( 13),( 15),( 18),( 19) and ( 21) are recast in the following system of matrix equations:

             dσ = C [dε -dε p ] dq = -Ddα dε p = ∇ σ G • Λ dα = ∇ q H • Λ ∇ σ F t • dσ + ∇ q F t • dq = 0. (54)
Combining Eqs.(54)-1 and (54)-3 one obtains:

dσ = C dε -∇ σ G • Λ . ( 55
)
The above equation, substituted into Eq.(54)-5 together with Eq.(54)-2,4, yields:

Λ = M -1 • ∇ σ F t • C • dε (56)
with:

M = ∇ σ F t • C • ∇ σ G + ∇ q F t • D • ∇ q H (57)
The substitution of Eq.(56) into Eq.(55) leads to the expression for the matrix of the elastoplastic tangent moduli C ep :

dσ = C ep dε C ep = C -C • ∇ σ G • M -1 • ∇ σ F t • C . ( 58
)
It is important to note that the M matrix and the matrix of the elastoplastic tangent modulus C ep are not necessarily symmetric. Both matrices would become symmetric only in the case of associative plasticity, namely:

∇ σ G ≡ ∇ σ F , ∇ q H ≡ ∇ q F . ( 59
)
5.2.2. Algorithmic tangent modulus. The algorithmic, or consistent [START_REF] Simo | Computational inelasticity[END_REF][START_REF] Zienkiewicz | The finite element method for solid and structural mechanics[END_REF], version of the elastoplastic tangent matrix, namely

C (h)
n+1 , is retrieved in similar way with its continuum version C ep (Eq.( 58)). For clarity we omit the index (h). One can rewrite system [START_REF] Stefanou | Three-dimensional Cosserat homogenization of masonry structures: Elasticity[END_REF] in the following algorithmic form:

               dσ n+1 = C [dε n+1 -dε p n+1 ] dq n+1 = -Ddα n+1 dε p n+1 = ∇ 2 σσ G n+1 • Λn+1 • dσ n+1 + ∇ σ G n+1 • d Λn+1 dα n+1 = ∇ 2 qq H n+1 • Λn+1 • dq n+1 + ∇ q H n+1 • d Λn+1 ∇ σ F t n+1 • dσ n+1 + ∇ q F t n+1 • dq n+1 = 0 , (60)
where the increment d Λn+1 is introduced. Eqs.(60)-1,3 are then combined, in order to give:

dσ n+1 = Θ n+1 dε n+1 -∇ σ G n+1 • d Λn+1 , (61)
where Θ n+1 is the matrix of the elastic algorithmic moduli [START_REF] Simo | Computational inelasticity[END_REF]. Its definition is (cf. Eq.( 37)):

Θ n+1 = C -1 + ∇ 2 σσ G n+1 • Λn+1 -1 . ( 62
)
In addition, Eqs.(60)-2,4 give:

dq n+1 = Ψ n+1 -∇ q H n+1 • d Λn+1 , (63)
where Ψ n+1 writes:

Ψ n+1 = D -1 + ∇ 2 qq H n+1 • Λn+1 -1 . (64)
Due to the similarity with Θ n+1 , matrix Ψ n+1 could be defined as the matrix of the hardening algorithmic moduli. Substituted into Eq.(60)-5, Eqs.( 61) and (63) hold:

d Λn+1 = M -1 n+1 • ∇ σ F t n+1 • Θ n+1 • dε n+1 , (65)
and consequently Eq.(61) becomes:

dσ n+1 = C n+1 dε n+1 C n+1 = Θ n+1 -Θ n+1 • ∇ σ G n+1 • M -1 n+1 • ∇ σ F t n+1 • Θ n+1 . (66)
In Eq.( 58), matrix M n+1 is defined as in Eq.( 57), with the exception that matrices C and D must be replaced here with their algorithmic versions Θ n+1 and Ψ n+1 . Notice, in conclusion, how the matrix of the algorithmic moduli (66) reduces to that of the continuum moduli (58), as the step increment is small. This shows that problem (60) is consistent with problem (54).

ELEMENT VALIDATION

In Section 4 and 5, the theory, the solution algorithm and the element implementation of multisurface plasticity for Cosserat materials have been presented in general manner. In this Section attention is focused on a specific finite element, formulated for micropolar plate models. The element is briefly introduced and its validation in the framework of multisurface plasticity is then carried out in two steps. Firstly, the accuracy of the CPP algorithm is tested locally, i.e. at the integration point level, through simple element tests. Secondly, the efficiency of the algorithm is assessed at the global (or element) level, based on more complex discretisations related to structural problems.

The micropolar plate element COSS8R

The element called COSS8R [START_REF] Godio | Dynamic finite element formulation for Cosserat elastic plates[END_REF] is a quadratic micropolar plate element (Figure 5). It possesses 8 nodes, all situated on the element sides. Each node is equipped with 6 degrees-of-freedom, i.e. 3 translations and 3 Cosserat rotations. In this element, the shape functions N and Φ, used respectively for the interpolation of the displacement and the rotational fields, are the same. The integration of the element stiffness matrix is made on 4 quadrature points (reduced integration), whereas the integration of the element mass matrix requires 9 quadrature points (full integration). This avoids the shear locking phenomenon (in statics) and the occurrence of massless degrees-offreedom (in dynamics). For further details, the reader is referred to [START_REF] Godio | Cosserat Elastoplastic Finite Elements for Masonry Structures[END_REF]. It is important to note that, in the extension to plasticity, the presence of multiple quadrature points requires the independent solution of the elastoplastic problem at each of these points.

Ω 2 Ω 1 Ω 3 U 1 U 2 U 3
Figure 5. The quadratic micropolar plate element COSS8R [START_REF] Godio | Dynamic finite element formulation for Cosserat elastic plates[END_REF].

membrane behaviour flexural behaviour For the sequel of this Section, it is useful to distinguish the variables describing the membrane response from the variables describing the flexural response of the plate element (Table I). Accordingly, the membrane elastic response is controlled by matrices A and D m , expressed as [START_REF] Papanastasiou | Numerical treatment of progressive localization in relation to borehole stability[END_REF][START_REF] Godio | Dynamic finite element formulation for Cosserat elastic plates[END_REF]:

u m = u 1 u 2 t u f = [u 3 ] ω c m = [ω 3 ] ω c f = ω 1 ω 2 t γ m = γ 11 γ 22 γ 12 γ 21 t γ f = γ 13 γ 31 γ 23 γ 32 t κ m = κ 31 κ 32 t κ f = κ 11 κ 22 κ 12 κ 21 t b m = b 1 b 2 t b f = b 3 l m = l 3 l f = l 1 l 2 t τ m =
τ m = Aγ m , µ m = D m κ m A =    K + G K -G 0 0 K -G K + G 0 0 0 0 G + G c G -G c 0 0 G -G c G + G c    , D m = 2Gl c 0 0 2Gl c ( 67
)
with K the compression modulus, G the shear modulus, G c the Cosserat shear modulus and l c the characteristic length of the micro-structured material. The flexural elastic response is controlled by matrices Q and D f , reading [START_REF] Altenbach | On the linear theory of micropolar plates[END_REF][START_REF] Godio | Dynamic finite element formulation for Cosserat elastic plates[END_REF]:

τ f = Qγ f , µ f = D f κ f Q =    0 0 0 0 0 α 4 0 0 0 0 0 0 0 0 0 α 4    , D f =    β 1 + β 2 + β 3 β 1 0 0 β 1 β 1 + β 2 + β 3 0 0 0 0 β 3 β 2 0 0 β 2 β 3    ( 68
)
with α 4 the material parameter related to the transverse shears, and β 1 , β 2 , β 3 those related to the torsions and the out-of-plane flexions of the plate.

Single element tests for multisurface plane plasticity

The scope of single element tests is twofold. On one hand, the algorithm accuracy is assessed not only pointwise, but also at the element level, i.e. after integration of the element's stiffness matrix. This guarantees the good response of the element and the convergence of the iterative scheme at the global level. On the other hand, with tests of this kind the performance of the element can be shown in special configurations, as for states of non-homogeneous deformation or in presence of multiple plastic surfaces. In this setting, the use of the SVD technique within the Newton-Raphson process is investigated.

Herein we carry out two series of tests based on two basic configurations (Figure 6): one involving uni-axial loading/unloading cycles, and one involving the use of multiple plastic surfaces for bi-axial stress states. All tests are displacement-controlled.

u 1 Uni-axial test u 1 u 2
Bi-axial test Figure 6. Configurations used for the single element tests in multisurface plasticity.

6.2.1. Uni-axial loading/unloanding cycles. The application of the displacement increment ∆u 1 leads to a homogeneous state of tension/compression/tension on the element (Figure 6-left). The uni-axial constitutive law is given by:

τ 11 = Eγ 11 , (69)
with E the elastic modulus. The plastic surface with linear isotropic hardening is given in the form:

|τ 11 | -τ y + q 1 ≤ 0 q 1 = hα 1 (70)
with τ y the uni-axial limit stress. Figure 7 illustrates the stress path followed during the test and the stress-strain relation curves. The evolution of the elastic domain during the elastoplastic loading phases (a-b) and (c-d) is apparent.

Bi-axial load tests for multisurface plasticity.

Five tests, indicated with B1 to B5, are carried out. In these tests, displacement increments are applied in x 1 and x 2 directions (Figure 6-right), leading the element to homogeneous bi-axial stress states. Different sets of plastic surfaces are considered for each test. In particular, we consider five plastic surfaces together with a linear isotropic hardening/softening function:

F 1 = τ 22 /τ y -1 + q 1 F 2 = τ 11 /τ y -1 + q 1 F 3 = (τ 11 + τ 22 )/(2τ y ) -1 + q 1 F 4 = (τ 11 /τ y ) 2 + (τ 22 /τ y ) 2 -1 + q 1 F 5 = (τ 11 /(1.4 τ y )) 2 + (τ 22 /(0.6 τ y )) 2 -1 + q 1 q 1 = hα 1 . (71)
Moreover, we consider a material with isotropic elastic constitutive law, given by Eq.( 67). All tests except Test B5 are carried out by using the same elastic properties (Figure 8), and the same number of equally spaced increments [START_REF] Simo | Consistent tangent operators for rate-independent elastoplasticity[END_REF]. This allows the convergence to be reached in Tests B1, B2 and to provide high level of accuracy in Test B3, B4. For Test B5, the total number of increments is doubled and the hardening modulus is 0.1h.

In Test B1 (Figure 8) surfaces F 1 , F 2 are used. These two intersecting plastic surfaces are activated during the monotonic bi-axial loading. For perfect plasticity (Figure 8-left), the stress path is limited by the boundary of the elastic domain. In case of isotropic hardening (Figure 8-right), both surfaces evolve towards the new stress state of the element. Notice that the symmetry of the problem is reflected in the symmetry of the solution, which is the same in x 1 -and x 2 -directions. Moreover, the same solution is reached at all integration points of the element. For brevity, only τ 11 as function of the total γ 11 , elastic γ e 11 and plastic deformation γ p 11 , and of the hardening parameter α is presented.

In Test B2 (Figure 9) surface F 3 is added to F 1 , F 2 and in Test B3 (Figure 10) surface F 4 and an identical surface denoted F 6 are employed. These two tests follow the examples of Figure 4, since they allow to consider redundant and coincident plastic surfaces respectively. If the SVD technique was not employed for the solution of the incremental elastoplastic problem (Section 4.2.1), the inversion of the Jacobian matrix would not be possible in these cases and the test would fail once the surfaces reached. This would hold for Test B2 in the case of perfect plasticity (N F = N σ ) and for Test B3 in both cases, since the same hardening variable is used for the two surfaces F 4 , F 6 . Nevertheless, these tests are handled by the element, thanks to the proposed solution strategy that involves the SVD technique. On one hand (Figure 9), the element is able to activate all the necessary surfaces. Since surface F 3 is redundant with F 1 and F 2 , the resulting deformations are actually the same with those of Test B1. On the other hand (Figure 10), both surfaces evolve in identical manner.

Test B4 (Figure 11) makes use of F 4 and F 6 . It consists of two successive monotonic load steps. First ∆u 1 > 0 and then ∆u 2 > 0 are applied. As a result of the form of the elastic domain, the second load step results in a loading step in direction x 2 and an unloading step in direction x 1 . This occurs both in case of hardening and of perfect plasticity, where the stress state is constrained to move along the boundary of the elastic domain. Also in this test, the SVD technique enables the element to follow exactly the stress path imposed, F 4 and F 6 being coincident.

In Test B5 (Figure 12) surfaces F 4 and F 5 are used. The test consists in five successive load steps, namely: a) ∆u 1 > 0, ∆u 2 = 0; b) ∆u 1 = 0, ∆u 2 > 0; c) ∆u 1 < 0, ∆u 2 = 0; d) ∆u 1 = 0, ∆u 2 < 0; e) ∆u 1 > 0, ∆u 2 = 0. Following these steps, the element carries out a complete loading/unloading cycle in x 1 and x 2 directions simultaneously. Accordingly, the stress path activates different plastic surfaces at each load step. Notice how the activation of the surfaces occurs precisely at their points of intersection. Moreover, in case of hardening plasticity the accumulated plastic strains induces the general expansion of the elastic domain. Hardening plasticity Figure 8. Element response to bi-axial test B1. Top: stress path in the (τ 11 -τ 22 )-plane. Bottom: stressstrain curve. In-plane elastic properties: K = 5MPa, G = 5MPa, Gc = 0MPa, lc = 0. Hardening modulus h = 10MPa, yield stress τy = 5MPa.

Structural benchmarks for micropolar plates

Structural benchmarks are made to test (and suggest) the use of the algorithm of multisurface plasticity for Cosserat materials and the developed plate element in cases of practical interest. Herein we carry out three benchmarks. In all the examples we consider 3D structural problems. The first benchmark is concerned with the development of in-plane strain localisation in hollow cylinders with micropolar properties. With the second benchmark we model the behaviour of micropolar square plates undergoing out-of-plane macroscopic plastic curvatures. The third benchmark aims at studying the load-bearing capacity of a shallow foundation resting upon a micropolar soil under plane strain conditions.

6.3.1. Strain localisation in hollow cylinders. We consider a hollow cylinder of height H = 60mm, radius R = 15mm and constant thickness t = 1mm. The cylinder is subjected to torsion (Figure 13).

In particular, the top and the bottom sides of the cylinder have all the degrees-of-freedom prevented except the axial rotation, denoted with Ω 3 , which is controlled. The lateral surface of the cylinder is free of stresses and couple-stresses. The cylinder is made of a micro-structured material. Its macroscopic behaviour is described by a Cosserat continuum with homogeneous isotropic centro-symmetric elastic properties (Eq.( 67)-( 68)). The in-plane response of the material follows here the tangential (θ) and longitudinal (z) directions of the cylinder. The out-of-plane direction is on the radial (r) direction (Figure 13). The elastic parameters of the material are: K = 4000MPa, G = 4000MPa, G c = 2000MPa, and α 4 = 1500MPa, β 1 = 2000MPa × mm, β 2 = 0MPa × mm, β 3 = 4000MPa × mm. The characteristic (or inner) length of the micro-structured material is l c = 2.5mm (Eq.( 67)).

The elastoplastic response of the material is described by means of the J2 (or Von-Mises) plasticity model formulated for Cosserat continua. The yield criterion for the generalised J2hardening-plasticity relies on the following single plastic surface [START_REF] De Borst | Simulation of strain localisation: a reappraisal of the Cosserat continuum[END_REF]:

F = 3J 2 -σ (γ p 2 ) ( 72
)
where σ represents the equivalent yield stress, J 2 is the second invariant of the deviatoric stress tensor and γ p 2 is the second invariant of the deviatoric plastic strain tensor. For the Cosserat continuum, J 2 and γ p 2 write respectively:

J 2 = 1 2 τ t m µ t m P τ m µ m , (73)
and:

γ p 2 = 2 3 γ t m κ t m P γ m κ m , (74)
with matrices P and P defined as:

P =        4 3 (a 1 + a 2 ) -2 3 (a 1 + a 2 ) 0 0 0 0 -2 3 (a 1 + a 2 ) 4 3 (a 1 + a 2 ) 0 0 0 0 0 0 2a 1 2a 2 0 0 0 0 2a 2 2a 1 0 0 0 0 0 0 2a 3 /l 2 c 0 0 0 0 0 0 2a 3 /l 2 c        P =        4 3 (b 1 + b 2 ) -2 3 (b 1 + b 2 ) 0 0 0 0 -2 3 (b 1 + b 2 ) 4 3 (b 1 + b 2 ) 0 0 0 0 0 0 2b 1 2b 2 0 0 0 0 2b 2 2b 1 0 0 0 0 0 0 2l 2 c b 3 0 0 0 0 0 0 2l 2 c b 3        . ( 75
)
Scalars a 1 , a 2 , a 3 and b 1 , b 2 , b 3 are the plastic parameters of the material. Following [START_REF] De Borst | Simulation of strain localisation: a reappraisal of the Cosserat continuum[END_REF], we use

a 1 = 1/4, a 2 = 1/4, a 3 = 1/2 and b 1 = 1/3, b 2 = 1/3, b 3 = 2/3.
Herein, strain-softening plasticity is simulated through a linear softening rule, with hardening modulus h = -400MPa and yield stress σ = 100MPa.

Simulations are carried out with three different finite element discretisations. The COSS8R being a flat element, it is suitable to choose a larger number along the tangential direction of the cylinder. Therefore, assemblies of 16 × 8, 24 × 16 and 32 × 24 elements are used. Analyses are run by updating the reference system at every load step. In this way large axial rotations are avoided. Moreover, an automatic step increment control technique is used throughout the analyses.

The resulting macroscopic behaviour of the cylinders is shown in Figure 14-left, in terms of normalized reaction moment M z and differential axial rotation ∆Ω z = Ω t z -Ω b z . The onset of plastic deformations (on the first branch) rapidly culminates with a softening branch, which is accompanied by strain localisation (Figure 14-right). The deformation pattern of the cylinders beyond the moment peak value is shown in Figure 15 for the different discretisations adopted. This result is expected since, due to the axial symmetry of the problem, the deformation is the same with that of other 2D structural problems [START_REF] De Borst | Simulation of strain localisation: a reappraisal of the Cosserat continuum[END_REF][START_REF] Sharbati | Computational aspects of the Cosserat finite element analysis of localization phenomena[END_REF]. Accordingly, the linear elastic branch is well captured even by the coarsest discretisation, and the softening branch is reached by adopting slightly finer discretisations. Strain localisation occurs exactly as in the shear layers studied by [START_REF] De Borst | Simulation of strain localisation: a reappraisal of the Cosserat continuum[END_REF]. The strain localisation shows then an apparent mesh-independency (Figure 14-right), i.e. the width of the localisation region does not depend on the number of elements falling inside this region. This is a fundamental feature of finite elements based on the Cosserat continuum [START_REF] De Borst | Fundamental issues in finite element analyses of localization of deformation[END_REF][START_REF] Papanastasiou | Numerical treatment of progressive localization in relation to borehole stability[END_REF][START_REF] Neff | a Numerical Solution Method for an Infinitesimal Elasto-Plastic Cosserat Model[END_REF].

Torsion of square plates.

As second example we consider a square plate, of span L = 10mm and thickness t = 1mm (Figure 16). The plate has two opposite free edges and two opposite clamped edges, with all the degrees-of-freedom prevented. A rotation Ω 2 is imposed at one of the clamped edges. The rotation produces a macroscopic deformation that involves the out-of-plane behaviour of the plate. We are here interested in the assessment of the ultimate (or limit) resulting moment M 2 supported by the plate, under multiple intersecting plastic surfaces.

The plate is made of a micro-structured material which has the same elastic micropolar properties with that composing the hollow cylinders considered in Section 6.3.1. However, criteria for the outof-plane behaviour of a Cosserat plate are absent in the literature. One criterion which is often used in practice is the Nielsen's criterion, formulated for reinforced concrete plates [START_REF] Makrodimopoulos | Lower bound limit analysis of cohesive-frictional materials using second-order cone programming[END_REF][START_REF] Krabbenhø Ft | Formulation and solution of some plasticity problems as conic programs[END_REF]. This criterion is usually expressed for conventional (or engineering) plate theories. Herein, it is formulated in order to take into account the out-of-plane statics of the micropolar plate model (Figure 1). It consists of six distinct and intersecting plastic surfaces:

F 1 = (µ 11 -µ 22 ) 2 /4 -(M + 21 -µ 21 )(M + 12 -µ 12 ) F 2 = (µ 11 -µ 22 ) 2 /4 -(M - 21 + µ 21 )(M - 12 + µ 12 ) F 3 = µ 21 -M + 21 F 4 = µ 12 -M + 12 F 5 = -µ 21 -M - 21 F 6 = -µ 12 -M - 12 . (76)
Parameters M + 21 , M - 21 , M + 12 , M - 12 are the (non-negative) out-of-plane yield flexions µ 12 and µ 21 of the reinforced concrete plate. M + and M -represent the ultimate moment leading respectively the upper and the lower side fiber of the plate to failure. For simplicity, these parameters are all considered equal to m = 5MPa × mm. For illustration purposes, simulations are carried out in the frame of perfect and associative plasticity. It is worth mentioning that the conventional torsional moment of the engineering theory of plates, M t , is here expressed in terms of the Cosserat torsions µ 11 and µ 22 by imposing: M t = (µ 11 -µ 22 ) /2. In this way the Nielsen's criterion may be extended to Cosserat plate models.

The plate is modeled by means of three different finite element discretisations. Figure 17 shows that, even with the coarsest 2 × 2 discretisation, the macroscopic response of the plate is captured accurately by the COSS8R element. At the beginning the plate behaves elastically, but as far as some of the plastic surfaces are activated (Eq.( 76)), plastic curvatures are developed until the formation of a macroscopic failure mechanism. This corresponds to the ultimate resulting moment M 2 . Figure 18 displays the number of activated surfaces once the limit moment is reached. It is found that, for the problem considered, only the first two plastic surfaces F 1 , F 2 are activated throughout the analyses. It is apparent how the algorithm of multisurface plasticity is able to select different sets of plastic surfaces at each integration point of the finite element. This confers high flexibility to the element and may increase the accuracy in problems characterised by states of non-homogeneous deformations.

6.3.3. Shallow foundation on micropolar soil. In the third benchmark we consider a shallow foundation supported by a homogeneous layer of soil (Figure 19). The foundation has a width B = 33.34m and is subjected to a distributed pressure T 2 . The layer of soil has dimensions L × H = 200m × 50m. Classical smooth conditions are applied to the boundaries of the layer. The geometry of the problem allows then to consider only half of the model. This is made by preventing the horizontal translations U 2 and the in-plane rotations Ω 3 along the symmetry plane. It is worth noting that plane strain conditions, for which the out-of-plane component τ 33 is non-null [START_REF] De Borst | Fundamental issues in finite element analyses of localization of deformation[END_REF], are formulated for this benchmark.

The foundation is considered as rigid and the interface with the soil perfectly cohesive. The soil is an isotropic material with micropolar elastic properties given by K = 34000kPa, G = 18000kPa and G c = 0.7G (Eq.( 67)). The elastoplastic response of the soil is described by means of the following set of plastic surfaces (Figure 20):

F 1 = 3J 2 + Ap -σ F 2 = 3J 2 -Cp -Cp c . ( 77
)
The first surface is the generalisation of the Drucker-Prager plasticity model for Cosserat continua [START_REF] De Borst | Fundamental issues in finite element analyses of localization of deformation[END_REF][START_REF] Papanastasiou | Numerical treatment of progressive localization in relation to borehole stability[END_REF][START_REF] Li | A consistent return mapping algorithm for pressure-dependent elastoplastic Cosserat continua and modelling of strain localisation[END_REF], for which p is the mean pressure and A = 0.39. The second surface introduces a generalised compression cap to the model, with p c = 2496kPa = 3σ the yield stress in compression and C = 0.81. No hardening is considered. The aim of the benchmark is twofold. On one hand, the variation of the load-bearing capacity of the foundation is assessed for soils with different internal lengths l c (scale effect). The effect of different yield criteria is also studied (Eq.( 77)). On the other hand, the ultimate load is computed in the limit of small lengths in order to obtain the solution for the Cauchy continuum.

Figure 21 shows the macroscopic response of the foundation in terms of normalized pressure T 2 and vertical displacement δ 2 . The scale effect is apparent on Figure 21-left: a larger internal length yields to a stiffer elastoplastic response and a higher load-bearing capacity. On the contrary, the activation of a compression cap leads to a consistent reduction of the ultimate load. The macroscopic response for the Cauchy continuum (Figure 21-right) is obtained by imposing an internal length which is very small as compared to the size of the foundation. Moreover, the conditions G c = 0 and a 3 = 0 allow respectively to retrieve the symmetry of the elastic response (Eq.( 67)) and the classical J2-plasticity model (Eq.( 72)). For validation purposes, the macroscopic response given by the COSS8R element is compared with that obtained by a finite element formulated for classical continua. To this end, a 6-noded Gauss-type mixed triangular element [START_REF] Krabbenhø Ft | Formulation and solution of some plasticity problems as conic programs[END_REF][START_REF] Krabbenhø Ft | An interior-point algorithm for elastoplasticity[END_REF] contained in the commercial code OptumG2 is used. The solution algorithm used in OptumG2 differs from that presented herein, since based on optimization methods [START_REF] Krabbenhø Ft | An interior-point algorithm for elastoplasticity[END_REF]. The comparison shows how the COSS8R element is able to cover the solution for the Cauchy continuum, with an error less then 1%.

CONCLUSIONS

The development of robust and efficient numerical procedures is a principal task prior to modeling materials with complex micro-structure. In the attempt of describing the macroscopic behaviour of such materials, this aspect is even more pronounced. Particularly important is in fact the determination of the overall response of the material, as function of the properties and the spatial arrangement of its micro-structure. Especially in the inelastic range, this is not a trivial task, since even when the macroscopic properties are established, one needs to resort to numerical tools for the assessment of the material response.

The use of the Cosserat continuum has gained much interest in this field, by virtue of its relatively simple formulation (limited number of variables and clear physical meaning) which yet allows to handle quite complex configurations of materials with pronounced micro-structure (masonry, soils, rocks, grid-works, etc.).

This paper was devoted to the development and validation of a numerical procedure for the analysis of micro-structured materials with macroscopic non-linear behaviour. For this purpose, the time-independent multisurface plasticity theory was extended to the Cosserat continuum. Multisurface plasticity theory allows the description of a large class of mechanisms, such as strain hardening and softening mechanisms, and the definition of both associative and non-associative properties for the material. Moreover, through this theory the elastic domain may be constructed by using a set of multiple plastic, intersecting surfaces of general form. When referred to the macroscopic description of the material, this may prove particularly useful, as each plastic surface may be related to distinct failure mechanisms at the micro-structure level.

The implementation of the multisurface plasticity theory into a Cosserat finite element formulation was presented in general manner. The implementation of the theory was based on a projection algorithm formulated for implicit time-integration schemes. The same algorithm may be encountered in more classical formulations, as the theory of classical (Cauchy) continua and the generalised theories of beams and shells. The paper focused on several numerical aspects concerning the element implementation, in relation with the iterative solution schemes used at the local (integration point) and global (element) level. A full Newton-Raphson method was used at both levels. At the local level, the projection algorithm required the use of a non conventional matrix-inversion technique (SVD). This technique was introduced in order to fully exploit all the salient features of the multisurface plasticity theory and to avoid some important limitations that were explored in details. It allowed the definition of plastic surfaces which are (hypothetically) unlimited in number and of the most general form. At the global level, the full rate of convergence was assured by the use of the elastoplastic stiffness tangent matrix of the Cosserat element. This matrix was herein derived explicitly in algorithmic form and for the very general multisurface and non-associative case.

A specific Cosserat finite element developed for micropolar plate models was finally presented. The element considered was already developed and tested in elasticity and dynamics. This paper was concerned with its extension and validation to the multisurface plasticity. Numerical tests assessed the (high) level of accuracy reached by the projection algorithm in case of multisurface Cosserat plasticity. Several benchmarks showed that the element provides accurate results under non-conventional loading and boundary conditions, involving strain softening, multiple plastic surfaces and strain localisation. Hardening plasticity Right: response for small internal lengths (lc = 0.00025m) and comparison with finite elements based on the Cauchy continuum [START_REF] Krabbenhø Ft | An interior-point algorithm for elastoplasticity[END_REF].

  vector µ: τ = τ 11 τ 22 τ 12 τ 21 τ 13 τ 31 τ 23 τ 32 t µ = µ 11 µ 22 µ 12 µ 21 µ 31 µ 32 t .

Figure 2 .

 2 Figure 2. Schematic representation of the elastic domain of a Cosserat continuum. According to the multisurface plasticity theory, it is a convex region in the generalised stress space (τ , µ; q) defined by multiple yield criteria.

Figure 3 .

 3 Figure 3. Illustration of the Closest-Point-Projection (CPP) algorithm adopted within the multisurface plasticity theory for Cosserat materials.

Figure 4 .

 4 Figure 4. Two possible configurations for which Jacobian matrix is singular for perfect plasticity.

11 e τ 11 vs. γ 11 p τ 11 vs. γ 11 - 1

 1111111 Stress path in the (τ 11 -τ 22 )-plane τ 11 vs. γ

Figure 7 .

 7 Figure 7. Element response to uni-axial cyclic loading. Elastic modulus E = 10MPa, hardening modulus h = 10MPa, yield stress τy = 5MPa.

Figure 9 .

 9 Figure 9. Element response to bi-axial test B2. Top: stress path in the (τ 11 -τ 22 )-plane. Bottom: stress-strain curve.

Figure 11 .

 11 Figure 11. Element response to bi-axial test B4. Top: stress path in the (τ 11 -τ 22 )-plane. Center and bottom: stress-strain curve.

Figure 12 .

 12 Figure 12. Element response to bi-axial test B5. Top: stress path in the (τ 11 -τ 22 )-plane. Center and bottom: stress-strain curve.

2 LFigure 16 . 2 Figure 17 . 2 Figure 18 .

 216217218 Figure 16. Square plate under torsion. Illustration of the boundary conditions and of the expected deformation pattern.

Figure 19 .

 19 Figure 19. Shallow foundation under distributed pressure supported by a soil with micropolar properties.Geometry and boundary conditions.

Figure 20 .

 20 Figure 20. Illustration of the set of plastic surfaces describing the elastoplastic response of the soil.

F 1 Figure 21 .

 121 Figure 21. Left: macroscopic response of the shallow foundation for soils with different internal lengths. Right: response for small internal lengths (lc = 0.00025m) and comparison with finite elements based on the Cauchy continuum[START_REF] Krabbenhø Ft | An interior-point algorithm for elastoplasticity[END_REF].

11 κ 22 κ 12 κ 21 κ 31 κ 32 t

  = γ 11 γ 22 γ 12 γ 21 γ 13 γ 31 γ 23 γ 32

			τ 22				μ 22
		μ 32	12			τ 32	μ 12
	τ 11	μ 31	τ 21	τ 11	μ 11	τ 31	τ 13	μ 21	μ 11
		τ 21	μ 31			μ 21	τ 23	τ 31
		τ 12	μ 32			μ 12	τ 32
			τ 22				μ 22
		in-plane behaviour

out-of-plane behaviour

Figure 1. Stresses and couple-stresses of a Cosserat plate.

The associated deformation measures of the plate are denoted with γ and κ and are:

γ t κ = κ

  τ 11 τ 22 τ 12 τ 21

t τ f = τ 13 τ 31 τ 23 τ 32 t µ m = µ 31 µ 32 t µ f = µ 11 µ 22 µ 12 µ 21 t Table I. Membrane and flexural static and kinematics of a Cosserat plate element.

  Element response to bi-axial test B3. Top: stress path in the (τ 11 -τ 22 )-plane. Bottom: stressstrain curve.

		8 0	2	4	6		8 8		0 8	2	4	6	8 8
		6	stress path				6		6		6
	τ 22 [MPa]	4					4	τ 22 [MPa]	4		4
		2					2		2		2
		0 0	2	4	6 F 4 =F 6		8 0		0 0	2	4	6 F 4 =F 6	8 0
			τ 11 [MPa]							τ 11 [MPa]
			Perfect plasticity					Hardening plasticity
		8 0.0	0.2	0.4	0.6	0.8 8		8 0.0	0.2	0.4	0.6	0.8 8
		6					6		6		6
	τ 11 [MPa]	4			τ 11 vs. γ 11 e	4	τ 11 [MPa]	4		4
		2			p τ 11 vs. γ 11 τ 11 vs. α τ 11 vs. γ 11	2		2		2
		0.0 0	0.2	0.4	0.6	0.8 0		0.0 0	0.2	0.4	0.6	0.8 0
			[%] γ 11 ,α							[%] γ 11 ,α
			Perfect plasticity					Hardening plasticity
	Figure 10.

r (k+1) -p ≈ r (k) + ∂ ζ r (k) ζ (k+1) -ζ (k) -p = 0. (33)Copyright c 0000 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (0000) Prepared using nmeauth.cls DOI: 10.1002/nme