## 1 Introduction

With the rise of additive manufacturing, topology optimization becomes a popular design methodology to determine the optimal distribution of materials in complex engineering structures(yunkang1998new, ; wang2003level, ; wang2004color, ; liu2016multi, ). Inevitable uncertainties in the additive manufacturing process and operating environment often undermine the performance of such topology designs. Classical deterministic design approaches often lead to unknowingly risky designs due to the underestimation of uncertainties, or inefficient and conservative designs that overcompensate for uncertainties. In the past decade, robust topology optimization (RTO) and reliability-based topology optimization (RBTO) are increasingly adopted as an enabling technology for topology design subject to uncertainty in aerospace, automotive, civil engineering, and additive manufacturing (chen2010level, ; chen2011new, ; asadpoure2011robust, ; guo2013robust, ; guo2015multi, ; zhang2016robust, ; jiang2017parametric, ). The former seeks for insensitive topology design via minimizing the propagation of input uncertainty, whereas the latter delivers reliable topology design by introducing probabilistic characterizations of response functions into the objective and/or constraints.

RTO and RBTO for realistic engineering applications confront two challenges: (1) the theoretically infinite-dimensional design vector; and (2) high-dimensional integration resulted from a large number of random variables. Both lead to the curse of dimensionality, which hinders or invalidates almost all RTO and RBTO methods. In RTO, the objective or constraint functions are usually expressed by first two moment properties, such as means and standard deviations, of certain stochastic responses, describing the objective robustness or feasibility robustness of a given topology. RBTO often contains probabilistic constraint functions, which restrict the probability of failure regarding certain failure mechanisms. Therefore, to solve a practical RTO or RBTO problem using gradient-based algorithm, an efficient and accurate method for statistical moments, reliability, and their sensitivity analysis for random responses are in demands.

The fundamental problem rooted in statistical moment or reliability analyses entails the evaluation of a high-dimensional integral in the entire support of random inputs or its unknown subdomain, respectively. In general, such an integral cannot be calculated analytically. Direct numerical quadrature can be applied, but it is computationally prohibitive when the number of random inputs exceeds three or four, especially when the evaluation of response function is carried out by expensive finite element analyses (FEA). Existing approaches for statistical moment and reliability analysis include the point estimate method (PEM) (huang2007, ), Taylor series expansion or perturbation method (huang2007, )

, tensor product quadrature (TPQ)

(lee2009, ), Neumann expansion method (yamazaki88, ), the first-order reliability method or FORM-based methods (kuschel1997, ; tu1999, ; du2004seq, ; chiralaksanakul2005, ; agarwal2006, ; liang2007, ), polynomial chaos expansion (PCE) (wang2006, ), statistically equivalent solution (grigoriu91, ), dimension-reduction method (xu04, ; xu05, ), and others (grigoriu02, ). Their topology sensitivities have relied mainly on two kinds of approaches: SIMP-based approaches (solid isotropic material with penalization) (amstutz2012topological, ) and topology-derivative-based approaches (novotny2003topological, ; allaire2005structural, ). The former is based on a fictitious density field representing a smooth transition between material and empty, which requires regularization procedures to get a clear topology. The latter introduces the topological derivative concept which defines the derivative of functionals whose variable is a geometrical domain with respect to singular topology perturbation. The topological derivative concept is mathematically rigorous and independent of the fictitious density field.Nonetheless, three major concerns arise when evaluating stochastic quantities and their sensitivity using existing approaches or techniques. First, when applied to large-scale topology optimization subjected to a large number of random inputs, many of those methods including Taylor series expansions, FORM-based methods, PEM, PCE, TPQ, and dimension-reduction methods, etc. begin to be inapplicable or inadequate. For example, although the Taylor series expansion, FORM-based methods, and PEM are inexpensive and simple, they deteriorate due to the lack of accuracy when the nonlinearity of a response function is high and/or when the input uncertainty is large. PCE approximates the random response via an infinite series of Hermite polynomials of Gaussian variables (or others) and is popular in stochastic mechanics in the last decades. However, when applied to high-dimensional systems, it is easily succumbed to the curse of dimensionality due to astronomically large numbers of terms or coefficients required to capture the interaction effects between random inputs. Rooted in the referential dimensional decomposition (RDD), the dimension-reduction approximates a high dimensional function via a set of low dimensional components, but it often results in sub-optimal estimations of the original function, and thus its stochastic moments and associated reliability. Second, to evaluate the topology sensitivity of stochastic quantities, many of the aforementioned methods may not be adequately efficient and accurate. Most of those methods rely on a fictitious density field, thus the sensitivities supplied are not the exact topology sensitivity. Furthermore, many of them resort to repetitive stochastic analyses especially for the sensitivity of reliability due to employed finite-difference techniques, which restrain their computational efficiency. Although Taylor series expansions, is able to perform stochastic sensitivities analysis economically, its accuracy is usually deteriorated by inherited errors from the associated stochastic analysis. Third, to the best of the author’s knowledge, in existing literature, there is no benchmark example that provides analytical or semi-analytical solutions for stochastic topology sensitivity analysis. A successful benchmark example certainly calls for analytical expressions of stress, strain, or other response functions in two domains - an original domain and a perforated domain - subject to the same loads and supports. These analytical expressions generally are not readily available even for simple domain and ordinary load cases. Furthermore verifying the performance of a certain method subject to high-dimensional random inputs often demands the benchmark example carrying on complex loads to accommodate a large number of random variables, which impede the implementation of analytical solution of stochastic topology sensitivity. These difficulties result in the lack of benchmark examples and make it impossible to verify the accuracy of existing and new algorithms, especially for high-dimensional cases.

