An Efficient Approach to Structural Optimization of Aircraft Structures ...

52nd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference
19th 4 - 7 April 2011, Denver, Colorado

AIAA 2011-2153

An efficient approach to structural optimization of aircraft structures considering parameter variation Aitor Baldomir1 and Santiago Hernández2 University of Coruña, 15071 La Coruña, Spain Félix Nieto3 University of Coruña, 15071 La Coruña, Spain and Jacobo Díaz4 University of Coruña, 15071 La Coruña, Spain

Nomenclature α β ε λ σ c F g i j p ∆p r S w X X* X*aprox

O

= = = = = = = = = = = = = = = = = =

magnitude of the change in design variables new design variable (Vanderplaats’ formulation) relative error buckling load factor normal stress positive constant (Vanderplaats’ formulation) objective function design constraint design variable index active constraint index parameter parameter increase inactive constraint index vector of feasible directions vertical displacement vector of design variables vector of optimum design variables vector of approximate optimum design variables

I. Introduction

ptimum design is currently very well established as a relevant design tool in some engineering fields as aerospace and car industries. The usual formulation of optimization problems was proposed by Schmit1 in 1960, exactly half a century ago. (1) F ( X, p) subject to (2) g j ( X, p)  0 j  1,..., m Where vectors X, p contain the set of design variables and fixed parameters, gj (j = 1,...,m) are the problem constraints and F is the objective function meaning the property of the design to be optimized. Consequently, vector X* of optimum values of design variables is X* = X (p), in other words it depends of the values designed for the fixed parameters of the problem. Similarly the optimum of the objective function F* is also F* = F (p). 1

Assistant Professor, Structural Mechanics Group, E.T.S.I. Caminos, Canales y Puertos. Professor, Structural Mechanics Group, E.T.S.I. Caminos, Canales y Puertos, Associate Fellow AIAA. 3 Assistant Professor, Structural Mechanics Group, E.T.S.I. Caminos, Canales y Puertos. 4 Assistant Professor, Structural Mechanics Group, E.T.S.I. Caminos, Canales y Puertos, Member AIAA. 1 American Institute of Aeronautics and Astronautics 2

Copyright © 2011 by Aitor Baldomir, Santiago Hernandez, Felix Nieto and Jacobo Diaz. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

It is not uncommon in real life problems that after finishing an optimization loop of a prototype, some of the parameters of the problem that were assumed to be constant must be modified due to manufacturing problems, design requirements or managerial decisions. The amount of changes may not be dramatic but variations of 20% 30% can occur to the authors’ experience. In such circumstance the obvious solution is to proceed again with an optimization process aimed to find out the best design for the new set of parameter values. Nevertheless this straight approach may be expensive in terms of time and computer cost, as real optimization problems in industry may require several hours of computer time and the use of expensive computational facilities. Also it seems unwise not take advantage of the existing information of the problem. Because of that, the values of the optimum solution for the previous values of the parameters could be a valuable starting point for the new solution. Therefore some techniques were formulated in the past to produce an approximate solution of the new optimal design using the previous one and a linear approximation by means of the first order derivatives of the design variables and the objective function with regards to the changing parameters. Amongst others, Sobieski2, Barthelemy3,4, Vanderplaats5 and Hernandez6 described different approaches for that purpose. After obtaining X* / p and F ( X* ) / p the new optimum could be approximate to

X* F ( X* ) (3) p and F ( X* )  p p p No activity has taken place in this topic in the past years but nevertheless its relevance is maintained by the fact that the size and complexity of the optimization problems being solved in industry is enlarging nowadays. In this abstract this optimization topic has been revisited and two main outputs can be mentioned: 1) Some inconveniences of the previous procedures have been eliminated. 2) The pre-existing techniques were only applied to simple structures. The advantages of the new technique are proven by comparing its performance with the results produced by the previous procedures in a pair of truss structures. In addition to that an aluminum panel similar to those used in the construction of aircraft fuselage is used to demonstrate the capabilities of the new approach. X  X* 

II. Sensitivity analysis with regards to fixed parameter using the feasible directions method This approach was presented by Vanderplaats and Yoshida7. In this procedure the varying parameter p was considered an additional design variable and the new vector of design variables was obtained as a combination of the former X* and a new direction S by

X  X*    S

i  1,..., n

(4)

and

X n i  p*    S n 1

(5)

The direction S was the output of an optimization problem formulated as

min F ( X* )T  S  c 

(6)

g j ( X* )T  S  0

(7)

subject to

jJ

 sn 1    0 p  0

sn 1    0 p  0

(8) (9)

In this formulation c is a positive term. As the aim is to minimize the objective function, β will be positive while T * F  X*   S would be negative or at least a small positive number. In the former case the objective function F(X ) would be decreased, in the latter would be increased the least as possible.

2 American Institute of Aeronautics and Astronautics

Figure 1.Design region and possible directions for vector S

III. Some enhancements of the approximate solution Constraint set tracking The formulation described in the previous paragraph allows that an active constraint at the optimum becomes passive when parameter p varies, but does not take in consideration if a new constraint turns up active. When this happen the direction S changes suddenly and the error of the approximate optimum with regards to the exact one can be large, because of that a procedure to improve the former formulation tracking constraints behavior has been explored. The main idea is to produce a linear expansion of the subset of passive constraints at the optimum point X*. n

