Florian Mitjana

Sonia Cafieri

Florian Bugarin

Stéphane Segonds

F Castanie

P Duysinx

Topological gradient in structural optimization under stress and buckling constraints

Keywords: topology optimization, topological gradient, buckling constraints, generalized eigenvalue problem

Introduction

Topology optimization plays a prominent role in structural engineering, where one seeks the optimal distribution of a given material in a given design space. Applications arise typically in industry and include aircraft [START_REF] Remouchamps | Application of a bi-level scheme including topology optimization to the design of an aircraft pylon[END_REF]; [START_REF] Zhu | Topology optimization in aircraft and aerospace structures design[END_REF]) and automotive [START_REF] Calvel | Conception d'organes automobiles par optimisation topologique[END_REF]) component design, and civil engineering [START_REF] Kingman | Applications of topology optimisation in structural engineering: high-rise buildings & steel components[END_REF]).

In the context of continuum structures, a mathematical optimization problem and an associated method to compute the optimal topology are defined according to two main approaches. The first one consists in associating a design variable to each element of the mesh defining the design domain. Design variables may be continuous, representing the value of the density associated to each element, or binary, to signify the presence/absence of each element (where the values 1 and 0 can be interpreted as extreme values for densities). Associated methods are sometimes referred to as density-based methods [START_REF] Bendsoe | Topology optimization: theory, methods, and applications[END_REF]; [START_REF] Deaton | A survey of structural and multidisciplinary continuum topology optimization: post 2000[END_REF]). The second approach relies on implicit functions to define the structural boundaries. Associated methods are Level-set methods [START_REF] Wang | A level set method for structural topology optimization[END_REF]; [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF]). In this paper, we consider an explicit parametrization of the design domain, that is discretized and meshed, and a binary variable is associated to each element. This yields to a combinatorial optimization problem. Topology optimization problems are characterized by constraints generally related to the inherent mechanical properties of the material and to the resistance of the structure under applied efforts in tension and compression. In this paper, we deal with stress and buckling constraints, and we aim at minimizing the structural mass.

Stress constraints are related to the integrity of the material under the applied external loads. Buckling constraints are related to the failure of a mechanical structure under instable phenomena due to compressive loads. From a mechanical point of view, these constraints play a crucial role in the behavior of the structure and have to be taken into account as soon as possible in the design process and so in the topological optimization [START_REF] Duysinx | Topology and generalized shape optimization: Why stress constraints are so important[END_REF]; [START_REF] Shanglong | Adaptive mesh refinement for topology optimization with discrete geometric components[END_REF]). However , these kind of constraints, while of most importance in real applications, are difficult to handle in an optimization process, especially due to their non-linear nature.

Considering stress constraints, another difficulty is related to the fact that one of such constraints is associated to each of the elements of the discretized domain, so that there are as many design variables as there are constraints.

Thus, the computational cost to handle such constraints may become critical for large-scale problems. To overcome this difficulty, aggregation techniques have been proposed, such as those based on the Kreisselmeier-Steinhauser function [START_REF] Kreisselmeier | Systematic control design by optimizing a vector performance index[END_REF]), the P-norm [START_REF] Duysinx | New developments in handling stress constraints in optimal material distribution[END_REF] or a unified aggregation and relaxation approach [START_REF] Verbart | A unified aggregation and relaxation approach for stress-constrained topology optimization[END_REF]). In this paper, we consider the formulation proposed by [START_REF] Verbart | A unified aggregation and relaxation approach for stress-constrained topology optimization[END_REF] and propose an efficient method to compute the derivatives of stress constraints.

Considering buckling phenomena in the context of topological optimization brings two kinds of difficulties: one related to the computational cost of the buckling analysis [START_REF] Bian | Large-scale buckling-constrained topology optimization based on assembly-free finite element analysis[END_REF]) and the other one to the poor convergence of some optimization processes [START_REF] Bruyneel | Discussion on some convergence problems in buckling optimisation[END_REF]; Gao et al. (2020)).

Buckling analysis of a structure is carried out through solving first a linear system to evaluate the structural displacements, and then a generalized eigenvalue problem implying two matrices. The first one is the stiffness matrix related to the structure, while the second one is the geometric stiffness matrix related to the deformations caused by the loading. The generalized eigenvalues represent the load factors of buckling, while the associated eigenvectors, related to displacements of the structure, are referred to as buckling modes. The smallest positive eigenvalue is called critical buckling load. For large scale problems, the buckling analysis can become highly memory and time consuming. Then, one challenge lies in an efficient method to solve the generalized eigenvalue problem.

Krylov sub-space methods [START_REF] Lehoucq | ARPACK users' guide: solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[END_REF]) and the Block Jacobi Conjugate Gradient method [START_REF] Ovtchinnikov | Computing several eigenpairs of hermitian problems by conjugate gradient iterations[END_REF]) lead to a very slow convergence due to the poor conditioning of the two involved matrices. Based on a shift-invert technique and a matrix factorization, [START_REF] Dunning | Levelset topology optimization with many linear buckling constraints using an efficient and robust eigensolver[END_REF] improved the rate of convergence of the Block Jacobi Conjugate Gradient method, but their aproach is highly memory consuming. The authors of [START_REF] Bian | Large-scale buckling-constrained topology optimization based on assembly-free finite element analysis[END_REF] proposed an assembly-free iterative solver to efficiently estimate the critical buckling load.

Their solver is based on the inverse iteration method, coupled with an efficient matrix-blocks computation that leads to low memory requirements.

In the context of structural optimization under buckling constraints, it is usual to compute only the critical buckling load as it is considered enough to assess the structural resistance to buckling [START_REF] Lindgaard | On compliance and buckling objective functions in topology optimization of snap-through problems[END_REF][START_REF] Luo | Structural topology optimization for maximum linear buckling loads by using a moving iso-surface threshold method[END_REF]; [START_REF] Bian | Large-scale buckling-constrained topology optimization based on assembly-free finite element analysis[END_REF]). However, as highlighted by [START_REF] Bruyneel | Discussion on some convergence problems in buckling optimisation[END_REF], as the topology of the structure changes during optimization, the buckling mode associated to the critical buckling load may also change, leading to a slow convergence or to divergence of the optimization process. Several works [START_REF] Bruyneel | Discussion on some convergence problems in buckling optimisation[END_REF]; [START_REF] Browne | A fast method for binary programming using first-order derivatives, with application to topology optimization with buckling constraints[END_REF]; [START_REF] Dunning | Levelset topology optimization with many linear buckling constraints using an efficient and robust eigensolver[END_REF]Gao et al. (2020)) showed that using a set of several eigenvalues in the optimization problem may overcome this difficulty. In this paper, building on the the inverse iteration method of [START_REF] Bian | Large-scale buckling-constrained topology optimization based on assembly-free finite element analysis[END_REF], we propose an extension of such a method that is able to efficiently compute multiple buckling modes.

Another difficulty of handling buckling constraints is related to the computation of their derivatives. By construction, the geometric stiffness matrix depends on the decision variables and on the stress of elements. In order to reduce the computational cost, the partial derivatives with respect to the stress are often neglected [START_REF] Neves | Generalized topology design of structures with a buckling load criterion[END_REF]; [START_REF] Bruyneel | Discussion on some convergence problems in buckling optimisation[END_REF]). However, this approach affects the convergence of the optimization process. Although an explicit formulation of the derivatives of the geometric stiffness matrix has been proposed by [START_REF] Browne | A fast method for binary programming using first-order derivatives, with application to topology optimization with buckling constraints[END_REF], their computation is highly computational expensive, thus untractable for large problems.

Most of topology optimization methods fall in the category of density-based methods. The most common one is the Solid Isotropic Material with Penalization (SIMP), based on a continuous relaxation of binary decision variables.

For problems under buckling constraints, the continuous relaxation implies ill-conditioning of the stiffness and the geometric stiffness matrices involved in the generalized eigenvalue problem, that may strongly impact the convergence of an iterative method (such as the inverse iteration method) for solving the generalized eigenvalue problem [START_REF] Suresh | Efficient generation of large-scale pareto-optimal topologies[END_REF]). Moreover, the continuous relaxation may introduce undesirable buckling modes, called spurious buckling modes [START_REF] Neves | Generalized topology design of structures with a buckling load criterion[END_REF]).

A recent approach, which combines the SIMP method with the Non-Uniform Rational Basis Spline (NURBS) hyper-surfaces framework [START_REF] Costa | A 2d topology optimisation algorithm in nurbs framework with geometric constraints[END_REF][START_REF] Costa | Nurbs hyper-surfaces for 3d topology optimization problems[END_REF])), allows one to describe the topology of a continuum structure through a geometrical entity. The number of design variables is reduced with respect to density-based methods, yielding an interesting advantage to address large-scale problems. This approach has recently been proved to efficiently handle constraints such as minimum and maximum length scale (Costa et al. (2019a,c)), structural displacements [START_REF] Rodriguez | Structural displacement requirement in a topology optimization algorithm based on isogeometric entities[END_REF]) or eigenfrequencies [START_REF] Costa | Eigen-frequencies and harmonic responses in topology optimisation: A cad-compatible algorithm[END_REF]); it has also been applied to multi-material topology optimization problems (Gao et al. (2020a,b)). Thanks to the geometrical properties of NURBS hyper-surfaces, it is possible to implicitly include in the problem formulation mechanical requirements such as buckling constraints. However, as this method is based on the SIMP approach, buckling analysis is still impacted by ill-conditioning of the stiffness and the geometric stiffness matrices.

