The Isotropic Material Design of In-Plane Loaded Elasto-Plastic Plates

This paper puts forward a new version of the Isotropic Material Design method for the optimum design of structures made of an elasto-plastic material within the Hencky-Nadai-Ilyushin theory. This method provides the optimal layouts of the moduli of isotropy to make the overall compliance minimal. Thus, the bulk and shear moduli are the only design variables, both assumed as non-negative fields. The trace of the Hooke tensor represents the unit cost of the design. The yield condition is assumed to be independent of the design variables, to make the design process as simple as possible. By eliminating the design variables, the optimum design problem is reduced to the pair of the two mutually dual Linear Constrained Problems (LCP). The solution to the LCP stress-based problem directly determines the layout of the optimal moduli. A numerical method has been developed to construct approximate solutions, which paves the way for constructing the final layouts of the elastic moduli. Selected illustrative solutions are reported, corresponding to various data concerning the yield limit and the cost of the design. The yield condition introduced in this paper results in bounding the values of the optimal moduli in the places of possible stress concentration, such as reentrant corners.


Introduction
The problem of designing structures made of a linear elastic material is one of the major topics of Free Material Design (FMD). Within this approach, all the elastic moduli of Hooke's tensor C are design variables. Usually, the aim is to minimize the compliance of the structure, while the unit cost is identified with the trace of Hooke's tensor or the sum of its eigenvalues (see [1,2], to mention the first papers on the topic). The additional assumption of isotropy reduces the number of design variables to two: the bulk modulus k and shear modulus µ (see [3,4], where this method, called there the Isotropic Material Design (IMD), was proposed). In the 3D setting, the eigenvalues of the Hooke tensor are: 3k, 2µ, 2µ, 2µ, 2µ, 2µ, hence tr C = 3k + 10µ. In the 2D setting, the eigenvalues of C are: 2k, 2µ, 2µ, hence tr C = 2k + 4µ. The present paper refers to those papers on FMD in which the Hooke tensor is subject only to the condition of positive semi-definiteness; in the case of the IMD method, this condition reduces to: k ≥ 0, µ ≥ 0. The upper bounds are absent to make the theory as simple as possible. Admitting the vanishing of moduli means working with the broadest possible class of the underlying microstructures. For instance, the hexagonal (or honeycomb in the plane) gridwork is characterized by a very small shear modulus, if ligaments are slender (see [5,6]). On the other hand, spiral microstructures are characterized by very small bulk modulus, which implies the effective Poisson ratio almost attaining its lower 2D limit equal −1 (see [7][8][9][10][11][12]). To encompass such a broad class of composites, it is necessary to admit the largest possible range of the bulk and shear moduli. Due to the simplicity of such modeling, it is possible to perform minimization over the moduli analytically, thus eliminating the design variables in the first step. Eventually, one arrives at two, mutually dual, linear constrained problems (LCP) in the meaning of The sign~means that the tensor is represented by the given matrix in the fixed Cartesian coordinate system. Both stress and strain tensors are symmetric. Let I represent the identity matrix, or I = diag [1,1,1]. The scalar product of two vectors u, v is defined by u · v = u x v x + u y v y + u z v z . The set of 2nd rank symmetric tensors will be denoted by E 2 s . The scalar product of σ, ε ∈ E 2 s is defined by: σ · ε = σ x ε x + 2σ xy ε xy + σ y ε y + 2σ xz ε xz + 2σ yz ε yz + σ z ε z .
The Euclidean norms of the vectors and tensors in E 2 s are denoted by u = √ u · u, σ = √ σ · σ. The trace of the tensor σ ∈ E 2 s is given by tr σ = σ x + σ y + σ z . The deviator of σ ∈ E 2 s is defined by: The Euclidean norm of the deviator reads: In the 2D case (d = 2), the tensors of stress and strain are represented by the matrices: The identity matrix is defined by I = diag [1,1]. The trace of the tensor σ ∈ E 2 s is given by tr σ = σ x + σ y . The deviator of σ ∈ E 2 s is defined by: dev σ = σ − 1 2 (tr σ) I (6) or dev σ ∼ 1 2 σ x − σ y σ xy σ yx The Euclidean norm of the deviator reads: For both cases of d = 2 or d = 3, the scalar product of the two tensors from E 2 s can be rewritten as: σ · ε = Tr σ · Tr ε + dev σ · dev ε (9) where: which is a modified trace of a tensor. According to the linear theory of the continuum media, the strain tensor is the symmetric part of the gradient of the displacement vector. In the case of d = 2, we define the operation: which determines the virtual strains corresponding to the virtual displacement field v. For a given function f (·) of argument σ ∈ E 2 s , one can define its polar by:

On the Hencky-Nadai-Ilyushin Theory of an Elasto-Plastic Body
Within the theory by Hencky-Nadai-Ilyushin (also called Hencky's theory, see [15]), the stress state σ is locally constrained by the plasticity condition: where Ω is the domain occupied by the body. The function F is assumed to be convex and continuous with respect to all stress components. Let us recall the HMH plasticity condition for isotropic metals proposed by Huber, Mises and Hencky (see [21]): where dev σ is given by Equation (4) and refers to the 3D setting, and σ 3D 0 is the plastic limit corresponding to the tensile test. Thus, the hydrostatic state of stress σ = pI cannot cause plastic yielding, irrespective of the sign of the pressure p.
The present paper deals with the optimum design of in-plane loaded, transversely homogeneous thin plates of thickness b; Ω will be its middle plane parameterized by the (x,y) system. In such a plate, the stress components σ z , σ xz , σ yz are negligible in comparison to other stresses. Substituting: σ z = 0, σ xz = 0, σ yz = 0 into Equation (14) leads to the HMH condition for the plane-stress problem: σ e f f ≤ σ 0 , σ e f f = σ 2 x − σ x σ y + σ 2 y + 3σ 2 xy .
Here, the stress resultants are involved, still denoted by σ x , σ xy , σ y of units N/m and σ 0 = bσ 3D 0 . It is worth noting that substitution σ z = 0, σ xz = 0, σ yz = 0 into Equation (4) does not lead to Equation (8). Indeed, Equation (15) now involves both stress invariants within the 2D setting, since: with Tr σ= (tr σ)/ √ 2 and dev σ defined by Equation (8). This shows that the function σ e f f (σ) is isotropic. The function polar to γ(σ) has the form: Thus, we see that its construction can be performed by inverting the coefficients in Equation (16). The simplicity of this construction follows from the orthogonality of the tensors Tr σ · I and dev σ and from Equation (9). Thus, the yield condition has the form of Equation (13) and F( σ) = γ(σ) − σ 0 . It is seen that in the considered case of the plane stress, the function γ(σ) has all the properties of a norm; in particular, it vanishes only if all stress components vanish.

Remark 1.
In the plane strain problem of structures made of the materials satisfying Equation (14), the function γ(σ) does not have the properties of a norm, since the condition γ(σ) = 0 implies σ = αI, α ∈ R. Thus, within the theory of elasto-plasticity, there is a vital difference between the plane stress and plane strain cases. The results of the present paper cannot be transferred to the plane strain case; it would require an independent analysis. The elastic energy stored in the plate, expressed in terms of the virtual stress field τ, is given by: the operations Tr(·) and dev(·) are defined by Equation (10) for d = 2 and Equation (6). The bulk modulus k(x, y) and the shear modulus µ(x, y) are determined like in the classical theory of in-plane loaded plates; their units are N/m. Any virtual stress field τ must satisfy the equilibrium equations, both local and along the loaded boundary of the domain, hence it should satisfy the virtual work equation: where ε(v) is given by Equation (11), while f (v) represents the virtual work of loads. If the body forces are neglected and the tractions of intensity g are applied along the part Γ 1 of the contour ∂Ω, then: where s is the natural parameter of the contour Γ 1 . The variational equation implies: -the local equations of equilibrium -the static boundary conditions τ x n x + τ yx n y = g x , τ xy n x + τ y n y = g y on the contour Γ 1 . Such statically admissible stress fields form the set Σ(Ω). In the problem considered, the stresses undergo the plasticity condition: where γ(τ) is given by Equation (16). The set of stress fields satisfying Equation (23) will be denoted by K(Ω). Thus, τ ∈ Σ(Ω) ∩ K(Ω). According to the results [15] concerning the Hencky-Nadai-Ilyushin theory, the unknown stress field σ is the minimizer of the problem: The stress field σ is accompanied with the displacement fields u x , u y such that: where the components of the so-called plastic strains ε p x (x, y), ε p y (x, y), ε p xy (x, y) are not kinematically compatible; they are not associated with any displacement field, i.e., there does not exist a vector field v such that Equation (11) holds. The pair (σ, ε p ) satisfies the variational inequality: representing the celebrated Hill's principle of maximal plastic work. Having Equation (25), the equilibrium Equations (21) and (22) and the plasticity condition Equation (23) can construct the displacement fields u x , u y in the elasto-plastic structure. While the plasticity condition cancels extrema of the stress fields (see [22]), admitting the plastic components of strain degenerates the layout of the displacement fields (see [23]).
Let us emphasize once again that the considered 2D problem is viewed as the plane stress problem of statics of a plate of constant thickness b. Thus, the intensities of the tractions, the elastic moduli, the plastic limit σ 0 as well as the stress components are measured in N/m. The virtual work and the compliance have the units Nm.

The Isotropic Material Design (IMD) Method within the Elasto-Plastic Range
The aim is to construct the strongest transversely homogeneous plate made of the isotropic material of non-negative bulk and shear moduli; just these moduli are the only design variables of the problem. The unit cost of the design is assumed as trace of the Hooke tensor. In the 2D case, the eigenvalues of the Hooke tensor are: 2k, 2µ, 2µ, hence the unit cost is equal to 2k + 4µ. The total cost is bounded by a constant Λ 0 : We shall assume in the sequel that the permissible stress σ 0 does not depend on the design variables (k, µ). Thus, the optimum design problem assumes the form: Y being the compliance of the optimal structure. Let us insert Equation (24) and perform minimization over the design variables (k, µ), by making use of the sets Σ(Ω), K(Ω) being independent of the design variables. The operation of minimization over the design variables can be performed by using the minimization result (see [24]): min Ω a 1 u 1 + a 2 u 2 dxdy over u 1 , u 2 such that : in which a 1 ≥ 0, a 2 ≥ 0 are given functions in the domain Ω. The equality above is attained for:û Upon assuming: we find: where: Assume that the problem (P) is solvable upon appropriate mathematical modification; let τ * be the minimizer. The optimal moduli are expressed by: where E 0 = Λ 0 /|Ω|. It is easy to note that Equation (27) is satisfied sharply. One can prove that the stress field in the optimal plate (in which the elastic moduli are determined by Equation (35)) coincides with the stress field τ * solving the problem (P). Thus, the method put forward makes it possible to form the safely designed least compliant plate structure in which the stress field satisfies both the equilibrium equations and the yield stress condition (Equation (23)).

The General Form of the Problem
The IMD method requires the construction of the problem dual to (P) (see Equation (33)). To this end, we first release the constraints τ ∈ Σ(Ω), and by treating the virtual displacements in Equation (19) as Lagrange multipliers, we rearrange (P) to the form: The operations: min and max can be interchanged (see [25]), which makes it possible to re-write Equation (36) as below: In the next step, we shall find the explicit form of Equation (38); its form will not involve the parameter σ 0 .

Construction of the Potential h(ε) and the Explicit Formulation of the Problem Dual to (P)
By using Equation (9) for the scalar product of two tensors from E 2 s , taking into account Equations (34) and (16) and remembering that d = 2, we rewrite the local problem (Equation (38)) in the form: Let us introduce the notation: and re-write Equation (39) in the form σ 0 h(ε) = σ 0 h 1 (a, c, b), where: Let us introduce a new notation: Equation (41) simplifies to the form: where B(0,1) is a unit ball: x 2 + y 2 + z 2 ≤ 1. We see that the parameter σ 0 is cancelled. Now, we introduce the spherical parameterization: Operation max over ϕ gives: which simplifies Equation (43) to the form: where: Let us introduce a division of the set E 2 s into the subdomains: Remark 2. The set D 0 coincides with the set: where ρ o (·) is the function polar to ρ(·). Indeed, the function ρ o (·) has the form: derived in [14], which confirms the above observation.
The division (Equation (48)) of E 2 s into subdomains can be shown in the plane of principal strains. Let us recall that: Now, we are ready to show the explicit formula for the potential h(ε) defined by Equation (38), see Figures 1 and 2.   Figure 1). Moreover, it is convex, of linear growth outside the central domain D0, vanishes at = ε0 and is non-negative. In conclusion, the problem dual to (P) (see Equation (33)) assumes the form: One can prove that this value coincides with Equation (33), and the duality gap is zero. The pair (P), (P*) constitutes the LCP problem in the meaning of Bouchitté and Fragalà [13].  The function h(ε) is continuous, i.e., it is continuously stitched along the lines Tr ε = 1, Tr ε = −1, dev ε = √ 2 (see Figure 1). Moreover, it is convex, of linear growth outside the central domain D 0 , vanishes at ε = 0 and is non-negative. In conclusion, the problem dual to (P) (see Equation (33)) assumes the form: One can prove that this value coincides with Equation (33), and the duality gap is zero. The pair (P), (P*) constitutes the LCP problem in the meaning of Bouchitté and Fragalà [13].
Note that the locking domain is just the domain D0 given by Equation (48) (see Remark 2). Moreover, one can prove that the displacement field ( ) , xy uu in the optimal structure (whose moduli are given by Equation (35)) is proportional to the maximizer  v of Equation (53): Thus, the optimization process introduces the bounds on strains, while the values of stresses follow the values of the optimal elastic moduli. □