g r ( X)  g r ( X* )   i 1

g r ( X* )  xi ; r  R xi

(10)

Imposing that the constraint becomes active g r ( X)  0

(11)

then n

g r ( X )  g r ( X* )   i 1

g r ( X* )  xi  0 r  R xi

(12)

recalling that

xi  xi*    si  xi*  xi

(13)

Substituting in the previous expression n

g r ( X* )   r   i 1

g r ( X* )  si  0; xi

r 

 g r ( X* ) g r ( X* )  si  xi i 1 n

rR

(14)

The lower value of r will be the one to turn up a passive constraint into an active one. Therefore the amount of parameter increment that will produce such event would be pact

pact   act  sn 1

(15)

If pact  p will not be necessary to update the set of active constraints. On the contrary if pact  p a step by step procedure is required along with an optimization procedure with different set of active constraints. The flowchart of this approach appears in figure 2.


Figure 2.Flowchart of the procedure

Evaluation of violated constraints The approximate solution can led to some amount of constraints violation and that circumstance decreases the quality of the solution. Therefore a screening of the values of the constraints is required in order for the procedure to be appropriate. In that regards an upper limit of the error  has been defined. If after finding the vector X any constraint is violated more than  the new values contained in X are rejected and instead a new direction S is evaluated. This strategy has proven very effective and provided a high level of accuracy in the approximate solutions as shown in the examples. In return it requires to solve in more occasions the problem defined by expression (6-9) but this is not an important drawback because the computer time needed by the procedure is very small.

IV. Application examples 1. Three bar truss This structure, shown in figure 3, has two design variables and is loaded with a force Q. The objective of the optimization is to minimize the volume subject to stress constraints and the vertical displacement of node 1. Therefore the problem can be written as

min F  (2 2  x1  x2 )  L

(16)

subject to xi  1 10 5 i  1, 2

(17)

g   1   máx

(18)

g w  w1  wmin

(19)

When max = 270 MPa the solution is x1  10.516 cm 2 ,

x2  5.443 cm3

and F *  3.520 cm3


L=1m Q = 360 KN E = 210000 MPa

Figure 3.Three bar truss

Let assume that different values for max are adopted, in other words max is selected as the variable parameter p. Table 1 shows both the values of the optimum obtained by a new optimization process and those corresponding to the approximate solution produced by the procedure described in paragraph 3 without tracking the set of active constraints. The numerical values are presented up to a maximum value of p  30% . It can be seen that up to an increment about 6% errors on constraints, design variables and objective function are small. But afterwards the error increases and this is because at p  6.25% the displacement constraints becomes active and therefore the procedure for finding direction S must take this fact in account. Δp (%) 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30

p=σmáx 270.0 275.4 280.8 286.2 291.6 297.0 302.4 307.8 313.2 318.6 324.0 329.4 334.8 340.2 345.6 351.0

x1*

x2*

10.516 10.310 10.111 9.920 9.646 9.375 9.119 8.876 8.646 8.427 8.219 8.021 7.833 7.653 7.481 7.317

5.443 5.337 5.234 5.135 5.301 5.493 5.674 5.846 6.009 6.163 6.310 6.450 6.583 6.710 6.832 6.948

x 1 * aprox x 2 * aprox ε x1 * (%) ----10.294 10.073 9.852 9.631 9.410 9.188 8.967 8.746 8.525 8.304 8.082 7.861 7.640 7.419 7.198

----5.365 5.287 5.209 5.131 5.052 4.974 4.896 4.818 4.740 4.661 4.583 4.505 4.427 4.349 4.270

-----0.1 -0.4 -0.7 -0.2 0.4 0.8 1.0 1.2 1.2 1.0 0.8 0.4 -0.2 -0.8 -1.6

ε x2 * (%)

ε F* (%)

ε gσ (%)

ε gw (%)

----0.5 1.0 1.4 -3.2 -8.0 -12.3 -16.2 -19.8 -23.1 -26.1 -28.9 -31.6 -34.0 -36.3 -38.5

----0.0 0.2 0.4 0.7 1.1 1.6 2.3 3.1 4.0 5.0 6.2 7.5 8.9 10.5 12.2

----0.0 0.2 0.4 0.6 1.0 1.5 2.0 2.6 3.4 4.2 5.1 6.2 7.3 8.6 10.0

----------------1.5 3.6 5.7 7.9 10.2 12.6 15.1 17.7 20.5 23.3 26.3 29.5

Table 1. Optimum X* and approximate X*aprox

Figure 4. Graphs of real and approximate optima

If the described procedure of tracking the set of passive constraints is carried out the numerical results are better than the previous ones. This can be observed in table 2 and graphical data appear in figure 5. 5 American Institute of Aeronautics and Astronautics

Δp (%)

p=σmáx

6.86 8 10 12 14 16 18 20 22 24 26 28 30

288.5 291.6 297.0 302.4 307.8 313.2 318.6 324.0 329.4 334.8 340.2 345.6 351.0

x1*

x2*

9.809 9.646 9.375 9.119 8.876 8.646 8.427 8.219 8.021 7.833 7.653 7.481 7.317

5.186 5.301 5.493 5.674 5.846 6.009 6.163 6.310 6.450 6.583 6.710 6.832 6.948