Another optimization approach, that is receiving an increasing attention [START_REF] Suresh | Stress-constrained topology optimization: a topological level-set approach[END_REF]; [START_REF] Deng | Multi-constrained topology optimization via the topological sensitivity[END_REF]), relies on topological gradients. Unlike the SIMP method, this kind of approach is not based on a problem relaxation, thus overcoming the problem of ill-conditioning of matrices in the generalized eigenvalue problem.

In this paper, a new algorithm, based on topological gradients, for solving topology optimization problems under stress and buckling constraints is presented. Our contributions are also related to the way stress and buckling constraints are handled. As regards to stress constraints, the formulation proposed by [START_REF] Verbart | A unified aggregation and relaxation approach for stress-constrained topology optimization[END_REF] is considered and a method to compute constraint derivatives by introducing an adjoint state is proposed. Concerning buckling analysis, the method based on the iterative inverse [START_REF] Bian | Assembly-free buckling analysis for topology optimization[END_REF]) is extended to compute a broader spectrum of eigenvalues. From the formulation of the buckling constraint of [START_REF] Browne | A fast method for binary programming using first-order derivatives, with application to topology optimization with buckling constraints[END_REF], we propose a method based on adjoint states to compute efficiently the associated buckling derivatives.

The article is organized as follows. In Section 2, the optimization problem formulation is presented. In Section 3, an extension of the iterative inverse method is proposed to perform the buckling analysis. Section 4 presents the optimization method proposed for the considered problem. Numerical results are presented in Section 5. Finally, Section 6 concludes the article.

Optimization problem modeling

In this section, we present the formulation of the topology optimization problem for continuous structures addressed in this paper. The design space, i.e., the design domain of the structure, is discretized in a regular grid of N cube elements in a three-dimensional space. Vertices of the cube elements are called nodes, and the total number of nodes is denoted by M . A binary variable x i is associated to each element i of the mesh meaning its presence or absence.

These decision variables will also be referred to as design variables.

We aim at optimizing the structural design in such a way that the mass of the structure is minimized. The structural mass is defined by the material density ρ multiplied by the sum of the design variables, each multiplied by the volume ω of the corresponding element. As in our model all the elements have equal volume, we can define the objective function, to be minimized, as follows:

ρω N i=1 x i (1)
As for the constraints, let us formulate the stress and the buckling constraints that characterize the problem under investigation. The evaluation of the structural stress and buckling first requires the solution of the following linear system:

Ku = f (2)
where K ∈ R 3M ×3M is the stiffness matrix, u ∈ R 3M is the displacement vector of the M nodes (each node has 3 possible displacements along the axes of the space, noted as x, y and z), f ∈ R 3M represents the loads applied on the structure. Note that one actually solves a system whose size is lower than 3M × 3M .

It is indeed easy to see that, for example, in order to keep fixed a portion of the structure, the associated nodal displacements have to be given a value of 0. More generally, boundary conditions on some nodal displacements must be imposed, thus reducing the size of (2). Once the vector u has been computed, the stress components {σ i x , σ i y , σ i z , σ i xy , σ i xz , σ i yz } of a structural element i can be evaluated in its midpoint as follows (the explicit dependence of each component from x and u being omitted for the sake of readability):

[σ i x , σ i y , σ i z , σ i xy , σ i xz , σ i yz ] = HB i u i (3)
where u i corresponds to the displacements of the nodes of the element i, H 140 is the Hooke's matrix computed from the problem data (material property) and B i is a matrix related to the geometry of the element i resulting from the discretization [START_REF] Cook | Concepts and applications of finite element analysis volume 4[END_REF]).

From equation (3), the Von Mises stress σ V M i of an element i can then be expressed, as a function of x = {x 1 , . . . , x N } and u, as follows:

σ V M i (x, u) = 1 √ 2 (σ i x -σ i y ) 2 + (σ i x -σ i z ) 2 + (σ i y -σ i z ) 2 + 6((σ i xy ) 2 + (σ i xz ) 2 + (σ i yz ) 2 ) (4)
Following [START_REF] Verbart | A unified aggregation and relaxation approach for stress-constrained topology optimization[END_REF], we finally formulate the stress constraint as:

G σ (x, u) := 1 P ln 1 N N i=1 e P gi(x,u) ≤ 0, where g i (x, u) = x i σ V M i (x, u) σ -1 ( 5
)
where σ is a constant inherent to the considered material, referred to as limit stress, and P represents an aggregation parameter [START_REF] Verbart | A unified aggregation and relaxation approach for stress-constrained topology optimization[END_REF]). When 145 P goes to infinity, this function converges to the maximal value of g i .

From the stress components in (3), the geometric stiffness matrix K σ ∈ R n×n can be assembled taking into account the blocks corresponding to each element i, through the following expression:

K σ = N i=1      Ωi G      σ i 0 0 0 σ i 0 0 0 σ i      G dxdydz      where σ i (x, u) =      σ i x σ i xy σ i xz σ i xy σ i y σ i yz σ i xy σ i yz σ i z      (6)
where Ω i represents the volume of element i and G is a matrix whose components can be computed from the discretization of the meshed domain [START_REF] Cook | Concepts and applications of finite element analysis volume 4[END_REF]), and [ ] represents the block assembly process. Thus, we are able to define the generalized eigenvalues problem related to buckling:

(K + λK σ )v = 0 (7)
where λ are scaling factors (critical loads) of the applied forces f, and the associated eigenvectors v represent the corresponding deformations of the structure (buckling modes).

To ensure that the structure resists to buckling, one needs to ensure that the critical buckling load λ cr is larger than a safety coefficient c s . Although using the critical buckling load λ cr is enough to assess the structural stability, several works [START_REF] Bruyneel | Discussion on some convergence problems in buckling optimisation[END_REF]; [START_REF] Browne | A fast method for binary programming using first-order derivatives, with application to topology optimization with buckling constraints[END_REF]; [START_REF] Dunning | Levelset topology optimization with many linear buckling constraints using an efficient and robust eigensolver[END_REF]; Gao et al. (2020) showed that using a set of several buckling modes allows one to avoid oscillations and slow convergence of the optimization process.

In [START_REF] Browne | A fast method for binary programming using first-order derivatives, with application to topology optimization with buckling constraints[END_REF], the authors use the positive semi-definite property of the stiffness matrix to formulate a unique buckling constraint which takes into account several buckling modes. Following [START_REF] Browne | A fast method for binary programming using first-order derivatives, with application to topology optimization with buckling constraints[END_REF], we formulate the buckling constraint as follows:

G λ := 1 m m j=1 v j (K + c s K σ )v j ≥ 0 (8)
where the vectors v j are the m buckling modes of the generalized eigenvalue 150 problem (7) and c s is the safety coefficient.

The problem of mass minimization under stress and buckling constraints can finally be summarized as the following nonlinear combinatorial optimization problem:

min x=(x1,...,x N ) ρω N i=1 x i s.t. 1 P ln 1 N N i=1 e P.gi(x,u) ≤ 0, where g i (x, u) = x i σ V M i (x, u) σ -1 1 m m j=1 v j (K + c s K σ )v j ≥ 0 Ku = f (K + λK σ )v = 0 x i ∈ {0, 1} ∀i = 1, . . . , N. (P)

Linear buckling analysis

In this section, we propose an extension of the inverse iteration method developed by [START_REF] Bian | Large-scale buckling-constrained topology optimization based on assembly-free finite element analysis[END_REF] to compute the m buckling eigenvectors associated to equation ( 8) (defining buckling constraints of the optimization problem (P)).

Let us first briefly recall the main elements of the method. For further details, please refer to [START_REF] Bian | Large-scale buckling-constrained topology optimization based on assembly-free finite element analysis[END_REF]. The method aims at avoiding assembly of the stiffness and of the geometric stiffness matrices, K and K σ , in order to reduce memory cost. The main operation, relying on a matrix-vector multiplication, is performed as a matrix-blocks computation [START_REF] Hughes | An element-by-element solution algorithm for problems of structural and solid mechanics[END_REF]. As an example, considering the stiffness matrix K and the structural displacement vector u, the matrix-vector multiplication Ku is performed as:

Ku = N i=1 K i u i (9)
where K i is the elementary stiffness matrix of element i, u i are the displacements of the nodes of element i and [ ] represents the block assembly process of the 155 K i u i vectors in the Ku vector. In the case of a regular mesh grid, meaning that all finite elements are identical (as it is in our case), the matrices K i are identical for all i. Thus, for a stiffness matrix-vector multiplication, only the stiffness matrix K i of an element needs to be stored, leading to low memory requirements [START_REF] Yadav | Large scale finite element analysis via assemblyfree deflated conjugate gradient[END_REF]).

160

The computational cost is also reduced through the application of the deflated conjugate gradient method [START_REF] Saad | A deflated version of the conjugate gradient algorithm[END_REF]) that appears to be particularly efficient to solve the (sparse) linear systems in the considered context.