This paper presents a novel framework for topology sensitivity analysis of statistical moments and reliability for complex engineering structures subject to a large number of random inputs. The framework, designed for dealing high-dimensional random inputs, is grounded on the polynomial dimensional decomposition (PDD), and thus it is capable of approximating the high-dimensional stochastic responses in an efficient and accurate manner. It also dovetails the deterministic topology derivatives with PDD and provides stochastic sensitivities in the exactly topological sense. For RTO, the proposed framework is endowed with analytical expressions for topology sensitivities of the first three stochastic moments. For RBTO, it supplies embedded Monte Carlo Simulation (MCS) and finite difference formulations to estimate topology sensitivities of failure probability. Furthermore, the evaluation of moments and/or reliability and their topology sensitivity is accomplished concurrently from only a single stochastic analysis. It is noteworthy that two benchmark examples, which provide analytical/semi-analytical topology sensitivity of moments and reliability, are developed for the calibration of stochastic topology sensitivity algorithms. The first example contains only two random variables but provides analytical expressions for moments, reliability, and their topology sensitivities. The second one accommodates 53 random variables, whereas the analytical expressions provided can be easily expanded to any positive number of random variables. The rest of this paper is organized as follows. Section 2 formally defines general RTO and RBTO problems, including a concomitant mathematical statement. Section 3 starts with a brief exposition of the polynomial dimensional decomposition and associated approximations, which result in explicit formulae for the first two moments and an embedded MCS formulation for the reliability of a generic stochastic response. Section 4 revisits the definition of topology derivative and describes the new framework of stochastic topology sensitivity analysis, which integrates PDD and deterministic topological derivative as well as numerical procedures for topology sensitivities of both stochastic moment and reliability. The calculation of PDD expansion coefficients is briefly described in Section 5. Section 6 presents three numerical examples. Two benchmark examples are developed to probe the accuracy and computational efforts of the proposed method. One three-dimensional bracket is used to demonstrate the feasibility of the new method for practical engineering applications. Finally, conclusions are drawn in Section 7.

## 2 Stochastic topology design problems

In the presence of uncertainties, the topology optimization problem often includes robust and/or probabilistic constraints. For RTO, both objective and constraint functions may involve the first two moment properties for the assessment of robustness. Whereas for RBTO, probabilistic functions are often embedded as constraints to restrict the failure probability and achieve a high confidence level on design. Nonetheless, a generic RTO and a generic RBTO problem are often formulated as the following mathematical programming problems

subject to | ||||

and

subject to | ||||

respectively, where is a bounded domain in which all admissible topology design are included; is an

-dimensional random input vector completely defined by a family of joint probability density functions

on the probability triple , where is the sample space; is the -field on ; is the probability measure associated with probability density ; is the th failure domain defined by response function ; expresses target failure probabilities; and are two non-negative, real-valued weights, satisfying , and are two non-zero, real-valued scaling factors; , , are non-negative, real-valued constants associated with the probabilities of constraint satisfaction; andare expectation operator and variance operator, respectively, with respect to the probability measure