x 1 * aprox x 2 * aprox ε x1 * (%) 9.758 9.515 9.229 8.944 8.658 8.373 8.087 7.802 7.516 7.230 6.945 6.659 6.374

5.175 5.347 5.549 5.751 5.953 6.155 6.357 6.558 6.760 6.962 7.164 7.366 7.568

-0.5 -1.4 -1.6 -1.9 -2.5 -3.2 -4.0 -5.1 -6.3 -7.7 -9.3 -11.0 -12.9

ε x2 * (%)

ε F* (%)

ε gσ (%)

ε gw (%)

-0.2 0.9 1.0 1.4 1.8 2.4 3.1 3.9 4.8 5.8 6.8 7.8 8.9

0.5 1.0 1.1 1.3 1.7 2.1 2.6 3.3 4.0 4.8 5.8 6.8 8.0

0.5 1.0 1.2 1.4 1.8 2.3 2.9 3.7 4.7 5.8 7.2 8.7 10.5

0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4

Table 2. Real and approximate optima tracking constraints behaviour

Figure 5. Real and approximate optima tracking constraints behaviour

Better approximation may be obtained by including the additional criterion limiting the numerical error in constraints. By doing that and defining a maximum error value of 2% the numerical results of table 3 and the graphs of figure 6 are obtained. Δp (%) 14 16 18 20 22 24 26 28 30

p=σmáx 307.8 313.2 318.6 324.0 329.4 334.8 340.2 345.6 351.0

x1*

x2*

8.876 8.646 8.427 8.219 8.021 7.833 7.653 7.481 7.317

5.846 6.009 6.163 6.310 6.450 6.583 6.710 6.832 6.948

x 1 * aprox x 2 * aprox ε x1 * (%) 8.658 8.418 8.193 7.968 7.743 7.518 7.293 7.068 6.843

5.953 6.19 6.349 6.508 6.667 6.826 6.986 7.145 7.304

-2.5 -2.6 -2.8 -3.1 -3.5 -4.0 -4.7 -5.5 -6.5

ε x2 * (%)

ε F* (%)

ε gσ (%)

ε gw (%)

1.8 3.0 3.0 3.1 3.4 3.7 4.1 4.6 5.1

1.7 1.5 1.6 1.8 2.0 2.3 2.7 3.2 3.7

1.8 1.7 1.9 2.1 2.4 2.8 3.4 4.0 4.8

0.4 -0.2 -0.2 -0.2 -0.2 -0.2 -0.2 -0.2 -0.2

Table 3. Real and approximate optima tracking constraints and limiting constraints error


Figure 6. Real and approximate optima tracking constraints and limiting constraints error

2. Ten bar truss

This well known structure contains ten design variables and stress and displacement constraints are included in the formulation of the problem. The objective function is the material volume and optimization can be written as min F   xi  li

(20)

subject to xi  0.1 250   i  250  w5  wmáx

i  1,...,10 i  1,...,10

(21) (22) (23)

Figure 7. Ten bar truss

Optimum values for a set of values of wmax (0.01 cm, 0.015 cm, 0.02 cm, 0.03 cm) appear in Table 4. xi (cm2) x1 x2 x3 x4 x5 x6 x7 x8 x9 x10

p =0.01

p =0.015

p =0.02

p =0.03

28.672 0.100 21.944 14.651 0.100 0.100 7.290 19.711 20.738 0.100

18.966 0.100 15.212 9.859 0.100 0.100 6.083 12.745 13.919 0.100

14.350 0.100 10.847 7.270 0.100 0.100 5.534 9.913 10.293 0.100

9.332 0.100 8.003 4.632 0.100 0.100 5.661 6.557 6.552 0.100

Table 4. Optimum design of different wmax

The allowable value of wmax has been chosen as varying parameter p and the search of the approximate optimum without tracking constraints has been worked out. Tables 5 to 7 show the numerical results. p+ Δp =0.011 Variable x1

p+ Δp =0.0165

p+ Δp =0.022

p+ Δp =0.033

X*

X*_aprox

X*

X*_aprox

X*

X*_aprox

X*

26.006

27.357

17.181

17.989

13.179

13.614

8.318

8.695

x2

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

x3

20.112

20.699

13.955

14.317

9.422

10.099

8.044

7.998

x4

13.354

13.363

8.990

8.916

6.544

6.531

4.121

3.982

X*_aprox

x5

0.100

0.450

0.100

0.321

0.100

0.347

0.100

0.212

x6

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

x7

7.018

6.023

5.778

5.180

5.542

5.223

5.720

5.654

x8

17.791

17.809

11.483

11.307

9.300

8.879

5.834

5.686

x9

18.898

18.916

12.705

12.585

9.259

9.248

5.833

5.633

x10

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

F Error F (%)

43820.549 -

44159.435 0.8

29860.409 -

29879.081 0.1

22920.984 -

23024.102 0.4

16384.567 -

16283.122 -0.6


Table 5. Real and approximate optima for p  10% p+ Δp =0.0125 Variable x1

p+ Δp =0.0188

p+ Δp =0.025

p+ Δp =0.0375

X*

X*_aprox

X*

X*_aprox

X*

X*_aprox

X*

22.829

25.384

15.191

16.177

11.573

12.511

7.938

7.740

x2

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

x3

17.919

18.832

12.010

12.423

8.301

8.976

8.062

