TOPOLOGICAL DERIVATIVE-BASED TOPOLOGY STRUCTURAL OPTIMIZATION UNDER DRUCKER-PRAGER-TYPE STRESS CONSTRAINTS

An algorithm for topology optimization of elastic structures under plane stress subject to the Drucker-Prager stress constraint is presented. The algorithm is based on the use of the topological derivative of the associated objective functional in conjunction with a levelset representation of the structure domain. In this context, a penalty functional is proposed to enforce the point-wise stress constraint and a closed formula for its topological derivative is derived. The resulting algorithm is of remarkably simple computational implementation. It does not require post-processing procedures of any kind and features only a minimal number of user-defined algorithmic parameters. This is in sharp contrast with current procedures for topological structural optimization with local stress constraints. The effectiveness and efficiency of the algorithm presented here are demonstrated by means of numerical examples. The examples show, in particular, that it can easily handle structural optimization problems with underlying materials featuring strong asymmetry in their tensile and compressive yield strengths.


INTRODUCTION
Over the last two decades or so, the development of algorithms for topology optimization of linear elastic load-bearing structures has attracted considerable attention in computational mechanics circles.As a result of the continuous research efforts in this direction a wide body of literature is currently available on this topic and various computational procedures are well established and can be applied to a range of practical problems of industrial interest [16,1,10].Many such procedures, almost invariably used in conjunction with finite element methods of structural analysis, are even available in off-the-shelf commercial software packages.
To date, most developments in this field have relied on so-called SIMP methods (solid isotropic material with penalization), where the physical black-and-white topology of the optimal structure, i.e. a topology consisting of either material (black) or empty space (white) at each point of the computational domain, is approximated by means of a fictitious density field displaying a smooth (grey) transition in the otherwise black-white interface (the boundary of the structure domain).Such methods have been widely applied with success to problems such as compliance minimization [10] but, despite its fundamental importance in engineering design, only a relatively small number of publications appear to deal with the incorporation of local (point-wise) stress constraints [2,3,11,14,18,21,27].This can be probably justified by the challenges resulting from the typically very large number of highly non-linear constraints involved as well as by the need for carefully designed stress relaxation procedures to address a side effect of the density field-based regularization of the original black-and-white problem [21].
More recently, a new class of methodologies for structural topology optimization has emerged based on the use of the topological derivative of the relevant objective functionals [29,12,25,5,24,26].The notion of topological derivative itself is a relatively new concept, introduced just over a decade ago [12,29].Further theoretical developments are reported, among others, in [23,4,28].An early application of this idea to topology compliance optimization, prior to its precise mathematical definition in a general context, is described in [17].The topological derivative concept extends the conventional notion of derivative to functionals whose variable is a geometrical domain subject to singular topology changes.In structural topology optimization for instance, it gives the exact sensitivity of the associated objective functionals to black-and-white-type topological perturbations such as the insertion of infinitesinal holes or inclusions of different material properties.Crucial here is the fact that the topological derivative of the objective functional contains fundamental information that accurately indicates descent directions associated with exact black-and-white-type topology changes, without the need for black-grey-white-type regularisation procedures.In this context, a topological derivative-based algorithm with a level-set representation of the structure domain has been proposed in [5] and shown to efficiently solve compliance minimization problems.More recently, following the ideas presented in [6] for the Laplace equation, this algorithm has been further developed in [8] to incorporate local stress constraints of the von Mises type by means of a penalty approach in plane stress problems.One striking feature of the algorithm of [8] is its simplicity of implementation.Once a suitable penalized objective functional has been defined and a closed formula for its topological derivative obtained, the locally stress-constrained topology optimization problem is treated algorithmically in exactly the same way as its unconstrained counterpart.In particular, in contrast to current methods of stress-constrained topology optimization, the topological derivative-based procedure does not require post-processing (e.g.procedures such as density filtering, ε-relaxation [21]) of any kind and only a minimal number of user-defined algorithmic parameters (e.g.penalty coefficient) are needed.This relative algorithmic simplicity is nothing but a natural consequence of the use of the topological derivative in defining the descent direction, which is based on the exact black-and-white definition of the topology optimization problem.In fairness to other methods of topology optimization, however, we should note that the striking algorithmic simplicity here comes at the cost of derivation of a closed formula for the topological derivative of the objective functional which may prove to be a laborious mathematical task.
Our main purpose in this paper is to extend the work reported in [8] to incorporate point-wise stress constraints of the Drucker-Prager type [15].In particular, we want to minimize the volume of the structure domain requiring at the same time the stress tensor at each point of the loaded optimized structure to be bound by a Drucker-Prager-type yield criterion.In this context, a suitable penalty functional for the enforcement of the Drucker-Prager constraint is proposed and a closed formula for its topological derivative is obtained.We recall that the Drucker-Prager yield criterion was originally conceived as a smooth approximation to the classical Mohr-Coulomb criterion for soils and geomaterials (refer for instance to [13]).Under plane stress (the case considered here) it may be used as a general model for materials with distinct tensile and compressive yield strengths, such as concrete, masonry and wood.The overall optimization algorithm is described in detail and numerical examples are presented to demonstrate its effectiveness and efficiency in the treatment of structural optimization under the present stress constraints.In particular, unlike stress-unconstrained optimization, the results here show that the obtained optimized structures are free from geometrical singularities that result in (highly undesirable) stress concentration.
The paper is organized as follows.Section 2 states the stress-constrained topology optimization problem and defines the penalized version to be solved by the algorithm.Section 3 presents a closed formula for the topological derivative of the corresponding penalized objective functional.The optimization algorithm is described in Section 4 and its application in numerical examples is presented in Section 5. Concluding remarks are drawn in Section 6.