. The evaluation of both and on certain random response demands statistical moment analysis (rosenblueth1981, ; hong1998, ; Kleiber92, ; rahman01p, ; evans1967, ; yamazaki88, ; grigoriu91, ; xu04, ; xu05, ; grigoriu02, ), which is not unduly difficult. In contrast, the evaluation of probabilistic constraint functions in RBTO, generally more complicated than and , is obtained from(3) |

which represents a failure probability from reliability analysis (cornell1969, ; madsen1986, ; hasofer1974, ; fiessler1979, ; hohenbichler1981, ; breitung1984, ; der1987, ; wu1990, ; wu1987, ; liu1991, ). The indicator function when and zero otherwise. For component-level RBTO, the failure domain, often adequately described by a single performance function , and component reliability analysis are relatively simple. In contrast, interdependent performance functions , are required for a system-level (series, parallel, or general) RBTO, leading to a highly complex failure domain and huge computational demand for system reliability analysis.

## 3 Polynomial dimensional decomposition method and uncertainty quantification

### 3.1 Polynomial dimensional decomposition

Consider a multivariate stochastic response of certain topology design subject to random input vector , representing any of the performance function in Eq. (2) or (2). Let be a Hilbert space of square-integrable functions with a probability measure supported on . Assuming independent components of , the PDD expansion of function generates a hierarchical representation(rahman2008, ; rahman2009, )

(4) |

of the original performance function, in terms of an infinite number of multivariate orthonormal basis (rahman2008, ; rahman2009, ) in , where is a -dimensional multi-index; contributes the constant component; for , commits all univariate component functions representing the individual contribution to from each single input variable; for , it brings in all bivariate component functions embodying the cooperative influence of any two input variables; and for , it admits -variate component functions quantifying the interactive effects of any input variables. For most performance functions in engineering applications, a truncated version of Eq. (4) is often accurate enough by retaining, at most, the interactive effects of input variables and th order polynomials,

(5) |

where

(6) |

and

(7) |

are referred to as expansion coefficients of PDD expansion (4)
or truncated PDD approximation (5). For and ,
Eq. (5) retains interactive effects among at most
input variables ,
and th order polynomial nonlinearity in , thus leading to
the so-called -variate, th-order PDD approximation. When
and , converges to in the mean-square
sense and engenders a sequence of hierarchical and convergent approximations
of . Based on the dimensional structure and nonlinearity of a
stochastic response, the truncation parameters and can be
chosen correspondingly. The higher the values of and permit
the higher the accuracy, but also endow the computational cost of
an th-order polynomial computational complexity (rahman2008, ; rahman2009, ).
Henceforth, the -variate, th-order PDD approximation will
be simply referred to as *truncated PDD approximation* in this
paper.

### 3.2 Stochastic moment analysis

For an arbitrary random response of certain topology design , let , if it exists, denote the raw moment of of order , where . Let denote the raw moment of of order , given an -variate, th-order PDD approximation of . The analytical expressions or explicit formulae for estimating the moments using PDD approximations are described as follows. Applying the expectation operator on and , the first moment or mean (rahman10, )

(8) |

of the -variate, th-order PDD approximation is simply the constant component in Eq. (5), whereas the second moment (rahman10, )

(9) |

is expressed as the sum of squares of all expansion coefficients of . It is straightforward that the estimation of the second moment evaluated by Eq. (9) approaches the exact second moment

(10) |

of when and . The mean-square convergence
of is ensured as its component functions will contain
all required bases of the corresponding Hilbert spaces. Furthermore,
the variance of is also mean-square
convergent.

### 3.3 Reliability analysis

The RBTO problem defined in Eq. (2) requires not only stochastic moment analysis but also evaluations of the probabilistic constraints, that is, the failure probability

(11) |

of a certain topology design with respect to certain failure set . In which, the indicator function when and zero otherwise. For component-level RBTO, the failure set is often adequately characterized by a single performance function as . Whereas for a system-level RBTO, it is usually described by multiple, interdependent performance functions , leading, for example, to and for series and parallel systems, respectively. Let or or be an approximate failure set as a result of -variate, th-order PDD approximations of or of . Then the embedded MCS estimate of the failure probability is