7.991

x4

11.775

11.431

7.853

7.071

5.743

5.423

3.938

3.007

X*_aprox

x5

0.100

0.976

0.100

0.100

0.100

0.718

0.100

0.381

x6

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

x7

6.640

4.122

5.531

2.139

5.572

4.757

5.745

5.644

x8

15.509

14.955

10.270

8.801

8.153

7.330

5.569

4.379

x9

16.664

16.184

11.093

9.974

8.128

7.681

5.569

4.253

x10

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

F Error F (%)

38827.291 -

38458.880 -0.9

26470.367 -

29435.616 11.2

20506.263 -

20132.633 -1.8

15931.823 -

14273.718 -10.4

Table 6. Real and approximate optima for p  25%

Table 7. Real and approximate optima for p  50%

It can be observed that the amount of error obtained is excessive. By applying the method of tracking the set of passive constraints the accuracy of the results is improved highly as depicted in tables 8 to 10. p+ Δp =0.011 Variable x1

X* 26.006

X*_aprox 27.357

p+ Δp =0.0165 X* 17.181

p+ Δp =0.022

X*_aprox 17.989

X* 13.179

p+ Δp =0.033

X*_aprox 13.614

X* 8.318

X*_aprox 8.695 0.100

x2

0.100

0.100

0.100

0.100

0.100

0.100

0.100

x3

20.112

20.699

13.955

14.317

9.422

10.099

8.044

7.998

x4

13.354

13.363

8.990

8.916

6.544

6.531

4.121

3.982

x5

0.100

0.450

0.100

0.321

0.100

0.347

0.100

0.212

x6

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

x7

7.018

6.023

5.778

5.180

5.542

5.223

5.720

5.654

x8

17.791

17.809

11.483

11.307

9.300

8.879

5.834

5.686

x9

18.898

18.916

12.705

12.585

9.259

9.248

5.833

5.633

x10

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

F Error F (%)

43820.549 -

44159.435 0.8

29860.409 -

29879.081 0.1

22920.984 -

23024.102 0.4

16384.567 -

16283.122 -0.6

Table 8. Real and approximate optima for p  10% with constraint tracking


p+ Δp =0.0125 Variable x1

p+ Δp =0.0188

p+ Δp =0.025

p+ Δp =0.0375

X*

X*_aprox

X*

X*_aprox

X*

X*_aprox

X*

22.829

25.627

15.191

16.850

11.573

12.511

7.938

7.933

x2

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

x3

17.919

19.058

12.010

13.141

8.301

8.976

8.062

8.060

x4

11.775

11.668

7.853

7.763

5.743

5.423

3.938

3.835

X*_aprox

x5

0.100

0.920

0.100

0.701

0.100

0.718

0.100

0.100

x6

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

x7

6.640

4.421

5.531

4.800

5.572

4.757

5.745

5.755

x8

15.509

15.298

10.270

9.727

8.153

7.330

5.569

5.426

x9

16.664

16.519

11.093

10.956

8.128

7.681

5.569

5.427

x10

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

F Error F (%)

38827.291 -

39190.345 0.9

26470.367 -

26940.484 1.8

20506.263 -

20132.633 -1.8

15931.823 -

15752.070 -1.1

Table 9. Real and approximate optima for p  25% with constraint tracking p+ Δp =0.015 Variable x1

p+ Δp =0.0225

p+ Δp =0.0375

p+ Δp =0.045

X*

X*_aprox

X*

X*_aprox

X*

X*_aprox

X*

18.966

23.521

12.880

14.871

7.938

7.932

7.938

7.933

x2

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

x3

15.212

16.873

9.219

11.106

8.062

8.057

8.062

8.060

x4

9.859

9.531

6.395

5.764

3.938

3.939

3.938

3.835

x5

0.100

1.546

0.100

1.363

0.100

0.100

0.100

0.100

x6

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

x7

6.083

3.781

5.547

4.069

5.745

5.718

5.745

5.755

x8

12.745

12.430

9.086

6.977

5.569

5.572

5.569

5.426

x9

13.919

13.496

9.049

8.128

5.569

5.535

5.569

5.427

x10

0.100

0.100

0.100

0.100

0.100

0.100

0.100

0.100

F Error F (%)

32684.166 -

33776.432 3.3

22473.673 -

21802.542 -3.0

15931.823 -

15898.559 -0.2

15931.823 -

15752.070 -1.1

X*_aprox

Table 10. Real and approximate optima for p  50% with constraint tracking

3. Aircraft fuselage panel Model definition

The next example presented is a curved panel corresponding to a cylindrical fuselage portion of an aircraft. The geometry of the model is defined in Figure 8 and it is composed by three structural elements: the skin or external surface of model, the stringers or longitudinal stiffeners and frames or transversal stiffeners.

Figure 8. Geometric model. Dimensions in millimetres


The separation between stiffeners is indicated in figure and their cross sectional area is a Z-shape geometry with internal dimensions provided in Figure 9.

Figure 9. Stiffener cross sections. Frames (left) and stringers (right). Dimensions in millimeters

