Bi-directional evolutionary topology optimization of geometrically nonlinear continuum structures with stress constraints

Journal Pre-proof

Bi-directional evolutionary topology optimization of geometrically nonlinear continuum structures with stress constraints Bin Xu , Yongsheng Han , Lei Zhao PII: DOI: Reference:

S0307-904X(19)30748-6 https://doi.org/10.1016/j.apm.2019.12.009 APM 13195

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

24 June 2019 1 December 2019 5 December 2019

Please cite this article as: Bin Xu , Yongsheng Han , Lei Zhao , Bi-directional evolutionary topology optimization of geometrically nonlinear continuum structures with stress constraints, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.12.009

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Elsevier Inc. All rights reserved.

Highlights 

The stress-constrained topology optimization model of geometrically nonlinear structure is first proposed.



The adopted p-norm global stress measure can well restrain the stress level and solve the stress concentration problem.



The compliance increases and the stress distribution becomes smoother with the strict stress constraint.



The maximal stress and compliance of geometry nonlinear stiffness design are smaller than those for linear design.

1

Revised Manuscript prepared for Applied Mathematical Modelling June 2019

Bi-directional evolutionary topology optimization of geometrically nonlinear continuum structures with stress constraints Bin Xu*, Yongsheng Han, and Lei Zhao Institute of Structural Health Monitoring and Control, School of Mechanics, Civil Engineering & Architecture, Northwestern Polytechnical University, Xi’an 710072, China

*Corresponding author Dr. Bin Xu Email: [email protected]

2

Abstract: This paper proposes a design method to maximize the stiffness of geometrically nonlinear continuum structures subject to volume fraction and maximum von Mises stress constraints. An extended bi-directional evolutionary structural optimization (BESO) method is adopted in this paper. BESO method based on discrete variables can effectively avoid the well-known singularity problem in density-based methods with low density elements. The maximum von Mises stress is approximated by the p-norm global stress. By introducing one Lagrange multiplier, the objective of the traditional stiffness design is augmented with p-norm stress. The stiffness and p-norm stress are considered simultaneously by the Lagrange multiplier method. A heuristic method for determining the Lagrange multiplier is proposed in order to effectively constrain the structural maximum von Mises stress. The sensitivity information for designing variable updates is derived in detail by adjoint method. As for the highly nonlinear stress behavior, the updated scheme takes advantages from two filters respectively of the sensitivity and topology variables to improve convergence. Moreover, the filtered sensitivity numbers are combined with their historical sensitivity information to further stabilize the optimization process. The effectiveness of the proposed method is demonstrated by several benchmark design problems.

Keywords: Topology optimization; Geometrically nonlinearity; Stress constraints; BESO method; Sensitivity analysis; Adjoint method

3

1. Introduction Topology optimization, as a powerful design method, enable designers to find the optimal structural topology or material layout for a required structural or material performance within a given design domain for specified objectives, constraints, and boundary conditions. Since the seminal paper by [1], topology optimization has undergone a remarkable development in both fields of academic research [2-4] and industrial application [5]. Over the past few decades, most contributions to topology optimization have focused on minimizing compliance (maximizing stiffness) or other global response functions constrained by volume fractions [6-8]. In engineering practice, strength is a basic design criterion. Since the groundbreaking paper [9], on this topic, continuous efforts such as stress-constrained topology optimization with different kind of criteria from von-Mises [10,11], with geometrical non-linearity [12] and considering thermal stress [13] have been made. However, stress-related topological optimization is far from mature compared with compliance-based topological optimization. Stress-related topology optimization imposes four main challenges [14-16]. The first issue is the singularity phenomenon of stress constraints in topological design, which results from the irregularity or degeneracy of the design space. To resolve this problem, many relaxation techniques have been proposed, such as the  -relaxed method [9,17] and qp approach [18,19]. The second challenge is that the local nature of stress, whether direct or adjoint method [9], will generate a large number of local stress constraints included in the optimization formulation, which leads to a prohibitive numerical computation burden. To address this issue, the active set strategy [9] has been proposed. Additionally, an alternative way to handle the large number of local stress constraints is the adoption of an aggregation function, such as the p-norm global stress and the Kresselmeier–Steinhauser (KS) functions [20]. Global stress measurement has computational efficiency in sensitivity assessment, but at the cost of losing sufficient control over local stress behavior [21]. The third noteworthy challenge is highly nonlinear stress behavior, where stress distributions are sensitive to even subtle topological changes, especially in critical areas of stress concentration [22]. To stabilize the convergence, a density filtering method was adopted by Le et al. [14] to achieve stable iteration and acceptable design solutions for the stress-constrained topology optimization problem. Another challenge is the 4

accuracy of stress assessments. For instance, Svärd [15] pointed out that inaccuracy in the evaluation of the stress field is related to the jagged nature of the optimal structure coming from the finite element mesh; Bruggi [16] compared optimal solutions achieved using stresses as main variables in mixed finite elements with conventional solutions based on displacement-based finite elements. The last three challenges are common issues in stress-related topology optimization problems. Compared to these three issues, the difference is that the singularity problems are only occurred in density-based approaches with intermediate densities. Unlike the previous method (  -relaxed method and qp approach), the method based on discrete variables, such as the evolutionary structural optimization (ESO) method [23] and its improved bi-directional ESO (BESO) approach [24] as significant branches of topology optimization have been widely adopted for both academic researches and engineering applications [25-27] and can avoid the singularity problem naturally. Nonlinear structural topology optimization has received increasing attention owing to the fact that nonlinear structures are among the most significant considerations in engineering applications. Early works on topology optimization were restricted to linear structural designs. In pursuing more realistic designs, continuous efforts have been conducted to extend topology optimization to nonlinear structural designs considering various sources of non-linearity, such as geometrical non-linearity [12,28-34], material non-linearity [35-40], and both geometrical and material non-linearities simultaneously [41-43]. The difficulty in topology optimization of nonlinear structures arises primarily from a large portion of computational time of the topology optimization. The computation time in nonlinear finite element analysis depends strongly on the size of the structures. Thus, the computation efficiency is critical especially for 3D problems. Moreover, in mathematical optimization methods such as SIMP the material density is penalized so that intermediate material densities invaluably appear in the topologies. As the topology develops, large displacements may cause the tangent stiffness matrix in low-density elements to become indefinite or even negative definite [44,45]. Then some additional scheme must be devised to circumvent the problem such as removing low-density elements [45] or relaxing the convergence criterion [44]. It seems that BESO method with discrete formulation has the potential to overcome these problems. 5

To our knowledge, most of the current literature on stress-constrained topology optimization is just for linear continuum structure. Notably, this is the first investigation of stress-constrained topology optimization of geometrically nonlinear continuum structures. By using an extended BESO method, the current concepts of stress-constrained topology optimization are extended to the design of geometrically nonlinear structures in this paper. The popular p-norm stress aggregation function is adopted to treat the vast number of local constraints imposed on all elements. An efficient sensitivity number formulation for stress-constrained geometrically nonlinear structures is derived from the adjoint method. As for the highly nonlinear stress behavior, the updated scheme takes advantages from two filters respectively of the sensitivity and topology variables to improve convergence. Moreover, the filtered sensitivity numbers are combined with their historical sensitivity information to further stabilize the optimization process. Several numerical examples demonstrate the validity of the presented method, in which the stress-constrained geometrically nonlinear designs are compared with traditional stiffness-based designs and stress-constrained linear designs to illustrate the merit of considering geometrically nonlinearity. It is observed that the proposed approach produces solutions that reduce stress concentration at the critical stress areas of geometrically nonlinear structures. The influence of varying stress constraints, prescribed target volume fractions, p-norm parameters and magnitude of external loads on the optimized results are investigated, and also demonstrate that the consideration of stress constraints in geometrically nonlinear structures is indispensable. The layout of this paper is organized as follows. Section 2 builds the stress-constrained topological design model of geometrically nonlinear structures including the geometrically nonlinear finite element analysis, effective stress models, p-norm stress measure and the definition of the topology optimization problem. Section 3 detailly derives the sensitivity numbers for the geometrically nonlinear structure under the global p-norm stress constraints. The extended BESO numerical techniques are proposed in Section 4. Section 5 validates the proposed method through a series of benchmark design problems. Finally, conclusions and future perspectives are drawn in Section 6.

