Reliability-based topology optimization design method for part structure considering bounded hybrid uncertainties

1. A reliability-based topology optimization design method for a part structure considering bounded hybrid uncertainties, comprising:

step 1) considering following uncertainties in manufacture and service processes of the part structure: regarding an amplitude and a direction of an external load with insufficient sample information as interval uncertainties, regarding a material property of the part structure with sufficient sample information as a bounded probabilistic uncertainty, and describing a bounded probabilistic uncertainty parameter as a random variable subject to generalized beta distribution;

step 2) discretizing a design domain of the part structure, comprising:

simplifying a force condition of the part structure into a two-dimensional plane stress state, retaining installation holes and removing structural details to improve calculating efficiency, placing the simplified part structure in a regular rectangular design domain, dividing the rectangular design domain into N x ×N y square elements, where N x and N y are numbers of divisions along x, y axes, respectively; and imparting, based on a solid isotropic material with penalization (SIMP) topology optimization framework, each element with a unique design variable ρ e ∈[0,1] (e=1,2, . . . , N x ·N y );

step 3) imposing physical constraints and geometric constraints on the discretized structure, comprising:

step 3.1) imposing the physical constraints including fixing or supporting and an external load according to a classical finite element method, which comprises setting the zero displacement at nodes with a constrained degree of freedom, and specifying the node where an external load is applied; and

step 3.2) the geometric constraints including specified holes in the structure and areas where materials are forcibly retained, setting design variables corresponding to elements in the holes as ρhd e≡0, setting design variables corresponding to elements in the areas where the materials are forcibly retained as ρ e ≡1, and keeping values of the design variables corresponding to elements in the holes and the areas where the materials are to be retained unchanged in subsequent optimization process;

step 4) establishing, by taking a space utilization rate of the design domain as an objective function and taking displacements of several key points of the part structure under a joint influence of interval and bounded probabilistic hybrid uncertainties as reliability constraint performance, a reliability-based topology optimization design model for the part structure, as shown in Eq.1:

min

ρ

V

⁡

(

ρ

)

=

∑

e

=

1

N

e

ρ

e

/

V

0

Eq

.

1

s

.

t

.

g

q

(

ρ

,

X

,

I

)

=

-

P

⁡

(

u

q

(

ρ

,

X

,

I

)

≤

u

qcri

)

+

P

q

≤

0

(

q

=

1

,

2

,

…

,

N

c

⁢

o

⁢

n

)

K

⁡

(

ρ

,

X

)

⁢

U

=

F

⁡

(

I

)

0

≤

ρ

min

≤

ρ

e

≤

1

⁢

(

e

=

1

,

2

,

…

,

N

e

)

X

=

(

X

1

,

X

2

,

…

,

X

m

)

T

,

I

=

(

f

1

,

f

2

,

…

,

f

n

,

α

1

,

α

2

,

…

,

α

n

)

T

where in Eq. 1, ρ=(ρ 1 , ρ 2 , . . . , ρ N x ·N y ) T is a design vector composed of design variables ρ e (e=1,2, . . . , N x ·N y ), and ρ min is a minimum allowable value for the design variables; a total number of elements is N e =N x ·N y ; a bounded probabilistic uncertainty vector X=(X 1 , X 2 , . . . , X m ) T comprises m uncertain material properties of the part structure; an interval uncertainty vector I=(f 1 , f 2 , . . . , f n , α 1 , α 2 , . . . , α n ) T comprises amplitudes f 1 , f 2 , . . . , f n and direction angles α 1 , α 2 , . . . , α n of n uncertain external loads on the part structure;

where V(ρ) is the space utilization rate of the design domain, corresponding to a total material usage of the part structure, and V 0 is a volume of the design domain; where g q (ρ,X,I) is a q th constraint function; u q (ρ,X,I) is a displacement at a q th key point regarded as a q th constraint performance5, denoted as u q for short, and u qcri is an allowable maximum value of u q (ρ, X, I); P(⋅) calculates an occurrence probability of an event in a bracket, P q is a reliability index of the q th constraint performance, and N con is a number of constraint functions; and

where in an equilibrium equation K(ρ,X)U=F(I) of the part structure, K(ρ,X) is a (2(N x +1)(N y +1))×(2(N x +1)(N y +1))-dimensional global stiffness matrix, which is affected by the design vector ρ and the bounded probabilistic uncertainty vector X; F(I) is a (2(N x +1)(N y +1))-dimensional nodal force vector affected by the interval uncertainty vector I; U is a (2(N x +1)(N y +1))-dimensional nodal displacement vector; and u q is extracted from U according to u q =L q T U, wherein in a (2(N x +1)(N y +1))-dimensional column vector L q , the element at a position corresponding to the q th key point is 1, and other elements are all 0;

step 5) calculating reliability of a constraint performance under bounded hybrid uncertainties, comprising:

step 5.1) searching a worst working condition of constraint performance u q , comprising:

step 5.1.1) setting X=μ X =(μ X 1 , μ X 2 , . . . , μ X m ) , where μ X 1 , μ X 2 , . . . , μ X m are means of uncertainties X 1 , X 2 , . . . , X m , the constraint performance u q is only affected by the interval uncertainty vector; and an uncertain external load F s (s=1,2, . . . , n) is rewritten as its components in horizontal and vertical directions such as F s =[f s cos α s , f s sin α s ] T ;

step 5.1.2) calculating, based on a linear elastic hypothesis in small deformation, displacement U s caused by F s according to Eq.2 based on e s x =[1 0] T , e s y =[0 1] T and F s , where e s x and e s y are respectively element nodal forces in the horizontal and vertical directions at an point where F s exerts;

U s =U s x +U s y =u s x f s cos α s +u s y f s sin α s =[u s x ,u s y ]·F s   Eq.2

where u s x =[u s x 0] T and u s y =[0 u s y ] T are respectively nodal displacement vectors calculated through the equilibrium equation of the part structure when only the element nodal force e s x or e s y acts;

step 5.1.3) calculating gradients of the constraint performance u q with respect to an amplitude and a direction of an uncertain external load according to Eq.3 and Eq.4:

∂

u

q

∂

f

s

=

L

q

T

⁢

∂

U

∂

f

s

=

L

q

T

⁢

∑

i

=

1

n

∂

U

i

∂

f

s

=

L

q

T

·

[

u

s

x

,

u

s

y

]

·

[

cos

⁢

α

s

,

sin

⁢

α

s

]

T

Eq

.

3

∂

u

q

∂

α

s

=

L

q

T

⁢

∂

U

∂

α

s

=

L

q

T

⁢

∑

i

=

1

n

∂

U

i

∂

α

s

=

L

q

T

·

[

u

s

x

,

u

s

y

]

·

[

-

f

s

⁢

sin

⁢

α

s

,

f

s

⁢

cos

⁢

α

s

]

T

;

Eq

.

4

step 5.1.4) solving, by utilizing results of Eq.3 and Eq.4, a worst working condition Ĩ q by a gradient search algorithm, and an external load of the worst working condition being F={tilde over (F)} q ;

step 5.2) restoring μ X to X under the worst working condition Ĩ q , so that constraint performance ũ q =u q (ρ,X,Ĩ q ) under the worst working condition is only manifested as a probabilistic type, and solving a performance fluctuation under the worst working condition to evaluate reliability, comprising:

step 5.2.1) calculating, according to Eq.5, a gradient of the constraint performance ũ q with respect to a bounded probabilistic uncertainty parameter X i (i=1 2, . . . , m) under the worst working condition:

∂

u

~

q

∂

X

i

=

L

q

T

⁢

∂

U

~

∂

X

i

=

-

L

q

T

⁢

K

-

1

⁢

∂

K

∂

X

i

⁢

U

~

=

-

L

q

T

⁢

K

-

1

(

∑

e

=

1

N

e

ρ

e

p

⁢

∂

k

e

∂

X

i

)

⁢

U

~

Eq

.

5

where a summation symbol is a combination operation of the element stiffness matrix defined by a finite element theory, k e is an element stiffness matrix, a nodal displacement vector Ũ under the worst working condition is obtained by solving a governing equation KŨ={tilde over (F)} q under the worst working condition; and p is a penalty factor;

step 5.2.2) searching, based on a result of Eq.5, two bounded probabilistic uncertainty vectors minimizing or maximizing u q (ρ,X,Ĩ q ) respectively according to Eq.6:

{

X

q

L

=

argmin

X

⁢

u

q

(

ρ

,

X

,

I

~

q

)

X

q

R

=

argmax

X

⁢

u

q

(

ρ

,

X

,

I

~

q

)

⁢

(

q

=

1

,

2

,

…

,

N

c

⁢

o

⁢

n

)

Eq

.

6

where corresponding to X q L and X q R , global stiffness matrix K is K q L and K q R , respectively, K q L is denoted as K q L =K(ρ,X q L ) and K q R is denoted as K q R =K(ρ,X q R ), K q L and K q R are uniformly marked as K q * (*=L,R), the nodal displacement vector Ũ is Ũ q L and Ũ q R , respectively, and Ũ q L and Ũ q R are uniformly marked as Ũ q * (*=L,R); and

step 5.2.3) denoting ũ q L =u q (ρ, X q L , Ĩ q ) and ũ q R =u q (ρ, X q R , Ĩ q ), defining [ũ q L , ũ q R ] as a performance fluctuation of the constraint performance u q under the worst working condition, and calculating a reliability {tilde over (R)} q of the constraint performance according to Eq .7:

R

~

q

=

1

2

⁢

tanh

⁢