A finite element model shown in Figure 10 having 8442 degrees of freedom was created using Abaqus/CAE v6.9-2 and it served to obtain the structural responses in order to verify the feasibility of the design. Skin was modeled by four node shell elements type S4R (Abaqus/CAE8) considering a thickness of 1.5 mm. For the stringers, three node bar elements B31 were used. The material considered was aluminum type AL2024-T3 for the skin and aluminum type AL7075-T6 for the stiffeners, according to C.T. Sun9, with mechanics properties presented in Table 11.

Figure 10. Finite element mesh and boundary conditions (left) and loads for buckling analysis (right)

Property E GPa

σu MPa

σY ρ MPa g/cm3

Material υ Aluminum 2024-T3 72 0.33 449 324 7075-T6 71 0.33 538 490 σu= tensile ultimate stress; σY= tensile yield stress

2.78 2.78

Table 11. Mechanical properties of the material

Boundary conditions are established along the edges of the panel. Firstly, all displacements and rotations are constrained for one of the short edges. Secondly, remaining edges are defined as follows: constraining the normal displacement with respect to the plane defined by the four panel corners and constraining all rotations. A shear buckling analysis of the panel was performed considering the load specified in Figure 10. A unit load was considered so that the eigenvalue obtained is directly the critical buckling load. Besides, a static analysis was performed with geometric nonlinearity including large-displacement formulation for a load case of internal pressure in the fuselage. In this case a maximum pressure of value 0.0621 MPa was considered, taking the data specified in the paper by M. Rouse et al.10.


Optimization problem and numerical results

Before obtaining the sensitivity analysis of the optimal solution, it was necessary to formulate and solve the initial optimization problem. The design variables considered were the skin thickness and the geometric dimensions and thickness of the cross section members of the stiffeners. Design variables are shown in Table 12 and appear in Figure 11. Variable Pe T1 T2 TA TB Te L1 L2 LA LB Le

Description

Min. Initial Max. Thickness of skin 0.4 1.5 3 Length of lower flange in frame cross section 5 25 50 Length of upper flange in frame cross section 5 15 50 Web height of frame cross section 20 60 100 Lip height of frame cross section 1 8 20 Thickness of members that compose the frame cross section 0.5 2 5 Length of lower flange in stringer cross section 5 22 50 Length of upper flange in stringer cross section 5 15 50 Web height of stringer cross section 10 25 50 Lip height of stringer cross section 1 6 20 Thickness of members that compose the stringer cross section 0.5 2 5

Table 12. Design variables. Initial value and lateral limits. Values in millimeters

Figure 11. Design variables. Frame (left), stringer (right) and skin (bottom)

Constraints were related to structure components geometry and mechanical responses as follows: - Upper and lower limits of design variables as indicated in Table 12. - Maximum value of strain for the internal load case ελmin, being λmin 125. - Maximum value of stiffeners displacement to a value of 5% of maximum panel displacement under buckling. This constraint prevents global panel buckling and guarantees that the buckling occurs between two stiffeners. Figure 12 shows the global buckling mode compared to the local buckling mode if this constraint is entered.


Figure 12. Global buckling of panel (left). Local buckling of panel (right)

- Slenderness of each part of stiffener cross section according to regulations of Eurocode 9 for aluminum structures11. Slenderness parameter β is defined in the document as b  t

(24)

where b is the length of the member and t is the thickness. This parameter must satisfy the constraint  2    3 according to the values given in Table 13. β1

W/ thermal treatment, Element WO/ welding

W/ thermal treatment welding, WO/ welding WO/ treatment

WO/ thermal treatment, welding

W/ thermal treatment, WO/ welding

β2 W/ thermal treatment welding, WO/ welding WO/ treatment


W/ thermal treatment, WO/ welding

β3 W/ thermal treatment welding, WO/ welding WO/ treatment


Outer

3ε

2.5ε

2ε

4.5ε

4ε

3ε

6ε

5ε

4ε

Inner

11ε

9ε

7ε

16ε

13ε

11ε

22ε

18ε

15ε

ε=(250/f0)1/2 with f0 in N/mm2

Table 13. Slenderness parameters extracted from Eurocode 9

Recalling that the material is AL7075-T6 aluminum type with thermal treatment Armao12 and using Table 11 the value ε = 0.714 is obtained. Eurocode 9 considers two types of members, Outer elements are those having a free edge and the other braced and inner elements are those with both ends braced. After these considerations and taking in account that welding was not used in panel manufacturing the constraints could be posed as follows: -

Web buckling constraint (inner element).

LA  11.423  Le  15.714 11.423    15.714   11.423  TA  15.714  Te

(26)

where LA, Le, TA and Te are the design variables specified in Figure 8.The buckling mode prevented is shown in Figure 13. -

Upper flange buckling constraint (inner element). L1  11.423  Le  15.714 11.423    15.714   11.423  T 1  15.714  Te 12 American Institute of Aeronautics and Astronautics

(27)

where L1, Le, T1 and Te are the design variables specified in Figure 8, and the buckling mode that we want to avoid is shown in Figure 13. -

Lip buckling constraints (outer element). LB  3.214  Le  3.571 3.214    3.571    3.214  TB  3.571  Te

(28)

In these expressions where LA, Le, TA and Te are design variables specified in Figure 11 and the buckling modes to be avoided are shown in Figure 13.

Figure 13. Buckling modes of web (left), upper flange (center), lip (right)

The objective function in this case is the structural weight per unit area of skin. Using the geometric dimensions of the structure and density of the material, provided in Table 11, it is possible to express the mass of each structural element as follows: -