THE TOPOLOGY OPTIMIZATION PROBLEM
Our purpose here is to find optimal topologies for two-dimensional elastic structures under plane stress condition loaded by a given system of mechanical loads with prescribed kinematical boundary conditions and subject to a point-wise constraint on the stress tensor.More specifically, we want to minimize the volume of the structure domain requiring at the same time the stress tensor at each point of the loaded optimized structure to be bound by a Drucker-Prager-type yield criterion.The corresponding optimization problem is mathematically stated in the following.

The constrained optimization problem
Let D ⊂ ℜ 2 be a bounded domain with Lipschitz boundary Γ defining the so-called hold-all domain (refer to Fig. 1).The domain of the sought optimal structure will be a subset of the hold-all domain.The boundary Γ is the union of three given non-overlapping subsets, Γ D , Γ N and Γ 0 .Displacements are prescribed on Γ D and non-zero and zero boundary tractions are prescribed respectively on Γ N and Γ 0 .In addition, we conveniently assume that the stress constraint is to be enforced on a given open subset D of D. Note that stress constraints cannot usually be enforced, for instance, in the surroundings of point supports or point loads and, hence, D ̸ = D in general.
Given a hold-all domain D and a stress constraint-enforcement subdomain D, the optimisation problem consists in finding a subdomain Ω ⊂ D (the optimal structure domain) that solves the following constrained minimization problem: with I Ω the objective functional subject to the elastic equilibrium equations, and a point-wise Drucker-Prager constraint on the stress tensor Σ: with Σ M the von Mises effective stress: and Σ d the stress deviator.The given scalar constants η and σ ⋆ in (4) are the Drucker-Prager yield criterion parameters [15,13] associated, respectively, with the Drucker-Prager cone angle and cohesion intersect.In (2,3), g is the prescribed boundary traction field on the given portion Γ N of the boundary and is assumed to belong to L 2 (Γ N ) 2 , n in (3) is the outward unit normal vector field on Γ and u Ω is the displacement field that solves the elastic equilibrium equations.The objective functional defined in (2) is well-suited for the minimization of the volume |Ω| of the structure subject to a point-wise stress constraint and has been used in [8] in conjunction with a von Mises stress constraint.The parameter β > 0 multiplying the compliance integral on the right hand side of (2) regularises the stress-constrained volume minimization problem which is otherwise ill-posed.The subscript Ω is used here to emphasise that the relevant quantities (e.g.I Ω , u Ω ) depend on the domain Ω -the design variable of problem (1).Throughout the paper, we assume (3) to hold in the weak sense and its solution, to be unique.The space V is the corresponding space of kinematicaly admissible displacement fields.The notation Σ(u Ω ) is used to emphasize that the stress tensor is a functional of the displacement field u Ω through the linear elastic constitutive equation: where e is the infinitesimal strain tensor, and C = 2µ Ω I + λ Ω (I ⊗ I), (9) with µ Ω and λ Ω denoting the Lamé coefficients and I and I the fourth-and second-order identity tensors respectively.The statement of the minimization problem is completed with the definition of a piece-wise constant Young's modulus field over D as follows: with That is, the original optimization problem, where the structure itself consists of the domain Ω of given elastic properties and the remaining part D \ Ω of the hold-all is empty (has no material), is approximated by means of the two-phase material distribution (10) over D where the empty region D \ Ω is occupied by a material (the soft phase) with Young's modulus, E soft , much lower than the given Young's modulus E hard of the structure material (the hard phase).Both phases share the same Poisson's ratio ν.The corresponding Lamé coefficients under plane stress read Figure 1.Sketch of the hold-all domain.