Construction of the Approximants of Statically Admissible Stresses
The optimal moduli , k   are determined by the solution to Equation (33). Therefore, it is thought appropriate to concentrate attention just on this problem and not on its dual form (Equation (53)). The aim of the present section is to show the numerical construction of sequences of sets ( ) h  approximating the set ( )  of statically admissible stresses, e.g., stresses equilibrating the given boundary traction load, hence satisfying the equilibrium Equations (21) and (22); index h symbolizes the mesh density parameter.
The description of the sequence of approximating sets ( ) h  needs specific notation, linked directly with the C++ programming syntax. The reader is asked to accept that the indices will start now from 0, not from 1. In particular, from now onward, the axes Note that the locking domain ρ o (ε(v(x, y))) ≤ 1 is just the domain D 0 given by Equation (48) (see Remark 2). Moreover, one can prove that the displacement field u x , u y in the optimal structure (whose moduli are given by Equation (35)) is proportional to the maximizer v * of Equation (53): Thus, the optimization process introduces the bounds on strains, while the values of stresses follow the values of the optimal elastic moduli.

Construction of the Approximants of Statically Admissible Stresses
The optimal moduli k * , µ * are determined by the solution to Equation (33). Therefore, it is thought appropriate to concentrate attention just on this problem and not on its dual form (Equation (53)). The aim of the present section is to show the numerical construction of sequences of sets Σ h (Ω) approximating the set Σ(Ω) of statically admissible stresses, e.g., stresses equilibrating the given boundary traction load, hence satisfying the equilibrium Equations (21) and (22); index h symbolizes the mesh density parameter.
The description of the sequence of approximating sets Σ h (Ω) needs specific notation, linked directly with the C++ programming syntax. The reader is asked to accept that the indices will start now from 0, not from 1. In particular, from now onward, the axes (x,y) will be denoted by (x 0 , x 1 ); consequently, we shall write f 0 , f 1 instead of f x , f y and f 00 , f 01 , f 10 , f 11 instead of f x , f xy , f yx , f y .
If Ω is a polygon, then the stress-based finite element method can be formulated as: find the interpolation σ h ∈ Σ h ⊂ Σ(Ω) of the statically admissible stress tensor field σ ∈ Σ(Ω), such that: where Dv represents the gradient of a vector field v and V h ⊂ V(Ω) is the finite elementwise subspace of functions υ h = (υ h0 , υ h1 ) : Ω → R 2 spanned by the polynomials of an appropriate degree. The P 1 (or Q 1 ) degree polynomials p = p(x) 11 x 0 x 1 ), p 00 , p 10 , p 01 (, p 11 ) ∈ R are adopted in this paper (see [26]). The finite element mesh in the domain Ω is composed of M 3-(or 4-node) finite elements Ω e ⊂ Ω covering the whole domain, provided it is a polygon. Let υ e h = υ h | Ω e be the truncation of υ h to the e-th element. Thus, the values υ h (x) = (υ h0 (x), υ h1 (x)) of the vector field υ h = (υ h0 , υ h1 ) ∈ V h truncated to the e-th element may be equivalently represented as two-dimensional vector: where υ e 0 , υ e 1 , . . . , υ e 2m , υ e 2m+1 are the unknown values of the scalar functions υ e hi (·), i = 0, 1 at three (or four) subsequent vertices of the triangle (quadrilateral) Ω e , while the polynomials ϕ e i : Ω e → R , i = 0, 1, . . . , m are the shape functions, which depend explicitly on the Cartesian co-ordinates z e i = z e i0 , z e i1 ∈ R 2 , i = 0, 1, . . . , m of the three (m = 2) or four (m = 3) vertices defining a triangular or quadrilateral finite element Ω e (see Figure 3). In the case considered, the formulae defining the shape functions in Equation (57) are relatively simple. However, even here it is thought appropriate to avoid using the functions ϕ e i (·) in Equation (57) and replace them by far more simple shape functions: φ i : ω → R, i = 0, 1, . . . , m , defined on the master element; in our problem, these functions are expressed by: for triangular reference (master) element and The implementation of the shape functions (Equations (58) and (59)) for an arbitrary Ω e element necessitates the introduction of a family of mappings F e = F e 0 , F e 1 : ω → Ω e , F e (ω) = Ω e , which link the master element ω with an arbitrary element Ω e such that ϕ e ∈ Ω e , ξ ∈ ω. This makes it possible to replace Equation (57) with a much simpler one: The geometric mapping F e : ω → Ω e is defined in a similar manner as the field υ e h has been constructed. Using the shape functions φ i and the Cartesian coordinates of nodes T of the finite element Ω e , we have the following simple relation:  The derivative of this mapping is a linear operator represented by the matrix: defined on master element ω (constant only for triangular element). On the basis of the easily calculated gradients ∇φ i (ξ) of the shape functions φ i (ξ), ξ ∈ ω, the gradients ∇ϕ e i (x), i = 0, 1, . . . , m of the shape functions ϕ e i (x), x = F e (ξ) ∈ Ω e are computed by: drawing upon the knowledge of the matrix (DF e (ξ)) −T being inverse-transpose to the matrix represented in Equation (62).
For the sake of simplicity, we assume that the load g applied to the boundary Γ 1 ⊂ Γ = ∂Ω of the design domain may have a different but constant value on selected sides of the polygon Ω, i.e.,∀x ∈ Γ g(x) = g = g 0 g 1 T = const, that is, we assume that a constant load is applied to the edge of any finite element, which is a fragment of the edge of the design domain Ω, possibly changing its value depending on the e-th number of the finite element Ω e . This allows us to assume that the vector g can be defined by three or four constants on each edge vector (see Figure 4): For the sake of simplicity, we assume that the load g applied to the boundary 1    =  of the design domain may have a different but constant value on selected sides of the polygon Ω, i.e., ( ) , that is, we assume that a constant load is applied to the edge of any finite element, which is a fragment of the edge of the design domain Ω, possibly changing its value depending on the e-th number of the finite element e  . This allows us to assume that the vector g can be defined by three or four constants on each edge vector (see Figure 4):  The calculation of the integral over the entire domain Ω and its boundary Γ (strictly Γ 1 ) in the variational Equation (56) can be reduced (as in classical, displacement-based FEM) to the calculation of the sum of the integrals over finite elements Ω e and their selected (i.e., loaded) boundaries Γ e , which coincide with the boundary Γ 1 ⊂ Γ ⊂ ∂Ω: Integration over Ω e and Γ e is shifted to the reference element ω and its boundaries ∂ω i , i = 0, 1, . . . , m. The left hand side is computed as follows: If the triangular element is used, the computation of the right hand side of Equation (56) is performed as below: while in the case of quadrilateral elements, the computation is performed in the way:  In each e-th finite element Ω e , the stress components (see Equation (71)) depend not only on ξ ∈ ω (i.e., x = F e (ξ) ∈ Ω e ) and on appropriately selected 3m + 3 indices i j , j = 0, . . . , 3m + 2 (from among all 3N indices {0, 1, . . . , 3N − 1}) defining local nodal stresses τ i j in e-th finite element, but additionally on s global parameters α k , k = 0, . . . , s − 1 defining the linear combinations of the s base vectors T k . In other words, upon constructing the solution (found only once) of linear, rectangular algebraic system B u T = Q u , one obtains a very simple approximation Σ α h of the statically admissible set of the stress fields Σ(Ω) determined by s global parameters α k ∈ R where in e-th finite element Ω e the following interpolations of the stress components hold:

Construction of the Approximate Solutions to the Problem (P) and Recovery of the Optimum Properties of the Initial Problem
The test fields τ ∈ Σ(Ω) of the problem (P) are interpolated by Equation (71) elementwise. These interpolations are x-dependent, which is underlined by now using the notation τ h (x, α).
Let us re-write Equation (23) in the form: According to the assumed stress field interpolation (Equation (71)), the discretized version of the problem (P) reads: find α * ∈ R s such that: Integration in Equation (73) is performed numerically on master element ω, i.e.: where here ξ = (ξ 0 , ξ 1 ) ∈ ω and w = w(ξ) are Gauss integration points and weights, respectively. In arbitrary element e and at arbitrary but fixed point ξ ∈ ω, the gradient: of the function R s α → ρ τ e h (ξ, α) ∈ R appearing in the mapping: can be computed by the rule: where: and: Equations (77)-(79) make it possible to calculate the quantity Π h given by Equation (76) and s components of its gradient for arbitrary design parameter α ∈ R s , i.e.: In arbitrary element e, at arbitrary point x ∈ Ω e and for arbitrary α ∈ R s , let us rewrite Equation (23) as: In arbitrary element e and at point x = F e (ξ) ∈ Ω e where ξ ∈ ω is arbitrary, the partial derivative of Equation (81) with respect to α k is equal to: where τ e hij = τ e hij (ξ, α) (i, j = 0, 1). For arbitrary p > 1, let us define the function: and write its derivative: In the algorithm for the numerical solution of the (P h ) problem proposed below, we assume that the yield condition in Equation (81) is satisfied at a finite number of points, i.e., at all Gaussian points. For this reason, we slightly modify the notation of the functional in Equation (81) and replace the lower index x symbolizing any point in Ω with a subscript denoting the successive Gaussian points counted in subsequent finite elements Ω e , e = 0, 1, . . . , M − 1, i.e.: where x g = F e (ξ) ∈ Ω e is the g-th image of the Gauss point ξ ∈ ω in the master element. The index g runs from 0 to G = m × M − 1, where m represents the number of Gauss points in ω. We will also omit the superscript e identifying the number of finite elements. Now, we are ready to formulate the algorithm for solving the (P h ) problem: Step 0. Find a solution T of the static problem B u T = Q u . From now, the design parameter is the vector α = [α 0 α 1 . . . α s−1 ] T ∈ R s .
Step 4. Starting with α = α 0 , apply any algorithm of the nonlinear mathematical programming to find the solution α * = argmin α∈R s f (α) ∈ R s of the unconstrained prob- where the function f (α) and its gradient ∇ f (α) are defined by Equations (86) and (87), respectively.
Step 5. If P k(α * ) < ε then STOP, otherwise calculate the new value of the penalty parameter as P = χ P and initialize design parameter α 0 = α * . Go to Step 4.
The approximants of the problem (P) (see Equation (33)) computed by the above algorithm will be denoted by Π * . The quantity Y * will represent approximants of the optimal compliance Y (see Equation (32)).