More precisely, [START_REF] Bian | Large-scale buckling-constrained topology optimization based on assembly-free finite element analysis[END_REF] used the deflated conjugate gradient method within an inverse iteration method [START_REF] Ipsen | Computing an eigenvector with inverse iteration[END_REF]; [START_REF] Golub | Inexact inverse iteration for generalized eigenvalue problems[END_REF]), that they adapted to the generalized eigenvalue problem (7) to determine the critical buckling load λ cr and its associated eigenvector v cr . The method consists in solving iteratively the following linear system starting from an initial approximation v 0 :

v k+1 = -K -1 K σ v k (10)
where the vector v k is normalized at each iteration. Once the sequence v k has converged to an eigenvector v, we can compute the associated eigenvalue λ by the Rayleigh quotient functional [START_REF] Timoshenko | Theory of elastic stability[END_REF]) related to the generalized eigenvalue problem (7):

R(v) = - v Kv v K σ v (11)
It is important to note that the eigenvalue λ corresponds to the eigenvalue with the smallest modulus of the generalized eigenvalue problem ( 7). If we assume that the geometric stiffness matrix is positive definite, then the eigenvalue λ always corresponds to the critical buckling load λ cr . However, as it is not necessarily the case, the eigenvalue λ might correspond either to the smallest 165 positive eigenvalue or to the smallest negative eigenvalue. Thus, in the context of structural optimization, this method can be applied only when the smallest eigenvalue has a positive modulus.

We propose to extend the inverse iteration method to compute a set of m positive eigenvalues of the generalized eigenvalue problem (7). The method is based on the minimization of the Rayleigh quotient functional of a set of vectors. First note that the minimization of the Rayleigh quotient of m vectors will converge to the m associated eigenvalues of smallest moduli. However, as noted above, the geometric stiffness matrix is not necessarily positive definite, so the associated eigenvalues are not always positive. To overcome this problem, let us consider the set Λ of the m > m eigenvalues of smallest moduli:

Λ := {λ j , ∀j = 1, . . . , m | |λ 1 | ≤ |λ 2 | ≤ • • • ≤ |λ m-1 | ≤ |λ m|}
Assuming that m positive eigenvalues exist for the considered problem (as is the case with most real-world applications in the field of structural optimization), and assuming that the cardinality of set Λ is large enough, we can extract the subset of m smallest positive eigenvalues. Thus, in order to determine the m eigenvectors associated with the m smallest positive eigenvalues, we propose to apply the inverse iteration method on a set V of m vectors of size 3M :

V ∈ R 3M × m := {v 1 , . . . , v m }
When solving the equation ( 10) for each vector of V, it is necessary to maintain the K-orthogonality of the vectors [START_REF] Dunning | Levelset topology optimization with many linear buckling constraints using an efficient and robust eigensolver[END_REF]), otherwise each component of V will converge to the eigenvector associated to the eigenvalue of smallest modulus. In order to maintain the K-orthogonality between the vectors belonging to V, we propose to use the Rayleigh-Ritz procedure [START_REF] Trefethen | Numerical linear algebra[END_REF]), consisting in solving the generalized eigenvalue problem, of dimension m, as follows:

V K σ Vw = λV KVw (12)
where the vectors v = Vw are called Ritz vectors (where w are the eigenvectors of ( 12)) and the eigenvalues λ are called Ritz values. Note that the Ritz vectors are K-orthogonal by construction. Moreover the system ( 12) is much simpler to solve than the initial problem (7) since m << 3M . This kind of problems can be solved by standard numerical linear algebra solvers as those implemented in the LAPACK library [START_REF] Anderson | LAPACK Users' guide[END_REF]).

The procedure that we propose to compute the set of m eigenvalues and eigenvectors of system ( 7) is summarized in Algorithm 1. Note that the Algorithm 1 stops when the m smallest positive eigenvalues have been determined.

Algorithm 1 Generalized inverse iteration method

Input: Number m of seeked eigenvectors, dimension m of the research space, and a tolerance

1: Initialize V ∈ R 3M × m = {v 1 , . . . , v m }, p = 1, l = 1 2: Solve V K σ Vw = λV KVw 3: V ← Vw 4: while l < m do 5: Compute z j = K σ v j , ∀j = 1, . . . , m 6:
Solve Ky j = z j , ∀j = 1, . . . , m 7: end while 17: end while Output: {λ j , v j | λ j ≥ 0}

Solve Y K σ Yw = λY KYw , where Y ∈ R 3M × m = {y 1 , . . . , y m } 8: V ← Yw 9: Compute λ j = - v j Kv j v j K σ v j , ∀j = 1, . . . , m 10: Compute r j = Kv j + λ j K σ v j 2 , ∀j = 1, . . . ,

Topological gradient-based optimization

In this section, we propose an algorithm for topology optimization of a problem under stress and buckling constraints, based on the concept of topological 180 gradient. After a brief recall of such a concept and of its use within an optimization algorithm, we focus the stress and buckling constraints of the optimization problem formulated in Section 2. In an unified framework, we present a computation of topological gradients for the stress constraints, and of gradients for the buckling constraints. We finally propose our algorithm.

T C (p i ) = lim r→0 C(Ω B ) -C(Ω) h(r) (13)
where h(r) is a positive scalar function such that lim r→0 h(r) = 0.

Topological gradients can then be used within topology optimization algorithms to define the principle of iterative material removal. Starting from a design domain that is initially considered full of material, and associated to a volume α equal to 1, a straightforward approach consists in fixing a threshold τ and removing the fraction of volume ∆α corresponding to the elements for which the topological gradient of the optimization criterion is less than τ . Thus, the quantity of removed material can be a different one at each iteration, as it depends on the (fixed) value of τ . Despite its simplicity, this kind of approach is not generally used in structural optimization, as it leads to a poor convergence of the optimization process [START_REF] Suresh | A 199-line matlab code for pareto-optimal tracing in topology optimization[END_REF]).

Another approach relies on fixing the fraction of volume ∆α and determining the threshold τ such that the amount of volume corresponding to the elements whose topological gradient is less than τ is equal to ∆α [START_REF] Deng | Multi-constrained topology optimization via the topological sensitivity[END_REF]).

Thus, unlike the previous approach, the value of τ is adapted throughout the optimization process. More specifically, given a domain Ω whose volume is α, and given a fraction ∆α of volume to remove, one has to determine a scalar τ , that can be assimilated to a cut-off value, in order to obtain a perturbed domain Ω τ whose volume is α -∆α. The value of τ can be determined from a topological gradient field T by identifying the elements whose associated topological gradient is larger than τ , while ensuring that the volume is equal to α -∆α.

Within an optimization algorithm, the value of the fraction of volume to remove is usually iteratively decreased, and iterations are run until a termination criterion is satisfied. The feasibility of the current topology for the addressed structure, i.e., the satisfaction of the constraints of the optimization problem, is checked at each iteration. This guarantees that a feasible, local albeit not global optimal solution, is always computed. Note that the value of the fraction of volume ∆α plays a crucial role, since it leads the entire optimization process.

Indeed, if its initial value is very small (for example, less than 1% of the initial volume) then the convergence will be slow, while if its value is large (for example, more than 30% of the volume) then the optimization may diverge [START_REF] Suresh | Efficient generation of large-scale pareto-optimal topologies[END_REF]).

An algorithmic scheme for structural optimization based on topological gradients is displayed in Figure 1. We compute a topological gradient, T σ , associated to the constraint G σ , and a gradient, T λ , associated to G λ . They both will be used to estimate the influence of a perturbation introduced on the structure whose topology is being optimized. The computation of T σ and T λ represents one of the contributions of this paper. Topological gradients have indeed been introduced in the literature for compliance constraints [START_REF] Feijoo | The topologicalshape sensitivity method in two-dimensional linear elasticity topology design[END_REF]) and for the P-norm of stress constraints [START_REF] Suresh | Stress-constrained topology optimization: a topological level-set approach[END_REF]), while constraints like G σ and G λ

have not been previously addressed in this context. Our approach to determine T σ and T λ is based on the adjoint state method [START_REF] Choi | Structural sensitivity analysis and optimization 1: linear systems[END_REF]; [START_REF] Deng | Topology optimization under thermo-elastic buckling[END_REF]; Allaire (2015)), as described below.

Stress constraints

Let us consider the stress constraint G σ . As it depends on the decision variables and on the structural displacements u, its derivative according to the decision variable x i is expressed as:

dG σ dx i = ∂G σ ∂x i + ∂G σ ∂u du dx i (14)
The structural displacements are obtained throught the equation Ku = f (see ( 2)). Thus by deriving according to x i , we obtain:

K du dx i + ∂K ∂x i u = 0 ⇔ du dx i = -K -1 ∂K ∂x i u (15)
Then, the derivative of G σ according to x i can be formulated as:

dG σ dx i = ∂G σ ∂x i - ∂G σ ∂u K -1 ∂K ∂x i u (16)
It is important to note that, to solve equation ( 16), we need to solve a linear system involving the stiffness matrix K. Then, computing the derivatives of G σ according to the decision variables x i , ∀i = 1, . . . , N , requires the solution of N linear systems, implying a potentially high computational cost. An alternative is represented by the adjoint method, that relies on the fact that the term -∂Gσ ∂u K -1 does not depend on the decision variables and remains constant regardless of the variable x i in (16). The adjoint method consists in defining an adjoint variable µ [START_REF] Choi | Structural sensitivity analysis and optimization 1: linear systems[END_REF]), solution of the following adjoint state:

Kµ = - ∂G σ ∂u ( 17
)
where K is the stiffness matrix, µ ∈ R 3M is the adjoint state and u corresponds to a vector of displacements. Then, the computation of the derivatives of G σ 250 according to x i , ∀i = 1, . . . , N involves only one additional linear system.

From the formulation of G σ in the equation ( 5), the term -∂Gσ ∂u can be expressed as:

- ∂G σ ∂u = - 1 P 1 1 N N i=1 e P.gi(x,u) 1 N N i=1 P x i σ e P.gi(x,u) ∂σ V M i ∂u (18)
where [ ] corresponds to the block assembly operator. From (4), the derivative of σ V M i with respect to a displacement u can be expressed as:

∂σ V M i ∂u = 1 σ V M i √ 2 (σ i x -σ i y )( ∂σ i x ∂u - ∂σ i y ∂u ) + (σ i x -σ i z )( ∂σ i x ∂u - ∂σ i z ∂u ) + (σ i y -σ i z )( ∂σ i y ∂u - ∂σ i z ∂u ) + 6(σ i xy ∂σ i xy ∂u + σ i xz ∂σ i xz ∂u + σ i yz ∂σ i yz ∂u ) (19)
Note that the derivatives of the stress components σ i x , σ i y , σ i z , σ i xy , σ i xz and σ i yz with respect to a displacement u are constant and correspond to :

∂σ i x ∂u = F 1,: ∂σ i y ∂u = F 2,: ∂σ i z ∂u = F 3,: ∂σ i xy ∂u = F 4,: ∂σ i xz ∂u = F 5,: ∂σ i yz ∂u = F 6,: (20)
where F = HB ∈ R 6×24 , with H and B matrices introduced in (3). In the case of a cube element (as those considered in our domain discretization), B i ∈ R 6×24

and the terms F i,: ∀i = 1, . . . , 6 correspond to to the rows of the F matrix.

Thus, the term -∂Gσ ∂u can be expressed as:

-∂Gσ ∂u = - 1 N i=1 e P.gi(x,u) N i=1            P x i σ e P.gi(x,u) 1 σ V M i √ 2            (σ i x -σ i y )(F 1,: -F 2,: )+ (σ i x -σ i z )(F 1,: -F 3,: )+ (σ i y -σ i z )(F 2,: -F 3,: )+ 6σ i xy F 4,: + 6σ i xz F 5,: + 6σ i yz F 6,:                       (21)
Then, the adjoint state µ can be determined by solving the adjoint system (17) by the deflated conjugate gradient method. Thus, the topological gradient T σ of an element x i is expressed as [START_REF] Feijoo | The topologicalshape sensitivity method in two-dimensional linear elasticity topology design[END_REF]):

T σ (x i ) = -
where ν is the Poisson ratio, σ(u(x i )) is the stress matrix, as noted in (6), according to the vector of displacements u of the element i and (µ(x i )) is the deformation matrix of the µ adjoint vector of the i element. The operator :

corresponds to the sum of the matrix components obtained from the Hadamard product. The tr operator is the trace operator of a matrix.

Buckling constraints

Let us now consider the buckling constraint G λ . First, observe that the derivative of G λ with respect to x i can be expressed as:

dG λ dx i = 1 m m i=1 2v i (K + c s K σ ) dv i dx i + v i ( dK dx i + c s dK σ dx i )v i (23)
The main difficulty in computing ( 23) lies in the term dK σ dx i . Indeed, the K σ matrix depends on the design variables x i and on the stress σ i of each structural element i, that is also depending on the design variables. A computation based on finite differences, as well as a direct computation using

dK σ dx i = ∂K σ ∂x i + ∂K σ ∂ σ d σ dx i , (24)
would be computationally expensive for large problems (the term ∂K σ ∂ σ d σ dx i requiring the solution of N linear systems involving the stiffness matrices). [START_REF] Browne | A fast method for binary programming using first-order derivatives, with application to topology optimization with buckling constraints[END_REF] proposed a direct method to compute (24), however this method is highly computational expensive, thus not tractable for large problems.

Then, we apply again the adjoint method.

Following the work of [START_REF] Deng | Topology optimization under thermo-elastic buckling[END_REF] on the derivative of one buckling eigenvalue, we propose to introduce in the expression of G λ the adjoint variables

γ i ∈ R ∀i = 1, . . . , m, β ∈ R 6N and w ∈ R 3M
, and consider the following expression:

G λ := 1 m m i=1 v i (K + c s K σ )v i + γ i (1 -v i Kv i ) + β ( σ -Yu) + w (f -Ku) (25)
where Y ∈ R 6N ×3M := N i=1 [(HB) i ] corresponds to the block assembly of the product of the Hooke's matrix H and the finite element deformations matrix B defined for each element i, and σ : Differentiating (25) with respect to x i gives the following expression:

= {σ 1 x , σ 1 y , σ 1 z , σ 1 xy , σ 1 xz , σ 1 yz , . . . , σ N x , σ N y , σ N z , σ N xy , σ N xz , σ N yz }
dG λ dx i = 1 m m i=1 2v i (K + c s K σ ) dv i dx i + v i ( dK dx i + c s ∂K σ ∂x i + c s ∂K σ ∂ σ d σ dx i )v i -γ i (2v i K dv i dx i + v i dK dx i v i ) + β ( d σ dx i - dY dx i u -Y du dx i )- w ( dK dx i u + K du dx i ) (26)
Eliminating the terms dv i dx i from equation ( 26) allows one to compute the adjoint variables γ i ∈ R∀i = 1, . . . , m:

dv i dx i (K + c s K σ )v i = γ i dv i dx i Kv i , ∀i = 1, . . . , m ⇔ dv i dx i ((1 -γ i )K + c s K σ )v i = 0, ∀i = 1, . . . , m (27)
Recalling that (K + λ i K σ )v i = 0, ∀i = 1, . . . , m, the adjoint variables γ i , ∀i = 1, . . . , m are chosen as:

γ i = 1 - c s λ i , ∀i = 1, . . . , m (28)
so that equation ( 26) becomes:

dG λ dx i = 1 m m i=1 v i ( dK dx i + c s ∂K σ ∂x i + c s ∂K σ ∂ σ d σ dx i )v i -γ i v i dK dx i v i ) + β ( d σ dx i - dY dx i u -Y du dx i ) -w ( dK dx i u + K du dx i ) (29)
The adjoint variable β is chosen in such a way that the term d σ dx i is eliminated from the equation (29) :

c s m m i=1 v i ( ∂K σ ∂ σ d σ dx i )v i + β d σ dx i = 0 (30)
After factorization and rearrangement, the adjoint variable β corresponds to [START_REF] Deng | Topology optimization under thermo-elastic buckling[END_REF]; [START_REF] Bian | Large-scale buckling-constrained topology optimization based on assembly-free finite element analysis[END_REF]):

β = - c s m m i=1 v i ∂K σ ∂ σ v i ( 31
)
where ∂K σ ∂ σ corresponds to:

∂K σ ∂ σ = N i=1 6 k=1 ∂K σ ∂ σik (32)
Equation ( 29) becomes:

dG λ dx i = 1 m m i=1 v i ( dK dx i + c s ∂K σ ∂x i )v i -γ i v i dK dx i v i ) -β ( dY dx i u + Y du dx i ) -w ( dK dx i u + K du dx i ) (33)
The last adjoint variable w is chosen such that the term du dx i is removed from the equation ( 33) :

β Y du dx i + w K du dx i = 0 (34)
It is then the solution of the following system:

Kw = -Y β (35)
Finally, dG λ dx i becomes:

dG λ dx i = 1 m m i=1 v i ( dK dx i +c s ∂K σ ∂x i )v i -γ i v i dK dx i v i ) -β dY dx i u-w dK dx i u (36)
Thus, by solving the single additional linear system (35), the derivative 275 (36) with respect to each design variable x i can be computed. This derivative corresponds to the expected variation of the value of the buckling constraint G λ in the case of removal of element x i from the structure. Recall that the topological gradient allows one to quantify the influence, on a given criterion, of a perturbation of the structure. [START_REF] Mirzendehdel | A pareto-optimal approach to multimaterial topology optimization[END_REF] observe that both the gradient and the topological-gradient fields provide equivalent information up to a scale factor, in the case of the compliance. We assume that such a relation also holds for buckling constraints. Thus, for the sake of simplicity of calculation, we consider equation ( 36) for the buckling constraint G λ with respect to x i , and we use the notation T λ (x i ) to refer to such derivative, similarly to the notation used for the stress constraint.

Algorithm

We present in the following our algorithm for structural topological optimization. It follows the general idea of topological gradient-based algorithms, that is to iteratively remove a fraction of the volume of the design domain on the basis of the values of topological gradients (see subsection 4.1). The specific character of the algorithm is in the way the gradients of stress and buckling constraints are computed and combined, as well as in the criteria considered to iteratively update the topology of the structure. In order to build a single topological gradient field T w , we follow the approach consisting in defining a weighted sum of the considered fields [START_REF] Suresh | An adaptive weighting strategy for multi-load topology optimization[END_REF]). For each element x i , we propose to combine the values of T σ and T λ through the following weighted sum:

T w (x i ) = T σ (x i ) + c s λ cr 2 T λ (x i ) (37)
where λ cr is the critical buckling load. The factor (c s /λ cr ) 2 allows us to control the weight of T λ . If the value of c s /λ cr is close to one, then the optimization is driven by both stress and buckling topological gradients; if λ 2 cr is much larger than c s , i.e., the structure features a high resistance to buckling, only T σ controls the optimization. Along the iterative process, the topology of the structure is updated following the approach presented in subsection 4.1, where the considered topological gradient is T w .

The topology is updated at each iteration in such a way that the problem constraints are satisfied. This ensures that a feasible solution is always computed. One can however observe [START_REF] Deng | Multi-constrained topology optimization via the topological sensitivity[END_REF]) that removing fractions ∆α of volume may lead to a divergence of the optimization process.

Hence, we propose a further condition to be satisfied at each iteration in order to update the topology by removing ∆α. This condition is based on the (relative) variation of the value of the stress constraint G σ between the one associated to the current structure, Ω, and the one associated to the perturbed one, Ω τ :

|G σ (Ω) -G σ (Ω τ )| |G σ (Ω)| ≤ ∆G σ (38)
where ∆G σ is an user-defined threshold. The smaller this threshold is, the smaller is the variation of the stress constraint that is accepted to remove a fraction ∆α.

Note that when an element is removed, it is no longer considered in the finite element analysis. This implies that, at each iteration, the size of the linear systems and of the generalized eigenvalues problem decrease, as well as the computing time for their solution. Moreover, no penalisation schemes are required for the stiffness matrix nor for the geometric stiffness matrix. At each iteration of the optimization process, some of the removed elements can be reintroduced to allow the process to escape local optima. This corresponds, from a mechanical point of view, to strengthen some parts of the current structure.

However, for a deleted element (i.e., x i = 0), the associated values of T σ and T λ are no longer available. Thus, we base their computation on the extrapolation of the values T σ and T λ associated to the elements that are actually present in the current structure. More precisely, let us consider the topological gradient T σ (a similar procedure applies to T λ ). A value of T σ is considered for each element of the meshed domain, assuming that such a value is zero for the elements which have been removed. Following [START_REF] Krishnakumar | Hinge-free compliant mechanism design via the topological level-set[END_REF]), the extrapolation process consists first in evaluating, for each node of the mesh, the mean of the topological gradients T σ of the elements connected to such a node, and assigning this value to the node. Then, for each element of the meshed domain, the value of T σ is recomputed by averaging the values assigned to its neighboring nodes. Note that this procedure leads to a smoothing of the topological gradients. This smoothing is usually referred to as a filtering scheme [START_REF] Bendsoe | Topology optimization: theory, methods, and applications[END_REF]) and allows one to prevent the checkerboard effect [START_REF] Sigmund | Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, meshdependencies and local minima[END_REF]).

320

The proposed algorithm is presented in Algorithm 2. Note that the procedure to compute multiple buckling modes, presented in Sect.III, is nested in Algorithm 2 (steps 3 and 7).

Algorithm 2 Topological optimization under stress and buckling constraints Input: ∆α, ∆α min , ∆G σ , σ

1: Initialisation : α = 1, Ω := {x i = 1, ∀i = 1, . . . , N } 2: Ω * ← Ω 3: Solve Ku = f and (K + λK σ )v = 0 for Ω 4: Compute T σ (x i ), T λ (x i ) et T w (x i ) ∀i = 1, . . . , N 5: while ∆α ≥ ∆α min do 6:
Determine τ such that volume(Ω τ ) = α -∆α 7:

Solve Ku = f and (K + λK σ )v = 0 for Ω τ 8:

if |G σ (Ω * ) -G σ (Ω τ )| |G σ (Ω * )| ≤ ∆G σ & (c s ≤ λ cr ) & σ V M i ≤ σ ; ∀i = 1, . . . , N then 9: α ← α -∆α 10: Compute T σ (x i ), T λ (x i ) et T w (x i ) ∀i = 1, . . . , N 11: Ω * ← Ω τ 12: else 13: ∆α ← ∆α 2 14:
end if 15: end while Output: Ω *

Numerical results

In this section, we first analyze the performance of Algorithm1 to compute a set of eigenvalues. Then, numerical results for topological optimization applications, obtained performing Algorithm 2, are illustrated and discussed.

In all numerical experiments in this section, the material considered is an aluminium with a Young modulus equals to 72 000 MPa and a Poisson ratio equals to 0.33. The density of aluminum is 0.00285 g.mm -3 . The value of the limit stress σ is 510 MPa. The safety coefficient c s is fixed to 1. Recall that we consider three-dimensional cubic elements to discretize the design domain of the considered structures.

From numerical experiments, we set the parameter P to 4 in the stress constraint (5), and the ∆G σ stress threshold is set to 0.01. The material fraction ∆α is initialized to 0.1. The threshold ∆α min is fixed at a value of 1 divided by the number of elements resulting from the discretization of the design domain.

This ratio corresponds to the volume occupied by an element in the initial structure. Thus, the solution algorithm stops when it is no longer possible to remove an element from the structure.

The proposed algorithms have been implemented in Fortran and parallelized by using OpenMP commands. In Algorithm 1 the tolerance is set to 10 -4 .

The codes were run on a 32-bit Linux OS computer with an Intel(R) Core(TM) i7-4710MQ CPU @ 2.50GHz and 16GB RAM.

Buckling analysis results

Let us first consider the computation of a set of buckling modes on the L-structure illustrated in Figure 2, through Algorithm 1. The hatched area corresponds to the foundation, i.e. clampled displacements. A load F is uniformly distributed over the red area in Figure 2 and is equal to 35,000 N. The structure is discretized with cubes of 1 mm in length, leading to 331 776 finite elements and 364 429 nodes.

We consider the computation of the first 3 buckling eigenvalues (i.e., m = 3 in Algorithm 1) associated to the considered structure. On the one hand, we focus on the importance of considering the dimension m of the set V large enough to compute the first 3 buckling eigenvalues. On the other hand, we 355 investigate the impact of the dimension m of V on the computational cost. Let us recall that, through Algorithm1, the eigenvalues are determined in ascending order according to their modulus. Thus, the first three eigenvalues obtained may not correspond to positive eigenvalues. In Table 1, we report the eigenvalues obtained for values of m varying from three to five. For the structure under study, we note that with m equals to three or four, we are not able to compute the first 3 positive eigenvalues. In the case of m equals to five, we observe that the three smallest positive eigenvalues correspond to the second, third and fifth eigenvalues of smallest modulus respectively. Thus, the dimension m of V must be greater than or equal to five in order to compute the first three positive eigenvalues. It is important to note that, along a process of structural optimization, the structural topology changes at each iteration.

Thus, the buckling eigenvalues may switch their positions along the process.

Therefore, the dimension m of V must be large enough to ensure that all the eigenvalues to be searched for are obtained at each iteration of the optimization process. However, if the dimension m of V is too large, then the buckling analysis may be highly time consuming.

