Topology optimization of structures considering local material uncertainties in additive manufacturing

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 360 (2020) 112786 www.elsevier.com/locate/cma

Topology optimization of structures considering local material uncertainties in additive manufacturing Baoshou Liua , Chao Jianga , Guangyao Lia , Xiaodong Huangb,a ,∗ a

Key Laboratory of Advanced Technology for Vehicle Body Design & Manufacture, Hunan University, Changsha, 410082, China b Faculty of Science, Engineering and Technology, Swinburne University of Technology, Hawthorn, VIC 3122, Australia Received 1 August 2019; received in revised form 10 December 2019; accepted 12 December 2019 Available online xxxx

Abstract Additive manufacturing provides more freedom for the design of structures but also exhibits prominent local uncertainties of material properties, which bring potential challenges. The performance of structures should depend not only on uncertain variations of material properties but also on the spatial occurrence frequency of the extreme material properties. This paper proposes a non-probabilistic reliability-based topology optimization algorithm by considering local material uncertainties in additive manufacturing. The whole design domain is divided into several uncertainty regions (URs), whose size is proportional to the spatial occurrence frequency of extreme material properties. Within each UR, these uncertain-but-bounded variations of materials are correlated by the multi-dimensional ellipsoid model. Then, the multi-ellipsoid model for all URs is established by considering the overall material uncertainties of the structure. Thereafter, a non-probabilistic reliability-based topology optimization (NRBTO) is proposed for minimizing structural volume against displacement constraints by considering material uncertainties during additive manufacturing. Several 2D and 3D examples are presented to illustrate the effectiveness of the proposed method. Compared with solutions resulting from the deterministic topology optimization (DTO), NRBTO provides conservative designs with a larger volume fraction due to material uncertainties. When smaller URs are assigned to indicate the high occurrence frequency of extreme material properties, the NRBTO design becomes even conservative. The extreme case is equivalent to the deterministic topology optimization using the lower bound material. c 2019 Elsevier B.V. All rights reserved. ⃝ Keywords: Additive manufacturing; Material uncertainties; Multi-ellipsoid model; Non-probabilistic reliability-based topology optimization (NRBTO)

1. Introduction Topology optimization of continuum structures seeks optimal material distribution with a given design domain to meet the prescribed constraints. Since the landmark paper of Bendsøe and Kikuchi [1] in 1988, several topology optimization methods have been developed to demonstrate their validity and applicability for the stiffness optimization problem for the deterministic conditions. These methods include the homogenization method [2,3], the ∗ Corresponding author at: Faculty of Science, Engineering and Technology, Swinburne University of Technology, Hawthorn, VIC 3122,

Australia. E-mail address: [email protected] (X. Huang). https://doi.org/10.1016/j.cma.2019.112786 c 2019 Elsevier B.V. All rights reserved. 0045-7825/⃝

2

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

solid isotropic material with penalization (SIMP) method [4,5], the level-set method (LSM) [6–8], the bi-directional evolutionary structural optimization (BESO) [9–12]. Deterministic conditions, such as loading, materials, and boundary, are commonly assumed in the traditional stiffness topology optimization problem. However, the deterministic optimized design may become less efficient or even break the specified constraint when some uncertainties are considered. Topology optimization considering aleatory uncertainties includes the reliability-based topology optimization (RBTO) [13–19] and the robust topology optimization (RTO) [20–26]. In RBTO, the uncertainties can be included by adopting the constraints of the reliability indexes or performance values. While RTO aims to minimize the linear sum of the mean and standard deviation of an objective function, so as to reduce the sensitivity of the objective function with respect to given uncertain parameters. Thus, the failure probability of an optimized design can be controlled within an acceptable level. Both probabilistic-based topology optimization techniques have been widely developed for various problems, such as traditional stiffness problem [27–30], geometrically nonlinear structures [31,32], dynamic optimization problems [32,33] and compliant mechanisms [34–38]. In this paper, we will consider topology optimization of structures with material uncertainties over the whole design domain. The material uncertainties of products are often caused by the manufacturing process, from the traditional casting to the modern additive manufacturing [39,40]. On account of heterogeneous properties of materials, Silva et al. [41] integrated the second-order perturbation approach to an RTO design for minimizing compliance under the hypothesis of uncertainties. Similarly, Silva et al. [42] extended the first-order perturbation to a stress-based RTO formulation by considering material uncertainties. Kang [43] presents a robust shape and RTO design for multi-material structures by considering the uncertain graded interface. Jung and Cho [31] present an RBTO design for three-dimensional geometrically nonlinear structures in the presence of the external loads and material property uncertainties by assuming the Gaussian normal distribution. Luo [16] presents an RBTO design for stress constrained topology optimization problem by considering material uncertainties. The design results of RBTO and RTO are more reliable than those of deterministic optimization and less conservative than those calculated by safety factor based methods. However, the material uncertainties caused by the manufacture process, e.g., selected laser melting (SLM), exhibit strong local variations and the complete statistical datum are hardly measured. Meanwhile, RTO design or RBTO design is extremely sensitive to the specific probability distribution of uncertainty parameters [44–46]. Therefore, the probabilistic RTO and RBTO techniques are difficult to deal with the random material uncertainties caused by the manufacturing process. Different from the probabilistic models, the non-probabilistic models can easily handle with uncertain-butbounded material uncertainties caused by the manufacturing process without the specific probability distribution. It is assumed that the lower and upper bounds of material uncertainties can be obtained from the incomplete measuring datum. Topology optimization considering the material uncertainties has been investigated by the researchers [47–53], where the whole structure was assumed to be composed of uniform material with uncertain properties. Xu [47] developed a non-probabilistic reliability optimization (NRBTO) on natural frequency of continuum structures with uncertain-but-bounded parameters based on BESO. Wang [51–53] presented a new NRBTO method for continuum structures with uncertain-but-bounded uncertainties. However, local material uncertainties with spatial variations in additive manufacturing are more prominent and unavoidable, which are typically caused by the temperature field and particle size in the selective laser melting (SLM). To our knowledge, NRBTO by considering the spatial variations of material properties has never been studied before. This paper will propose a new NRBTO method by considering local material uncertainties in additive manufacturing. The effects of local material uncertainties for the design of structures should embody the variation ranges of material properties and the spatial occurrence frequency of the extreme material properties. When additive manufacturing achieves a high overall uniformity of the material, the extreme properties of the material may incidentally occur over a large region. On the contrary, a small region should be assumed for the high-frequency occurrence for the extreme properties of the material. In this paper, the whole design domain will be divided into 2D or 3D uncertain regions (URs), whose size is proportional to the spatial occurrence frequency of the extreme material properties. Within each UR, the uncertain-but-bounded variations of the materials are correlated by the multi-dimensional ellipsoid model [47–50,54–56]. Then, the multi-ellipsoid model [57–59] for all URs is established for material uncertainties of the whole structure. Therefore, the multi-ellipsoid model is used to integrate dependent and independent uncertainties without the specific probability distribution in all URs considering local material uncertainties. Finally, an NRBTO algorithm is proposed for minimizing structural volume against the displacement constraints.

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

3