6

2. Nonlinear analysis and optimization problem statement 2.1. Geometrically nonlinear analysis For many industrial applications, the maximum stiffness of a structure is pursued. Consider a typical nonlinear structure subject to an applied load which increases monotonously with displacement up to a maximum value, F * . The corresponding nonlinear load-displacement curve is depicted in Fig. 1. To maximize the structural stiffness, the natural choice of the optimization objective is to minimize the displacement U* or the end compliance F *T U * in the deflected configuration [30]. However, the minimization of the end compliance may result in degenerated structures which can only support the maximum load they are designed for. To avoid this problem and make sure that the structure is stable for any load up to the maximum design load, one may minimize the complementary work WC shown as the area bounded by OAB in Fig. 1. 2.1.1. Displacement-strain conversion matrix Element design variables and quadrangle elements are adopted to discretize the design domain. Independent displacement field U can be interpolated by node displacement as

U  N de

(1)

where N and d e denote the shape function and elemental node displacement vector, respectively. From the definition of Green strain [46], it can be written as linear and nonlinear parts, as follows

E  EL  E N

(2)

Eq. (2) can be expressed in matrix form as follows  E L  B L 0d e , B L 0  LN   1 E N  2 B N 0d e , B N 0  AG

7

(3)

where L is differential operator. B L 0 and B N 0 denote the linear and nonlinear transformation matrixes between node displacement and element strain, respectively.

A and G denote the derivative matrix of nodal displacement with respect to coordinates and the derivative matrix of the shape function with respect to coordinates , respectively. Substitute Eq. (3) to Eq. (2), and the Green strain can be further written as

E  Β0d e

(4)

B0  BL0  1 2 B N 0  BL0  1 2 AG

(5)

where B0 denotes the transformation matrix between node displacement d e and element strain E . Because A contains the derivative of the displacement, the relationship between Green strain vector E and node displacement vector d e is nonlinear. From Eqs. (4-5), the relationship between the increment of Green strain  E and the increment of node displacement  de is as follows

 E   E L   E N  B L 0 de 

1 1  A  Gde  AG de  B L 0 de  AG de 2 2

(6)

Eq. (6) can be rewritten as

 E  B0 de

(7)

B0  B L 0  B N 0

(8)

where B0 is defined as the transformation matrix between  E and  de . 2.1.2. Residual of the equilibrium equation The second type of Piola-Kirchhoff stress vector is denoted as S . The physical force and surface load are denoted by p 0 and q 0 respectively, and the virtual work equation within the element can be expressed as



e0

 ET SdV    uT p0 dV    uT q0 dA e0

Ae 0

8

(9)

where the integral domains e0 and Ae 0 are the volume of the element in its initial configuration and the corresponding surface area, respectively. It can be obtained from Eq. (1) that

 U  N de

(10)

Substitute Eq. (7) and Eq. (10) to Eq. (9), and consider the arbitrariness of  de , Eq. (9) can be further written as



e0

B0T SdV   NT p0 dV   NT q0 dA e0

Ae 0

(11)

The element node displacement vector d e is connected with the total displacement vector d of the finite element by the selection matrix c , i.e.

 de  c d

(12)

Then the equilibrium equation of the element can be rewritten as cT  B0T SdV  cT  NT p0 dV  cT  NT q0 dA e0

e0

Ae 0

(13)

By integrating the equilibrium equations of all the elements, the equilibrium equation of the whole finite element system can be obtained. If the integral over the whole area of the object V0 and over the entire boundary A0t means summation, then the equilibrium equation of the system can be simplified as



V0

B0T SdV   NT p0 dV   NT q0 dA V0

A0 t

(14)

Eq. (14) can be rewritten as t R  d    B0T SdV  F  Fi n  F0 V0

(15)

where R can be defined as the residual of the equilibrium equation. F int represents the internal force vector, and F denotes the external force vector, i.e. F   NT p0 dV   NT q0 dA V0

A0 t

(16)

2.1.3. Tangential stiffness matrix In many applications, the equilibrium equation needs to be solved iteratively by Newton-Raphson method. To this end, it is necessary to establish the equilibrium 9

equation to represent the structure in an incremental form. The corresponding stiffness matrix is the tangential stiffness matrix. According to the definition of tangential stiffness matrix and assume that external load is independent of node displacement, then we have

 R    B0T  S  dV     B0T  SdV  K T  d V0

V0

(17)

where K T represents the tangential stiffness matrix. The first integral of Eq. (17) can be written as

 B V0

T 0

 S  dV   B0T DB0 dV  d  K M  d V0

(18)

where K M is the tangential stiffness matrix related to the constitutive tensor, which can be expressed as

KM  KL  K N

(19)

K L   B L 0T DB L 0 dV

(20)

where V0

K N   B L 0T DB L 0 dV   B L 0T DB L 0 dV   B L 0T DB L 0 dV V0

V0

V0

(21)

where K L is the linear stiffness matrix under the assumption of small displacement, and K N is caused by large displacement. The second integral of Eq. (17) can be written as

   B  SdV   T 0

V0

V0

GT  AT SdV

(22)

Based on matrix transformation, Eq. (22) can be expressed as

  G MGdV  d  K  d T

S

V0

(23)

where K S is the tangential stiffness matrix caused by the stress state. And matrix

M is composed of the components of the stress vector S . Then the tangential stiffness matrix of the system can be written as

KT  K L  K N  K S

10

(24)

2.1.4. Solution of equation The Newton-Raphson method is used to solve the nonlinear finite element model in this paper. Newton-Raphson method can be drawn as the following flow chart. where superscript n represents the physical quantity corresponding to the nth iteration step, and K Tn  K T  d n  . The global tangential stiffness matrix K T can be obtained by solving the nonlinear finite element equations iteratively as shown in Fig. 2.

K T  K Tlast

(25)

where K Tlast denotes the global tangential stiffness matrix obtained in the last iterative step. In practical analysis, d 0  0 can be usually taken. The magnitude of

 d or  R may vary according to the convergence accuracy requirements. d 0  0 and  d  0.001 are adopted in this paper. After the last iterative step, the stress and strain of each element can be obtained according to the node displacement. 2.2. Effective stress model The design domain is divided into nel quadrilateral elements, and the nel-dimensional design variable vector is represented by x  ( x1 , x2 ,

, xnel )T . By

using BESO method [8], the elemental pseudo densities can be defined as

xi  0 or 1 ,

i  1, 2,

, nel

(26)

where the one and zero represent solid and void elements, respectively. The effective constitutive matrix is linked to the design variable xi and can be expressed as

Di  xi D0

(27)

where D0 denote the constitutive matrix of solid element. In practical application, in order to avoid re-meshing and the singularity of the overall stiffness, the cavity

11

element adopts very low stiffness. The stress at the center of element is adopted as the element effective stress by

σi  (Di Eilast ) / xi  D0Eilast

(28)

where Elast which is obtained in the last iterative step of nonlinear element analysis, i denotes the strain vector of i th element. last Elast  Blast i i ui

where

Blast i

and

ulast i

respectively denote the transformation

(29) matrix

of

displacement-strain and the node displacement vector of i th element which are also obtained in the last iterative step of nonlinear element analysis. On the basis of the definition of Eq. (28), it can be found that the element stress is independent of its corresponded design variable. 2.3. Global p-norm stress measure On basis of the last nonlinear finite element iteration solution, the stress of each element can be calculated according to the node displacement. The p-norm stress measure can be calculated as 1/ p

 PN

 nel p      vm ,i   i 1 

(30)

where  vm,i and p denote the von Mises stress at the center of the i th element and the stress norm parameter, respectively. The p-norm approaches the maximum stress

 max when p   and gives on the other hand the average stress when p  1 . From Eq. (30), seemingly, a larger value of p is desirable as it gives a more precise approximation of the maximum stress from any point of the design domain. However, the global p-norm stress measure problem becomes ill-conditioned when

