Optim Eng (2009) 10: 109–132 DOI 10.1007/s11081-008-9044-4

Application of an iterative global approximation technique to structural optimizations Luca Lanzi · Alessandro Airoldi · Clive Chirwa

Received: 13 April 2007 / Accepted: 4 March 2008 / Published online: 5 April 2008 © Springer Science+Business Media, LLC 2008

Abstract An iterative global approximation technique based on the Kriging method is proposed. The technique is validated through analytical test cases and then applied to solve two practical optimization problems: the optimization of aluminium-foam filled absorbers against crashworthiness requirements and the optimization of composite stiffened panels against buckling and strength constraints. The absorbers of the first application consist of two co-axial aluminium alloy tubes filled with lightweight aluminium foam. They were optimized to collapse at a controlled force level and to be the lightest possible. Explicit Finite element analyses were performed to evaluate the structural behavior of the absorbers in the sample points used to build the approximation. In the second application stiffened panels were optimized against buckling and strength constraints. The Tsai-Wu criterion was used to estimate first-ply failures as strength limit of the structure. Non-linear Riks analyses were performed with ABAQUS/Standard to evaluate the shell behavior in the sample points used to build the response surfaces. Basing on the obtained results the proposed iterative procedure seems a promising alternative to the classic a-priori building of response surface allowing better accuracy and saving of sample points. Keywords Global approximation · Kriging method · Optimization · Aluminium foam · Composite panels · Crashworthiness · Buckling L. Lanzi Airbus, New Filton House, Filton, Bristol BS99 7AR, UK e-mail: [email protected] A. Airoldi () Dipartimento di Ingegneria Aerospaziale, Politecnico di Milano, Via La Masa 34, Milano 20156, Italy e-mail: [email protected] C. Chirwa Engineering and Design Department, The University of Bolton, Deane Road, Bolton BL3 5AB, UK e-mail: [email protected]

110

L. Lanzi et al.

1 Introduction Structural design and optimization involving non-linear phenomena such as postbuckling and/or crashworthiness are nowadays greatly limited due to the high CPU time and hardware resources they require. As a matter of fact, non-linearities are becoming of primary interest in structural design both in aerospace and automotive fields. This trend is proved to be true by the increasing efforts of the automotive industries at improving vehicle safety as well as by the meaningful efforts of aerospace companies to enhance the understanding of bird-strike phenomena, ditching, micro-drills of satellite, helicopter vertical crashes and so on. The response of structures under crash conditions is highly non-linear due to the experienced load peaks and the meaningful inertial forces which produce significant deformations, changing the boundary conditions resulting from contacts, failure of materials and/or structural joints. An appealing possibility to handle crashworthiness problems is indeed the formalization of optimization problems, such as trying to minimize the ratio between the absorbed energy and the absorber weight constraining the maximum force level. In any case, since crashworthiness requirements has to assure the safety of the occupants within a certain range of impact conditions, an actual design philosophy should account for several other requirements related to the occupant survivability as well as to the post impact integrity. Another example of the primary role played by the geometrical non-linearity is relevant to the buckling and post-buckling of composite shells. In this respect, it is well known and commonly accepted in the aeronautical field that an optimal design of aluminium alloy shells could be achieved allowing these structures to carry loads in the post-buckling range. The use of very efficient aluminium alloy structures is made possible by validated design procedure and by large database of experimental results. The weight reduction obtained by designing structures able to undergo loads higher than the first buckling one, could suggest the extension of this design philosophy also to composite materials. In both these fields of application, as well as in many others related to the structural design, the need to account for non-linear aspects with adequate accuracy levels requires numerical computations, such as detailed non-linear finite element analyses, with high computational time and costs. These aspects become of primary importance when optimization strategies want to be used solving a design problem. A promising strategy to be used in the preliminary design and sizing could be based on fast response surfaces capable to identify the most appealing design solutions containing as much as possible the number of finite element analyses. In fact expensive computations might be reduced in number and partially replaced using meta-modeling techniques, i.e. global approximation methods. A global approximation is a mathematical model aimed at capturing and approximating the behavior of an unknown function over the all domain of interest only knowing its values in a limited number of sample points. The challenge is then to define a reliable and fast global approximation using a set of sample points as limited as possible. Nowadays a variety of well consolidated global approximation methods exists and was successfully used in practical applications: Polynomial regressions (Mayers and

Application of an iterative global approximation technique

111

Montgomery 1995) Multivariate Adaptive Regression Splines, Least Square and/or Moving Least Squares (Levin 1998), Radial Basis Functions (Hussain et al. 2002), Artificial Neural Networks (Rafiq et al. 2001) and Kriging models (Chung and Alonso 2002) are but a few examples of the most common approaches. The choice between these approximation methods should, indeed, take into consideration not only their performance but also other factors related to their flexibility, their theoretical complexity and the effort required to build, use and update the approximation (Jin et al. 2001; Carpenter and Barthelemy 1993). In any case, the allocation of the sample points used to build the approximation turns out to be very important for the final performances of the approximation. In literature there are many criteria and schemes to allocate a-priori the sample points in a regular and convex domain of interest: Box-Behnken, Koshal, D-Optimality and Latincubes are but a few examples (Mayers and Montgomery 1995). Similar criteria have been formalized and proposed to account for irregular domains of interest as well as to deal with discrete variables (Lanzi et al. 2004). Despite of this wide number of criteria, it is practically impossible to say which is the best allocation method and which is the minimum number of sample points required to reach a desired accuracy when all the sample points are allocated a-priori. A promising solution, appeared only recently (Martin and Simpson 2002; Jurecka et al. 2005; Weiyu and Batilly 2002), is based on the idea of building global approximations through an iterative process avoiding the a-priori allocation of the sample points. The present work moves from this emerging idea to discuss and formalize a basic iterative global approximation procedure. Its theoretical aspects and background hypotheses will be discussed in Sect. 2. The proposed procedure will then be applied to solve two practical problems: the first one relevant to the crashworthiness performances of aluminium foam filled absorbers; the second one relevant to the post-buckling behavior of composite stiffened panels. 2 Assessment of iterative global approximations From a general standpoint, an iterative technique of building response surfaces is based on the idea that the accuracy of the approximation could be increased by adding a limited number of sample points ad hoc allocated within the domain of interest. Accordingly a method for building iterative global approximations starts with a first guess approximation, built using a limited number of sample points, and evolves through a series of new refined approximations obtained by adding few new sample points to the starting ones. A basic hypothesis in defining an iterative global approximation is that, at each new iteration the settlement of the new sample points can be driven only by information coming out from the values of the sample points available in the previous iteration. In fact, the behavior of the function to be approximated is completely unknown except than in the currently available sample points. Accordingly, a decision criterion must be defined to allocate the new sample points iteration by iteration. A first approach, to be used when the global approximation is part of an optimization procedure, could involve a decision criterion to resample the area that appear the most promising in terms of the objective function. The success of this approach de-