(12) |

where is the sample size, is the th realization
of , and ,
equal to *one* when and
*zero* otherwise, is an approximation of the indicator function
.

Note that the stochastic moment analysis and reliability analysis for RTO and RBTO are quite similar to the ones in a general robust design optimization (RDO) and reliability-based design optimization (RBDO) (ren2013robust, ; rahman2014novel, ; ren2015reliability, ) except that the former is affiliated with certain topology designs . However, topology sensitivity analysis of moments and reliability is distinct from sensitivity analysis in RDO and RBDO due to the disparate topology change associated, and is elaborated in the next section.

## 4 Stochastic topology sensitivity analysis

To evaluate the topology sensitivity of a stochastic response, a new framework is proposed here which dovetails PDD and deterministic topological derivative. It relies fundamentally on the topology derivative (td5, ; td6, ; td4, ; td1, ; td2, ; td3, ; L4, ; L6, ; L7, ) of a deterministic objective function . The new method provides a closed formula and an embedded MCS formulation for the topological derivative of stochastic moments and reliability, respectively. Before presenting the new framework itself, a brief revisit on the idea of topological derivative appears to be necessary and should be convenient to those not yet familiar with the concept.

### 4.1 Topological derivative - revisit

Pioneered by Schumacher(shumacher1995topologieoptimierung, ), Sokolowski and Zochowski (sokolowski1999topological, ; sokolowski2009topological, ), and Garreau et al. (garreau2001topological, ), the topological derivative measures the change of a performance functional when an infinitesimal hole is introduced in the reference domain in which a boundary-value problem is defined. For a given reference domain , a point , and a hole with the radius of , a translated and rescaled hole can be defined by and the perforated domain is as shown in Fig. 1.

For a small , if admits the topological asymptotic expansion

(13) |

then is called the topological derivative at point and is applicable to general boundary value problems including the linear elastic system

(14) |

where is the elastic tensor, and denote Dirichlet boundary and Neumann boundary of , respectively. The topological asymptotic expansion (13) contains two performance functions and . The former is related to the reference domain and evaluated by solving (14), whereas the latter is affiliated with the perforated domain and the associated boundary value problem

(15) |

where the Neumann type condition is prescribed on , i.e., the boundary with the opposite normal vector. Comparing Eq. (14) and Eq. (15), it concludes that on and on . Moreover, it is proved that , where is reminder of higher order compared to , is the solution of the following external problem (garreau2001topological, )

(16) |

as . Solutions for various cases of isotropic elasticity are summarized in Table 1, for more details and an easy solution utilizing the Eshelby tensor, refer to Appendix A.

Both and admit a general class of performance functions. Consider the compliance of the structure as the performance functional, , which can be augmented by a Lagrange multiplier to introduce the governing equation as follows,

(17) |

by noticing being the solution of Eq. (14) in advance, where can be any kinematically admissible field that meets appropriate smoothness requirements. Similarly for the perforated domain,

(18) |

The change of compliance after perforation

(19) |

employing on and on . Integrate the second term of the above equation by parts twice and the third term one time, meanwhile applying divergence theorem,

(20) |

noticing on , on , on , and is always the normal of the current integration surface during the above derivation. Take as the displacement solution of the following adjoint problem

(21) |

and apply the solution on of the external problem for the three-dimensional case, we have

(22) |

identifying for the three-dimensional case, where is the second-order unit tensor, is the fourth-order identity tensor, and is the stress solution at of the adjoint problem. Further calculations lead to

(23) |

noticing for this case. Therefore the corresponding topological derivative has a concrete form

(24) |

where the fourth-order tensor . The evaluation of requires the stress solution at from both the original problem and the adjoint problem. In the case that , the latter becomes self-adjoint and only the solution of Eq. (14) is needed. The expressions of for various cases are summarized in Table 1.

### 4.2 Topology sensitivity of stochastic moments

Let be a response function of the linear system (14) subject to random input , which can be uncertain loads, geometry, or material properties. For a point , taking topology derivative of th moments of the response function and applying the Lebesgue dominated convergence theorem, which permits the interchange of the differential and integral operators, yields

(25) |

that is, the topology derivative is obtained from the expectation of a product comprised of the response function and its topology derivative.