We now investigate the impact of the dimension m of V on the computational cost. For different values of m, we report in Table 2 the computing time, the number of iterations, and the number of Kv and K σ v products evaluated to converge to the seeked buckling eigenvalues. We note that the product K σ v is on average 7 times more expensive than the product Kv. The operation K σ v leads moreover to a larger memory footprint than Kv, since each geometric stiffness matrix element depends on its own stress matrix (6) (while K is constituted by equal blocks). Through the proposed method, the number of K σ v operations is kept low, thus yielding to important savings in time and memory consumption. This crucial point about memory consumption already led the authors in [START_REF] Yadav | Large scale finite element analysis via assemblyfree deflated conjugate gradient[END_REF] to develop an alternative approach to what is commonly done in standard commercial software (such as ANSYS - [START_REF] Stolarski | Engineering analysis with ANSYS software[END_REF], where sparse direct solvers may require huge memory consumption. In such a case, the authors recall that, as an example, for a matrix with a size of 1 million, 1GB of memory is needed to store the stiffness matrix and 10 to 20 GB additional memory for the factorization.

The method that we propose follows a similar approach than the one in [START_REF] Yadav | Large scale finite element analysis via assemblyfree deflated conjugate gradient[END_REF] for efficient memory handling, with the additional advantage to be able to compute multiple eigenvalues.

Comparing the results obtained for m equals to 5 and 6, we observe a reduction of the computing time by a factor of 1.56, although the number of Kv products remains stable. The number of iterations and the number of K σ v products decrease when the dimension of the search space is increased. Nevertheless, for larger values of m, computing time increases due to a larger number of Kv products, although the number of iterations to converge decreases.

We conclude that the choice of the value of m of V in Algorithm 1 should try to meet the right trade-off between computing time saving and the need to capture every proper buckling value required throughout the optimization process.

Topological optimization results

In this section, we present numerical results to validate our topology optimization approach and to illustrate its efficiency on real mechanical structures.

We consider two different mechanical structures: a beam under a compressive load, and a real application arising in the aeronautic field. The first application additionally allows us to highlight the interest of considering multiple buckling eigenvectors in the optimization process, while the second shows the importance of buckling constraints on the optimization process.

Beam under a compressive load

The first numerical example of topological optimization is a structure under a compressive load illustrated in Figure 3. The dimensions of the structure are l, equals to 240mm, h, equals to 42mm and e, equals to 6mm. The load F is 18,360 N and is applied uniformly over the red area in Figure 3. The hatched area corresponds to the fixed area. The structure is discretized with cubes with 415 a length of 2mm, leading to 52 920 elements. We report the results obtained when solving problem (P) by applying the proposed topological optimization algorithm. Table 3 gives, in terms of the considered number of buckling modes, the optimal mass (i.e., the value of the objective function of (P)), the final critical buckling load, the final maximum Von Mises stresses max Figures 4 and 5 illustrate the evolution, with respect to the number of computed buckling eigenvectors, of the structural mass and the critical buckling load respectively. As regards the structural mass, when the optimization is run taking into account only the first buckling eigenvector, the optimized structure has a mass of 178.50 g (Figure 6(a)). A reduction of 17% is observed when the first four buckling modes are considered (Figure 7(a)). From Figure 5, the evolution of the critical buckling load when only one buckling mode is considered, appears erratic. This is characteristic of buckling mode-switching behaviour.

A smoother evolution of the critical buckling load is observed when multiple buckling modes are considered. This is related to a reduced mode-switching behaviour along the optimization process.

Whatever the number of buckling modes considered, the proposed algorithm always converges towards an optimized structure. In the specific case of a beam in compression, the optimized structures, and therefore their masses, change according to the number of buckling modes considered in the constraint (8). From a mechanical point of view, this influence of buckling modes can be explained by comparing the two obtained structures (Figure 6(b) and Figure 7(b)). Indeed, in the case of a beam in compression, its resistance to buckling is characterized by its moment of inertia with respect to the y and z axes. For the case of four buckling modes, we observe that the material is mostly placed in the corners of the design space and the legs are connected each other by stiffeners in order to counter the effect of the first buckling mode, generating a better buckling resistance compared to the structure obtained with only one buckling eigenvec-tor. Note that this solution is very interesting from a point of view of structural engineering, that would add these stiffeners in a post-processing phase if they were not automatically generated by the optimization algorithm. In addition, when four buckling modes are considered, the critical buckling load decreases less rapidly than in the case of one buckling mode, allowing a better layout of the material within the structure and generating a lighter structure.

Gooseneck -aeronautical structure

The second numerical example is an aeronautical application (illustrated in In this example, numerical experiments show that considering more than one buckling mode does not modify the final optimal solution. This can also be explained from a mechanical point of view, considering the dimensions of the structure and the direction of the applied loads. Thus, taking into account that the optimization process for this structure may be associated to high computing times, we limit the optimization to one buckling mode. We compare the optimal structure that is obtained with and without buckling constraint, to highlight the important role of the buckling constraint on the the optimal structural design.

Table 4 reports the optimal structural mass, the final critical buckling loads, and the maximal Von Mises stresses, obtained in the two cases corresponding to the optimization carried out with and without considering buckling constraints. Optimization without buckling constraint leads to a structure (Figure 11) 480 with a mass of 199.52 g, which does not resist to buckling as the critical load is equal to 0.743. Taking into account the buckling constraint in the optimization problem leads to an optimized structure (Figure 12) with a mass of 210.05 g, corresponding to an increase of 5.2%.

In Figure 10, we illustrate the evolution of the critical buckling load during 485 the optimization. It can be seen that in the first iterations the values of the critical buckling load are very high. For the initial structure, the buckling analysis leads to a critical buckling load of 277.97. Recall that in the optimization al-gorithm the global topological gradient T w is a weighted sum of the topological gradients T σ and T λ . In particular, the topological gradient associated with the buckling constraint, T λ , is multiplied by the inverse of the square of the critical load (see (37)). Thus, for this application, at the initial iterations only the topological gradient T σ related to the stress constraints drives the optimization process. This explains the similarity of the structures obtained in Figures 11 and12.

The major difference is in the buckling load: its value associated to the later structure is greater than one, meaning that this structure is not subject

Conclusion

We address the topological optimization of continuum structures, whose mass is to be minimized, subject to stress and to buckling constraints.

We propose, for the buckling analysis, a variant of the inverse-iterations algorithm that is devised to compute a set of eigenvalues. It features low memory requirements, thus allowing one to handle large-scale problems, that are not uncommon in the context of topology optimization, specially when 3D discretizations are considered.

To optimize the topology of structures under stress and buckling constraints, we develop an algorithm that follows a topological gradients-based approach. We use aggregation formulations for the two considered constraints, and apply the adjoint method to efficiently compute their gradients. The proposed algorithm allows one to handle optimization of structures subject to constraints that are attracting an increasing attention in the structural engineering community. It relies on our proposed buckling analysis and gradient computation, that are embedded in a simple scheme to iteratively compute the optimized topology for the targeted structure. Numerical results on mechanical structures, including a real application arising in the aeronautical field, show the effectiveness of the proposed optimization algorithm.

The proposed topology optimization algorithm features a flexibility allowing one to easily implement further constraints, e.g., the compliance. We consider this as a perspective for future work. Moreover, the inclusion of additive manufacturing constraints, such as the trajectory of the deposition of the material in the case of the Wire Arc Additive Manufacturing method [START_REF] Radel | Skeleton arc additive manufacturing with closed loop control[END_REF]), will be investigated.

1854. 1 .

 1 Topological gradient frameworkTopology optimization problems rely on identifying optimal material distributions. Algorithms for topology optimization of mechanical structures are, thus, basically based on a principle of iterative material removal, starting from an initial domain.The notion of topological gradient, introduced by[START_REF] Eschenauer | Bubble method for topology and shape optimization of structures[END_REF], quantifies the influence, on an optimization criterion, of a perturbation of the design domain. From the point of view of a mechanical structure, this can be seen as the insertion of a hole in the structure, corresponding to an elementary amount of material that would be removed from the initial domain.Let us consider a spatial domain Ω and let Ω B ⊂ Ω be the domain obtained perturbing Ω by the introduction of a ball B(p i , r) of radius r centered in a point p i ∈ Ω. Given an optimization criterion C : Ω → R, the topological gradient of C in p i is defined as[START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF] Céa et al. (2000);[START_REF] Novotny | Topological sensitivity analysis for three-dimensional linear elasticity problem[END_REF]):

Figure 1 :

 1 Figure 1: Optimization algorithm based on topological gradients

  is the global stress component vector. Note that the two 270 formulations of G λ in (8) and (25) are equivalent. Indeed, the eigenvectors v i are K-orthoganal vectors, the term σ -Y u corresponds to the computation of the stress components from equation (3) and the term f -Ku corresponds to equilibrium equation (2).

Figure 2 :

 2 Figure 2: L-shaped structure

Figure 3 :

 3 Figure 3: Beam under a compressive load: Design space

  Finally, comparing the plots of the Von Mises stress in Figures 6(d) and 7(e), one can note that the stress values are more homogeneous for the structure optimized by considering 4 buckling modes than the values obtained considering only one buckling mode. When six buckling modes are considered, the mass of the obtained optimal structure has a value of 149.44 g, that is very close to that obtained in the case of four modes. Moreover, one can note that the distributions of the Von Mises stress in Figures 7(e) and 8(f) are also very close. Comparing the structures obtained in the two cases, displayed respectively in Figures 7(a) and 8(a), one can remark that they differ with respect to the points where the loads are applied. These structures correspond to different local optima.

Figure 4 :

 4 Figure 4: Beam under a compressive load: Evolution of the mass according to the number of buckling modes considered

Figure 9 )

 9 Figure 9), referred to as gooseneck. Several structures of this type appear on the landing gear doors of an aircraft. The area in red corresponds to a load of -12 128 N along the x axis, -35 N along the y axis, and 3 456 N along the z axis. The hatched area represents the fixed area. The structure is discretized by cubic finite elements with a length of 4mm, leading to 31 670 elements.

Figure 9 :

 9 Figure 9: Gooseneck: Design space of the structure

  to buckling failure. The first and the second buckling modes for the structures obtained respectively without buckling constraints and with such constraints are illustrated in Figure 11(c)-(d)-(e) and Figure 12(c)-(d).

Figure 10 :Figure 12 :

 1012 Figure 10: Gooseneck structure: Evolution of the critical buckling load

  m

	11:	while r p ≤ do
	12:	p = p + 1
	13:	if λ p ≥ 0 then
	14:	l = l + 1
	15:	end if
	16:

Table 1 :

 1 L-shaped structure: Eigenvalues obtained for different dimensions of the set V

	Dimension m	λ 1	λ 2	λ 3	λ 4	λ 5
	of the set V
	3	-2.222 2.677 7.038	NA	NA
	4	-2.222 2.677 7.038 -7.577	NA
	5	-2.222 2.677 7.038 -7.577 11.352

Table 2 :

 2 L-shaped structure: Computing time according to the size of the eigenvalue search space

	Dimension m	Computing time	Number of	Number of	Computing time	Number of	Computing time
	of the set V		iterations	products Kv	of one product Kv	products Kσv	of one product Kσv
	5	49min 3s	33	27 777		506
	6	31min 21s	18	28 099		328
	8	38min 47s	17	39 117		428
					0.135s		0.98s
	9	44min 21s	15	41 043		422
	10	49min 51s	15	45 089		463