The penalized optimization problem
The presence of the point-wise stress constraint (4) makes it difficult to treat the above constrained optimization problem directly.This issue has been recently discussed in some detail by Le et al. [21] in the context of SIMP methods for structural optimization [10].To tackle the problem here we follow a radically different approach proposed in [8].It relies on a topological derivative-based algorithm in conjunction with an approximation of the original constrained problem by means of a penalty regularization of the point-wise stress constraint.The penalized problem is obtained in the following.
Before defining the corresponding penalty functional it is convenient in the present case to re-phrase the stress constraint (4) in terms of normalized quantities.To this end we define the normalized stress tensor: and the normalized cohesion intersect-related parameter of the Drucker-Prager yield surface: By rewriting (4) as squaring both sides and making use of the above definitions, we obtain after a straightforward manipulation an equivalent statement of the Drucker-Prager stress constraint in terms of normalized stresses: where Alternatively, by taking the elastic law (7,9) into account, ( 16) can be expressed as where with µ and λ the normalized Lamé coefficients: and With the above at hand, we now proceed to define the penalized objective function.Then, let Φ : ℜ + → ℜ + be a nondecreasing function of class C 2 .To allow a proper justification in the analysis, we further assume that the derivatives Φ ′ and Φ ′′ are bounded.The penalty functional is defined as According to (4), in the original optimization problem the stress is constrained in Ω ∩ Dthe portion of the elastic structure subject to the stress constraint.In the penalized version ( 22) -where a soft phase has been introduced to mimic the void region -the constraint must be imposed over the entire D. With the above penalty function, we define a corresponding penalized objective functional as where the scalar α > 0 is a given penalty coefficient.The original constrained optimization problem ( 1)-( 4) with point-wise constraints can then be approximated by the following penalized optimization problem: Problem ( 24) provides a good approximation to (1)-( 4) so long as (a) the penalty coefficient α is sufficiently large; and In particular, in the present paper we shall adopt a function Φ of the following functional format: where p ≥ 1 is a given real parameter and Φ p : ℜ + → ℜ + is defined as With this choice, the penalized problem (24) to be solved here reads explicitly

TOPOLOGICAL DERIVATIVES
The unconstrained minimization problem (27) will be solved in this paper by the algorithm described in Section 4, which relies fundamentally on the concept of topological derivative.This section provides a closed formula for the topological derivative of the penalized objective functional of ( 27) to be used in the algorithm.Before presenting the closed formula itself, a brief discussion on the relatively recent concept of topological derivative appears to be convenient and should be helpful to those not yet familiar with the idea.

The topological derivative concept
The notion of topological derivative extends the conventional definition of derivative to functionals whose variable is a geometrical domain subjected to singular topology changes.The idea can be introduced by considering a generic functional G(Ω) of a given domain Ω and assuming that Ω is subject to topology changes consisting, say, of the introduction of a circular hole of radius ε centered at an arbitrary point x ∈ Ω.The resulting topologically changed domain, denoted Ω ε ( x), is the set defined as (refer to Fig. 3) Figure 3.An example of topological domain perturbation.
where B ε ( x) denotes the closure of the domain of the inserted hole.The topological derivative of the functional G exists if its value G(Ω ε ) for the topologically perturbed domain Ω ε can be expressed as a sum of the functional G(Ω) evaluated for the original domain Ω, a term f (ε) D T G( x) that varies linearly with a function f (ε) and a remainder of the form o(f (ε)).The function f : ℜ + → ℜ + must be such that f (ε) → 0 when ε → 0 + and the remainder o(f (ε)) vanishes faster than f (ε), with respect to ε, namely: The right hand side of ( 29) is named the topological asymptotic expansion of G and the field D T G : Ω → ℜ is the topological derivative of the functional G evaluated at the original domain Ω for the considered type of topological perturbation (the introduction of a circular hole).The topological derivative D T G itself can be expressed as The analogy between (29,31) and the corresponding expressions for a conventional derivative should be noted.