For simplicity, we denote by , and construct an -variate, th-order PDD approximation as

(26) |

Replacing and of Eq. (25) by their -variate, th-order PDD approximations and , respectively, we have

(27) |

For , employing the zero mean property and orthonormal property of the PDD basis yields analytical formulation for topology sensitivity of first three moments

(28) |

(29) |

(30) |

(31) |

which requires expectations of various products of three random orthonormal polynomials (ren2013robust, ). However, if

follows classical distributions such as Gaussian, Exponential, and Uniform distribution, then the expectations are easily determined from the properties of univariate Hermite, Laguerre, and Legendre polynomials

(busbridge1948, ; ren13, ; rahman2014novel, ). For general distributions, numerical integration methods will apply.### 4.3 Topology sensitivity of reliability

Using PDD to approximate the performance function , the Monte Carlo estimate for topology sensitivity of failure probability is

(32) |

where is the sample size; is the th realization of ; and and are the indicator functions for failure domains and , respectively. The PDD approximation of the response function of the current topology design is , while at perturbed design , it is . When takes finite values, Equation (32) leads to a finite-difference approximation

(33) |

of the topology derivative for reliability. Calculating the topology derivative of reliability requires , which is simply obtained from

(34) |

requiring no additional PDD expansion or FEA. This Monte Carlo estimation involves only two PDD approximations, Eq. (5) for the response function itself and Eq. (26) for its deterministic topology derivative, both of which are generated from the same stochastic analysis. Therefore it requires no additional computational cost once the stochastic analysis is done, facilitating a novel and highly efficient sensitivity analysis approach for RBTO.

## 5 Calculation of PDD Coefficients

The expansion coefficients in Eq. (5) and Eq. (26) are defined by -dimensional integrations and etc. For large , direct numerical integration is often prohibitive, especially when FEA is involved in the Gauss point evaluation. Instead, we will use the dimension-reduction method (xu04, ; rahman2004, ; xu05, ), which entails multiple low-dimensional integrations as an effective replacement of a single -dimensional integration.

Let , which is commonly adopted as the mean of , be a reference point, and represent an -variate referential dimensional decomposition (RDD) component function of , where and . Given a positive integer , when in the above -dimensional integration is replaced with its -variate RDD approximation, the coefficients are estimated from(xu04, )

(35) |

(36) |

requiring evaluation of, at most -dimensional integrations. For each integration involved, the Gauss quadrature rule applies. For engineering problems, the evaluation of Gauss points often relies on FEA. For instance, each FEA with realized at certain gauss point supplies response function value for that Gauss point. Whereas to approximate the coefficients for the topology sensitivity or in section 4.2, each FEA provides stress results for Eq. 24 and further produces values at the corresponding Gauss point. The reduced integration is significantly more efficient than performing one -dimensional integration owing to a much fewer number of Gauss points required by the former, particularly when . On the other hand, it facilitates the calculation of coefficients approaching their exact value as . In addition, the same set of Gauss points thus the same set of FEAs will be reused for the evaluation of coefficients in Eq. (26), rendering a significantly efficient framework for stochastic topology sensitivity analysis.

## 6 Numerical Examples

In this section, two new benchmark examples are developed for the analytical or semi-analytical solution of moments and reliability and their topology sensitivities. The first one involves two random variables and renders analytical expression for stochastic quantities and their topology sensitivities of compliance. The second one contains 53 random variables to test the accuracy and efficiency of the proposed method for high dimensional problems by developing corresponding analytical and semi-analytical solutions. The third example is a three-dimensional bracket, whose topology has already been optimized, illustrating a practical application of the proposed method. In all examples, orthonormal polynomials and associated Gauss quadrature rules consistent with the probability distributions of input variables, including classical forms, if they exist, were employed. No unit for length, force, and Young’s modulus is specified in all examples for simplicity while permitting any consistent unit system for the results.

### 6.1 A round disk subject to a uniform pressure

Assuming the plane stress state, consider a round disk subject to a uniform pressure as shown in Fig. 2, where is the polar coordinate system with its origin locating at the center of the disk. The Young’s module and pressure are random variables. The Poisson’s ratio , and is deterministic.

Assume follows inverse uniform distribution on with the probability density function (PDF)

(37) |

and uniform distribution on . For this particular problem, the exact compliance is readily available, it is