Mass of frames in kg/mm2

MT = 0.006496  Te  (2  TB+TA+T1- 2  Te+T2 - 2  Te) -

Mass of stringers in kg/mm2

ML = 0.01624  Le  (LB+ LA+ L1 - Le+ L2 - 2  Le) -

(29)

(30)

Mass of skin in kg/mm2

MP = 2.3548 Pe

(31)

The objective function expressed in kg/m2 is: F = 1 10 6

( MT  ML  MP )  1.189  ( MT  ML  MP ) 1450  580

(32)

To solve this problem the VisualDoc 6.213 code was used and structural responses of the model were obtained using the code Abaqus/CAE v6.9-2. A script written in Python was created using VisualScript 6.2 software that communicates the optimization code and analysis code. The iterative process converged after eight iterations and the evolution of the objective function is shown in Figure 14. This optimization process required a time of 3 hours using a computer with a Pentium 4 processor (3.4 GHz) and 2 GB of RAM. The optimum design variables are shown in the last column of Table 14. Moreover, a drawing of the comparison between optimal and initial design of the stiffeners is provided in Figures 15 and 16. The active constraints at the optimum design are: 1. 2. 3. 4. 5.

Load buckling factor constraint Maximum displacement of the stiffeners for buckling analysis Maximum strain in the skin for internal pressure Web buckling in frames and stringers Lip buckling in stringers


6. Lower limit of the variable T1

Figure 14. Objective function history Variable Initial Optimum Pe 1.5 1.7 T1 25 5.0 T2 15 9.7 TA 60 22.2 TB 8 2.7 Te 2 1.4 L1 22 5.2 L2 15 21.1 LA 25 21.7 LB 6 4.9 Le 2 1.4

Table 14. Optimum design. Dimensions in millimeters

Figure 15. Initial design of the frames (left in black color), optimum design (right in green color) and a comparison between both designs (center)


Figure 16. Initial design of the stiffeners (left in black color), optimum design (right in green color) and a comparison between both designs (center)

Figure 17 shows the buckling mode of the structure for the optimal values of design variables. The buckling mode for the optimal solution appears for a value of λ=125.02 and buckling instability occurs locally between two stiffeners, verifying that there is no global buckling.

Figure 17. Buckling mode for optimum design Sensitivity analysis and linear approximation of the optimal solution

The optimal solution to the above problem is the starting point for implementing the methodology to obtain approximate results of the optimal solution when varying a fixed parameter. Two fixed parameters were considered the minimum value of buckling load factor λmin =125 and the radius of curvature of the panel, r=1500 mm. The input data to obtain the approximate optimum design using the previous methodology are the set of active constraints, the gradients of the constraints and gradient of the objective function. In this case all of them are a byproduct of the optimization method itself, and therefore it is not necessary to obtain them again. To update all active constraints and implement the flow diagram of Figure 2, it will be necessary to calculate gradients in the new approximate designs. These values were obtained by finite differences.

p=λmin (∆p=-10%) The first approximate optimal solution was performed considering a 10% of reduction in the parameter p=λmin, from a value of λmin=125 to λmin=112.5. Initially, an approximate optimum solution was produced obtaining only once the feasible direction S for the total variation of parameter. The problem (6-9) was formulated and solved for this example, obtaining the S vector, which contains the 12 directions components associated to each design variable when the parameter p is modified. The resulting solution was: S= [-.003, 0, .3108, .0643, -.5599, .0041, .0411, -.0324, -.0452, -.0075, -.0021, -.7622] Knowing the parameter variation and applying the expressions (7-8), the new approximate optimum solution can be determined. Table 15 shows the result for p=λmin=112.5 (∆p=-10%) and also provides the optimal solution of the original optimization problem for the same parameter value.


Variable

p=λmin=112.5 (∆p=-10 %)

p=λmin=93.75 (∆p=-25 %)

X*

Pe 1.62 T1 5.01 T2 9.58 TA 22.55 TB 2.84 Te 1.44 L1 5.02 L2 21.06 LA 21.13 LB 4.78 Le 1.34 6.20 Mass (kg/m2) X*=optimum solution

 X

 X act

X*

 X

 X act

1.64 5.00 14.79 23.23 -6.51 1.48 5.87 20.55 21.00 4.75 1.35 6.12

1.64 5.00 13.23 21.67 1.00 1.38 5.09 20.64 21.35 4.70 1.33 6.22

1.53 5.00 9.49 21.32 1.87 1.37 6.49 18.43 21.48 4.02 1.37 5.86

1.58 5.00 22.44 24.84 -20.26 1.57 6.89 19.77 19.85 4.59 0.54 4.99

1.55 5.00 13.51 22.26 1.00 1.42 5.00 19.86 20.58 4.51 1.27 5.87

 =approximate solution without updating active constraints X  =approximate solution updating active constraints X act

Table 15. Optimum solution and approximate optimum design for p=λmin. Approximated results for ∆p=-10 % and ∆p=-25 %