p increases. In the process of topology optimization, the larger the value p , as the algorithm only for peak stress, and not for other peak stress, it will produce serious oscillation. Therefore, an appropriate norm parameter value must be selected so that the smoothness of the problem does not deteriorate and the maximum stress can be 12

approximated sufficiently. A detailed discussion on the property of p-norm global measure can refer to [47]. The element von Mises stress, in general, is calculated on the centroid as

 vm,i   σTi Vσi 

1/2

(31)

where σ i and V denote the stress vector of the i th element and the constant stress coefficient matrix, respectively. 2.4. Topology optimization model With regard to the geometrically nonlinear topology optimization problem, the mathematical model for maximizing stiffness with both constraints on volume fraction and the maximum von Mises stress is described as follows Find: x   x1 , x2 ,

, xnel 

1 n  Minimize: f  x   W c  lim   FT  U h  U h 1   n   2 h 1 

Subject to:

R  F  Fi n t  0

(32a) (32b)

nel

V (x)   xi vi  Vreq

(32c)

i 1

max *  vm   vm

xi  0 or 1 ,

(32d) i  1, 2,

, nel

(32e)

max * where  vm and  vm are the maximum stress and the value of stress constraints,

respectively. vi , V  x  and Vreq are the volume of the i th element, the total material volume and the required material volume, respectively. U is the displacement vector, h is the incremental number of the load vector and n is the total number of load increments. The size of the increment is determined by F = F n

(33)

From Eqs. (32a) and (33), a higher value of n is desirable as it gives a more precise approximation of the complementary work WC. However, as n increases, the computational time consumed in nonlinear finite element analysis will increase so the total calculation time is extended. 3. Sensitivity analysis 13

During optimization, the sensitivity of design objectives to design variables needs to be analyzed to update the structural topology. In BESO framework, the traditional method to consider additional constraints in addition to volume constraints is to extend the design objective function by introducing Lagrange multiplier [48,49]. Following the same strategy, we represent the modified design objective as max * f1  x  = f  x      vm   vm 

(34)

with   0 . When the stress constraint is not violated, the augmented design objective is equivalent to the original one with   0 . Otherwise, an appropriate value  is found through the iterative process until the stress constraint is satisfied, thus making the design objective a compromise between compliance and maximum von Mises stress. The minimization of f1  x  is equivalent to the minimization of max f1*  x  = f  x      vm

(35)

The equivalence of the optimization procedure using two objectives f1  x  and

f1*  x  is guaranteed on the condition that the  is iteratively determined upon max * and  vm .  vm

In order to avoid the local nature of stress, the design objective is further modified max and replaced by the maximum stress  vm with

f 2  x  = f  x      PN

(36)

The derivative of the modified design objective f 2  x  with respect to xi equals f 2  x  f  x   =   PN xi xi xi

(37)

The first term in Eq. (37) can be easily derived using the adjoint method. The sensitivity of the complementary work with respect to the design variable xi is 1 n f  x   U U h 1    lim    FhT  FhT1   h   xi xi   n   2 h 1  xi

(38)

An adjoint equation is introduced by adding a series of vectors of Lagrangian multipliers  h into the objective function as

14

f  x   lim n 

1 n   FhT  FhT1   Uh  Uh1   hT  R h  R h1  2 h1 

(39)

where Rh and R h1 are the residual forces in load increment steps h and h  1 . Thus

R h  R h1  Fh  Fhint  Fh1  Fhint1  0

(40)

Because both Rh and R h1 are equal to zero, the modified objective function in Eq. (39) is the same as the original objective function in Eq. (32a). The sensitivity of the modified objective function is f  x   U U h 1  1 n   lim   FhT  FhT1   h   xi xi  n  2 h 1   xi   R U h R h 1 U h 1   R h  R h 1     hT  h    xi  U h xi U h 1 xi  

(41)

Note that the derivative of  h is not included in Eq. (41) because it would be multiplied by

 R h  R h1 

which is zero. It is assumed that there is a linear

force-displacement relationship in a small load increment. Therefore,

R h R h 1   K th U h U h 1

(42)

where K th is the tangential stiffness matrix in the h th step. By substituting the relationship Eq. (42) into Eq. (41), the sensitivity of the modified objective function can be rewritten as f  x   U U h 1  T   R h  R h 1   1 n   lim   FhT  FhT1   hT K th   h     h xi xi  xi n  2 h 1   xi 

(43)

In order to eliminate the unknowns Uh xi  Uh1 xi ,  h can be chosen as

K thh  Fh  Fh1

(44)

This equation defines the adjoint system. From the assumption of linear force-displacement relationship in a small increment, the increment of the force can be approximately expressed by

K th  U h  U h 1   Fh  Fh 1

(45)

Comparing Eq. (44) with Eq. (45),  h can be obtained as

h  Uh  Uh1

15

(46)

By substituting  h into Eq. (43) and utilizing Eq. (40), the sensitivity of the objective function can be expressed as

f  x   Fint Fint  1 n   lim   UTh  UTh1   h  h1  xi xi  n  2 h 1  xi

(47)

In the evolutionary structural optimization method, a structure can be optimized by removing and adding elements. That is to say that, the element itself is treated as the design variable. Thus, when one element is totally removed from the system, the variation of the objective function is 1 n  f  x    lim    UTh  UTh 1  Fhint  Fhint1   n   2 h 1 

(48)

Form Eq. (15), the variation of internal force due to removing one element can be written as

Fhint  cT Fhe

(49)

Substituting Eq. (49) into Eq. (48), the variation of the objective function can be expressed as

f  x  1 n   f  x    lim    UTh  UTh 1  CeT Fhe  CeT Fhe1   xi n   2 h 1    lim  E  E n

n 

h 1

e h

e h 1

  E

(50)

e n

where Ene is the total strain energy of the removed element. Eq. (50) means that the decrease of the total external work due to removing one element is equal to the total strain energy of the element in its final deformed state and irrelevant to the size of displacement intervals. Next, the second term of Eq. (37) is derived. According to Eq. (30), the derivative of p-norm stress  PN with respect to the i th design variable xi can be expressed as

 nel p 1  vm. j   PN   1PN p   vm  ,j xi xi   j 1

(51)

By differentiating Eq. (31) with respect to the stress vector, we can obtain

 vm, j σ j



1 2 1 T 1 T σ j Vσ j  2σTj V   vm  , jσ j V 2

(52)

With the help of Eq. (52), the derivative of the von Mises stress of the j th

16

element with respect to the i th design variable can be obtained by the chain rule

 vm, j xi



 vm, j σ j σ j

xi

1 T   vm , jσ j V

σ j xi

(53)

In consideration of Eq. (28), it is known that the constitutive matrix D0 of solid element and the last elemental strain-displacement transformation matrix B last i obtained from nonlinear finite element analysis are independent with the corresponding design variables, the derivation in Eq. (53) can be further written as

 vm, j xi  vm, j xi



1 T vm, j j



1 T vm , j j

σ VD0B

last j

ulast j

(54)

xi

Ulast σ VD0B L j xi last j

(55)

where the matrix L j gathers the nodal displacement vector of the j th element from the global displacement vector satisfying

ulast  L j Ulast j

(56)

By differentiating the last state equilibrium equation

K tlast Ulast  F

(57)

K tlast last U last U  K tlast 0 xi xi

(58)

One obtains

Therefore, Eq. (55) can be approximately expressed as  vm, j xi

1 T last t   vm , j σ j VD0 B j L j  K last 

1

K tlast last U xi

(59)

Substituting Eq. (37) into Eq. (29), therefore Eq. (29) can be further rewritten as t nel  t 1 K last  PN 1 p  p 2 T last   PN   vm, j σ j VD0B j L j   K last  Ulast xi xi  j 1 

(60)

With the solution of the following adjoint problem p 2 last K tlast μ    vm , j  D0 B j L j  Vσ j nel

T

(61)