112

L. Lanzi et al.

Fig. 1 Schema of an iterative global approximation procedure

pends on the quality of the initial approximation. In fact, if the initial approximation is accurate enough, the process will lead to the identification of the global optimum in few iterations. On the contrary, if the quality of the initial approximation is not good enough, namely if the approximation model does not mirror with enough accuracy the behavior of the unknown function, the process might lead to the identification of a local optimum ignoring at all the presence of other and even better solutions. These limitations could be possibly overcome re-formulating the optimization process as a multi-objective one where objective functions and approximation error would be simultaneously minimized. Genetic and Evolutionary multi-objective algorithms (Deb 2005; Coello and Christiansen 2000) where generations of designs can be used to update the approximation models whilst searching for the optimal solutions seems particularly suitable for this purpose. A second approach, adopted in the following of this work, consists in the resampling the area where the approximation errors are higher. Clearly, this second approach needs an error estimator capable to foresee the local accuracy of the approximation knowing the values of the function to be approximated only in the current sample points. If a reliable error estimator exists, such an approach will ideally achieve a good accuracy overall the domain of interest correctly mirroring the behavior of the unknown function. As a result, when this approach is used within an optimization procedure, it will try to approximate the unknown function with the same accuracy level overall the domain increasing the chances to correctly identify global and local optima requiring as a drawback an higher number of sample points. The main advantage of this second approach is its intrinsic capability to optimize a-priori the accuracy of the responses surfaces with no regards on the later use of the response. Hence, this second approach allows to maximize the accuracy of the approximation given a minimum number of sample points without knowing if the approximation will be used later as objectives or constraints within a particular optimization process. The same global approximations can be used several times within different optimization runs without requiring any new sample point. In any case, an iterative global approximation method will present a flux diagram similar to the one of Fig. 1. The updating schemes could differ in terms of the adopted convergence and allocation criteria which could be problem dependent. Within an approach based on the error estimation, the updating of the global approximation could be stopped when the error estimator gives a result lower than a pre-

Application of an iterative global approximation technique

113

defined value for a certain number of consecutive iterations or, in a simpler way, when a maximum number of sample points has been allocated. However, the need of an error estimator must be accounted for in order to select a suitable approximation method among those nowadays available. The Kriging method (Chung and Alonso 2002; Jin et al. 2001; Weiyu and Batilly 2002) has been chosen in this work since it intrinsically defines an error estimator: the Mean Squared Error (MSE). The basic idea of the Kriging method is presented in, while Sect. 2.2 is focused on the formalization of the updating procedure and on the decision criterion used to select the new sample points. 2.1 Kriging method The Kriging method tries to approximate a given function y( x ) as the sum of a polyˆ where βˆ is the vector of the regression coefficients, and nomial regression f ( x , β), a stochastic process z( x ) that expresses a localized deviation from the polynomial regression: ˆ + z( y( ˆ x ) = f ( x , β) x ).

(1)

ˆ = F · βˆ where the matrix F The polynomial regression can be expressed as f ( x , β) contains the values of each polynomial function in the point x . Once given the N sample points xS , the stochastic process is built in such a way that: ˆ z( x ) = c( x )T · R−1 · ε(β)

(2)

where RN xN is the correlation matrix of the sample points based on a generic correˆ expresses the errors of the regression model f ( ˆ x , β) lation function r( x 1 , x2 ); ε(β) with respect to the known function values in the sample points, namely: ˆ = y − F · βˆ ε(β)

(3)

denoting with y the column vector of the values of the function to be approximated in the sample points. Finally, c( x ) is the correlation vector between any current point x and the N sample points xS . It is computed as: c( x ) = [r( x , xS1 ), r( x , xS2 ), . . . , r( x , xSm )].

(4)

Basing on the results presented in previous works (Weiyu and Batilly 2002; Sacks et al. 1989), it was decided to use a zero order approximation and a Gaussian correlation function RN xN , defined as follow: r( x 1 , x2 ) = e−

K

1 2 2 k=1 θk ·(xk −xk )

.

(5)

The Kriging approximation is completely defined when the regression coefficients βˆ and the correlation parameters θk , which determine the characteristics of the approximation itself, are assigned. Solving a classical minimum square estimation, the computation of the regression coefficients βˆ should be carried out as follow: βˆ = (FT · R−1 · F)−1 · FT · R−1 · F

(6)

114

L. Lanzi et al.

while the K values of θk can be identified by satisfying the so called ‘likelihood’ criterion, namely minimizing the function: 1

ψ(θk ) = det(R)− M ·

ˆ ˆ T · R · ε(β) ε(β) . M

(7)