Table 3 :

 3 Beam under a compressive load: Mass, Critical buckling load, Maximal Von Mises stress and Computing time

	Number of considered	Mass (g)	Critical buckling	Maximal Von Mises	Computing time
	buckling modes		load	stress
	1	178.50	1.009	299.25	3h 45 min
	4	148.17	1,07	304.94	5h 38 min
	6	149.44	1.0001	265.04	8h 55min

Table 4 :

 4 Gooseneck structure: Optimization results with and without buckling constraint

	Number of considered	Mass (g)	Critical buckling	Maximal Von	Computational time
	buckling eigenvectors		load	Mises stress
	0	199.52	0.743	497.44	30min
	1	210.05	1.07	469.18	7h 23min

+ ν σ(u(x i )) : (µ(x i )) + 1 -3ν 1 -ν

tr(σ(u(x i )))tr( (µ(x i )))(22)

Acknowledgement

The authors gratefully acknowledge the support by AVANTIS Project and the French National Association of Reasearch and Technology (ANRT).


domain: chim phys sde sdv

text:
Fabio Pietrucci
email: fabio.pietrucci@upmc.fr

Mauro Boero
email: mauro.boero@ipcms.unistra.fr

Wanda Andreoni
email: wanda.andreoni@epfl.ch

How natural materials remove heavy metals from water: mechanistic insights from molecular dynamics simulations

come

How natural materials remove heavy metals from water:

mechanistic insights from molecular dynamics simulations

INTRODUCTION

Lead, cadmium, mercury and arsenic are common pollutants of fresh and waste waters and constitute a threat to human and animal health as well as to agriculture [1][2][3] . Concern has increased dramatically over the years and led to several different strategies to reduce water contamination by heavy metals, including adsorption onto active carbon, membrane filtration and chemical precipitation [START_REF] Sharma | Heavy Metals in Water: Presence, Removal and Safety[END_REF] . An especially interesting strategy is based on inexpensive and natural materials obtained from waste produced in industrial or agricultural processes, e.g., coffee residues such as spent-coffee-grains (SCG), dried tea leaves, chestnut shells, orange or lemon peels [see e.g., 5-22]. Despite obvious economic and environmental advantages, this technology has not yet been exploited on a large scale. A vast quantity of reports on the "performance" of natural materials is scattered in the scientific literature, from which it appears, however, that understanding of the chemistry involved in the metal capture has rarely been the research focus. Attempts to identify specific responsible components have mainly relied on FTIR or XPS, namely on the comparison of the spectra of the dried material before and after immersion in polluted water [see, e.g., Refs 6, 11, 17, 23, 24]. The interpretation, however, is not easy. Even in cases when one can distinguish the response of different functional groups, one cannot ignore the fact that the same group is present in many different constituents. On the other hand, the complexity of the systems to be investigated explains the lack of simulations at the molecular level which could provide useful insights into the relevant mechanisms. The scope of the present paper is to explore the potential of state-of-the-art computer experiments for research in this field and to start by studying, via atomistic models and molecular dynamics (MD) simulations, the interaction of polysaccharides -common components of vegetables and fruit peels -with heavy metals immersed in water. More precisely, our investigation aims at clarifying how adsorption/desorption processes take place, namely at characterizing the driving mechanisms, the nature of the ion adsorption sites and the energetics. As a prototype of highly toxic metals we have chosen lead, which is mainly detected as Pb 2+ in aqueous environments. As reference polysaccharide(s) we have selected hemicellulose(s) because of its (their) special structural characteristics: amorphous, partially or fully soluble in water and heterogeneous, namely consisting of chains of diverse monosaccharides. Moreover, hemicellulose was found amongst the dominant components of SCG and coffee silverskin. 15,19,24

METHODS

Here we give a brief account of the models and methodologies used in our simulations, both classical and ab initio MD. A detailed description is given in the Electronic Supplementary Information (ESI).

Classical MD

Among the major groups of hemicellulose constituents, mannans are detected in larger percentage in SCG. They consist of β (1-4) linked D-mannose substituted every about 100 residues, at O6, with a single D-galactose (Gal). On this basis, we construct an oligomer model with the sequence (Man20GlcA10Man20)Gal(Man20GlcA10Man20)Gal(Man20GlcA10Man20)Gal(Man20GlcA 10Man20) where Man = D-Mannose, Gal = D-Galactose, GlcA = D-Glucuronic acid. This is the system simulated using classical molecular dynamics, solvated in explicit water (31782 molecules) and with 20 Pb 2+ ions that neutralize the negative charges associated with 40 carboxylate (COO) anions. Other functional groups are present in the organic that contain oxygen and can thus also act as adsorption centers for Pb, namely 570 hydroxy (C-OH) and 406 ether (-COC-) groups. As force-fields we use GLYCAM_06 25 for the organic subsystem, TIP3P 26 for water and a simple Lennard-Jones potential for the lead ions 27 .

Ab initio MD