However, the approximate design is invalid because the design variable TB has a negative value, not an appropriate value for a geometric dimension. This fact shows that 10% parameter variation cannot be done by solving in only one step the optimization problem (6-9). It will be necessary to update all active constraints and regain the vector S, according to the new methodology presented in the flowchart of Figure 2. Applying this procedure, the previous approximation is valid up to a parameter value p=122.72 (∆p=-1.8%), as at this value the constraint of minimal value for the variable TB is activated. This for this amount of parameter variation, a new feasible direction problem was performed giving the following values for the components of direction S: S= [-.0036, 0, .2464, -.0661, 0, -.0042, .0215, -.0326, -.0246, -.0141, -.0043, -.9657] Then, it was verified that no other constraint becomes active for the remainder of the parameter variation interval and the values of the approximate solution is given in Table 15. The results are very accurate, with an error of the objective function of only 0.25%.

p=λmin (∆p=-25%) The next variation on the parameter p=λmin is a reduction of 25%, from λmin=125 to λmin=93.75. Again, if the approximation of the optimal solution is done with a single evaluation of vector S, obtained for the point corresponding to λmin=125, the results would lead to a nonsense design with negative geometrical values. Thus, the approximation of the optimal solution was obtained by updating the design constraints, to find that for a parameter value λmin=108.62 (∆p=-13.1%) another constraint becomes active, which in this case is the lower limit of the variable SL1. Once again, it was necessary to obtain the vector S at this new point, resulting: S= [-.0039, 0, .2770, -.0844, 0, -.0054, 0, -.0333, -.0462, -.0077, -.0022, -.9554] With these data the approximate optimum solution could be obtained for p=λmin=93.75 (∆p=-25%) and the results are shown in Table 15. Looking at the numerical values can be seen that the error in the objective function is almost zero. Furthermore the percentage of violation of the constraints is below 2%. The approximate optimum solutions presented below have been obtained with the same computer and the computing time has not exceeded four minutes in any case.

p=r (∆p=-5%) The other parameter to which this methodology was applied is the radius of curvature of the panel. The first variation on this parameter is a reduction of 5%, from r=1500 mm to r=1425 mm. For the parameter c a value of c=1 turn out to be efficient. The result obtained is: S= [1.2e-4, 0, .0473, .0191, -.135, 1.23e-3, 7.42e-3, -.0156, 5.82e-3, 1.3e-3, 3.74e-4, -.99] This approximate solution obtained without updating active constraints is not valid because one of design variables is negative as shown in Table 16.


p=r=1425 (∆p=-5 %)

Variable

X*

Pe 1.68 T1 5.00 T2 6.35 TA 20.00 TB 1.00 Te 1.27 L1 5.00 L2 19.41 LA 22.59 LB 4.91 Le 1.44 6.29 Mass (kg/m2) X*=optimum solution

 X

 X act

1.70 5.00 13.27 23.62 -7.53 1.50 5.75 19.90 22.19 4.98 1.41 6.32

1.71 5.00 8.45 20.01 1.00 1.27 5.00 18.36 23.84 4.56 1.39 6.36

p=r=1350 (∆p=-10 %)

 X

 X act

1.71 5.00 16.86 25.07 -17.73 1.60 6.32 18.72 22.63 5.07 1.44 6.19

1.70 5.00 8.07 20.00 1.00 1.28 5.00 18.48 24.81 4.41 1.23 6.21

X* 1.73 5.00 5.00 20.00 1.00 1.27 5.00 18.45 18.52 4.21 1.18 6.09

 =approximate solution without updating active constraints X  =approximate solution updating active constraints X act

Table 16. Optimum solution and approximate optimum design for p=r. Approximated results for ∆p=-5 % and ∆p=-10 %

To overcome this error and to get a change of parameter value ∆p=-5%, it was necessary to solve the approximate problem four consecutive times due to activation of the lower limits of the variables TA, TB and L1. Table 17 shows the results of the direction S and the approximate design variables for the four parameter variations that alter the set of active constraints. The comparison between the optimum design and the approximate design for ∆p=-5% is provided in Table 16. Var. 1 2 3 4 5 6 7 8 9 10 11 12

∆p=-12.3 (-0.82%)  S X 1.20E-04 1.692 -1.17E-12 5.000 4.73E-02 10.278 1.91E-02 22.414 -1.35E-01 1.003 1.23E-03 1.427 7.42E-03 5.284 -1.56E-02 20.889 5.82E-03 21.819 1.30E-03 4.893 3.74E-04 1.389 -9.90E-01 -----

∆p=-20.6 (-1.37%)  S X 2.22E-04 1.693 -7.26E-04 4.994 -5.91E-02 9.783 -8.38E-03 22.344 -4.66E-12 1.003 -5.39E-04 1.422 -3.37E-02 5.002 -5.39E-02 20.437 2.51E-02 22.029 1.25E-02 4.998 1.62E-03 1.402 -9.97E-01 -----

∆p=-33.8 (-2.25%)  S X 1.83E-04 1.696 -9.69E-13 4.994 -4.88E-02 9.132 -1.75E-01 20.013 -9.69E-13 1.003 -1.12E-02 1.273 -1.31E-06 5.002 -8.70E-02 19.276 4.14E-02 22.581 -2.74E-03 4.961 2.66E-03 1.438 -9.79E-01 -----

∆p=-75 (-5%)  S X 2.74E-04 1.707 2.22E-03 5.001 -1.64E-02 8.448 2.24E-08 20.013 -9.90E-05 0.999 1.42E-09 1.273 3.06E-08 5.002 -2.19E-02 18.363 3.03E-02 23.843 -9.68E-03 4.558 -1.13E-03 1.391 -1.00E+00 -----

