This paper deals with the applications of stochastic spectral methods for structural topology optimization in the presence of uncertainties. A non-intru- sive polynomial chaos expansion is integrated into a topology optimization algo- rithm to calculate low-order statistical moments of the mechanical-mathematical model response. This procedure, known as robust topology optimization, can op- timize the mean of the compliance while simultaneously minimizing its standard deviation. In order to address possible variabilities in the loads applied to the mechanical system of interest, magnitude and direction of the external forces are assumed to be uncertain. In this probabilistic framework, forces are described as a random eld or a set of random variables. Representation of the random ob- jects and propagation of load uncertainties through the model are eciently done through Karhunen{Lo eve and polynomial chaos expansions. We take advantage of using polygonal elements, which have been shown to be e ective in suppressing checkerboard patterns and reducing mesh dependency in the solution of topology optimization problems. Accuracy and applicability of the proposed methodology are demonstrated by means of several topology optimization examples. The ob- tained results, which are in excellent agreement with reference solutions computed via Monte Carlo method, show that load uncertainties play an important role in optimal design of structural systems, so that they must be taken into account to ensure a reliable optimization process. Keywords topology optimization stochastic spectral approach polynomial chaosKarhunen{Lo eve expansion robust optimization polygonal nite element N. Cuellar A. Pereira I. F. M. Menezes Ponti cal Catholic University of Rio de Janeiro { PUC{Rio, Department of Mechanical Engi- neering, Rua Marqu^ es de S~ ao Vicente, 225, Rio de Janeiro, 22453-900, RJ, Brazil E-mail: ncuellar@tecgraf.puc-rio.br E-mail: anderson@tecgraf.puc-rio.br E-mail: ivan@puc-rio.br A. Cunha Jr Universidade do Estado do Rio de Janeiro { UERJ, Nucleus of Modeling and Experimentation with Computers { NUMERICO, Rua S~ ao Francisco Xavier, 524, Rio de Janeiro, 20550-900, RJ, Brazil E-mail: americo@ime.uerj.br 2 N. Cuellar et al. 1 Introduction Due to new requirements of design associated with the most modern engineering applications, mechanical systems with very complex geometrical con gurations are becoming increasingly common. In this context, some of the most promising design approaches are based on topology optimization (TO), which seeks to nd the best layout for a system, by optimizing the material distribution in a prede ned design domain [35, 8]. The growing popularity of TO solutions is demonstrated by their wide range of application in various elds such as structural mechanics [58, 59, 45, 43, 17, 33], composite and multi-materials [62, 67], nanotecnology [36], uid mechanics [39, 19], uid-structure interaction [1], medicine [38], etc. In general, the physical systems underlying TO applications are subjected to a series of uncertainties (e.g. unknown loads, geometrical imperfections, uctuations in physical properties, etc) so that, usually, their response is not well predicted by the traditional (deterministic) tools of engineering analysis. For this reason, there is a consensus among computational engineering experts that uncertainties e ects must be incorporated into any computational predictive model [49]. The modeling and quanti cation of uncertainties is necessary in order to predict a possible range of variability for the mathematical model response, and to conduct applied tasks, such as analysis and design, in a robust way [5, 51]. Notice, however, that the majority of the works in TO area, currently available in the literature, are limited to deterministic analyses. The need for robust design and analysis of uncertainties in topological opti- mization applications naturally induces the search for computationally ecient frameworks for TO. For this purpose, TO literature started to take uncertainty quanti cation (UQ) into account over the last decade, as can be seen in several pa- pers addressing the two issues [65, 28, 16, 66, 42, 44, 68, 20, 27, 60, 3, 14, 25, 29, 61]. Some of these works are based on classical techniques for stochastic computa- tion like Monte Carlo (MC) method [14, 20] or series expansion [16, 68, 27, 3, 25] which, despite of being very simple in conceptual terms, they are limited by the high computational cost, the former, or very small range of applicability, the latter. These limitations open space for spectral-based approaches [29, 60, 44, 65, 42, 28], that use state-of-the-art tools for representing and propagating uncertainties in computational models, like Karhunen{Lo eve (KL) and generalized polynomial chaos (gPC) expansions. A recent work by Keshavarzzadeh et. al [28] presents a non-intrusive gPC strategy to propagate uncertainties in topology optimization problems. They use non-intrusive polynomial chaos expansion to evaluate low- order statistics of compliance and volume and the uncertainties are considered in the applied loads and also in the geometry of the problems. The classical formulation for topology optimization, which corresponds to min- imize the structural compliance, is commonly carried out on uniform grids con- sisting of Lagrangian-type nite elements (e.g., linear quads). However, this choice of discretization, together with density methods, su er from the well-known nu- merical instabilities, such as the checkerboard patterns. Unstructured "Voronoi" meshes, generated from an initial set of random points, have been shown to be e ective in suppressing checkerboard patterns [58]. Moreover, compared to stan- dard Lagrangian-type uniform grids, polygonal elements are more versatile in dis- cretizing complex domains and in reducing mesh dependency in the solutions of topology optimization [58, 2]. The geometrical exibility of the polygonal nite Title Suppressed Due to Excessive Length 3 elements also make them very attractive for adaptive mesh re nement schemes in topology optimization problems [37, 26]. The computational code used here was developed based on PolyTop [59], a MATLAB code for solving topology op- timization problems using either structured or unstructured polygonal meshes in arbitrary two-dimensional domains. The modular structure of PolyTop , where the analysis routine and optimization algorithm are separated from the choice of topol- ogy optimization formulation, together with a non-intrusive way of computing the statistical measures, allowed us to implement a robust topology optimization code in a very straightforward way, with only a few modi cations in the original PolyTop code. The aim of this paper is to present a computationally ecient and accurate non-intrusive robust topology optimization approach using polygonal elements. For this purpose, the PolyTop framework by Talischi et al. [59], which employs polygonal nite elements in TO, is combined with a consistent methodology for stochastic analysis that uses a non-intrusive gPC strategy to propagate parame- ters uncertainties through the computational model. This combination generates a framework that is computationally ecient for stochastic simulations, and robust to numerical instabilities typical of TO problems, such as checkerboards, one-node connections, and mesh dependency. The novel approach is used to solve TO prob- lems that seek to minimize an objective function based on the low-order statistical moments of the compliance function of a structure, subjected to uncertainties on the external load, and satisfying volume constraints. The remainder of this paper is organized as follows. Stochastic spectral meth- ods are introduced in section 2, together with the mathematical formulation and basic steps to obtain KL and gPC expansions. In section 3, the TO problem is brie y described, as well as stochastic procedure to propagate uncertainties within TO, which is called robust topology optimization (RTO). In section 4, numerical examples are presented, and the proposed methodology is compared with a refer- ence result obtained with MC method. Finally, some remarks and suggestions for future work are presented in section 5. 2 Stochastic spectral methods 2.1 Preliminary de nitions and notation Consider a probability space ( ;F;P), whereis the sample space, Fa- eld over, andP:F ! [0;1] denotes the probability measure. It is assumed that the distribution PX(dx) of any real-valued random variable Xin this probability space admits a density x7!pX(x) with respect to dx. The set of values where this density is not zero is dubbed the support of X, being denoted by SuppX. In this probabilistic setting, any realization of random variable Xis denoted byX() for2, and the mathematical expectation operator is de ned by EfXg=Z Rx PX(dx); (1) so that the mean value and standard deviation of Xare given by X=EfXgand X= (En X2o EfXg2)1=2, respectively. The random variable Xis said to be of 4 N. Cuellar et al. second-order if En jXj2o =Z Rjxj2PX(dx)<+1: (2) The space of all second-order random variables in ( ;F;P), denoted by L2(;P), is a Hilbert space [12] equipped with the inner product h;i:L2(;P)L2(;P)! Rsuch that hX;Yi=EfXYg=Z Z R2xy PX;Y(dx;dy ); (3) wherePX;Y(dx;dy ) is the joint distribution of the random variables X;Y2L2(;P). This inner product induces a norm kk:L2(;P)!Rwhere kXk= hX;Xi1=2=Z Rjxj2PX(dx)1=2 : (4) Further ahead it will also be helpful to consider L2(D), the set of all real-valued square integrable functions de ned on the spatial domain DRd; d1. This set of functions is also a Hilbert space [12], with inner product h;iD:L2(D)L2(D)!R de ned byD ; 0E D=Z D (x) 0(x)dx; (5) for ; 02L2(D). 2.2 Karhunen{Lo eve expansion The KL expansion [24, 63] is one of the most widely used and powerful techniques for analysis and synthesis of random elds, providing a denumerable represen- tation, in terms of the spectral decomposition of the correlation function, for a random eld parametrized by a nondenumerable index [51, 53]. Let the map ( x;)2D7!U(x;)2Rbe an arbitrary real-valued random eld, indexed by the spatial coordinate vector x2D  Rd; d1, denoted in an abbreviated way as U(x) orU. By construction, for a xed x2D,U(x;) is a real-valued random variable, and U(;), for a xed 2, is a function of x, dubbed realization of the random eld. The correlation of U(x) is the function KU(;):DD! Rde ned for any pair of vectors xandx0by means of KU(x;x0) =En U(x)U(x0)o : (6) Suppose that random eld U(x) is second-ordered and mean-square continuous, properties respectively de ned by En jU(x)j2o <+1;8x2D; (7) and limx0!x U(x0)U(x) 2 = 0: (8) Under these assumptions, the linear integral operator KU (x) =Z DKU(x;x0) (x0)dx0(9) Title Suppressed Due to Excessive Length 5 de nes a Hilbert-Schmidt operator [53, 51], which has denumerable family of eigen- pairs(n; n) +1 n=1such that Z DKU(x;x0) n(x0)dx0=n n(x); x2D; (10) wherenare the eigenvalues and nthe corresponding eigenfunctions of the op- erator de ned by Eq.(9). Besides that, the sequence of eigenvalues is such thatP+1 n=1n<+1and12n! 0; and the family of functions f ng+1 n=1de nes an orthonormal Hilbertian basis in L2(D), i.e. h m; niD=mn; (11) where Kronecker delta is such that mn= 1 ifm=nandmn= 0 form6=n. Therefore, applying two standard results of functional analysis [15], the theo- rems of Hilbertian basis and orthogonal projection, it is possible to show that the random eld U(x) admits a decomposition U(x) =U(x) ++1X n=1p n n(x)n; (12) wherefng+1 n=1is a family of random variables de ned by n=1pnhUU; niD; (13) which are centered (zero mean) and mutually uncorrelated i.e. n= 0; andEfmng=mn: (14) A nite dimensional approximation for U(x), denoted by Ukl(x), is con- structed by truncation of the series in Eq.(12) i.e. Ukl(x) =U(x) +klX n=1p n n(x)n; (15) where the integer klis chosen such that energy (kl) =Pkl n=1nP+1 n=1n ; (16) for a heuristically chosen threshold (e.g. = 90%). In practice, as a closed formula for nis not available in general, energy (kl) is estimated using a nite (but large) number of eigenvalues, instead of an in nite quantity. This procedure is justi ed in light of the eigenvalues decreasing property. One of the main diculties to apply KL expansion to discrete random elds is the determination of the eigenvalues and corresponding eigenfunctions of the correlation function. Analytical solutions for the Fredholm integral equation in (10) are almost never available. However, for some special cases, such as exponential and Gaussian autocovariance functions, an analytical solution can be obtained by converting the integral equation into a di erential equation through successive derivatives [24, 63]. 6 N. Cuellar et al. Several numerical methods can be used to solve the eigenvalue problem of Eq.(10), such as the direct method, projection methods, among others [4, 9]. In this study, the direct method is employed to transform the Fredholm integral equa- tion into a nite dimensional eigenvalue problem, whose the solution provides an approximation for the desired eigenvalues/eigenvectors of the in nite dimensional problem. In this numerical procedure, a set of Mrealizations of the random eld U(x) and its mean function U(x) are numerically generated1and grouped into the matrices U=2 66664U1(x1)U2(x1)::: UM(x1) U1(x2)U2(x2)::: UM(x2) ............ U1(xn)U2(xn)::: UM(xn)3 77775;=2 66664U(x1)U(x1)::: U(x1) U(x2)U(x2)::: U(x2) ............ U(xn)U(xn)::: U(xn)3 77775;(17) which are used to de ne the zero mean matrix ^U=U. Then, the correlation matrix is estimated with the aid of K^U=1 M^U^UT; (18) and the discrete eigenvalue problem K^U=; (19) is solved to obtain the matrices and, which present approximations for the rstMeigenfunctions/eigenvalues on the columns/main diagonal. For further details on theoretical and practical aspects of KL expansion the reader is encouraged to see [34, 6, 40, 41]. 2.3 Generalized polynomial chaos expansion The gPC expansion is a theoretical tool used to construct representations for random elds, with a denumerable or nondenumerable set of index, in terms of a denumerable collection of random variables weigthed by deterministic coecients [51, 53]. It was introduced in the engineering community by R. Ghanem [24, 55, 22] as a tool to compute approximate responses for problems involving random elds with unknown distribution and since the early 2000s, especially after the work of Xiu and Karniadakis [64], it has been used in many applications of computational stochastic mechanics, a trend that should increase [56]. For the sake of theoretical development, consider a second-ordered random variableU:!Rwhich can be written in terms of a (possibly in nity) set of independent random variables 1();2();3();. Collecting these independent variables into the random vector () = (1();2();3();), dubbed the germ, it is possible to rewrite the original variable in the parametric form U=U(). 1These numerical realizations are de ned in a computational mesh x1;x2;  ;xn. Title Suppressed Due to Excessive Length 7 In this context of second-ordered random variables, gPC expansion theory says that such random variable Uadmits a spectral representation U() =+1X n=0un n(); (20) where the set of orthonormal polynomials f 0; 1; 2;gis a basis for L2(;P), and the deterministic coecients u0;u1;u2;correspond to the coordinates of U in this in nite dimensional base. By de nition 0() = 1, and due to the orthogo- nality propertyh m; ni=mn, one has U=EfUg (21) =Ef 0Ug =+1X n=0unh 0; ni =u0; i.e., the mean of Uis equal to the coordinate u0of gPC expansion. Besides that, it is easy to see that Ef ng= 0; n1. Thus, it can also be shown that standard deviation of Ucan be written as 2 U=En U2o EfUg2(22) =E8 >< >:0 @+1X n=1un n()1 A29 >= >; =+1X n=1u2 n: The gPC expansion has attracted the attention of many researchers due to its rapid convergence property and its capability to estimate statistical moments, which allows this method to eciently reduce computational e ort in highly non- linear engineering design applications [34, 41]. Of course, for any purpose of numerical computation, it is necessary to parametrize the random variable Uwith a nite number rvof independent random variables, i.e.,= (1;2;;rv), and restrict up to ppcthe order of the polynomials in the basis, so that the series in Eq.(20) is truncated with pc+1 = (rv+ppc)!=(rv!ppc!) terms, giving rise to the approximation Upc() =pcX n=0un n(); (23) which is mean-square convergent to Uwhenppc;rv!+1. In practice, rvis obtained from the stochastic modeling and ppcis speci ed such that a compromise between accuracy and eciency can be established. In this work, the gPC polynomial order is selected such that the low-order statistics estimations become invariant for increasing values of ppc. This spectral expansion can be easily extended to a second-ordered random eldU:D!R, considering this eld as a family of random variables 8 N. Cuellar et al. U(x;)2L2(;P), parameterized by the index x, and by letting the determin- istic coecients of the gPC expansion depend on x, i.e. U(x;) =+1X n=0un(x) n(); (24) where the coordinates are now a deterministic function of x. As in the case of a random variable, the truncation of the series in Eq.(24) results in the approxima- tion Upc(x;) =pcX n=0un(x) n(): (25) Basically, two approaches for implementing the gPC expansion are available in the literature, one called the intrusive, and another one dubbed non-intrusive. The intrusive approach is based on a stochastic version of Galerkin method [63], and consists in modifying the deterministic model in order to take into account the uncertainty propagation. The main disadvantage of this formalism is related to the level of diculty in modifying the code associated with the deterministic model, particularly if commercial software is being used, due to code access restrictions [51, 23]. On the other hand, the non-intrusive approach requires no modi cation in the deterministic model and therefore, can be treated as a black box [51, 23]. In this case, a probabilistic collocation approach, based on a sparse grid method [21], can be used to estimate the coecients (coordinates) of the expansion. In this work, a di erent non-intrusive approach based on linear regression is employed, where the random eld (variable) of interest is evaluated in a nite set of gqpossible realizations of the germ , and thus the coecients of the expansion are obtained through u= ( T )1 TU, the solution of the mean-square problem uU; (26) where2 664 0(1) pc(1) ......... 0 gq  pc gq3 775 |{z } 2 66664u0 u1 ... upc3 77775 |{z} u2 664U(x;1) ... U(x;gq)3 775: |{z} U(27) For further information about the basic aspects of gPC expansion the reader is encouraged to see the