The topological derivative of the penalized objective functional
In the minimization problem (27) the hold-all domain is split as the union of a subset Ω occupied by the hard phase and its complement D \ Ω occupied by the soft phase.In this case, it is appropriate to consider topological perturbations consisting of the introduction a circular inclusion of domain B ε ( x) made of hard phase material if the perturbation point x lies in the soft phase domain and made of soft phase material if x lies in the hard phase domain.The corresponding perturbed structural domain Ω ε ( x), i.e. the domain of the hard phase after the introduction of the inclusion, reads The topological derivative of the unconstrained objective functional ( 27) is given by the sum of topological derivatives of each term on the right hand side of (27) with respect to the class of topological perturbations defined by (32).The first term D T |Ω| above is trivial.Here we have The topological derivative D T K Ω of the compliance functional is known.It has been used in the context of structural optimization with topological derivative-based algorithms (refer, for instance, to [4,19] for a detailed derivation).Its closed formula is where the scalar ρ is the fourth-order tensor T is the polarization tensor given by with γ the elastic modulus contrast and the constants a and b given by The derivation of the topological derivative D T J Ω of the penalty functional (22) for the Drucker-Prager stress constraint is rather involved.For the sake of clarity we limit ourselves to presenting only the final formula here.The closed formula for D T J Ω reads where with χ D the characteristic function of D: The functions ζ 1 , ζ 2 and Ψ ρ are given by and with where σ I and σ II are the eigenvalues of σ(u Ω ).The vector field v Ω in (41) is the solution of the adjoint equation (48) Formula ( 41) is valid for all x ∈ D \ ∂ D \ ∂Ω.

THE TOPOLOGY DESIGN/OPTIMIZATION ALGORITHM
The numerical solution of the penalized minimization problem ( 27) is undertaken here by the algorithm proposed in [5] in conjunction with a finite element approximation of the elastic boundary value problem (3) and the adjoint equation (48).The algorithm relies essentially on an optimality criterion based on the topological derivative of the objective function and on a level-set representation of the structure domain.It was proven very successful in the context of unconstrained structural optimization and optimization in problems of flow through porous media [5], in structural optimization under a von Mises stress constraint [8] and in the topology optimization of elastic microstructures [7].
With the level-set representation, the current structure domain Ω is characterized by a level-set function ψ ∈ L 2 (D) as and its complement as (50)

Topological derivative-based local optimality condition
The establishment of a local optimality condition based on the topological derivative field is straightforward.Indeed, note that for any given structure domain Ω, a negative (positive) value of the topological derivative D T I α Ω (x) at an arbitrary point x ∈ D indicates that the introduction of an infinitesimal circular inclusion centered at that point produces a perturbed domain whose objective functional value is smaller (greater) than that of the original domain.From this observation we have that a sufficient condition of local optimality under the considered class of topological perturbations is that That is, (51) implies in the present context that the introduction of an infinitesimal circular inclusion at any point of D can only cause an increase in the value of the objective functional and, hence, Ω is indeed a locally optimum structure domain.
The algorithm described below is based on an alternative sufficient condition of local optimality, particularly convenient for use in conjunction with the level-set representation of the structure domain.The alternative condition can be established by first defining the scalar function, and then observing that, exclusively in terms of the function g and the level-set ψ, condition (51) is equivalent to Now, we note that (53) holds if the function g is a strictly positive scalar multiple of ψ, i.e.
or, equivalently, where θ is the angle between the functions g and ψ in L 2 (D).Hence, (54) or ( 55) are also sufficient conditions of local optimality.In particular, (55) will be used in the algorithm described below.