Case Studies and Discussion
In the analysis of plate structures loaded in the plane, deforming within the linear elastic range, it is impossible to prevent singularities of stresses around critical points or along some lines. These points are reentrant corners, places where the load is concentrated or where the boundary conditions change abruptly and the structure loses its support. One can achieve better control over the stress level if the structure is not supported and the load is self-equilibrated; however, such problems are usually not practical. The stress-based LCP problem of the IMD method within the elastic range (the specific case of problem (P), Equation (33), with the yield condition being neglected) also suffers from the drawback of the possible appearance of stress singularities. Thus, according to Equation (35), the optimal moduli blow up at these places. To be more precise, the bulk modulus becomes infinite where the trace of stress tensor is singular; the shear modulus blows up where the norm of the stress deviator tends to infinity. In the plane stress problem considered, the HMH condition assumes the form of Equation (15). That is why introduction of the yield condition (Equation (23)) alleviates all components of stress. Thus, one can expect that the condition (Equation (23)) in the IMD setting should bring about cutting all extremes of all components of the stress field solving the auxiliary problem (P), hence making regular all layouts of the optimal elastic moduli.
Two optimum design problems are considered: -Designing the material layout within the rectangular cantilever plate (of the in-plane dimensions 2L by 4L, see Figure 5a) subjected to a lateral constant traction of intensity g x : Examples 1a, 1b, 1c; - The optimum design of the L-shaped plate, see Figure 5b, subjected to the vertical shearing traction along one vertical side: Examples 2, 3.
Within the purely elastic IMD method, the optimum cantilever plate suffers singular layouts of the moduli around the left and right ends of the support. Due to the linear elastic approach, the moduli are proportional to the magnitude of the load, while the shape of the layout is load-independent. The plastic version of the IMD introduces an essential change: the optimal layout of the moduli does depend upon the ratio: g x /σ 0 , hence the layout of the optimal moduli becomes dependent on the magnitude of the load. We also have control over the size of plastic zones. One of the aims of the present paper is to analyze sequences of the optimal designs corresponding to various values of the ratio g x /σ 0 .
The optimal designs of the rectangular cantilever plate have been constructed by using the special software based on the numerical scheme outlined in Section 5. Two kinds of finite elements are used: the triangular (T) and quadrilateral (Q) described in Section 5. Both FE meshes are regular.
The same software has been used to design the optimal moduli within: the L-shaped plate of sharp corners (see Figure 5b,c) and within the L-shaped plate with the reentrant corner being slightly rounded (see Figure 5e). The plate in Figure 5b is meshed by quadrilateral finite elements; the mesh is regular. The meshes for the plate in Figure 5c The value of the referential modulus E 0 (this is not Young's modulus, its units are N/m) appearing in the isoperimetric condition (Equation (27)) is assumed now as E 0 = 2k 0 + 4µ 0 , where: are characteristic bulk and shear stiffnesses of the plate of thickness b, made of the referential homogeneous material with moduli E, ν. The values of the remaining parameters appearing in the penalty function algorithm are adopted as follows: The last quantity ftol is a parameter used in the gradient-oriented frprmn( . . . ) procedure in C++ (see [27] implementing the Fletcher-Reeves-Polak-Ribiere algorithm of the minimization of functions without constraints). Numerical integration has been performed for the master element on the basis of the rules of integration with one and four Gauss points for triangular (T) and quadrilateral (Q) finite elements, respectively. All the data are now given, and the results are ready to be replicated.

-
Designing the material layout within the rectangular cantilever plate (of the in-plane dimensions 2L by 4L, see Figure 5a) subjected to a lateral constant traction of intensity gx: Examples 1a, 1b, 1c; - The optimum design of the L-shaped plate, see Figure 5b, subjected to the vertical shearing traction along one vertical side: Examples 2, 3.
(e) (f) Case 1a. The lateral horizontal traction of intensity g x = 0.01 · σ 0 applied to the left edge (see Figure 5a).
The optimum design problem (Equation (28)) has been solved by applying the numerical method outlined in Sections 5 and 6. The two regular FEM meshes composed of 34 × 69 = 2346 quadrilateral and 68 × 69 = 4692 triangular finite elements were used. It has occurred that for sufficiently dense FEM meshes, the results obtained for triangular and quadrilateral elements are practically identical (see Figures 6 and 7). For this reason, the next results of optimal distributions of elastic moduli will be presented for a mesh spanned only by quadrilateral or only by triangular finite elements. The optimum design problem (Equation (28)) has been solved by applying the numerical method outlined in Sections 5 and 6. The two regular FEM meshes composed of 34 × 69 = 2346 quadrilateral and 68 × 69 = 4692 triangular finite elements were used. It has occurred that for sufficiently dense FEM meshes, the results obtained for triangular and quadrilateral elements are practically identical (see Figures 6 and 7). For this reason, the next results of optimal distributions of elastic moduli will be presented for a mesh spanned only by quadrilateral or only by triangular finite elements.    The optimal layouts of the moduli k * , µ * have been constructed by Equation (35), and the moduli E * , ν * are computed by: (see Figures 6 and 7). Because the traction is small, Equation ( The optimal layouts of the moduli k * , µ * , E * , ν * have been constructed (see Figure 8)  Zero (or numerically close to zero) values of the optimal moduli k * and µ * mean in practice the need to cut off these sub-areas from the entire Ω domain. In Figure 9 the same as in Figure 8, the optimal distributions of elastic moduli are shown with a clearly visible modification of the optimal shape of Ω consisting of cutting off the right upper corner of the cantilever at all those points where both optimal values of k * and µ * are equal to zero or are numerically close to zero. However, the correct cutting off of the material inside the Zero (or numerically close to zero) values of the optimal moduli k * and µ * mean in practice the need to cut off these sub-areas from the entire Ω domain. In Figure 9 the same as in Figure 8, the optimal distributions of elastic moduli are shown with a clearly visible modification of the optimal shape of Ω consisting of cutting off the right upper corner of the cantilever at all those points where both optimal values of k * and µ * are equal to zero or are numerically close to zero. However, the correct cutting off of the material inside the design domain cannot be easy programmed. For this reason, in the further examples, the empty domain within the design domain will not be cut off. The optimal layouts of the moduli k * , µ * , E * , ν * have been constructed (see Figure 10 Having constructed the optimal designs for three subsequent, increasing magnitudes of the lateral traction, one can discuss the influence of the parameter g x /σ 0 on the final solutions. Along with the increase in the lateral load g x , one can observe that those zones of the design domain Ω expand, in which the optimal moduli k * and µ * assume high or moderate values; those zones are shown in orange and red. This is very visible while comparing the layouts of the moduli k * and µ * in the vicinity of the lower vertices of the design domain, along the lower horizontal edge and both the vertical edges. In the case of the small load, which does not induce the plastic zones within the design domain, making the mesh denser causes the shrinking of the zones of high values of the optimal elastic moduli (see Figure 6) to several finite elements (whose dimensions are smaller and smaller if the mesh is made denser) around the corners. Just in these elements, the values of the optimal moduli grow up, thus making the cost condition satisfied. These values tend to infinity along with making the mesh denser and denser. By the introduction of the plastic limit within the whole design domain, we ban the mentioned tendency to accumulate the high values of the optimal moduli around some points; the zones of high values of the moduli become broader along with the expansion of plastic zones. This tendency is easy to verify by comparing the optimal layouts of the elastic moduli shown in Figures 6-10. The plastic zones are places where the γ function attains the upper bound-see places in yellow in Figure 8, Figure 9, where the plot of the function γ becomes flat. Let us note that the intensity of the load can be increased only up to a certain limit; if this limit is exceeded, the problem (P h ) ceases to be solvable. Moreover, it is worth stressing that in each case (presented in Figures 6-10) the optimal Poisson ratio assumes the values from the whole admissible range: (−1,1); in particular, the auxetic zones (with negative Poisson ratio) appear in all cases, where necessary. In the case of the appearance of optimal plastic zones, the shape of the sub-domains in which the Poisson ratio remains negative changes slightly, always keeping the full range of its extremely small negative values. This can be partially explained by recalling the well-known properties of auxetic materials, in particular those related to the influence of negative Poisson's ratio on the values of the stress concentration factor in the design of body components subjected to stress: "When the Poisson's ratio becomes negative, stress concentration factors are reduced in some situations and unchanged or increased in others."-see [28]. The results of many studies suggest that very often (but not always) a negative Poisson's ratio gives the lowest possible (i.e., the most desirable) value of the stress concentration factor, which can be, in an analogous way, justified by our numerical results of optimal distributions of elastic moduli minimizing the compliance of the elasto-plastic body with a simultaneous demand to meet the Mises plasticity condition at all points within the design domain Ω. However, the study does not analyze the impact of the optimal auxetic sub-domains on the values of the stress concentration factors. Many very interesting results on this subject can be found, e.g., in the monograph [8]. Having constructed the optimal designs for three subsequent, increasing magnitudes of the lateral traction, one can discuss the influence of the parameter gx/σ0 on the final solutions. Along with the increase in the lateral load gx, one can observe that those zones of the design domain Ω expand, in which the optimal moduli k * and µ * assume high or moderate values; those zones are shown in orange and red. This is very visible while  The design problem has been solved with the use of regular and irregular FE meshes composed of 2523 quadrilateral or 5833 triangular finite elements, respectively. The Lshaped cantilever is loaded with the vertical tangent traction of intensity g y = 0.1 · σ 0 applied to the right lower vertical edge. The optimal layouts of the moduli k * , µ * ,E * , ν * have been constructed (see Figure 11) of the stress concentration factor, which can be, in an analogous way, justified by our numerical results of optimal distributions of elastic moduli minimizing the compliance of the elasto-plastic body with a simultaneous demand to meet the Mises plasticity condition at all points within the design domain Ω. However, the study does not analyze the impact of the optimal auxetic sub-domains on the values of the stress concentration factors. Many very interesting results on this subject can be found, e.g., in the monograph [8].

Example 2.
Optimum design of the L -shaped cantilever plate (see Figure 5b-d).
The design problem has been solved with the use of regular and irregular FE meshes composed of 2523 quadrilateral or 5833 triangular finite elements, respectively. The Lshaped cantilever is loaded with the vertical tangent traction of intensity gy = 0.1 • σ0 applied to the right lower vertical edge. The optimal layouts of the moduli   The plate is covered with an irregular mesh of 5803 triangular finite elements. The cantilever is loaded with the vertical tangent traction of intensity g y = 0.1 · σ 0 applied to the right lower vertical edge. The optimal layouts of the moduli k * , µ * , E * , ν * have been constructed (see Figure 12) All the remarks concerning the interpretation of the results concerning the rectangular cantilever apply here. Moreover, by making the reentrant corner curve smoothly, we alleviate the stress concentration, thus making the optimal Young modulus and Poisson's ratio layouts much more regular (see Figures 11 and 12).

Conclusions
The hitherto existing works on topology optimization enhanced with local stress constraints have been formulated within the elastic range: on the stress components, being associated with the displacement field, the local constraints are imposed; they can concern All the remarks concerning the interpretation of the results concerning the rectangular cantilever apply here. Moreover, by making the reentrant corner curve smoothly, we alleviate the stress concentration, thus making the optimal Young modulus and Poisson's ratio layouts much more regular (see Figures 11 and 12).

Conclusions
The hitherto existing works on topology optimization enhanced with local stress constraints have been formulated within the elastic range: on the stress components, being associated with the displacement field, the local constraints are imposed; they can concern all the components of stresses (see [29]) or the effective stress (see, e.g., [30]). In the present paper, another formulation of the topology optimization problem is set forth: the Hencky-Nadai-Ilyushin elasto-plastic theory is adopted in which the stress state is not linked directly with the displacement field. Thus, the optimal structure (here: an in-plane loaded plate) works within the elasto-plastic range. Consequently, the optimal design does depend upon the ratio: intensity of the load/yield stress. One of the aims of this paper is to analyze the variation of the design for a given load, if the yield stress level varies. It turns out that the approximants Π * of the optimal compliance calculated for subsequent values of the plasticity limit and fixed intensity of the traction load decrease with the increasing value of yield stresses σ 0 (see Figure 13a). If for an assumed intensity of the traction load the yield stress is taken too small, it is not possible to attain the minimum Π * of the mapping Π h , which means that an optimal solution does not exist. Similar conclusions hold in the case of increasing the load g x for the assumed constant value of the yield stress σ 0 (see Figure 13b). modulus k * /b, (c) shear modulus µ * /b, (d) Young's modulus E * /b, (e) Poisson's ratio ν * in the case of FEM mesh composed of triangular elements.
All the remarks concerning the interpretation of the results concerning the rectangular cantilever apply here. Moreover, by making the reentrant corner curve smoothly, we alleviate the stress concentration, thus making the optimal Young modulus and Poisson's ratio layouts much more regular (see Figures 11 and 12).

Conclusions
The hitherto existing works on topology optimization enhanced with local stress constraints have been formulated within the elastic range: on the stress components, being associated with the displacement field, the local constraints are imposed; they can concern all the components of stresses (see [29]) or the effective stress (see, e.g., [30]). In the present paper, another formulation of the topology optimization problem is set forth: the Hencky-Nadai-Ilyushin elasto-plastic theory is adopted in which the stress state is not linked directly with the displacement field. Thus, the optimal structure (here: an in-plane loaded plate) works within the elasto-plastic range. Consequently, the optimal design does depend upon the ratio: intensity of the load/yield stress. One of the aims of this paper is to analyze the variation of the design for a given load, if the yield stress level varies. It turns out that the approximants П * of the optimal compliance calculated for subsequent values of the plasticity limit and fixed intensity of the traction load decrease with the increasing value of yield stresses σ0 (see Figure 13a). If for an assumed intensity of the traction load the yield stress is taken too small, it is not possible to attain the minimum П * of the mapping Пh, which means that an optimal solution does not exist. Similar conclusions hold in the case of increasing the load gx for the assumed constant value of the yield stress σ0 (see Figure 13b).
The research planned will concern the design of the underlying microstructures exhibiting the given effective yield limit, characterized by the effective moduli predicted by the IMD method.  Data Availability Statement: Data sharing is not applicable to this article.

Conflicts of Interest:
The authors declare no conflict of interest.