j 1

The sensitivity numbers can be further simplified as t   PN 1 p T  K last   PN μ  Ulast  xi  xi 

17

(62)

Recalling the effective constitutive matrix Eq. (27), the derivative of the tangent stiffness matrix with respect to the i th design variables equals K tlast nel T k j 0  Lj L j  LTi k i  Li xi xi j 1

(63)

0 where k i and k i  are the stiffness matrix of the i th element and the elemental

stiffness corresponding to solid element, which can be respectively written as k i  xi 

B 

D0 Bilast d i

(64)

k i   

B 

D0 B last i d i

(65)

i

0

i

last T i

last T i

where i represents the domain of the i th element. After substituting Eq. (63) into Eq. (62), the sensitivity can be eventually evaluated as

 PN   1PN p μTi k i 0ulast i xi

(66)

where μ i denotes the vector of the adjoint nodal values of the i th element. Finally, substituting Eq. (40) and Eq. (66) into Eq. (37), the sensitivity of modified objective function is f 2  x  0 =  Ene     1PN p μTi k i ulast i xi

(67)

The iterative determination scheme for the Lagrange multiplier  is provided in detail in the following section. 4. Numerical techniques Within BESO framework, the sensitivity numbers denoting the relative ranking of elemental contribution of the design objective are used to guide material removal and addition [22]. Based on Eq. (67), the sensitivity number for the modified objective function is defined as

 i   xi

f 2  x  0 = xi Ene     1PN p μTi k i ulast i xi





(68)

where xi denotes an indicator to ensure that sensitivity numbers are evaluated only

18

for solid elements ( xi  1 ) and are enforced to zero for void elements ( xi  0 ). 4.1. Update scheme of design variables The update scheme of design variables in the framework of BESO method has been proposed by [24] and detailly summarized by [8,21]. This is briefly stated in this subsection. l The target volume fraction V   at the current optimization iteration ( l th) is

calculated by



V    max Vreq , 1  cer V  l

l 1



(69)

where cer denotes the evolutionary ratio. Sensitivity numbers are firstly smoothed by [50]

  

nel

i

w 

sen j j 1 ij nel sen j 1 ij

(70)

w

where wijsen denotes weight factor

wijsen  max 0, rsen    i, j 

(71)

where rsen and   i, j  are the sensitivity filter radius and the element center-to-center distance, respectively. Same as Ref. [21], the sensitivity numbers are averaged with their historical information of last two design iterations l 

i 

 il    il 1   il 2 3

,

for l  2

(72)

All the elements are divided into a vector based on their sensitivity. Then two th th threshold parameters  del and  add are used to update the topology design

variables for material removal and addition, respectively [8]

19

xi l 1

 0    th , x  l   1 i del i  th   1  add   i , xi l   0  l  otherwise  xi

(73)

th th The parameters  del and  add are obtained from a iterative algorithm [21,24].

An additional filtering on the topological variables using the same scheme as Eq. (60) is proposed

  

nel

xi

den j j 1 ij nel den j 1 ij

w x

(74)

w

where the linear weight factor wijden is defined as

wijden  max 0, rden    i, j 

(75)

with density filter radius rden . 4.2. Determination of the Lagrange multiplier The validity of the above variable updating scheme essentially depends on the Lagrange multiplier. An empirical scheme [22] is developed for the determination of the Lagrange multipliers so that the maximum von Mises stress could be effectively constrained through the controlling of the aggregated p-norm stress. The Lagrange multiplier is initially set to zero, reducing the augmented design objective function to the traditional compliance. In order to determine the Lagrange multiplier, experimental updates and additional finite element analysis of design variables are required. The Lagrange multiplier is increased iteratively in an trial exponential manner if the maximum von Mises stress  vm of the trial test violates * the critical von Mises stress  vm . trial *   k 1    k   2k , if  vm   vm

(76)

k where k and    denote iteration number and the Lagrange multiplier at the k th

update iteration, respectively. trial Once the trial stress  vm is below the critical stress, stop the operation of Eq.

20

(76). The target Lagrange multiplier can be obtained by





k 

trial * ,  target  min  vm   vm

k  1, 2,

, kmax

(77)

where kmax denotes the prescribed maximum iteration number. In such case, the trial * algorithm finds  vm is closest to  vm among all iterations, and the corresponding

Lagrange multiplier value is adopted. 4.3. Overall update algorithm The sensitivity and density filter radius are defined by satisfying rden  rsen  he , where he is the typical element length. The required number of optimization iterations to achieve the target prescribed volume fraction from a full solid design can be estimated according to the evolutionary ratio cer . As shown in Fig. 3, the flow chart of the proposed algorithm procedures is concluded as: 1. Define design domain, load case, boundary conditions and finite element mesh. 2. Initialize the topological design full of solid elements. 3. Prescribed BESO parameters: Vreq , cer , carmax , rsen , rden . 4. Estimate the number of optimization iterations upon Vreq and cer . 5. Loop: (a) Calculate the target volume fraction in the current iteration as Eq. (69); (b) Perform nonlinear finite element analysis iteratively; (c) Determine 

iteratively according to Eqs. (76) and Eq. (77);

(d) Evaluate sensitivity numbers according to Eq. (68); (e) Filter and average sensitivity numbers by Eq. (70) and Eq. (72); (f) Update design variables according to the scheme in Eq. (73); (g) Filter topology variables upon Eq. (74); (h) Convert the filtered topology into discrete design; 6. End.

21

5. Numerical examples The effectiveness of the method proposed in this paper is certified through a series of 2D benchmark examples in this section. Quadrilateral elements are used in three 2D cases. Young’s modulus and Poisson’s ratio of the solid material are set to be 30 106 and 0.3, respectively. All the parameters involved in the extended BESO

method shown in section 4 remain unchanged in the following examples. The evolutionary rate cer determining the relative percentage of material to be removed at each design iteration is set to be cer  1% . The maximal admission ratio corresponding to the allowed maximal percentage of recovered material during per iteration is set to be carmax  0.5% . The stress norm parameter are set to be p  6 . The evolutionary design process starts with an initial design filled with solid materials. 5.1. L-bracket design Figure 4 shows the design domain of the well-known L-bracket benchmark with characteristic dimensionless dimensions. Elements of length he  1 are set to discretize the L-bracket (200 elements along the left vertical boundary). The top end of the L-bracket is clamped, and the vertical load F =16 is distributed over eight nodes (red dots as shown in Fig. 4) to avoid stress concentration. The volume fraction is set to be 50%. The filter radiuses for sensitivity numbers and topology variables are set to be thrice the typical finite element size rsen  3he and rden  6he , respectively. Figure 5 shows the comparison of geometrically nonlinear topological designs of the L-bracket benchmark under different stress constraints. Fig. 5(b-d) shows the * corresponding design results with  vm =8.5, 7.7, 7.2, respectively. The conventional

geometrically nonlinear stiffness design without stress constraint is shown in Fig. 5(a). It can be observed from Fig. 5 that a stricter stress constraint yields a higher compliance, which is physically reasonable from the constitution of the design objective. The stress distribution of the optimal topologies under stress constraint becomes smoother. This shows that the proposed method can effectively solve the stress concentration problem of geometrically nonlinear structures. By comparing

22

with geometrically nonlinear stiffness design without stress constraint as shown in Fig. 5(a), it can be found that the maximum von Mises stress is controlled, however with * * an apparent compliance increase up to 1.627% (  vm =8.5), 6.770% (  vm =7.7), * 12.501% (  vm =7.2) (stiffness decrease).