{

P

·

(

u

qcri

-

u

q

C

)

·

[

1

+

(

P

·

(

u

qcri

-

u

q

C

)

)

γ

]

}

+

1

2

Eq

.

7

where u q C =(ũ q R +ũ q L )/2 is a midpoint of the performance fluctuation under the worst working condition; P is a multiplier calculated by P=1(u q W −ε u ), where ε u is a small constant for adjusting reliability at a boundary position of the performance fluctuation under the worst working condition; u q W =(ũ q R −ũ q L )/2 is a radius of the performance fluctuation under the worst working condition; and γ∈{2i|i∈N + } is a regulatory factor; and

establishing the constraint functions in Eq.1 as Eq.8:

g q (ρ, X,I )=− {tilde over (R)} q +P q ≤0( q= 1, 2, . . . , N con )   Eq.8

step 6) calculating gradients of objective and constraint functions with respect to the design variables:

step 6.1) calculating the gradient of the objective function through Eq.9:

∂

V

⁡

(

ρ

)

∂

ρ

e

=

1

V

0

⁢

(

e

=

1

,

2

,

…

,

N

e

)

Eq

.

9

step 6.2) solving the gradient of a constraint function as follows:

step 6.2.1) writing a gradient expression of g q (ρ,X,I) according to a chain rule, as shown in Eq.10:

∂

g

q

(

ρ

,

X

,

I

)

∂

ρ

e

=

∂

g

q

(

ρ

,

X

,

I

)

∂

u

~

q

L

⁢

∂

u

~

q

L

∂

ρ

e

+

∂

g

q

(

ρ

,

X

,

I

)

∂

u

~

q

R

⁢

∂

u

~

q

R

∂

ρ

e

⁢

(

e

=

1

,

2

,

…

,

N

e

)

Eq

.

10

step 6.2.2) denoting a function in a bracket of tanh(⋅) in Eq.7 as R(u qcri ), and calculating gradient terms ∂g q (ρ,X,I)/∂ũ q L and ∂g q (ρ,X,I)/∂ũ q R in Eq.10 according to Eq.11:

∂

g

q

(

ρ

,

X

,

I

)

∂

u

~

q

=

sec

⁢

h

2

(

R

⁡

(

u

qcri

)

)

2

·

∂

R

⁡

(

u

qcri

)

∂

u

~

q

*

⁢

(*

=

L

,

R

)

Eq

.

11

where a gradient term ∂R(u qcri )/∂ũ q * is as follows:

∂

R

⁡

(

u

qcri

)

∂

u

~

q

L

=

u

qcri

-

u

q

C

2

⁢

P

2

-

P

2

+

(

u

qcri

-

u

q

C

)

γ

+

1

2

⁢

P

2

-

P

⁡

(

γ

+

1

)

2

⁢

(

u

qcri

-

u

q

C

)

γ

Eq

.

12

∂

R

⁡

(

u

qcri

)

∂

u

~

q

R

=

-

u

qcri

-

u

q

C

2

⁢

P

2

-

P

2

-

(

u

q

⁢

c

⁢

r

⁢

i

-

u

q

C

)

γ

+

1

2

⁢

P

2

-

P

⁡

(

γ

+

1

)

2

⁢

(

u

q

⁢

c

⁢

r

⁢

i

-

u

q

C

)

γ

Eq

.

13

step 6.2.3) giving gradient terms ∂ũ q L /∂ρ e in Eq.10 in a uniform form according to the SIMP framework:

∂

u

~

q

*

∂

ρ

e

=

-

L

q

T

(

K

q

*

)

-

1

⁢

∂

K

q

*

∂

ρ

e

⁢

U

~

q

*

⁢

(*

=

L

,

R

)

Eq

.

14

where K q * and Ũ q * (*=L,P) are defined in 5.2.2; and a gradient term ∂K q * /∂ρ e is calculated according to Eq.15:

∂

K

q

*

∂

ρ

e

=

p

⁢

ρ

e

p

-

1

⁢

〈

k

e

⁢

q

*

〉

⁢

(*

=

L

,

R

)

Eq

.

15

where k eq * is an element stiffness matrix extracted from K q * , and k eq * is a square matrix reconstructed by performing a combined operation on elements in k eq * according to the element stiffness matrix, and is consistent with K q * in dimensionality; and

step 6.2.4) substituting all the gradient terms in Eq.11 to Eq.15 into Eq.10 to obtain a gradient of the constraint function g q (ρ,X,I); and

step 7) updating the design variables by using a moving asymptote algorithm based on the gradients of the objective and constraint functions with respect to the design variables, checking a difference value between an objective function value in a current iteration and an objective function value in a previous iteration, wherein the difference value for a first iteration is defined as an objective function value, and when the difference value is less than a convergence threshold, outputting the updated design variables, or otherwise, repeating the steps 5) to 7).