The algorithm
The algorithm itself aims to generate a sequence {ψ i } of level-set functions (a sequence of structural domains {Ω i }) that will produce for some iteration n a domain Ω n such that (55) is satisfied to within a given small numerical tolerance ϵ θ > 0: The iterative procedure starts with the choice of an initial guess for the optimal structure domain, i.e. with the choice of a starting level-set function ψ 0 ∈ L 2 (D).For simplicity, the function ψ 0 is chosen as a unit vector of L 2 (D).With S denoting the set of unit vectors of L 2 (D), the algorithm is explicitly given by ψ 0 ∈ S, where i denotes a generic iteration number and κ i ∈ [0, 1] is a step size determined by a linesearch performed at each iteration in order to decrease the value of the objective functional I α Ω i .Note that the right hand side of (57) 2 is a convex combination between ψ i−1 and g i−1 up to a positive multiplicative constant and that, by construction of the iteration formula, we have The iterative process is stopped when for some iteration i the optimality condition (56) is satisfied to the desired degree of accuracy, i.e. if If at some iteration i, the line-search step size κ i is found to be smaller than a given numerical tolerance ϵ κ > 0: i.e. if the topology is effectively no longer changing, and at the same time the optimality condition (59) is not met, then a uniform mesh refinement of the hold-all domain D is carried out and the iterative procedure is continued.The goal of this coarse-to-fine procedure is twofold: on one hand to reduce the computer cost, and on the other hand to be less subject to getting stuck in bad local minima.Of course, as the continuous topology optimization problem has no global minimum in general, the mesh dependency of the obtained domain cannot be completely avoided.Note that a uniform mesh refinement will enrich the discretized space of level-set functions and will in general improve the accuracy of the elastic solutions.At this point we should add that more sophisticated approaches could be adopted, for example, by generating at each iteration a finite element mesh having element edges matching exactly the phase interface defined by the corresponding level-set function.Also, mesh density could be defined according to a suitable error estimator, leading to potential savings in computing time and ensuring the error is bounded to a desired level in the attained optimal structure domain.We emphasize, however, that our main purpose here is to demonstrate the robustness of the topological derivative-based iteration (57).Hence, we choose to rely on a simpler approach to avoid any potential lack of robustness being masked by procedures of a more peripheral nature.
In the computation of D T I α Ω according to expression (33) the stresses and the topological derivatives are first computed within the finite elements (at Gauss points) and then extrapolated to nodes.The final discretized version of the field D T I α Ω used in the iterations is generated by the finite element shape functions with smoothed nodal values obtained in a standard fashion.The level-set functions ψ and the discretized field D T I α Ω are generated by the same shape functions used in the finite element approximation of the direct and adjoint boundary value problems (3) and (48).The material properties E hard or E soft are assigned to nodes of the mesh depending on whether they are at points with ψ < 0 (hard phase) or ψ > 0 (soft phase).In this way, elements crossed by the hard-soft phase interface (defined by ψ = 0) will have Young's moduli between the values E hard and E soft , obtained by a standard interpolation of the nodal Young's moduli using the element shape functions.Obviously, according to the above procedure, the resolution of the optimal structure domain depends directly on the fineness of the adopted mesh.

NUMERICAL EXAMPLES
The effectiveness of the algorithm described above is demonstrated in this section by means of numerical examples.In order to avoid numerical ill-conditioning of the optimization problem we use in all examples, without loss of generality, a normalized version of the objective functional of (27) defined as with the normalizing factors V 0 and K 0 being respectively the area and the compliance functional of the initial guess Ω 0 for the optimum structure domain, here taken as Ω 0 = D.In all the examples, we adopt the Young's modulus contrast E soft /E hard = 10 −3 .