The detailed topological evolution with von mises stress distribution for the case of Fig. 5(d) is shown in Fig. 6. As can be observed, the stress concentration appears initially at the L-bracket corner when the design domain is full of solid material. The maximal von Mises stress reduces from the initial 7.4703 to 6.7375 during the first design iteration through the material removal at the L-bracket corner. The maximal von Mises stress has been effectively reduced for 12.795% from 7.4703 to 6.5145 through the material removal at the L-bracket corner and smoothed boundaries of the topological structure of L-bracket. However, the compliance has been increased for 36.435% from 5.8474×10-4 to 7.9779×10-4. Figure 7 shows the comparison of topological designs and von Mises stress distribution of the L-bracket under different cases. By comparing Fig. 7(a) and (b), it can be seen that the maximal stress and stiffness for geometrically nonlinear stiffness design without stress constraint are respectively smaller and larger than that for linear stiffness design without stress constraint. This indicates that geometrically nonlinear structures can absorb more strain energy. Comparing Fig. 7(b) and (c), the maximal stress and compliance for geometrically nonlinear stress design are respectively smaller and larger than that for geometrically nonlinear stiffness design without stress constraint. Stress design pursues maximum stress minimization at the expense of structural stiffness. Comparing Fig. 7(b) and (d), the maximal von Mises stress and compliance for geometrically nonlinear stiffness design with stress constraint are respectively smaller and larger than that for geometrically nonlinear stiffness design without stress constraint. Stress-constrained design seek to control the maximum stress of the structure at the expense of structural stiffness. Comparing with linear stiffness design without stress constraints, thin parts of the topological structures for stiffness design with geometrically nonlinearity and stress design are disappeared. As can be observed, the stress concentrations are eliminated for stress design and stress-constrained geometrically nonlinear stiffness design. The corner of the L-bracket of the topological structure for stress design and geometrically

23

nonlinear stiffness design with stress constraint are smoother. This shows that the proposed method can effectively solve the stress concentration problem of geometrically nonlinear structures under stress constraints. Figure 8 shows the sensitivity distribution (3th iteration) of the L-bracket benchmark under different cases. The sensitivity distributions around the corner of the L-bracket for stress design and stress-constrained stiffness design are smoother than that for stiffness design without stress constraints. Figure 9 shows the comparison of the influence of volume fraction on the topological structure and von Mises stress distribution of the L-bracket benchmark under the case of geometrically nonlinear stiffness design with stress constraint * =7.2. With the increase of volume fraction, the objective function decrease  vm

gradually. The maximum von Mises stress of the topological structures are less than the prescribed stress constraint. 5.2. MBB beam with one pre-existing crack notch The MBB beam with one pre-existing crack notch is illustrated in Fig. 10. The beam is discretized with 360 by 120 square elements of length he  1. Vertical load

F =34 is distributed over 17 nodes to avoid stress concentration. The crack notch of length 30 at the middle of the bottom edge. The volume fraction is set 50%. The filter radiuses for sensitivity numbers and topology variables are set to be thrice the typical finite element size rsen  3he and rden  4he , respectively. Figure 11 shows the comparison of geometrically nonlinear topological designs of the MBB beam under different stress constraints. Fig. 11(b-d) shows the * corresponding design results with  vm =6.5, 6.0, 5.8, respectively. The conventional

geometrically nonlinear stiffness design without stress constraint is shown in Fig. 11(a). Similar to the first example, a stricter stress constraint yields a higher compliance. By comparing with geometrically nonlinear stiffness design without stress constraint, it can be found that the maximum von Mises stress is controlled, * however with an apparent compliance increase up to 1.183% (  vm =6.5), 7.070%

24

* * (  vm =6.0), 10.958% (  vm =5.8) (stiffness decrease).

The detailed topological evolution with von mises stress distribution for the case corresponding to Fig. 11(d) is shown in Fig. 12. Fig. 13 shows the comparison of topological designs and von Mises stress distribution of the MBB beam under different cases. By comparing Fig. 13(a) and (b), the maximal stress and compliance for geometrically nonlinear stiffness design without stress constraint are respectively smaller and larger than that for linear stiffness design without stress constraint. Comparing Fig. 13(b) with (c), the maximal stress and compliance for geometrically nonlinear stress design are respectively smaller and larger than that for geometrically nonlinear stiffness design without stress constraint. Stress design pursues maximum stress minimization at the cost of the decrease of structural stiffness. Similarly comparing Fig. 13(b) with (d), the maximal stress and compliance for stress-constrained geometrically nonlinear stiffness design are respectively smaller and larger than that for geometrically nonlinear stiffness design without stress constraint. Also, in some degree, stress-constrained design can control the maximum stress of the structure at the expense of structural stiffness. Comparing with linear stiffness design without stress constraints, similar to the first example, thin parts of the topological structures for stiffness design with geometrically nonlinearity and stress design are disappeared. The stress concentration is eliminated for stress design and stress-constrained geometrically nonlinear stiffness design. The corner of the crack notch of the topological structure for stress design and stress-constrained geometrically nonlinear stiffness design are smoother. This shows that the proposed method can effectively solve the stress concentration problem with stress requirement the objective or the constraint function. Figure 14 shows the evolution histories of the MBB beam with one pre-existing crack notch under different cases. Comparing with stiffness design without stress constraints, the curves of the topological structures for stress-constrained geometrically nonlinear stiffness design and geometrically nonlinear stress design are more oscillation. This is due to the geometrically nonlinearity and the highly stress

25

nonlinear behaviour. As can be observed, BESO starts from the full design and gradually decreases the total volume of the structure using an evolutionary rate 1%. The curves converge to almost constant values when the volume fraction remains constant at the final stage. Figure 15 shows the comparison of the influence of p-norm on the topological structure and von Mises stress distribution of the MBB beam under the case of * geometrically nonlinear stiffness design with stress constraint  vm =5.8. With the

increase of p-norm, the objective function also increase gradually. The maximum von Mises stress of the topological structures are less than the prescribed stress constraint. The comparison of the influence of external force on the topological structure and * von Mises stress distribution of the MBB beam with stress constraint  vm =5.8 are

shown in Fig. 16. With the increase of external force, the objective function also increase gradually. The maximum von Mises stress of the topological structures satisfies the prescribed stress constraint. 5.3. U-shaped structure design The design domain of the U-shaped structure is shown in Fig. 17 with characteristic dimensionless dimensions. The U-shaped structure is divided into 200 elements along the bottom horizontal boundary and 160 elements along the left vertical fixed boundary. The vertical load F =12 is distributed over six nodes (red dots as shown in Fig. 17) to avoid stress concentration. The volume fraction is set to be 50%. The filter radiuses for sensitivity numbers and topology variables are set to be twice the typical finite element size rsen  2he and rden  4he , respectively. The comparison of geometrically nonlinear topological designs of the U-shaped structure under different stress constraints are shown in Fig. 18. Fig. 18(b-d) shows * the corresponding design results with  vm =7.0, 6.6, 6.2, respectively. The

conventional geometrically nonlinear stiffness design without stress constraint is shown in Fig. 18(a). Similar to the first two numerical examples, a stricter stress constraint results in a higher compliance. By comparing with geometrically nonlinear stiffness design without stress constraint as shown in Fig. 18(a), it can be seen that the maximum von Mises stress is controlled, however with an apparent compliance

26

* * * increase up to 9.649% (  vm =7.0), 12.652% (  vm =6.6), 16.640% (  vm =6.2) (stiffness

decrease). Figure 19 shows the detailed topological evolution with von mises stress distribution for the case corresponding to Fig. 18(d). It can be observed that the stress concentration appears initially at two corners of the U-shaped structure when the design domain is full of solid material. The maximal von Mises stress has been effectively reduced for 8.244% from 5.5338 to 5.0776 through the material removal at two corners of the U-shaped structure and smoothed boundaries of the topological structure of the U-shaped structure. Figure 20 shows the comparison of topological designs and von Mises stress distribution of the U-shaped structure under different cases. The maximal stress and stiffness for geometrically nonlinear stiffness design without stress constraint are respectively smaller and larger than that for linear stiffness design without stress constraint. Geometrically nonlinear structures can absorb more strain energy. Meanwhile, the maximal stress and compliance for geometrically nonlinear stress design are respectively smaller and larger than that for geometrically nonlinear stiffness design without stress constraint. Stress design pursues maximum stress minimization at the expense of structural stiffness. The maximal stress and compliance for stress-constrained geometrically nonlinear stiffness design are respectively smaller and larger than that for geometrically nonlinear stiffness design without stress constraint. To a certain extent, stress-constrained design can control the maximum stress of the structure at the expense of structural stiffness. Figure 21 shows the evolution histories of the U-shaped structure under different cases. Similar to the second example, comparing with stiffness design without stress constraints, the curves of the topological structures for stress-constrained geometrically nonlinear stiffness design and geometrically nonlinear stress design are more oscillation because of the geometrically nonlinearity and the highly stress nonlinear behaviour. 6. Conclusions The main contribution of this work is to extend the BESO method for stiffness maximization design of geometrically nonlinear continuum structures with both volume fraction and maximum von Mises stress constraints. BESO method based on 27