The atomistic model simulated with ab initio MD corresponds to a cut from the one considered above. It consists of 4 organic chains, namely Man3, Man3, Man4, Man1GlcA8 (saturated with H at the ends) (842 atoms) containing 8 carboxylates (COO), 53 hydroxy (C-OH and 34 (-COC-) ether groups, 4 Pb atoms and 149 water molecules. The system is neutral and -as verified -after equilibration the Pb atoms were ionized, with Bader charges in the range 1.5-1.7. By "ab initio" MD we denote Car-Parrinello MD that is based on density-functional-theory (DFT), and, in particular, with gradient-corrected exchange-correlation functionals 28 and semi-empirical corrections for van-der-Waals interactions. 29 .

Metadynamics

Calculation of activated events implies the use of acceleration procedures to empower MD. We have chosen metadynamics 30,[START_REF]Handbook of Materials Modeling -Methods: Theory and Modeling[END_REF] (MTD) as enhanced configurational sampling method in both classical and ab initio approaches. The path collective variables (CV) [START_REF] Branduardi | [END_REF] introduced in Ref. 33 are adopted since these are uniquely suited for the study of reactions in complex systems such as aqueous solutions because all atoms in the system (metal ions, organic and solvent) can in principle be allowed to participate in the reaction under study. In this case, in addition to the Pbs, oxygens of the organic are included as well as those of all water molecules. To inspect the results of MTD simulations, we plot the free-energy surface (FES) in the CV space reconstructed from the bias potential. [START_REF] Pietrucci | [END_REF] We further analyze it in terms of the Pb-O coordination number (see the ESI) relative to the first shell, with a cutoff of 3.3 Å.

We remark that the values estimated for the free-energy barriers (FEBs) are affected by a statistical error of at least 5 kcal/mol, with a tendency to overestimation because the FES is estimated from the bias potential accumulated until the "forward" (or the "backward") transition occurs, without attaining the regime of multiple forward and backward transitions.

RESULTS and DISCUSSION

The key questions we tackle are: How easy is the metal capture? Do all types of oxygen-containing groups in the model act as adsorption sites and if yes to what extent? How easy is the metal ion detachment from the organic? How do water molecules participate in these reactions? First, three distinct unbiased MD simulations were run at room temperature, starting from configurations with all metal ions fully solvated and far away from the organic. As shown in Figure 1 (and Figures S1 andS2), adsorption occurs spontaneously within tens of nanoseconds. Moreover, as illustrated in the examples of Figure 1, four-five oxygen atoms succeed in binding to a given Pb ion within 100 ns and all three types of functional groups containing oxygen contribute to the metal capture to some extent although, as expected from their larger partial negative charge, carboxylates are the principal chelating agents. As a second step, the desorption process was investigated via four independent MD-MTD simulations -also at room temperature -starting with Pb 2+ located in different environments. We employ the S path coordinate (see SI) as the main indicator of the "degree of adsorption": the region with S<1.2 includes configurations with Pb 2+ tightly bound to the organic (Figure 2A), while the region with S>1.6 includes configurations with Pb 2+ dissociated from the organic and fully solvated by water molecules (Figure 2D). The region 1.2<S<1.6 (see Figures 2B and2C) includes transition pathways between A and D. Figure 3 shows the FES reconstructed from two (a,b) out of four (a-d; see Figure S3) simulations. As expected, the bound state corresponds to a deep free-energy minimum, while the transition region towards the dissociated state features an almost monotonic growth of the free energy. The relative stability of the dissociated state cannot be quantified from MTD simulations due to the absence of reverse (association) transitions. However, as we have reported above, unbiased MD clearly indicates the lack of significant barriers for association. On the contrary, all four MTD simulations reveal that desorption involves sizable barriers: The estimated values from the bias accumulated until full dissociation are 32 (a), 45 (b), 56 (c) and 32 kcal/mol (d).

In order to understand the link between the free-energy profile and the structural evolution of the system, we exploit the correlation between the path coordinate S and the Pb 2+ -O coordination pattern. The latter can be further dissected into the coordination numbers of Pb 2+ with oxygen atoms belonging to its different potential adsorption species: water or the organic via carboxylate-(-COO -), hydroxy-(C-OH) and ether groups (-COC-) (Figures 3 andS3). In the bound state (S<1.2), carboxylate groups are the main chelating agents (coordination between 4 and 6) in agreement with the results of the simulations of the capture process. Hydroxy groups also make an important contribution (coordination between 1 and 4) while ether groups are seldom in contact. Moreover, water molecules (up to 4) coordinate the cation, indicating that the binding pocket is not dry. The loss of contacts with organic groups during the unbinding is systematically compensated by that of water molecules, until full solvation. To a good approximation, Pb 2+ is constantly coordinated to 7-9 oxygen atoms, in agreement with the coordination number computed from the pair distribution function of the solvated ion (Figure S4). Analysis of the MD trajectories of solvated Pb 2+ indicates an average lifetime of 189 ± 27 ps for water molecules in the first shell (with standard deviation 186 ps, i.e. close to the average as expected for Poisson statistics of rare water unbinding events). For comparison with pure water, the lifetime of a tagged molecule in the first shell is 8.3 ± 1.6 ps (with standard deviation 7.5 ps). We remark that such a long lifetime for the first solvation shell hampers the brute-force sampling of the dynamics of Pb 2+-in water and justifies our MTD protocol based on CVs that include the solvent. Turning on to discuss the results of ab initio MD, some more clarifications of our model and protocol are necessary. DFT calculations are expected to provide a better description of the relevant intermolecular interactions (Pb-water; Pb-organic, water-organic, water-water etc.). Still, when comparison is made with classical MD, one should care about the effects of the unavoidable limitations in the model size and in the duration of these simulations (20-40 times faster than in classical MD). First, we notice that, although offering several possibilities for Pb adsorption, the size limitations do not allow more than one hydration shell to complete after detachment from the organic. However, we can characterize in detail the early steps of the desorption and of the reverse reaction. MTD was applied, at room temperature, to two out of four Pb ions having different local environments. The FEBs for desorption (Figures 4 andS5) estimated from the bias accumulated until the most dissociated state are 103 and 59 kcal/mol. On the other hand, the estimated FEBs relative to the removal of the metal ions from solvated configurations -end points of the unbinding simulations -and re-capture by the organic are 30 and 10 kcal/mol (Figures 5 andS6). Inspection of the atomic trajectories suggests that these barriers correspond to the displacement of water molecules from the path of approach of the metal ion to the organic.

The slow relaxation of the solvent is expected to cause a more severe overestimation of the FEBs relative to classical MD. In order to quantify this statement, we have repeated the four classical MD simulations in Figure 3 with the MTD setup adjusted (SI) to that of ab initio simulations. The values thus estimated from the bias potentials (Figure S7) are indeed higher by a factor of at least 1.6. Once applied to our estimates from DFT-MD, this result explains at large the seeming quantitative discrepancy with classical MD. From comparison of the results of all our simulations, a coherent picture emerges of the equilibrium configurations as well as of binding and unbinding processes.

(i)

The tendency of the Pb ion to bind to diverse oxygen-based functional groups -that we observed in the extensive classical MD simulations -is also observed in the ab initio simulations of the capture in Figure 4. These observations agree with several FTIR and XPS data on bio-organic materials with adsorbed Pb (see, e.g., Refs 11, 35). However, recent XPS data on SCGs were interpreted as demonstrating an exclusive binding of Pb to carboxylates 18 . This argument was further emphasized as support for the application of the Langmuir model in its simplest form (single-site adsorption), assuming that all adsorbents are identical. However, the spectra resolution does not seem to be adequate to support these conclusions. [START_REF] Lopez | O1s spectra are deconvoluted into two (instead of typical three) components and ascribed to oxygen atoms involved in either single (C-O) or double (C=O) bonds. Only the peak position of the latter appears to be sensitive to the interaction with the metal ions, whereas the height and surface areas of both contributions are affected[END_REF] Vibrational spectra calculated from the ab initio MD trajectories before and after injection of the metal ions show that the major changes induced by Pb are blue-shifts of the stretching modes of carbon-oxygen bonds in carboxylates (Figure S8). These results agree with the interpretation of infrared spectra. 6,37 (ii)

The adsorption locations are never dry, namely water molecules participate in the Pborganic bound states. This is reminiscent of other biological environments, such as active sites in protein-ligand complexes or at protein-protein interfaces [38][39][40] . (iii) The gradual replacement of the organic residues with water molecules in the first coordination shell accompanies the ion removal from the organic and vice versa in the adsorption process, with the bound state corresponding to a coordination of 8-9 in classical MD (Figures 3 andS3) and 8 or 6 in ab initio MD (Figures 4 andS5). The reduced availability of binding sites in the latter case explains the difference.

(iv)

The uptake from the organic is easy, fast and irreversible, as can be expected for an efficient capture.

CONCLUSIONS

Several bioorganic materials such as SCG are known to be efficient agents for the capture of toxic metals from water but knowledge of their functioning at the molecular level is missing. Motivated by the importance of the related issues we have used modeling and advanced molecular dynamics protocols to gain new insights into the dynamics of the uptake of lead ions from aqueous solutions by hemicellulose -a most relevant component of SCG -as well as their release back to water.

In conclusion, our simulations were able to predict -from scratch -that the polymer can indeed act as a very efficient capture agent for lead ions, namely in a fast and irreversible manner: By offering multiple and effective adsorption sites (oxygen-based functional groups), components like hemicellulose -as our simulations show -can easily capture the metal ions at room temperature whereas overcoming a range of sizable barriers is required for their desorption under the same conditions. Moreover, the twofold role of water emerges, namely not only as a competitor of the polymer for the uptake of the metal ions, but also as its "partner" in the bound metal-polymer states, thus satisfying the eagerness of lead for high coordination with oxygen atoms.

From a methodological viewpoint: The need of a metadynamics protocol like ours, namely including all solvent molecules as potential participants in the events, is proven. It becomes also clear that improving on quantitative assessments will have to mainly rely on methodological advances in alleviating the time-scale problem, [START_REF]Handbook of Materials Modeling -Methods: Theory and Modeling[END_REF][START_REF]Handbook of Materials Modeling -Methods: Theory and Modeling[END_REF] given the slow reorganization of the hydrogenbond network to accommodate changes in coordination with the metal ion.

Although our study has characterized model systems, the knowledge we have gained can be extended to more complex scenarios. We believe that the deepening of our understanding of the physics and chemistry underlying water pollution by heavy metals or organic waste is essential for progress on the effective utilization and the design of removal agents. In this challenging task, application of state-of-the-art simulation methods on simplified but still realistic models could play a central role, as is currently the case for protein science 42,43 and drug discovery 44 .

Figure 1 :

 1 Figure 1: Unbiased classical MD: Examples of unbiased spontaneous binding of two out of 20 Pb ions in solution. The coordination number corresponds to the nearest neighbors (see text).

Figure 2 :

 2 Figure 2: Classical MD-MTD. Snapshots of configurations with Pb (blue) strongly bound to the organic (A), along the unbinding path (B, C) and fully solvated in water (D). Note that the light blue background represents the explicit solvent of ~32000 molecules, of which only those in the first coordination shell of one metal atom are shown.

Figure 3 .

 3 Figure 3. Upper panels: Free energy surface for the unbinding of Pb reconstructed from two out of four (see Figure S3) independent classical metadynamics simulations in different environments. The shaded area (S>1.8) represents the region where the bias profile is not representative of the free energy, due to lack of reverse transition (binding). Note in (a) the location of the configurations of Figure 2. Lower panels: Evolution of coordination number of Pb and O atoms of water molecules and different groups in the organic (see legend) as a function of the S path coordinate. Vertical bars represent standard deviations.

Figure 4 .

 4 Figure 4. Upper panel: Free energy surface for the unbinding of Pb reconstructed from one out of two (Figure S5) independent ab initio metadynamics simulations in different environments. The shaded area (S>1.7) represents the region where the bias profile is not representative of the free energy, due to lack of reverse transition (binding). Lower panel: Evolution of the coordination numbers of Pb and O belonging to water and to different organic functional groups, as a function of the S path coordinate. The vertical bars represent standard deviations of coordination numbers.

Figure 5 .

 5 Figure 5. Upper panel: Free energy surface for the binding of Pb reconstructed from one out of two ab initio metadynamics simulations (see Figure S6 for both) started from the end (unbound state, S>1.8) of the simulation in

Figure 4 .

 4 Figure 4. In the shaded area -S<1.2 -the bias is not representative of the free-energy. Lower panel: Evolution of the coordination of Pb with oxygen atoms belonging to water and to different organic functional groups as function of the S path coordinate.

Acknowledgments

The research presented here was made possible thanks to the support from PRACE (WASHME; PRA12-3086) at CINECA. Other computer grants are also acknowledged: from Grand Equipement National de Calcul Intensif (GENCI) (under allocation DARI-A0060906092; A008090609; A0050910143; A0070811069) and ISCRA at CINECA. We wish to thank Roberto Cingolani, Athanassia Athanassiou and Despina Fragouli for stimulating discussions and Massimiliano Guarrasi and Piero Lanucara for technical help at CINECA.

Footnotes

Electronic Supporting Information available: Additional information on structural models and computational methods, details of the metadynamics protocols, additional illustrations of the simulation results (Figures S1-S8) discussed in the main text.

The authors contributed equally to this work.

The authors declare no competing financial interest.