Wall under shear load
The first example consists of wall under shear load (see Fig. 4).The hold-all domain is a rectangle of size 2×1 clamped at its bottom edge.The loading consists of a unit uniformly distributed horizontal force g = (1, 0) applied along a central portion of length 0.2 of the top edge of the hold-all domain.The material parameters E hard = 1.0, ν = 0.3 and σ = 1 are used.For the penalty coefficient and compliance weighting factor we choose α = 25 with p = 32 and β = 1/4.The optimization procedure is carried out for three different values of η.Firstly we use η = 0, corresponding to a von Mises stress constraint and then adopt η = 0.4 and η = −0.4.The positive η corresponds to a standard Drucker-Prager material with yield strength greater in compression than in tension.The negative value η = −0.4models a material with yield strength in tension than in compression.An initial uniform mesh containing 6400 linear triangles and 3321 nodes was adopted to discretize the hold-all domain.During the optimization procedure, one step of uniform mesh refinement of the hold-all domain was required in all cases to achieve convergence with a tolerance ϵ θ = 1 • .Convergence was attained in 26 iterations for the von Mises constraint case (η = 0) and 39 iterations in the other two cases (η = 0.4 and η = −0.4).The final mesh contains contains 25600 elements and 13041 nodes.The optimal topologies obtained are shown in Fig. 5.As one would expect, a symmetric structure is obtained under the von Mises constraint.The optimal domains for the other two cases are flipped images of each other and, as expected, under the conventional Drucker-Prager constraint (η = 0.4) the member under compression (on the right) is bulkier than the member under tensile stresses (on the left).

L-bracket
Now we turn our attention to a classical structural optimization problem containing a geometrical singularity -the L-bracket problem subject to stress constraints.The hold-all domain and loading are illustrated in Fig. 6.This problem has been studied by a number of authors and various strategies have been proposed for the treatment of the von Mises stress constraint, exclusively in the context of SIMP-based methods of structural optimization (refer to [21] and references therein).The solution of this optimization problem (with a slightly different loading condition to that of Fig. 6) under a von Mises constraint by a topological derivative-based approach has been recently proposed in [8].Here we show the application of the topological-derivative approach to the case of Drucker-Prager-type constraints.The lengths of the horizontal and vertical branches of the L-bracket are respectively 2 m and 2.5 m measured along their centre lines.Both have identical width of 1 m.The structure is clamped at the top edge and a point load g = −(0, 40) KN/m is applied to the corner of the right tip.The elastic properties of the structure material are E hard = 12500 MPa and ν = 0.2.The Drucker-Prager yield criterion parameters are set as η = −0.3703and σ ⋆ = 63.85MPa.These are chosen so that the Drucker-Prager yield surface matches the compressive and tensile uniaxial yield strengths [13] of a natural wood, given respectively by f c = 46.6MPa and f t = 101.4MPa.The stress constraint is not enforced in the white region of radius R = 0.15 m directly under the point of load application (shown in Fig. 6).The initial (non-uniform) mesh discretizing the hold-all domain has 14236 three-noded triangular elements and 7323 nodes with a higher density of elements around the reentrant corner that gives rise to the stress singularity.Figure 7 shows the optimum structures obtained without and with the enforcement of the Drucker-Prager stress constraint.In the stress-constrained case, the penalty coefficient adopted in the penalized objec-tive functional was α = 10 4 with p = 32.In both cases we set β = 1/3.As in the previous example, one step of uniform mesh refinement is performed in both cases to achieve convergence.The final mesh here has a total of 58240 elements and 29532 nodes.The convergence tolerance adopted for both unconstrained and stress-constrained problems is ϵ θ = 1 • with a total number of iterations required for convergence being 39 and 62, respectively.The evolution of the objective functional, volume fraction and angle θ throughout the iterations of the optimization algorithm is shown in Fig. 8(a-c).We remark here that the adopted tolerances are quite stringent and a converged design for practical purposes is in fact obtained for the constrained case with the (quite satisfactory) initial mesh at iteration 39.This is where a sharp variation in θ and I α Ω is depicted in Figs.8(a) and 8(c), corresponding to the mesh refinement step.Figure 8(d) shows the history of the worst stress ratio in the structure: whose maximum admissible value is 1.It should be observed that in the stress constrained case shown in Fig. 7(b) the reentrant corner has been rounded by the algorithm.The corresponding worst stress ratio in this case (shown in Fig. 8(d)) is 1.0184 for the converged structural domain -very close to its saturation value of 1.In the unconstrained case, on the other hand, the worst stress ratio blows up when minimizing the compliance due to the geometrical singularity of the reentrant corner.It is worth noting here that the rounding off of the reentrant corner in the stressconstrained problem has been achieved by the present algorithm in a most natural manner without any added post-processing techniques.This is a mere consequence of the use of the exact formula (33) for the topological derivative of the objective functional.This formula gives the exact sensitivity of the penalized objective functional with respect to the considered black-and-white-type topological changes.The only approximation here is the use of a penalty term to enforce the required stress constraint.In SIMP-based methodologies on the other hand, some of the exact information on the sensitivity to black-and-white-type topological changes (i.e.first order terms of the topological asymptotic expansion of the objective functional) is inevitably lost with the introduction of the regularized density field that approximates the sharp black-white transition.The enforcement of stress constraints with such methods poses a more significant challenge and requires, for instance, the use of post-processing techniques to retrieve stresses.In this context, many such procedures have been proposed and used with success in a number of stress-constrained problems (a recent overview is provided in [21]).