discrete variables can effectively avoid the well-known singularity problem in density-based methods with low density elements. The maximal stress is approximately evaluated using the conventional p-norm global stress. In order to consider the additional stress constraints, the traditional compliance design objective function is extended by introducing a Lagrange multiplier and using the p-norm stress measurement method. For highly nonlinear stress behavior, the optimization process is stabilized by filtering topological variables and sensitivity numbers. The optimal results of several 2D benchmark design problems clearly demonstrated the effectiveness and practicability of the proposed method. The main conclusions from this research are as follows: (1) The adopted p-norm global stress measure can well restrain the stress level and solve the stress concentration problem of geometrically nonlinear structures. (2) The compliance of topological structure increases and the stress distribution becomes smoother with the strict stress constraint. (3) The compliance and the maximal von Mises stress of geometrically nonlinear stiffness design are smaller than those for linear design. (4) The stress-constrained design is seek to control the maximum stress of the structure at the expense of structural stiffness. (5) The evolution histories of stress-constrained geometrically nonlinear stiffness design are more oscillation due to the geometrically nonlinearity and the highly stress nonlinear behavior. The work presented in this paper can be seen as a useful step towards a more realistic and comprehensive framework for integrating stress-constrained stiffness design (pursue the balance between structural strength and stiffness) and the topological design of geometrically nonlinear structures. However, additional researches are needed to include other stress-constrained nonlinear problems, such as stress-constrained simultaneous geometrical and material nonlinearity, fatigue, damage and fracture etc. Finally, the proposed method is also expected to be extended and applied for the design of 3D benchmark design problems and multi-material micro-structures to restrain the stress level at the lower scale within the context of multiscale modeling and nonlinear design.

28

Acknowledgements This work was sponsored by the National Natural Science Foundation of China (11872311). Disclosure statement No potential conflict of interest was reported by the authors. References [1] M.P. Bendsøe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Eng. 71(2) (1988) 197–224. doi: 10.1016/0045-7825(88)90086-2. [2] M.P. Bendsøe, O. Sigmund, Topology optimization: theory, methods and applications, (2003) Springer, Berlin. [3] J.D. Deaton, R.V. Grandhi, A survey of structural and multidisciplinary continuum topology optimization: post 2000, Struct. Multidiscip. Optim. 49(1) (2014) 1–38. doi: 10.1007/s00158-013-0956-z. [4] X. Huang, Y.M. Xie, Evolutionary topology optimization of continuum structures:

methods

and

applications,

(2010)

Wiley,

Chichester.

doi:

10.1002/9780470689486. [5] J. Zhu, W. Zhang, L. Xia, Topology optimization in aircraft and aerospace structures design, Arch. Comput. Methods Eng. 23(4) (2016) 595-622. doi:10.1007/s11831-015-9151-2. [6] Z. Li, T. Shi, Q. Xia, Eliminate localized eigenmodes in level set based topology optimization for the maximization of the first eigenfrequency of vibration, Adv. Eng. Softw. 107 (2017) 59–70. doi: 10.1016/j.advengsoft.2016.12.001. [7] P. Wei, M.Y. Wang, Piecewise constant level set method for structural topology optimization, Int. J. Numer. Methods Eng. 78(4) (2009) 379–402. doi: 10.1002/nme.2478. [8] L. Xia, Q. Xia, X. Huang, Y.M. Xie, Bi-directional evolutionary structural optimization on advanced structures and materials: a comprehensive review, Arch. Comput. Methods Eng. 25(2) (2018b) 437–478. doi: 10.1007/s11831-016-9203-2. [9] P. Duysinx, M.P. Bendsøe, Topology optimization of continuum structures with 29

local stress constraints, Int. J. Numer. Methods Eng. 43 (8) (1998) 1453–1478. doi: 10.1002/(SICI)1097-0207(19981230)43:8<1453::AID-NME480>3.0.CO;2-2. [10] S.H. Jeong, S.H. Park, D.H. Choi, G.H. Yoon, Topology optimization considering static failure theories for ductile and brittle materials, Comput. Struct. 110 (2012) 116-132.

doi: 10.1016/j.compstruc.2012.07.007.

[11] Y. Luo, Z. Kang, Topology optimization of continuum structures with drucker-prager yield stress constraints, Comput. Struct. 90 (2012) 65-75. doi: 10.1016/j.compstruc.2011.10.008. [12] S. J. Moon, G. H. Yoon, A newly developed qp-relaxation method for element connectivity parameterization to achieve stress-based topology optimization for geometrically nonlinear structures, Comput. Meth. Appl. Mech. Eng. 265 (2013) 226-241. doi: 10.1016/j.cma.2013.07.001. [13] A. Takezawa, G.H. Yoon, S.H. Jeong, M. Kobashi, M. Kitamura, Structural topology optimization with strength and heat conduction constraints, Comput. Meth. Appl. Mech. Eng. 276 (2014) 341-361. doi:

10.1016/j.cma.2014.04.003.

[14] C. Le, J. Norato, T. Bruns, C. Ha, D. Tortorelli, Stress-based topology optimization for continua, Struct. Multidiscip. Optim. 41 (4) (2010) 605–620. doi: 10.1007/s00158-009-0440-y. [15] H. Svärd, Interior value extrapolation: a new method for stress evaluation during topology optimization, Struct. Multidiscip. Optim. 51 (3) (2015) 613–629. doi: 10.1007/s00158-014-1171-2. [16] M. Bruggi, Topology optimization with mixed finite elements on regular grids, Comput.

Methods

Appl.

Mech.

Eng.

305

(2016)

133–153.

doi:

10.1016/j.cma.2016.03.010. [17] G. Cheng, X. Guo, Epsilon-relaxed approach in structural topology optimization, Struct. Multidiscip. Optim. 13

(4) (1997) 258–266. doi:

10.1007/BF01197454. [18] M. Bruggi, On an alternative approach to stress constraints relaxation in topology optimization, Struct. Multidiscip. Optim.36 (2008) 125–141. doi: 10.1007/s00158-007-0203-6. [19] M. Bruggi, P. Duysinx, Topology optimization for minimum weight with compliance and stress constraints, Struct. Multidiscip. Optim. 46 (3) (2012) 369–384. doi: 10.1007/s00158-012-0759-7. [20] R. Yang, C. Chen, Stress-based topology optimization, Struct. Multidiscip. 30

Optim. 12 (2) (1996) 98–105. doi: 10.1007/BF01196941. [21] L. Xia, L. Zhang, Q. Xia, T.L. Shi, Stress-based topology optimization using bi-directional evolutionary structural optimization method, Comput. Methods Appl. Mech. Engrg. 333 (2018) 356-370. doi: 10.1016/j.cma.2018.01.035. [22] F. Zhao, L. Xia, W.X. Lai, Q. Xia, T.L. Shi, Evolutionary topology optimization of continuum structures with stress constraints, Struct. Multidiscip. Optim. 59 (2019) 647–658. doi: 10.1007/s00158-018-2090-4. [23] Y.M. Xie, G.P. Steven, A simple evolutionary procedure for structural optimization,

Comput.

Struct.

49

(1993)

885–896.

doi:

10.1016/0045-7949(93)90035-c. [24] X. Huang, Y.M. Xie, Convergent and mesh-independent solutions for bi-directional evolutionary structural optimization method, Finite Elem. Anal. Des. 43 (2007) 1039–1049. doi: 10.1016/j.finel.2007.06.006. [25] F. Fritzen, L. Xia, M. Leuschner, P. Breitkopf, Topology optimization of multiscale elastoviscoplastic structures, Inter. J. Numer. Methods Eng. 106 (6) (2016) 430–453. doi: 10.1002/nme.5122. [26] L. Xia, F. Fritzen, P. Breitkopf, Evolutionary topology optimization of elastoplastic structures, Struct. Multidiscip. Optim. 55 (2) (2017) 569–581. doi: 10.1007/s00158-016-1523-1. [27] L. Xia, D. Da, J. Yvonnet, Topology optimization for maximizing the fracture resistance of quasi-brittle composites, Comput. Methods Appl. Mech. Engrg. 332 (2018) 234–254. doi: 10.1016/j.cma.2017.12.021. [28] T. Bruns, D. Tortorelli, An element removal and reintroduction strategy for the topology optimization of structures and compliant mechanisms, Int. J. Numer. Meth. Eng. 57(10) (2003) 1413–1430. doi: 10.1002/nme.783. [29] T. Buhl, C. Pedersen, O. Sigmund, Stiffness design of geometrically nonlinear structures using topology optimization, Struct. Multidiscip. Optimiz. 19(2) (2000) 93–104. doi: 10.1007/s001580050089. [30] H. Gea, J. Luo, Topology optimization of structures with geometrical nonlinearities,

Comput.

Struct.

79(20–21)

(2001)

1977–1985.

doi:

10.1016/S0045-7949(01)00117-1. [31] Y. Luo, M. Wang, Z. Kang, Topology optimization of geometrically nonlinear structures based on an additive hyperelasticity technique, Comput. Methods Appl. Mech. Eng. 286 (2015) 422–441. doi: 10.1016/j.cma.2014.12.023. 31

[32] C. Pedersen, T. Buhl, O. Sigmund, Topology synthesis of large-displacement compliant mechanisms, Int. J. Numer. Meth. Eng. 50(12) (2001) 2683–2705. doi: 10.1002/nme.148. [33] F. Wang, B. Lazarov, O. Sigmund, J. Jensen, Interpolation scheme for fictitious domain techniques and topology optimization of finite strain elastic problems, Comput. Methods Appl. Mech. Eng. 276 (2014a) 453–472. doi: 10.1016/j.cma.2014.03.021. [34] G. Yoon, Y. Kim, Element connectivity parameterization for topology optimization of geometrically nonlinear structures, Int. J. Solids Struct. 42(7) (2005) 1983–2009. doi: 10.1016/j.ijsolstr.2004.09.005. [35] M. Bendsøe, J. Guedes, S. Plaxton, J. Taylor, Optimization of structure and material properties for solids composed of softening material, Int. J. Solids Struct. 33(12) (1996) 1799–1813. doi: 10.1016/0020-7683(95)00121-2. [36] K. Maute, S. Schwarz, E. Ramm, Adaptive topology optimization of elastoplastic

structures,

Struct.

Optimiz.

15(2)

(1998)

81–91.

doi:

10.1007/BF01278493. [37] S. Schwarz, K. Maute, E. Ramm, Topology and shape optimization for elastoplastic structural response, Comput. Methods Appl. Mech. Eng. 190(15–17) (2001) 2135–2155. doi: 10.1016/S0045-7825(00)00227-9. [38] G. Yoon, Y. Kim, Topology optimization of material-nonlinear continuum structures by the element connectivity parameterization, Int. J. Numer. Meth. Eng. 69(10) (2007) 2196–2218. doi: 10.1002/nme.1843. [39] K. Yuge, N. Kikuchi, Optimization of a frame structure subjected to a plastic deformation, Struct. Optimiz. 10(3–4) (1995) 197–208. doi: 10.1007/BF01742592. [40] K. Yuge, N. Iwai, N. Kikuchi, Optimization of 2-D structures subjected to nonlinear deformations using the homogenization method, Struct. Optimiz. 17(4) (1999) 286–299. doi: 10.1007/BF01207005. [41] X. Huang, Y.M. Xie, Topology optimization of nonlinear structures under displacement

loading,

Eng.

Struct.

30(7)

(2008)

2057–2068.

doi:

10.1016/j.engstruct.2008.01.009. [42] X. Huang, Y.M. Xie, G. Lu, Topology optimization of energy-absorbing structures,

Int.

J.

Crash

worthiness

12(6)

(2007)

663–675.

doi:

10.1080/13588260701497862. [43] D. Jung, H. Gea, Topology optimization of nonlinear structures, Finite Elem. 32

Anal. Des. 40(11) (2004) 1417–1427. doi: 10.1016/j.finel.2003.08.011. [44] T. Buhl, C.B.W. Pedersen, O. Sigmund, Stiffness design of geometrically nonlinear structures using topology optimization, Struct. Multidiscip. Optimiz. 19 (2000) 93–104. doi: 10.1007/s001580050089. [45] T.E. Burns, D.A. Tortorelli, An element removal and reintroduction strategy for the topology optimization of structures and compliant mechanisms, Int. J. N um. Methods Eng. 57 (2003) 1413–1430. doi: 10.1002/nme.783. [46] M.A. Crisfield, Non-linear finite element analysis of solids and structures, (1991) New York: Wiley. [47] P. Duysinx, O. Sigmund, New development in handling stress constraints in optimal

material

distribution,

in:

Proc.

7th

AIAA/USAF/NASA/ISSMO

Symposium on Multidisciplinary Analysis and Optimization. A Collection of Technical Papers (Held in St. Louis, Missouri), Vol. 3, 1998, pp. 1501–1509. [48] X. Huang, Y.M. Xie, Evolutionary topology optimization of continuum structures with an additional displacement constraint, Struct. Multidiscip. Optim. 40 (2010) 409–416. doi: 10.1007/s00158-009-0382-4. [49] Z.H. Zuo, Y.M. Xie, X. Huang, Evolutionary topology optimization of structures with multiple displacement and frequency constraints, Adv. Struct. Eng. 15(2) (2012) 359–372. doi: 10.1260/1369-4332.15.2.359. [50] O. Sigmund, A 99 line topology optimization code written in Matlab, Struct. Multidiscip. Optim. 21 (2) (2001) 120–127. doi: 10.1007/s001580050176.

33

Figures:

Force

F*

O

B

A WC

Displacement

U*

Fig. 1. Load-deflection curves in nonlinear finite element analysis.

34

Start

0 Input the information of initial approximate solution d ,

convergence precision  d or  R , and external load P

The tangential stiffness matrix of the structure is obtained by n calculating  R d d dn  K T ,  n  0,1, 2,...

 

n n n n Solve the equations K T d   R d , and obtain the d

d n 1  d n  d n , and calculate R  d n1 

No

n Convergence criterion d

or R  dn1 

dn   d

P   R satisfied? Yes End

Fig. 2. Flow chart of iterative solution procedure of nonlinear equations.

35

Start

End

Define design domain, load case, boundary and FE mesh

Visualization Yes

Set the initial topology full of solid material

No

Convergence criterion and Vreq satisfied?

Define BESO parameters: Vreq ,

cer , carmax , rsen , rden .

Convert the filtered topology into discrete design

Estimate the number of design iterations

Filter topology variables upon (74)

Calculate the target volume of the current iteration upon (69)

Update topology variables upon (73)

Perform nonlinear FEA iteratively

Filter and average sensitivity numbers by (70) (72)

Determine  iteratively according to (76) (77)

Evaluate sensitivity numbers according to (68)

Fig. 3. Flow chart of the stress-constrained geometrically nonlinear stiffness design problem.

36

80

200

F

80

200 Fig. 4. Design domain, boundary conditions and load of L-bracket benchmark.

(a) W c =7.0914×10-4 (b) W c =7.2068×10-4 (c) W c =7.5715×10-4 (d) W c =7.9779×10-4 max max max max =8.9515 =8.4619 =7.5813 =6.5145  vm  vm  vm  vm Fig. 5. Design results and von Mises stress distribution of the L-bracket. (a) * Geometrically nonlinear stiffness design without stress constraint; (b)  vm =8.5; (c) * * =7.7; (d)  vm =7.2.  vm