The paper is organized as follows, Section 2 introduces the non-probabilistic uncertainty analysis for the multiellipsoid model. In Section 3, topology optimization formulations are derived. Section 4 describes the topology optimization algorithm. Some numerical examples are presented in Section 5 to demonstrate the effectiveness of the proposed method. Finally, conclusions are drawn. 2. Non-probabilistic reliability analysis using convex models 2.1. The multi-dimensional ellipsoid model The interval and ellipsoid models are two most commonly used in the field of non-probabilistic uncertain analysis. The interval model implicitly presumes that all uncertain parameters are uncorrelated and thus vary independently. However, local material properties in additive manufacturing are correlated within a certain spatial range (e.g. the proposed UR). Therefore, we will use the multi-dimensional ellipsoid model in this paper. Here, we assume that the bounds of material uncertainties within each UR can be measured and the variation of material property, yi , in point i is expressed by the interval as [ ] yi = yiL , yiR (1) where yiL and yiR represent the lower bound and upper bound of yi , respectively. Thus, the midpoint yiC and radius yiW of the interval can be calculated according to its lower and upper bounds as. ⎧ y L + yiR ⎪ ⎨ yiC = i 2 L (2) R ⎪ ⎩ y W = yi − yi i 2 The material properties over the whole structure can be expressed by the following vector. y = {y1 , . . . , yi , . . . , yn }T

(3)

where n is the total number of points inside each section of the structure. In order to overcome the ill-conditioned covariance matrix and ensure numerical precision, uncertain parameter yi in y space can be transformed to a regularized variable δi in a dimensionless vector δ space as δi =

yi − yic yiw

(4)

In the convex model theory, all possible values are assumed to be bounded in the multidimensional ellipsoid model and this model in δ space can be expressed by { } δ ∈ Ωδ = δ|δ T W δ δ ≤ 1 (5) where W δ is the correlation coefficient matrix in δ space, which is a real symmetric positive-definite matrix. These data of the single ellipsoid can be objectively extracted from outputs of instrument measurements or determined from the manufacturing tolerant specifications. Then the eigenvalue decomposition is performed for W δ . ΦT W δ Φ = Λ and ΦT Φ = I

(6)

where Φ is an orthogonal matrix composed of normalized eigenvectors and Λ is a diagonal matrix composed of the eigenvalues of W δ . I is an identity matrix. Thus, the correlation coefficient matrix W δ can be calculated according to the correlation analysis [60] as. ⎡ ⎤ b y1 y1 b y1 y2 · · · b y1 yn ⎢b y y b y y · · · b y yn ⎥ cov(yi , y j ) 2 2 2 ⎥ ⎢ 2 1 √ Wδ = ⎢ (7) ⎥ and b yi y j = √ .. .. .. .. ⎣ ⎦ cov(yi , yi ) cov(y j , y j ) . . . . b yn y1 b yn y2 · · · b yn yn where the correlation coefficient, b yi y j , is a dimensionless quantity and represents the degree of linear correlation between uncertain parameters yi and y j . In general, 0 ≤ b yi y j ≤ 1. When i = j, b yi y j = 1. If yi is uncorrelated to y j , b yi y j = 0. When uncertain parameters are uncorrelated to each other, W δ becomes an identity matrix and

4

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

Fig. 1. Two-dimensional ellipsoid models for (a) uncorrelated uncertainty parameters; (b) correlated uncertainty parameters; (c) a two-dimensional interval model.

two-dimensional ellipsoid model is shown in Fig. 1(a). Fig. 1(b) shows the ellipsoid model when two uncertainty parameters are correlated. When W δ is a null matrix, the two-dimensional ellipsoid model will be degenerated into a two-dimensional box as shown in Fig. 1(c). And the linear transformation between the dimensionless vector δ and the normalized vector q is introduced as q = Λ1/2 ΦT δ The normalized uncertainties q in space Ωq can be expressed by a hyper-sphere with unit radius. { } q ∈ Ωq = q|q T q ≤ 1

(8)

(9)