Table 17. Feasible directions S and approximate optimum design for each parameter variation that modifies the set of active constraints, from ∆p=0 to ∆p=-5%

The error in the objective function is only 1.1% and all constraints are satisfied. This example shows that preexisting methods may not be effective when the set of active constraints is not updated, even for small parameter variations. However, using the linear approximation presented in (14) to update the set of active constraints, the results improve substantially.

p=r (∆p=-10%) An interval of a 10% of parameter variation from r=1500 mm to r=1350 mm was also studied. In this case it was required to update the active constraints six times to approximate the optimal solution: four times from 0% to 5% and two additional times from -5% to -10%. Table 18 shows the result for the vector S when active constraints are updated and parameter variation ranges from -5% to -10%. The update of vector S was done not only when the set of active constraints changes, but also when the rate of violation of any constraint exceeded 3%. Table 16 shows the approximate solution compared to the optimum solution for a 10% reduction in radius of curvature. 17 American Institute of Aeronautics and Astronautics

The result presents an error of 2.05% in the objective function and all constraints are satisfied. ∆p=-115 (-7.67%)  S Var. X 1 -2.20E-04 1.698 2 1.50E-03 5.000 3 -7.55E-03 8.143 4 1.69E-03 20.000 5 4.23E-10 1.000 6 1.07E-04 1.277 7 1.36E-07 5.000 8 3.20E-03 18.493 9 -4.45E-03 23.663 10 8.49E-03 4.901 11 -4.96E-04 1.371 12 -9.98E-01 -----

∆p=-150 (-10%)  S X -2.07E-05 1.698 1.31E-03 5.000 -2.05E-03 8.070 2.12E-03 20.000 2.95E-14 1.000 1.18E-04 1.281 1.07E-02 5.000 -4.75E-04 18.476 3.23E-02 24.806 -1.40E-02 4.407 -4.00E-03 1.229 -9.98E-01 -----

Table 18. Feasible directions S and approximate optimum design for each parameter variation that modifies the set of active constraints, from ∆p=-5% to ∆p=-10%

V. Concluding remarks The following conclusions can be extracted from this work: a) Values of fixed parameters may change along the design process and this circumstance must be taken in account in optimization methodologies. b) Sensitivity analysis with regards to fixed parameters is a valuable technique that may avoid repeating costly optimization loops in presence of changes in parameter values. c) Although some methodologies already existed, an enhanced technique that tracks constraint behaviour has been defined and implemented. d) Three examples have been presented to test the efficiency of the proposed approach with regards to previous ones and the numerical results show some advantages for the new technique. e) The computer time required for the method presented to produce a very accurate approximate solution is about 1-3% of the computer time needed for a new optimization procedure.

Acknowledgments The research leading to these results has received funding from the Spanish Ministry of Science and Innovation under Project DPI2010-16238. Also, the Dutch National Aerospace Laboratory (NRL) is acknowledged for providing the finite element model of the panel.

References 1

Schmit, L.A., “Structural Design by Systematic Synthesis”, Second Conference on Electronic Computation, ASCE, Pittsburg, PA, pp. 105-132, 1960. 2 Sobieszczansky-Sobieski, J.; Barthelemy, J.F.M. and Riley, K.M., “Sensitivity of Optimum Solutions of Problem Parameters”, AIAA J., Vol. 20, No. 9, pp. 1291-1299, 1982. 3 Barthelemy, J.F.M. and Sobieszczansky-Sobieski, J., “Optimum Sensitivity Derivatives Nonlinear Programming”, AIAA J., Vol. 21, No. 6, pp. 913-915, 1983. 4 Barthelemy, J.F.M. and Sobieszczansky-Sobieski, J., “Extrapolation of Optimum Design Based on Sensitivity Derivatives”, AIAA J. Vol. 21, No. 5, pp 797-799, 1983. 5 Vanderplaats, G.N., “Numerical Optimization Techniques for Engineering Design”, Vanderplaats 3rd edition, 2001. Vanderplaats Research & Development Inc., Colorado Springs, USA. 6 Hernández, S.: “Métodos de Diseño Óptimo de Estructuras”, Colegio Oficial de Ingenieros de Caminos, Canales y Puertos, Madrid, 1990.


7 Vanderplaats, G.N. and Yoshida, H., “Efficient Calculation of Optimum Design Sensitivity”, AIAA J., Vol. 23, No. 11, pp. 1798-1803, 1985. 8 Abaqus/CAE v6.9-2, Inc. Abaqus/CAE User´s Manual, 2008. 9 Sun, C. T., “Mechanics of Aircraft Structures”, Wiley, 2nd Edition, June 2006. 10 Rouse, M; Ambur D. R.; Bodine J. and Dopker B., “Evaluation of a Composite Sandwich Fuselage Side Panel With Damage and Subjected to Internal Pressure”, NASA Technical Memorandum 11 Eurocode 9, Design of all aluminum structures, Part 1-1: General Common rules, BS EN 1999-1-1, 2007. 12 Armao, F., “On Aluminum Welding, MetalForming Magazine”, June 2000. 13 VisualDoc 6.2, “User’s Manual, 2009”. Vanderplaats Research & Development, Inc. Colorado Springs, CO. (http://www.vrand.com)