37

W c =5.8474×10-4 max =7.4703  vm Iteration 1

W c =5.7948×10-4 max =6.7375  vm Iteration 2

W c =6.4228×10-4 max =6.8427  vm Iteration 10

W c =5.9692×10-4 max =5.7669  vm Iteration 17

W c =6.3853×10-4 W c =6.8202×10-4 W c =7.1404×10-4 W c =7.9779×10-4 max max max max =5.6370 =7.0537 =5.4045 =6.5145  vm  vm  vm  vm Iteration 28 Iteration 36 Iteration 43 Iteration 68 Fig. 6. The topological evolution with stress field for the case of geometrically * nonlinear stiffness design with stress constraint  vm =7.2.

(a) W c =7.1448×10-4 (b) W c =7.0914×10-4 (c) W c =9.4923×10-4 (d) W c =7.9779×10-4 max max max max =9.6499 =8.9515 =6.2665 =6.5145  vm  vm  vm  vm Fig. 7. Topological designs and von Mises stress distribution of the L-bracket benchmark under different cases. (a) Linear stiffness design without stress constraint; (b) Geometrically nonlinear stiffness design without stress constraint; (c) Geometrically nonlinear stress design; (d) Geometrically nonlinear stiffness design * with stress constraint  vm =7.2.

38

(a) (b) (c) (d) Fig. 8. Sensitivity distribution (3th iteration) of the L-bracket benchmark under different cases. (a) Linear stiffness design without stress constraint; (b) Geometrically nonlinear stiffness design without stress constraint; (c) Geometrically nonlinear stress * design; (d) Geometrically nonlinear stiffness design with stress constraint  vm =7.2.

(a) W c =7.9779×10-4 (b) W c =7.2901×10-4 (c) W c =6.0848×10-4 max max max =6.5145 =5.7595 =7.0467  vm  vm  vm Fig. 9. Comparison of the influence of volume fraction on topological structure and von Mises stress distribution of the L-bracket benchmark under the case of * geometrically nonlinear stiffness design with stress constraint  vm =7.2.

(a) Vreq =50%; (b) Vreq =60%; (c) Vreq =70%.

39

F

120 30 360 Fig. 10. MBB beam design with one pre-existing crack notch.

(a) W c =2.5022×10-4

max =7.2902 (b) W c =2.5318×10-4  vm

max =6.4547  vm

max max (c) W c =2.6791×10-4  vm =5.9822 (d) W c =2.7764×10-4  vm =4.5759 Fig. 11. Design results and von Mises stress distribution of the MBB beam for different maximum von Mises stress constraints. (a) Geometrically nonlinear stiffness * * * design without stress constraint; (b)  vm =6.5; (c)  vm =6.0; (d)  vm =5.8.

40

max W c =1.7963×10-4  vm =4.2215 Iteration 1

max W c =1.7018×10-4  vm =5.3760 Iteration 2

max W c =2.0044×10-4  vm =3.8802 Iteration 11

max W c =1.8818×10-4  vm =4.9384 Iteration 15

max W c =1.9972×10-4  vm =4.6861 Iteration 30

max W c =2.3949×10-4  vm =3.9670 Iteration 42

max max W c =2.8229×10-4  vm W c =2.7764×10-4  vm =4.5887 =4.5759 Iteration 72 Iteration 100 Fig. 12. The topological evolution with stress field for the case of geometrically * nonlinear stiffness design with stress constraint  vm =5.8.

41

(a) W c =2.4271×10-4

max =7.7037 (b) W c =2.5022×10-4  vm

max =7.2902  vm

max max (c) W c =3.8641×10-4  vm =4.0414 (d) W c =2.7764×10-4  vm =4.5759 Fig. 13. Topological designs and von Mises stress distribution of the MBB beam

under different cases. (a) Linear stiffness design without stress constraint; (b) Geometrically nonlinear stiffness design without stress constraint; (c) Geometrically nonlinear stress design; (d) Geometrically nonlinear stiffness design with stress * constraint  vm =5.8.

42

(a) W c =2.4271×10-4

max =7.7037  vm

(b) W c =2.5022×10-4

max =7.2902  vm

max max (c) W c =3.8641×10-4  vm =4.0414 (d) W c =2.7764×10-4  vm =4.5759 Fig. 14. The evolution histories of the MBB beam under different cases. (a) Linear

stiffness design without stress constraint; (b) Geometrically nonlinear stiffness design without stress constraint; (c) Geometrically nonlinear stress design; (d) Geometrically * nonlinear stiffness design with stress constraint  vm =5.8.

43

(a) W c =2.2658×10-4

max =5.7957  vm

(b) W c =2.7764×10-4

max =4.5759  vm

max max (c) W c =2.9773×10-4  vm =5.7732 (d) W c =3.0665×10-4  vm =5.6184 Fig. 15. Comparison of the influence of p-norm on topological structure and von Mises stress distribution of the MBB beam under the case of geometrically nonlinear * stiffness design with stress constraint  vm 5.8. (a) p=4; (b) p=6; (c) p=8; (d) p=10.

(a) W c =0.6173×10-4

max =3.6089  vm

(b) W c =1.3889×10-4

max =5.4133  vm

max (c) W c =2.7764×10-4  vm =4.5759 Fig. 16. Comparison of the influence of external force on topological structure and von Mises stress distribution of the MBB beam under the case of geometrically * nonlinear stiffness design with stress constraint  vm 5.8. (a) F=25; (b) F=34.

44

60

160

F 100

F 60

200 Fig. 17. Design domain, boundary conditions and load of U-shaped structure.

(a) W c =4.1652×10-4 (b) W c =4.5671×10-4 (c) W c =4.6922×10-4 (d) W c =4.8583×10-4 max max max max =7.4699 =6.7317 =6.5997 =5.0776  vm  vm  vm  vm Fig. 18. Design results and von Mises stress distribution of the U-shaped structure for different maximum von Mises stress constraints. (a) Geometrically nonlinear stiffness * * * design without stress constraint; (b)  vm =7.0; (c)  vm =6.6; (d)  vm =6.2.

45

W c =2.6420×10-4 max =5.5338  vm

Iteration 1

W c =4.3107×10-4 max =5.2183  vm

W c =2.6836×10-4 max =5.5635  vm

W c =2.9022×10-4 max =5.7829  vm

W c =4.1206×10-4 max =5.5519  vm

Iteration 6

Iteration 20

Iteration 42

W c =5.3323×10-4 max =5.4902  vm

W c =4.8942×10-4 max =5.7330  vm

W c =4.8583×10-4 max =5.0776  vm

Iteration 50 Iteration 65 Iteration 80 Iteration 98 Fig. 19. The topological evolution with stress field plot for the case of geometrically * nonlinear stiffness design with stress constraint  vm =6.2.

(a) W c =4.2199×10-4 (b) W c =4.1652×10-4 (c) W c =8.4923×10-4 (d) W c =4.8583×10-4 max max max max =7.7158 =7.4699 =5.1665 =5.0776  vm  vm  vm  vm Fig. 20. Topological designs and von Mises stress distribution of the U-shaped structure under different cases. (a) Linear stiffness design without stress constraint; (b) Geometrically nonlinear stiffness design without stress constraint; (c) Geometrically nonlinear stress design; (d) Geometrically nonlinear stiffness design with stress * constraint  vm =6.2. 46

(a) W c =4.2199×10-4

max =7.7158  vm

(b) W c =4.1652×10-4

max =7.4699  vm

max max (c) W c =8.4923×10-4  vm =5.1665 (d) W c =4.8583×10-4  vm =5.0776 Fig. 21. The evolution histories of the U-shaped structure under different cases. (a) Linear stiffness design without stress constraint; (b) Geometrically nonlinear stiffness design without stress constraint; (c) Geometrically nonlinear stress design; (d) * Geometrically nonlinear stiffness design with stress constraint  vm =6.2.

47