Bridge design
This last example considers the design of a bridge.The hold-all domain is a rectangle 180 m long and 60 m high illustrated in Fig. 9.The bridge is assumed clamped at the two bottom supports of equal length a = 9 m.A uniformly distributed traffic load g = −(0, 400) KN/m 2 is applied to the edge of the dark strip of height h = 3 m indicated in Fig. 9 that represents the road and will remain unchanged throughout the optimization process.The strip is positioned at a distance c = 27 m from the top of the hold-all domain.The material properties are E hard = 27500 MPa and ν = 0.2.For the purpose of comparison, the optimization procedure is carried for two cases: (a) No stress constraints (α = 0), and (b) The Drucker-Prager stress constraint with yield strength parameters η = 0.417 and σ ⋆ = 5.05 MPa.These parameters are obtained from the Drucker-Prager biaxial fit model [13] to match a tensile and compressive yield strength of f c = 30.5MPa and f t = 2.75 MPa respectively.For the stress-constrained case we adopt the penalty coefficient α = 10 3 and in both cases we choose β = 1/10 and the convergence tolerance ϵ θ = 1 • .The stress constraint is not enforced within the white region of size 15 × 15m adjacent to the bottom supports.Due to symmetry, only half of the hold-all domain is discretized.The initial (uniform) mesh has 4800 elements and 2501 nodes.In both cases, two steps of uniform mesh refinement are performed leading to a final mesh of 76800 and 38801 nodes.Figure 10 shows the optimized topologies obtained for the two cases.The total number of iteration required for convergence was 16 and 13, respectively, for the unconstrained and constrained cases.Note that the unconstrained optimization results in the well-known tie-arch bridge design.In this design some structural members are under tensile and others under compressive dominant stresses.The stress-constrained optimization with the Drucker-Prager criterion, on the other hand, results in a radically different design where all members are subject to compressive dominant stresses.Such designs are typical in practice for materials whose compressive strength is much higher than their tensile strength (such as concrete).Its automatic generation here clearly demonstrates the success of the proposed topology optimization procedure.

CONCLUSION
This paper has extended the result derived in [8] to incorporate the Drucker-Prager stress constraint within a topological derivative-based algorithm for topology optimization of elastic structures.In this context a penalty functional has been proposed to enforce the point-wise Drucker-Prager constraint and a closed formula for its topological derivative has been presented.The overall algorithm, which uses the topological derivative to indicate the descent direction in conjunction with a level-set representation of the structure domain, is of remarkably simple computational implementation.In particular, it does not feature postprocessing procedures (such as filtering or relaxation) of any kind and only a minimal number of user-defined algorithmic parameters are needed.This is in sharp contrast with current methodologies of topological structural optimization with local stress constraints.Numerical examples have demonstrated the effectiveness and efficiency of the algorithm in the solution of topology optimization problems under the considered class of constraints.The algorithm was shown, for instance, to efficiently handle topology optimization with materials displaying strong asymmetry in their tensile and compressive uniaxial yield strengths.From a practical standpoint, we believe this fact to be particularly relevant in that it opens the possibility for the efficient automatic design/optimization of structures made of a much wider range of materials than that for which stress-constrained topology optimization has been mainly used so far.Extensions of the present approach to the three-dimensional elasticity system is currently under investigation.

Figure 4 .
Figure 4. Wall under shear load.Initial guess and boundary conditions.

4 Figure 5 .
Figure 5. Wall under shear load.Obtained design for different values of η.

Figure 9 .
Figure 9. Bridge design.Initial guess and boundary conditions.

Figure 10 .
Figure 10.Bridge design.Obtained design for the unconstrained and constrained cases.