It is obvious that the normalized uncertainties set Ωq is a hyper-sphere with a unit radius. With the combination of Eqs. (4) and (8), the transformation between the original uncertainties y and the normalized uncertainties q can be expressed by )−1 ] ( ) [( y = diag y W · Λ1/2 ΦT δ · q + yC = S · q + yC (10) where diag( y W ) is the diagonalization matrix of y W in the uncertainties. Meanwhile, the transformation vector S is ( ) ( )−1 S = diag y W · Λ1/2 ΦT δ (11) 2.2. Non-probabilistic reliability index of the multi-ellipsoid model Considering two interval uncertainties in two-dimensional ellipsoid model, the ellipse model should be established from experimental samples of the uncertainties, for which y1 and y2 in the original y space are shown in Fig. 2(a) and q1 and q2 in the normalized q space are shown in Fig. 2(b). In the structural reliability analysis, the failure and safety regions are mathematically defined by the performance function, G ( y) < 0 and G ( y) > 0, respectively, as shown in Fig. 2(a). Based on the transformation in Eq. (11), the normalized performance function g (q) in the q space defines the failure region (g (q) < 0) and the safe region (g (q) > 0) as shown in Fig. 2(b). g (q) = 0 is the limit-state surface. In this paper, the design domain will be divided into N S URs and the local material variation within each UR is described by an ellipsoid set as explained above. The multi-ellipsoid model is used to tackle with independent and dependent local material uncertainties with N S URs, simultaneously. Meanwhile, the multi-ellipsoid model will avoid the extremely conservative design from the multidimensional parallelepiped model [61–63], similar to the dimensional interval model. Thus, the material uncertainties over the whole design domain can be expressed by { }T Y = y1T , . . . , y Tj , . . . , y TN S (12)

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

5

Fig. 2. (a) Original domain of the multi-dimensional ellipse model; (b) Normalized domain of the multi-dimensional ellipse model; (c) Normalized domain of the multi-ellipsoid model, where q1 and q2 are correlated and q3 is independent.

where y j is the vector of the uncertain variables in the jth UR. Based on the transformation vector S in the ellipsoid model, the normalized uncertainties over the whole design domain for the multi-ellipsoid model can be expressed by }T { (13) Q = q 1T , . . . , q Tj , . . . , q TN S where q j is the vector of the normalized uncertainties in the jth UR. Fig. 2(c) schematically illustrates a threedimensional Q space in the multi-ellipsoid model for the case of two ellipsoid sets, where one set consists of two uncertainties (q1 , q2 ) while the other one only consists of q3 . The solid cylinder in Fig. 2(c) represents the normalized space of uncertainties and the gray surface denotes the limit-state surface g ( Q) = 0. In the normalized Q space, the length of the vector Q is defined as √ √ √ ∥ Q∥ = max( q 1T q 1 , q 2T q 2 , . . . , q TN S q N S ) (14) where ∥ Q∥ represents the Euclidean norm of the multi-ellipsoid model. As shown in Fig. 2(c), the solid cylinder represents the bounds of the normalized uncertainties (∥ Q∥ = 1) and it can be proportionally expanded to the

6

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

dotted cylinder. Point C is the tangency point between the limit-state surface and the dotted cylinder. It is noted that point C has minimal distance to the origin on the basis of length measure for all points on the limit-state surface. Therefore, point C is the most-probable-point. The non-probabilistic reliability is defined as the maximum tolerance degree for uncertainties in the safe region [64], therefore, it can be quantified by the distance between the origin and point C. Thus, the quantified reliability index η can defined as the distance from the origin to point C in the normalized space of the multi-ellipsoid model ⎧ (√ { )} √ √ ⎨ η = sgn (g (0)) · min max q 1T q 1 , q 2T q 2 , . . . , q TN S q N S q (15) ⎩ s.t. g ( Q) = 0 where the sign function sgn (g (0)) aims at defining a negative reliability index when the limit-state function is negative at the origin in the normalized space. It is obvious that the greater the non-probabilistic reliability index η means that the structure will allow for the greater variations of uncertain parameters. Particularly, η = 1 implies that structure is critical for the reference parameter uncertainties and a greater η offers a specified safety margin. 2.3. Performance measure approach (PMA) The reliability analysis methods mainly include simulation methods [65–67], such as Monte Carlo simulation (MCS), Importance Sampling (IS) and Subset Simulation (SS), and the MMP-based (most-probable point) analysis including reliability index approach (RIA) and performance measure approach (PMA). As for simulation methods, the number of samples will determine the accuracy, efficiency, and robustness of the methodology. However, considering local material uncertainties, uncertainties exhibit strong local variations and the complete statistical datum are hardly measured. Therefore, it will lack enough samples to calculate an accurate joint probability density function (PDF) distribution, which makes the simulation methods to lose their accuracy and efficiency. For the MMP-based analysis, RIA may yield singularity in many RBTO applications though it is more efficient in evaluating the violated probabilistic constraint [68]. Meanwhile, sensitivity analysis of RIA is especially difficult for a min–max problem shown in Eq. (15), especially when a larger number of design variables are involved and the high nonlinearity problem may lead to numerical instability and poor convergence [57–59]. Instead, the popular performance measure approach (PMA) can be used to overcome numerical difficulties and achieve better convergence. The basic idea of the performance measure approach is to transform the min–max problem in Eq. (15) into an equivalent optimization problem as follows. { α = min g ( Q) q (16) T s.t. q j q j ≤ η2 ( j = 1, 2, . . . , N S ) where g ( Q) is the performance function and η is the specified non-probabilistic reliability index. In contrast to the direct method relying on sensitivity analysis of reliability index, the performance measure approach checks the concerned performance value at the most-portable-point in the normalized space with regard to the specified non-probabilistic reliability index. 3. Topology optimization problems and sensitivity analysis 3.1. Topology optimization problems Lightweight design of structures satisfying one or multiple displacement constraints [69] is often required in many engineering practices. It is a typical deterministic topology optimization (DTO) problem which can be expressed mathematically as the follows ∑N ⎧ ρe Ve ⎪ ⎪ min : V f = ∑e=1 ⎪ ⎪ N ρ ⎪ ⎪ e=1 Ve ⎧ ⎨ ⎪ ⎪ Ku = f DTO (17) ⎨ ⎪ ⎪ ∗ ⎪ s.t. : ⎪ u k − u k ≤ 0 (k = 1, 2, . . . , Nc ) ⎪ ⎪ ⎪ ⎪ ⎩ ⎩ ρe ∈ {0, 1} (e = 1, 2, . . . , N )

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

7

Fig. 3. Flow chart of NRBTO where the most-probable-point is searched in the inner loop and optimized topology is obtained in the outer loop.

where ρ = {ρ1 , ρ2 , . . . , ρe , . . . , ρ N } is the vector of the design variables and ρe denotes the relative density of the eth element. N is the total number of elements within the whole design domain. Ve denotes the volume of the eth element. u k and u ∗k are the nodal displacement at point k and its corresponding upper constraint. Nc is the total number of specified displacement constraints. K, u and f denote the global stiffness matrix and the displacement and force vectors, respectively. The design variable, ρe = 0 for void and 1 for solid. When uncertain material properties of a structure are considered, the displacement field of the structure becomes uncertain too. Thus, the displacement constraints are replaced by the performance constraints defined by the specified non-probabilistic reliability index, η. The corresponding non-probabilistic reliability-based topology optimization (NRBTO) problem can be stated as follows ∑N ⎧ ⎪ e=1 ρe Ve ⎪ min : V = ⎪ ∑N f ⎪ ρ ⎪ ⎪ e=1 Ve ⎧ ⎨ ⎪ ⎪ Ku = f NRBTO (18) ⎨ ⎪ ⎪ ⎪ s.t. : αk (ρ) ≥ 0 (k = 1, 2, . . . , Nc ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎩ ρe ∈ {0, 1} (e = 1, 2, . . . , N ) where αk (ρ) is the solution of a minimization problem given in Eq. (16). Meanwhile, the kth performance function is gk (ρ, Q) = u ∗k − u k (ρ, Q). The above nested problem will be divided into two sequential optimization loops as shown in Fig. 3. The inner loop searches the most-probable-point and calculates the performance constraints for the outer loop. The topology optimization in the outer loop seeks the optimized design satisfying the specified constraints. To mathematically solve the topology optimization problems in Eqs. (17) and (18), the discrete design variables should be relaxed to be continuous as 0 < ρmin ≤ ρe ≤ 1. Here, ρmin is a small value, e.g. 10−3 , to avoid the

8

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

Fig. 4. (a) A structure has meshed with square elements in the interior region and triangular elements at boundaries, and the resulting displacement at the loading point is 9.202. (b) The same structure is analyzed with the fixed square mesh (ρe = 1 for interior elements and 0.5 for boundaries elements) and the predicted displacement at the loading point is 9.203 for the ersatz material model and 10.663 for the SIMP model with p = 3.

singularity. Meanwhile, a pure 0/1 design results in a zig-zag boundary, which is less practical in applications. Ideally, the final design contains a layer of intermediate elements along the boundary of the structure to represent a smooth design. The material property of intermediate elements can be interpolated by the ersatz material model [7,8] as E (ρe ) = ρe E 0

(19)

where E 0 denotes Young’s modulus of the solid material. It is noted that we use the ersatz material model rather than the well-known solid isotropic material penalization (SIMP) model [4,5], which assumes the power-law relationship p between the Young’s modulus and the relative density of the element as E (ρe ) = ρe E 0 . Since the displacement will be used as the constraint in this paper, it is necessary to check the calculation of the displacement for a smooth design using the fixed mesh. Therefore a simple numerical test is conducted here. A two-bar smooth structure is composed of the material with Young’s modulus E 0 = 1 and Poisson’s ratio v = 0.3 under plane stress as shown in Fig. 1(a). The concentrated horizontal force F = 1 is applied at the center of the bottom as shown and the resulting displacement at the loading point is 9.202. When this structure is analyzed by using the fixed mesh as shown in Fig. 4(b). The relative density of boundary elements is assigned with 0.5 so that both models denote the same volume fraction. When the ersatz material model is used for those boundary intermediate elements, the predicted displacement at the loading point is 9.203 with only 0.01% error. However, the SIMP model with p = 3 gives the displacement at the same point, 10.663 with a 15.88% error. Therefore, we use the ersatz material model rather than the SIMP model in the paper. Consequently, the global stiffness matrix K can be expressed by K=

N ∑

ρe k0

(20)

e=1

where k0 denotes the elemental stiffness matrix for a solid element. 3.2. Sensitivity analysis In order to solve the minimum problem shown in Eq. (16) in the inner loop, the sensitivity of gk (ρ, q) with respect to the uncertainties Q should be calculated. It is noted that the design variables ρ keep constant in the inner loop, ∂u k (ρ, Q) ∂gk (ρ, Q) =− ∂Q ∂Q

(21)

Q) In order to calculate ∂u k∂(ρ, , a virtual load fk , whose components are equal to unity at the degree of freedom Q of u k and zero for others, is applied for the structure. Thus,

u k (ρ, Q) = fkT u

(22)

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

The sensitivity

∂u k (ρ, Q) ∂Q

9

of the kth displacement constraint with respect to Q can be expressed by

∂fT ∂u k (ρ, Q) ∂u ∂u = fkT + k u = fkT ∂Q ∂Q ∂Q ∂Q

(23)

The derivative of the equilibrium equation, Ku = f, results in ∂u ∂K u+K =0 ∂Q ∂Q

(24)

Substituting Eq. (24) into Eq. (23) gives ∂K ∂K ∂u k (ρ, Q) = −(fkT K−1 ) u = −ukT u ∂Q ∂Q ∂Q

(25)

where uk is the solution of the adjoint system, Kuk = fk . In order to solve the topology optimization problem, the sensitivity of the structural volume with respect to the design variable ρe can be easily calculated as ∂V Ve = ∑N ∂ρe e=1 Ve

(26)

Similar to the calculation of ρe can be expressed by

∂u k (ρ, Q) , ∂q

the sensitivity

∂u k (ρ, Q) ∂ρe

∂u k (ρ, Q) ∂K = −ukT u ∂ρe ∂ρe

of the kth displacement constraint with respect to

(27)

In the NRBTO problem, the sensitivity of the constraint, αk (ρ), with regard to the design variable ρe can be expressed by ∂αk (ρ) ∂gk (ρ, Q) ∂u k (ρ, Q) ∂K = =− = ukT u ∂ρe ∂ρe ∂ρe ∂ρe

(28)

To damping the variations of design variables, all sensitivities will be averaged with their historical values as BESO [10] ( )l−1 ( )l ( )l−1 ( )l ∂u k ∂u k ∂αk (ρ) ∂αk (ρ) + + ∂ρe ∂ρe ∂ρe ∂ρe ∂u k ∂αk (ρ) = and = (29) ∂ρe 2 ∂ρe 2 where l is the current iteration number and l > 1. It is noted that, once a solution is stably convergent, the sensitivities in the last two iterations should be equal. It means that this average scheme has no effect on the calculation of sensitivities of the convergent solution. 4. Topology optimization algorithm 4.1. Update scheme Topology optimization will iteratively update the design variables as the flow chart shown in Fig. 5. In order to consider the given constraints, the objective function in Eq. (18) can be modified by introducing the Lagrange multipliers as min : f = V f −

Nc ∑

Λk αk +

∑

( ) ς j g j − g ∗j

(30)

k=1

where Λk ≥ 0 and ς j ≥ 0 are the Lagrange multipliers for the corresponding constraint functions. g j and g ∗j denote the constraints resulting from the 0/1 constraints of the design variables. These constraints include the upper and lower bounds of the design variables, the relative density filter constraint and the floating projection constraint. Since all constraints defined in g j will be strictly enforced in the late update scheme, the last term in Eq. (30) can be ignored.

10

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

Fig. 5. Design domain and boundary conditions of the cantilever.

Consequently, the objective function can be equivalently expressed by the bounded Lagrange multipliers ( ) Nc Nc ∑ ∑ min : fˆ = 1 − λk V f − λk αk k=1

(31)

k=1

where 0 ≤ λk < 1 is the bounded Lagrangian multiplier for the constraint αk (ρ) ≥ 0. The sensitivity of the modified objective function with respect to the design variables can be expressed by ( ) Nc Nc ∑ ∑ ∂Vf ∂ fˆ ∂αk = 1− λk − λk (32) ∂ρe ∂ρe ∂ρe k=1 k=1 Thus, the optimality condition can be expressed by ) ( Nc Nc ∑ ∑ ∂αk /∂ρe λk − λk A ≡ 1− =0 ∂ V f /∂ρe k=1 k=1

(33)

To achieve an optimal solution, the design variable ρe can be firstly updated according to the following equation. ρe = ρel−1 (1 − A)

(34)

The negative sign before A denotes the search direction for the minimization problem. In order to enforce the upper and lower bounds of the design variables, the design variables will be modified by ⎧ ⎨ρmin ρˆe = ρe ⎩ 1

when ρe ≤ ρmin when ρmin < ρe < 1 when ρe ≥ 1

(35)

In order to avoid numerical instability such as checkerboard pattern and mesh-dependent behavior [70], the following filter constraint is adopted for further modifying the design variable as ∑N wen ρˆn (36a) ρ˜e = ∑n=1 N n=1 wen where the subscript n denotes the number of the element and the wen is the linear weight factor as { r − ren i f ren < rmin wen = min 0 i f ren ≥ rmin

(36b)

where the ren is the distance between the center of the eth element and the nth element. rmin is the specified filter radius. Since the ersatz material model with the linear interpolation scheme is used, the resulting design may contain a lot of intermediate elements. Therefore, the following floating projection constraint is also adopted to simulate the

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

11

original 0/1 constraint in the problem. [ ] tanh (β · th) + tanh β · (ρ˜e − th) l (37) ρe = tanh (β · th) + tanh [β · (1 − th)] where ∑ N βl controls ∑ N the steepness of the projection and th is the threshold value which is determined according to ρ = e=1 e e=1 ρ˜e . It can be seen that the threshold, th, varies during the whole optimization procedure. We should point out that the above projection is used as one of the constraints as defined in Eq. (30), which is fundamentally different from the one used in the traditional three-phase SIMP method, where the projection converts variables from the filter phase to the physical phase. In the optimization, β initially assigns a small value (e.g. 10-6 in this paper) until the solution is convergent. Then, β increases with ∆β (e.g., 0.5) until another convergent solution is achieved. β is stopped to increase when the solution achieves a clear topology. In each iteration, the moving limit γ may also be applied as ρmin ≤ ρel−1 (1 − γ ) ≤ ρel ≤ ρel−1 (1 + γ ) ≤ 1. In this paper, γ is set as 0.02 to limit the variation of the design variables in each iteration. 4.2. Calculation of Lagrange multipliers Since the Lagrange multipliers are bounded with 0 ≤ λk < 1, their values can be efficiently determined according to the corresponding constraints. For example, αk in the next iteration can be approximately estimated by αkl = αkl−1 +

N ∑ ) ∂αkl ( l ρe − ρel−1 l ∂ρe e=1

(38)

where αkl−1 is the unbalanced αk in the last iteration. It can be further expressed by αkl = αkl−1 +

N ∑ ) ∂αkl ∂ρel ∂ ρ˜e ∂ ρˆe ( ρe − ρel−1 l ∂ρe ∂ ρ˜e ∂ ρˆe ∂ρe e=1

(39)

Note that ρe is an implicit function of the Lagrange multipliers as shown in Eqs. (33) and (34). In order to ensure the constraint αk ≥ 0, the Lagrange multiplier λk can be updated according to the Newton–Raphson procedure as αk m−1 λm − (40) k = λk ∂αk /∂λk where m is the inner iteration number initial value λ0k will be set as 0.5 in this paper. The Lagrange ⏐ mand the ⏐ m−1 multipliers are obtained when ∆λk = ⏐λk − λk ⏐ is small enough, e.g., 10−6 . 4.3. Representation of an optimized topology In order to represent the optimized design with a smooth boundary, elemental densities, ρel , can be linearly interpolated into the whole design domain, ρ l (x, y). Then, a level-set function, φ l (x, y) = ρ l (x, y) − th, is constructed for representing the structural boundary by the zero level-set ⎧ l ⎨φ (x, y) > 0 for solid region φ l (x, y) = 0 for boundary (41) ⎩ l φ (x, y) < 0 for void region Obviously, the represented smooth topology is accurate when the optimized solution only contains boundary intermediate elements. Therefore, a topological error is estimated by the volume of the intermediate elements below the threshold, th, dividing by the total volume of the structure. When the topological error is smaller than a certain value, such as 3% in this paper, the solution can be assumed to be very close to a smooth design. Then, the steepness of the floating projection, β, will then keep as the constant until the solution is strictly convergent. 5. Results and discussion In this section, some numerical examples are presented to demonstrate the effectiveness of the proposed NRBTO formulation. As explained above, the whole design domain will be divided by several URs. For simplicity, the

12

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

Fig. 6. DTO design with a uniform material by the ersatz material model and SIMP model (left: expressed by elements; right: expressed by the level-set function). (a) The ersatz material model and V f = 47.08%; (b) SIMP model ( p = 3) and V f = 50.25%; (c) Distribution of vertical deformations of the current DTO design (left) and the SIMP DTO design (right).

Young’s modulus of the material, E 0 , in all URs varies between [8 MPa, 12 MPa] and Poisson’s ratio is v = 0.3. The correlation coefficient matrix W δ is assumed as the identity matrix except for being clearly specified. The filter radius for all examples is rmin = 2 mm. Generally, a large size of URs can be recommended for additive manufacturing, which can achieve the uniform material over a large region. The size of URs depends on the printing parameters such as particle size, laser power and scan speed. In the following examples, various size of URs are arbitrarily selected to illustrate the effects of various local material uncertainties on the design of structures. 5.1. A cantilevers with a single displacement constraint The first example considers the design problem for a cantilever shown in Fig. 5 with dimensions 240 mm in length and 80 mm in height. The design domain is discretized into 240 × 80 four-node plane-stress elements. A concentrate force, F = 0.1 N, is acted at the right-bottom corner. The maximum downward displacement at the loading point is specified as u ∗ = 1.8 mm. The objective is to find the optimal topology that minimizes the total material volume under this displacement constraint. The non-probabilistic reliability indexes are specified as η = 1.0 or 2.0 in the NRBTO problem. Generally, η = 1.0 will just satisfy the reliability required in practical engineering. A greater η implies that an optimized structure can tolerate a greater variation of material properties. In this example, the whole design domain is divided into 12 × 4, 24 × 8 and 30 × 10 URs, respectively, to represent the different occurrence frequency of the extreme material properties caused by the additive manufacturing. When we assume that the structure is composed of a uniform material, the optimized design resulting from DTO is shown in Fig. 6(a) with the volume fraction, 47.08%. In order to compare the ersatz material model with the SIMP model, the DTO design with u ∗ = 1.8 mm is also conducted by using SIMP with p = 3 and its optimized result is shown in Fig. 6(b). The final volume fraction becomes 50.25%, which is much higher than that of the current solution. Meanwhile, the smooth boundaries of both designs are extracted and simulated the commercial COMSOL Multiphysics software as shown in Fig. 6(c). It indicates that the vertical displacement at the loading point for the SIMP design is 1.717 mm, which is much lower than the constraint value, 1.8 mm. The corresponding displacement for the current design using the ersatz material model is 1.795 mm, which is very close to the specified constraint value.

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

13

Fig. 7. Optimal topologies (left: expressed by elements; right: expressed by the level-set function) (a) NRBTO design with η = 1.0 and 12 × 4 URs, V f = 48.13%; (b) NRBTO design with η = 2.0 and 12 × 4 URs, V f = 49.38%.

Fig. 8. Iteration histories of the volume fraction, topologies and displacement constraint for DTO and NRBTO with η = 1.0 and 12 × 4 URs.

Considering the material uncertainties with 12 × 4 URs, the optimized design resulting from the NRBTO with η = 1.0 is presented in Fig. 7(a) with the volume fraction, 48.13%. The increase of the volume fraction comes from the effect of the material uncertainties. For safety reason, the optimized design resulting from the NRBTO with η = 2.0 is given in Fig. 7(b) with the volume fraction, 49.38%, which is even less than that of the DTO design using SIMP. Fig. 8 shows the iteration histories of the volume fraction, topologies and displacement constraint for DTO and NRBTO with η = 1.0. It indicates that the displacement constraints for both cases are stably satisfied at the last stage of the optimization. Meanwhile, the volume fractions are stably convergent to an optimal value and the final designs have clear topologies. With consideration for the serve variation of material properties in additive manufacturing, the design domain can be divided into 24 × 8 URs. The NRBTO results are shown in Fig. 9 and the final volume fractions are 48.86% and 50.93% for η = 1.0 and 2.0, respectively. When the design domain is divided into 30 × 10 URs, the NRBTO

14

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

Fig. 9. Optimized design of NRBTO with 24 × 8 URs (left: expressed by elements; right: expressed by the level-set function): (a) η = 1.0, V f = 48.86%; (b) η = 2.0, V f = 50.93%.

Fig. 10. Optimized designs of NRBTO with 30 × 10 URs (left: expressed by elements; right: expressed by the level-set function): (a) η = 1.0, V f = 49.21%; (b) η = 2.0, V f = 51.81%.

Fig. 11. Optimized design (left: expressed by elements; right: expressed by the level-set function) of NRBTO with η = 1.0 and 240 × 80 URs, V f = 67.53%, which is equivalent to the DTO design using the weakest material, E 0 = 8 MPa.

designs shown in Fig. 10 have the final volume fractions, 49.21% and 51.81% for η = 1.0 and 2.0, respectively. It can be seen that increasing the reliability index η, or the total number of URs leads to a more conservative design with a larger volume fraction. This is understandable because the increase of the reliability index denotes the increase of the safety factor, and the increase of the total number of URs indicates a higher occurrence frequency of extreme material properties. When the design domain is divided into 240 × 80 URs, it means that each UR (or ellipsoid set) only contains one element. The NRBTO design for η = 1.0 is shown in Fig. 11 with the volume fraction, V f = 67.53%. Actually, this design is equivalent to the DTO design using the weakest material, E 0 = 8 MPa.

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

15

Fig. 12. Design domain and boundary conditions of the long beam.

Fig. 13. Optimal results (left: expressed by elements; right: expressed by the level-set function): (a) DTO design, V f = 46.63%; (b) NRBTO design with η = 1.0 and 12 × 3 URs, V f = 47.37%; (c) NRBTO design with η = 2.0 and 12 × 3 URs, V f = 48.74%; (d) NRBTO design with η = 1.0 and 24 × 6 URs, V f = 48.41%; (e) NRBTO design with η = 2.0 and 24 × 6 URs, V f = 49.77%.

Therefore, it is the most conservative case for η = 1.0. Compared with Figs. 7, 9, 10, the change of topologies in Fig. 11 is obvious. Although topologies of NRBTO designs change substantially with the consideration of a smaller size for URs, the change pattern in topologies is not different from those of the RTO design in reference [41], which have complex topologies with many thin members while using a small correlation length. This difference may be attributed to the difference between the non-probabilistic model and the probabilistic model under the specific probability distribution.

16

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

Fig. 14. NRBTO results with the correlation coefficient 0.6 (left: expressed by elements; right: expressed by the level-set function): (a) NRBTO design with η = 1.0 and 12 × 3 URs, V f = 47.37%; (b) NRBTO design with η = 2.0 and 12 × 3 URs, V f = 49.05%.

Fig. 15. Design domain and boundary conditions of the simply-supported beam.

5.2. A long cantilever with two displacement constraints As illustrated in Fig. 12, the design domain of a long cantilever has dimensions of 240 mm × 60 mm and is discretized into 240 × 60 four-node plane-stress elements. Two concentrated forces (F1 = 0.1 N and F2 = 0.6 N) are applied on the point A and B at the bottom edge. The vertical displacement constraints are enforced at the loading points, which are u A ≤ u ∗A = 3.8 mm and u B ≤ u ∗B = 4.2 mm, respectively. The design domain is divided with 12 × 3 and 24 × 6 URs, respectively, which denote the different severity of material variations. The proposed NRBTO problem will be solved under the specific reliability indexes η = 1.0 and 2.0, respectively. For comparison, the DTO design is shown in Fig. 13(a) with the final volume fraction, 46.63. The NRBTO designs for 12 × 3 URs are shown in Fig. 13(b) and (c) with the volume fractions, 47.37% and 48.74%, respectively. The NRBTO designs for 24 × 6 URs are presented in Fig. 13(d) and (e) with the volume fractions, 48.41% and 49.77%, respectively. Compared to the DTO design, the NRBTO solutions have different topologies with larger volume fractions. Similar to the previous example, increasing the reliability index η, or the total number of URs leads to a more conservative design with a larger volume fraction. As described in Section 4, all Lagrange multipliers are limited as 0 ≤ λk < 1. λk = 0 indicates that the corresponding constraint is inactive. Therefore, the final Lagrange multipliers are shown in Table 1, which well indicates both displacement constraints are active. It should be noted that the correlation coefficients in W δ may also affect the optimized solution. For instance, b yi y j = 0.6(when i ̸= j) and b yi y j = 1(when i = j) are assumed. The NRBTO designs for 12 × 3 URs are shown in Fig. 14(a) for the specified reliability index η = 1.0 and Fig. 14(b) for η = 2.0, respectively. The corresponding volume fractions are 47.47% and 49.05%. Compared to the solutions using the identity matrix shown in Fig. 13(b) and (c), numerical examples indicate that the increasing correlation among design variables causes the final solution having a larger volume fraction and more thin members.

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

17

Fig. 16. Optimal results (left: expressed by elements; right: expressed by the level-set function): (a) DTO design, V f = 39.44%; (b) NRBTO design with η = 1.0 and 15 × 5 URs, V f = 40.41%; (c) NRBTO design with η = 2.0 and 15 × 5 URs, V f = 42.65%; (d) NRBTO design with η = 1.0 and 27 × 9 URs, V f = 40.94%; (e) NRBTO design with η = 2.0 and 27 × 9 URs, V f = 43.93%. Table 1 Lagrange multipliers of displacement constraints for the final design. η DTO design

λA

λB

0.1753

0.0045

η = 1.0

12 × 3 URs 24 × 6 URs

0.1814 0.1798

0.0052 0.0068

η = 2.0

12 × 3 URs 24 × 6 URs

0.1966 0.1896

0.0065 0.0068

5.3. A simply-supported beam with three displacement constraints As shown in Fig. 15, a simply-supported beam undertakes three external loads at its bottom. The design domain with dimensions of 270 mm×90 mm is discretized into 270 × 90 four-node plane-stress elements. The concentrated forces, F1 = F2 = F2 = 0.5 N, are applied to the 1/3, 1/2 and 2/3 of the bottom edge (A, B and C) as shown in Fig. 15. The vertical displacements at A, B and C are restricted as u A = u C ≤ u ∗1 = 1.4 mm and u B ≤ u ∗2 = 1.6 mm.

18

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786 Table 2 Lagrange multipliers of displacement constraints for the final design. η

λA

DTO design

λB

λC

0.1351

0.0608

0.1351

η = 1.0

15 × 5 URs 27 × 9 URs

0.1432 0.1438

0.0617 0.062

0.1432 0.1438

η = 2.0

15 × 5 URs 27 × 9 URs

0.15 0.1507

0.0622 0.0641

0.15 0.1507

Fig. 17. Design domain and boundary conditions of the 3D cantilever.

The proposed NRBTO will be solved under the specific reliability indexes η = 1.0 and 2.0, respectively. Two types of the URs, 15 × 5 and 27 × 9, are used. Fig. 16 shows the optimized results for DTO and NRBTO using different reliability index and URs. It clearly indicates that NRBTO by considering material uncertainties provides conservative designs with a larger volume fraction, which changes the optimized topology as well. The Lagrange multipliers for final designs are listed in Table 2. It shows that all of those three displacement constraints are active. 5.4. A 3D cantilever with a single displacement constraint The developed NRBTO can be equally extended for 3D cases such as a 3D cantilever shown in Fig. 17. The design domain is discretized by 60 × 40 × 10 brick elements with size 1 mm×1 mm×1 mm. A vertical concentrated force, F1 = 5 N is applied at the center of the right surface. The maximum allowable vertical displacement at the load is u ∗A = 3 mm. The optimized design resulting from DTO is shown in Fig. 18(a), which is expressed by finite elements (left) or the level-set function (right). The final volume fraction is V f = 12.71%. To account for the uneven severity of the material, the whole design domain is divided into 6 × 4 × 1 volume URs. Fig. 18(b) and (c) show the optimized results for NRBTO using different reliability indexes, η = 1.0 and 2.0. Similar to 2D cases, NRBTO provides conservative designs by considering the uncertainties of the material, where the final volume fraction increases from 12.71% to 13.22% for η = 1.0 and 13.86% for η = 2.0. Fig. 18(d) shows the optimized design for the extreme case of η = 1.0 with 60 × 40 × 10 URs. The optimized design is equivalent to the DTO design using the weakest material, E 0 = 8 MPa. The final volume fraction increases to 16.36%. 6. Conclusion This paper aims to create a reliable topological design by considering local material uncertainties, which often occur in additive manufacturing. These uncertainties are qualified by the variation range and occurrence frequency

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

19

Fig. 18. Optimal results (left: expressed by elements; right: expressed by the level-set function): (a) DTO design, V f = 12.71%; (b) NRBTO design with η = 1.0 and 6 × 4 × 1 URs, V f = 13.22%; (c) NRBTO design with η = 2.0 and 6 × 4 × 1 URs, V f = 13.86%; (d) NRBTO design with η = 1.0 and 60 × 40 × 10 URs, V f = 16.36%.

20

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

of extreme material properties. The design domain is divided into 2D or 3D uncertainty regions (URs), whose size is proportional to the occurrence frequency of the extreme material properties. The multi-ellipsoid model is employed to integrate these uncertain-but-bounded material uncertainties in all URs. Then, the non-probabilistic reliabilitybased topology optimization (NRBTO) problem is formulated as minimizing structural volume against single or multiple displacement constraints. In the proposed NRBTO algorithm, the ersatz material model is used to avoid the possible overestimation of the displacements and the floating projection technique pushes the design variables to solid or void for a clear topology. Some 2D and 3D numerical examples have demonstrated the applicability and the effectiveness of the proposed method by considering the local material uncertainties. Compared to the result of the deterministic design, a larger volume fraction is achieved to ensure a more reliable design due to the consideration of local material uncertainties. The extreme case is equivalent to the deterministic topology optimization using the lower bound material. The future research on local stress constrained topology optimization with local material uncertainties is highly recommended. Acknowledgments Authors wish to acknowledge the financial support from the Australian Research Council (FT130101094), Hunan Provincial Innovation Foundation for Postgraduate (CX20190278) and the China Scholarship Council. References [1] M.P. Bendsøe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Engrg. 71 (1988) 197–224. [2] K. Suzuki, N. Kikuchi, A homogenization method for shape and topology optimization, Comput. Methods Appl. Mech. Engrg. 93 (1991) 291–318. [3] S. Nishiwaki, M.I. Frecker, S. Min, N. Kikuchi, Topology optimization of compliant mechanisms using the homogenization method, Int. J. Numer. Methods Eng. 42 (1998) 535–559. [4] M.P. Bendsøe, Optimal shape design as a material distribution problem, Struct. Optim. 1 (1989) 193–202. [5] M.P. Bendsøe, O. Sigmund, Topology Optimization: Theory, Methods and Applications, Springer-Verlag, Berlin, 2003. [6] J.A. Sethian, A. Wiegmann, Structural boundary design via level set and immersed interface methods, J. Comput. Phys. 163 (2000) 489–528. [7] M.Y. Wang, X. Wang, D. Guo, A level set method for structural topology optimization, Comput. Methods Appl. Mech. Engrg. 192 (2003) 227–246. [8] G. Allaire, F. Jouve, A.M. Toader, Structural optimization using sensitivity analysis and a level-set method, J. Comput. Phys. 1 (2004) 363–393. [9] O.M. Querin, G.P. Steven, Y.M. Xie, Evolutionary structural optimization (ESO) using a bi-directional algorithm, Eng. Comput. 15 (1998) 1031–1048. [10] X. Huang, Y.M. Xie, Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method, Finite Elem. Anal. Des. 43 (14) (2007) 1039–1049. [11] X. Huang, Y.M. Xie, Evolutionary Topology Optimization of Continuum Structures: Methods and Applications, John Wiley & Sons, Chichester, 2010. [12] X. Huang, Y.H. Zuo, Y.M. Xie, Evolutionary topological optimization of vibrating continuum structures for natural frequencies, Comput. Struct. 88 (2010) 357–364. [13] G.W. Melenka, J.S. Schofield, M.R. Dawson, J.P. Carey, Evaluation of dimensional accuracy and material properties of the makerbot 3d desktop printer, Rapid Prototyp. J. 21 (2013) 618–627. [14] J.S. Yu, S.R. Kim, J.Y. Park, S.Y. Han, Reliability-based topology optimization based on bidirectional evolutionary structural optimization, Proc. World Acad. Sci. Eng. Technol. 19 (2010) 529–538. [15] C. Kim, S. Wang, K.R. Bae, H. Moon, K.K. Choi, Reliability-based topology optimization with uncertainties, J. Mech. Sci. Technol. 20 (2006) 494–504. [16] Y. Luo, M. Zhou, M.Y. Wang, Z. Deng, Reliability-based topology optimization for continuum structures with local failure constraints, Comput. Struct. 143 (2014) 73–84. [17] M. Jalalpour, J.K. Guest, T. Igusa, Reliability-based topology optimization of trusses with stochastic stiffness, Struct. Saf. 43 (2013) 41–49. [18] M. Jalalpour, M. Tootkaboni, An efficient approach to reliability-based topology optimization for continua under material uncertainty, Struct. Multidiscip. Optim. 53 (2016) 759–772. [19] S. Bobby, A. Suksuwan, S.M.J. Spence, A. Kareem, Reliability-based topology optimization of uncertain building systems subject to stochastic excitation, Struct. Saf. 66 (2017) 1–16. [20] J. Zhao, C. Wang, Robust structural topology optimization under random field loading uncertainty, Struct. Multidiscip. Optim. 50 (2014) 517–522. [21] J.K. Guest, T. Igusa, Structural optimization under uncertain loads and nodal locations, Comput. Methods Appl. Mech. Engrg. 198 (2008) 116–124.

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

21

[22] A. Asadpoure, M. Tookkaboni, J.K. Guest, Robust topology optimization of structures with uncertainties in stiffness – Application to truss structures, Comput. Struct. 89 (2011) 1131–1141. [23] X. Guo, W. Zhang, L. Zhang, Robust structural topology optimization considering boundary uncertainties, Comput. Methods Appl. Mech. Engrg. 253 (2013) 356–368. [24] J. Wu, J. Gao, Z. Luo, T. Brown, Robust topology optimization for structures under interval uncertainty, Adv. Eng. Softw. 99 (2016) 36–48. [25] J.L. Wu, Z. Luo, H. Li, N. Zhang, Level-set topology optimization for mechanical metamaterials under hybrid uncertainties, Comput. Methods Appl. Mech. Engrg. 319 (2017) 414–441. [26] J.T. Liu, H.C. Gea, Robust topology optimization under multiple independent unknown-but-bounded loads, Comput. Methods Appl. Mech. Engrg. 329 (2018) 464–479. [27] G. Kharmanda, N. Olhoff, A. Mohamed, M. Lemaire, Reliability-based topology optimization, Struct. Multidiscip. Optim. 26 (2004) 295–307. [28] S.K. Chen, W. Chen, Sanghoon Lee, Level set based robust shape and topology optimization under random field uncertainties, Struct. Multidiscip. Optim. 41 (2010) 507–524. [29] J.S. Yu, S.R. Kim, J.Y. Park, S.Y. Han, Reliability-based topology optimization based on bidirectional evolutionary structural optimization, Proc. World Acad. Sci. Eng. Technol. 19 (2010) 529–538. [30] Y.S. Eom, K.S. Yoo, J.Y. Park, S.Y. Han, Reliability-based topology optimization using a standard response surface method for three-dimensional structures, Struct. Multidiscip. Optim. 43 (2011) 287–295. [31] H.S. Jung, S. Cho, Reliability-based topology optimization of geometrically nonlinear structures with loading and material uncertainties, Finite Elem. Anal. Des. 41 (2004) 311–331. [32] M.A. Valdebenito, G.I. Schuëller, Efficient strategies for reliability-based optimization involving non-linear, dynamical structures, Comput. Struct. 89 (2011) 1797–1811. [33] Z.C. He, Y. Wu, Eric Li, Topology optimization of structure for dynamic properties considering hybrid uncertain parameters, Struct. Multidiscip. Optim. 57 (2017) 1–14. [34] M. Schevenels, B.S. Lazarov, O. Sigmund, Robust topology optimization accounting for spatially varying manufacturing errors, Comput. Methods Appl. Mech. Engrg. 200 (2011) 3613–3627. [35] B.S. Lazarov, M. Schevenels, O. Sigmund, Topology optimization considering material and geometric uncertainties using stochastic collocation methods, Struct. Multidiscip. Optim. 46 (2012) 597–612. [36] M. Jansen, G. Lombaert, M. Diehl, B.S. Lazarov, O. Sigmund, M. Schevenels, Robust topology optimization accounting for misplacement of material, Struct. Multidiscip. Optim. 47 (2013) 317–333. [37] K. Maute, D.M. Fragopol, Reliability-based design of MEMS mechanisms by topology optimization, Comput. Struct. 81 (2003) 813–824. [38] K.H. Cho, J.Y. Park, M.G. lm, S.Y. Han, Reliability-based topology optimization of electro-thermal-compliant mechanisms with a new material mixing method, Int. J. Precis. Eng. Manuf. 13 (2012) 693–699. [39] L. Stabile, M. Scungio, G. Buonanno, F. Arpino, G. Ficco, Airborne particle emission of a commercial 3D printer: the effect of filament material and printing temperature, Indoor Air (2016). [40] G.W. Melenka, J.S. Schofield, M.R. Dawson, J.P. Carey, Evaluation of dimensional accuracy and material properties of the makerbot 3d desktop printer, Rapid Prototyp. J. 21 (2013) 618–627. [41] G.A. da Silva, E.L. Cardoso, Topology optimization of continuum structures subjected to uncertainties in material properties, Internat. J. Numer. Methods Engrg. 106 (2016) 192–212. [42] G.A. da Silva, E.L. Cardoso, Stress-based topology optimization of continuum structures under uncertainties, Comput. Methods Appl. Mech. Engrg. 313 (2017) 647–672. [43] Z. Kang, C. Wu, Y. Luo, M. Li, Robust topology optimization of multi-material structures considering uncertain graded interface, Compos. Struct. 208 (2019) 395–406. [44] I. Elishakoff, Essay on uncertainties in elastic and viscoelastic structures: from A.M Freudenthal’s criticisms to modern convex modelling, Comput. Struct. 17 (1995) 871–895. [45] Y. Ben-Haim, I. Elishakoff, Convex Models of Uncertainty in Applied Mechanics, Elsever, Amsterdam, 2013. [46] Y. Ben-Haim, A non-probabilistic concept of reliability, Struct. Saf. 14 (1994) 227–245. [47] B. Xu, L. Zhao, Y.M. Xie, J. Jiang, Topology optimization of continuum structures with uncertain-but-bounded parameters for maximum non-probabilistic reliability of frequency requirement, J. Vib. Control 23 (2015) 2557–2566. [48] M. Xia, D. Gu, G. Yu, D. Dai, H. Chen, Q. Shi, Influence of hatch spacing on heat and mass transfer, thermodynamics and laser processability during additive manufacturing of Inconel 718 alloy, Int. J. Mach. Tools Manuf. 109 (2016) 147–157. [49] Y. Ben-Haim, Convex models of uncertainty in radial pulse buckling of shells, J. Appl. Mech. 60 (3) (1993) 683–688. [50] H.E. Lindberg, Convex models for uncertain imperfection control in multimode dynamic buckling, J. Appl. Mech. 59 (4) (1992) 937–945. [51] L. Wang, D. Liu, Y. Yang, X. Wang, Z. Qiu, A novel method of non-probabilistic reliability-based topology optimization corresponding to continuum structures with unknown but bounded uncertainties, Comput. Methods Appl. Mech. Engrg. 326 (2017) 573–595. [52] L. Wang, J. Liang, D. Liu, W. Chen, A novel reliability-based topology optimization framework for the concurrent design of solid and truss-like material structures with unknown-but-bounded uncertainties, Int. J. Numer. Methods Eng. 119 (2019) 239–260. [53] L. Wang, D. Liu, Y. Yang, J. Hu, Novel methodology of Non-probabilistic Reliability-based Topology Optimization (NRBTO) for multi-material layout design via interval and convex mixed uncertainties, Comput. Methods Appl. Mech. Engrg. 346 (2019) 550–573. [54] Y. Luo, Z. Kang, A. Li, Structural reliability assessment based on probability and convex set mixed model, Comput. Struct. 87 (2009) 1408–1415.

22

B. Liu, C. Jiang, G. Li et al. / Computer Methods in Applied Mechanics and Engineering 360 (2020) 112786

[55] Z. Kang, Y. Luo, Reliability-based structural optimization with probability and convex set hybrid models, Struct. Multidiscip. Optim. 42 (2010) 89–102. [56] L. Wang, J. Liang, D. Wu, A non-probabilistic reliability-based topology optimization (NRBTO) method of continuum structures with convex uncertainties, Struct. Multidiscip. Optim. 58 (2018) 2601–2620. [57] Y. Luo, Z. Kang, Z. Luo, A. Li, Continuum topology optimization with non-probabilistic reliability constraints based on multi-ellipsoid convex model, Struct. Multidiscip. Optim. 39 (2009) 297–310. [58] Y. Luo, Z. Kang, Z. Yue, Maximal stiffness design of two-material structures by topology optimization with non-probabilistic reliability, AIAA J. 50 (50) (2012) 1993–2003. [59] Z. Kang, Y. Luo, Non-probabilistic reliability-based topology optimization of geometrically nonlinear structures using convex models, Comput. Methods Appl. Mech. Engrg. 198 (2009) 3228–3238. [60] C. Jiang, X. Han, G.Y. Lu, J. Liu, Z. Zhang, Y.C. Bai, Correlation analysis of non-probabilistic convex model and corresponding structural reliability technique, Comput. Methods Appl. Mech. Engrg. 200 (2011) 2528–2546. [61] C. Jiang, Q.F. Zhang, X. Han, J. Liu, D.A. Hu, Multidimensional parallelepiped model-a new type of non-probabilistic convex model for structural uncertainty analysis, Internat. J. Numer. Methods Engrg. 103 (2015) 31–59. [62] B.Y. Ni, C. Jiang, X. Han, An improved multidimensional parallelepiped non-probabilistic model for structural uncertainty analysis, Appl. Math. Model. 40 (2016) 4727–4745. [63] J. Zheng, Z. Luo, C. Jiang, B. Ni, J. Wu, Non-probabilistic reliability-based topology optimization with multidimensional parallelepiped convex model, Struct. Multidiscip. Optim. 1 (2017) 1–17. [64] I. Elishakoff, Discussion on: a non-probabilistic concept of reliability, Struct. Saf. 17 (1995) 195–199. [65] S.K. Au, J.L. Beck, A new adaptive importance sampling scheme for reliability calculation, 21 (1999) 135–158. [66] B. Echard, N. Gayton, M. Lemaire, AK-MCS: An active learning reliability method combining Kriging and Monte Carlo simulation, Struct. Saf. 33 (2011) 145–154. [67] M. Norouzi, E. Nikolaidis, Integrating subset simulation with probabilistic re-analysis to estimate reliability of dynamic systems, Struct. Multidiscip. Optim. 48 (2013) 533–548. [68] J. Tu, K.K. Choi, Y.H. Park, A new study on reliability-based design optimization, J. Mech. Des. 121 (1999) 557–564. [69] L. Yin, W. Yang, Optimality criteria method for topology optimization under multiple constraints, Comput. Struct. 79 (2001) 1839–1850. [70] O. Sigmund, J. Petersson, Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboard, mesh-dependencies and local minina, Struct. Optim. 16 (1998) 68–75.