More detailed descriptions on how regression coefficients βˆ and the correlation parameters θk can be computed are reported in (Weiyu and Batilly 2002; Sacks et al. 1989). One of the main advantages of the Kriging approximation is its intrinsic possibility to estimate the Mean Squared Error in any point x within the domain of interest as:      F 0 FT · x )] · MSE(x) = σ 2 · 1 − [FT , cT ( . (8) c( x) F R The availability of such an error estimator makes the Kriging approximation an ideal candidate to be used within an iterative global approximation schema. 2.2 Updating schema Exploiting the availability of the MSE coming directly from the Kriging method, the iterative updating procedure could be formalized so as to allocate new sample points where the MSE is the maximum. Following the scheme of Fig. 1, the analysis of the global approximation will consist in searching for the areas that exhibit the higher expected errors. The identified areas will then be re-sampled with new sample points. The searching of the new sample points could be performed and implemented in different ways basing on the problem needs and on the experience of the user. As far as these aspects are concerned, this work is mainly aimed at developing a simple but general procedure, underlining few future enhancement possibilities. A first hypothesis herein considered is that a single sample point could be computed at each iteration. This basically means that a serial process could be used to evaluate the sample points iteration by iteration, instead of a parallel one. This hypothesis, however, is not so restrictive and might be removed in different ways if a parallel cluster is available for the evaluation of the sample points as shown in (Martin and Simpson 2002). It seems thus natural to add the new sample point as closed as possible to the maximum MSE inside the overall domain of interest. Accordingly, the core task of the updating procedure would consist in the searching of the points with the maximum MSE computed by means of (8). Methods such as Genetic Algorithm (GA), stochastic and evolutionary techniques as well as direct search strategies seem the most promising in the identification of the regions with the maximum MSE. In fact these techniques will reduce the possibilities to be trapped in local MSE minima and will work properly even in presence of nonsmooth and discontinuous MSE. Besides, if it is true that these methods will require a higher number of evaluations of the approximation models with respect of gradient based methods, it is also true that a single evaluation of the approximation model is

Application of an iterative global approximation technique

115

usually negligible in terms of computational time. Accordingly whichever will be the strategy used for the identification of the maximum MSE points, the robustness of the search algorithm should be preferred to its efficiency with respect to the required evaluations of the approximation model. For this reason a Genetic Algorithm (GA) (Khoo and Chen 2001; Goldberg 1989) was used in this work to identify the sample points to be added in the updating of the approximation. In this respect, the choice of a robust even if quite simple GA with respect to other algorithms was mainly driven by the characteristics of the real word applications presented in the following paragraphs where the presence of high non-linear behaviors (due to crash performances for example) might have lead to non-smooth MSE especially in the first iterations of the iterative approximation. The DNA codification, the operators which define the GA, i.e. selection, crossover and mutation, as well as the convergence parameters could be eventually tuned for the specific problem basing on the number of design variables, the dimensions of the domain and the availability of CPU time and resources. Limiting to this work, the same GA parameters and operators were used for all the applications. In particular, a binary DNA codification with 8 bits for each design variable was used to potentially explore a full factorial grid of 256 nodes along each design dimension. An elitist selection was used to preserve the best optimal solution trough the iterations with cross-over and mutation probabilities fixed to 0.6 and 0.05, respectively. The searching for the optimal solution was stopped after a maximum number of 350 iterations or after 50 consecutives iterations without improvement of the objective function, namely the expected MSE. The number of members to be used in each iteration, denoted with N, was selected according to the number M of design variables defining the domain of interest as follow: N = 5 · 2M .

(9)

The above formalized procedure can handle one response surface at a time. As a matter of fact, two or more global approximations can be built simultaneously, as it is required in many structural applications, if the sample points of the different surfaces can be evaluated by the same structural analysis. This is for instance the case of a structural optimization problem in which the objective function and/or different constraints can be computed by the same finite element analysis for a given point in the domain of interest. More generally speaking, this happens every time two or more quantities to be approximated share the same domain of interest and can be evaluated by the same computation. In these cases the iterative procedure can be extended introducing a least square concept. Accordingly, a new MSE function can be defined as the weighted sum of the error estimators of each approximation. Once again, a single sample point is added at each new iteration where the value of this new defined fitness function is the maximum one. The whole iterative approximation procedure, as well as the genetic algorithm and the Kriging approximation algorithm were implemented in Matlab (Matlab 2002).

116

L. Lanzi et al.

3 Benchmark test cases The iterative procedure has been validated considering several analytical test cases based on known analytical functions. These benchmarks were carried out to compare the accuracy levels reached by the iterative procedure with respect to those reached by a standard a-priori allocation of the sample points. The accuracy evaluation is provided in terms of Maximum Absolute Error (MAE) and the Overall Mean Squared Error (MSE) defined as follow: MAE = max{abs(F ( x ) − f ( x )},

(10)

N 1  · MSE = (F ( xi ) − f ( xi ))2 N

(11)

x∈D

i

where F ( x ) are the exact values of the function to be approximated and f ( x ) are its approximated values. Aiming at a better assessment of the iterative procedure, two examples are reported in the following. Both the examples were defined on a 2-D domain of interest and an accurate evaluation of the Maximum Absolute Error and of the Average Mean Squared Error of the final obtained approximation has been performed on a reference regular grid of 150 × 150 points. 3.1 First test case The first test case compares the accuracy level of the proposed iterative procedure to the one obtained using a full factorial grid. The analytical function used to define this test case is defined in the unary domain 0 ≤ x ≤ 1, 0 ≤ y ≤ 1 as:    b·π F (x, y) = a · r(x, y) + sin √ · r(x, y) · [c · sin(d · ϑ(x, y))] (12) 2  where r(x, y) = x 2 + y 2 , θ (x, y) = arctg(x/y), a = 0.8, b = 2.4, c = 1.55 and d = 1.3. The Kriging approximation built by a regular grid of 5 × 5 sample points a-priori allocated is reported in Fig. 2 together with the analytical (exact) surface. The iterative process was initialized using four sample points located in the corners of the domain of interest, namely (0, 0), (0, 1), (1, 0) and (1, 1). It was stopped when the total number or 30 sample points was reached. The evolution of the response surface through the iterations and the allocated sample points are reported in Fig. 3. Table 1 summarizes the accuracies of the approximations. The results relevant to the MSE show that the iterative procedure allowed reaching a higher accuracy level with a saving in the number of sample points. In fact, the accuracy reached by the iterative procedure with 17 sample points overcomes the one obtained by the a-priori allocation of 25 points thus showing a 30% saving of sample points. Similarly, when the iterative procedure reached the same number of the sample points used for the a-priori built approximation, the iterative method led to an average absolute error of about two thirds and a MSE less than half those obtained by a regular grid of 5 × 5 points a-priori allocated.

Application of an iterative global approximation technique

117

Fig. 2 First test case: analytical surface of the first test function (a) and Kriging approximation based on a regular grid of 5 × 5 sample points (b)

Fig. 3 First test case: response surface evolution of the first test function through the iterations

As far as the convergence behavior of the updating schema is concerned, as discussed in Sect. 2.1, it is driven by the expected maximum MSE provided by the Kriging model. Figure 4 shows the Maximum Expected MSE iteration by iteration: this is the expected MSE in correspondence of which a new sample point will be

118

L. Lanzi et al.

Table 1 First test case: accuracy of the approximations Points

Max. abs. error

MSE

5 × 5 sample points allocated a-priori

25

0.3353

0.0934

Iterative procedure

17

0.2919

0.0602

25

0.2249

0.0423

Fig. 4 First test case: convergence curves of the iterative procedure for the first test function

allocated to update the approximation in the next iteration. Figure 4 reports also the reference MSE and the MAE evaluated on the grid of 150 × 150 points. The error graph of Fig. 4 underlines that: – the Maximum Expected MSE coming out from the Kriging Approximation according to (8), is in general greater than the reference MSE evaluated overall the domain of interest; – the iterative procedure exhibits a convergent behavior not only with respect to the average response, namely with respect to the MSE, but also in terms of local accuracy as proved by the decreasing trend of the Maximum Absolute Error; – however, the behavior of the approximation in terms of Maximum Expected MSE exhibits a non-monotonic trend characterized by irregular peaks. The peaks reflect sudden changes in the approximation model during the updating process. It is worth noting how the convergence of the approximation becomes more regular as the number of sample points increase, as well as the frequency of the peaks is reduced. 3.2 Second test case The second test case is focused on the use of the iterative procedure to built more than one function simultaneously. The iterative procedure is applied to approximate two different test functions defined in the bi-dimensional domain −2 ≤ x ≤ 2, −2 ≤ y ≤ 2: F1 (x, y) = a · (y − x 2 )2 + b · (1 − x)2 ,

(13)

Application of an iterative global approximation technique

119

F2 (x, y) = c · y 2 + d · x 2 − e · x · y

(14)

where a = 0.5, b = 0.005, c = 0.5, d = 3 and e = 1.5. The surfaces obtained by the iterative procedure are compared to the analytical ones in Fig. 5. Once again, the errors obtained by the iterative procedure are compared with those obtained by an a-priori allocation of 5 × 5 sample points in Table 2. The iterative procedure seems the most effective in terms of accuracy and reliability. In fact, using the iterative procedure, 19 sample points obtained a better accuracy reached by 5 × 5 sample points a-priori allocated. Similarly, if the number of sample points of the iterative procedure is increased up to 25, the errors become significantly lower than those experienced by the a-priori allocation. Figure 6a–c compares the error levels expected by the Kriging approximation with the true ones computed on a grid of 150 × 150 points. In particular, Fig. 6a reports the results relevant to F1 (x, y) while Fig. 6b reports those relevant to F2 (x, y). Also in this second case, the expected errors returned by the Kriging approximations are higher with respect to the true overall MSE proving the reliability of the error estimator.

Fig. 5 Second test case: analytical surface of F1 (x, y) (a), approximation of F1 (x, y) obtained by the iterative procedure with 19 sample points (b), analytical surface of F2 (x, y) (c) and approximation of F2 (x, y) obtained by the iterative procedure with 19 sample points (d) Table 2 Second test case: accuracy of the approximations for function F1 (x, y) and for function F2 (x, y) Sample points

5 × 5 sample points

F1 (x, y)

F2 (x, y)

Max. abs. error

MSE

Max. abs. error

MSE

25

0.0394

0.0179

7.92 × 10−4

3.72 × 10−4

19

0.0022

5.49 × 10−4

7.53 × 10−4

3.04 × 10−4

25

3.89 × 10−4

1.67 × 10−4

6.01 × 10−4

2.66 × 10−4

Allocated a-priori iterative procedure

120

L. Lanzi et al.

Fig. 6 Second test case: convergence curves and approximation errors on the first function (a) and on the second function (b)

3.3 Remarks The results of both test cases show that the iterative procedure performs better than an a-priori allocation of the sample points in terms of global accuracy (Mean Squared Error overall the domain of interest) and of local accuracy (Maximum Absolute Error). Another advantage of the iterative procedure is that it does not require any a priori estimation of the number of sample points needed to reach the desired accuracy level. More in general, the updating process should be stopped when the desired level of maximum MSE is reached. However, since the convergence behavior of the iterative procedure is general non-monotonic, a convergence criterion based on the maximum MSE should be satisfied for few consecutive iterations before the updating procedure could be stopped.

4 Applications Some theoretical aspects of iterative approximation techniques and their applications to some benchmark cases have been presented in the above paragraphs. In the following the proposed iterative procedure will be used to deal with real world structural optimization requiring quite expensive and detailed finite element analyses. Generally speaking, in order to apply an iterative global approximation technique to solve practical optimizations two distinct steps should be identified: in the first one the global approximations of the constraints and/or of the objective functions are built; in the second one the obtained approximations are used within an optimization algorithm to identify the best configurations. In the first phase the domain of interest will be analysed in a number of sample points as much limited as possible by means

Application of an iterative global approximation technique

121

Fig. 7 Typical configuration of the lightweight aluminium foam filled absorbers

of detailed finite element analyses used to iteratively build the global approximations accordingly to the scheme presented in Fig. 1. The computational efforts due to this first phase are directly related to the number of sample points to be analysed, i.e. the number of finite element analyses required to build the approximations. The second phase, namely the search for the optimal solutions, exploits the obtained response surfaces. From a general standpoint, the computational costs of the second phase should be significantly lower than those of the first phase because the structural responses are given by the already available approximations and no more expensive finite element analyses are required. In the following, two applications are presented: the optimization of the crashworthiness performances of aluminium foam filled absorbers and the optimization of composite stiffened panels accounting for post-buckling and strength constraints. 4.1 Optimization of aluminium foam tubes against crash requirements The following paragraph is focused on the application of the iterative procedure to the optimal design of energy absorbers for crashworthiness applications: aluminium tubes filled with aluminium lightweight foam. The absorber configuration consists of two co-axial alloy tubes filled by lightweight aluminium foam. A typical configuration of the absorbers is reported in Fig. 7. The method of using foam fillers to reinforce tubes began from the research activities of Thornthon et al. (Thornthon 1980) on polyurethane foam and has been more recently extended exploiting the properties of aluminium foams (Seitzberger et al. 1997). The crushing of aluminium foam proceeds through the formation, multiplication and propagation of localized failures through the specimen (Blazy et al. 2004). This process maintains the average stress level at an almost constant plateau for a wide range of average volumetric strain in the specimen. Thereafter the foam densification increases the stress up to significantly higher level. This typical behavior is

122

L. Lanzi et al.

Fig. 8 Stress-strain (volumetric) curve of the lightweight aluminium foam

Table 3 Design variables and domain of interest of the aluminium foam filled absorbers Minimum

Maximum

X1 —radius of the inner aluminium tube [mm]

20.0

40.0

X2 —thickness of the foam filler [mm]

10.0

40.0

X3 —thickness of the inner and outer aluminium tubes [mm]

0.50

2.00

reported in Fig. 8 where the stress vs. volumetric strain curve of aluminium foam is reported. The crushing phase corresponding to the almost constant stress level indicates that foam fillers can be used to reinforce aluminium tubes so to reduce the load oscillations usually experienced by hollow aluminium tubes and to obtain load-shortening curves characterized by more regular plateau. Furthermore, it is known that the synergy between the foam and the co-axial tubes increases the overall efficiency of the absorbing system in terms of their specifying absorbed energy, as discussed in previous works (Chirwa et al. 2005). 4.1.1 Domain of interest In the application presented in this work, the height of the absorber was fixed at 200 mm while the radius of the inner aluminium tube, the thickness of the foam filler and the thickness of the aluminium tubes were assumed as design variables, accordingly to Table 3. As far as the impact conditions are concerned, the absorbers were crushed with a constant impact velocity of 8 m/s up to a shortening of 150 mm corresponding to a stroke efficiency of 75% with respect to the their whole height. 4.1.2 Finite element model Finite element analyses were used to evaluate the structural behavior in the sample points during the iterative building of the global approximations. Further finite ele-

Application of an iterative global approximation technique

123

Fig. 9 Examples of load-stroke curves obtained by Ls-Dyna

ment analyses were performed at the end of each optimization to verify the approximation accuracy and the reliability of the obtained results. Finite element analyses were performed by using Ls-Dyna (Hallquist 2001). The outer and inner aluminium tubes were modeled by means of 4-nodes shell elements with an elastic-plastic constitutive law and a Cowper-Symonds law to account for the strain rate sensitivity of the plastic response. The aluminium foam was modeled by using constant stress solid elements and an incremental constitutive law based on the stress-volumetric strain curve for foam materials (Ls-Dyna MAT. 63). The total number of elements in the model depended on the values of the design variables. In particular, the finite element models were built so that the shell elements on the inner tubes have size of 3×3 mm about, and at least 5 bricks were used through the foam thickness constraining the aspect ratio of all the elements to be close to one. 4.1.3 Global approximations Three different global approximations were simultaneously built: mean force, maximum force and absorbed energy. The iterative procedure was started by using 4 different sample points in the corner of the domain of interest and was stopped after a total number of 35 sample points was reached. Figure 9 reports three different load-shortening curves obtained by Ls-Dyna during the iterative building of the approximations, while Fig. 10 shows the final position of the sample points evaluated by the iterative procedure. The maximum errors expected by the approximations, computed according to (8), are 14 kN on the maximum force, 8 kN on the mean force and 2.23 kJ on the absorbed energy. The total CPU time required to build the approximation, i.e. to run the 35 LsDyna analyses, was of 160 hours about on a common PC Pentium IV, 1500 MHz, 512 Mb of RAM memory. Neglecting the time required by the iterative procedure to re-build and analyse the approximation to find the sample point for the next iteration, an average time of 4 hours and 46 minutes for each Ls-Dyna run was required.

124

L. Lanzi et al.

Fig. 10 Final positions of the 35 allocated sample points

4.1.4 Optimization results In real world applications, once the maximum stroke has been fixed, absorbing devices should be designed to collapse at a controlled force level. The best absorber is the one that guarantees the desired force level with the minimum weight and with the lower ratio between maximum force and mean force. Accordingly, the optimization problem used to drive the absorber design can be formalized as follow: ⎧ min(Weight( x )), ⎪ ⎪ ⎨ X LOWER ≤ F UPPER , F¯MEAN x ) ≤ F¯MEAN MEAN ( ⎪ ⎪ ⎩ x ) ≤ F¯MAX . FMAX (

(15)

In order to show the capabilities of the proposed approximation technique, the opLOWER = 80 kN and timization problem was solved by fixing F¯MAX = 100 kN, F¯MEAN UPPER F¯MEAN = 85 kN. The searching for the optimum was performed by using a SQP algorithm, because of the continuity of the design variables as well as of the objective and of the constraint functions. During the optimization process, the evaluation of the constraints, i.e. FMAX ( x ) and FMEAN ( x ), was demanded to the global approximations while the evaluation of the objective function, namely the weight of the absorber, was provided analytically according to the following equation:  = π · [2 · ρAL · (2 · x1 + x2 ) · x3 + ρF OAM · (x32 + 2 · x1 · x3 )] · h W (X)

(16)

being ρAL and ρF OAM the densities of the aluminium alloy and of the foam filler, respectively, and denoting with h the total height of the absorbers.

Application of an iterative global approximation technique

125

Table 4 Characteristics of the optimal absorber configuration Constraint values

Approximation values

Ls-dyna Values

Abs. err.

% Err

Weight [kg]

–

0.543

0.543

–

–

Mean force [kN]

80–85

80.55

82.70

2.156

2.61

Maximum force [kN]

< 100

87.05

96.50

9.456

9.79

Absorbed energy [kJ]

–

12.63

12.89

0.2621

2.03

Limiting to the considered example, the optimization process required 16 iterations and 106 evaluations of the approximation to identify the optimal absorber configuration. Comparing these results with the 35 finite element computations required to build the approximation, it can be observed that the effort required to build the approximation was significantly lower than the one that would have been required to directly apply the optimization procedure using the finite element analyses. The identified configuration is characterized by a radius of the inner aluminium tube of 27.5 mm, a foam filler of 11 mm and the inner and outer tubes 1.85 mm thick. In order to validate the results obtained by the global approximations, a further LsDyna analysis was performed. Table 4 reports the main characteristics of the identified configuration as predicted by the global approximations and by the Ls-Dyna analysis. Results of Table 4 underline that the absolute errors between the approximations and the values foreseen by Ls-Dyne are within the expected tolerance limits. In fact, the error on the maximum force is 9.456 kN against the maximum expected value of 14 kN; the error on the mean force and on the absorbed energy are equal to 2.156 kN and 0.2621 kJ, against the maximum expected values of 8 kN and 2.23 kJ, respectively. Figure 11 shows the behavior of the optimal absorber configuration in terms of its load-stroke curve and deformed shape. As proved in Fig. 11, the optimized absorber exhibits a stable and regular crash front assuring good absorption capabilities with an almost constant force level avoiding initial peaks loads. In this respect, the initial optimization requirements seem completely satisfied proving the good formulation of the optimization problem. 4.2 Optimization of composite stiffened panels against buckling and strength constraints The iterative procedure was then applied to optimize composite stiffened panels against buckling and strength constraints. The optimization and the design of composite stiffened shells under post-buckling and strength constraints is widely discussed in aerospace literature and it represents, probably, one of the emerging challenges for the extensive industrial application of composite materials in primary aircraft structures. An interesting general discussion of the engineering aspects of the problem can be found in the review provided by Nemeth (Nemeth and Starnes 1998) as well as other more recent works (Arbocz and

126

L. Lanzi et al.

Fig. 11 Load-stroke curve and deformation of the optimal absorber configuration

Table 5 Composite stiffened panels: material properties used in the finite element models

Young’s modulus E 11 [N/mm2 ]

58615

Young’s modulus E 22 [N/mm2 ]

58615

Shear modulus G12 = G23 = G13 [N/mm2 ]

3064

Poisson ratio ν 12

0.05

Density ρ [kg/m3 ]

1510

Nominal Ply thickness [mm]

0.32

Starnes 2002). Optimization strategy have been recognized as a promising approach to this problem since the late nineties, as in the works by Wiggenraad et al. and Nagendra et al. (Wiggenraad et al. 1998; Nagendra et al. 1996) while surface response techniques have been more recently applied by Rikards et al. (Rikards et al. 2004). However, the efforts of this paper are devoted more in the assessment of an iterative global approximation procedure rather than that in the investigation of composite structures under post-buckling requirements. The provided references are suggested for a better understanding of optimization problems relevant to the design of composite stiffened panels. In the considered optimization problem, after a limit load was defined, the composite panel was designed to exhibit a first buckling beyond a given load level that could be lower than the limit load, exploiting the post-buckling response of the panel. A stable post-buckling range, without failure, was also required until the ultimate load, defined as 1.5 times the limit one. 4.2.1 Domain of interest The baseline configuration of the composite stiffened panels is shown in Fig. 12. It consists of flat panels with two T shaped stringers located on the free edges and having a web thickness double of the foot one. The panel was supposed to be manufactured laminating woven fabric plies of carbon fibre reinforced plastic (CFRP) with an almost equal behavior in 0◦ and 90◦ directions, as reported in Table 5. The structure was loaded under compression and the lateral edges are free. The length of the panel was fixed at 650 mm and the distance between the stringer was

Application of an iterative global approximation technique

127

Fig. 12 Baseline configuration of the composite stiffened panels

Table 6 Design variables and domain of interest of the composite stiffened panels Minimum

Maximum

X1 —semi-foot of the T shaped stringers [mm]

15.0

25.0

X2 —height of the stringer webs [mm]

15.0

35.0

X3 —panel skin thickness [mm]

1.0

3.5

X4 —stringer foot thickness [mm]

1.0

3.0

0.30

0.70

0.30

0.70

X5 —thickness percentage of the 45◦ oriented layers in the stringers X6 —thickness percentage of the 45◦ oriented layers in the skin

fixed at 250 mm. During the optimization 6 design variables were considered according to Table 6. The optimization problem took in account both sizing aspect, such as the foot and the web of stringer, as well as the lamination sequences. With respect to the lamination sequences, a homogenization approach was followed. Thus, the percentage of 45◦ oriented layers was considered as a continuous design variable while the relative position of the layers in the laminate has been fixed so that 45◦ layers remained in the middle while 0◦ layers were stacked on the laminate surfaces. As far as the stacking sequences of the stringer are concerned, they were supposed to be manufactured by two co-cured L-shaped beams. Accordingly, their foot has three homogenized layers [0◦ /45◦ /0◦ ] oriented while the web has [0◦ /45◦ /0◦ ]S homogenized layers. The overall thicknesses of the foot and the web of the stringer were considered as design variables. 4.2.2 Finite element model Also in this second application, finite element analyses were used to build the response surfaces needed to evaluate the constraints during the optimization process. A further finite element analysis was performed at the end of the optimization to verify the approximation accuracy and the reliability of the obtained results. Finite element analyses were performed by using ABAQUS/Standard (Abaqus Theory and Users’ manuals, 1998). Following the results obtained in previous works (Lanzi 2004), full non-linear analyses were performed by using the arch-length method proposed by Riks (1972, 1979). Stiffened panels were modeled by 4-nodes

128

L. Lanzi et al.

shells, reduced integration, with six degrees of freedom at each node. Three integration points along the thickness for each homogenized ply were used. After a preliminary sensitivity study on the element size and basing on the experience gained in previous works (Lanzi 2004), the dimensions of shell elements were assumed of about 8 × 8 mm. Aiming at the evaluation of the strength capabilities of the investigated structure, a first ply failure approach was used. Even if this kind of approach might be considered conservative to completely exploit the post-buckling capabilities of composite stiffened panels, it is commonly used in real structure design. The Tsai-Wu criterion (Jones 1999) was adopted to localize first ply failures experienced by the panel skin and stringers in post-buckling. Accordingly to this criterion a failure happens if: 

1 1 − C T σ¯ 11 σ¯ 11



 · σ11 +

1 1 − C T σ¯ 22 σ¯ 22

+ 2 · F12 · σ11 · σ22 > 1

 · σ22 +

2 σ11 T σ¯ 11

C · σ¯ 11

+

2 σ22 T σ¯ 22

C · σ¯ 22

+

2 σ12 2 σ¯ 12

(17)

where σij are the current stress components as experienced by the structure while σ¯ ij are the material strength properties and the interaction term F12 is assumed equal to −0.5 for the considered woven fabric. 4.2.3 Global approximations In order to account both first buckling and first-ply failure loads, two global approximations were simultaneously built: critical load and first-ply failure load. The iterative procedure was started by using 10 points inside the domain of interest and stopped when a total number of 150 sample points had been allocated. The maximum errors expected by the approximations, computed according to (8), are 2 kN on the first buckling load and 20 kN on first-ply failure one. The total time required to build the approximation, namely to allocate and compute 150 sample points, was of 23 hours about on a common PC Pentium IV, 1500 MHz, 512 Mb of RAM memory, corresponding to an average time of about 9 minutes for each sample point. Eight minutes were used by ABAQUS/Standard to run the nonlinear analysis with the modified Riks’ method, half minute was used to post-process the load curve and the stress values so to extract the first-ply failure index. The remaining time, half minute about, was required by the iterative process to re-build and analyse the obtained approximations. 4.2.4 Optimization results As it usually happens in the design of aeronautical stiffened panels, the optimization problem was formulated to minimize the weight of the panel and to guarantee minimum buckling and first-ply failure loads. Supposing that the structure can overcome the first buckling load in a restricted number of operating flight conditions, the minimum first buckling load was fixed at the 0.65 of the limit load. Similarly, the structure was required to carry loads up to the ultimate one without first-ply failures

Application of an iterative global approximation technique

129

Table 7 Characteristics of the optimized stiffened panel configuration Constraint values

Approximation values

ABAQUS/Standard Values

Abs. err.

Weight [kg]

–

1.040

1.040

–

First buckling load [kN]

45

45

45.2

0.2

First-ply failure load [kJ]

100

100

105.2

5.2

localized in the stringers. The ultimate load value is fixed at 1.5 times the limit load. Accordingly, the optimization problem was formulated as follow: ⎧ min(Weight( x )), ⎪ ⎪ ⎨ X (18) x ) ≥ F¯BUCKLING ≈ 0.65 · F¯LIMIT , FBUCKLING ( ⎪ ⎪ ⎩ x ) ≥ F¯PLY−FAILURE ≈ 1.50 · F¯LIMIT . FPLY−FAILURE ( In order to show the capabilities of the proposed approach, the optimization problem was solved fixing F¯BUCKLING = 45 kN and F¯PLY−FAILURE = 100 kN corresponding to a limit load of about 67 kN. Since the design variables, the objective function as well as the constraint equations were continuous inside the domain of interest, the searching for the optimum x ) and was performed by using a SQP algorithm. Constraint values, i.e. FBUCKLING ( x ), were evaluated by the global approximations while the objective FPLY−FAILURE ( function, namely the weight of the stiffened panel, was provided analytically according to the following equation:  = ρ · [4 · x1 · (x3 + x4 ) + 4 · x2 · x4 + (d − 2 · x1 ) · x3 ] · L W (X)

(19)

being ρ the density of the composite fabric and denoting with L the free length of the panels, i.e. 650 mm, and with d the distance between the web of the stringers, i.e. 250 mm. Limiting to this case, the optimization algorithm required 40 iterations and 301 function evaluations to identify the optimal panel configuration. The identified configuration is characterized by stringers having a foot 30 mm width and 1.2 mm thick corresponding to webs 30.4 mm height and 2.4 mm thick. The percentage of 45◦ layers in the stringer is of 0.3, lying on the lower bound of the domain of interest as expected by common engineering feeling. The panel skin is 3 mm thick with a fraction of 45◦ oriented layers equal to 0.67. The results obtained by the global approximations were validated by a further finite element analysis. Table 7 compares the main characteristics of the identified configuration as foreseen by the global approximations against those computed by the ABAQUS/Standard analysis. Results underline again how the errors between the approximations and the finite element analysis are within the expected tolerance limits. Figure 13 shows the behavior of the optimal configuration in terms of its loadshortening curve and its deformed shape. As required by the formulation of the optimization problem, the optimized panel configuration exhibits a local skin buckling

130

L. Lanzi et al.

Fig. 13 Load-shortening curve and first-ply failure of the optimized stiffed panel

between the stringers and a wide large and stable post-buckling range. At about 80 kN the webs of the stringers experience local buckling phenomena. However, the equilibrium path remains stable up to first ply failures occurring in the stringer webs at a load of about 105 kN.

5 Conclusive remarks This work investigates and formalizes an iterative procedure to build and update global approximations exploiting the availability of the MSE in the Kriging method as error estimator. The iterative procedure was validated considering several analytical test cases. Two of them have been reported and discussed showing that the iterative procedure performs better than an a-priori allocation of the sample points in terms of global accuracy (Mean Squared Error overall the domain of interest) and of local accuracy (Maximum Absolute Error). The validation tests allowed also to critically review the non-monotonic convergence behavior of the procedure. An advantage of the iterative procedure is that it does not require any a priori estimation of the number of sample points needed to reach the desired accuracy level. More in general, the updating process can be stopped when the desired level of maximum MSE is reached and maintained for few consecutive iterations. The iterative approximation has then been used to solve two different optimizations problems relevant to real world applications: the optimization of aluminium foam filled absorbers against crashworthiness requirements and the optimization of composite stiffened panels accounting for their post-buckling behavior and strength. The absorbers of the first application were optimized to collapse at a controlled force level and to be the lightest possible. Finite element analyses were performed with Ls-Dyna to evaluate the structural behavior of the absorbers in the sample points used to build the approximation. Three different design variables were considered and the iterative procedure allowed containing the total number of sample points to 35 reaching an accuracy level adequate to perform reliable optimization searches. In fact, the optimized absorber configuration satisfies all the imposed constraints and

Application of an iterative global approximation technique

131

exhibits a stable and regular crash front assuring good absorption capabilities with an almost constant force level avoiding initial peaks loads. In the second application stiffened panels were optimized against buckling and strength constraints. The Tsai-Wu criterion was used to estimate fist-ply failures as strength limit of the structure. Non-linear Riks analyses were performed with ABAQUS/Standard to evaluate the panel behavior in the sample points. The optimized panel configuration exhibits a local skin buckling between the stringers and a wide large stable post-buckling range as required by the formulation of the optimization problem. First ply failures occurred in the stringer webs as expected by common engineering feeling. Concluding, the obtained results seem to prove how the proposed procedure could be a promising alternative to the classic a-priori building of response surface allowing better accuracy and further saving of sample points so to potentially reduce the total computational efforts required by optimization involving high-fidelity computational analyses. Furthermore, the feasibility and flexibility of the proposed global approximation techniques could be exploited to express and solve different optimization problems, changing objectives, constraint functions and values among the structural responses available as global approximations. In fact, once built, the same global approximations can be used several times within different optimization runs without requiring any new finite element analyses.

References Abaqus (1998) Abaqus® theory and users’ manuals. Hibbitt, Karlsson & Sorensen, Pawtucket Arbocz J, Starnes JH Jr (2002) Future directions and challenges in shell stability analysis. Thin-Walled Struct 40:729–754 Blazy JS, Marie-Louise A, Forest S, Chastel Y, Pineau A, Awadec A, Grolleron C, Moussy F (2004) Deformation and fracture of aluminium foams under proportional and non proportional multi-axial loading: statistical analysis and size effect. Mech Sci 46:217–244 Carpenter WC, Barthelemy JFM (1993) A comparison of polynomial approximations and artificial neural nets as response surfaces. Struct Optim 4:166–174 Chirwa EC, Lehmhus D, Mao M, Chen T, Lanzi L (2005) Mechanics of lightweight aluminium foam wrapped in carbon fibre reinforced composites. In Alves M, Jones N (eds) Impact loading of lightweight structures. WIT Press Chung HS, Alonso JJ (2002) In: 9th AIAA/ISSMO symposium on multidisciplinary analysis and optimization, Atlanta, USA, 4–6 September 2002. AIAA, Washington, pp 1–21 Coello CA, Christiansen AD (2000) Multiobjective optimization of trusses using genetic algorithms. Comput Struct 75:647–660 Deb K (2005) A population-based algorithm-generator for real-parameter optimization. Soft Comput 9:236–253 Goldberg DE (1989) Genetic algorithms in search, optimization and machine learning. London, AddisonWesley Hallquist JO (2001) LS-DYNA keyword users manual, version 960, vols I and II. Livermore Software Technology Corporation Hussain MF, Burton RR, Joshi SB (2002) Metamodeling: radial basis functions, versus polynomials. Eur J Oper Res 138:142–154 Jin R, Chen W, Simpson TW (2001) Comparative studies of metamodelling techniques under multiple modeling criteria. Int J Struct Multidiscipl Optim 23:1–13 Jones RM (1999) Mechanics of composite materials. Tailor & Francis, London

132

L. Lanzi et al.

Jurecka F, Bletzinger KU, Ganser M Update schemes for global approximations in robust design optimization. In: Proceedings of the 6th world congress of structural and multidisciplinary optimization, Rio de Janeiro, Brazil, 30 May–03 June 2005, pp 1–10 Khoo LP, Chen CH (2001) Integration of response surface methodology with genetic algorithms. Int J Adv Manuf Technol 18:483–489 Lanzi L (2004) A numerical and experimental investigation on composite stiffened panels into postbuckling. Int J Thin-Walled Struct 42:1645–1664 Lanzi L, Bisagni C, Ricci S (2004) Neural network systems to reproduce crash behavior of structural components. Comput Struct 82:93–108 Levin D (1998) The approximation power of moving least-squares. Int J Math Comput 67:1517–1531 Martin JD, Simpson TW (2002) Use of adaptive metamodeling for design optimization. In: Proceedings of the 9th AIAA/ISSMO symposium on multidisciplinary analysis and optimization, Atlanta, 4–6 September 2002. AIAA, Washington, pp 1–9 Matlab (2002) Matlab® user manuals, Matlab 6.3, version 5. The Math Works, Natick Mayers RH, Montgomery DC (1995) Response surface methodology: process and product optimization using designed experiments. Wiley, New York Nagendra S, Jestin D, Gürdal Z, Haftka RT, Watson LT (1996) Improved genetic algorithm for the design of stiffened composite panels. Comput Struct 58:543–555 Nemeth MP, Starnes JH Jr The NASA monographs on shell stability design recommendations: a review and suggested improvements. NASA TP-1998-206290 Rafiq MY, Bugmann G, Easterbrook DJ (2001) Neural network design for engineering applications. Comput Struct 21(17):1541–1552 Rikards R, Abramovich H, Auzins J, Korjakins A, Ozolinsh O, Kalnins K, Green T (2004) Surrogate models for optimum design of stiffened composite shells. Compos Struct 63:243–251 Riks E (1972) The application of Newton’s method to the problem of elastic stability. J Appl Mech 39:1060–1065 Riks E (1979) An incremental approach to the solution of snapping and buckling problems. Int J Solid Struct 15:529–551 Sacks J, Welch WJ, Mitchell TJ, Wynn HP (1989) Design and analysis of computer experiments. Stat Sci 4:409–435 Seitzberger M, Rammerstorfer FG, Degischer HP, Gradinger R (1997) Crushing of axially compressed steel tubes filled with aluminium foam. Acta Mech 125:93–105 Thornthon PH (1980) Energy absorption by foam-filled structures. SAE paper 800081 Weiyu L, Batilly SM (2002) Gradient-enhanced response surface approximations using Krgiging models. In: 9th AIAA/ISSMO symposium on multidisciplinary analysis and optimization, Atlanta, USA, 4–6 September 2002. AIAA, Washington, pp 1–11 Wiggenraad JFM, Arendsen P, da Silva Pereira JM (1998) Design optimization of stiffened composite panels with buckling and damage tolerance constraints. AIAA-98-1750, pp 420–430